Research Papers: General

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

[+] Author Affiliations
Daqing Piao

Oklahoma State University, School of Electrical and Computer Engineering, Stillwater, Oklahoma 74078, United States

Randall L. Barbour

SUNY Downstate Medical Center, Department of Pathology, Brooklyn, New York 11203, United States

NIRx Medical Technologies LLC, Glen Head, New York 11545, United States

Harry L. Graber

NIRx Medical Technologies LLC, Glen Head, New York 11545, United States

Daniel C. Lee

University of Oklahoma College of Medicine, Department of Surgery, Oklahoma City, Oklahoma 73104, United States

J. Biomed. Opt. 20(10), 105005 (Oct 14, 2015). doi:10.1117/1.JBO.20.10.105005
History: Received June 23, 2015; Accepted August 28, 2015
Text Size: A A A

Open Access Open Access

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.01mm1 and μs=1.0mm1 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.

Figures in this Article

Near-infrared spectroscopy (NIRS)1 is becoming an increasingly important modality for a number of investigational and clinical needs, including noninvasive functional study of neurophysiological connectivity,2 bedside monitoring of cerebral hemodynamics,3 and evaluating responses to cancer treatment.4 Most NIRS systems for long duration monitoring operate under the continuous-wave (CW) or steady-state condition.5 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.6 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±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.7 at source–detector spacings >4cm yielded DPFs of 6.26 (±14.1%) for an adult head, 5.51 (±18%) for an adult leg, 4.99 (±9%) for the head of a newborn infant, and 4.16 (±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.6 suggested a wavelength-, medium-, and geometry-dependent constant8 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.9 and Boas et al.10 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 head11 (e.g., prefrontal or parietal), and even different ages12 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.13 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.14 Whereas estimates of DPFs from complex geometries are readily available with the use of discrete numerical methods [e.g., finite-element method (FEM)1517], 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.1822 A particularly complete set of solutions in time domain and frequency domain in a number of geometries were presented by Arridge et al.18 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.22 developed the corresponding solutions of time-resolved DE using the more widely used extrapolated boundary condition (EBC)21,23,24 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,25 frequency domain26 and time domain,27 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. 252627 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. 252627 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.

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.2933 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 χ=(x,y,z) to a detector at χ=(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 χ=(ρ,ϕ,z) to a detector at χ=(ρ,ϕ,z). The sphere domain will be modeled in the spherical coordinates, where we are concerned about photon pathlength from a source at χ=(r,θ,ϕ) to a detector at χ=(r,θ,ϕ).

Graphic Jump Location
Fig. 1
F1 :

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

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 χ=(0,0,0), and the detector at χ=(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 χ=(0,0,0), and the detector at χ=(0,d,0). The physical source at χ=(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 Ra=1/μs, into the medium. The “real” isotropic source χreal thus has the coordinates of (0,0,Ra). The effect of the medium-applicator interface on photon diffusion as modeled by the EBC21,23,24 sets zero the photon fluence rate on an imaginary boundary located at a distance of Rb=2AD away from the physical boundary, where A is a coefficient23 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 χimag', which has the coordinates of (0,0,Ra+2Rb). The distances from the detector to the “real” isotropic source χreal and the “image” source χimag are denoted by lr and li, respectively. The notations of lr=|χχreal| and li=|χχ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=|χχ| 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.

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 R0, then a detector χ on the medium-applicator interface locates at (R0,ϕ,z), and a directional source χ on the medium-applicator interface locates at (R0,ϕ,z). The geometric symmetry requires that the directional source χ be modeled by a “real” isotropic source χreal positioned along the radial direction of χ and inwardly at a radial distance of Ra=1/μs, at (R0Ra,ϕ,z). The extrapolated zero-boundary is concentric with and at a radial distance of Rb=2AD outward from the physical boundary. Apparently, as the radius R0 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.25

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.

Spherical Geometry

A spherical geometry is illustrated in Fig. 1(d). The radius of the spherical applicator is R0, then a detector χ on the medium-applicator interface locates at (R0,θ,ϕ), and a directional source χ on the medium-applicator interface locates at (R0,θ,ϕ). The geometric symmetry requires that the directional source χ be modeled by a “real” isotropic source rreal positioned along the radial direction of χ and inwardly at a radial distance of Ra=1/μs, [i.e., at (R0Ra,θ,ϕ)]. The extrapolated zero boundary is concentric with and at a radial distance of Rb=2AD outward from the physical boundary. Apparently, as the radius R0 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.

Steady-State Photon Diffusion in an Infinite Homogeneous Domain

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

2Ψ(χ)μaDΨ(χ)=S(χ)D,(1)
where Ψ is the photon fluence rate and S is the source term. For an infinite medium and a point source at χ with an intensity S as represented by S·δ(χ), Eq. (1) has the well-known closed-form solution of Ref. 25Display Formula
Ψ(χ,χ)=S4πD1|χχ|eμeff|χχ|,(2)
where μeff=μa/D=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 1): Display Formula
Ψ(χ,χ)=SD(μeff)l=0il(μeffr<)kl(μeffr>)m=llYlm*(θ,ϕ)Ylm(θ,ϕ),(3)
where il and kl 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 Ylm is the spherical harmonics function. The numerical identity between Eq. (3) and the well-known form of Eq. (2) is justified in Appendix 2.

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

We follow the notations summarized by Fantini et al.34 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 Display Formula

Ψ=ΨrealΨimag=S4πDlreμefflrS4πDlieμeffli,(4)
where Display Formula
lr=d1+(Ra/d)2,Ra=1/μs,(5)
Display Formula
li=d1+[(Ra+2Rb)/d]2,Rb=2AD.(6)

The principles to find the solutions to steady-state photon diffusion associated with a homogeneous infinite-length cylinder25 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.

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 R0, the physical source is located at (R0,θ,ϕ) and the detector is located at (R0,θ,ϕ), both on the physical boundary. The equivalent “real” isotropic source locates at (R0Ra,θ,ϕ), and the extrapolated boundary locates at a radial position of ρr>=R0+2AD 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<=R0Ra and the detector locates at r>=R0+Rb, is Display Formula

Ψreal|extr=1D(μeff)l=0S·il[μeff(R0Ra)]kl[μeff(R0+Rb)]m=llYlm*(θ,ϕ)Ylm(θ,ϕ),(7)
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<=R0+Rb, is Display Formula
Ψimag|extr=1D(μeff)l=0Sl*il[μeff(R0+Rb)]kl[μeffr>]m=llYlm*(θ,ϕ)Ylm(θ,ϕ),(8)
where the Sl* terms associated with each “l”’th order are not the same. Based on the essence of “image-source,”35,36 the two unknown terms Sl* and r> associated with the l’th order “image” source (the kl component) can be expressed by a single unknown term Sl associated with the same order “real” source (the il component), that is Display Formula
Sl*kl(μeffr>)=Slil[μeff(R0Ra)].(9)

Applying Eq. (9) to the EBC of Ψreal|extrΨimag|extr=0, we have Display Formula

Sl=Skl[μeff(R0+Rb)]il[μeff(R0+Rb)]l=0,1,2,(10)

Now for the “real” source but evaluating at the physical boundary, the source still locates at r<=R0Ra but the detector locates at r>=R0. For the “image” source also evaluating at the physical boundary, the detector locates at r<=R0 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 Display Formula

Ψ=Ψreal|physΨimag|phys=SD(μeff)l=0il[μeff(R0Ra)]kl[μeffR0]m=llYlm*(θ,ϕ)Ylm(θ,ϕ)SD(μeff)l=0il[μeff(R0Ra)]il[μeffR0]kl[μeff(R0+Rb)]il[μeff(R0+Rb)]m=llYlm*(θ,ϕ)Ylm(θ,ϕ)=SD(μeff)l=0il[μeff(R0Ra)]kl[μeffR0]m=llYlm*(θ,ϕ)Ylm(θ,ϕ)1il[μeffR0]kl[μeffR0]kl[μeff(R0+Rb)]il[μeff(R0+Rb)].(11)

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

Ψ=SD1(R0Ra)R0l=0Il+1/2[μeff(R0Ra)]Kl+1/2(μeffR0)2l+14πPl(cosγ)1Il+1/2[μeffR0]Kl+1/2[μeffR0]Kl+1/2[μeff(R0+Rb)]Il+1/2[μeff(R0+Rb)],(12)
where Il and Kl are respectively the l’th order modified Bessel function of the first and the second kinds, Pl 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 1).

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.01mm1 and μs=0.25mm1, the middle traces μa=0.01mm1 and μs=0.50mm1, and the bottom traces μa=0.01mm1 and μs=1.0mm1. 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(d2*Ψ) 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.0mm1 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 study37 dedicated to experimental examination of the analytical results of photon fluence rate specific to infinite-length cylinder,25 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.

Graphic Jump Location
Fig. 2
F2 :

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.01mm1 with μs=0.25mm1 (top traces), μs=0.5mm1 (middle traces), and μs=1.0mm1 (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.0mm1, but nevertheless the numerical computations marked by the solid lines clearly follow the traces indicated by the FEM results.

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 (R0Ra,θ,ϕ), but the image of the “real” isotropic source with respect to this “hypothetical” semi-infinite boundary will be at (R0+Ra+2Rb,θ,ϕ). 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 Display Formula

Ψimagsemi|phys=SD(μeff)l=0il[μeffR0]kl[μeff(R0+Ra+2Rb)]m=llYlm*(θ,ϕ)Ylm(θ,ϕ).(13)

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 Display Formula

Ψimag|phys=SD(μeff)l=0il[μeffR0]kl[μeff(R0+Ra+2Rb)]·ηl·m=llYlm*(θ,ϕ)Ylm(θ,ϕ),(14)
where Display Formula
ηl=kl[μeff(R0+Rb)]il[μeff(R0+Rb)]il[μeff(R0Ra)]kl[μeff(R0+Ra+2Rb)].(15)

Note that R0Ra, Rb and if the argument of the modified Bessel functions is much greater than 1 (i.e., μeffR010 that corresponds to R057.8mm for μa=0.01mm1 and μs=1.0mm1) the modified spherical Bessel functions in Eq. (15) can be simplified by their asymptotic expressions,38 then Eq. (15) becomes Display Formula

ηl=R0+Ra+2RbR0Ra.(16)

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 Display Formula

Ψimag|phys=Ψimagsemi|phys·R0+Ra+2RbR0Ra.(17)

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

Ψ=Ψreal|physΨimag|phys=Ψreal|physΨimagsemi|phys·R0+Ra+2RbR0Ra.(18)

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 (R0+Ra+2Rb)/(R0Ra)1, and the Ψimagsemi|phys of Eq. (18) will become the Ψimag in Eq. (4) because the detector point locating at (R0,θ,ϕ) 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 Display Formula

Ψ=S4πDeμefflrlrS4πDeμeffliliR0+Ra+2RbR0Ra.(19)

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): Display Formula

lr=d1+Δr;whereΔr=Ra2d2RaR0,(20)
Display Formula
li=d1+Δi;whereΔi=(Ra+2Rb)2d2+Ra+2RbR0.(21)

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 μeffR0=8.66 will indicate that the approximate formula for the infinite cylinder and the sphere is generally valid when the radius satisfies μeffR010, that is, the radius (diameter) is more than 10 (20) times of the effective attenuation pathlength 1/μeff.

Table Grahic Jump Location
Table 1Analytical (large-radius approximation) representations of the photon diffusion.

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·DPF·d], where “Att” is the ratio of measured light over the incident light. We use the definition of DPF provided by Fantini et al.9 or Arridge et al.18 that becomes DPF={[ln(S/Ψ)]/μ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] Display Formula
    DPFinf=1dμa[ln(SΨinf)]=μa(μeff)=3μs2μa.(22)
  • Semi-infinite geometry [identical to Eq. (9) of Ref. 9] Display Formula
    DPFsemi=1dμa[ln(SΨsemi)]=1dexp(μefflr)exp(μeffli)exp(μefflr)lrexp(μeffli)liμa(μeff)3μs2μad3μsμad3μsμa+1.(23)
  • Infinite-length cylindrical geometry, evaluated along the azimuthal direction Display Formula
    DPFazi=1dμa[ln(SΨazi)]=1dexp(μefflr)[exp(μeffli)]R0+Ra+2RbR0Raexp(μefflr)lrexp(μeffli)liR0+Ra+2RbR0Raμa(μeff)=3μs2μad3μsμa[ηaziΔiΔr][ηazi1](d3μsμa+1)[ηaziΔiΔr][ηazi1],(24)
    Display Formula
    whereηazi=R0+Ra+2RbR0Ra.(25)
  • Spherical geometry Display Formula
    DPFsph=1dμa[ln(SΨsph)]=1dexp(μefflr)[exp(μeffli)]R0+Ra+2RbR0Raexp(μefflr)lrexp(μeffli)liR0+Ra+2RbR0Raμa(μeff)=3μs2μad3μsμa[ηsphΔiΔr][ηsph1](d3μsμa+1)[ηsphΔiΔr][ηsph1],(26)
    Display Formula
    whereηsph=R0+Ra+2RbR0Ra.(27)

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

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),39,40 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 R0=25mm and R0=50mm 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 100mm×100mm 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).

