1. Introduction

The development of diagnostic imaging techniques commonly entails the use of “phantoms”, which try to mimic human tissue as closely as possible. They are conceived in such a way that their optical properties are similar to those of tissues and, therefore, constitute a precious tool which provides a well-defined environment to study light propagation. Quoting a recent comprehensive review by Pogue and Patterson [1], these phantoms “are used for a number of purposes, including: (i) initially testing system designs, (ii) optimizing signal to noise in existing systems, (iii) performing routine quality control, (iv) comparing performance between systems.”

Our group is conducting research in optoacoustic (OA) imaging [2] where phantoms are widely used based on hydrogels (such as gelatin, agar, polyvinylalcohol) or PVC plastisol and other resins (see [1,3–5] and references therein). We opted for gelatin phantoms because they are relatively easy to produce and have the advantage of being adjustable to various purposes: (i) the phantom's optical properties can be customized by incorporating scattering agents (like TiO₂ or intralipid) or absorbing elements (like India ink, or any dye of interest), and (ii) the phantom's mechanical/elastic properties can be modified by varying the fraction of gelatin. Elastic properties especially play a valuable role: indeed, the combination of OA imaging with ultrasound elastography holds promise for reduced clutter and increased imaging depth in the framework of clinical imaging [6,7].

Our gelatin phantoms are produced in a two-step approach. First, a phantom precursor stock is prepared by dissolving gelatin in water and by incorporating titanium dioxide (TiO₂) as the scattering component. During the second step, different phantoms can be developed from this precursor stock by adding various amounts of supplementary scattering agents, absorbing structures, and acoustic scatterers (such as flour) to obtain more realistic optical and acoustic properties for combined OA and ultrasound imaging.

Accurate knowledge of the optical properties of the phantom precursor stock is essential within this process. Determining the optical properties of such phantoms or tissues is not an obvious task and has been subject to intensive research over the past decades (see e.g. the review by Kim and Wilson [8]). In this paper, we employ a series of methods by which we were able to fully characterize the optical properties of the precursor stock (which is simply referred to as “phantom” or “sample” throughout this manuscript) at a wavelength of 780 nm. This wavelength is particularly appropriate since it lies within the diagnostic window where, owing to low optical absorption, OA imaging can achieve several centimeters of imaging depth [9].

The way this paper is organized is outlined in the diagram of Fig. 1 . After describing how the gelatin phantom was produced, we describe how the optical properties of the phantom were determined at 780 nm, as follows:

 

Fig. 1 Diagram summarizing the different stages of the work presented in this paper: the optical properties (refractive index, transport mean free path, absorption coefficient, scattering coefficient, and scattering law) of the phantom were experimentally determined by the indicated set of experiments. Monte Carlo simulations were performed ex-post to validate the experimental results.

Download Full Size | PDF

In the last part we present Monte Carlo (MC) simulations performed to validate the experimental results: we simulated the TOF experiments by providing the experimentally determined optical properties of the phantom as inputs to the MC program.

2. Gelatin phantom preparation

2.1. Ingredients and preparation of the phantom

The phantom precursor mixture consists of 16.7% weight fraction of gelatin dissolved in water and 1 gram per liter (g L⁻¹) of TiO₂ for optical scattering. The TiO₂ concentration was chosen according to the data in [4], which suggest a reduced scattering coefficient of 0.4 mm⁻¹ at 1 g L⁻¹. This is below the range 0.5-1.3 mm⁻¹ typically found in human breast tissue [10,11], but we anticipate the need for further additives, including e.g., flour for acoustic contrast, which will modify the scattering properties.

The TiO₂ particles were obtained from Millenium Chemicals (type Tiona RL11A). They consist of 99% rutile, with a refractive index of 2.4884 at 780 nm [12]. Their size distribution lies between 180 nm and 400 nm diameter, with a clear peak at 250 nm (weight percentage, determined by disc centrifugation).

The phantom in its precursor form is prepared as follows:

2.2. Refractive index

