On the geometry dependence of differential pathlength factor for near-infrared spectroscopy. I. Steady-state with homogeneous medium

Daqing Piao; Randall L. Barbour; Harry L. Graber; Daniel C. Lee M.D.

doi:10.1117/1.JBO.20.10.105005

14 October 2015 On the geometry dependence of differential pathlength factor for near-infrared spectroscopy. I. Steady-state with homogeneous medium

Daqing Piao, Randall L. Barbour, Harry L. Graber, Daniel C. Lee M.D.

Author Affiliations +

Funded by: National Institutes of Health

Journal of Biomedical Optics, Vol. 20, Issue 10, 105005 (October 2015). https://doi.org/10.1117/1.JBO.20.10.105005

Abstract

This work analytically examines some dependences of the differential pathlength factor (DPF) for steady-state photon diffusion in a homogeneous medium on the shape, dimension, and absorption and reduced scattering coefficients of the medium. The medium geometries considered include a semi-infinite geometry, an infinite-length cylinder evaluated along the azimuthal direction, and a sphere. Steady-state photon fluence rate in the cylinder and sphere geometries is represented by a form involving the physical source, its image with respect to the associated extrapolated half-plane, and a radius-dependent term, leading to simplified formula for estimating the DPFs. With the source-detector distance and medium optical properties held fixed across all three geometries, and equal radii for the cylinder and sphere, the DPF is the greatest in the semi-infinite and the smallest in the sphere geometry. When compared to the results from finite-element method, the DPFs analytically estimated for 10 to 25 mm source–detector separations on a sphere of 50 mm radius with μ_a=0.01 mm−¹and μ_s′=1.0 mm−¹ are on average less than 5% different. The approximation for sphere, generally valid for a diameter ≥20 times of the effective attenuation pathlength, may be useful for rapid estimation of DPFs in near-infrared spectroscopy of an infant head and for short source–detector separation.

1. Introduction

Near-infrared spectroscopy (NIRS)¹ is becoming an increasingly important modality for a number of investigational and clinical needs, including noninvasive functional study of neurophysiological connectivity,² bedside monitoring of cerebral hemodynamics,³ and evaluating responses to cancer treatment.⁴ Most NIRS systems for long duration monitoring operate under the continuous-wave (CW) or steady-state condition.⁵ Quantitation of the changes of chromophores, most importantly oxyhemoglobin and deoxyhemoglobin, within the highly diffusing tissue typically relies upon a priori estimates of the effective pathlength of the photons that propagate from a source to a detector. Delpy et al.⁶ first suggested utilizing the ratio of the effective pathlength of the photon trajectory over fixed separations of source–detector pairs as the differential pathlength factor (DPF) in the routine for quantitative NIRS. A DPF value of $5.3 \pm 0.3$ across a rat head was experimentally obtained by using time-of-flight techniques. Although the DPF was observed to be approximately constant for source–detector spacings greater than 2.5 cm, frequency-domain NIRS studies by Duncan et al.⁷ at source–detector spacings $> 4 cm$ yielded DPFs of 6.26 ( $\pm 14.1 %$ ) for an adult head, 5.51 ( $\pm 18 %$ ) for an adult leg, 4.99 ( $\pm 9 %$ ) for the head of a newborn infant, and 4.16 ( $\pm 18.8 %$ ) for an adult arm. These experimentally obtained values indicate that the DPF may depend on the dimensions and shape of the tissue structure under investigation, providing the same source–detector separations (fixed physical configuration) and tissue optical properties (likely subject condition).

In terms of the analytical form of DPF, Delpy et al.⁶ suggested a wavelength-, medium-, and geometry-dependent constant⁸ in addition to the term straightforwardly derived from the solution to photon diffusion in an unbounded homogeneous medium. In the case of a homogeneous semi-infinite medium, explicit analytical formation of the DPF was given by reports including Fantini et al.⁹ and Boas et al.¹⁰ A semi-infinite geometry is a reasonable representation of the measurement geometry for a set of optodes separated by a few centimeters when externally placed on large tissue areas. However, in cases involving measurement of the brain, different head sizes (infant to adult), different global shapes of the head (Caucasian to Asian), different regions of the head¹¹ (e.g., prefrontal or parietal), and even different ages¹² may result in different DPFs. Therefore, an analytical form of DPF that indicates the effect of the size and shape of the nonsemi-infinite geometry, even at the simplest case of a homogeneous medium that is the basis for analysis of more realistic geometries, may be informative for the development of better strategies to the quantitation of the chromophore changes when the wavelength dependency of DPF is also taken into account.¹³ In addition, the more accurate a priori estimation of the DPF for the subsets of optode configurations with shorter source–detector separations, the more reliable the determination of the shallow-layer contamination to the deeper-layer signal.¹⁴ Whereas estimates of DPFs from complex geometries are readily available with the use of discrete numerical methods [e.g., finite-element method (FEM)¹⁵^–¹⁷], access to analytical estimators of DPFs that can be quickly implemented with reasonable tolerance can hold considerable practical and clinical value, especially in situations where near immediate or real-time quantitates of DPFs are required (e.g., acute care environments) to allow timely diagnosis and intervention. Mainly lacking in support of this, an analytically rendered understanding of the DPF in geometries more complex than the semi-infinite case, specifically how the shape and dimension of a spherically curved geometry resembling an infant head affect the DPF, is relevant to NIRS applications for better clinical care.

Mathematical treatments of analytical solutions to the photon diffusion equation (DE) in a medium bounded externally by a circular cylindrical applicator or a spherical applicator have been investigated over a course of nearly two decades in many comprehensive studies.¹⁸^–²² A particularly complete set of solutions in time domain and frequency domain in a number of geometries were presented by Arridge et al.¹⁸ in 1992; the DPFs could be obtained from the explicit solutions derived for the mean time of flight in a variety of geometries using the zero boundary condition. In 2001, Sassaroli et al.²² developed the corresponding solutions of time-resolved DE using the more widely used extrapolated boundary condition (EBC)²¹^,²³^,²⁴ for curved geometries including infinite cylinder and sphere; the DPFs for infinite cylinder and sphere can also be derived or computed from the results that were in excellent agreement with Monte Carlo (MC) simulations. While these studies as represented by Refs. 18 and 22 have theorized the computational availability of accurate time-domain DPFs for curved geometries of interest, the computational complexity arising from applying high-precision arithmetic library on infinite-series-summations involving various special functions renders them out of reach for most members of the NIRS community. Solutions for the CW reflectance and subsequently the DPF that preferably involve only the regular functions and do not impose lengthy summations, while perhaps less accurate, may be much more convenient for most clinically oriented CW NIRS applications in rapidly estimating the DPF for timely physiological assessment. This present study aims to develop the first steps toward this unmet need.

In our recent approach to a new routine of solutions of DE for an infinite-length cylinder domain as demonstrated in a series of work applicable to CW,²⁵ frequency domain²⁶ and time domain,²⁷ the solutions for DE were presented using a pair of the modified Bessel functions of the first and the second kinds. The basic form of the DE solution for the unbounded homogeneous medium geometry shown in Refs. 25 26.–27 may be similar to the one used by Ref. 28; however, our method of utilizing the EBC in the case of the medium bounded within an infinite-length cylinder was new in terms of leading to a more intuitive interpretation of the solution. The solution obtained for the infinite-length cylinder domain is presented to consist of two parts: the first part is associated with the “real” isotropic source that is the common approach of treating the light normally incident to the medium from the source fiber or channel, and the second part accounts for the contribution of the “image” source as a result of the medium boundary that is represented by the “real” source term scaled by a factor determined by the radius of the cylinder. This solution is subsequently formatted to a form similar to that for the semi-infinite geometry except a shape-curvature-associated term that makes the solution for the cylinder domain approaching the one for semi-infinite geometry as the radius of the cylinder domain increases without bound. Such treatment allowed obtaining approximating solutions of CW DE in cylinder geometries of relatively large radius using functions more convenient to evaluate than special functions and without involving the infinite series of summations. Additionally, the analytical methodology explored in Refs. 25 26.–27 applies to two opposite medium geometries, one for the medium enclosed by the optode applicator as with the arm or leg, and the other for the medium enclosing the optode applicator as with the prostate via transrectal probing.

In continuation of our demonstrated approaches to the homogeneous infinite cylinder domain, we develop in this present work the solution of steady-state DE in a homogeneous spherical domain (limited to the external applicator case although the approach is extendable to the internal applicator case) to a form that is also composed of two parts: the first part is associated with the “real” isotropic source, and the second part is the contribution of the “image” source term that is represented by the “real” source term scaled by a factor depending upon the radius of the sphere. The format of the solution for the spherical domain, or the subsequent approximation, shall also be similar to that for the semi-infinite geometry except for a shape-curvature-associated term, which nonetheless must be different from its counterpart for the infinite cylinder, to make the solution for the sphere domain approaching the one for semi-infinite geometry, as the radius of the sphere increases without bound.

With the solutions of steady-state DE associated with the geometries of interest becoming available, the DPF can then be evaluated by calculating the partial derivative of the attenuation of the photon fluence with respect to the absorption and dividing that by the source–detector separation. The usefulness of the quantitative patterns of analytically approximated DPF associated with curved geometries can then be examined by using numerical simulation, which in this work is limited to FEM method.

2. Geometries Considered and Analytic Preparations for Steady-State Photon Diffusion in Infinite and Semi-Infinite Geometries

We consider a scattering domain with an absorption coefficient of $μ_{a}$ and a diffusion coefficient of $D = 1 / (3 μ_{s}^{'})$ where $μ_{s}^{'}$ is the reduced scattering coefficient. It must be noted that it is yet to reach a consensus regarding whether the diffusion coefficient in steady-state measurement shall include the absorption coefficient.²⁹^–³³ We argue that, as the DPF is unbiasedly defined by the use of mean time of flight through time-dependent DE, the definition of DPF in steady-state may agree with the one in the time-resolved condition only if the diffusion coefficient is not dependent upon the absorption coefficient. Therefore, this study uses the diffusion coefficient as $D = 1 / (3 μ_{s}^{'})$ .

This study intends to analytically compare the DPF in three geometries of the bounded medium that include a semi-infinite geometry, an infinite-length cylinder geometry, and a spherical geometry. Since the solutions for bounded geometries originate from the solution for the infinite medium geometry, it becomes imperative to establish a general analytical routine to the studied geometries as are shown in Fig. 1, and from which to evaluate the DPF of steady-state measurement. For obvious conveniences, the infinite and semi-infinite domains will be represented using the Cartesian coordinates, where we are concerned about photon pathlength from a source at ${\vec{χ}}^{'} = (x^{'}, y^{'}, z^{'})$ to a detector at $\vec{χ} = (x, y, z)$ . The infinite-length cylinder domain will be represented with the polar cylindrical coordinates, where we are concerned about photon pathlength from a source at ${\vec{χ}}^{'} = (ρ^{'}, ϕ^{'}, z^{'})$ to a detector at $\vec{χ} = (ρ, ϕ, z)$ . The sphere domain will be modeled in the spherical coordinates, where we are concerned about photon pathlength from a source at ${\vec{χ}}^{'} = (r^{'}, θ^{'}, ϕ^{'})$ to a detector at $\vec{χ} = (r, θ, ϕ)$ .