Graphic Jump Location
Fig. 3
F3 :

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.01mm1, μs=1.0mm1. (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.011mm1 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.

Graphic Jump Location
Fig. 4
F4 :

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.01mm1, μs=1.0mm1, refractive index of the medium is 1.33 (the same for the following figures), and R0=25mm 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.01mm1, μs=1.0mm1, R0=25mm, and R0=50mm 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 R0, 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 μeffR0=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.

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.12 At 832 nm, the DPFs measured by using phase-resolved NIRS from the head of neonates were 4.67±0.65, while the DPFs of adults were 5.86±0.98,7 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.01mm1, with the μs changing from 0.5mm1 to 1.5mm1. 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.0mm1, with the μa changing from 0.01mm1 to 0.1mm1. 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.5mm1 to 1.5mm1. 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.

Graphic Jump Location
Fig. 5
F5 :

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.

Graphic Jump Location
Fig. 6
F6 :

The dependence of DPF on the absorption coefficient from 0.01 to 0.1mm1 at a fixed reduced scattering coefficient of 1.0mm1. 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.

Graphic Jump Location
Fig. 7
F7 :

The dependence of DPF on the reduced scattering coefficient from 0.5 to 1.5mm1 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.8 Ijichi et al.41 measured DPFs over a considerably narrow postconceptional range of age from 30 to 42 weeks by phase-resolved measurements. The DPF values 4.31±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.43mm1 to 0.79mm1 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.41 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.01mm1, μs=1.0mm1. 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.41 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. 678, 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.

Graphic Jump Location
Fig. 8
F8 :

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.01mm1 and μs=1.0mm1.

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 calculations42 and the analytical solution of the diffusion equation incorporating the zero-boundary condition.18 However, there are experimental studies of the adult head43 suggesting that the DPF increases at short source–detector separation (2cm). 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)44 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.45 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.

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 μeffR010, 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.01mm1 and μs=1.0mm1 is on average 5% different compared to the value obtained from FEM.