For the evaluation of TOF measurements (Section 3) and for the MC simulations (Section 6), it is necessary to know the refractive index n of the phantom at the wavelength of$\lambda =780nm$. After having prepared a phantom sample, we measured at the room temperature of 22°C both its refractive index $n_D\approx 1.363$ at the sodium d-line wavelength ${\lambda }_D=589.29nm$ and its corresponding Abbe number $\upsilon =1/(n_F-n_C)=56.72$, where $n_F$ and $n_C$ are the refractive indices at the wavelengths ${\lambda }_F=486.13nm$ and ${\lambda }_C=656.27nm$. These two values were then used to determine the two parameters A and B of the two-term form of the Cauchy extrapolation formula $n(\lambda )=A+B/{\lambda }^2$. This law gives the refractive index of the phantom at 780 nm: $n=1.359\pm 0.002$. The measurements were done with an Abbe refractometer (Type 1T by ATAGO). Since the refractrometer is calibrated with water, we use the error estimate of the IAPWS formulation for pure water in the wavelength range 700-1100 nm [13].

3. Determination of the transport mean free path and absorption length

3.1. Principle of the measurements

Our technique is a variant of position and time resolved diffuse reflectometry [14], but we place photon injection and detection deep inside the scattering medium, in order to minimize the influence of the boundaries [15–17]. The principle is illustrated in Fig. 2. Photons generated by a pulsed picosecond diode laser (LDH-P-785, PicoQuant) are delivered into the phantom through a multimode optical fiber (50 mm core diameter) and injected at a point S well inside the scattering medium. A second fiber collects the scattered light at a point D at a distance r from S and guides it to a single photon detector (id100-50, ID Quantique). The time of flight (TOF) of the individual photons from S to D is measured by a fast timer (SPC 530 PicoQuant) triggered by the laser pulse. The resulting TOF distribution is characteristic of the optical properties of the material through which the light propagates. Thus, by fitting an appropriate analytical model to the recorded distribution curve, the optical properties of the phantom can be extracted.

3.2. Diffusion model

We employ the photon diffusion model [18], which is an approximation to describe the light transport in strongly scattering media (i.e., where ${\mu }_s>>{\mu }_a$) and over large distances such that $r>>{\ell }^{*}$. According to this approximation, the probability $P(r,t)d{r}^{3}dt$ to detect a photon in a volume element $d{r}^3$ and time interval dt obeys the diffusion equation:

where $D_p$is the photon diffusion coefficient and ${\alpha }_p$ is the absorption rate. $D_p$ and ${\alpha }_p$are related to the phantom’s transport mean free path ${\ell }^{*}$and absorption length ${\ell }_a$ respectively:

Here $c_l$ is the light speed in the medium. (In the biomedical photonics community one often uses the reduced scattering coefficient ${{\mu }^{\prime }}_s=1/{\ell }^{*}$and${\mu }_a=1/{\ell }_a$. Note that we assume ${\mu }_s>>{\mu }_a$.) The Green’s solution of Eq. (1) is

This simple solution applies in infinite space and for a δ-source in space and time: $S(r,t)=\delta (r)\delta (t)$. To approximate the former condition, we employ a large sample container, so that photons are injected and detected far from boundaries. To regard S(r) as a point source, we require that$ds\ll {\ell }^{*}$, where ds is the characteristic size of the source. This condition is easily met by using an optical fiber as an illuminator.

The time dependence of the measured distribution M(t) results from the convolution of P(t) with the instrument response function IRF(t) that includes the finite width and shape of the laser pulse and the timing jitter of the detection: $M=IRF\otimes P$. To determine the least square estimates of the parameters D_p and α_p, we apply the standard iterative reconvolution technique, using the Levenberg – Marquardt algorithm provided by Matlab. A third parameter is needed to account for the transition from the ballistic photon transport to diffusion that occurs within approximately a distance of 2${\ell }^{*}$ from the injection point. We include in the fitting procedure a (manually) variable shift τ_d of the measured TOF-distribution with respect to the measured IRF. This shift does not play a significant role if $r>>{\ell }^{*}$. Since ${\ell }^{*}$is not known a priori, we must carry out a series of measurements at varying source-detector separations. These estimated values should converge to the “true” values as r increases, so that most of the photon paths are long enough to assure the validity of the diffusion approximation.

3.3. Setup and measurement procedure

