**=**

*μ*_{a}**0.01 mm**−

^{1 }and

**=**

*μ*_{s}′**1.0 mm**−

**are on average less than 5% different. The approximation for sphere, generally valid for a diameter ≥**

^{1}**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)^{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\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.^{7} at source–detector spacings $>4\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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.^{6} suggested a wavelength-, medium-, and geometry-dependent constant^{8} 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 head^{11} (e.g., prefrontal or parietal), and even different ages^{12} 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)^{15}16.^{–}^{17}], 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.^{18}19.20.21.^{–}^{22} 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 domain^{26} 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. 2526.–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. 2526.–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 ${\mu}_{a}$ and a diffusion coefficient of $D\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}1/(3{\mu}_{s}^{\prime})$ where ${\mu}_{s}^{\prime}$ 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.^{29}30.31.32.^{–}^{33} 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{\mu}_{s}^{\prime})$.

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 ${\overrightarrow{\chi}}^{\prime}=({x}^{\prime},{y}^{\prime},{z}^{\prime})$ to a detector at $\overrightarrow{\chi}=(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 ${\overrightarrow{\chi}}^{\prime}=({\rho}^{\prime},{\varphi}^{\prime},{z}^{\prime})$ to a detector at $\overrightarrow{\chi}=(\rho ,\varphi ,z)$. The sphere domain will be modeled in the spherical coordinates, where we are concerned about photon pathlength from a source at ${\overrightarrow{\chi}}^{\prime}=({r}^{\prime},{\theta}^{\prime},{\varphi}^{\prime})$ to a detector at $\overrightarrow{\chi}=(r,\theta ,\varphi )$.

## 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 ${\overrightarrow{\chi}}^{\prime}=(\mathrm{0,0},0)$, and the detector at $\overrightarrow{\chi}=(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 ${\overrightarrow{\chi}}^{\prime}=(\mathrm{0,0},0)$, and the detector at $\overrightarrow{\chi}=(0,d,0)$. The physical source at ${\overrightarrow{\chi}}^{\prime}=(\mathrm{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/{\mu}_{s}^{\prime}$, into the medium. The “real” isotropic source ${\overrightarrow{\chi}}_{\text{real}}^{\prime}$ thus has the coordinates of $(\mathrm{0,0},-{R}_{a})$. The effect of the medium-applicator interface on photon diffusion as modeled by the EBC^{21}^{,}^{23}^{,}^{24} sets zero the photon fluence rate on an imaginary boundary located at a distance of ${R}_{b}=2AD$ away from the physical boundary, where $A$ is a coefficient^{23} 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 ${\overrightarrow{\chi}}_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}}^{\text{'}}$, which has the coordinates of $(\mathrm{0,0},{R}_{a}+2{R}_{b})$. The distances from the detector to the “real” isotropic source ${\overrightarrow{\chi}}_{\text{real}}^{\prime}$ and the “image” source ${\overrightarrow{\chi}}_{\mathrm{imag}}^{\prime}$ are denoted by ${l}_{r}$ and ${l}_{i}$, respectively. The notations of ${l}_{r}=|\overrightarrow{\chi}-{\overrightarrow{\chi}}_{\text{real}}^{\prime}|$ and ${l}_{i}=|\overrightarrow{\chi}-{\overrightarrow{\chi}}_{\mathrm{imag}}^{\prime}|$ 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=|\overrightarrow{\chi}-{\overrightarrow{\chi}}^{\prime}|$ 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 $\overrightarrow{\chi}$ on the medium-applicator interface locates at $({R}_{0},\varphi ,z)$, and a directional source ${\overrightarrow{\chi}}^{\prime}$ on the medium-applicator interface locates at $({R}_{0},{\varphi}^{\prime},{z}^{\prime})$. The geometric symmetry requires that the directional source ${\overrightarrow{\chi}}^{\prime}$ be modeled by a “real” isotropic source ${\overrightarrow{\chi}}_{\text{real}}^{\prime}$ positioned along the radial direction of ${\overrightarrow{\chi}}^{\prime}$ and inwardly at a radial distance of ${R}_{a}=1/{\mu}_{s}^{\prime}$, at $({R}_{0}-{R}_{a},\varphi \prime ,z\prime )$. The extrapolated zero-boundary is concentric with and at a radial distance of ${R}_{b}=2AD$ 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.^{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.

## 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 $\overrightarrow{\chi}$ on the medium-applicator interface locates at $({R}_{0},\theta ,\varphi )$, and a directional source ${\overrightarrow{\chi}}^{\prime}$ on the medium-applicator interface locates at $({R}_{0},{\theta}^{\prime},{\varphi}^{\prime})$. The geometric symmetry requires that the directional source ${\overrightarrow{\chi}}^{\prime}$ be modeled by a “real” isotropic source ${\overrightarrow{r}}_{\text{real}}^{\prime}$ positioned along the radial direction of ${\overrightarrow{\chi}}^{\prime}$ and inwardly at a radial distance of ${R}_{a}=1/{\mu}_{s}^{\prime}$, [i.e., at $({R}_{0}-{R}_{a},{\theta}^{\prime},{\varphi}^{\prime})$]. The extrapolated zero boundary is concentric with and at a radial distance of ${R}_{b}=2AD$ 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

## (1)

$${\nabla}^{2}\mathrm{\Psi}(\overrightarrow{\chi})-\frac{{\mu}_{a}}{D}\mathrm{\Psi}(\overrightarrow{\chi})=-\frac{S(\overrightarrow{\chi})}{D},$$## (2)

$$\mathrm{\Psi}(\stackrel{\rightharpoonup}{\chi},{\overrightarrow{\chi}}^{\prime})=\frac{S}{4\pi D}\frac{1}{|\overrightarrow{\chi}-{\overrightarrow{\chi}}^{\prime}|}{e}^{-{\mu}_{\mathrm{eff}}|\overrightarrow{\chi}-{\overrightarrow{\chi}}^{\prime}|},$$## (3)

$$\mathrm{\Psi}(\stackrel{\rightharpoonup}{\chi},{\overrightarrow{\chi}}^{\prime})=\frac{S}{D}({\mu}_{\mathrm{eff}})\sum _{l=0}^{\infty}{i}_{l}({\mu}_{\mathrm{eff}}{r}_{<}){k}_{l}({\mu}_{\mathrm{eff}}{r}_{>})\sum _{m=-l}^{l}{Y}_{lm}^{*}({\theta}^{\prime},{\varphi}^{\prime}){Y}_{lm}(\theta ,\varphi ),$$## 2.5.

### 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

## (4)

$$\mathrm{\Psi}={\mathrm{\Psi}}_{\text{real}}-{\mathrm{\Psi}}_{\text{imag}}=\frac{S}{4\pi D{l}_{r}}{e}^{-{\mu}_{\mathrm{eff}}\hspace{0.17em}{l}_{r}}-\frac{S}{4\pi D{l}_{i}}{e}^{-{\mu}_{\mathrm{eff}}\hspace{0.17em}{l}_{i}},$$## 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^{25} 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},{\theta}^{\prime},{\varphi}^{\prime})$ and the detector is located at $({R}_{0},\theta ,\varphi )$, both on the physical boundary. The equivalent “real” isotropic source locates at $({R}_{0}-{R}_{a},{\theta}^{\prime},{\varphi}^{\prime})$, and the extrapolated boundary locates at a radial position of ${\rho}_{r>}={R}_{0}+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}_{<}={R}_{0}-{R}_{a}$ and the detector locates at ${r}_{>}={R}_{0}+{R}_{b}$, is

## (7)

$${{\mathrm{\Psi}}_{\text{real}}|}_{\mathrm{extr}}=\frac{1}{D}({\mu}_{\mathrm{eff}})\sum _{l=0}^{\infty}S\xb7{i}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}-{R}_{a})]{k}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{b})]\sum _{m=-l}^{l}{Y}_{lm}^{*}({\theta}^{\prime},{\varphi}^{\prime}){Y}_{lm}(\theta ,\varphi ),$$## (8)