Considering a source at x of (r,θ,ϕ) and a detector at x of (r,θ,ϕ) in spherical coordinates, the equation for the Green’s function of Eq. (1) is Display Formula

2G(χ,χ)μeff2G(χ,χ)=δ(χχ).(28)

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

δ(χχ)=1r2δ(rr)δ(cosθcosθ)δ(ϕϕ),(29)
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)] Display Formula
δ(χχ)=1r2δ(rr)l=0m=llYlm*(θ,ϕ)Ylm(θ,ϕ).(30)

Substituting Eqs. (29) and (30) into Eq. (28) and expanding 2 in spherical coordinates lead to Display Formula

1r2r2[rG(χ,χ)]+1r2sinθθ(sinθG(χ,χ)θ)+1r2sin2θ2G(χ,χ)ϕ2μeff2G(χ,χ)=1r2δ(rr)l=0m=llYlm*(θ,ϕ)Ylm(θ,ϕ).(31)

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)] Display Formula

G(χ,χ)=l=0m=llgl(r,r)Ylm*(θ,ϕ)Ylm(θ,ϕ),(32)
where gl(r,r) is the radial Green’s function. Substituting Eq. (32) into Eq. (31) leads to Display Formula
{d2dr2+2rddr[μeff2+l(l+1)r2]}gl=1r2δ(rr),(33)
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) Display Formula
gl=(μeff)il(μeffr<)kl(μeffr>),(34)
where il and kl 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 Display Formula
G(χ,χ)=(μeff)l=0il(μeffr<)kl(μeffr>)m=llYlm*(θ,ϕ)Ylm(θ,ϕ).(35)

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 Display Formula

