Closed-form solution of the steady-state photon diffusion equation in the presence of absorbing inclusions

Rosario Esposito; Fabrizio Martelli; Sergio De Nicola

doi:10.1364/OL.39.000826

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 $φ (r, r_{0})$ in a homogeneous scattering medium of absorption coefficient $μ_{a}$ and in the presence of an absorptive defect with contrast $Δ μ_{a} (r)$ is

D \nabla^{2} φ (r, r_{0}) - (μ_{a} + Δ μ_{a} (r)) φ (r, r_{0}) = - \frac{δ (r - r_{0})}{4 π},

where

D = 1 / (3 μ_{s}^{'})

is the diffusion coefficient, with

μ_{s}^{'}

the reduced scattering coefficient of the host medium. We assume an isotropic point-like source located at

r_{0}

, and we consider

Δ μ_{a} (r) = Δ μ_{a}

inside

V

and

Δ μ_{a} (r) = 0

outside

V

. Equation (1) can be solved iteratively, and this leads to a scheme of successive approximations for the average intensity [16]

φ (r, r_{0}) = φ_{0} (r, r_{0}) + \sum_{n = 1}^{\infty} φ_{n} (r, r_{0}) .

The zero-order term

φ_{0} (r, r_{0})

is the average intensity in the background medium in the absence of the defect

(Δ μ_{a} (r) = 0)

, and the

n

th order term

φ_{n} (r, r_{0})

can be written recursively according to the equation

φ_{n} (r, r_{0}) = - Δ μ_{a} \int_{V} d r' G (r - r') φ_{n - 1} (r', r_{0}),

where

G (r - r')

is the Green’s function of Eq. (1) with

(Δ μ_{a} (r) = 0)

. Note that

G (r - r_{0})

is the fluence at point

r

due to a source located at

r_{0}

, and is related to the zero-order term

φ_{0} (r, r_{0})

by the simple relation

G (r - r_{0}) = 4 π φ_{0} (r, r_{0})

. In the following we will show that, under suitable approximations, we can derive a closed-form expression of the diffuse intensity series given by Eq. (2). In order to do this, it is necessary to have a general way of expressing the

n

th order term in the expansion series for the intensity. The first-order term is

φ_{1} (r, r_{0}) = - \frac{Δ μ_{a}}{4 π} \int_{V} d r' G (r - r') G (r' - r_{0}),

where the integration is performed over the volume

V

of the absorptive inclusion. If the linear size of the volume

V

is much smaller than the distance

| r_{c} - r_{0} |

between the center

r_{c}

of the defect and the source point

r_{0}

, then we can approximate

G (r' - r_{0})

with

G (r_{c} - r_{0})

in the integral of Eq. (4), and we find

φ_{1} (r, r_{0}) = - \frac{Δ μ_{a} V}{4 π} G (r_{c} - r_{0}) \bar{G} (r),

where we have defined the average value of the Green’s function over the inclusion volume

\bar{G} (r) = \frac{1}{V} \int_{V} d r' G (r - r') .

The second-order term is

φ_{2} (r, r_{0}) = \frac{Δ μ_{a}^{2} V^{2}}{4 π} \bar{G} (r_{c}) G (r_{c} - r_{0}) \bar{G} (r),

where the product of two average values of the Green’s function appears. Proceeding in this way, we find that the general expression of the

n

th order term is

φ_{n} (r, r_{0}) = \frac{{(- Δ μ_{a} V \bar{G} (r_{c}))}^{n}}{4 π \bar{G} (r_{c})} G (r_{c} - r_{0}) \bar{G} (r) .

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

Δ φ (r, r_{0}; Δ μ_{a}) = (- Δ μ_{a} V + \frac{Δ μ_{a}^{2} V^{2} \bar{G} (r_{c})}{1 + Δ μ_{a} V \bar{G} (r_{c})}) \frac{G (r_{c} - r_{0}) \bar{G} (r)}{4 π},

where we have defined the perturbation of the average intensity as

Δ φ (r, r_{0}; Δ μ_{a}) = \sum_{n = 1}^{\infty} φ_{n} (r, r_{0})

.

Equation (9) gives a simple closed-form expression for the correction to the unperturbed intensity $φ_{0} (r, 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 $Δ μ_{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 $Δ μ_{a} \to \infty$ , the limit value of the correction to the perturbed intensity reduces to

Δ φ (r, r_{0}; Δ μ_{a} ≫ μ_{a}) = - \frac{G (r_{c} - r_{0}) \bar{G} (r)}{4 π \bar{G} (r_{c})} .

If the linear size of the volume

V

of the inclusion is smaller than the distance between the center

r_{c}

of the defect and the measurement point

r

, as is often the case, we can further approximate the average value

\bar{G} (r)

in Eq. (9) with

G (r - r_{c})

. In the following, we specialize the developed general formalism to the case of an absorptive spherical inclusion of radius

R

embedded in an infinitely extended slab of thickness

d

, as illustrated in Fig. 1. The source, for the case of a pencil light beam impinging the medium at

(x_{0}, y_{0}, 0)

, can be placed at

r_{0} = (x_{0}, y_{0}, z_{0} = 1 / (μ_{s}^{'} + μ_{a}))

[17], and the inclusion is located at

r_{c} = (x_{c}, y_{c}, z_{c})

. The transmitted flux is calculated at

r_{m} = (x_{m}, y_{m}, z_{m} = d)

, while in a reflectance configuration, we have

r_{m} = (x_{m}, y_{m}, z_{m} = 0)

. The Green’s function

G (r - r_{0})

of the homogeneous absorbing slab of Fig. 1 can be expressed as an infinite sum of dipoles located along the

z

axis that represents pairs of negative and positive point-like sources, namely

G (r - r_{0}) = \frac{1}{4 π D} \sum_{n = - \infty}^{\infty} \frac{\exp (- \sqrt{\frac{μ_{a}}{D}} ρ_{+})}{ρ_{+}} - \frac{\exp (- \sqrt{\frac{μ_{a}}{D}} ρ_{-})}{ρ_{-}} .

We have defined the quantities

ρ_{\pm}^{2} = ρ^{2} + {(z - z_{\pm})}^{2}

,

z_{+} = 2 n d_{e} + z_{0}

, and

z_{-} = 2 n d_{e} - 2 z_{e} - z_{0}

where the extrapolated length is

d_{e} = d + 2 z_{e}

. The square of the radial distance from the

z

axis is

ρ^{2} = {(x - x_{0})}^{2} + {(y - y_{0})}^{2}

. The average intensity given by Eq. (10) vanishes at extrapolated flat surfaces outside the slab at a distance

z_{e} = 2 A D

from the boundaries, where the coefficient

A

accounts for the refractive index mismatch between the slab and the external [17]. Taking into account the extrapolated boundary conditions, the outward photon flux measured at the position

r_{m}

on the surface of the slab is given by

I (r_{m}, r_{0}) = (2 π / A) φ (r_{m}, r_{0})

[17,18].

Fig. 1. Schematic of the heterogeneous slab with the absorptive inclusion.

Download Full Size | PDF

The calculation of the perturbed intensity $φ (r, r_{0}; Δ μ_{a})$ through the slab requires the calculation of the average value $\bar{G} (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 $r_{c} = (x_{c}, y_{c}, z_{c})$ . Although it is possible to formulate a general expression for $\bar{G} (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 ≪ | r_{c} - r_{m} |$ . Accordingly, we have the following simple expression for the average value at the center of the defect

\bar{G} (r_{c}) = \frac{1}{V μ_{a}} [1 - e^{- R α} (1 + R α)],

where

α = \sqrt{μ_{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 $| Δ I / I_{0} | = | (I - I_{0}) / I_{0} |$ , where $I_{0}$ is the unperturbed transmittance. In this regard, Fig. 2 shows the dependence of $| Δ I / I_{0} |$ for the transmittance through a slab versus the absorption contrast $Δ μ_{a}$ of the defect for three values of the volume $V$ : 100, 300, and $500 {mm}^{3}$ . The perturbation of the flux $Δ I$ in the slab geometry is calculated as $Δ I (r_{m}, r_{0}) = (2 π / A) Δ φ (r_{m}, r_{0})$ , where the perturbed fluence $Δ φ$ 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 ( $ρ = 0$ ) in panel (a), and at radial distance $ρ = 5 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 $| Δ 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 {mm}^{3}$ ( $R \leq 5 mm$ ). In particular, excellent agreement between the theory and the MC results curve is obtained for volumes lower than $300 {mm}^{3}$ . Indeed, the discrepancies (solid blue curve) are lower than 3% for volumes increasing from 50 to $300 {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 \geq 300 {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.

Fig. 2. Relative contrast $| Δ I / I_{0} |$ for the transmittance through a slab plotted against the absorption contrast $Δ μ_{a}$ between a spherical absorptive inclusion and the background medium.

Download Full Size | PDF

Fig. 3. Relative contrast as a function of the inclusion volume $V$ of the black defect. The red dotted curve is the model prediction and the solid black curve is the MC simulation. The solid blue curve shows the relative discrepancies. The slab geometry and the optical parameters are the same as in Fig. 2(a).

Download Full Size | PDF

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]

Closed-form solution of the steady-state photon diffusion equation in the presence of absorbing inclusions

Abstract

References

Cited By

Figures (3)

Equations (12)

Optics Letters