## Abstract

Over the past decade, developments towards near-infrared (NIR) optical tomography involve the recovery of interior optical maps from boundary measurements using the first principles of light propagation models. The refractive-index mismatch parameter in the boundary condition of the light propagation model, namely the diffusion equation, can significantly impact model prediction of measurements and therefore image recovery. In this contribution, the influence of refractive-index mismatch parameter between predictions and referenced measurements of fluorescence-enhanced frequency-domain photon migration (FDPM) are established; its greater influence on emission over excitation predictions are demonstrated, and the methods to accurately determine refractive index mismatch parameter from basic principles are reviewed.

© 2002 Optical Society of America

## 1. Introduction

Near-infrared (NIR) optical tomography used for recovery of endogenous- or exogenous-contrasted tissue regions depends upon accurate exterior surface measurements of diffusely propagated light, an accurate model describing the transport of light, as well as an optimization approach to recover the interior optical properties. The diffuse propagation of light in random media can be predicted by the diffusion approximation of the radiative transport equation [1–3] under the conditions of predominant scattering over absorption. Various boundary conditions are used to represent the light propagation at the surface of the medium [4–10]. Many early works fail to identify the significant impact of refractive index-mismatch at the boundary when determining optical properties for spectroscopy or imaging applications [11,12]. More recently, several investigators have highlighted the significance of this index-mismatch parameter for demonstrating optical tomography in phantom studies [13,14]. Bartlett *et al*. [13] considered the refractive index-mismatch at the boundary interface to be an unknown parameter and employed error-minimization methods to determine the relative refractive index, *n _{rel}* of the medium. In another method, Pogue

*et al*. [14] suitably chose the index-mismatch parameter to minimize the model mismatch between experiments and simulations. In both cases, the index-mismatch parameter in the boundary condition was not accounted from basic principles of diffuse reflection of unpolarized light. Instead, the relative refractive index was used as an added parameter estimate.

Since the relative refractive index can be determined directly, other investigators have derived expressions from basic principles to determine the refractive index mismatch parameter at boundary surfaces [4,15,16]. However, the refractive index that is generally employed involves the air-phantom interface; whereas most phantom studies involve a solution held in a container made of plastic or any other material. Consequently, the relative refractive index for a container-phantom interface, which differs from the air-phantom interface, must be employed to mimic the actual phantom. Reflections through the container may need to be considered depending on the material of the container [8,17]. Without proper consideration of this index mismatch, erroneous model predictions result in greater model-mismatch between measurement and simulations as well as contribute to the ill posedness of the inverse imaging problem.

In this paper, the methods for accurately determining the refractive index mismatch at interfaces are reviewed and demonstrated in a phantom study using measurements of fluorescence-enhanced frequency-domain photon migration (FDPM). The significance of the proper choice of index-mismatch parameter at the boundaries in the fluorescence measurements is also highlighted.

## 2. Background: Boundary conditions for the diffusion equation

Fluorescence-enhanced optical imaging involves the use of exogenous fluorescing contrast agents to increase the optical contrast between the normal and diseased tissues [18,19]. At the Photon Migration Laboratory (PML), the measurement technique employed is the frequency-domain using fluorescing agents, and hence the technique is termed as fluorescence-enhanced FDPM.

The propagation of excitation light and the generation and subsequent propagation of fluorescent light in random media are given by a set of coupled diffusion equations, which predict the excitation and the emission fluence [20,22].

Three boundary conditions, namely i) partial current boundary condition [4,8] ii) extrapolated boundary condition [5,6,10] and iii) zero-boundary condition [6,7] are commonly employed. However, the partial current boundary condition, which states that the photon leaving the tissue surface never returns, is the most rigorous and widely used of the three [4,8]. At excitation (subscript x) and emission (subscript m) wavelengths (785 and 830 nm, respectively), the partial current or Robin boundary condition is given by

Here, *D _{x,m}*(

*r*) is the optical diffusion coefficient at position

*r*(cm),

*Φ*(

_{x,m}*r*,

*ω*) is the complex fluence describing the wave located at position

*r*and modulated at frequency,

*ω*, and

*γ*is the index-mismatch parameter, which is a function of the effective refractive index (

*R*) at the boundary surface.

_{eff}There are at least two methods reported in literature for calculating the effective refractive index, *R _{eff}*, which is based upon the physical principle that the diffuse flux normal to the interface is equal to the reflected part of the flux directed in the opposite direction [4, 8, 9, 15, 16].

The first method employs an empirical relationship between the relative refractive index with respect to air (*n _{medium}*/

*n*) and the internal reflection or the effective refractive index,