The phantom is contained in a soda-lime glass aquarium ($180mm\times 240mm\times 210mm$). The setup was designed to perform measurements at different source-detector separations without deforming either the phantom or the fibers (see Fig. 3 ). The multimode fibers (core diameter of 50 mm), both on the laser and the detector side, are guided through hypodermic needles (diameter 0.8 mm) in order to reach the inner parts of the gelatin phantom with minimal damage and better accuracy. Up to eight such optodes can be positioned in the orifices of a holder-plate and fixed perpendicularly to the phantom’s top surface. The heads of the needles are placed on the holder so that their tips are deep inside the phantom. This allows both the injection and the collection of the light far from the phantom-air boundary (typically at a depth of 40 mm beneath the top surface). To avoid distortions caused by the insertion of the needles, the phantom was poured in its liquid state into the aquarium where the holder and fibers had been positioned beforehand; the phantom then solidified with the needles already in place. The accuracy of the positioning of the fiber tips was inspected by taking photographs prior and after the solidification of the gelatin: we found a maximum deviation of $0.6\text{\hspace{0.17em}mm}$ from the intended separation.

 

Fig. 3 Detail of the setup used for the TOF measurements. (a) The gelatin phantom (5 L) is contained in a soda-lime glass aquarium. (b) and (c) In the holder-plate, the optodes can be placed at eight positions. Each pair of optode/orifice on the holder-plate corresponds to a different source-detector separation r. The tips of the needles are inside the phantom, at a depth of 40 mm beneath the phantom/air interface.

Download Full Size | PDF

Each pair of fiber/orifice on the holder plate corresponds to a different source-detector separation r. Thus, measurements at different source-detector separations are carried out by successively selecting different pairs of fibers and connecting each fiber to the laser and the single photon counting module. Prior to each measurement the IRF corresponding to the chosen pair of optodes was determined by measuring the time-resolved scattering from a sheet of white paper immersed in water. With the present laser and photon detector the width of the IRF is small when compared with the width of the TOF distribution M(t). Therefore, the precise shape of the IRF does not play a crucial role. What matters, however, is the timing, i.e., the shift between IRF(t) and M(t). Strictly speaking, for a precise determination of the IRF, one should place the fiber tips close together so that the photon path length is shorter than $1{\ell }^{*}$. However, this exact positioning is not easy to achieve (for example, we must insert filters to adjust the photon-count rate), and, therefore, in the data evaluation we shall rely on the adjustable parameter${\tau }_d$to compensate the resulting inaccuracy of timing. An example of the best fit of the measured data by the diffusion model is shown in Fig. 4 .

 

Fig. 4 Example of the fit of the measured data (thin solid line in black) with the diffusion model (thicker line in red). The corresponding IRF is shown in green dots. For comparison, we also include the fit with the corresponding MC simulations (see Section 6) with the same parameters ${\ell }^{*}$ and ${\ell }_a$ (thicker line in blue).

Download Full Size | PDF

3.4. Results and discussion

In Fig. 5 we compile the results of five series of measurements on five different phantoms. First of all we observe that ${\ell }^{*}$ is rather large, larger than the transport mean free path in human breast (typically around 1 mm [10,11]). Correspondingly large is the influence of inaccuracies and limitations of the diffusion approximation, which results in the data scatter at small distances r. From an analysis of the influence of primary experimental errors on the least square estimates of the parameters ${\ell }^{*}$ and ${\ell }_a$ we found that the overall variance is composed of three contributions:

 

Fig. 5 The values of the transport mean free path (top) and absorption length (bottom) obtained from the TOF measurements with five different phantoms. The vertical dashed line is a “guide to the eye” indicating roughly the validity limit of the diffusion approximation.

Download Full Size | PDF

Nevertheless, the estimates of ${\ell }^{*}$ of the five phantoms converge with a reasonable consistency to the values ${\ell }^{*}=4.5\pm 0.5mm$ and ${\ell }_a=190\pm 15mm$, as r increases beyond approximately six times${\ell }^{*}$.

The transport mean free path of almost 5 mm (reduced scattering coefficient ${{\mu }^{\prime }}_s$ of only 0.2 mm⁻¹) is substantially larger than the typical values reported for elastomer based phantoms with the same TiO₂ concentrations (DMSO 0.9 mm [19], PVCP 2.5 mm [4], Polyurethane 1.25 mm [20]). On the other hand, much larger values for ${\ell }^{*}$ have been reported for a comparable gelatin based system: Pogue et al. [21] measured ${{\mu }^{\prime }}_s(c)$ of TiO₂ in 7% gelatin (at $\lambda =690\text{\hspace{0.17em}nm}$) and found a linear concentration dependence ${{\mu }^{\prime }}_s(c)\approx 0.14\cdot c$ (${{\mu }^{\prime }}_s$ in mm⁻¹, c in g L⁻¹) up to c = 40 g L⁻¹. Their data suggest that ${\ell }^{*}$ is as large as 7 mm for 1 g L⁻¹ TiO₂ (unfortunately, the particle size is not specified).