$${\mathrm{\Psi}}_{\mathrm{imag}}{|}_{\mathrm{extr}}=\frac{1}{D}({\mu}_{\mathrm{eff}})\sum _{l=0}^{\infty}{S}_{l}^{*}{i}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{b})]{k}_{l}[{\mu}_{\mathrm{eff}}{r}_{>}]\sum _{m=-l}^{l}{Y}_{lm}^{*}({\theta}^{\prime},{\varphi}^{\prime}){Y}_{lm}(\theta ,\varphi ),$$^{35}

^{,}

^{36}the two unknown terms ${S}_{l}^{*}$ and ${r}_{>}\text{\hspace{0.17em}}\text{\hspace{0.17em}}$ 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

## (9)

$${S}_{l}^{*}{k}_{l}({\mu}_{\mathrm{eff}}{r}_{>})={S}_{l}{i}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}-{R}_{a})].$$Applying Eq. (9) to the EBC of ${\mathrm{\Psi}}_{\text{real}}{|}_{\mathrm{extr}}-{\mathrm{\Psi}}_{\mathrm{imag}}{|}_{\mathrm{extr}}=0$, we have

## (10)

$${S}_{l}=S\frac{{k}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{b})]}{{i}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{b})]}\phantom{\rule[-0.0ex]{1em}{0.0ex}}l=\mathrm{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 $\text{\hspace{0.17em}}\text{\hspace{0.17em}}{r}_{>}={R}_{0}$. For the “image” source also evaluating at the physical boundary, the detector locates at $\text{\hspace{0.17em}}\text{\hspace{0.17em}}{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

## (11)

$$\mathrm{\Psi}={\mathrm{\Psi}}_{\text{real}}{|}_{\mathrm{phys}}-{\mathrm{\Psi}}_{\mathrm{imag}}{|}_{\mathrm{phys}}=\frac{S}{D}({\mu}_{\mathrm{eff}})\sum _{l=0}^{\infty}{i}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}-{R}_{a})]{k}_{l}[{\mu}_{\mathrm{eff}}{R}_{0}]\sum _{m=-l}^{l}{Y}_{lm}^{*}({\theta}^{\prime},{\varphi}^{\prime}){Y}_{lm}(\theta ,\varphi )\phantom{\rule{0ex}{0ex}}-\frac{S}{D}({\mu}_{\mathrm{eff}})\sum _{l=0}^{\infty}{i}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}-{R}_{a})]{i}_{l}[{\mu}_{\mathrm{eff}}{R}_{0}]\frac{{k}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{b})]}{{i}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{b})]}\sum _{m=-l}^{l}{Y}_{lm}^{*}({\theta}^{\prime},{\varphi}^{\prime}){Y}_{lm}(\theta ,\varphi )=\frac{S}{D}({\mu}_{\mathrm{eff}})\sum _{l=0}^{\infty}{i}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}-{R}_{a})]{k}_{l}[{\mu}_{\mathrm{eff}}{R}_{0}]\sum _{m=-l}^{l}{Y}_{lm}^{*}({\theta}^{\prime},{\varphi}^{\prime}){Y}_{lm}(\theta ,\varphi )\u27e81-\frac{{i}_{l}[{\mu}_{\mathrm{eff}}{R}_{0}]}{{k}_{l}[{\mu}_{\mathrm{eff}}{R}_{0}]}\frac{{k}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{b})]}{{i}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{b})]}\u27e9.$$Equation (11) is convertible to the following form that may be more convenient for numerical evaluation:

## (12)

$$\mathrm{\Psi}=\frac{S}{D}\sqrt{\frac{1}{({R}_{0}-{R}_{a}){R}_{0}}}\sum _{l=0}^{\infty}{I}_{l+1/2}[{\mu}_{\mathrm{eff}}({R}_{0}-{R}_{a})]{K}_{l+1/2}({\mu}_{\mathrm{eff}}{R}_{0})\frac{2l+1}{4\pi}{P}_{l}(\mathrm{cos}\text{\hspace{0.17em}}\gamma )\phantom{\rule{0ex}{0ex}}\u27e81-\frac{{I}_{l+1/2}[{\mu}_{\mathrm{eff}}{R}_{0}]}{{K}_{l+1/2}[{\mu}_{\mathrm{eff}}{R}_{0}]}\frac{{K}_{l+1/2}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{b})]}{{I}_{l+1/2}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{b})]}\u27e9,$$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 ${\mu}_{a}=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and ${\mu}_{s}^{\prime}=0.25\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, the middle traces ${\mu}_{a}=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and ${\mu}_{s}^{\prime}=0.50\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, and the bottom traces ${\mu}_{a}=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and ${\mu}_{s}^{\prime}=1.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{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 $\mathrm{ln}({d}^{2}*\mathrm{\Psi})$ 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 ${\mu}_{s}^{\prime}=1.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{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^{37} 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.

## 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},{\theta}^{\prime},{\varphi}^{\prime})$, 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},{\theta}^{\prime},{\varphi}^{\prime})$. 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

## (13)

$${\mathrm{\Psi}}_{\mathrm{imag}}^{\mathrm{semi}}{|}_{\mathrm{phys}}=\frac{S}{D}({\mu}_{\mathrm{eff}})\sum _{l=0}^{\infty}{i}_{l}[{\mu}_{\mathrm{eff}}{R}_{0}]{k}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{a}+2{R}_{b})]\sum _{m=-l}^{l}{Y}_{lm}^{*}({\theta}^{\prime},{\varphi}^{\prime}){Y}_{lm}(\theta ,\varphi ).$$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

## (14)