_{air}*R*of uniformly diffuse radiation. Eq. (3) represents the relationship developed from experimental measurements of

_{eff}*R*as a function of the relative refractive index,

_{eff}*n3*[15,16].

However, this empirical calibration is valid only for interfaces where at least one of the refractive indices is equal to one, or in other words, one of the surfaces is air.

In the second approach, *R _{eff}* is determined directly from Fresnel’s reflections,

*R*and

_{j}*R*and is in turn related to

_{Φ}*n*[4].

_{rel}where

and *R _{Fresnel}* is the Fresnel reflection coefficient for light incident upon the boundary at an angle of incidence, Θ[4].

The value of *R _{eff}* varies for different interfaces, consequently varying the

*γ*term in the partial current boundary condition. Table 1 provides the

*R*value, calculated using Eq. (4), (5) and (6), for various interfaces used in phantom tomography experiments. In order to simulate the tissue’s optical properties in phantom studies, diluted Intralipid solutions were used in this study. Intralipid, which is a suspension of lipids and nutrients in water, has similar optical properties as tissue in the NIR wavelength region [23]. Since dilute solutions of Intralipid (1% by volume) were used, its refractive index was approximated to the refractive index of water (

_{eff}*n*= 1.33).

Herein, we study the effect of incorporating different *R _{eff}* values for the boundary conditions of the coupled diffusion equations, and highlight their influence on the accuracy of the predicted excitation and emission fluence measured at the interface.

## 3. Materials and methods

Phantom measurements for optical tomography were conducted using 4×8×8 cm^{3} black acrylic box filled with 1% Intralipid solution (Pharmacia and Upjohn Co., Clayton NC) containing, 0.01μM indocyanine green (ICG) (Aldrich Chemical Co., Milwaukee WI). The black acrylic container has minimum transmission and refraction at the boundaries owing to its color; hence the reflections through the container are neglected during the estimation of the index-mismatch parameter. A 1×1×1 cm^{3} closed cuvette filled with 1% Intralipid containing 1μM ICG was positioned within the phantom (<1,4,3>cm from origin) in order to represent the target contrasted by 100 fold greater concentration of fluorescent dye. Figure 1a is a photograph of the actual acrylic phantom coupled to fibers for delivery and collection of light and Figure 1b represents a grid of the source-detector locations on the 3D phantom surfaces. The details of the optical properties of the homogeneous background and the target are given in Table 2.

#### 3.1 FDPM measurements

Millimeter diameter, multimode optical fibers (numerical aperture=0.39) were used as point sources and detectors, with 25 optical fibers spaced 1 cm apart from each other and 2 cm apart from the surface boundaries on either side of the 8×8 cm^{2} surface. Single-pixel (single point source and point detector) FDPM measurements were performed using 50 mW, 785 nm laser diode (Thorlabs Inc, Newton NJ) biased at 70 mA and modulated at a frequency of 100 MHz. The excitation and emission reflectances/transmittances (i.e. source and detectors are on the same surface or different surfaces, respectively) were collected at the phantom surface using appropriate interference filters and detected using modulated photomultiplier tubes (R928 model, Hamamatsu, Japan). The details of the instrumentation used herein are explained elsewhere [24]. All measurements were referenced to a measurement conducted at a specified “reference” position in order to account for the unknown source strength on the phantom. In addition, the phase shift and amplitude attenuation owing to instrument function was measured by the method of multiplexing, and then subtracted from the final measurements in order to obtain a final signal representative of light propagation through the phantom. [24,25]. There were 50 fiber positions, of which four were chosen as source fibers The ‘x’ symbols in Figure 1b represent the source locations on one side of the two 8×8 cm^{2} surfaces, each surface containing two sources at similar locations. For each of the four sources, one of the adjacent detector positions that is closest to the source (i.e. 1 cm away from the source in our case) was chosen as the “reference” position. Each measurement at all positions on the phantom was then reported relative to that made at the “reference” position. Hence, all measurements were referenced at a given wavelength for a given source location. In other words, referenced emission measurements consisted of relative measurements between a point and a “reference” location at the emission wavelength; while referenced excitation measurements consisted of relative measurements between a point and a “reference” location at the excitation wavelength.

#### 3.2 Data Analysis

The measurements are comprised of the phase-shift (θ) and amplitude attenuation (*I _{ac}*) of the diffuse excitation and emission reflectance at the surface of the phantom. The reflectance includes contributions from the fluence as well as flux at the interface as shown in Eq. (7). However, it is to be noted that the fluence or the flux at the surface of the phantom cannot be measured individually, but only as a combined reflectance term. The diffuse reflectance,