4. Scattering coefficient and asymmetry

4.1. Extinction measurements

The extinction coefficient ${\mu }_e$ of the phantom was determined at 780 nm using a Perkin Elmer LAMBDA^TM 750 absorption spectrometer. With the absorption coefficient ${\mu }_a=1/{\ell }_a$ as small as 0.005 mm⁻¹, we expect that most of the extinction is due to scattering, i.e., ${\mu }_e={\mu }_s$. If the concentration of the TiO₂ solution is low enough that multiple scattering can be neglected, then the transmission follows Beer-Lambert’s law$I=I_0\exp [-d{\mu }_s(c)]$, where d is the thickness of the sample, $I_0$ is the intensity of the illuminating light and I the intensity of the transmitted light. Furthermore, ${\mu }_s(c)$ depends on the concentration c of the TiO₂ particles. In the simplest case, i.e., when concentration dependent clustering of titania particles can be neglected, ${\mu }_s$is proportional to the TiO₂ concentration c in the sample:

where ${\sigma }_{s\_TiO2}$ is the scattering cross-section. Thus, the absorbance $A=-{\log }_{10}(I/I_0)$ is expected to be linear in d and c: $A={\log }_{10}(e){\sigma }_{s}cd$. The experimental data shown in Fig. 6 show that A is indeed linear in both, the sample thickness d and concentration c, which indicates that we indeed measured in the single scattering regime and that the cross-section of the individual scattering entities does not change with the concentration. Moreover, the linear regime characteristic of the Beer-Lambert’s law extends up to the concentration c = 1 g L⁻¹, which we used to prepare the phantom.

 

Fig. 6 The absorbance A of TiO₂ solutions (containing 16.7% weight percentage of gelatin) as the function of concentration c (left: sample thickness d = 0.1mm; right: d = 1mm).

Download Full Size | PDF

The scattering coefficient at a TiO₂ concentration of 1 g L⁻¹ as obtained from the linear fit to the data in Fig. 6 is${\mu }_s(c=1g{L}^{-1})=2.86\pm 0.04m{m}^{-1}$. The error estimate represents the standard statistical error of the least square parameter estimation.

4.2. Asymmetry parameter g

Knowing the scattering coefficient μ_s, we can determine the phantom’s asymmetry parameter $g=〈\cos (\theta )〉$, i.e., the average of the cosine of the scattering angle θ, which is related to the transport mean free path and the scattering coefficient according to

Using ${\ell }^{*}$ = 4.5 ± 0.5 mm, and ${\mu }_s=2.86$ ± 0.04 mm⁻¹, we obtain $g=0.93\pm 0.01$. This large asymmetry indicates substantial forward scattering, which is much more pronounced than one would expect from Mie scattering on TiO₂ particles. The largest g one can get with isotropic spherical particles of the given refractive index ratio would be about 0.7 (regardless of the size of the particles!). Moreover, finite element calculations of Thiele and French [22] indicate that the effect of non-spherical shape and anisotropy is not sufficiently large to explain the large asymmetry. Apparently, the scattering entities are not bare isolated TiO₂ particles. We speculate that the effective scatterers are some large (though sparse) TiO₂-gelatin clusters formed by micro-phase separation during the solidification of the system.

5. Assessing the scattering law

The ultimate characterization of a phantom tissue would consist of the complete determination of the scattering law $p(\theta ,\varphi )$ (also called “phase function”) characterizing the angular distribution of scattered photons in a single scattering event. With samples exhibiting strongly asymmetric, nearly forward, scattering, this is a rather difficult task, which requires specialized instrumentation. Here we report preliminary results obtained with a commercial goniometer originally intended for dynamic light scattering measurement (ALV GmbH Langen). The available hardware does not allow a precise determination of $p(\theta )$ at small angles, but is sufficient to check the consistency of the results obtained so far.

5.1. Goniometric measurements