$${\mathrm{\Psi}}_{\mathrm{imag}}{|}_{\mathrm{phys}}=\frac{S}{D}({\mu}_{\mathrm{eff}})\sum _{l=0}^{\infty}{i}_{l}[{\mu}_{\mathrm{eff}}{R}_{0}]{k}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{a}+2{R}_{b})]\xb7{\eta}_{l}\xb7\sum _{m=-l}^{l}{Y}_{lm}^{*}({\theta}^{\prime},{\varphi}^{\prime}){Y}_{lm}(\theta ,\varphi ),$$## (15)

$${\eta}_{l}=\frac{{k}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{b})]}{{i}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{b})]}\frac{{i}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}-{R}_{a})]}{{k}_{l}[{\mu}_{\mathrm{eff}}({R}_{0}+{R}_{a}+2{R}_{b})]}.$$Note that ${R}_{0}\gg {R}_{a}$, ${R}_{b}$ and if the argument of the modified Bessel functions is much greater than 1 (i.e., ${\mu}_{\mathrm{eff}}{R}_{0}\ge 10$ that corresponds to ${R}_{0}\ge 57.8\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ for ${\mu}_{a}=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and ${\mu}_{s}^{\prime}=1.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$) the modified spherical Bessel functions in Eq. (15) can be simplified by their asymptotic expressions,^{38} then Eq. (15) becomes

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

## (17)

$${\mathrm{\Psi}}_{\mathrm{imag}}{|}_{\mathrm{phys}}={\mathrm{\Psi}}_{\mathrm{imag}}^{\mathrm{semi}}{|}_{\mathrm{phys}}\xb7\frac{{R}_{0}+{R}_{a}+2{R}_{b}}{{R}_{0}-{R}_{a}}.$$Hence, for a sphere of large diameter, Eq. (12) approximates to

## (18)