*R*(

*r*,

*ω*), is a function of the distance,

*r*, between the detected area and the incident beam. It also depends on the numerical aperture of the detector and contributions arising from the refractive-index mismatch at the interface. For FDPM reflectance measurements,

*R*(

*r*,

*ω*) is given by the following empirical equation [26].

where *A*, *B* are constants that vary depending on the numerical aperture of the detector and the *R _{eff}* value at the interface. The terms

*Φ*(

*r*,

*ω*) and

*∂Φ*(

*r*,

*ω*)/

*∂n*are the fluence and the flux at the surface boundaries respectively. Similar expressions can be written for time-dependent as well as CW reflectance measurements. For the case of partial current boundary condition (see Eq. (1)), one observes that the fluence is directly proportional to the flux at the boundaries and the proportionality constant is a function of the diffusion coefficient (

*D*) and index-mismatch parameter (

_{x,m}*γ*). Upon rearranging Eq. (1) and substituting for the flux term in Eq. (7),, one can see that the reflectance

*R*(

*r*,

*ω*), is directly proportional to the fluence alone, with a new proportionality constant. Upon referencing the measured diffuse reflectance,

*R*(

*r*,

*ω*) to that measured at a reference point

*R*(

*r*,

_{ref}*ω*), the new proportionality constant cancels out leaving behind the relationship:

where the suffix *ref* refers to the location of the reference measurement.

The constants *A* and *B* are a function of the numerical aperture of the collecting and illuminating fibers and the refractive index of the phantom’s interior and surfaces. In the case of absolute measurements, constants *A* and *B* need to be accounted for. However, upon referencing, the influence of constants *A* and *B* are cancelled as shown in Eq. (8). Thus, by conducting referenced FDPM measurements, we can directly compare referenced measurements to predictions without the need to empirically characterize the coefficients *A*, or *B* in Eq. (7). We also observe that upon using referenced measurement and the partial current boundary condition (where flux is proportional to the fluence), the sensitivity of these measurements on the flux and fluence terms is eliminated.

In terms of phase shift and amplitude attenuation, fluence is given by Eq. (9)

and Eq. (8) can be expressed in terms of *I _{ac}* and

*θ*as shown in Eq. (10).

where, *I _{ac}*/

*I*is the ac ratio and

_{ac,ref}*Δθ*is the relative phase shift (

*Δθ*=

*θ*-

*θ*) at each location

_{ref}*r*, and frequency

*ω*.

#### 3.3 Model prediction

For predictions of the referenced measured reflectance given by equation (8) and (9), the finite element method was employed. A finite element mesh was generated using Gambit software (Fluent. Inc, New Hampshire) for the geometry in Figure 1, and was discretized into 56423 tetrahedral elements and 11649 nodes. The Galerkin approximation of the finite-element method was employed to solve the forward problem of the coupled diffusion equations [27,28]. The details of this approach are explained in detail elsewhere [27]. The partial current boundary condition given in Eqs. (1) through (6) was employed with varying *R _{eff}* values to determine the index-mismatch parameter value (

*γ*) for the different cases as discussed in the next section. The aim of this study was to determine which value of the index-mismatch parameter simulates the actual measurements and how to estimate it from basic principles.

The excitation and emission fluence data obtained at each node of the phantom was referenced against the measurement at the reference position at the same wavelength. Both the experimental and simulated referenced fluence is represented in terms of ac ratio, *I _{ac}*/

*I*and relative phase shift,

_{ac,ref}*Δθ*for comparison studies.

## 4. Results and discussion

Solutions of the coupled diffusion equations with three different values of *γ* in the partial current boundary conditions were obtained to predict the referenced measurements of ac ratio, *I _{ac}*/

*I*and relative phase shift,

_{ac,ref}*Δθ*at the excitation and emission wavelengths. The three cases were: Case 1 in which the value of

*R*for an air-phantom interface was chosen for all six interfaces of the phantom (

_{eff}*R*= 0.4311); Case 2 in which the value of

_{eff}*R*for the phantom-acrylic interface was chosen for all six interfaces (

_{eff}*R*= 0.0222); and Case 3 in which the values of

_{eff}*R*corresponded to the actual phantom interfaces. In the latter case, five of the interfaces were phantom-acrylic and one was an air-phantom interface.

_{eff}#### 4.1 Model prediction for referenced measurements

Figures 2 and 3 represent the ratio of ac intensity, *I _{ac}*/

*I*and the relative phase-shift,

