## Abstract

We have developed a theoretical model for photon migration through scattering media in the presence of an absorbing inhomogeneity. A closed-form solution for the average diffuse intensity has been obtained through an iterative approximation scheme of the steady-state diffusion equation. The model describes absorbing defects in a wide range of values. Comparisons with the results of Monte Carlo simulations show that the error of the model is lower than 3% for size inclusion lower than 4 mm and absorption contrast up to the threshold value of the “black defect.” The proposed model provides a tractable mathematical basis for diffuse optical and photoacoustic tomographic reconstruction techniques.

© 2014 Optical Society of America

A physical understanding of photon migration through highly scattering media is required to study biological tissues by means of optical methods for medical diagnostic applications. Light propagation in tissue is well described by the diffusion equation (DE), and this has led to the development of several analytical solutions of this equation that provide a formal mathematical basis to separate the effects of scattering from absorption in tissue. Some models can even describe the effects of local changes in tissue optical properties. In optical tomography, these changes can be related to composition, vascularization, blood oxygen saturation, and structural properties of normal tissue, benign lesions, and tumors as well [1]. Furthermore, in quantitative photoacoustic tomography, changes in scattering and absorption properties are combined with the ultrasound propagation to estimate the concentration of chromophores, such as hemoglobin and melanin, inside tissues [2,3]. The DE is quite versatile and it is also adopted to detect concentrations and lifetimes of fluorescent probes that are used to improve tumor contrast [4,5]. The underlying mathematical model of the measured fluorescence probes signal is formally similar to that employed for describing the diffuse intensity in turbid media in the presence of absorptive inhomogeneities. Indeed, the detection of absorbing inclusions is an issue widely investigated in diffuse optical imaging.

Many approaches have attempted to provide forward solvers describing the effects of absorbing inclusions using a perturbation scheme applied to the DE [6–12] or a random walk scheme [13,14]. In particular, the Born approximation has been widely used to describe absorbing defects for imaging applications to biological tissues, although the accuracy of the method is limited to weak perturbations of a few percent. To overcome this limitation, several iterative schemes have been developed to improve the accuracy of the calculated diffuse intensity. These approaches have been based on the calculation of high-order terms in the expansion series of the intensity [10–12]. However, in many cases, the use of these approximations has been severely limited by the long computation time [10]. On the other hand, computationally faster methods using a heuristic perturbative approach based on higher-order perturbation theory and the Padé approximants were developed [11,12]. These forward solvers show a negligible computation time for the calculation of the average intensity, although the complexity of the theory behind such solvers may discourage their use. Moreover, all these theories show deficiencies when used with high or totally absorbing (“black”) defects where a saturation effect of contrast is achieved. Situations like this can be usually found in biological tissues such as in the case of vessels, which are structures with a strong localized absorption and produce a marked relative contrast in optical and photoacoustic tomography. For this evident argument, it will be desired to provide improved forward solvers that are able to overcome the limitations of the existing ones.

In this Letter, we derive a closed-form expression for the average intensity of the steady-state DE for a turbid medium in the presence of an absorbing defect through an iterative approximation scheme. The main novelty of the proposed solution compared to conventional approaches for the intensity calculations is twofold. (i) A closed-form expression of the diffuse intensity which is amenable both to weak and strong absorbing inhomogeneity up to the case of a “black defect.” (ii) The simplicity of the analytical expression of the intensity. This second aspect makes it easy to understand the role of the physical parameters of the proposed forward model and also simplifies its use in the forward and inverse problem. Furthermore, the method calculates the intensity in a computation time much smaller than previous approaches. Typically it takes about 0.1 ms using a homebuilt routine package written in MATLAB [15] and an Intel(R) Core(TM)i7-2600 CPU@3.40 GHz processor with 8 GB memory. This computation time is at least one order of magnitude lower than the Fortran routine for the fourth-order method based on the Padé approximants described in Ref. [11].

The steady-state DE for the average diffuse intensity $\phi (\mathbf{r},{\mathbf{r}}_{0})$ in a homogeneous scattering medium of absorption coefficient ${\mu}_{a}$ and in the presence of an absorptive defect with contrast $\mathrm{\Delta}{\mu}_{a}(\mathbf{r})$ is