The geometry of the scattering experiment is sketched in Fig. 7 . The light source was a 780 nm CW diode laser (Laser Components). The scattered photons were collected with a multimode fiber receiver focused at the center of the sample and counted with a single photon detector. The sample is immersed in a toluene filled cylindrical vat, to minimize geometrical errors resulting from refractive index mismatch. A polarizer-analyzer pair selects the SS-component (perpendicular to the scattering plane) of the phase function.

 

Fig. 7 Top view of the light scattering experiment. The light emitted by the laser is vertically polarized. We measured the light scattered by TiO₂ suspensions (containing 16.7% weight percentage of gelatin) of different concentrations at various angles 40°≤θ≤140°.

Download Full Size | PDF

An accurate determination of the single scattering law of the strongly scattering sample would require an extrapolation to a sample cell whose dimensions are smaller than $1/{\mu }_s$, i.e., smaller than 0.4 mm. This is not practicable with the given setup and even less with a heterogeneous sample. However, the observed linearity of ${\mu }_s(c)$ and ${{\mu }^{\prime }}_s(c)$ (see [21]) indicates that the structure of the scattering entities does not depend on the TiO₂ concentration. This suggests a simpler alternative: We carried out measurements on TiO₂ dispersions (each containing 16.7% gelatin) with decreasing concentrations, and from this series we extrapolated the trend of the scattering law. Examples of the extrapolation procedure are indicated in the inset in Fig. 8 .

 

Fig. 8 Estimating the scattering law of the TiO₂ scatterers immersed in gelatin: the pvgHG scattering law model is fitted to the extrapolated data (represented by “+” symbols in the main figure) corresponding to infinitesimal concentrations of TiO₂, i.e., where the single scattering regime applies. The fitting is achieved by modifying the set of parameters {g, ν}. The extrapolation procedure which gives the values corresponding to infinitesimal concentrations is indicated in the inset figure.

Download Full Size | PDF

For each TiO₂-gelatin dispersion, the scattered intensity $I_{meas}(\theta )$ was measured at different angles${40}^{\text{o}}\le \theta \le 140°$. (The apparatus is not optimized for forward scattering and therefore we restrict the measurements to scattering angles larger than 40°). The measured signal $I_{meas}(\theta )$ is normalized, so that $I_n(\theta )$ is obtained as follows:

Here $I_{blank}(\theta )$ is the blank signal, measured with water in place of the TiO₂-gelatin sample and $I_{tol}(\theta )$ is the reference signal measured with toluene. Since $I(\theta )\propto {\mu }_{s}p(\theta )$, the normalized signal can be expressed as:

${\mu }_{s\_tol}$and $p_{tol}(\theta )$ are the scattering coefficient and the scattering law of the toluene, respectively and ${\mu }_s$(c) is the concentration-dependent scattering coefficient of the phantom as previously determined. Toluene scatters as an almost perfectly isotropic dipole scatterer (Rayleigh scattering), thus:$p_{tol}(\theta )=3/8\pi$. (The theoretical background of polarized Rayleigh scattering is discussed for example in [23].) Thus, the scattering law of the TiO₂ suspension can be determined according to

Moreover, ${\mu }_{s\_tol}$is identical to the so-called Rayleigh ratio $R_{tol}$of the toluene. From the data provided by Moreels et al. [24], we determine the Rayleigh extrapolation formula $R_{tol}(\lambda )=1.65\times {10}^8/{\lambda }^4$ (R_tol in m⁻¹ and λ in nm), which leads to ${\mu }_{s\_tol}=4.46\times {10}^{-4}{\text{m}}^{\text{-1}}$. In Fig. 8 the extrapolated values are indicated with symbols “+”. Note, that since ${\mu }_{s\_tol}$, $p_{tol}(\theta )$ and ${\mu }_s$ are known, an adequate model should match not only the shape but also the magnitude of the experimental curve without adjusting a normalization pre-factor.

5.2. Modeling of the scattering law

As the first guess, an obvious model would be the Mie scattering law corresponding to nearly spherical TiO₂ particles with a certain size distribution. However, as already mentioned, disregarding the size of the particles, the largest value of g that can be obtained with TiO₂ dispersed in the gelatin gel is 0.7. This is much smaller than the previously determined g = 0.93. Thus, we need another scattering law to model the proposed clustering of the TiO₂ in the gelatin. To account for the interference between the individual scatterers forming a cluster, we resort to the Rayleigh-Debye approximation. (Which is also suited to describe light scattering in tissue in a more realistic fashion than Mie scattering. See Rička et al. [23].) In the Rayleigh-Debye approximation one combines polarized scattering from elementary Rayleigh scatterers with a structure factor that characterizes their mutual arrangement. We employ the polarized version of the generalized Henyey-Greenstein scattering law (pvgHG) (discussed in more details in Refs. [23,25]), which assumes the following “fractal” structure factor F:

Here q is the scattering vector,$|q{|}^2=2k^2[1-\cos (\theta )]$, where k is the wavenumber in the medium and θ the scattering angle. The pvgHG structure factor has two modifiable parameters: the correlation length L characterizes the extent of correlation between the elementary scatterers, i.e., roughly the size of the proposed clusters. The power law exponent ν is related to the shape of the correlation, i.e., roughly the shape of the clusters. The case $\nu =0$ corresponds to the well-known Henyey-Greenstein phase function widely used in the field of tissue scattering. (Note, however, that we also include polarization. For more on these topics see Section 7.3 in [23].) The correlation length L is not yet known, but we already have an estimate of the asymmetry parameter g. Therefore, for the convenience of the data analysis, we numerically convert the parameters $\{L,\nu \}$ into $\{g,\nu \}$.

Thus, the pvgHG model is fitted to the experimental data by manually modifying the set of parameters$\{g,\nu \}$. (Recall that$\{g,\nu \}$are the only adjustable parameters; there is no adjustable pre-factor.) The examples in Fig. 8 indicate the sensitivity to the parameter g. The model with g = 0.90 and ν = 0 (curve 1) fails to give an accurate description of the data but with g = 0.93 (curve 2) is more satisfactory. A further improvement is achieved by varying the parameter ν. A reasonably good fit is obtained with $\{g=0.93,\nu =-0.3\}$ (curve 3).

From g = 0.93 we estimate a correlation length L in the order of 10 μm. The exact value depends on the parameter ν, but the data quality does not allow its unambiguous determination. Nevertheless, the scattering law is consistent with the previous results from Sections 3 and 4.

6. MC simulations

In the last stage of the investigation, MC simulations were performed to assess the robustness of the experimental results. We employ a simulation tool that we recently developed [25] to study polarized light propagation in complex geometries. While polarization is not important in the diffusion limit, we shall take advantage of the following features of our software:

i) A user interface allows creation of complex shapes in the simulation space through which the ray tracing algorithm follows the paths of photons across the interfaces and through the different regions. ii) The user can choose the optical properties for different regions of the sample, including the interface properties (e.g. Fresnel reflections) and scattering laws. (Currently Mie and pvgHG are implemented.) These features allow transposing the TOF experiment into a realistic simulation (compare Fig. 9 and Fig. 2). The optical properties of the phantom and its container that were given as inputs to the simulations are listed in Table 1 .

 

Fig. 9 Screenshot from the MC simulation: the simulation tool was used to reproduce the TOF measurements presented in Section 3 (depicted in Fig. 2). The colored disc visualizes the number of photons detected by the individual detector rings.

Download Full Size | PDF

Table 1. The optical properties (at λ = 780 nm) assigned to the three different materials used in the experiment: the quantities listed here correspond to the experimental results and were provided as input parameters to the MC simulation program

View Table | View all tables in this article

To simplify the simulation design and to speed up the computation we permit ourselves certain simplifications of photon launching and detection, which should not affect the results in the diffusion limit. Specifically, we neglect the presence of the needles, because they are made of non-absorbing material and thus create only a minor disturbance.

Photons are simply launched downwards at a position corresponding to the tip of the source fiber, 40mm beneath the phantom’s top surface. The individual photons are launched with randomly chosen initial polarization (the so-called unpolarized light) to mimic the depolarization of the laser light occurring in a multimode fiber.

The detection is simulated as a concentric array of detectors (r resolution of 1mm), which register incoming photons without absorbing them (see Fig. 9). We used two versions of this detection scheme. In the first version, a photon is registered when it crosses the horizontal plane of the receiving fiber tip within a numerical aperture NA = 0.87. The fibers used during the measurements have NA = 0.2, but increasing the NA reduces the simulation time. The second detection scheme registers a photon when a scattering event occurs within a narrow slab (thickness 1 mm) adjacent to the detector plane, again into the aperture of 0.87. (To compare the resulting TOF distributions for different scattering laws, they must be normalized by dividing with μ_s.) This second detection scheme provides approximately a six fold improvement of counting statistics but is only appropriate for distances within the validity of the diffusion approximation.