_{ac,ref}*Δθ*at the excitation and emission wavelengths, respectively at each of the 50 detectors in response to one source, located at <0,5,5> cm in the phantom The symbols (black hollow circles) denote measurements while the colored lines (connecting the symbols in Figure 2) represent predictions made using one of the three cases described above. Table 3 lists the average of the relative errors in ac ratio and average of the absolute errors in relative phase shift (in degrees) between the experimental and simulated data. From Figures 2 and 3, and the errors listed in Table 3, we observe that the model mismatch at both the excitation and emission wavelength is smallest in case 3, where the

*R*was calculated for each individual surface appropriately. Case 1 does not display a good model match in both the excitation and emission data in comparison to cases 2 or 3. The reason may be due to the approximation of air-tissue interface made for all the phantom interfaces, which do not correspond to the actual physical system. However, one may observe that the mean relative error in excitation ac data for case 1 is smaller compared to that obtained in cases 2 and 3 for the specific data set presented. Upon evaluating data obtained at the other three source locations, we found that this is not a consistent observation and underscores the lesser sensitivity of excitation measurements over emission measurements.

_{eff}At the emission wavelength, the degree of model mismatch is clearly evident in the ac ratio data, *I _{ac}*/

*I*as shown in Figure 3(a). The reason for the large impact may be owing to the larger attenuation of emission light by as much as 3 or 4 orders of magnitude over the measured excitation signal. Hence, even small inaccuracies of model predictions due to erroneous boundary condition parameters may constitute a greater proportion of the model-mismatch errors at emission wavelengths. Consequently, the errors arising from the model mismatch could prevent convergence of the inverse imaging problem to an accurate solution.

_{ac,ref}In case 2, which is a closer approximation to the actual physical system, the errors in *I _{ac}*/

*I*and

_{ac,ref}*Δθ*are comparable to the errors given by case 3 (representing the actual system) at excitation wavelength. Whereas, at emission wavelength, errors in ac ratio and relative phase shift are significantly higher for case 2 compared to the errors for case 3, which is also accounted by the weak emission signal in comparison to the excitation signal. All these comparisons were performed from data obtained for a particular source at <0,5,5> cm, located on one of the 8×8 cm

^{2}surface boundary. In brief, the choice of suitable value for the index-mismatch parameter would greatly improve the model mismatch between the measurements and the forward model, especially in the emission measurements of the fluorescence-enhanced optical tomography.

As a general trend in Figures (2) and (3), we also observe that there was greater mismatch for all of the three cases in which the detectors were located in transillumation geometry (i.e., when the detectors were located on the x=4 cm plane and the source was located on the x=0 cm plane), or alternatively detectors 1 to 25 in the x-axis of Figures (2) and (3). Thus the errors from the measurements on the transillumination side are enhanced due to greater source-detector distances in the system as well as the heterogeneity location. In this case, the heterogeneity was closer to the source and detectors on the x=0 cm plane relative to the detectors on the x=4 cm plane; thus weakening the measured signal at the detectors away from the source or the heterogeneity.

## 5. Conclusion

In summary, the choice of the appropriate index-mismatch parameter values at the boundary interface is vital to reduce the model mismatch of the forward model, especially in fluorescence-enhanced optical tomography. Determination of this index-mismatch parameter from basic principles would be more appropriate to minimize the error propagation in the coupled-diffusion equation, in turn improving the image reconstructions to a considerable extent. The reconstructions performed using the above stated 3D phantom and case 3 for the index-mismatch parameter are presented elsewhere [29]. The model mismatch in terms of ac ratio and relative phase shift, between measurements and predictions was prominent at emission wavelength than excitation wavelength, possibly due to the stronger excitation signal over the emission signal (3-4 orders of magnitude), thus masking the effect of the index-mismatch parameter.

## 6. Acknowledgements

This research was supported in parts by NIH (R01 CA88082) and (R01 CA67176).

## References and links

**1. **K. M. Case and P. F. Zweifel, *Linear Transport Theory* (Addison-Wesley, Massachusetts, 1967).

**2. **J. J. Duderstadt and L. J. Hamilton, *Nuclear Reactor Analysis* (Wiley, New York, 1976).

**3. **A. Ishimaru, “Diffusion of light in turbid media,” Appl. Opt. **28**, 2210–2215 (1989). [CrossRef] [PubMed]

**4. **R. C. Haskell, L. O. Scassand, T-T. Tsay, T-C. Feng, M. S. Mc Adams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A. **11**, 2727–2741 (1994). [CrossRef]

**5. **M. S. Patterson, B. Chance, and B. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. **28**, 2331–2336 (1989). [CrossRef] [PubMed]