Fig. 1

Four cases of medium geometry. (a) Infinite geometry, where an isotropic source is at ${\vec{χ}}^{'} = (0,0, 0)$ and a detector is at $\vec{χ} = (0, d, 0)$ . (b) Semi-infinite geometry wherein the source and the detector locate at ${\vec{χ}}^{'} = (0,0, 0)$ and $\vec{χ} = (0, d, 0)$ , respectively. The equivalent isotropic source locates inwardly with respect to the interface at $(0,0, - R_{a})$ , and the image of it with respect to the boundary locates at ( $0, 0, R_{a} + 2 R_{b}$ ). (c) Infinite-length cylindrical geometry, where the source and the detector locate at $(R_{0}, ϕ^{'}, z^{'})$ and $(R_{0}, ϕ, z)$ , respectively. The equivalent isotropic source locates inwardly with respect to the curvature of the interface at $(R_{0} - R_{a}, ϕ^{'}, z^{'})$ , and the image of it with respect to the boundary locates at a to-be-determined radial position of $(ρ^{?}, ϕ^{'}, z^{'})$ , which will reach $(R_{0} + R_{a} + 2 R_{b}, ϕ^{'}, z^{'})$ as the radius of the cylinder increases toward infinity. (d) Spherical geometry, where the source and the detector locate at $(R_{0}, θ^{'}, ϕ^{'})$ and $(R_{0}, θ, ϕ)$ , respectively. The equivalent isotropic source locates inwardly with respect to the curvature of the interface at $(R_{0} - R_{a}, θ^{'}, ϕ^{'})$ , and the image of it with respect to the boundary locates at a to-be-determined radial position of $(r^{?}, θ^{'}, ϕ^{'})$ , which will reach $(R_{0} + R_{a} + 2 R_{b}, θ^{'}, ϕ^{'})$ as the radius of the cylinder increases toward infinity.

2.1.

Infinite Geometry and Semi-Infinite Geometry