The second-order term is

According to Eq. (8), the sum of the geometrical series at the righthand side of the Eq. (2) can be easily performed, and we obtain

Equation (9) gives a simple closed-form expression for the correction to the unperturbed intensity ${\phi}_{0}(\mathbf{r},{\mathbf{r}}_{0})$ due to an absorptive inclusion, and it is written as the sum of two terms. The first term is proportional to the product $\mathrm{\Delta}{\mu}_{a}V$ of the absorption contrast and its volume, and it is the leading contribution of the first-order perturbation approximation, while the second term is the contribution of high-order terms in the perturbation expansion. Equation (9) also shows that, in case of a “black defect,” i.e., when $\mathrm{\Delta}{\mu}_{a}\to \infty $, the limit value of the correction to the perturbed intensity reduces to

The calculation of the perturbed intensity $\phi (\mathbf{r},{\mathbf{r}}_{\mathbf{0}};\mathrm{\Delta}{\mu}_{a})$ through the slab requires the calculation of the average value $\overline{G}(\mathbf{r})$ of the Green’s function. According to the definition given by Eq. (6), the average value is obtained by integrating the Green’s function over the volume $V$ of the spherical region of the defect centered at ${\mathbf{r}}_{\mathbf{c}}=({x}_{c},{y}_{c},{z}_{c})$. Although it is possible to formulate a general expression for $\overline{G}(\mathbf{r})$ which takes into account the contributions of both the positive and negative sources, it turns out that the average value at the center of the defect is essentially determined by the positive zero-order term, corresponding to the solution for a homogeneous infinite medium geometry, because the higher-order terms in Eq. (10) give a negligible contribution for $R\ll |{\mathbf{r}}_{\mathbf{c}}-{\mathbf{r}}_{\mathbf{m}}|$. Accordingly, we have the following simple expression for the average value at the center of the defect

where $\alpha =\sqrt{{\mu}_{a}/D}$.The accuracy of the model for the slab geometry is investigated by means of comparisons with the results of Monte Carlo (MC) simulations [17] for the relative contrast $|\mathrm{\Delta}I/{I}_{0}|=|(I-{I}_{0})/{I}_{0}|$, where ${I}_{0}$ is the unperturbed transmittance. In this regard, Fig. 2 shows the dependence of $|\mathrm{\Delta}I/{I}_{0}|$ for the transmittance through a slab versus the absorption contrast $\mathrm{\Delta}{\mu}_{a}$ of the defect for three values of the volume $V$: 100, 300, and $500\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$. The perturbation of the flux $\mathrm{\Delta}I$ in the slab geometry is calculated as $\mathrm{\Delta}I({\mathbf{r}}_{\mathbf{m}},{\mathbf{r}}_{\mathbf{0}})=(2\pi /A)\mathrm{\Delta}\phi ({\mathbf{r}}_{\mathbf{m}},{\mathbf{r}}_{\mathbf{0}})$, where the perturbed fluence $\mathrm{\Delta}\phi $ is obtained with Eq. (9). The inclusion is assumed to be located midway between the source and the detector at ${z}_{c}=d/2$ in coaxial configuration ($\rho =0$) in panel (a), and at radial distance $\rho =5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ in panels (b) and (c). The thickness of the slab $d$ and the optical properties are reported in the inset of the figure. The discrepancies between the presented solver and MC simulations are within 11% for the all situations considered in the figure. Remarkably, the theory also accounts for the saturation of the relative intensity for increasing values of the absorption contrast when the inclusion behaves like a black defect, and this is observed for a wide range of optical and geometrical properties that cover tissue optics applications. Figure 3 shows the relative contrast $|\mathrm{\Delta}I/{I}_{0}|$ in the case of a black defect versus the inclusion volume: the red dotted curve is the model and the solid black curve is the MC simulation. The results are obtained for the same slab and the same optical parameters used for case (a) of Fig. 2. As can be noticed, the theoretical predictions are in good agreement with MC simulations for volumes of the black defect ranging from 50 to $500\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ ($R\le 5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$). In particular, excellent agreement between the theory and the MC results curve is obtained for volumes lower than $300\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$. Indeed, the discrepancies (solid blue curve) are lower than 3% for volumes increasing from 50 to $300\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$. This is an expected result because the theoretical model is derived within the small volume approximation. Conversely, the theory overestimates the MC simulations when $V\ge 300{\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}}^{3}$, with discrepancies that tend to increase up to a value within 11% in the considered range. In this case, the accuracy is limited essentially by the larger values of the volume of the inclusion.