$$\mathrm{\Psi}={\mathrm{\Psi}}_{\text{real}}{|}_{\mathrm{phys}}-{\mathrm{\Psi}}_{\mathrm{imag}}{|}_{\mathrm{phys}}={\mathrm{\Psi}}_{\mathrm{real}}{|}_{\mathrm{phys}}-{\mathrm{\Psi}}_{\mathrm{imag}}^{\mathrm{semi}}{|}_{\mathrm{phys}}\xb7\frac{{R}_{0}+{R}_{a}+2{R}_{b}}{{R}_{0}-{R}_{a}}.$$The ${\mathrm{\Psi}}_{\mathrm{real}}{|}_{\mathrm{phys}}$ of Eq. (18) is essentially the same as the ${\mathrm{\Psi}}_{\text{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 ${\mathrm{\Psi}}_{\mathrm{imag}}^{\mathrm{semi}}{|}_{\mathrm{phys}}$ of Eq. (18) will become the ${\mathrm{\Psi}}_{\mathrm{imag}}$ in Eq. (4) because the detector point locating at $({R}_{0},\theta ,\varphi )$ 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

## (19)

$$\mathrm{\Psi}=\frac{S}{4\pi D}\frac{{e}^{-{\mu}_{\mathrm{eff}}{l}_{r}}}{{l}_{r}}-\frac{S}{4\pi D}\frac{{e}^{-{\mu}_{\mathrm{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):

## (20)

$${l}_{r}=d\sqrt{1+{\mathrm{\Delta}}_{r}};\phantom{\rule[-0.0ex]{1em}{0.0ex}}\text{where}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{\Delta}}_{r}=\frac{{R}_{a}^{2}}{{d}^{2}}-\frac{{R}_{a}}{{R}_{0}},$$## (21)

$${l}_{i}=d\sqrt{1+{\mathrm{\Delta}}_{i}};\phantom{\rule[-0.0ex]{1em}{0.0ex}}\mathrm{where}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{\Delta}}_{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 ${\mu}_{\mathrm{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 ${\mu}_{\mathrm{eff}}{R}_{0}\ge 10$, that is, the radius (diameter) is more than 10 (20) times of the effective attenuation pathlength $1/{\mu}_{\mathrm{eff}}$.

## Table 1

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

Geometry | Fluence rate | Distance parameters |
---|---|---|

Infinite | ${\mathrm{\Psi}}_{\mathrm{inf}}=\frac{S}{4\pi D}\frac{{e}^{-{\mu}_{\mathrm{eff}}d}}{d}$ | $d$ |

Semi-infinite | ${\mathrm{\Psi}}_{\mathrm{semi}}=\frac{S}{4\pi D}\frac{{e}^{-{\mu}_{\mathrm{eff}}{l}_{r}}}{{l}_{r}}-\frac{S}{4\pi D}\frac{{e}^{-{\mu}_{\mathrm{eff}}{l}_{i}}}{{l}_{i}}$ | ${l}_{r}=d\sqrt{1+{\mathrm{\Delta}}_{r}}$@ ${R}_{0}\to \infty $${l}_{i}=d\sqrt{1+{\mathrm{\Delta}}_{i}}$ @ ${R}_{0}\to \infty $ |

Large infinite cylinder, azimuthal direction | ${\mathrm{\Psi}}_{\mathrm{azi}}=\frac{S}{4\pi D}\frac{{e}^{-{\mu}_{\mathrm{eff}}{l}_{r}}}{{l}_{r}}-\frac{S}{4\pi D}\frac{{e}^{-{\mu}_{\mathrm{eff}}{l}_{i}}}{{l}_{i}}\sqrt{\frac{{R}_{0}+{R}_{a}+2{R}_{b}}{{R}_{0}-{R}_{a}}}$ | ${l}_{r}=d\sqrt{1+{\mathrm{\Delta}}_{r}}$${l}_{i}=d\sqrt{1+{\mathrm{\Delta}}_{i}}$ |

Large sphere | ${\mathrm{\Psi}}_{\mathrm{sph}}=\frac{S}{4\pi D}\frac{{e}^{-{\mu}_{\mathrm{eff}}{l}_{r}}}{{l}_{r}}-\frac{S}{4\pi D}\frac{{e}^{-{\mu}_{\mathrm{eff}}{l}_{i}}}{{l}_{i}}\frac{{R}_{0}+{R}_{a}+2{R}_{b}}{{R}_{0}-{R}_{a}}$ | ${l}_{r}=d\sqrt{1+{\mathrm{\Delta}}_{r}}$${l}_{i}=d\sqrt{1+{\mathrm{\Delta}}_{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 $\mathrm{Att}=\mathrm{exp}[-{\mu}_{a}\xb7\mathrm{DPF}\xb7d]$, 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 $\mathrm{DPF}=\{\partial [\mathrm{ln}(\mathrm{S}/\mathrm{\Psi})]/\partial {\mu}_{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]

• Semi-infinite geometry [identical to Eq. (9) of Ref. 9]

## (23)

$${\mathrm{DPF}}_{\mathrm{semi}}=\frac{1}{d}\frac{\partial}{\partial {\mu}_{a}}\left[\mathrm{ln}\left(\frac{S}{{\mathrm{\Psi}}_{\mathrm{semi}}}\right)\right]=\frac{1}{d}\frac{\mathrm{exp}(-{\mu}_{\mathrm{eff}}{l}_{r})-\mathrm{exp}(-{\mu}_{\mathrm{ef}f}{l}_{i})}{\frac{\mathrm{exp}(-{\mu}_{\mathrm{eff}}{l}_{r})}{{l}_{r}}-\frac{\mathrm{exp}(-{\mu}_{\mathrm{eff}}{l}_{i})}{{l}_{i}}}\frac{\partial}{\partial {\mu}_{a}}({\mu}_{\mathrm{eff}})\approx \frac{\sqrt{3{\mu}_{s}^{\prime}}}{2\sqrt{{\mu}_{a}}}\frac{d\sqrt{3{\mu}_{s}^{\prime}{\mu}_{a}}}{d\sqrt{3{\mu}_{s}^{\prime}{\mu}_{a}}+1}.$$• Infinite-length cylindrical geometry, evaluated along the azimuthal direction

## (24)

$${\mathrm{DPF}}_{\mathrm{azi}}=\frac{1}{d}\frac{\partial}{\partial {\mu}_{a}}\left[\mathrm{ln}\left(\frac{S}{{\mathrm{\Psi}}_{\mathrm{azi}}}\right)\right]=\frac{1}{d}\frac{\mathrm{exp}(-{\mu}_{\mathrm{eff}}{l}_{r})-[\mathrm{exp}(-{\mu}_{\mathrm{eff}}{l}_{i})]\sqrt{\frac{{R}_{0}+{R}_{a}+2{R}_{b}}{{R}_{0}-{R}_{a}}}}{\frac{\mathrm{exp}(-{\mu}_{\mathrm{eff}}{l}_{r})}{{l}_{r}}-\frac{\mathrm{exp}(-{\mu}_{\mathrm{eff}}{l}_{i})}{{l}_{i}}\sqrt{\frac{{R}_{0}+{R}_{a}+2{R}_{b}}{{R}_{0}-{R}_{a}}}}\frac{\partial}{\partial {\mu}_{a}}({\mu}_{\mathrm{eff}})\phantom{\rule{0ex}{0ex}}=\frac{\sqrt{3{\mu}_{s}^{\prime}}}{2\sqrt{{\mu}_{a}}}\frac{d\sqrt{3{\mu}_{s}^{\prime}{\mu}_{a}}[{\eta}_{\mathrm{azi}}{\mathrm{\Delta}}_{i}-{\mathrm{\Delta}}_{r}]-[{\eta}_{\mathrm{azi}}-1]}{(d\sqrt{3{\mu}_{s}^{\prime}{\mu}_{a}}+1)[{\eta}_{\mathrm{azi}}{\mathrm{\Delta}}_{i}-{\mathrm{\Delta}}_{r}]-[{\eta}_{\mathrm{azi}}-1]},$$• Spherical geometry

## (26)

$${\mathrm{DPF}}_{\mathrm{sph}}=\frac{1}{d}\frac{\partial}{\partial {\mu}_{a}}\left[\mathrm{ln}\left(\frac{S}{{\mathrm{\Psi}}_{\mathrm{sph}}}\right)\right]=\frac{1}{d}\frac{\mathrm{exp}(-{\mu}_{\mathrm{eff}}{l}_{r})-[\mathrm{exp}(-{\mu}_{\mathrm{eff}}{l}_{i})]\frac{{R}_{0}+{R}_{a}+2{R}_{b}}{{R}_{0}-{R}_{a}}}{\frac{\mathrm{exp}(-{\mu}_{\mathrm{eff}}{l}_{r})}{{l}_{r}}-\frac{\mathrm{exp}(-{\mu}_{\mathrm{eff}}{l}_{i})}{{l}_{i}}\frac{{R}_{0}+{R}_{a}+2{R}_{b}}{{R}_{0}-{R}_{a}}}\frac{\partial}{\partial {\mu}_{a}}({\mu}_{\mathrm{eff}})\phantom{\rule{0ex}{0ex}}=\frac{\sqrt{3{\mu}_{s}^{\prime}}}{2\sqrt{{\mu}_{a}}}\frac{d\sqrt{3{\mu}_{s}^{\prime}{\mu}_{a}}[{\eta}_{\mathrm{sph}}{\mathrm{\Delta}}_{i}-{\mathrm{\Delta}}_{r}]-[{\eta}_{\mathrm{sph}}-1]}{(d\sqrt{3{\mu}_{s}^{\prime}{\mu}_{a}}+1)[{\eta}_{\mathrm{sph}}{\mathrm{\Delta}}_{i}-{\mathrm{\Delta}}_{r}]-[{\eta}_{\mathrm{sph}}-1]},$$

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 ${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),^{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 ${R}_{0}=25\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ and ${R}_{0}=50\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ with consistent node density of 1 mm that was comparable to the reduced scattering pathlength $(1/{\mu}_{s}^{\prime})$. The mesh for the semi-infinite geometry was generated for a $100\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}\times 100\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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).

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 ${\mathrm{\Psi}}_{1}$ measured at the 50 detector positions were forward calculated at the following setting of the optical properties: ${\mu}_{a}=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, ${\mu}_{s}^{\prime}=1.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$. (2) The second set of the photon fluence data ${\mathrm{\Psi}}_{2}$ measured at the same set of 50 detector positions were forward calculated with the ${\mu}_{a}$ at $0.011\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and other parameters remaining unchanged from the values of the previous step. (3) The DPF was calculated as $\{[\mathrm{ln}({\mathrm{\Psi}}_{1}/{\mathrm{\Psi}}_{2})]/\mathrm{\Delta}{\mu}_{a}\}/d$. Note that NIRFAST by default includes ${\mu}_{a}$ in the calculation of the diffusion coefficient. To make the FEM evaluation as consistent with the analytical approaches as would be available, the ${\mu}_{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.

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: ${\mu}_{a}=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, ${\mu}_{s}^{\prime}=1.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, ${R}_{0}=25\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, and ${R}_{0}=50\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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 ${\mu}_{\mathrm{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.^{12} 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$,^{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 ${\mu}_{a}=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, with the ${\mu}_{s}^{\prime}$ changing from $0.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ to $1.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$. The results are shown in Fig. 6, where it is observed that the three folds of increase of ${\mu}_{s}^{\prime}$ at the specific fixed value of ${\mu}_{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 ${\mu}_{s}^{\prime}$. The second case is to fix the reduced scattering coefficient at ${\mu}_{s}^{\prime}=1.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, with the ${\mu}_{a}$ changing from $0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ to $0.1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$. The results are shown in Fig. 7, where it is observed that the 10 folds of increase of ${\mu}_{a}$ at the given specific fixed ${\mu}_{s}^{\prime}$ 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 ${\mu}_{a}$ increases. The third case is to fix the ratio of the absorption coefficient to the reduced scattering coefficient to $1/100$, with the ${\mu}_{s}^{\prime}\text{\hspace{0.17em}}\text{\hspace{0.17em}}$ changing from $0.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ to $1.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{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 ${\mu}_{s}^{\prime}$ 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 ${\mu}_{s}^{\prime}$ at a fixed ${\mu}_{a}$ resulted in nearly 100% changes of DPF as shown in Fig. 6.

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\pm 0.42\text{\hspace{0.17em}}\text{\hspace{0.17em}}$ 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 ${\mu}_{s}^{\prime}$ at 835 nm changed from $0.43\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ to $0.79\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ for postconceptional age ranging from 30 to 42 weeks. The cerebral blood volume, which affected the global ${\mu}_{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 ${\mu}_{\mathrm{a}}=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, ${\mu}_{\mathrm{s}}^{\prime}=1.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{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.^{41} were experimentally determined by using the phase of the intensity modulated light, where the ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ were estimated by fitting the absolute DPF to a semi-infinite model. Therefore, when the sets of ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ were test plotted into a figure similar to Figs. 6Fig. 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 ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ of the neonatal brains could have been different.

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^{42} and the analytical solution of the diffusion equation incorporating the zero-boundary condition.^{18} However, there are experimental studies of the adult head^{43} suggesting that the DPF increases at short source–detector separation ($\sim 2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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)^{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 ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ 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 ${\mu}_{\mathrm{eff}}$ ${R}_{0}\ge 10$, that is, the radius (diameter) is more than 10 (20) times of the effective attenuation pathlength, $1/{\mu}_{\mathrm{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 ${\mu}_{a}=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and ${\mu}_{s}^{\prime}=1.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ is on average 5% different compared to the value obtained from FEM.

## Appendices

## Appendix A

Considering a source at ${\overrightarrow{x}}^{\prime}$ of $({r}^{\prime},{\theta}^{\prime},{\varphi}^{\prime})$ and a detector at $\overrightarrow{x}$ of $(r,\theta ,\varphi )$ in spherical coordinates, the equation for the Green’s function of Eq. (1) is

## (28)

$${\nabla}^{2}G(\overrightarrow{\chi},{\overrightarrow{\chi}}^{\prime})-{\mu}_{\mathrm{eff}}^{2}G(\overrightarrow{\chi},{\overrightarrow{\chi}}^{\prime})=-\delta (\overrightarrow{\chi}-{\overrightarrow{\chi}}^{\prime}).$$The Dirac delta function in Eq. (28) expressed in spherical coordinates (Ref. 38, pages 120, 427) is

## (29)

$$\delta (\overrightarrow{\chi}-{\overrightarrow{\chi}}^{\prime})=\frac{1}{{r}^{2}}\delta (r-{r}^{\prime})\delta (\mathrm{cos}\text{\hspace{0.17em}}\theta -\mathrm{cos}\text{\hspace{0.17em}}{\theta}^{\prime})\delta (\varphi -{\varphi}^{\prime}),$$## (30)

$$\delta (\overrightarrow{\chi}-{\overrightarrow{\chi}}^{\prime})=\frac{1}{{r}^{2}}\delta (r-{r}^{\prime})\sum _{l=0}^{\infty}\sum _{m=-l}^{l}{Y}_{lm}^{*}({\theta}^{\prime},{\varphi}^{\prime}){Y}_{lm}(\theta ,\varphi ).$$Substituting Eqs. (29) and (30) into Eq. (28) and expanding ${\nabla}^{2}$ in spherical coordinates lead to

## (31)

$$\frac{1}{r}\frac{{\partial}^{2}}{\partial {r}^{2}}[rG(\overrightarrow{\chi},{\overrightarrow{\chi}}^{\prime})]+\frac{1}{{r}^{2}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\theta}\frac{\partial}{\partial \theta}\left(\mathrm{sin}\text{\hspace{0.17em}}\theta \frac{\partial G(\overrightarrow{\chi},{\overrightarrow{\chi}}^{\prime})}{\partial \theta}\right)+\frac{1}{{r}^{2}\text{\hspace{0.17em}}{\mathrm{sin}}^{2}\text{\hspace{0.17em}}\theta}\frac{{\partial}^{2}G(\overrightarrow{\chi},{\overrightarrow{\chi}}^{\prime})}{\partial {\varphi}^{2}}-{\mu}_{\mathrm{eff}}^{2}G(\overrightarrow{\chi},{\overrightarrow{\chi}}^{\prime})=-\frac{1}{{r}^{2}}\delta (r-{r}^{\prime})\sum _{l=0}^{\infty}\sum _{m=-l}^{l}{Y}_{lm}^{*}({\theta}^{\prime},{\varphi}^{\prime}){Y}_{lm}(\theta ,\varphi ).$$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)]

## (32)

$$G(\stackrel{\rightharpoonup}{\chi},{\overrightarrow{\chi}}^{\prime})=\sum _{l=0}^{\infty}\sum _{m=-l}^{l}{g}_{l}(r,{r}^{\prime}){Y}_{lm}^{*}({\theta}^{\prime},{\varphi}^{\prime}){Y}_{lm}(\theta ,\varphi ),$$## (33)

$$\{\frac{{d}^{2}}{d{r}^{2}}+\frac{2}{r}\frac{d}{dr}-[{\mu}_{\mathrm{eff}}^{2}+\frac{l(l+1)}{{r}^{2}}\left]\right\}{g}_{l}=-\frac{1}{{r}^{2}}\delta (r-{r}^{\prime}),$$## (34)

$${g}_{l}=({\mu}_{\mathrm{eff}}){i}_{l}({\mu}_{\mathrm{eff}}{r}_{<}){k}_{l}({\mu}_{\mathrm{eff}}{r}_{>}),$$## (35)

$$G(\stackrel{\rightharpoonup}{\chi},{\overrightarrow{\chi}}^{\prime})=({\mu}_{\mathrm{eff}})\sum _{l=0}^{\infty}{i}_{l}({\mu}_{\mathrm{eff}}{r}_{<}){k}_{l}({\mu}_{\mathrm{eff}}{r}_{>})\sum _{m=-l}^{l}{Y}_{lm}^{*}({\theta}^{\prime},{\varphi}^{\prime}){Y}_{lm}(\theta ,\varphi ).$$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

## (36)

$$\mathrm{\Psi}(\stackrel{\rightharpoonup}{\chi},{\overrightarrow{\chi}}^{\prime})=\frac{S}{D}({\mu}_{\mathrm{eff}})\sum _{l=0}^{\infty}{i}_{l}({\mu}_{\mathrm{eff}}{r}_{<}){k}_{l}({\mu}_{\mathrm{eff}}{r}_{>})\sum _{m=-l}^{l}{Y}_{lm}^{*}({\theta}^{\prime},{\varphi}^{\prime}){Y}_{lm}(\theta ,\varphi ).$$## 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:

## (37)

$${i}_{l}({\mu}_{\mathrm{eff}}{r}_{<})=\sqrt{\frac{\pi}{2{\mu}_{\mathrm{eff}}{r}_{<}}}{I}_{l+1/2}({\mu}_{\mathrm{eff}}{r}_{<}),$$## (38)

$${k}_{l}({\mu}_{\mathrm{eff}}{r}_{<})=\sqrt{\frac{2}{\pi {\mu}_{\mathrm{eff}}{r}_{>}}}{K}_{l+1/2}({\mu}_{\mathrm{eff}}{r}_{>}),$$## (39)

$$\mathrm{cos}\text{\hspace{0.17em}}\gamma =\mathrm{cos}\text{\hspace{0.17em}}\theta \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}{\theta}^{\prime}+\mathrm{sin}\text{\hspace{0.17em}}\theta \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\theta}^{\prime}\text{\hspace{0.17em}}\mathrm{cos}(\varphi -{\varphi}^{\prime}),$$## (40)

$${P}_{l}(\mathrm{cos}\text{\hspace{0.17em}}\gamma )=\frac{4\pi}{2l+1}\sum _{m=-l}^{l}{Y}_{lm}^{*}({\theta}^{\prime},{\varphi}^{\prime}){Y}_{lm}(\theta ,\varphi ),$$## (41)

$$\mathrm{\Psi}(\stackrel{\rightharpoonup}{\chi},{\overrightarrow{\chi}}^{\prime})=\frac{S}{4\pi D}\sum _{l=0}^{\infty}\sqrt{\frac{1}{{r}_{<}\xb7{r}_{>}}}{I}_{l+1/2}({\mu}_{\mathrm{eff}}{r}_{<}){K}_{l+1/2}({\mu}_{\mathrm{eff}}{r}_{>})(2l+1){P}_{l}(\mathrm{cos}\text{\hspace{0.17em}}\gamma ).$$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 ${\mu}_{a}=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and ${\mu}_{s}^{\prime}=1.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{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},\theta ,\varphi )$ and $({R}_{0}+d,\theta ,\varphi )$, are presented in Fig. 9. The results are presented in terms of $\mathrm{ln}(d*\mathrm{\Psi})$ versus $d$, as it leads to a convenient linear relationship as evident from Eq. (2). At a radial coordinate of ${R}_{0}=50\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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.

## 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.