A homogeneous medium of infinite geometry is illustrated in Fig. 1(a), wherein the source is regarded as being isotropic. In the infinite geometry, it is convenient to set the source at ${\vec{χ}}^{'} = (0,0, 0)$ , and the detector at $\vec{χ} = (0, d, 0)$ , where $d$ is the source–detector separation.

A homogeneous medium of semi-infinite geometry with the physical source and the detector located on the medium boundary is depicted in Fig. 1(b). In the semi-infinite geometry, it becomes convenient to set the physical source at ${\vec{χ}}^{'} = (0,0, 0)$ , and the detector at $\vec{χ} = (0, d, 0)$ . The physical source at ${\vec{χ}}^{'} = (0,0, 0)$ launches photon into the medium at an initial direction orthogonal to the medium-applicator interface, and the resulted photon diffusion is treated as being excited by an equivalent “real” isotropic source located at one step of reduced scattering, represented by $R_{a} = 1 / μ_{s}^{'}$ , into the medium. The “real” isotropic source ${\vec{χ}}_{real}^{'}$ thus has the coordinates of $(0,0, - R_{a})$ . The effect of the medium-applicator interface on photon diffusion as modeled by the EBC²¹^,²³^,²⁴ sets zero the photon fluence rate on an imaginary boundary located at a distance of $R_{b} = 2 A D$ away from the physical boundary, where $A$ is a coefficient²³ related to the refractive index differences across the physical boundary. This EBC is conventionally accommodated by setting a sink or a negative “image” of the “real” isotropic source, with respect to the extrapolated boundary. The “image” source for the semi-infinite geometry has the opposite strength of the “real” isotropic source and locates at ${\vec{χ}}_{i m a g}^{'}$ , which has the coordinates of $(0,0, R_{a} + 2 R_{b})$ . The distances from the detector to the “real” isotropic source ${\vec{χ}}_{real}^{'}$ and the “image” source ${\vec{χ}}_{imag}^{'}$ are denoted by $l_{r}$ and $l_{i}$ , respectively. The notations of $l_{r} = | \vec{χ} - {\vec{χ}}_{real}^{'} |$ and $l_{i} = | \vec{χ} - {\vec{χ}}_{imag}^{'} |$ also apply to other medium geometries studied in this work when involving a boundary. In all studied geometries, the straight distance between the physical source and the detector $d = | \vec{χ} - {\vec{χ}}^{'} |$ is referred to as the “line-of-sight” source–detector distance, and hereafter specified as source–detector distance. It is with respect to this $d$ , which is to remain identical across different medium geometries considered in this work, that the DPF is to be expressed, facilitating the evaluation of the effect of shape and dimension of the nonsemi-infinite-domains on DPF in comparison to the semi-infinite geometry on DPF.

2.2.

Infinite-Length Cylindrical Geometry

The homogeneous medium bounded in an infinite-length cylindrical geometry is illustrated in Fig. 1(c). The radius of the cylinder is $R_{0}$ , then a detector $\vec{χ}$ on the medium-applicator interface locates at $(R_{0}, ϕ, z)$ , and a directional source ${\vec{χ}}^{'}$ on the medium-applicator interface locates at $(R_{0}, ϕ^{'}, z^{'})$ . The geometric symmetry requires that the directional source ${\vec{χ}}^{'}$ be modeled by a “real” isotropic source ${\vec{χ}}_{real}^{'}$ positioned along the radial direction of ${\vec{χ}}^{'}$ and inwardly at a radial distance of $R_{a} = 1 / μ_{s}^{'}$ , at $(R_{0} - R_{a}, ϕ', z')$ . The extrapolated zero-boundary is concentric with and at a radial distance of $R_{b} = 2 A D$ outward from the physical boundary. Apparently, as the radius $R_{0}$ of the infinite-length cylinder reaches infinity, the photon fluence rate associated with a source–detector pair on the infinite-length cylindrical medium-applicator interface shall reach that on the semi-infinite medium interface for the same source–detector distance. This characteristic has served for “checking” the results developed for infinite-length cylindrical geometry.²⁵

If one considers reaching the curved geometry from a semi-infinite geometry that extends to infinities at two orthogonal dimensions of the direction, bending on one dimension of the direction will make the effect of an infinite-length cylinder, and concurrent bending on both dimensions of the direction will reach the effect of a sphere. From such an intuitive perspective, we could expect that the effect of the one-dimensional (1-D) bending of the semi-infinite domain on photon diffusion will be seen when evaluating the photon diffusion along the azimuthal direction of the resulted infinite-length cylinder geometry, and that effect on photon diffusion when compared to semi-infinite geometry will be less than the double-dimensional bending of the semi-infinite domain to reaching the sphere for the same source–detector distances. Consequently, we evaluate the infinite-cylinder domain in the azimuthal direction only.

2.3.

Spherical Geometry

A spherical geometry is illustrated in Fig. 1(d). The radius of the spherical applicator is $R_{0}$ , then a detector $\vec{χ}$ on the medium-applicator interface locates at $(R_{0}, θ, ϕ)$ , and a directional source ${\vec{χ}}^{'}$ on the medium-applicator interface locates at $(R_{0}, θ^{'}, ϕ^{'})$ . The geometric symmetry requires that the directional source ${\vec{χ}}^{'}$ be modeled by a “real” isotropic source ${\vec{r}}_{real}^{'}$ positioned along the radial direction of ${\vec{χ}}^{'}$ and inwardly at a radial distance of $R_{a} = 1 / μ_{s}^{'}$ , [i.e., at $(R_{0} - R_{a}, θ^{'}, ϕ^{'})$ ]. The extrapolated zero boundary is concentric with and at a radial distance of $R_{b} = 2 A D$ outward from the physical boundary. Apparently, as the radius $R_{0}$ of the spherical geometry reaches infinity, the photon fluence rate associated with a source–detector pair on the spherical medium-applicator interface shall reach that on the semi-infinite medium interface for the same source–detector distance. It shall be noted, however, that the rate at which the photon fluence approaches the value of the semi-infinite geometry will be different in spherical geometry comparing to the infinite-length cylindrical geometry, given the same radius.

2.4.

Steady-State Photon Diffusion in an Infinite Homogeneous Domain

The equation of steady-state photon diffusion in a homogeneous medium is expressed by

Eq. (1)

\nabla^{2} Ψ (\vec{χ}) - \frac{μ_{a}}{D} Ψ (\vec{χ}) = - \frac{S (\vec{χ})}{D},

where

Ψ

is the photon fluence rate and

S

is the source term. For an infinite medium and a point source at

{\vec{χ}}^{'}

with an intensity

S

as represented by

S \cdot δ ({\vec{χ}}^{'})

, Eq. (1) has the well-known closed-form solution of Ref. 25

Eq. (2)

Ψ (\overset{⇀}{χ}, {\vec{χ}}^{'}) = \frac{S}{4 π D} \frac{1}{| \vec{χ} - {\vec{χ}}^{'} |} e^{- μ_{eff} | \vec{χ} - {\vec{χ}}^{'} |},

where

μ_{eff} = \sqrt{μ_{a} / D} = \sqrt{3 μ_{s}^{'} μ_{a}}

is the effective attenuation coefficient. The solution to Eq. (1) also can be written in the form of the following eigenfunction expansion (for derivation, see Appendix A):

Eq. (3)

Ψ (\overset{⇀}{χ}, {\vec{χ}}^{'}) = \frac{S}{D} (μ_{eff}) \sum_{l = 0}^{\infty} i_{l} (μ_{eff} r_{<}) k_{l} (μ_{eff} r_{>}) \sum_{m = - l}^{l} Y_{l m}^{*} (θ^{'}, ϕ^{'}) Y_{l m} (θ, ϕ),

where

i_{l}

and

k_{l}

are respectively the

l

’th order modified spherical Bessel function of the first and the second kinds,

r_{<}

and

r_{>}

are respectively the smaller and greater radial coordinates of the source and detector, and

Y_{l m}

is the spherical harmonics function. The numerical identity between Eq. (3) and the well-known form of Eq. (2) is justified in Appendix B.

2.5.

Steady-State Photon Diffusion in a Homogeneous Semi-Infinite Medium

We follow the notations summarized by Fantini et al.³⁴ for the semi-infinite geometry having a directional source and an isotropic detector located on a planar boundary, as illustrated in Fig. 1(b). As stated in Sec. 2.1, the directional source is modeled as an isotropic source placed one reduced scattering distance into the medium. Then, based on the extrapolated boundary and the “image” source approach, the steady-state photon fluence rate reaching the detector located on the physical boundary is determined by the equivalent “real” isotropic source and the “image” of it with respect to the extrapolated boundary. That gives

Eq. (4)

Ψ = Ψ_{real} - Ψ_{imag} = \frac{S}{4 π D l_{r}} e^{- μ_{eff} l_{r}} - \frac{S}{4 π D l_{i}} e^{- μ_{eff} l_{i}},

where

Eq. (5)

l_{r} = d \sqrt{1 + {(R_{a} / d)}^{2}}, R_{a} = 1 / μ_{s}^{'},

Eq. (6)

l_{i} = d \sqrt{1 + {[(R_{a} + 2 R_{b}) / d]}^{2}}, R_{b} = 2 A D .

3. Steady-State Photon Diffusion in a Homogeneous Medium Bounded in a Spherical Geometry

The principles to find the solutions to steady-state photon diffusion associated with a homogeneous infinite-length cylinder²⁵ will not be duplicated here. The solutions associated with infinite-cylinder geometry will however be presented when needed for the comparison with the solutions obtained for other geometries.

3.1.

Solution of Steady-State Photon Diffusion Associated with a Homogeneous Medium of Spherical Geometry

In surface imaging of a medium of spherical geometry with a radius of $R_{0}$ , the physical source is located at $(R_{0}, θ^{'}, ϕ^{'})$ and the detector is located at $(R_{0}, θ, ϕ)$ , both on the physical boundary. The equivalent “real” isotropic source locates at $(R_{0} - R_{a}, θ^{'}, ϕ^{'})$ , and the extrapolated boundary locates at a radial position of $ρ_{r >} = R_{0} + 2 A D$ outside the physical boundary. The geometric symmetry determines that the image of the “real” isotropic source with respect to the extrapolated boundary must locate along the radial direction of the “real” isotropic or the physical source. This image source and the “real” isotropic source collectively set zero the photon fluence rate on the extrapolated boundary.

Based on Eq. (3), the photon fluence rate associated with the “real” isotropic source and evaluating on the extrapolated boundary, for which the source locates at $r_{<} = R_{0} - R_{a}$ and the detector locates at $r_{>} = R_{0} + R_{b}$ , is

Eq. (7)

{Ψ_{real} |}_{extr} = \frac{1}{D} (μ_{eff}) \sum_{l = 0}^{\infty} S \cdot i_{l} [μ_{eff} (R_{0} - R_{a})] k_{l} [μ_{eff} (R_{0} + R_{b})] \sum_{m = - l}^{l} Y_{l m}^{*} (θ^{'}, ϕ^{'}) Y_{l m} (θ, ϕ),

where the notation “

left | right

” indicates evaluating the “left” as a source on the “right” as a boundary. Note that the source term “

S

” is placed inside the summation. This basically indicates that any “

l

”’th order of the “real” source has the same intensity of “

S

,”,as compared to different intensity for each order of the “image” source as will be shown in the following equation. Similarly, the photon fluence rate associated with the image source and evaluating on the extrapolated boundary, for which the source locates at radial position of a to-be-determined

r_{>}

and the detector locates at

r_{<} = R_{0} + R_{b}

, is

Eq. (8)

Ψ_{imag} |_{extr} = \frac{1}{D} (μ_{eff}) \sum_{l = 0}^{\infty} S_{l}^{*} i_{l} [μ_{eff} (R_{0} + R_{b})] k_{l} [μ_{eff} r_{>}] \sum_{m = - l}^{l} Y_{l m}^{*} (θ^{'}, ϕ^{'}) Y_{l m} (θ, ϕ),

where the

S_{l}^{*}

terms associated with each “

l

”’th order are not the same. Based on the essence of “image-source,”³⁵^,³⁶ the two unknown terms

S_{l}^{*}

and

r_{>}

associated with the

l

’th order “image” source (the

k_{l}

component) can be expressed by a single unknown term

S_{l}

associated with the same order “real” source (the

i_{l}

component), that is

Eq. (9)

S_{l}^{*} k_{l} (μ_{eff} r_{>}) = S_{l} i_{l} [μ_{eff} (R_{0} - R_{a})] .

Applying Eq. (9) to the EBC of $Ψ_{real} |_{extr} - Ψ_{imag} |_{extr} = 0$ , we have

Eq. (10)

S_{l} = S \frac{k_{l} [μ_{eff} (R_{0} + R_{b})]}{i_{l} [μ_{eff} (R_{0} + R_{b})]} l = 0,1, 2, \dots

Now for the “real” source but evaluating at the physical boundary, the source still locates at $r_{<} = R_{0} - R_{a}$ but the detector locates at $r_{>} = R_{0}$ . For the “image” source also evaluating at the physical boundary, the detector locates at $r_{<} = R_{0}$ and the source terms are now known through Eqs. (9) and (10). Collectively, the photon fluence rate sensed by a detector at the physical boundary becomes

Eq. (11)

Ψ = Ψ_{real} |_{phys} - Ψ_{imag} |_{phys} = \frac{S}{D} (μ_{eff}) \sum_{l = 0}^{\infty} i_{l} [μ_{eff} (R_{0} - R_{a})] k_{l} [μ_{eff} R_{0}] \sum_{m = - l}^{l} Y_{l m}^{*} (θ^{'}, ϕ^{'}) Y_{l m} (θ, ϕ) - \frac{S}{D} (μ_{eff}) \sum_{l = 0}^{\infty} i_{l} [μ_{eff} (R_{0} - R_{a})] i_{l} [μ_{eff} R_{0}] \frac{k_{l} [μ_{eff} (R_{0} + R_{b})]}{i_{l} [μ_{eff} (R_{0} + R_{b})]} \sum_{m = - l}^{l} Y_{l m}^{*} (θ^{'}, ϕ^{'}) Y_{l m} (θ, ϕ) = \frac{S}{D} (μ_{eff}) \sum_{l = 0}^{\infty} i_{l} [μ_{eff} (R_{0} - R_{a})] k_{l} [μ_{eff} R_{0}] \sum_{m = - l}^{l} Y_{l m}^{*} (θ^{'}, ϕ^{'}) Y_{l m} (θ, ϕ) ⟨ 1 - \frac{i_{l} [μ_{eff} R_{0}]}{k_{l} [μ_{eff} R_{0}]} \frac{k_{l} [μ_{eff} (R_{0} + R_{b})]}{i_{l} [μ_{eff} (R_{0} + R_{b})]} ⟩ .

Equation (11) is convertible to the following form that may be more convenient for numerical evaluation:

Eq. (12)

Ψ = \frac{S}{D} \sqrt{\frac{1}{(R_{0} - R_{a}) R_{0}}} \sum_{l = 0}^{\infty} I_{l + 1 / 2} [μ_{eff} (R_{0} - R_{a})] K_{l + 1 / 2} (μ_{eff} R_{0}) \frac{2 l + 1}{4 π} P_{l} (\cos γ) ⟨ 1 - \frac{I_{l + 1 / 2} [μ_{eff} R_{0}]}{K_{l + 1 / 2} [μ_{eff} R_{0}]} \frac{K_{l + 1 / 2} [μ_{eff} (R_{0} + R_{b})]}{I_{l + 1 / 2} [μ_{eff} (R_{0} + R_{b})]} ⟩,

where

I_{l}

and

K_{l}

are respectively the

l

’th order modified Bessel function of the first and the second kinds,

P_{l}

is the

l

’th order Legendre polynomial, and

γ

is the angle formed between the position vector of the source and that of the detector (see also Appendix A).

The value determined by Eq. (11) or Eq. (12) is compared against FEM results (configuration of FEM will be detailed in subsequent sections), as shown in Fig. 2, for two radii of 25 mm and 50 mm at three sets of optical properties: the top traces correspond to $μ_{a} = 0.01 {mm}^{- 1}$ and $μ_{s}^{'} = 0.25 {mm}^{- 1}$ , the middle traces $μ_{a} = 0.01 {mm}^{- 1}$ and $μ_{s}^{'} = 0.50 {mm}^{- 1}$ , and the bottom traces $μ_{a} = 0.01 {mm}^{- 1}$ and $μ_{s}^{'} = 1.0 {mm}^{- 1}$ . A total of 51 detector positions are placed at a source–detector distance of 1 to 50 mm at a step of 1 mm from a single source. At the radius of 25 mm, the farthest detector point is positioned at the opposite to the source point, and at the radius of 50 mm, the farthest detector point is positioned at 60 deg with respect to the source point. The results are presented in terms of $\ln (d^{2} * Ψ)$ versus $d$ , as it leads to a convenient near-linear relationship as the source–detector distance becomes longer than approximately five times of the reduced scattering pathlength. The oscillation observed of the analytical result corresponding to $μ_{s}^{'} = 1.0 {mm}^{- 1}$ is due to arithmetic overflow in computing the modified Bessel functions. It is shown that at source–detector separations less than approximately five times of the reduced scattering pathlength (i.e., in the subdiffusion regime) there is a discrepancy between the analytically evaluated photon fluence rate and its FEM counterpart. As observed in our previous study³⁷ dedicated to experimental examination of the analytical results of photon fluence rate specific to infinite-length cylinder,²⁵ the discrepancy between the diffusion-based analytical quantification and the FEM simulation in the subdiffusion regime was largely attributed to the difference between a spatially impulsive source treated in the analytical solution and a spatially distributed Gaussian source profile utilized in the FEM tool available for this work. The spatially impulsive source treated in the analytical solution was physically analogous to the single-point launching of forward directional photons into the medium, whereas the spatially distributed Gaussian source profile implemented in the FEM had the size effect manifested increasingly at shorter distances from the source. Additionally, the treatment to boundary effect (e.g., the probability of a boundary-reaching photon returning to the medium) in the analytical model was based on EBC, whereas a Robin-type boundary condition was programmed in the FEM; the subtle difference in the boundary conditions could have been amplified at the subdiffusion regime. In the diffusion regime, however, the analytical computations marked by the solid lines and the FEM results given in discrete markers clearly match, corroborating that Eq. (11) is the solution to the steady-state photon diffusion in a homogeneous medium bounded in a sphere.

Fig. 2

Comparison of the value determined by Eq. (11) or (12) against finite-element method (FEM) results. The parameters are chosen as: radii of 25 mm and 50 mm, optical properties of $μ_{a} = 0.01 {mm}^{- 1}$ with $μ_{s}^{'} = 0.25 {mm}^{- 1}$ (top traces), $μ_{s}^{'} = 0.5 {mm}^{- 1}$ (middle traces), and $μ_{s}^{'} = 1.0 {mm}^{- 1}$ (bottom traces). A total of 51 detector positions are placed at a source–detector distance of 1 to 50 mm at a step of 1 mm from the single source position. The results given by FEM are relatively smooth for all sets of parameters considered in contrast to the variation of the numerical computation based on Eq. (12) for the 50 mm radius at $μ_{s}^{'} = 1.0 {mm}^{- 1}$ , but nevertheless the numerical computations marked by the solid lines clearly follow the traces indicated by the FEM results.

3.2.

Large-Diameter Approximation of the Solution for a Spherical Domain

As shown in Fig. 1(d), if a plane tangential to the sphere at the physical source position is considered as an imaginary planar boundary to a semi-infinite geometry of the medium to which the medium of the spherical geometry will approach as the radius of the sphere increases without bound, then the “real” isotropic source in this semi-infinite geometry still locates at $(R_{0} - R_{a}, θ^{'}, ϕ^{'})$ , but the image of the “real” isotropic source with respect to this “hypothetical” semi-infinite boundary will be at $(R_{0} + R_{a} + 2 R_{b}, θ^{'}, ϕ^{'})$ . According to Eq. (3), the photon fluence rate sensed by a detector on the spherical boundary due to the image of the “real” isotropic source associated with this “hypothetical” semi-infinite boundary is

Eq. (13)

Ψ_{imag}^{semi} |_{phys} = \frac{S}{D} (μ_{eff}) \sum_{l = 0}^{\infty} i_{l} [μ_{eff} R_{0}] k_{l} [μ_{eff} (R_{0} + R_{a} + 2 R_{b})] \sum_{m = - l}^{l} Y_{l m}^{*} (θ^{'}, ϕ^{'}) Y_{l m} (θ, ϕ) .

The photon fluence rate sensed by a detector on the spherical boundary due to the image of the “real” isotropic source associated with the spherical boundary, as seen in Eq. (11), can be rewritten to

Eq. (14)

Ψ_{imag} |_{phys} = \frac{S}{D} (μ_{eff}) \sum_{l = 0}^{\infty} i_{l} [μ_{eff} R_{0}] k_{l} [μ_{eff} (R_{0} + R_{a} + 2 R_{b})] \cdot η_{l} \cdot \sum_{m = - l}^{l} Y_{l m}^{*} (θ^{'}, ϕ^{'}) Y_{l m} (θ, ϕ),

where

Eq. (15)

η_{l} = \frac{k_{l} [μ_{eff} (R_{0} + R_{b})]}{i_{l} [μ_{eff} (R_{0} + R_{b})]} \frac{i_{l} [μ_{eff} (R_{0} - R_{a})]}{k_{l} [μ_{eff} (R_{0} + R_{a} + 2 R_{b})]} .

Note that $R_{0} ≫ R_{a}$ , $R_{b}$ and if the argument of the modified Bessel functions is much greater than 1 (i.e., $μ_{eff} R_{0} \geq 10$ that corresponds to $R_{0} \geq 57.8 mm$ for $μ_{a} = 0.01 {mm}^{- 1}$ and $μ_{s}^{'} = 1.0 {mm}^{- 1}$ ) the modified spherical Bessel functions in Eq. (15) can be simplified by their asymptotic expressions,³⁸ then Eq. (15) becomes

Eq. (16)

η_{l} = \frac{R_{0} + R_{a} + 2 R_{b}}{R_{0} - R_{a}} .

Comparing to its counterpart for infinite cylinder [Eq. (3.1.9) in Ref. 25], Eq. (16) for sphere geometry is different only in the power of the entire fraction, that is, one or one-half—an intuitive analogy of the double-sided bending or single-sided bending from a semi-infinite geometry. Substituting Eq. (16) into Eq. (14) and comparing to Eq. (13) we have

Eq. (17)

Ψ_{imag} |_{phys} = Ψ_{imag}^{semi} |_{phys} \cdot \frac{R_{0} + R_{a} + 2 R_{b}}{R_{0} - R_{a}} .

Hence, for a sphere of large diameter, Eq. (12) approximates to

Eq. (18)

Ψ = Ψ_{real} |_{phys} - Ψ_{imag} |_{phys} = Ψ_{real} |_{phys} - Ψ_{imag}^{semi} |_{phys} \cdot \frac{R_{0} + R_{a} + 2 R_{b}}{R_{0} - R_{a}} .

The $Ψ_{real} |_{phys}$ of Eq. (18) is essentially the same as the $Ψ_{real}$ in Eq. (4). As the radius of the sphere reaches infinity without bound, we have $(R_{0} + R_{a} + 2 R_{b}) / (R_{0} - R_{a}) \to 1$ , and the $Ψ_{imag}^{semi} |_{phys}$ of Eq. (18) will become the $Ψ_{imag}$ in Eq. (4) because the detector point locating at $(R_{0}, θ, ϕ)$ will reach the “hypothetical” semi-infinite boundary. By using the expression of the photon fluence rate given in Eq. (2), we can rewrite Eq. (18) as

Eq. (19)

Ψ = \frac{S}{4 π D} \frac{e^{- μ_{eff} l_{r}}}{l_{r}} - \frac{S}{4 π D} \frac{e^{- μ_{eff} l_{i}}}{l_{i}} \frac{R_{0} + R_{a} + 2 R_{b}}{R_{0} - R_{a}} .

For the purpose of assessing the general patterns of the steady-state photon fluence rate versus the source–detector distance, the following expressions of the involved distances will be implemented in Eq. (19):

Eq. (20)

l_{r} = d \sqrt{1 + Δ_{r}}; where Δ_{r} = \frac{R_{a}^{2}}{d^{2}} - \frac{R_{a}}{R_{0}},

Eq. (21)

l_{i} = d \sqrt{1 + Δ_{i}}; where Δ_{i} = \frac{{(R_{a} + 2 R_{b})}^{2}}{d^{2}} + \frac{R_{a} + 2 R_{b}}{R_{0}} .

The formula for the photon fluence rates associated with the four geometries of consideration, including infinite medium, semi-infinite medium or half-plane, large-radius infinite-cylinder when evaluated along the azimuthal direction, and large-radius sphere is tabulated in Table 1. The numerical results shown in the later sections performed for conditions including $μ_{eff} R_{0} = 8.66$ will indicate that the approximate formula for the infinite cylinder and the sphere is generally valid when the radius satisfies $μ_{eff} R_{0} \geq 10$ , that is, the radius (diameter) is more than 10 (20) times of the effective attenuation pathlength $1 / μ_{eff}$ .

Table 1

Analytical (large-radius approximation) representations of the photon diffusion.

Geometry	Fluence rate	Distance parameters
Infinite	$Ψ_{\inf} = \frac{S}{4 π D} \frac{e^{- μ_{eff} d}}{d}$	$d$
Semi-infinite	$Ψ_{semi} = \frac{S}{4 π D} \frac{e^{- μ_{eff} l_{r}}}{l_{r}} - \frac{S}{4 π D} \frac{e^{- μ_{eff} l_{i}}}{l_{i}}$	$l_{r} = d \sqrt{1 + Δ_{r}}$ @ $R_{0} \to \infty$ $l_{i} = d \sqrt{1 + Δ_{i}}$ @ $R_{0} \to \infty$
Large infinite cylinder, azimuthal direction	$Ψ_{azi} = \frac{S}{4 π D} \frac{e^{- μ_{eff} l_{r}}}{l_{r}} - \frac{S}{4 π D} \frac{e^{- μ_{eff} l_{i}}}{l_{i}} \sqrt{\frac{R_{0} + R_{a} + 2 R_{b}}{R_{0} - R_{a}}}$	$l_{r} = d \sqrt{1 + Δ_{r}}$ $l_{i} = d \sqrt{1 + Δ_{i}}$
Large sphere	$Ψ_{sph} = \frac{S}{4 π D} \frac{e^{- μ_{eff} l_{r}}}{l_{r}} - \frac{S}{4 π D} \frac{e^{- μ_{eff} l_{i}}}{l_{i}} \frac{R_{0} + R_{a} + 2 R_{b}}{R_{0} - R_{a}}$	$l_{r} = d \sqrt{1 + Δ_{r}}$ $l_{i} = d \sqrt{1 + Δ_{i}}$

4. Differential Pathlength Factors Associated with Homogeneous Semi-Infinite, Infinite-Length Cylinder When Evaluated along the Azimuthal Direction and Spherical Geometries

The DPF was originally introduced in Ref. 6 to relate the changes of the light to the changes of the absorption coefficient, at a given source–detector pair, by the modified Beer–Lambert law of $Att = \exp [- μ_{a} \cdot DPF \cdot d]$ , where “Att” is the ratio of measured light over the incident light. We use the definition of DPF provided by Fantini et al.⁹ or Arridge et al.¹⁸ that becomes $DPF = {\partial [\ln (S / Ψ)] / \partial μ_{a}} / d$ using the set of equations in Table 1. The DPF for infinite medium geometry is obtained by direct derivation operation, whereas the DPFs for the semi-infinite, cylindrical when evaluated along the azimuthal direction, and spherical geometries are obtained by involving Taylor series expansion. The results are:

• Infinite geometry [identical to Eq. (4) of Ref. 9]
Eq. (22)
${DPF}_{\inf} = \frac{1}{d} \frac{\partial}{\partial μ_{a}} [\ln (\frac{S}{Ψ_{\inf}})] = \frac{\partial}{\partial μ_{a}} (μ_{eff}) = \frac{\sqrt{3 μ_{s}^{'}}}{2 \sqrt{μ_{a}}} .$
• Semi-infinite geometry [identical to Eq. (9) of Ref. 9]
Eq. (23)
${DPF}_{semi} = \frac{1}{d} \frac{\partial}{\partial μ_{a}} [\ln (\frac{S}{Ψ_{semi}})] = \frac{1}{d} \frac{\exp (- μ_{eff} l_{r}) - \exp (- μ_{ef f} l_{i})}{\frac{\exp (- μ_{eff} l_{r})}{l_{r}} - \frac{\exp (- μ_{eff} l_{i})}{l_{i}}} \frac{\partial}{\partial μ_{a}} (μ_{eff}) \approx \frac{\sqrt{3 μ_{s}^{'}}}{2 \sqrt{μ_{a}}} \frac{d \sqrt{3 μ_{s}^{'} μ_{a}}}{d \sqrt{3 μ_{s}^{'} μ_{a}} + 1} .$
• Infinite-length cylindrical geometry, evaluated along the azimuthal direction
Eq. (24)
${DPF}_{azi} = \frac{1}{d} \frac{\partial}{\partial μ_{a}} [\ln (\frac{S}{Ψ_{azi}})] = \frac{1}{d} \frac{\exp (- μ_{eff} l_{r}) - [\exp (- μ_{eff} l_{i})] \sqrt{\frac{R_{0} + R_{a} + 2 R_{b}}{R_{0} - R_{a}}}}{\frac{\exp (- μ_{eff} l_{r})}{l_{r}} - \frac{\exp (- μ_{eff} l_{i})}{l_{i}} \sqrt{\frac{R_{0} + R_{a} + 2 R_{b}}{R_{0} - R_{a}}}} \frac{\partial}{\partial μ_{a}} (μ_{eff}) = \frac{\sqrt{3 μ_{s}^{'}}}{2 \sqrt{μ_{a}}} \frac{d \sqrt{3 μ_{s}^{'} μ_{a}} [η_{azi} Δ_{i} - Δ_{r}] - [η_{azi} - 1]}{(d \sqrt{3 μ_{s}^{'} μ_{a}} + 1) [η_{azi} Δ_{i} - Δ_{r}] - [η_{azi} - 1]},$
Eq. (25)
$where η_{azi} = \sqrt{\frac{R_{0} + R_{a} + 2 R_{b}}{R_{0} - R_{a}}} .$
• Spherical geometry
Eq. (26)
${DPF}_{sph} = \frac{1}{d} \frac{\partial}{\partial μ_{a}} [\ln (\frac{S}{Ψ_{sph}})] = \frac{1}{d} \frac{\exp (- μ_{eff} l_{r}) - [\exp (- μ_{eff} l_{i})] \frac{R_{0} + R_{a} + 2 R_{b}}{R_{0} - R_{a}}}{\frac{\exp (- μ_{eff} l_{r})}{l_{r}} - \frac{\exp (- μ_{eff} l_{i})}{l_{i}} \frac{R_{0} + R_{a} + 2 R_{b}}{R_{0} - R_{a}}} \frac{\partial}{\partial μ_{a}} (μ_{eff}) = \frac{\sqrt{3 μ_{s}^{'}}}{2 \sqrt{μ_{a}}} \frac{d \sqrt{3 μ_{s}^{'} μ_{a}} [η_{sph} Δ_{i} - Δ_{r}] - [η_{sph} - 1]}{(d \sqrt{3 μ_{s}^{'} μ_{a}} + 1) [η_{sph} Δ_{i} - Δ_{r}] - [η_{sph} - 1]},$
Eq. (27)
$where η_{sph} = \frac{R_{0} + R_{a} + 2 R_{b}}{R_{0} - R_{a}} .$

Note that Eq. (23) is the same DPF for semi-infinite medium as was used by Fantini et al.,⁹ and both Eqs. (24) and (26) will reach Eq. (23) as the radius $R_{0}$ becomes infinity.

5. Numerical Evaluation of the Differential Pathlength Factors Associated with Homogeneous Semi-Infinite, Infinite Cylinder Evaluated along the Azimuthal Direction and Spherical Geometries

The patterns of DPF among the three specific configurations associated respectively with semi-infinite geometry, infinite-length cylinder geometry when evaluated along the azimuthal direction, and sphere were examined by using near-infrared frequency-domain absorption and scattering tomography (NIRFAST),³⁹^,⁴⁰ an FEM platform for solving photon diffusion problems. A half-plane for semi-infinite geometry and an infinite-length cylinder cannot be produced by FEM, but could be reasonably represented in FEM if the length dimension remained substantially greater than the source–detector spacings. The cylinder and the sphere domains were examined at two relatively small radii of $R_{0} = 25 mm$ and $R_{0} = 50 mm$ with consistent node density of 1 mm that was comparable to the reduced scattering pathlength $(1 / μ_{s}^{'})$ . The mesh for the semi-infinite geometry was generated for a $100 mm \times 100 mm$ surface area and 50 mm depth, the mesh for the infinite-length cylinder geometry at a 100 mm length with the aforementioned 25 mm or 50 mm radius, and the mesh for the sphere at the same set of radius. Figure 3 displays the meshes (cylindrical and spherical meshes at 25 mm radius only).

Fig. 3

FEM meshes for the configurations to be compared. In the semi-infinite geometry, one source is set in the middle-sagittal plane and at 3 cm away from one edge, and a total of 50 detector points at an interval of 1 mm are placed in the middle-sagittal plane to cover a source–detector range from 1 mm to 50 mm. In the sphere geometry and the cylinder geometry evaluating along the azimuthal direction, one source is at the middle azimuthal plane, and a total of 50 detector points are placed in the middle azimuthal plane to have a range of the source–detector distance (“line-of-sight” or the chord distance) from 1 mm to 50 mm with an interval of 1 mm. In these two configurations, the detector with 50 mm distance from the source is set at the direct opposite of the source point.

In the configuration for the semi-infinite geometry, one source point was set in the middle-sagittal plane and at 30 mm away from one edge of the domain. A total of 50 detector points at an interval of 1 mm were placed in the middle-sagittal plane to cover a source–detector range from 1 mm to 50 mm. In the configurations for the sphere geometry and the cylinder geometry but evaluating along the azimuthal direction, one source point was set at the middle azimuthal plane, and a total of 50 detector points were placed in the middle azimuthal plane to have a range of the source–detector distance (“line-of-sight” distance or the chord distance) from 1 mm to 50 mm at an interval of 1 mm.

The DPFs were calculated using the forward solver of NIRFAST by the following steps: (1) the first set of the photon fluence data $Ψ_{1}$ measured at the 50 detector positions were forward calculated at the following setting of the optical properties: $μ_{a} = 0.01 {mm}^{- 1}$ , $μ_{s}^{'} = 1.0 {mm}^{- 1}$ . (2) The second set of the photon fluence data $Ψ_{2}$ measured at the same set of 50 detector positions were forward calculated with the $μ_{a}$ at $0.011 {mm}^{- 1}$ and other parameters remaining unchanged from the values of the previous step. (3) The DPF was calculated as ${[\ln (Ψ_{1} / Ψ_{2})] / Δ μ_{a}} / d$ . Note that NIRFAST by default includes $μ_{a}$ in the calculation of the diffusion coefficient. To make the FEM evaluation as consistent with the analytical approaches as would be available, the $μ_{a}$ has been excluded from the calculation of the diffusion coefficient when using NIRFAST through the workspace open to the users. The DPFs for the source–detector distances ranging from 1 mm to 50 mm in the three specified configurations are shown in the left column of Fig. 4.

Fig. 4

Comparison of the differential pathlength factors (DPFs) for the source–detector separation between 1 mm and 50 mm at three configurations, evaluated using FEM (left column) and the set of Eq. (27) (right column). The set properties are $μ_{a} = 0.01 {mm}^{- 1}$ , $μ_{s}^{'} = 1.0 {mm}^{- 1}$ , refractive index of the medium is 1.33 (the same for the following figures), and $R_{0} = 25 mm$ and 50 mm. The inset at the top row compares the DPF of the sphere calculated using FEM and using Eq. (25).

The DPFs evaluated according to the set of Eqs. (23), (24), and (26) are given in the right column of Fig. 4. The equations were evaluated with the following parameters: $μ_{a} = 0.01 {mm}^{- 1}$ , $μ_{s}^{'} = 1.0 {mm}^{- 1}$ , $R_{0} = 25 mm$ , and $R_{0} = 50 mm$ over the source–detector distance from 3 mm to 50 mm. From both columns of Fig. 4, one may observe the following qualitative patterns of DPF at the same source–detector separation: (1) the DPF evaluated on a semi-infinite medium boundary is greater than the DPF evaluated along the azimuthal direction of the infinite-length cylindrical medium boundary; (2) the latter DPF of (1) is greater than the DPF evaluated on a spherical boundary, given the same radial dimension for the cylinder and the sphere. It is also obvious from Eqs. (23), (24), and (26) that as $R_{0} \to \infty$ , Eqs. (24) and (26) will all reach Eq. (23), whereas at different rates. The inset at the top row illustrates how the DPF of the sphere calculated using Eq. (26), which is an approximation, fairs with the one calculated with FEM. The analytically approximated DPF for the spherical geometry of 50 mm radius at the given set of optical properties (resulting in $μ_{eff}$ $R_{0} = 8.66$ ) reached an average error of not more than 5% over the 10 mm to 25 mm source–detector distances when compared with the FEM results performed for the same spherical geometry. If the analytic DPF of the semi-infinite geometry were to be compared to, the DPF approximated for the spherical geometry of 50 mm radius at the set of optical properties as used for Fig. 4 has shown a maximum difference of about 9% at a source–detector distance of 9 mm as shown in Fig. 5. The relative difference of the DPF approximated for the spherical geometry with respect to the analytic semi-infinite one keeps reducing as the radius of the sphere increases—less than 5% at a radius of 110 mm for the given set of optical properties.

6. Discussions

The dependence of DPF on age has been well documented. In 1996, Duncan et al. showed that the older the subject was the larger the DPF became.¹² At 832 nm, the DPFs measured by using phase-resolved NIRS from the head of neonates were $4.67 \pm 0.65$ , while the DPFs of adults were $5.86 \pm 0.98$ ,⁷ at the same source–detector separation of 4.3 cm. The increase of DPF from infant to adult may be attributed to either a change of the optical properties of the brain, or the change of the size of the head, or a combination of both. For the appreciation of the effect of the changes of optical properties on DPF, we evaluated the DPFs of a semi-infinite geometry and a spherical geometry of 50 mm radius, using the set of Eqs. (23), (24), and (26), at three different combinations of optical properties. The first case is to fix the absorption coefficient at $μ_{a} = 0.01 {mm}^{- 1}$ , with the $μ_{s}^{'}$ changing from $0.5 {mm}^{- 1}$ to $1.5 {mm}^{- 1}$ . The results are shown in Fig. 6, where it is observed that the three folds of increase of $μ_{s}^{'}$ at the specific fixed value of $μ_{a}$ nearly doubled the DPFs for both the semi-infinite geometry (the upper or blue profile) and the spherical geometry (the lower or green profile). It is also observed from Fig. 6 that the absolute difference of the DPFs between the semi-infinite and spherical geometry seems to be insensitive to the $μ_{s}^{'}$ . The second case is to fix the reduced scattering coefficient at $μ_{s}^{'} = 1.0 {mm}^{- 1}$ , with the $μ_{a}$ changing from $0.01 {mm}^{- 1}$ to $0.1 {mm}^{- 1}$ . The results are shown in Fig. 7, where it is observed that the 10 folds of increase of $μ_{a}$ at the given specific fixed $μ_{s}^{'}$ has a nearly 50% reduction of the DPFs for both the semi-infinite geometry (the upper or blue profile) and the spherical geometry (the lower or green profile). It is also observed from Fig. 7 that the absolute difference of the DPFs between the semi-infinite and spherical geometry diminishes as the $μ_{a}$ increases. The third case is to fix the ratio of the absorption coefficient to the reduced scattering coefficient to $1 / 100$ , with the $μ_{s}^{'}$ changing from $0.5 {mm}^{- 1}$ to $1.5 {mm}^{- 1}$ . The results are shown in Fig. 8, where it is observed that if the absorption-to-reduced scattering ratio is kept the same, the three folds of increase of $μ_{s}^{'}$ has increased the DPFs approximately 20% for both the semi-infinite geometry (the upper or blue profile) and the spherical geometry. In comparison, the same change of $μ_{s}^{'}$ at a fixed $μ_{a}$ resulted in nearly 100% changes of DPF as shown in Fig. 6.

Fig. 5

The effect of the radius (in mm) of the spherical geometry on the DPF evaluated over 3 mm to 50 mm range, when compared to semi-infinite geometry for the same source-detector distances, are shown. The horizontal red dashed line marks the 5% error with respect to the semi-infinite case.

Fig. 6

The dependence of DPF on the absorption coefficient from 0.01 to $0.1 {mm}^{- 1}$ at a fixed reduced scattering coefficient of $1.0 {mm}^{- 1}$ . The upper blue profile corresponds to the DPF of a semi-infinite geometry, and the lower green profile a spherical geometry of 50 mm in radius, for source–detector distances from 1 mm to 50 mm in both geometries.

Fig. 7

The dependence of DPF on the reduced scattering coefficient from 0.5 to $1.5 {mm}^{- 1}$ at a fixed absorption-to-reduced scattering ratio of $1 / 100$ . The upper blue profile corresponds to the DPF of a semi-infinite geometry, and the lower green profile a spherical geometry of 50 mm in radius, for source–detector distances from 1 mm to 50 mm in both geometries.

The age dependence of the DPFs reported by a number of studies was summarized in Ref. 8. It was shown that DPF increased an amount of more than 1 from neonatal to 25 years of age, and an amount of slightly less than 1 from 25 years of age to 50 years of age.⁸ Ijichi et al.⁴¹ measured DPFs over a considerably narrow postconceptional range of age from 30 to 42 weeks by phase-resolved measurements. The DPF values $4.31 \pm 0.42$ at 835 nm and at source–detector separation of 2.9 cm were in good agreement with those reported by Duncan et al. at 832 nm. The optical properties of the neonatal brains had shown significant changes over the narrow range of the postconceptional ages. The average $μ_{s}^{'}$ at 835 nm changed from $0.43 {mm}^{- 1}$ to $0.79 {mm}^{- 1}$ for postconceptional age ranging from 30 to 42 weeks. The cerebral blood volume, which affected the global $μ_{a}$ , nearly doubled over the same period of postconceptional age. However, the DPF values had shown little if any change with the postconceptional age. The qualitative evaluations shown in Fig. 8 may provide some insights to the observations of Ijichi et al.⁴¹ on the no-or-little change of DPF over the range of postconceptional age of the subjects despite the significant changes of the scattering and absorption properties of the neonatal brain. If the near doubling of the reduced scattering coefficient were nearly counterbalanced by the near doubling of the absorption coefficient as a result of the increase of the cerebral hemoglobin content, the absorption-to-reduced scattering ratio could have remained relatively consistent for the subject. If that were the case, the variation of the DPFs of the subjects over the range of reduced scattering change could be small and not easily distinguishable from the inherent intersubject variations.

This work intends to evaluate how the shape as well as size of the diffusing medium (bounded in a regular geometry like a sphere) will affect the DPF. This initial line of evaluation is based on treating the medium as homogeneous. If one considers reaching the curved geometry from a semi-infinite geometry that extends to infinities at two dimensions of the direction, bending on one dimension of the direction will make the effect of an infinite-length cylinder, and bending on both dimensions of the direction will reach the effect of a spherical domain. From such an intuitive perspective, we could expect that the effect to photon diffusion due to a spherical curvature (as is the result of double-dimensional bending of the semi-infinite domain) likely will double that due to the azimuthal direction of an infinite-length cylinder geometry (as is the result of 1-D bending of the semi-infinite domain). The expressions of the photon fluence rate as given in the set of Eqs. (23), (24), and (26), which are resulted from the approximations of their respective general solutions, do show a quasidoubling effect of the spherical geometry when compared to the azimuthal direction of the infinite-length cylinder geometry. Such doubling effect of the spherical geometry compared to the azimuthal direction of the cylinder geometry with respect to the semi-infinite geometry is observable in Fig. 4, wherein the trend shown by the FEM becomes more pronounced at the radius of 50 mm compared to the radius of 25 mm. It is also evident that, as the radial dimension of either the spherical or cylindrical geometry increases, the deviation of the DPFs of the two curved geometries from that of the semi-infinite geometry diminishes. Therefore, in the cases of adult head that easily exceeds a local radius of 100 mm, the convention of using the DPFs estimated for a semi-infinite geometry is appropriate for source–detector distance greater than 2.5 cm; whereas for spherical geometries with sizes comparable to infant head, adoption of the DPFs corresponding to semi-infinite geometry should be used with caution.

With the development of denser optodes in fNIRS as well as the attention to neonatal applications, however, the noticeable changes of the DPF at smaller source–detector distances and at a spherical geometry comparing to the semi-infinite geometry may become too significant not to be accounted for when inferring clinical markers for intervention. As shown in Fig. 8, for example, the DPF corresponding to a 2 cm source–detector separation in a sphere of 50 mm radius as is comparable to the size of a neonatal head, can be approximately 5% smaller than that associated with a semi-infinite geometry at the optical properties of the medium set at $μ_{a} = 0.01 {mm}^{- 1}$ , $μ_{s}^{'} = 1.0 {mm}^{- 1}$ . A 5% reduction of the DPF translates to 5% more changes of absorption coefficient at the same level of fNIRS signal change. If the shape dependency of the DPF is not accounted for, there will be increased probability of underestimating the absorption changes by assuming the DPFs of semi-infinite medium, and an amount of 5% underestimation of the change of absorption coefficient may be a concern if the baseline parameter to be monitored is close to a clinically important threshold. We note that the DPFs reported by Ijichi et al.⁴¹ were experimentally determined by using the phase of the intensity modulated light, where the $μ_{a}$ and $μ_{s}^{'}$ were estimated by fitting the absolute DPF to a semi-infinite model. Therefore, when the sets of $μ_{a}$ and $μ_{s}^{'}$ were test plotted into a figure similar to Figs. 6 Fig. 7–8, unsurprisingly, they always appeared to agree with the semi-infinite model. Were the phase-resolved measurements fitted using a model that had accounted for the smaller head size and shape, the resulted absolute changes of the $μ_{a}$ and $μ_{s}^{'}$ of the neonatal brains could have been different.

Fig. 8

The relative value of the DPF approximated for the spherical geometry with respect to the analytic DPF of the semi-infinite geometry at $μ_{a} = 0.01 {mm}^{- 1}$ and $μ_{s}^{'} = 1.0 {mm}^{- 1}$ .

This study is limited as only the simplest homogeneous domain has been considered. Irrespective of the geometry dependence, the DPF is shown to decrease at short source–detector separation down to about 5 mm. This pattern of DPF associated with homogeneous domain is consistent with previous reports developed using MC calculations⁴² and the analytical solution of the diffusion equation incorporating the zero-boundary condition.¹⁸ However, there are experimental studies of the adult head⁴³ suggesting that the DPF increases at short source–detector separation ( $\sim 2 cm$ ). This trend of DPF increasing at short source–detector separation obtained from actual tissue was justified as relating to the inhomogeneity of tissue domain with which the use of partial differential pathlength factor (PDPF)⁴⁴ was articulated. An immediate improvement of the modeling of the tissue structure to accommodate PDPF is a multilayer geometry accounting for the effect of the blood content of cerebrospinal fluid on the optical properties of brain as indicated by Robertson et al.⁴⁵ It, nevertheless, could be expected that, the qualitative pattern of the geometry dependence of PDPF, specifically how it changes over different geometries at fixed optical properties and source–detector separation, may remain indifferent from that of DPF. Future works to extend the presented analytical approach to multilayer geometry are under planning. Furthering the analytical understanding of the geometry effect to the DPF and PDPF may eventually lead to quantitative, rapid, and accurate estimation of the $μ_{a}$ and $μ_{s}^{'}$ that would be particularly influential in the clinical evaluation and management of disorders associated with early brain development and acute-care cerebral monitoring.

7. Conclusion

This work analytically examined the DPF of steady-state photon diffusion in a homogeneous medium associated with a few regularly shaped geometries, including a half-plane, an infinite-length cylinder, and a sphere. With the analytical procedures established for steady-state photon diffusion in an infinite-length cylinder domain that can be extended to deriving DPF, the solution to steady-state photon diffusion in a spherical domain is developed by applying the “image” source method with the EBC and involving the modified spherical Bessel functions of the first and the second kinds. The solution is converted into a format employing the physical source and its image with respect to its associated semi-infinite geometry and a radius-dependent term accounting for the dimension of the sphere. The solutions can be simplified at the general condition of $μ_{eff}$ $R_{0} \geq 10$ , that is, the radius (diameter) is more than 10 (20) times of the effective attenuation pathlength, $1 / μ_{eff}$ , a form that involve regular functions without lengthy summations. With these analytical preparations, the DPFs of steady-state photon diffusion in homogeneous medium associated with the three bounded geometries, including the semi-infinite geometry, along the azimuthal direction of the infinite-length cylinder geometry, and the spherical geometry, are compared for the same line-of-sight source–detector distance. The patterns of DPF estimated for the two curved geometries with a radius of 25 mm or 50 mm agree in the global patterns with those given by FEM. Quantitatively, the DPF approximated for a source–detector separation range of 10 to 25 mm on a spherical medium of 50 mm radius with $μ_{a} = 0.01 {mm}^{- 1}$ and $μ_{s}^{'} = 1.0 {mm}^{- 1}$ is on average 5% different compared to the value obtained from FEM.

Appendices

Appendix A

Considering a source at ${\vec{x}}^{'}$ of $(r^{'}, θ^{'}, ϕ^{'})$ and a detector at $\vec{x}$ of $(r, θ, ϕ)$ in spherical coordinates, the equation for the Green’s function of Eq. (1) is

Eq. (28)

\nabla^{2} G (\vec{χ}, {\vec{χ}}^{'}) - μ_{eff}^{2} G (\vec{χ}, {\vec{χ}}^{'}) = - δ (\vec{χ} - {\vec{χ}}^{'}) .

The Dirac delta function in Eq. (28) expressed in spherical coordinates (Ref. 38, pages 120, 427) is

Eq. (29)

δ (\vec{χ} - {\vec{χ}}^{'}) = \frac{1}{r^{2}} δ (r - r^{'}) δ (\cos θ - \cos θ^{'}) δ (ϕ - ϕ^{'}),

where the delta functions involving both

θ

and

ϕ

can be decomposed to the spherical harmonics as [Ref. 38, page 108, Eq. (3.56), or page 120, Eq. (3.117)]

Eq. (30)

δ (\vec{χ} - {\vec{χ}}^{'}) = \frac{1}{r^{2}} δ (r - r^{'}) \sum_{l = 0}^{\infty} \sum_{m = - l}^{l} Y_{l m}^{*} (θ^{'}, ϕ^{'}) Y_{l m} (θ, ϕ) .

Substituting Eqs. (29) and (30) into Eq. (28) and expanding $\nabla^{2}$ in spherical coordinates lead to

Eq. (31)

\frac{1}{r} \frac{\partial^{2}}{\partial r^{2}} [r G (\vec{χ}, {\vec{χ}}^{'})] + \frac{1}{r^{2} \sin θ} \frac{\partial}{\partial θ} (\sin θ \frac{\partial G (\vec{χ}, {\vec{χ}}^{'})}{\partial θ}) + \frac{1}{r^{2} \sin^{2} θ} \frac{\partial^{2} G (\vec{χ}, {\vec{χ}}^{'})}{\partial ϕ^{2}} - μ_{eff}^{2} G (\vec{χ}, {\vec{χ}}^{'}) = - \frac{1}{r^{2}} δ (r - r^{'}) \sum_{l = 0}^{\infty} \sum_{m = - l}^{l} Y_{l m}^{*} (θ^{'}, ϕ^{'}) Y_{l m} (θ, ϕ) .

The Green’s function is expanded to a form similar to the right-hand side of Eq. (30) as [Ref. 38, page 427, Eq. (9.95)]

Eq. (32)

G (\overset{⇀}{χ}, {\vec{χ}}^{'}) = \sum_{l = 0}^{\infty} \sum_{m = - l}^{l} g_{l} (r, r^{'}) Y_{l m}^{*} (θ^{'}, ϕ^{'}) Y_{l m} (θ, ϕ),

where

g_{l} (r, r^{'})

is the radial Green’s function. Substituting Eq. (32) into Eq. (31) leads to

Eq. (33)

{\frac{d^{2}}{d r^{2}} + \frac{2}{r} \frac{d}{d r} - [μ_{eff}^{2} + \frac{l (l + 1)}{r^{2}}]} g_{l} = - \frac{1}{r^{2}} δ (r - r^{'}),

which is the differential equation for the modified spherical Bessel function. The solution that satisfies the finiteness at the origin and the infinity and assures the correct discontinuity in slope of the Dirac delta function is (approach similar to that shown in Ref. 38 on page 427)

Eq. (34)

g_{l} = (μ_{eff}) i_{l} (μ_{eff} r_{<}) k_{l} (μ_{eff} r_{>}),

where

i_{l}

and

k_{l}

are respectively the

l

’th order modified spherical Bessel function of the first and the second kinds, and

r_{<}

and

r_{>}

are respectively the smaller and greater radial coordinates of the source and detector. Thus, the Green’s function is

Eq. (35)

G (\overset{⇀}{χ}, {\vec{χ}}^{'}) = (μ_{eff}) \sum_{l = 0}^{\infty} i_{l} (μ_{eff} r_{<}) k_{l} (μ_{eff} r_{>}) \sum_{m = - l}^{l} Y_{l m}^{*} (θ^{'}, ϕ^{'}) Y_{l m} (θ, ϕ) .

Convolving the Green’s function with the source term in Eq. (1), assuming a point source, renders the solution, expressed using the modified spherical Bessel functions, to the steady-state photon diffusion in an infinite homogeneous medium as

Eq. (36)

Ψ (\overset{⇀}{χ}, {\vec{χ}}^{'}) = \frac{S}{D} (μ_{eff}) \sum_{l = 0}^{\infty} i_{l} (μ_{eff} r_{<}) k_{l} (μ_{eff} r_{>}) \sum_{m = - l}^{l} Y_{l m}^{*} (θ^{'}, ϕ^{'}) Y_{l m} (θ, ϕ) .

Appendix B

The numerical identity between Eq. (3) and the well-known form of Eq. (2) is worth being justified. With the use of the following definitions and identities:

Eq. (37)

i_{l} (μ_{eff} r_{<}) = \sqrt{\frac{π}{2 μ_{eff} r_{<}}} I_{l + 1 / 2} (μ_{eff} r_{<}),

Eq. (38)

k_{l} (μ_{eff} r_{<}) = \sqrt{\frac{2}{π μ_{eff} r_{>}}} K_{l + 1 / 2} (μ_{eff} r_{>}),

Eq. (39)

\cos γ = \cos θ \cos θ^{'} + \sin θ \sin θ^{'} \cos (ϕ - ϕ^{'}),

Eq. (40)

P_{l} (\cos γ) = \frac{4 π}{2 l + 1} \sum_{m = - l}^{l} Y_{l m}^{*} (θ^{'}, ϕ^{'}) Y_{l m} (θ, ϕ),

where

I_{l}

and

K_{l}

are respectively the

l

’th order modified Bessel function of the first and the second kinds,

γ

is the angle formed between the position vector of the source and that of the detector, and

P_{l}

is the

l

’th order Legendre polynomial, Eq. (3) converts into a form that is more convenient for numerical evaluation:

Eq. (41)

Ψ (\overset{⇀}{χ}, {\vec{χ}}^{'}) = \frac{S}{4 π D} \sum_{l = 0}^{\infty} \sqrt{\frac{1}{r_{<} \cdot r_{>}}} I_{l + 1 / 2} (μ_{eff} r_{<}) K_{l + 1 / 2} (μ_{eff} r_{>}) (2 l + 1) P_{l} (\cos γ) .

The value determined by Eq. (41) is compared to that by Eq. (2) for a number of configurations of the source–detector pair, at common optical properties of $μ_{a} = 0.01 {mm}^{- 1}$ and $μ_{s}^{'} = 1.0 {mm}^{- 1}$ and different upper limits of $l$ for the summation. The results for the case of source–detector pair aligning with the origin, that is, the source and the detector are respectively at $(R_{0}, θ, ϕ)$ and $(R_{0} + d, θ, ϕ)$ , are presented in Fig. 9. The results are presented in terms of $\ln (d * Ψ)$ versus $d$ , as it leads to a convenient linear relationship as evident from Eq. (2). At a radial coordinate of $R_{0} = 50 mm$ , it takes the summation of 30 terms with Eq. (41) to reach an error less than 0.1% at source–detector distance greater than 10 mm, when compared to the value given by Eq. (2). At a larger radial coordinate of $R_{0} = 200 mm$ , it takes summation of 100 terms with Eq. (41) to reach an error less than 0.1% at source–detector distance greater than 10 mm, when compared to the value given by Eq. (2). In the cases when the source–detector pair is not aligned with the origin (i.e., the source and the detector are separated in either azimuthal or elevational directions), the values given by Eq. (41) as a function of the source–detector distance oscillate with respect to the expected correct values (due to the limit of the numerical arithmetic in the computational software) with respect to the accurate value at the summation up to 200 terms. We note that should Eq. (41) be numerically nonidentical to Eq. (2), the numerical evaluation based on Eq. (41) for any source–detector configuration would have systematically deviated from that based on Eq. (2). With the comparison shown in Fig. 9, an upper limit of $l = 250$ is implemented in the numerical evaluations of the solutions involving Eq. (3) or Eq. (41) when needed in the rest of the study.

Fig. 9

Numerical comparison of the values given by Eqs. (41) and (2), at a common set of optical properties of $μ_{a} = 0.01 {mm}^{- 1}$ and $μ_{s}^{'} = 1.0 {mm}^{- 1}$ , and different upper limit of $l$ for the summation. The source and the detector are respectively at $(R_{0}, θ, ϕ)$ and $(R_{0} + d, θ, ϕ)$ . The results are presented in terms of the linear relationship expected for $\ln (d * Ψ)$ versus $d$ . At the radial coordinates of $R_{0} = 50 mm$ and $R_{0} = 200 mm$ , it takes the summation of 30 terms and 100 terms, respectively, with Eq. (41) to reach an error less than 0.1% at source–detector distance greater than 10 mm when compared to the value given by Eq. (2).

Acknowledgments

This work was supported in part by the National Institutes of Health (Grant No. 1R21NS067278) and a research initiative fund provided by the University of Oklahoma College of Medicine and the Oklahoma City VA Medical Center.

References

1.

F. F. Jöbsis, “Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science, 198 (4323), 1264 –1267 (1977). http://dx.doi.org/10.1126/science.929199 SCIEAS 0036-8075 Google Scholar

2.

B. Molavi et al., “Analyzing the resting state functional connectivity in the human language system using near infrared spectroscopy,” Front. Hum. Neurosci., 7 921 (2014). http://dx.doi.org/10.3389/fnhum.2013.00921 Google Scholar

3.

C. G. Favilla et al., “Optical bedside monitoring of cerebral blood flow in acute ischemic stroke patients during head-of-bed manipulation,” Stroke, 45 (5), 1269 –1274 (2014). http://dx.doi.org/10.1161/STROKEAHA.113.004116 SJCCA7 0039-2499 Google Scholar

4.

A. E. Cerussi et al., “Diffuse optical spectroscopic imaging correlates with final pathological response in breast cancer neoadjuvant chemotherapy,” Philos. Trans. A Math Phys. Eng. Sci., 369 (1955), 4512 –4530 (2011). http://dx.doi.org/10.1098/rsta.2011.0279 Google Scholar

5.

C. H. Schmitz et al., “Instrumentation and calibration protocol for imaging dynamic features in dense-scattering media by optical tomography,” Appl. Opt., 39 6466 –6486 (2000). http://dx.doi.org/10.1364/AO.39.006466 APOPAI 0003-6935 Google Scholar

6.

D. T. Delpy et al., “Estimation of optical pathlength through tissue from direct time of flight measurement,” Phys. Med. Biol., 33 (12), 1433 –1442 (1988). http://dx.doi.org/10.1088/0031-9155/33/12/008 Google Scholar

7.

A. Duncan et al., “Optical pathlength measurements on adult head, calf and forearm and the head of the newborn infant using phase resolved optical spectroscopy,” Phys. Med. Biol., 40 (2), 295 –304 (1995). http://dx.doi.org/10.1088/0031-9155/40/2/007 Google Scholar

8.

F. Scholkmann and M. Wolf, “General equation for the differential pathlength factor of the frontal human head depending on wavelength and age,” J. Biomed. Opt., 18 (10), 105004 (2013). http://dx.doi.org/10.1117/1.JBO.18.10.105004 Google Scholar

9.

S. Fantini et al., “Non-invasive optical mapping of the piglet brain in realtime,” Opt. Express, 4 (8), 308 –314 (1999). http://dx.doi.org/10.1364/OE.4.000308 OPEXFF 1094-4087 Google Scholar

10.

D. A. Boas et al., “The accuracy of near infrared spectroscopy and imaging during focal changes in cerebral hemodynamics,” Neuroimage, 13 (1), 76 –90 (2001). http://dx.doi.org/10.1006/nimg.2000.0674 NEIMEF 1053-8119 Google Scholar

11.

Y. Hoshi, “Towards the next generation of near-infrared spectroscopy,” Philos. Trans. A Math Phys. Eng. Sci., 369 (1955), 4425 –4439 (2011). http://dx.doi.org/10.1098/rsta.2011.0262 Google Scholar

12.

A. Duncan et al., “Measurement of cranial optical path length as a function of age using phase resolved near infrared spectroscopy,” Pediatr. Res., 39 (5), 889 –894 (1996). http://dx.doi.org/10.1203/00006450-199605000-00025 Google Scholar

13.

T. Talukdar, J. H. Moore and S. G. Diamond, “Continuous correction of differential path length factor in near-infrared spectroscopy,” J. Biomed. Opt., 18 (5), 056001 (2013). http://dx.doi.org/10.1117/1.JBO.18.5.056001 JBOPFO 1083-3668 Google Scholar

14.

L. Gagnon et al., “Further improvement in reducing superficial contamination in NIRS using double short separation measurements,” Neuroimage, 85 (pt 1), 127 –135 (2014). http://dx.doi.org/10.1016/j.neuroimage.2013.01.073 NEIMEF 1053-8119 Google Scholar

15.

S. R. Arridge et al., “A finite element approach for modeling photon transport in tissue,” Med. Phys., 20 (2 Pt 1), 299 –309 (1993). http://dx.doi.org/10.1118/1.597069 Google Scholar

16.

M. Schweiger, S. R. Arridge and D. T. Delpy, “Application of the finite-element method for the forward and inverse models in optical tomography,” J. Math. Imaging Vis., 3 (3), 263 –283 (1993). http://dx.doi.org/10.1007/BF01248356 Google Scholar

17.

X. Deulin and J. P. L'Huillier, “Finite element approach to photon propagation modeling in semi-infinite homogeneous and multilayered tissue structures,” Eur. Phys. J. Appl. Phys., 33 133 –146 (2006). http://dx.doi.org/10.1051/epjap:2006010 EPAPFV 1286-0042 Google Scholar

18.

S. R. Arridge, M. Cope and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol., 37 (7), 1531 –1560 (1992). http://dx.doi.org/10.1088/0031-9155/37/7/005 Google Scholar

19.

A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. I. Homogeneous case,” Opt. Express, 18 9456 –9473 (2010). http://dx.doi.org/10.1364/OE.18.009456 OPEXFF 1094-4087 Google Scholar

20.

A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. II. Layered case,” Opt. Express, 18 9266 –9279 (2010). http://dx.doi.org/10.1364/OE.18.009266 OPEXFF 1094-4087 Google Scholar

21.

A. Da Silva et al., “Optical calibration protocol for an x-ray and optical multimodality tomography system dedicated to small-animal examination,” Appl. Opt., 48 D151 –D162 (2009). http://dx.doi.org/10.1364/AO.48.00D151 APOPAI 0003-6935 Google Scholar

22.

A. Sassaroli et al., “Performance of fitting procedures in curved geometry for retrieval of the optical properties of tissue from time-resolved measurements,” Appl. Opt., 40 185 –197 (2001). http://dx.doi.org/10.1364/AO.40.000185 APOPAI 0003-6935 Google Scholar

23.

R. C. Haskell et al., “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A, 11 2727 –2741 (1994). http://dx.doi.org/10.1364/JOSAA.11.002727 JOAOD6 0740-3232 Google Scholar

24.

D. Contini, F. Martelli and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt., 36 4587 –4599 (1997). http://dx.doi.org/10.1364/AO.36.004587 APOPAI 0003-6935 Google Scholar

25.

A. Zhang et al., “Photon diffusion in a homogeneous medium bounded externally or internally by an infinitely long circular cylindrical applicator. I. Steady-state theory,” J. Opt. Soc. Am. A, 27 648 –662 (2010). http://dx.doi.org/10.1364/JOSAA.27.000648 JOAOD6 0740-3232 Google Scholar

26.

A. Zhang and D. Piao, “Photon diffusion in a homogeneous medium bounded externally or internally by an infinitely long circular cylindrical applicator. IV. Frequency-domain analysis,” J. Opt. Soc. Am. A, 29 1445 –1458 (2012). http://dx.doi.org/10.1364/JOSAA.29.001445 JOAOD6 0740-3232 Google Scholar

27.

D. Piao, “Photon diffusion in a homogeneous medium bounded externally or internally by an infinitely long circular cylindrical applicator. VI. Time-domain analysis,” J. Opt. Soc. Am. A, 31 2232 –2243 (2014). http://dx.doi.org/10.1364/JOSAA.31.002232 JOAOD6 0740-3232 Google Scholar

28.

S. A. Walker, D. A. Boas and E. Gratton, “Photon density waves scattered from cylindrical inhomogeneities: theory and experiments,” Appl. Opt., 37 (10), 1935 –1944 (1998). http://dx.doi.org/10.1364/AO.37.001935 Google Scholar

29.

K. Furutsu and Y. Yamada, “Diffusion approximation for a dissipative random medium and the applications,” Phys. Rev. E, 50 3634 (1994). http://dx.doi.org/10.1103/PhysRevE.50.3634 Google Scholar

30.

R. Aronson and N. Corngold, “Photon diffusion coefficient in an absorbing medium,” J. Optical Soc. Am. A, 16 1066 –1071 (1999). http://dx.doi.org/10.1364/JOSAA.16.001066 Google Scholar

31.

R. Aronson, “Radiative transfer implies a modified reciprocity relation,” J. Opt. Soc. Am. A, 14 486 –490 (1997). http://dx.doi.org/10.1364/JOSAA.14.000486 Google Scholar

32.

R. Elaloufi, R. Carminati and J.-J. Greffet, “Definition of the diffusion coefficient in scattering and absorbing media,” J. Opt. Soc. Am. A, 20 678 (2003). http://dx.doi.org/10.1364/JOSAA.20.000678 JOAOD6 0740-3232 Google Scholar

33.

R. Pierrat, J. J. Greffet and R. Carminati, “Photon diffusion coefficient in scattering and absorbing media,” J. Opt. Soc. Am. A, 23 1106 (2006). http://dx.doi.org/10.1364/JOSAA.23.001106 JOAOD6 0740-3232 Google Scholar

34.

S. Fantini, M. A. Franceschini and E. Gratton, “Semiinfinite-geometry boundary problem for light migration in highly scattering media: a frequency-domain study in the diffusion approximation,” J. Opt. Soc. Am. B, 11 2128 –2138 (1994). http://dx.doi.org/10.1364/JOSAB.11.002128 JOBPDE 0740-3224 Google Scholar

35.

K. S. Fine and C. F. Driscoll, “The finite length diocotron mode,” Phys. Plasmas, 5 601 –607 (1998). http://dx.doi.org/10.1063/1.872752 PHPAEN 1070-664X Google Scholar

36.

S. T. Cui, “Electrostatic potential in cylindrical dielectric media using the image charge method,” Mol. Phys., 104 2993 –3001 (2006). http://dx.doi.org/10.1080/00268970600926647 MOPHAM 0026-8976 Google Scholar

37.

A. Zhang et al., “Photon diffusion in a homogeneous medium bounded externally or internally by an infinitely long circular cylindrical applicator. II. Quantitative examinations of the steady-state theory,” J. Opt. Soc. Am. A., 28 66 –75 2011). http://dx.doi.org/10.1203/01.PDR.0000175638.98041.0E Google Scholar

38.

J. D. Jackson, Classical Electrodynamics, 3rd edJohn Wiley & Sons, Inc., Hoboken, New Jersey (1998). Google Scholar

39.

H. Dehghani et al., “Near infrared optical tomography using NIRFAST: algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng., 25 711 –732 (2009). http://dx.doi.org/10.1002/cnm.1162 CANMER 0748-8025 Google Scholar

40.

M. Jermyn et al., “Fast segmentation and high-quality three-dimensional volume mesh creation from medical images for diffuse optical tomography,” J. Biomed. Opt., 18 (8), 086007 (2013). http://dx.doi.org/10.1117/1.JBO.18.8.086007 JBOPFO 1083-3668 Google Scholar

41.

S. Ijichi et al., “Developmental changes of optical properties in neonates determined by near-infrared time-resolved spectroscopy,” Pediatr. Res., 58 (3), 568 –573 2005). http://dx.doi.org/10.1203/01.PDR.0000175638.98041.0E Google Scholar

42.

P. van der Zee et al., “The effect of optode positioning on optical pathlength in near infrared spectroscopy of brain,” Adv. Exp. Med. Biol., 277 79 –84 (1990). http://dx.doi.org/10.1007/978-1-4684-8181-5_11 Google Scholar

43.

P. van der Zee et al., “Experimentally measured optical pathlengths for the adult head, calf and forearm and the head of the newborn infant as a function of inter optode spacing,” Adv. Exp. Med. Biol., 316 143 –153 (1992). http://dx.doi.org/10.1007/978-1-4615-3404-4_17 Google Scholar

44.

M. Hiraoka et al., “A Monte Carlo investigation of optical pathlength in inhomogeneous tissue and its application to near-infrared spectroscopy,” Phys. Med. Biol., 38 (12), 1859 –1876 (1993). http://dx.doi.org/10.1088/0031-9155/38/12/011 Google Scholar

45.

J. L. Robertson et al., “Effect of blood in the cerebrospinal fluid on the accuracy of cerebral oxygenation measured by near infrared spectroscopy,” Adv. Exp. Med. Biol., 812 233 –240 (2014). http://dx.doi.org/10.1007/978-1-4939-0620-8_31 Google Scholar

Biography

Biographies for the authors are not available.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Daqing Piao, Randall L. Barbour, Harry L. Graber, and Daniel C. Lee M.D. "On the geometry dependence of differential pathlength factor for near-infrared spectroscopy. I. Steady-state with homogeneous medium," Journal of Biomedical Optics 20(10), 105005 (14 October 2015). https://doi.org/10.1117/1.JBO.20.10.105005

Published: 14 October 2015

Access the abstract

JOURNAL ARTICLE
14 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 29 scholarly publications.

Explore citations on Lens.org

KEYWORDS

Sensors

Spherical lenses

Optical spheres

Diffusion

Finite element methods

Near infrared spectroscopy

Scattering

1.

Introduction

2.

Geometries Considered and Analytic Preparations for Steady-State Photon Diffusion in Infinite and Semi-Infinite Geometries

Fig. 1

2.1.

Infinite Geometry and Semi-Infinite Geometry

2.2.

Infinite-Length Cylindrical Geometry

2.3.

Spherical Geometry

2.4.

Steady-State Photon Diffusion in an Infinite Homogeneous Domain

Eq. (1)

Eq. (2)

Eq. (3)

2.5.

Steady-State Photon Diffusion in a Homogeneous Semi-Infinite Medium

Eq. (4)

Eq. (5)

Eq. (6)

3.

Steady-State Photon Diffusion in a Homogeneous Medium Bounded in a Spherical Geometry

3.1.

Solution of Steady-State Photon Diffusion Associated with a Homogeneous Medium of Spherical Geometry

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

Eq. (11)

Eq. (12)

Fig. 2

3.2.

Large-Diameter Approximation of the Solution for a Spherical Domain

Eq. (13)

Eq. (14)

Eq. (15)

Eq. (16)

Eq. (17)

Eq. (18)

Eq. (19)

Eq. (20)

Eq. (21)

Table 1

4.

Differential Pathlength Factors Associated with Homogeneous Semi-Infinite, Infinite-Length Cylinder When Evaluated along the Azimuthal Direction and Spherical Geometries

Eq. (22)

Eq. (23)

Eq. (24)

Eq. (25)

Eq. (26)

Eq. (27)

5.

Numerical Evaluation of the Differential Pathlength Factors Associated with Homogeneous Semi-Infinite, Infinite Cylinder Evaluated along the Azimuthal Direction and Spherical Geometries

Fig. 3

Fig. 4

6.

Discussions

Fig. 5

Fig. 6

Fig. 7

Fig. 8

7.

Conclusion

Appendices

Appendix A

Eq. (28)

Eq. (29)

Eq. (30)

Eq. (31)

Eq. (32)

Eq. (33)

Eq. (34)

Eq. (35)

Eq. (36)

Appendix B

Eq. (37)

Eq. (38)

Eq. (39)

Eq. (40)