In conclusion, we have developed a simple theoretical approach for photon migration through scattering media in the presence of an absorptive inhomogeneity. A closed-form expression for the average diffuse intensity is derived in the framework of an iterative approximation scheme of the steady-state DE. This method addresses the case of a single inclusion but it should be remarked that this limitation, which is also present in current analytical and perturbation schemes, can be released somewhat in the presence of multiple defects whose relative distances are large enough to make negligible their interactions. The slab geometry was used to test the closed-form solution given by Eq. (9) in transmission configuration. Comparable accuracy was also obtained in the case of reflectance geometry. It is worth pointing out that this solution can be applied to both homogeneous or inhomogeneous geometries. Indeed, there exists a variety of available analytical solutions of the DE such as, for instance, the parallelepiped, the sphere, the cylinder, and several layered geometries [17]. Furthermore, Eq. (9) is applicable also in case no explicit or analytical expressions of the Green’s function are available.

## References

**1. **T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, Rep. Prog. Phys. **73**, 076701 (2010). [CrossRef]

**2. **B. Cox, S. Arridge, K. Kostli, and P. Beard, Appl. Opt. **45**, 1866 (2006).

**3. **A. Q. Bauer, R. E. Nothdurft, T. N. Erpelding, L. V. Wang, and J. P. Culver, J. Biomed. Opt. **16**, 096016 (2011). [CrossRef]

**4. **T. Correia, N. Ducros, C. D’Andrea, M. Schweiger, and S. Arridge, Opt. Lett. **38**, 1903 (2013). [CrossRef]

**5. **V. Ntziachristos and R. Weissleder, Opt. Lett. **26**, 893 (2001). [CrossRef]

**6. **M. R. Ostermeyer and S. L. Jacques, J. Opt. Soc. Am. A **14**, 255 (1997). [CrossRef]

**7. **S. Carraresi, T. S. M. Shatir, F. Martelli, and G. Zaccanti, Appl. Opt. **40**, 4622 (2001). [CrossRef]

**8. **S. De Nicola, R. Esposito, and M. Lepore, Phys. Rev. E **68**, 21901 (2003). [CrossRef]

**9. **S. Feng, F. A. Zeng, and B. Chance, Appl. Opt. **34**, 3826 (1995). [CrossRef]

**10. **D. Grosenick, A. Kummrow, R. Macdonald, P. M. Schlag, and H. Rinneberg, Phys. Rev. E **76**, 061908 (2007). [CrossRef]

**11. **A. Sassaroli, F. Martelli, and S. Fantini, Appl. Opt. **48**, D62 (2009). [CrossRef]

**12. **A. Sassaroli, F. Martelli, and S. Fantini, J. Opt. Soc. Am. A **27**, 1723 (2010). [CrossRef]

**13. **A. G. Gandjbakhche, R. F. Bonner, R. Nossal, and G. H. Weiss, Appl. Opt. **35**, 1767 (1996). [CrossRef]

**14. **V. Chernomordik, D. Hattery, A. H. Gandjbakhche, A. Pifferi, P. Taroni, A. Torricelli, G. Valentini, and R. Cubeddu, Opt. Lett. **25**, 951 (2000). [CrossRef]

**15. **G. Started, *The Language of Technical Computing* (Prentice Hall, 1997), Vol. L, p. 186.

**16. **G. Arfken, *Mathematical Methods for Physicists*, 2nd ed. (Elsevier, 1971), Vol. 39.

**17. **F. Martelli, S. Del Bianco, A. Ismaelli, and G. Zaccanti, *Light Propagation through Biological Tissue and Other Diffusive Media: Theory, Solutions, and Software* (SPIE, 2010).

**18. **R. Esposito, S. De Nicola, M. Lepore, and P. L. Indovina, J. Opt. Soc. Am. A **23**, 1937 (2006). [CrossRef]