All the simulations presented in this paper have been performed on a Linux PC equipped with an Intel Pentium D computer processor with a clock speed of 3.40 GHz and a cache size of 2048 KB. The path generation routine written in C + + has been compiled with a GNU compiler (gcc version 4.3.2).

6.1. Testing the diffusion limit, influence of scattering law

In a first series of simulations we tested whether the setup indeed fulfills the conditions for the applicability of the diffusion approximation. The results for two optode distances r are shown in Fig. 10 . The transport mean free path and absorption length are set to the previously determined values${\ell }_{input}^{*}=4.5\text{\hspace{0.17em}mm}$,${\ell }_{a\_input}=190\text{\hspace{0.17em}mm}$, but we vary the scattering law (see inset in the bottom of Fig. 10) and its parameters. Obviously, at $r=5\text{\hspace{0.17em}mm}$ the diffusion approximation fails to describe the simulated data (which is not surprising, see e.g [26].), but for 36 mm distance (8${\ell }^{*}$) we are already well within the diffusion regime. Moreover, at $r=36mm$ there is no way to distinguish the scattering law used in the simulation, since ${\ell }^{*}$ and ${\ell }_a$ are the only parameters that matter in the diffusion limit. However, the influence of the scattering law becomes apparent at an optode separation of about 1${\ell }^{*}$, despite the fact that all simulations are done with the same ${\ell }^{*}$ and ${\ell }_a$. Apparently, the asymmetry g is important on the length scale ${\ell }^{*}$ and on the time scale that corresponds to a mean photon path length of 1${\ell }^{*}$ . The asymmetry g as well as the form of the scattering law affect angularly resolved reflectance experiments [27,28].

 

Fig. 10 Examples of the MC simulations where we varied the scattering law and its parameters. The curves in the inset of the bottom figure are polarization averaged. The data were generated in an overnight simulation including 60 million photons for each scattering law (total number of scattering events: 47 × 10⁹). The simulated curves show the raw, non-normalized and unsmoothed data with a time-resolution of 2 ps. At small source-detector separations (top), the limitation of the diffusion model (red line) is clearly visible but the differences between the scattering laws are hardly discernable in the semi-log plot. These differences are significant at short times, as shown in the linear plot in the inset of the top figure. (The small bump at 25 ps is probably a commensurability artifact due to the finite resolution of the ring detector.) At large source-detector separations (bottom), the applicability of the diffusion model becomes apparent. At both distances, the diffusion model is shifted by τ_d ≈50ps.

Download Full Size | PDF

6.2. Analyzing the simulated data

In Fig. 4 the outcome of the simulation (pvgHG law {$g=0.93,\nu =-0.3$}) is compared with the measured data obtained at the optode separation $r=36\text{\hspace{0.17em}mm}$. The simulation is convoluted with the IRF of the experiment, i.e., we regard the simulation as a model function M(t), an alternative to the diffusion model. Both models match the measurement within the statistical error. However, there is a slight mismatch between the diffusion model and the simulation model with the same parameters. This could indicate that at $r=36\text{\hspace{0.17em}mm}$ we are not yet within the range of the validity of the diffusion approximation.

To further elucidate this issue, we analyzed the simulation data by fitting the diffusion model. Note that thereby we regard the simulation data as equivalent to a measurement, the only difference from a real measurement is that we have now a perfect estimate of the IRF, namely a δ-peak positioned at t = 0. The fit parameters are the transport mean free path ${\ell }^{*}$, absorption length ${\ell }_a$ and the shift ${\tau }_d$ that accounts for the transition from ballistic to diffusive regime of photon propagation. Examples of the results for four optode distances are summarized in Table 2 . The estimates of the parameter uncertainty are derived from visual inspection of the quality of the fit. (Statistical analysis in terms of χ² would not be appropriate, since we expect a systematic discrepancy between the diffusion model and the data for short photon paths.)

Table 2. Parameters obtained by fitting the diffusion model to the simulated data with${\ell }_{input}^{*}=4.5mm,{\ell }_{a\_input}=190mm$

View Table | View all tables in this article