Ψ(χ,χ)=SD(μeff)l=0il(μeffr<)kl(μeffr>)m=llYlm*(θ,ϕ)Ylm(θ,ϕ).(36)

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: Display Formula

il(μeffr<)=π2μeffr<Il+1/2(μeffr<),(37)
Display Formula
kl(μeffr<)=2πμeffr>Kl+1/2(μeffr>),(38)
Display Formula
cosγ=cosθcosθ+sinθsinθcos(ϕϕ),(39)
Display Formula
Pl(cosγ)=4π2l+1m=llYlm*(θ,ϕ)Ylm(θ,ϕ),(40)
where Il and Kl 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 Pl is the l’th order Legendre polynomial, Eq. (3) converts into a form that is more convenient for numerical evaluation: Display Formula
Ψ(χ,χ)=S4πDl=01r<·r>Il+1/2(μeffr<)Kl+1/2(μeffr>)(2l+1)Pl(cosγ).(41)

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.01mm1 and μs=1.0mm1 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 (R0,θ,ϕ) and (R0+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 R0=50mm, 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 R0=200mm, 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.

Graphic Jump Location
Fig. 9
F9 :

Numerical comparison of the values given by Eqs. (41) and (2), at a common set of optical properties of μa=0.01mm1 and μs=1.0mm1, and different upper limit of l for the summation. The source and the detector are respectively at (R0,θ,ϕ) and (R0+d,θ,ϕ). The results are presented in terms of the linear relationship expected for ln(d*Ψ) versus d. At the radial coordinates of R0=50mm and R0=200mm, 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).

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.

Jöbsis  F. F., “Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science. 198, (4323 ), 1264 –1267 (1977). 0036-8075 CrossRef
Molavi  B.  et al., “Analyzing the resting state functional connectivity in the human language system using near infrared spectroscopy,” Front. Hum. Neurosci.. 7, , 921  (2014).CrossRef
Favilla  C. G.  et al., “Optical bedside monito