**6. **A. H. Hielscher, S. L. Jacques, L. Wang, and F. K. Tittel, “The influence of boundary conditions on the accuracy of diffusion theory in time-resolved reflectance spectroscopy of biological tissues,” Phys. Med. Biol. **40**, 1957–1975 (1995). [CrossRef] [PubMed]

**7. **T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties *in vivo*,” Med. Phys. **9**, 879–888 (1992). [CrossRef]

**8. **M. Keijzer, W. M. Star, and P. R. Storchi, “Optical diffusion in layered media,” Appl. Opt. **27**, 1820–1824 (1988). [CrossRef] [PubMed]

**9. **R. Aronson, “Boundary conditions for diffusion of light,” J. Opt. Soc. Am. A. **12**, 2532–2539 (1995). [CrossRef]

**10. **R. Aronson, “Extrapolation distance for diffusion of light,” in *Photon Migration and Imaging in Random Media and Tissues*, B. Chance. and R. Alfano, Proc. Soc. Photo-Opt. Instrum. Eng. **1888**, 297–305 (1993).

**11. **M. G. Nichols, E. L. Hull, and T. H. Foster, “Design and testing of a white-light, steady-state diffuse reflectance spectrometer for determination of optical properties of highly scattering systems,” Appl. Opt. **36**, 93–104 (1997). [CrossRef] [PubMed]

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

**13. **M. A. Bartlett and H. Jiang, “Effect of refractive index on the measurement of optical properties in turbid media,” Appl. Opt. **40**, 1735–1741 (2001). [CrossRef]

**14. **B. W. Pogue, S. Geimer, T. O. McBride, S. Jiang, U. L. Osterberg, and K. D. Paulsen, “Three-dimensional simulation of near-infrared diffusion in tissue: boundary condition and geometry for finite-element image reconstruction,” Appl. Opt. **40**, 588–600 (2001). [CrossRef]

**15. **R. A. J. Groenhuis, H. A. Ferwerda, and J. J. Ten Bosch, “Scattering and absorption of turbid material determined from reflection measurements. 1. Theory,” Appl. Opt. **22**, 2456–2462 (1983). [CrossRef] [PubMed]

**16. **W. G. Egan and T. W. Hilgeman, *Optical properties of inhomogeneous materials* (Academic, New York, 1979).

**17. **D. J. Durian, “Influence of boundary reflection and refraction on diffusive photon transport,” Phys. Rev. E **50**, 857–866 (1994). [CrossRef]

**18. **E. M. Sevick-Muraca and D. Y. Paithankar, “Fluorescence imaging system and measurement,” U. S. Patent No. 5,865,754 (2 February 1999).

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

**20. **E. M. Sevick and C. L. Burch, “Origin of phosphorescence signals reemitted from tissues.” Opt. Lett. **19**, 1928–1930 (1994). [CrossRef]

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

**22. **C. L. Hutchinson, J. R. Lakowicz, and E. M. Sevick-Muraca, “Fluorescence life-time based sensing in tissues: A computational study,” Biophys. J. **68**, 1574–1582 (1995). [CrossRef] [PubMed]

**23. **H. J. van Staveren, C. J. M. Moes, J. van Marle, S. A. Prahl, and M. J. C. van Gemert, “Light scattering in Intralipid-10% in the wavelength range of 400-1100 nm,” Appl. Opt. **30**, 4507–4514 (1991). [CrossRef] [PubMed]

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

**25. **D. J. Hawrysz, *Bayesian approach to the inverse problem in contrast-enhanced, three-dimensional, biomedical optical imaging using frequency domain photon migration*, PhD Thesis, Purdue University, May 2001.

**26. **A. Kienle and M. S. Patterson, “Improved solution of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,” J. Opt. Soc. Am. A. **14**, 246–254 (1997). [CrossRef]

**27. **R. Roy and E. M. Sevick-Muraca, “Truncated Newton’s optimization scheme for absorption and fluorescence optical tomography: Part I Theory and Formulation” Opt. Express **4**, 353–371 (1999). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-4-10-353. [CrossRef] [PubMed]

**28. **O. C. Zeinkiewicz and R. L. Taylor. *The finite element methods in engineering science* (McGraw-Hill, New York, 1989).

**29. **M. J. Eppstein, D. J. Hawrysz, A. Godavarty, and E. M. Sevick-Muraca, “Three-dimensional, near-infrared fluorescence tomography with Bayesian methodologies for image reconstruction from sparse and noisy data sets,” Proc. of Natl. Acad. Science2002 (accepted). [CrossRef]