The optimal shift ${\tau }_d=50\text{\hspace{0.17em}ps}$, corresponding to 11 mm of the photon path, is consistently found for all optode separations. Thus, the ballistic regime of photon propagation extends to about 2${\ell }^{*}$. Note, that the uncertainty of the parameter ${\tau }_d$ increases with the optode distance. This shows that the shift ceases to be important as the spacing enters the range of validity of the diffusion approximation. The parameter ${\ell }^{*}$ converges within the given uncertainty to the input value ${\ell }_{input}^{*}=4.5\text{\hspace{0.17em}mm}$. However, the parameter ${\ell }_a$ converges to a value which is significantly smaller than the input value ${\ell }_{a\_input}=190\text{\hspace{0.17em}mm}$. This indicates the increasing influence of the boundaries on this parameter. Apparently, the optode depth of 40 mm below the air interface is too small, causing a cut-off of long photon paths. Thus, the ${\ell }_a$ measured with the real phantom is slightly overestimated (which we take into account in the final error estimate in Section 7, Conclusions). However, the resulting error of about 10 mm remains well within the experimental error estimate.

7. Conclusions

The optical parameters of our phantom tissue precursor consisting of 1 g L⁻¹ TiO₂ in 16.7% gelatin are summarized as follows:

The transport mean free path ${\ell }^{*}$of almost 5 mm (reduced scattering coefficient ${{\mu }^{\prime }}_s$ of only 0.2 mm⁻¹) is substantially larger than initially expected for the given concentration of TiO₂ particles. This large ${\ell }^{*}$reflects the large asymmetry of underlying scattering law. From the measured ${{\mu }^{\prime }}_s$and ${\mu }_s$we obtain g = 0.93( ± 0.01), which is much larger than one would expect for Mie scattering from bare TiO₂ particles (g = 0.5-0.7). Indeed, our preliminary measurements with a simple goniometer indicate that the angular dependence of scattering is not compatible with the TiO₂ particles, but can be reasonably modeled with a generalized Henyey-Greenstein structure factor (Eq. (9) with g = 0.93 and shape parameter $\nu =-0.3$. We speculate that the scattering is due to sparse microstructures in the order of 10 μm in size.

These findings must be taken into account in the next steps of the phantom development. When designing a tissue phantom, one should not blindly expect the optical properties of the product to be the sum of the optical properties of its components. The optical properties depend on the microstructure of the system, as it results from interactions between the constituents, and whose formation depends on the physico-chemical parameters such as pH and ionic strength. Certainly, the microstructure of the gelatin-TiO₂ phantoms would deserve a thorough investigation. (This, however, would require a dedicated small angle light scattering setup.)

The parameters ${\ell }^{*}$and ${\ell }_a$ (which are the relevant parameters in many practical situations) were determined using a novel version of the fiber-optic TOF technique in which we immerse the injection and detection optodes deep inside the scattering medium. This geometry minimizes the influence of the interfaces on the TOF distribution and allows the analysis of the data in terms of the simple free space diffusion approximation.

The robustness of the approach was tested using our new simulation tool, which allows not only realistic modeling of the experimental geometry but also provides a simple way to vary the scattering law and its parameters. In the diffusion regime, the parameters obtained from simulated and experimentally determined data agree well within their error estimates. The simple diffusion approximation can be reliably applied down to optode separations of less than 6${\ell }^{*}$, if one corrects for a small time shift ${\tau }_d\approx 2{\ell }^{*}/c_l$ (in our case about 50 ps) that accounts for the transition from ballistic to diffusive regime of light propagation.

Moreover, our simulations confirm that the above four parameters, namely n, ${\ell }^{*}$, ${\ell }_a$ and μ_e, represent a necessary and sufficient set of optical parameters of a realistic tissue phantom. Usually n, ${\ell }^{*}$, ${\ell }_a$ are the most relevant parameters. The fourth parameter μ_e (or rather the asymmetry parameter g that corresponds to a certain μ_e when ${\ell }^{*}$ is fixed) plays a role only at length and time scales that involve a small number of scattering events. (This may be the case, for example, with multiphoton imaging of the skin.) In any case, the exact shape of the scattering law seems to play a negligible role in most practical situations (although there may be exceptions [27]). This certainly simplifies the design and characterization of the phantom tissue.

These insights combined with the experimental and computational tools described here create a powerful platform on which to build realistic tissue phantoms. Such optically, acoustically, and geometrically realistic phantoms will serve as essential tools in the development of important new biomedical imaging techniques.