Novel data types for frequency-domain diffuse optical spectroscopy and imaging of tissues: characterization of sensitivity and contrast-to-noise ratio for absorption perturbations

Angelo Sassaroli; Giles Blaney; Sergio Fantini

doi:10.1364/BOE.485651

1. Introduction

In near-infrared spectroscopy (NIRS) there are three domains of investigations: continuous-wave (CW), frequency-domain (FD) and time-domain (TD). The latter two domains fall into the common definition of time-resolved domains. One of the advantages of these two domains is the possibility to measure the absolute values of the optical properties of tissues, namely the absorption coefficient ${\mu _a}$ and the reduced scattering coefficient $\mu {^{\prime}_s}$. Moreover, these two domains not only incorporate CW as data type, but also offer alternative data types that can be used for deep imaging of spontaneous or induced hemodynamic changes. In fact, it is well known that CW despite having the highest contrast-to-noise ratio (CNR) than any other data type, it is mostly sensitive to absorption changes occurring close to source and detector. Therefore, in an in vivo experiment where one is interested in probing/monitoring deep tissues, CW data will be likely polluted by confounds due to superficial hemodynamic changes. On the contrary, other data types both in TD and in FD are characterized by lower CNR but a different distribution of sensitivity. In particular, data types such as the mean time of flight in TD or the phase in FD, have a higher deep-to-superficial sensitivity ratio than CW data. In TD two “families” of data types have been used for imaging of deep hemodynamic changes: a) time gates [1] and b) higher order moments [2–4], specifically the mean time of flight ($\langle{t}\rangle$), and the variance $\left(\operatorname{Var}(t)=\left\langle t^2\right\rangle-\langle t\rangle^2\right)$ These two higher order moments have been used also for absolute recovery of the optical properties both from single-distance (SD) [5] or two distances [6]. Recently the use of spatial gradients in TD was also proposed as a method to increase the sensitivity of higher order moments to deep tissues [7,8]. Between the two higher order moments, several studies have shown the advantages of variance [7,9,10]. A comprehensive review of TD instrumentation and methods for functional imaging can be found in the work of Torricelli et al., [1].

The standard data types used in FD are AC amplitude, CW intensity and phase data ($\phi $). For typical modulation frequencies used in FD ($f \approx 100\; MHz$), AC and CW have similar sensitivity and CNR features to absorption changes. On the contrary it is well known that for typical optical properties of tissue (and typical modulation frequencies) phase data is substantially equivalent to mean time of flight: $\phi \approx \omega\langle t\rangle$ [11]. A note of caution is that this approximation becomes a poor one for lower values of the absorption coefficient and/or for higher values of the modulation frequency. Both AC and phase data have been used for measuring the absolute values of the optical properties and for the detection of hemodynamic changes in tissues. In particular, these data types have been regularly used in FD breast imaging [12], and sporadically for functional imaging of the brain [13]. In the last few years there has been a renewed interest in the use of phase data in functional studies of the brain. In particular, the advantages of incorporating phase data for image reconstruction in reflectance tomography has been demonstrated both theoretically [14] and in vivo during a protocol of brain activation [15]. Despite these successes of using phase data in functional FD NIRS, it is well known that its major limitation is the lower CNR than intensity data. The main question that we are trying to address in this work is: can we do any better than phase in FD NIRS? We think that it is possible, and the “recipe” we are proposing is to define new FD data types by using combinations of the three standard data types, that show higher sensitivity to deeper regions or better CNR features than phase. The guiding principle we are proposing is based on the observation that in FD one has the characteristic function (${{\rm X}_t}(\omega )$) of the photons’ flight time ($t$). The moments’ expansion of the real and imaginary part of the characteristic function is straightforward: the real part depends on the even moments, while the imaginary part depends on the odd moments, where the scaling factor of each moment $\langle{t^k}\rangle$ is ${\omega ^k}$. Therefore, we can imagine that at least for a range of interest of the optical properties, a few terms of the expansion will be appropriate to describe the real and imaginary parts. Linear combinations of the real and imaginary parts with phase data are used to enhance the role of some higher order moments or to reduce the effect of negatively interfering term(s) in a data type.

In the following section the theory of the novel FD data types is introduced (2.1). The calculation of Jacobians, and the definitions of sensitivity, depth-sensitivity, and deep-to-superficial sensitivity ratio, are provided in the following sections (2.2) and (2.3), respectively. In section (2.3) we define the sensitivity not only for the typical single-distance (SD) measurements, but also for the spatial gradients of FD data types, which we have previously named as dual-slope (DS) [16]. Next section (2.4) is dedicated to contrast-to-noise ratio (CNR) considerations. The results include the following sections. In section (3.1) the new data types and their fourth order moment expansion are compared. In section (3.2) the sensitivity of the higher order moments up to $\langle{t^3}\rangle$ are compared with the sensitivity of phase and CW data, for both SD and DS source-detector arrangements. The results of these section offer some guidelines for selecting six new FD data types that are compared with phase data in terms of sensitivity, in section (3.3) and CNR, in section (3.4). A comprehensive summary of the results of this investigation with some critical considerations will be provided in the discussion and conclusion.

2. Theory

2.1 Definition of new data types

In this section, the term “Green’s function” means (unless otherwise specified) the Green’s function of the output flux, which is obtained by applying Fick’s law to the Green’s function of the diffusion equation (DE) for the fluence rate. The latter Green’s function is obtained for a source type $\delta ({\mathbf r} )\delta (t )$ in TD and $\delta ({\mathbf r} )$ in CW and FD. It is known that given the time-domain (TD) Green’s function (${\mathrm{{\cal G}}_{TD}}(t )$, with dimensions $[{{\mathrm{{\cal G}}_{TD}}} ]= {L^{ - 2}}{t^{ - 1}}$), one can obtain the Green’s function for the complex reflectance $\tilde{R}(\omega )$ (i.e., the FD Green’s function: $\tilde{R}(\omega )= {\mathrm{{\cal G}}_{FD}}(\omega );\; [{{\mathrm{{\cal G}}_{FD}}} ]= {L^{ - 2}}$) by taking the Fourier transform of ${\mathrm{{\cal G}}_{TD}}(t )$. The k^th order moment of the photon’s arrival time (from the point source to a given point at the boundary of the medium) is defined by:

(1)$$\langle{t^k}\rangle = \mathop \smallint \nolimits_0^\infty {t^k} {\mathrm{\hat{{\cal G}}}_{TD}}(t )dt$$

where:

(2)$${\mathrm{\hat{{\cal G}}}_{TD}}(t )= \frac{{{\mathrm{{\cal G}}_{TD}}(t )}}{{\mathop \smallint \nolimits_0^\infty {\mathrm{{\cal G}}_{TD}}(t )dt}}$$

is the normalized TD Green’s function, which is the probability density (per unit time) for a detected photon to be detected at time t (i.e., $[{{{\mathrm{\hat{{\cal G}}}}_{TD}}} ]= {t^{ - 1}}$). The denominator of Eq. (2) is also called DC (for direct current), (i.e., the continuous wave (CW) Green’s function for the fluence rate $[{{\mathrm{{\cal G}}_{CW}}} ];\; [{{\mathrm{{\cal G}}_{CW}}} ]= {L^{ - 2}}$):

(3)$$DC = \mathop \smallint \nolimits_0^\infty {\mathrm{{\cal G}}_{TD}}(t )dt$$

We note that “CW” and “DC” have the same meaning and are used interchangeably. As a guiding principle to generate new FD data types, we note that if we consider the arrival time t of detected photons as a random variable (having ${\mathrm{\hat{{\cal G}}}_{TD}}(t )$ as its probability density function) its characteristic function ${{\rm X}_t}(\omega )$ is given by [17]:

(4)$${{\rm X}_t}(\omega )= {\langle e^{i\omega t}\rangle} = \frac{{\tilde{R}(\omega )}}{{DC}} = \frac{{\tilde{R}(\omega )}}{{\tilde{R}({\omega = 0} )}} = \left[ {\mathop \sum \nolimits_{k = 0}^\infty \frac{{{{({i\omega } )}^k}}}{{k!}}{\langle t^{k}\rangle}} \right]$$

where $\tilde{R}(\omega )= AC{e^{i\phi }}$ is the complex reflectance (AC is the intensity and $\phi $ the phase of the detected wave), $\; \omega = 2\pi f$ is the angular frequency ($f$ is the modulation frequency). We note that the characteristic function, as the moment generating function in CW (appendix A), depends on the entire series of infinite moments, with decreasing magnitude. However, the characteristic function can be divided into real and imaginary parts, which depend only on the even and odd moments, respectively. Therefore, by skipping every other term of the expansion we can expect that the first term of the real and imaginary parts has a stronger role in a data type, at least for some range of the optical properties and source-detector distances. Therefore, from Eq. (4) we can define the following two data types:

(5)$$1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]= 1 - \frac{{AC}}{{DC}}cos(\phi )= \mathop \sum \nolimits_{k = 0}^\infty \frac{{{{({ - 1} )}^k}{\omega ^{2k}}}}{{({2k)!} }}{\langle t^{2k}\rangle} \approx \frac{{{\omega ^2}}}{{2!}}{\langle t^{2}\rangle} - \frac{{{\omega ^4}}}{{4!}}{\langle t^{4}\rangle} + \mathrm{{\cal O}}({{\omega^6}{\langle t^{6}\rangle}} )$$

(6)$$\mathrm{\Im }[{{{\rm X}_t}(\omega )} ]= \frac{{AC}}{{DC}}sin(\phi )= \mathop \sum \nolimits_{k = 0}^{\infty} \frac{(-1)^k \omega^{2 k+1}}{(2 k+1) !}\left\langle t^{2 k+1}\right\rangle \approx \omega \langle t \rangle - \frac{{{\omega ^3}}}{{3!}}{\langle t^{3} \rangle} + \mathrm{{\cal O}}({{\omega^5}{\langle t^{5}\rangle}} )$$

Where $\mathrm{\Re }$ and $\mathrm{\Im }$ are the real and imaginary part of a complex number, respectively, and $\mathrm{{\cal O}}$ (i.e., “big O”) indicates the order of the next term in the expansion (therefore the order of the infinitesimal remainder when $\omega \to 0$). The data type defined in Eq. (5) contains two leading terms that are negatively interfering since they have opposite signs. However, in some situations (for example for higher values of the absorption coefficient of the medium or short source-detector separations) it may be possible that the contribution of the first term is much greater than the second. From Eq. (5) and (6) we can define different combinations such as:

(7)$$\mathrm{\Im }({{{\rm X}_t}(\omega )} )- \mathrm{\Re }({{{\rm X}_t}(\omega )} )+ 1 \approx \omega \langle t \rangle + \frac{{{\omega ^2}}}{{2!}}{\langle t^{2}\rangle} - \frac{{{\omega ^3}}}{{3!}}{\langle t^{3}\rangle} - \frac{{{\omega ^4}}}{{4!}}{\langle t^{4}\rangle} + \mathrm{{\cal O}}({{\omega^5}{\langle t^{5}\rangle}} )$$

where the first and second scaled moments are the leading terms and add together (positively interfering terms), while the third and the fourth moments are the negatively interfering terms (because of their negative sign). If we assume that $\phi \approx \omega\langle t\rangle$, we can define other data types such as:

(8)$$\phi - \mathrm{\Re }\left( {\frac{{\tilde{R}}}{{DC}}} \right) + 1 \approx \omega \langle t \rangle + {\; }\frac{{{\omega ^2}}}{{2!}}{\langle t^{2}\rangle} - \frac{{{\omega ^4}}}{{4!}}{\langle t^{4}\rangle} + \mathrm{{\cal O}}({{\omega^6}{\langle t^{6}\rangle}} )$$

which likely improves the sensitivity to deeper regions with respect to the previous data type, due to the absence of the negatively interfering term $\langle{t^3}\rangle$. Also, we can consider:

(9)$$\phi - \mathrm{\Im }\left( {\frac{{\tilde{R}}}{{DC}}} \right) \approx \frac{{{\omega ^3}}}{{3!}}{\langle t^{3}\rangle} + \mathrm{{\cal O}}({{\omega^5}{\langle t^{5}\rangle}} )$$

where the leading term should be more sensitive to deeper tissue regions than $\phi $ or $\langle{t}\rangle$.

We can think also of more elaborated combination such as:

(10)$$1 - \mathrm{\Re }\left( {\frac{{\tilde{R}}}{{DC}}} \right) - \frac{1}{2}{\phi ^2} \approx \frac{{{\omega ^2}}}{{2!}}Var(t )- \frac{{{\omega ^4}}}{{4!}}\langle{t^4}\rangle + \mathrm{{\cal O}}({{\omega^6}\langle{t^6}\rangle} )$$

In actuality, a further expansion of the phase which includes also third order moments, leads to different expansion of those data types where the phase is added or subtracted (see appendix B for a full list of data types, their fourth order expansion, and their Jacobians). We note that the expansion of FD data types in moments is intrinsic to the domain, therefore it is valid in any geometry and any distribution of the optical properties. The expansion of the FD data types in appendix B is valid also within the radiative transport equation (RTE). However, in this work we will focus on diffusive media characterized by the absorption and reduced scattering coefficients. In particular, we will choose the semi-infinite homogeneous medium geometry, and in appendix C we provide the expressions of the FD and TD Green’s function and moments up to the fourth order in this geometry. In summary, the novel data types differ between themselves and with respect to phase because of different moment expansions, therefore it is reasonable that they also have different spatial sensitivity properties to localized or global absorption perturbations.

2.2 Calculation of the Jacobians of the data types

The Jacobians of a data type with respect to a localized change of the absorption coefficient is defined by (${J_{Y,i}} = \frac{{\partial Y}}{{\partial {\mu _{ai}}}}$) where Y is a data type and ${\mu _{ai}}$ is a local value of the absorption coefficient. Also, the symbol ${J_Y}$ will be used for a uniform change of the absorption coefficient in the whole medium (i.e., $\; {J_Y} = \frac{{\partial Y}}{{\partial {\mu _a}}}$). If we divide a medium into N voxels, we have that ${J_Y} = \mathop \sum \nolimits_{i = 1}^N {J_{Y,i}}$. A voxel is a region of the medium where both the optical properties and the photon fluence rate can be considered constant. For the practical calculations of this work, we considered cubic voxels with 1 mm sides. For AC and CW we will consider the data types $Y = ln({{r^2}AC} )$ (where $r$ is the source-detector separation [18]) and similar definition for DC. By these definitions, all the data types considered in this work are unitless and the Jacobians have units of length ($[{{J_{Y,i}}} ]= L$), therefore we unify the notation by naming “generalized partial pathlengths” the opposite of the Jacobians (${\langle l\rangle _{Y,i}} \equiv{-} {J_{Y,i}}$). We note that the only Jacobian that corresponds to a physical partial or total pathlength of detected photons is the CW case (i.e., $Y = ln({{r^2}DC} )$). This observation also justifies the definition of generalized partial pathlength as the opposite of a Jacobian. The reason is that an increase in the absorption coefficient causes a decrease of CW and therefore a negative Jacobian. All the Jacobians of the novel FD data types (appendix B) can be derived from the FD one:

(11)$${\langle\tilde{l}\rangle_{FD,i}}({{{\mathbf r}_s},{{\mathbf r}_i},{{\mathbf r}_d},\omega } )={-} \frac{1}{{{\mathrm{{\cal G}}_{FD}}({{{\mathbf r}_s},{{\mathbf r}_d},\omega } )}}\frac{{\partial {\mathrm{{\cal G}}_{FD}}({{{\mathbf r}_s},{{\mathbf r}_d},\omega } )}}{{\partial {\mu _{ai}}}} ={-} {J_{FD,i}}({{{\mathbf r}_s},{{\mathbf r}_i},{{\mathbf r}_d},\omega } )$$

where ${{\mathbf r}_s}$, ${{\mathbf r}_d}$, ${{\mathbf r}_i}$ are the location of the source, the detector, and the focal absorption change, respectively. Note the tilde symbol at the top of the symbol (Eq. (11)) which is used for a complex quantity. From now on we will often omit the arguments of the Jacobians (and the pathlengths). We note that ${J_{FD,i}}$ is a concise expression that incorporates both generalized partial pathlengths of AC and phase [16]:

(12)$${\langle l \rangle_{AC,i}} = \mathrm{\Re }[{{\; }{{\langle\tilde{l}\rangle}_{FD,i}}} ]$$

(13)$${\langle l \rangle_{\phi ,i}} = \mathrm{\Im }[{{{\langle \tilde{l}\rangle}_{FD,i}}} ]$$

Also, the CW Jacobian ${J_{CW,i}}$ is a particular case of Eq. (11) for $\omega = 0$ (${J_{CW,i}} = {\; }{J_{FD,i}}({\omega = 0} )$). In a previous work [19] we have proposed an analytical formula for the calculation of the partial pathlengths in FD:

(14)$${\langle\tilde{l}\rangle_{FD,i}} \cong \frac{{{\mathrm{{\cal G}}_{FD,\mathrm{\Phi }}}({{{\mathbf r}_s},{{\mathbf r}_i},\omega } ){\mathrm{{\cal G}}_{FD}}\; ({{{\mathbf r}_i},{{\mathbf r}_d},\omega } )}}{{{\mathrm{{\cal G}}_{FD}}\; ({{{\mathbf r}_s},{{\mathbf r}_d},\omega } )}}d{V_i}$$

where ${\mathrm{{\cal G}}_{FD,\mathrm{\Phi }}}({{{\mathbf r}_s},{{\mathbf r}_i},\omega } )$, ${\mathrm{{\cal G}}_{FD}}\; ({{{\mathbf r}_i},{{\mathbf r}_d},\omega } )$, and ${\mathrm{{\cal G}}_{FD}}({{{\mathbf r}_s},{{\mathbf r}_d},\omega } )$ are the FD Green’s function for fluence rate calculated at ${{\mathbf r}_i}$ when photons are emitted from ${{\mathbf r}_s}$, the FD Green’s function for the output flux measured at the detector’s point $\; {{\mathbf r}_d}$ when photons are emitted from ${{\mathbf r}_i}$, and the FD Green’s function for the output flux measured at the detector’s point $\; {{\mathbf r}_d}$ when photons are emitted from ${{\mathbf r}_s}$, respectively; $d{V_i}$ is the infinitesimal volume of the voxel i (centered at ${{\mathbf r}_i}$). For the calculation of the partial pathlength in a finite region comprising many voxels, the differential expression in Eq. (14) must be integrated over the finite region of interest. The expressions of ${\mathrm{{\cal G}}_{FD,\mathrm{\Phi }}}({{{\mathbf r}_s},{{\mathbf r}_i},\omega } )$ and ${\mathrm{{\cal G}}_{FD}}\; ({{{\mathbf r}_s},{{\mathbf r}_d},\omega } )$ used for this work are given in appendix C. We also note that a lower case “$ l $” (such as: ${\langle l \rangle_{Y,i}}$) is used to indicate a generalized partial pathlength of the data type Y, while an upper case “$ L $” is used to indicate a generalized total pathlength (such as: ${J_Y} ={-} {\langle L \rangle_Y}$). The total generalized pathlength ${\langle L \rangle_Y}$ can be obtained in two ways: a) by directly applying Eq. (11) for a global uniform change of ${\mu _a}$; b) by integration of Eq. (14) over the whole medium.

2.3 Definition of single-distance (SD), single-slope (SS), and dual-slope (DS) sensitivities

We define sensitivity to localized absorption perturbations in terms of a normalized Jacobian. For example, the single-distance sensitivity (SD) of a data type Y, involves data collected with one source and one detector and is defined as ${S_{SDY,i}}$:

(15)$${S_{SDY,i}} = \frac{{{J_{Y,i}}({{{\mathbf r}_s},{{\mathbf r}_i},{{\mathbf r}_d},\omega } )}}{{{J_Y}({{{\mathbf r}_s},{{\mathbf r}_d},\omega } )}} = \frac{{{\langle l \rangle_{Y,i}}({{{\mathbf r}_s},{{\mathbf r}_i},{{\mathbf r}_d},\omega } )}}{{{\langle L \rangle_Y}({{{\mathbf r}_s},{{\mathbf r}_d},\omega } )}}$$

The single-slope sensitivity (SS) involves data collected with one source and two detectors; therefore, we are interested in calculating the normalized Jacobian of $SSY = \frac{{Y({{{\mathbf r}_s},{{\mathbf r}_{d2}},\omega } )- Y({{{\mathbf r}_s},{{\mathbf r}_{d1}},\omega } )}}{{|{{{\mathbf r}_{d2}} - {{\mathbf r}_{d1}}} |}}$ (i.e., the single slope of $Y$) which is defined as:

(16)$${S_{SSY,i}} = \frac{{{J_{Y,i}}({{{\mathbf r}_s},{{\mathbf r}_i},{{\mathbf r}_{d2}},\omega } )- {J_{Y,i}}({{{\mathbf r}_s},{{\mathbf r}_i},{{\mathbf r}_{d1}},\omega } )}}{{{J_Y}({{{\mathbf r}_s},{{\mathbf r}_{d2}},\omega } )- {J_Y}({{{\mathbf r}_s},{{\mathbf r}_{d1}},\omega } )}} = \frac{{{\langle l \rangle_{Y,i}}({{{\mathbf r}_s},{{\mathbf r}_i},{{\mathbf r}_{d2}},\omega } )- {\langle l \rangle_{Y,i}}({{{\mathbf r}_s},{{\mathbf r}_i},{{\mathbf r}_{d1}},\omega } )}}{{{\langle L \rangle_Y}({{{\mathbf r}_s},{{\mathbf r}_{d2}},\omega } )- {\langle L \rangle_Y}({{{\mathbf r}_s},{{\mathbf r}_{d1}},\omega } )}}$$

For AC data type ${S_{SSY,i}}$ is the sensitivity of $ln\left( {\frac{{r_2^2}}{{r_1^2}}\frac{{AC({{{\mathbf r}_s},{{\mathbf r}_{d2}},\omega } )}}{{AC({{{\mathbf r}_s},{{\mathbf r}_{d1}},\omega } )}}} \right)/|{{{\mathbf r}_{d2}} - {{\mathbf r}_{d1}}} |$, where $r_i^2 = {|{{{\mathbf r}_{di}} - {{\mathbf r}_s}} |^2}$. Similar definition is valid for the sensitivity of DC data. Finally, the dual-slope sensitivity (DS) of a data type Y [16], involves data collected with two sources and two detectors:

(17)$${S_{DSY,i}} = \frac{{{\langle l \rangle_{Y,i}}\left( {{{\mathbf r}_{s1}},{{\mathbf r}_i},{{\mathbf r}_{d2}},\omega } \right) - {\langle l \rangle_{Y,i}}\left( {{{\mathbf r}_{s1}},{{\mathbf r}_i},{{\mathbf r}_{d1}},\omega } \right) + {\langle l \rangle_{Y,i}}\left( {{{\mathbf r}_{s2}},{{\mathbf r}_i},{{\mathbf r}_{d1}},\omega } \right) - {\langle l \rangle_{Y,i}}\left( {{{\mathbf r}_{s2}},{{\mathbf r}_i},{{\mathbf r}_{d2}},\omega } \right)}}{{{\langle L \rangle_Y}\left( {{{\mathbf r}_{s1}},{{\mathbf r}_{d2}},\omega } \right) - {\langle L \rangle_Y}\left( {{{\mathbf r}_{s1}},{{\mathbf r}_{d1}},\omega } \right) + {\langle L \rangle_Y}\left( {{{\mathbf r}_{s2}},{{\mathbf r}_{d1}},\omega } \right) - {\langle L \rangle_Y}\left( {{{\mathbf r}_{s2}},{{\mathbf r}_{d2}},\omega } \right)}}$$

${S_{DSY,i}}$ is the sensitivity of the sum of two slopes, which involves two sources and two detectors. For $\textrm{AC and DC }$ it is the sensitivity of the slopes of two natural logarithms (as the one defined above) where each slope involves one source and two detectors. We note that the difference of the Jacobians for the two sources in Eq. (17) is taken in reverse order, because of the requirement that the far detector for one source must be the close detector for the other one [16]. For the symmetric DS arrangement used in this work, Eq. (17) simplifies into the average of the sensitivities of the two single slopes:

(18)$${S_{DSY,i}} = \frac{1}{2}({{S_{SSY,i}}({{{\mathbf r}_{s1}}} )+ {S_{SSY,i}}({{{\mathbf r}_{s2}}} )} )$$

where ${S_{SSY,i}}({{{\mathbf r}_{sj}}} )$ is Eq. (16) calculated for the source ${{\mathbf r}_{sj}}\; \; \; \; j = 1,2$.

It will be of interest to calculate also the depth sensitivity which can be defined for SD, SS, and DS arrangements of data types. Here we provide only the definition for SD:

(19)$${S_{SDY,z}} = {S_{SDY}}({{{\mathbf r}_s},z,{\; }{{\mathbf r}_d},\omega } )= \frac{{\frac{{\partial Y({{{\mathbf r}_s},\textrm{z},\; \; {{\mathbf r}_d},\omega } )}}{{\partial {\mu _{az}}}}}}{{\frac{{\partial Y({{{\mathbf r}_s},{{\mathbf r}_d},\omega } )}}{{\partial {\mu _a}}}}} = \frac{{{\langle l \rangle_{Y,z}}({{{\mathbf r}_s},z,{\; }{{\mathbf r}_d},\omega } )}}{{{\langle L \rangle_Y}({{{\mathbf r}_s},{{\mathbf r}_d},\omega } )}}$$

The Jacobian ${J_{Y,z}}$ (${J_{Y,z}} ={-} {\langle l\rangle_{Y,z}}$) is calculated with respect to the variation of the absorption coefficient in a “layer” at depth $\textrm{z}$ and having “infinitesimal” thickness. For the numerical calculations of this work, we used a finite layer having thickness of 1 mm and sides of 10 cm in both x and y directions, centered in the midpoint between source and detector. All the previous definitions of sensitivity apply also to the moments. In appendix D the general form of the moments’ Jacobians are provided for the SD case. From Eq. (19) by integration, we can define the deep-to-superficial sensitivity ratio (DSSR) which is the sensitivity to a uniform absorption change occurring at $z \ge {z_{top}}$ divided by the sensitivity to a uniform absorption change at $z \le {z_{top}}$, where ${z_{top}}$ is the thickness of a top layer:

(20)$$DSS{R_Y}({{{\mathbf r}_s},{\; }{{\mathbf r}_d},\omega } )= \frac{{\mathop \smallint \nolimits_{{z_{top}}}^\infty {S_{SDY}}({{{\mathbf r}_s},z,{\; }{{\mathbf r}_d},\omega } )dz}}{{\mathop \smallint \nolimits_0^{{z_{top}}} {S_{SDY}}({{{\mathbf r}_s},z,{\; }{{\mathbf r}_d},\omega } )dz}}$$

Having in mind the application of brain imaging in adults, we have chosen as the thickness of top layer ${z_{top}} = 15.5\; mm$, which is indicative of a typical thickness of the extracerebral layer. Finally, we note that all sensitivities are normalized in the whole medium to “1”.

2.4 Contrast-to-noise ratio considerations

The sensitivities allow for a comparison of noiseless data types in the space of all the parameters they depend on, the main ones being the source-detector separation(s), the optical properties, and the modulation frequency. However, they only provide information on the spatial distribution of the sensitivity of a given data type. They do not consider the magnitude of the change of a given data type as a result of a local absorption perturbation (contrast) in relation to the measurement error in that data type (noise). This information, which is of crucial importance in actual measurements, is obtained by the contrast-to-noise ratio (CNR). If one had a precise model of noise for the data types, one would not need to consider the sensitivity. In fact, the CNR, not only incorporates the sensitivity (note that the spatially normalized CNR is the sensitivity; see Eq. (21) below) but it tells also how feasible the detection of a perturbation with a data type is. In the lack of a precise noise model, we resort to phantom measurements, which are carried out only in a small subset of the space of possible parameters. In this work we have used a commercially available FD instrument (ISS Imagent V2, Champaign, IL) for the collection of data on highly scattering phantoms. The instrument operates at a modulation frequency $f = 140\; \textrm{MHz}$. The errors on the different data types were assessed directly as the standard deviations over 20-s acquisition on two phantoms. The first was a solid phantom with optical properties: ${\mu _a} = 0.0085\; \textrm{m}{\textrm{m}^{ - 1}}\; ,\mu {^{\prime}_s} = 1.05\; \textrm{m}{\textrm{m}^{ - 1}}$. The second one was a liquid phantom with optical properties: ${\mu _a} = 0.012\; \textrm{m}{\textrm{m}^{ - 1}}\; ,\mu {^{\prime}_s} = 0.56\; \textrm{m}{\textrm{m}^{ - 1}}$. We have used an optical probe with a linear arrangement of sources and detectors with the short and long source-detector separations of 25 and 35 mm, respectively. The data collection and analysis were done as follows. The data were acquired at sampling rate of 9.8 Hz. After detrending, the optical properties were determined at each time point by using an iterative fitting procedure [20]. We note that the precise error on the retrieved optical properties was not the focus of this work, therefore we did not investigate all sources of errors, especially the bias induced by the error in the source-detector distance and the error due to optical coupling between probe and phantom (the latter is smaller for the DS case). However, we point out that the precision of the optical properties in the 20 seconds acquisition was better than one part per million ($\frac{{{\sigma _\mu }}}{\mu }\; \le {10^{ - 6}}$, were $\mu $ is either the absorption or the reduced scattering coefficient). Therefore, the values of the optical properties were truncated to a typical significant digit that is usually found when all the sources of errors are considered. For the calculation of the CNR we considered a small variation of the absorption coefficient ($\Delta {\mu _{ai}}$) in a sub-region of a semi-infinite homogeneous medium. The variation of the data type (i.e., the contrast $\Delta Y$) was calculated by using a first order perturbation theory of the DE [19]. Therefore, the focal CNR for a data type Y was defined as:

(21)$$\textrm{CN}{\textrm{R}_Y} ={-} \frac{{\Delta Y}}{{{\sigma _Y}}} = \frac{{{\langle l \rangle_{Y,i}}\Delta {\mu _{ai}}}}{{{\sigma _Y}}}$$

which refers to the SD arrangement. By the definition provided by Eq. (21) we achieve two objectives: 1) we retain the information on the sign of the change of a data type for a localized increase of the absorption coefficient; 2) $\textrm{CN}{\textrm{R}_Y} > 0$ in most regions of the medium except in some subsets of a superficial region. Similar definitions apply to the SS and DS arrangements. We considered two types of absorption perturbations: a) a box-shaped region of size $16\; \textrm{mm}\; \times 16\; \textrm{mm} \times 2\; \textrm{mm}$ (realistic perturbation); b) a layered region of thickness $1\; \textrm{mm}$ (for the calculation of the depth-$CNR$). The $ \textrm{CNR}$ was calculated either in the plane $xz\; ({y = 0} )$ (case a) or at different depths (case b). The box-shaped region represents a realistic size of perturbation that may occur during brain activation. In both cases we considered $\Delta {\mu _{ai}} = 0.0033\; \textrm{m}{\textrm{m}^{ - 1}}$ which also reflects a realistic absorption change measured during brain activation [21].

3. Results

3.1 Comparison of FD data types and their fourth order moment expansion

In figure (1) thirteen data types are plotted (including AC, DC and phase) and their approximation up to the fourth order generalized moment (appendix B) as the function of source-detector separation. For the phase ($\phi $) we report not only the third order approximation but also $\omega\langle{t}\rangle$. The first term of the expansion of two data types, namely $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]- \frac{1}{2}{\phi ^2}$ (Eq. (30)) and $1 - \frac{{AC}}{{DC}}$ (Eq. (33)), is the scaled variance of the photon’s arrival time, therefore we report the fourth order approximation and $\frac{{{\omega ^2}}}{2}Var(t )$. The plots refer to typical values of the optical properties of tissues ${\mu _a} = 0.01\; m{m^{ - 1}},\; \mu {^{\prime}_s} = 1\; m{m^{ - 1}}$. In all the results we chose the inner and outer refractive indices as ${n_{in}} = 1.4,\; {n_{out}} = 1$, respectively, and a modulation frequency $f = 140\; MHz$. From figure (1) we can see a good agreement between data types and their approximations at shorter source-detector separations, with worse agreement at larger separations. This is expected since the higher order moments depend on the skewness of the TD Green’s function which becomes more pronounced at larger source-detector separation.

The skewness of the TD Green’s function becomes more evident also at smaller values of the absorption coefficient and higher scattering values. Therefore, the agreement between data types and any truncated expansion of the moments critically depends on three parameters: 1) the source-detector separation, 2) the optical properties, 3) the modulation frequency. It is expected to improve at shorter source-detector separations, higher values of the absorption coefficient, lower values of the reduced scattering coefficient, and lower modulation frequencies. For the last reason, all the results obtained in this work will apply also for modulation frequencies $f < 140\; MHz$. The least accurate agreement is found for high values of the reduced scattering coefficient and low values of the absorption coefficient. From Fig. 1 it turns out that $ \omega\langle t\rangle$ is already a good approximation of phase (discrepancy of less than ${0.8^ \circ }$ for $r \le 30\; mm$). About the two data types that depend on variance (panels (h), (k)), $1 - \frac{{AC}}{{DC}}$ (panel (k)) is a better approximation of variance (discrepancy of less than 8% for $r \le 30\; mm$). Similar agreement between the other data type (panel (h)) and variance is found only at very short source-detector separation ($r \le 13\; mm$).

Fig. 1. Data types (“true”) and their approximation (“app”) up to the fourth order of the generalized moments as the function of the source-detector separation (r). The optical properties of the medium and the inner and outer refractive indices are ${\mu _a} = 0.01\; m{m^{ - 1}},\; \mu {^{\prime}_s} = 1m{m^{ - 1}},$ ${n_{in}} = 1.4,\; {n_{out}} = 1$, respectively. We note that the order of the data types (a)…(l) is the same as that in appendix B. The modulation frequency used is $f = 140\; MHz$.

Download Full Size | PDF

3.2 Sensitivity of higher order moments

The study on the sensitivity was carried out for nine different combinations of the optical properties: ${\mu _a} \in \{{0.005,\; 0.01,\; 0.015} \}\; m{m^{ - 1}}$ and $\mu {^{\prime}_s} \in \{{0.5,\; 1,\; 1.5} \}\; m{m^{ - 1}}$ and for both SD and DS arrangements. Three source-detector separations were considered: 15, 25, 35 mm, which were used to define two DS sets, the first one with source-detector separations (25, 35) mm, and the second with (15, 25) mm. We show the plots only for the typical case ${\mu _a} = 0.01\; m{m^{ - 1}}$ and $\mu {^{\prime}_s} = 1\; m{m^{ - 1}}$, and for the SD (DS) we show only the results at 35 mm (25, 35 mm). A more comprehensive summary of all the results will be reported in the discussion. For the TD case we will show the sensitivity only for a “point” perturbation (a cube having $1\; mm$ sides), while for the FD case we will also show the depth sensitivity. These two extreme cases will give us some ideas about the sensitivities of more realistic spatial perturbations that occur during functional studies of brain activation.

Figure 2 shows a schematic of the geometry and the reference system used for the calculation of the sensitivities and CNR for both SD and DS source-detector arrangements. The DS is a linear symmetric arrangement with $|{{{\mathbf r}_{S1}} - {{\mathbf r}_{D1}}} |= |{{{\mathbf r}_{S2}} - {{\mathbf r}_{D2}}} |$ and $|{{{\mathbf r}_{S1}} - {{\mathbf r}_{D2}}} |= |{{{\mathbf r}_{S2}} - {{\mathbf r}_{D1}}} |$. In this section we compare the sensitivity of the moments $\langle{t^k}\rangle\; \; k = 1,\; 2,\; 3$, phase data, CW, variance ($Var(t )$), and one of the two pseudo-variances ($ \left(P \operatorname{var}(t)=\left\langle t^2\right\rangle-\frac{2}{3}\langle t\rangle^2\right) $; appendix B). For the SD case (Fig. 3) we can see some expected features of the sensitivity of the moments compared with those of phase and CW data. In particular, we point out the following: 1) CW has the highest sensitivity at depths less than about $10\; mm$. The CW data become increasingly less sensitive than other moments deeper in the medium; 2) the sensitivity of phase and $\langle t \rangle$ are very similar throughout the medium and they are the highest at an intermediate range of depths ($10\; mm\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }15\; mm$); 3) deeper in the medium ($z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel> \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }15\; mm$) $Var(t )$ becomes the most sensitive followed by $Pvar(t )$, $\langle{t^3}\rangle$ and $\langle{t^2}\rangle$. For the DS arrangement (Fig. 4) we note that: 1) the behavior of CW data and $Var(t )$ are similar to the SD case; 2) the intermediate range where $\langle t \rangle$, $\phi $ show the highest sensitivity (with the exception of $Var(t )$) extends to $z \approx 20\; mm$; 3) the hierarchy of sensitivities that was found in the SD case for $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel> \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }15\; mm$, is established deeper in the tissue. For example, $Pvar(t )$ has a higher sensitivity than $\langle t \rangle$ and $\phi $ only when $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel> \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }23\; mm$. For $\langle{t^3}\rangle$ and $\langle{t^2}\rangle$ this happens for $z > 26\; mm$.

Fig. 2. Schematic of the reference systems and source-detector arrangements for single-distance SD (left) and dual-slope DS (right) cases. ${S_i}$ are light sources while ${D_i}$ are collection points (detectors).

Download Full Size | PDF

Fig. 3. Sensitivity profiles of $\langle{t^k}\rangle\; \; k = 1,2,3$, phase data ($\phi $), CW, variance ($Var(t )$) and pseudo-variance $Pvar(t )$ for the SD arrangement (${S_{SD}},\; $ Eq. (15)) in a homogeneous semi-infinite medium geometry with ${\mu _a} = 0.01\; m{m^{ - 1}},\; \mu {^{\prime}_s} = 1\; m{m^{ - 1}},\; {n_{in}} = 1.4,\; {n_{out}} = 1,$ (where ${n_{in}}$ and ${n_{out}}$ are the internal and external refractive indices, respectively) The positions of the source ${S_1}$ and detector ${D_1}$ is $x = 0,\; x = 35\; mm$, respectively. The perturbation is a cubic defect with side $1\; mm$. Each subplot represents the sensitivity profiles along the x axis but at a different depth (title of the subplot).

Download Full Size | PDF

Fig. 4. Sensitivity profiles for the DS arrangement (${S_{DS}},\; $ Eq. (17)) for the same medium of Fig. 3. The positions of the sources ${S_1},{S_2}$ are $x = 0,\; x = 60\; mm$, respectively; the positions of the detectors ${D_1},{D_2}$ are $x = 25,\; x = 35\; mm$, respectively.

Download Full Size | PDF

We also note that the SD and DS configurations have different distributions of sensitivities, more spread for the SD case and more focal for the DS case. Finally, we note that except for SD CW data, both the SD and DS source-detector arrangements are characterized by regions of negative sensitivity which occur in a superficial layer of the medium, especially in a central region between source and detector for the SD case and between the two detectors for the DS case. The only exception for the DS case is CW data, having the region of negative sensitivity between each source and the closer detector.

These results agree with our intuitive understanding of the sensitivity of higher order moments for absorbing perturbations: the higher the moment the stronger sensitivity to deeper regions of diffusive media. They also give some indications about how to “sift” through the list of data types (appendix B) to select those which might have sensitivity or CNR features similar to or better than phase in deeper region of tissues. Since in FD a single higher order moment cannot be extracted, one method is to choose data types where the positive interference (i.e., terms adding together) between some higher order moments is much stronger than the negative interference (i.e., terms with opposite signs) with others. This can be assessed by the total generalized pathlength of a data type (${\langle L\rangle_Y} ={-} {J_Y}$) which gives an indication of the contrast for a uniform absorption change. For example, for a typical value of the optical properties of tissue (${\mu _a} = 0.01\; m{m^{ - 1}},\; \mu {^{\prime}_s} = 1\; m{m^{ - 1}}$) and at a source-detector distance $d = 35\; mm$, the data type $\mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ (Eq. (25)) has ${L_{\mathrm{\Im }[{{{\rm X}_t}(\omega )} ]}} = 9.64\; mm$. In the same situation, the data type $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$ has ${L_{1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]}} \approx 39.2\; mm$. The data type $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ is devised to eliminate the negative interference between the first and second term of the expansion of $\mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$. In fact, we obtain that ${L_{\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]}} \approx 33.05\; mm$. We also note that by adding together (or multiplying) two higher order moments (positive interference) we obtain a data type with intermediate sensitivity. By using the criterion of maximizing the total generalized pathlength, we have selected among the data types in appendix B six data types (with the highest ${\langle L\rangle_Y}$ in the list) that will be compared with phase.

This was done by having in mind applications of deep imaging of tissue (e.g., the adult brain) where the targeted absorption changes occur at a depth of 10-20 mm. It may be possible that other novel data types in appendix B can be used for imaging of tissue at an intermediate range of depths. The selected data types (and phase) are identified in appendix B by Eqs. (23, 24, 26,27, 28, 31, 33).

3.3 Sensitivity of FD novel data types

Figure 5 shows the sensitivity profiles for the SD arrangement of the selected FD data types compared with phase $\phi $. They were generated for the same medium and defect of Fig. 3. The data types $1 - \frac{{AC}}{{DC}}$ and $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ in agreement with their moment expansion (appendix B), have almost the same relationship with respect to phase as $Var(t )$ and $Pvar(t )$ of Fig. 3. We notice that for depths $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }15\; mm$ $\phi + \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ and $\mathrm{\Im }[{{{\rm X}_t}(\omega )} ]- \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$ have slightly higher sensitivity than other data types. This feature is more relevant for different combinations of the optical properties (data not shown). Another important feature of all the data types is the slight negative sensitivity which occurs in the central region (between source and detector) in the superficial region. Ideally the most favorable data types not only should have the highest sensitivity at depths $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel> \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }10\; mm$, but they should also have minimal negative sensitivity to a superficial region. Figure 6 shows the depth sensitivity for the SD arrangement and the medium of Fig. 3. We recall that in this case a layered perturbation (1 mm thick) at different depths ($z$) was considered. For this choice of the optical properties the only two data types that show more enhanced different features from phase are $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ and $1 - \frac{{AC}}{{DC}}$, which have an increased depth sensitivity at $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel> \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }13\; mm$. This happens for any combination of the optical properties considered in this study, but especially at lower $\mu {^{\prime}_s}$ or ${\mu _a}$ (results not shown).

Fig. 5. Sensitivity profiles of some new data types and phase for the SD arrangement (${S_{SD}},\; $ Eq. (15) for the medium and source-detector arrangement of Fig. 3.

Download Full Size | PDF

Fig. 6. Depth sensitivity profiles of some new data types and phase for the SD arrangement (${S_{SD}},\; $ Eq. (19)) for the medium and source-detector of Fig. 3.

Download Full Size | PDF

As we can see, the modulations of positive and negative sensitivity that were present in Fig. 5 close to the boundary, are missing in Fig. 6 and none of the data types shows a negative depth sensitivity. This happens because the region of negative sensitivity is overcompensated by the spikes of positive sensitivity which occurs at the locations of source and detector (Fig. 5). Therefore, we can imagine that for a realistic size of absorption perturbations which occur for example during brain activation, the negative region of sensitivity might play a minor role. About the deep-to-superficial sensitivity ratio (DSSR; Eq. (20)), four data types have larger values than phase: $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$, $\; \phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$, $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi $, and $1 - \frac{{AC}}{{DC}}$, which show an increase of this metric of 12%, 27%, 6% and 114%, respectively. Figure 7 shows the sensitivity profiles for the DS arrangement for the medium of Fig. 4. In this arrangement some data types feature larger oscillations (than for the SD case) between negative and positive sensitivity in the superficial region. This is due to the nonlinear behavior of the total Jacobians of some data types with respect to the source-detector separation (negative convexity, also known as concavity), which causes the denominator of Eq. (17) to be smaller than that of other data types. Especially this is true for $\phi + \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ and $\mathrm{\Im }[{{{\rm X}_t}(\omega )} ]- \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$. These two data types have also the highest sensitivity in a central region between the detectors down to a depth $z \approx 20\; mm$, where they are overcome by $1 - \frac{{AC}}{{DC}}$.

Fig. 7. Sensitivity profiles of some new data types and phase for the DS arrangement (${S_{DS}},\; $ Eq. (17)) for the medium and source-detector arrangement of Fig. 4.

Download Full Size | PDF

The sensitivity hierarchy that was present in the SD arrangement with $1 - \frac{{AC}}{{DC}}$ and $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ with highest sensitivity for $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel> \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }20\; mm$ is valid only for the first data type and is established deeper in tissue for the second one (data not shown). It is worth noticing that $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ is the data type showing the least oscillations between positive and negative sensitivity in the superficial region (Fig. 7 at $z = 5.5\; mm$). This is an important feature of this data type that can be useful when the limited number of SD or DS channels does not allow for a tomographic image reconstruction and only effective hemodynamic changes (which are assumed to occur uniformly in tissue) are measured. In this case, regions of negative sensitivity can cause misleading artifacts, e.g., a falsely measured effective decrease in absorption, when a true increase is present deeper in tissues. Figure 8 shows the depth sensitivity for the DS arrangement. We can see that a distinctive feature between data types is whether they have a negative sensitivity to a superficial region ($\phi $, $\; \phi + \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$, $\mathrm{\Im }[{{{\rm X}_t}(\omega )} ]- \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$, $1 - \frac{{AC}}{{DC}}$) or don’t ($\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$, $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$, $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi $). We note also that the maximum of phase data $\phi $ occurs deeper in the medium than the last three data types, however phase also retains some negative sensitivity in the first 5 mm of the medium. On the contrary the maxima of $\phi + \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$, $\mathrm{\Im }[{{{\rm X}_t}(\omega )} ]- \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$ and $1 - \frac{{AC}}{{DC}}$ occur deeper in tissue than phase data $\phi $, but these data types have a deeper region of negative sensitivity that stretches down to $z \approx 10\; mm$. Therefore, for these data types in the DS arrangement the region of negative sensitivity in the superficial region may play a bigger role (than for the SD case), even for realistic sizes of absorption perturbations occurring during brain activation. About the DSSR, three data types have higher values than phase:$\; \phi + \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$, $\mathrm{\Im }[{{{\rm X}_t}(\omega )} ]- \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ 1$, and $1 - \frac{{AC}}{{DC}}$, with an increase of 100%, 40% and 157%, respectively. However, these values are biased by the large negative sensitivity of these data types in the superficial region. A different definition of DSSR which considers the integral of the absolute sensitivity (see Eq. (20)), would show that only $\mathrm{\Im }[{{{\rm X}_t}(\omega )} ]- \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ 1$, and $1 - \frac{{AC}}{{DC}}$ have larger DSSR than phase, with an increase of 3% and 136% respectively.

Fig. 8. Depth sensitivity profiles of some new data types and phase for the DS arrangement for the medium of Fig. 4.

Download Full Size | PDF

3.4 Contrast-to-noise ratio of FD novel data types

In Fig. 9 are shown the standard deviations of the data types in the SD (left) and DS (right) arrangements for both solid (phantom 1) and liquid (phantom 2) phantoms having the optical properties:$\; {\mu _{a1}} = 0.0085\; m{m^{ - 1}}$, $\mu {^{\prime}_{s1}} = 1.05\; m{m^{ - 1}}$, and ${\mu _{a2}} = 0.012\; m{m^{ - 1}}$, $\mu {^{\prime}_{s2}} = 0.56\; m{m^{ - 1}}$, respectively. We note that beside the different units of the two standard deviations, ${\sigma _{DS}}$ is about ten times smaller than ${\sigma _{SD}}$ and this is expected because $|{{{\mathbf r}_{S1}} - {{\mathbf r}_{D2}}} |- |{{{\mathbf r}_{S1}} - {{\mathbf r}_{D1}}} |= 10\; mm$ (Fig. 2).

Fig. 9. Standard deviations of the FD data types for the SD (left; for d = 35 mm) and DS (right) arrangement, for the solid phantom (phantom 1) and the liquid phantom (phantom 2). The optical properties of the phantoms are: ${\mu _{a1}} = 0.0085\; m{m^{ - 1}}$, $\mu {^{\prime}_{s1}} = 1.05\; m{m^{ - 1}}$ for phantom1, and ${\mu _{a2}} = 0.012\; m{m^{ - 1}}$ $\mu {^{\prime}_{s2}} = 0.56\; m{m^{ - 1}}$ for phantom 2.

Download Full Size | PDF

In the SD arrangement the results can be summarized by considering the depth-CNR (Fig. 10). The data types that show a slightly higher depth-CNR (at all depths) than phase are: $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$, $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi $. The data type $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ shows a smaller depth-CNR than phase for $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }15\; mm$ but larger depth-CNR deeper in the medium. For the phantom 2 (data not shown) only $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi $ has a slightly higher depth-CNR than phase, and $\mathrm{\Im }({{{\rm X}_t}(\omega )} )- \mathrm{\Re }({{{\rm X}_t}(\omega )} )$ shows the same behavior of phase. These conclusions can also be drawn by scanning the realistic perturbation (data not shown). One point worth noting is that the negative sensitivity to a superficial region is not shown in Fig. 10, but it is retained for the realistic perturbation (as in Fig. 5; data not shown). The advantage of some novel data types is better shown in the DS arrangement. In Figs. 11 and 12 are shown the CNR for the phantom 1 and 2, respectively, and for the realistic perturbation.

Fig. 10. Depth-CNR in the SD arrangement for the medium and source-detector of Fig. 3 when the noise of phantom 1 is considered.

Download Full Size | PDF

Fig. 11. CNR in the DS arrangement of Fig. 4, for phantom 1 and the realistic perturbation (box-shaped region of size $16\; \textrm{mm}\; \times 16\; \textrm{mm} \times 2\; \textrm{mm}$ and $\Delta {\mu _{ai}} = 0.0033\; \textrm{m}{\textrm{m}^{ - 1}}$.

Download Full Size | PDF

Fig. 12. CNR in the DS arrangement of Fig. 4, for phantom 2 and the realistic perturbation (box-shaped region of size $16\; \textrm{mm}\; \times 16\; \textrm{mm} \times 2\; \textrm{mm}$ and $\Delta {\mu _{ai}} = 0.0033\; \textrm{m}{\textrm{m}^{ - 1}}$.

Download Full Size | PDF

First, we notice how all the data types, except for $\mathrm{\Im }({{{\rm X}_t}(\omega )} )- \mathrm{\Re }({{{\rm X}_t}(\omega )} )$ and $\phi + \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ have smaller oscillations between positive and negative sensitivity in the top region ($z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }6\; mm$) than phase data, with the data type $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ having the smallest oscillations.

Also, in the range $10\; mm\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }20\; mm$ three data types ($\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ],\; 1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$, $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi $) have more spread CNR features (but lower peak values) than phase data. This behavior leads (by integration) to a higher CNR in a layer located in the aforementioned range of depths, as it is shown in Figs. 13, 14. The pros and cons of having more spread CNR features will be discussed later. Finally, we can see that for a realistic perturbation the CNR of the data types is larger than 1 down to a depth $z \approx 20\; mm$. These results are consistent with depth-CNR shown in Fig. 13 and 14, for phantom 1 and 2, respectively. In both cases we can see that $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ],\; 1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$, $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi $ are the only data types not showing any negative CNR to a top layer. At the same time these data types have higher CNR than phase at least to a depth of $z \approx 20\; mm$. If we consider phantom 1, the increase of CNR in the whole medium with respect to phase is 19% for the data type $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$, 64% for $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$), and 11% for $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi $. If we consider the increase of CNR with respect to phase at depths $z \ge 15\; mm$ we have: 11% for $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$, 35% for $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$), and 7% for $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi $. For phantom 2, the increase of CNR in the whole medium with respect to phase is 18% for the data type $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$, 27% for $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$), and 16% for $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi $. If we consider the increase of CNR with respect to phase at depths $z \ge 15\; mm$ we have: 6% for $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$, 7% for $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$), and 10% for $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi $.

Fig. 13. Depth-CNR in the DS arrangement for the medium and source-detector of Fig. 4 when the noise of phantom 1 is considered.

Download Full Size | PDF

Fig. 14. Depth-CNR in the DS arrangement for the medium and source-detector of Fig. 4 when the noise of phantom 2 is considered.

Download Full Size | PDF

4. Discussion

In this work, we have proposed combinations of the standard FD data types (AC, DC and phase) to obtain novel data types that show better sensitivity or CNR than phase to deeper tissue perturbations in the absorption coefficient. The guiding principle was based on the observation that in FD one has the characteristic function of the photons’ arrival time, which has a straightforward moments’ expansion. We carried out a detailed study on the sensitivity features of the moments with the purpose of using these results to define novel FD data types. The study was carried out for nine different combinations of the optical properties and at three source-detector separations (defining two DS sets). The study confirmed that higher order moments have a higher sensitivity than phase and $\langle t \rangle$ in deeper regions of diffusive media $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel> \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }15 - 20\; mm$. More specifically, the investigated moments showing higher sensitivity than phase and $\langle t \rangle$ are$\; Var(t )$, $Pvar(t )$, $\langle{t^3}\rangle$ and $\langle{t^2}\rangle$ (in descending order). However, the depth where this property is established depends on the optical properties and the source-detector separations for both SD and DS cases. For the lowest values of the optical properties considered (${\mu _a} = 0.005\; m{m^{ - 1}},\; \mu {^{\prime}_s} = 0.5\; m{m^{ - 1}}$) and at a source-detector separation $d = 35\; mm$, this property is valid for $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel> \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }25\; mm$. On the contrary, for the largest value of the optical properties considered (${\mu _a} = 0.015\; m{m^{ - 1}},\; \mu {^{\prime}_s} = 1.5\; m{m^{ - 1}}$) the property is valid for $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel> \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }15\; mm$. Also, the depth where this property is valid is shifted upwards (towards smaller depths) when we consider $d = 25\; mm$. Similar conclusions are drawn for the DS case, however except for $Var(t )$ and $Pvar(t )$, $\langle{t^3}\rangle$ and $\langle{t^2}\rangle$ become more sensitive than phase deeper in the medium than for the SD case.

Based on these results, the way we combined the real and imaginary part of the characteristic function with phase data was such to obtain a higher order moment ${\omega ^k}\langle{t^k}\rangle\; k \ge 2$ (e.g., $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$) or a generalized third order moment (e.g., $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$) as the first term of the expansion. Other data types were devised to obtain the first and second order moments as first two positively interfering terms of the expansion. The sensitivity of the selected FD data types for a “point” defect is as follows.

1) The sensitivity relationship of the data type $1 - \frac{{AC}}{{DC}}$ to phase is similar to $Var(t )$ and phase in the previous study of the moments. This is true for both SD and DS cases.
2) Similarly, the data type $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ relates to phase as $Pvar(t )$ did in the previous study but only in the SD arrangement. In the DS arrangement the higher sensitivity of this data type with respect to phase is reached in deeper regions of the medium. One important property of this data type is that in both SD and DS cases it shows always the least negative sensitivity to a superficial region. Another important property is the increase in the SD deep-to-superficial sensitivity ratio (DSSR) with respect to phase by 41% and 27% at the source-detector separation of 25 and 35 mm, respectively.
3) The data types $\mathrm{\Im }({{{\rm X}_t}(\omega )} )- \mathrm{\Re }({{{\rm X}_t}(\omega )} )+ 1$ and $\phi + \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ have similar sensitivity features and they show the highest sensitivity down to a depth $z \approx 11 - 18\; mm$ (depending on the optical properties and source-detector separation) for both SD and DS cases. However, in the DS case this property is true only in a focal region between the two detectors. This behavior is typical for all the data types investigated: those having higher sensitivity in a central region between the detectors, have smaller sensitivity in the outer region (i.e., between each source and the closer detector), and vice versa. These data types are also characterized by a nonlinear global Jacobian with respect to the source-detector separation (negative convexity; data not shown) especially at higher values of the reduced scattering coefficient. This means that the global Jacobians at two source-detector separations can be very similar and this may cause misleading results, especially when one wants to retrieve effective hemodynamic changes (i.e., by assuming homogeneous absorption changes).
4) The data types $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$ and $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi $ have also similar sensitivity features. In the SD arrangement their sensitivities overcome that one of phase deeper in the medium $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel> \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }15 - 20\; mm$. This property depends to a lesser extent on the optical properties than on the source-detector separation. The global Jacobian of the data type $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$ also shows a nonlinear behavior with respect to the source-detector separation but to a lesser extent than the previous two data types (data not shown). The data type $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi $ is characterized by the highest global Jacobian in all the cases studied.

For depth-sensitivity we note the following:

1) For the SD case the data types that consistently and more markedly shows higher sensitivity than phase deeper in the medium is $1 - \frac{{AC}}{{DC}}$. This happens in the range $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel> \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }10 - $18 mm depending on the optical properties and source-detector separations. Other three data types $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$ and $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi $ and especially $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ shows slightly higher sensitivity deeper in the medium, and this property is more evident at lower values of the reduced scattering coefficient and for a source-detector separation d = 25 mm. The data types $\mathrm{\Im }({{{\rm X}_t}(\omega )} )- \mathrm{\Re }({{{\rm X}_t}(\omega )} )+ 1$ and $\phi + \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ have similar depth sensitivity of phase in most cases studied. The only exception is at lower values of the reduced scattering coefficient where they show higher sensitivity than phase for $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }15\; mm$. Finally, we point out that the SD sensitivities are positive in all the case studied.
2) For the DS case we note that the only data type showing always positive sensitivity at all depths is $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$. All other data types, in particular $\mathrm{\Im }({{{\rm X}_t}(\omega )} )- \mathrm{\Re }({{{\rm X}_t}(\omega )} )+ 1$ and $\phi + \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ and phase are characterized by a negative sensitivity in the superficial region $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }5 - 10\; mm$. For higher values of the reduced scattering coefficient these data types show wider oscillations between negative (in the superficial region) and positive values (deeper in the medium). This is a consequence of the nonlinearity of the total Jacobian with respect to the source-detector separation, which causes the denominator of Eq. (17) to be rather small. The data type $1 - \frac{{AC}}{{DC}}$ has similar behavior with respect to phase as the SD case. However, it is more affected by the region of negative sensitivity in the superficial region. Other data types show a higher sensitivity than phase for $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }10\; mm$.

In summary, the sensitivity study points out that there is not a single data type that prevails in all possible situations considered, and that careful consideration of the optical properties, source-detector arrangement (SD versus DS, including the distances used) and also location and size of a defect could lead to different choice of the data type used. Nevertheless, we have shown that for typical ranges of the optical properties and for both SD and DS arrangements some of the novel FD data types can have higher sensitivity than phase data. Also, we pointed out that in the DS arrangement some data types have higher sensitivity in a central region (between the two detectors) but lower in the outer region. Of course, this is a desirable property if one knows the exact location of a targeted hemodynamic change. Since, this is seldom the case, one must arrange source-detector arrays with a wide coverage of a region. Successful image reconstruction (with least artifacts) relies on the overlapping sensitivities of different source-detector units, and this is favored by data types having more spread sensitivity features.

The study on CNR has shown that:

1) For realistic perturbation and both SD and DS case, all the data types considered (except $1 - \frac{{AC}}{{DC}}$) are characterized by $CNR\; \ge 1$ for $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }20\; mm$.
2) Three data types $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]$, $1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi $ and especially $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ have higher CNR than phase at least down to $z\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel< \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }20\; mm$. If we look at Fig. 9 the latter data type has the lowest standard deviation (at least in the case of phantom 2). This means that $\phi $ and $\mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ are highly correlated and the individual errors do not add independently. This is true also for the other data type where phase is added or subtracted. These three data types show some promise for FD NIRS.

One final note is about the accuracy of the CNR results obtained by using first order perturbation theory of the DE, as we did in this work. We have repeated some calculations by using Monte Carlo (MC) simulations with a statistic of fifty thousand detected photons in each detector. We were able to replicate the comparison between the DS CNR of $\phi $ and $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ (Fig. 13). If we neglect the region $z \le 2\; mm$ for the case of $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$, and $z \le 6\; mm$ for the case of $\phi $, the average discrepancy between MC and DE calculations is 6% and 1.6%, respectively. If we consider $z \ge 6\; mm$ also for the first data type, the discrepancy between DE and MC calculations is 2.3%. These comparisons show that calculations carried out with first order perturbation theory of the DE are reliable, and that even more accurate calculation methods confirm the results shown.

5. Conclusions and future studies

We have shown that in FD some data types can be defined that for typical values of the optical properties of tissue show higher sensitivity and/or CNR for deeper absorption perturbations. In future investigation we intend to 1) extend the CNR study on more optical phantoms; 2) use the data types for image reconstruction in phantoms and in vivo; 3) consider optimal source-detector arrangements where the novel data types can have best performance; 4) define other data types that can be considered a “compromise” between the higher CNR of intensity data $ln({{r^2}AC} )$ and the higher sensitivity to deeper regions of the data types used in this study ($Y$), such as $Yln({{r^2}AC} )$. Ultimately, the goal of this study is to realize non-invasive optical measurements of biological tissues that provide selective sensitivity to deeper tissue. This is especially relevant for muscle and brain studies but has broader implications to depth-resolved measurements in other tissues as well.

Appendix A

Given the arrival time of a photon in a detector as a random variable, the moment generating function is given by:

(22)$$\frac{{{\mathrm{{\cal G}}_{CW}}({{\mu_a} + \Delta {\mu_a}} )}}{{{\mathrm{{\cal G}}_{CW}}({{\mu_a}} )}} = \langle{e^{ - \Delta {\mu _a}vt}}\rangle = \mathop \sum \nolimits_{k = 0}^\infty \frac{{{{({ - v\Delta {\mu_a}} )}^k}}}{{k!}}\langle{t^k}\rangle = \mathop \sum \nolimits_{k = 0}^\infty {({ - 1} )^k}\frac{{{{({\Delta {\mu_a}} )}^k}}}{{k!}}\langle{L^k}\rangle$$

In Eq. (23) $\langle{L^k}\rangle$ is the path length moment of k^th order. We recognize a formal similarity with the characteristic function, which is obtained when $- \Delta {\mu _a}v$ is replaced by $i\omega $.

Appendix B

We provide a full list of the data types considered in this work, together with the traditional $AC$, $DC$ and phase $\phi $. It is also provided a generalized fourth order moment approximation of each data type. By generalized moment of order n, we mean a moment $\langle{t^{k_1}\rangle}^{{m_1}}\langle{t^{k_2}\rangle}^{{m_2}}$ where ${k_1}{m_1} + \; {k_2}{m_2} = n$. We also note that $Var(t )= \langle{t^2}\rangle - \langle{t}\rangle^2$ is the variance of the arrival time of detected photons. One note is about the expansion of phase $\phi $ and $\frac{{AC}}{{DC}}$ beyond the first and second order, respectively. It is done by considering their expression in terms of $\mathrm{\Re }[{{{\rm X}_t}(\omega )} ]= x$ and $\mathrm{\Im }[{{{\rm X}_t}(\omega )} ]= y$ and carrying out a two-variable Taylor expansion around the point ($x,y$) = (1,0). The expansion in terms of $x,y$ is converted back into generalized moments of increasing order. In the following formulas we define $P {var}(t)=\left\langle t^2\right\rangle-\frac{2}{3}\langle t\rangle^2$, $Pva{r_1}(t )= \langle{t^2}\rangle - \frac{1}{2}\langle{t}\rangle^2$ (pseudo-variances) which are two recurring terms in the expansion. We remind that ${{\rm X}_t}(\omega )= \frac{{\tilde{R}(\omega )}}{{DC}} = \frac{{AC(\omega )}}{{DC}}{e^{i\phi (\omega )}}$.

(23)$$\phi \approx \omega \langle t\rangle - \frac{{{\omega ^3}}}{{3!}}\langle{t^3}\rangle + \frac{{{\omega ^3}}}{2}\langle t\rangle Pvar(t )+ \mathrm{{\cal O}}({{\omega^5}\langle{t^5}\rangle} )$$

(24)$$1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]= 1 - \frac{{AC}}{{DC}}cos(\phi )= \frac{{{\omega ^2}}}{2}\langle{t^2}\rangle - \frac{{{\omega ^4}}}{{4!}}\langle{t^4}\rangle + \mathrm{{\cal O}}({{\omega^6}\langle{t^6}\rangle} )$$

(25)$$\mathrm{\Im }[{{{\rm X}_t}(\omega )} ]= \frac{{AC}}{{DC}}sin(\phi )= \omega \langle t \rangle - \frac{{{\omega ^3}}}{{3!}}\langle{t^3}\rangle + \mathrm{{\cal O}}({{\omega^5}\langle{t^5}\rangle} )$$

(26)$$\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]= \frac{{{\omega ^3}}}{2}\langle t\rangle Pvar(t )+ \mathrm{{\cal O}}({{\omega^5}\langle{t^5}\rangle} )$$

(27)$$\phi + \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]= 2\omega \langle t \rangle + \frac{{{\omega ^3}}}{2}\langle t\rangle Pvar(t )- \frac{{{\omega ^3}}}{3}\langle{t^3}\rangle + \mathrm{{\cal O}}({{\omega^5}\langle{t^5}\rangle} )$$

(28)$$1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi = \omega \langle t \rangle + \frac{{{\omega ^2}}}{2}\langle{t^2}\rangle - \frac{{{\omega ^3}}}{{3!}}\langle{t^3}\rangle + \frac{{{\omega ^3}}}{2}\langle t\rangle Pvar(t )- \frac{{{\omega ^4}}}{{4!}}\langle{t^4}\rangle + \mathrm{{\cal O}}({{\omega^5}\langle{t^5}\rangle} )$$

(29)$$\mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi - 1 = \omega \langle t \rangle - \frac{{{\omega ^2}}}{2}\langle{t^2}\rangle - \frac{{{\omega ^3}}}{{3!}}\langle{t^3}\rangle + \frac{{{\omega ^3}}}{2}\langle t\rangle Pvar(t )+ \frac{{{\omega ^4}}}{{4!}}\langle{t^4}\rangle + \mathrm{{\cal O}}({{\omega^5}\langle{t^5}\rangle} )$$

(30)$$ 1-\mathfrak{R}\left[X_t(\omega)\right]-\frac{1}{2} \phi^2=\frac{\omega^2}{2} \operatorname{Var}(t)+\omega^4\left[-\frac{1}{4 !}\left\langle t^4\right\rangle+\frac{1}{3 !}\langle t\rangle\left\langle t^3\right\rangle-\frac{1}{2}\langle t\rangle^2 {Pvar}(t)\right]+\mathcal{O}\left(\omega^6\left\langle t^6\right\rangle\right) $$

(31)$$\mathrm{\Im }({{{\rm X}_t}(\omega )} )- \mathrm{\Re }({{{\rm X}_t}(\omega )} )+ 1 \approx \omega \langle t \rangle + \frac{{{\omega ^2}}}{{2!}}\langle{t^2}\rangle - \frac{{{\omega ^3}}}{{3!}}\langle{t^3}\rangle - \frac{{{\omega ^4}}}{{4!}}\langle{t^4}\rangle + \mathrm{{\cal O}}({{\omega^5}\langle{t^5}\rangle} )$$

(32)$$\mathrm{\Im }({{{\rm X}_t}(\omega )} )+ \mathrm{\Re }({{{\rm X}_t}(\omega )} )- 1 \approx \omega \langle t \rangle - \frac{{{\omega ^2}}}{{2!}}\langle{t^2}\rangle - \frac{{{\omega ^3}}}{{3!}}\langle{t^3}\rangle + \frac{{{\omega ^4}}}{{4!}}\langle{t^4}\rangle + \mathrm{{\cal O}}({{\omega^5}\langle{t^5}\rangle} )$$

(33)$$1 - \frac{{AC}}{{DC}} = \frac{{{\omega ^2}}}{2}Var(t )+ {\omega ^4}\left[ { - \frac{1}{{4!}}\langle{t^4}\rangle +\frac{1}{3 !}\langle t\rangle\left\langle t^3\right\rangle-\frac{1}{4}\langle t\rangle^2 Pvar_{1}(t )} \right] + \mathrm{{\cal O}}({{\omega^6}\langle{t^6}\rangle} )$$

(34)$$DC = \mathop \smallint \nolimits_0^\infty {\mathrm{{\cal G}}_{TD}}(t )dt = DC{t^0} = DC\mathop \smallint \nolimits_0^\infty {t^0}{\mathrm{\hat{{\cal G}}}_{TD}}(t )dt$$

(35)$$AC = DC\sqrt {{\mathrm{\Re }^2}[{{{\rm X}_t}(\omega )} ]+ {\mathrm{\Im }^2}[{{{\rm X}_t}(\omega )} ]} $$

These data types, except $AC$ and $DC$, are unitless. $AC$ data type does not have any fourth order approximation because we provide the approximation of $AC/DC$ (Eq. (33)). The Jacobians of the data types are provided below. Here only the global Jacobians are provided, but the expressions are valid also for local Jacobians by substitution of the generalized global pathlengths with the generalize partial pathlengths. Similar to the Jacobians for $AC$ and $\phi $, also the Jacobians of the two data types in (Eq. (24)) and (Eq. (25)) can be unified if we define a complex Jacobian $\langle\tilde{L}\rangle_{X_t}=X_t(\omega)\left[\langle\tilde{L}\rangle_{F D}-\langle L\rangle_{C W}\right]$.

(36)$$ \langle L\rangle_\phi=\mathfrak{J}\left\lceil\langle\tilde{L}\rangle_{F D}\right\rceil $$

(37)$$ \langle L\rangle_{1-\Re\left[X_t(\omega)\right]}=-\Re\left[\langle\tilde{L}\rangle_{X_t}\right] $$

(38)$$ \langle L\rangle_{\Im\left[X_t(\omega)\right]}=\Im\left[\langle\tilde{L}\rangle_{X_t}\right] $$

(39)$$ \langle L\rangle_{\phi-\Im\left[X_t(\omega)\right]}=\Im\left[\langle\tilde{L}\rangle_{F D}-\langle\tilde{L}\rangle_{X_t}\right] $$

(40)$$ \langle L\rangle_{\phi+\Im\left[X_t(\omega)\right]}=\mathfrak{\Im}\left[\langle\tilde{L}\rangle_{F D}+\langle\tilde{L}\rangle_{X_t}\right] $$

(41)$$ \langle L\rangle_{1-\Re\left[X_t(\omega)\right]+\phi}=\mathfrak{I}\left[\langle\tilde{L}\rangle_{F D}\right]-\Re\left[\langle\tilde{L}\rangle_{X_t}\right] $$

(42)$$ \langle L\rangle_{\Re\left[X_t(\omega)\right]+\phi-1}=\mathfrak{I}\left[\langle\tilde{L}\rangle_{F D}\right]+\Re\left[\langle\tilde{L}\rangle_{X_t}\right] $$

(43)$${\langle L\rangle_{1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]- \frac{1}{2}{\phi ^2}}} ={-} 2\varphi \mathrm{\Im }[{{\langle{\tilde{L}\rangle}_{FD}}} ]- \mathrm{\Re }[{{\langle \tilde{L}\rangle_{{{\rm X}_t}}}} ]$$

(44)$$ \langle L\rangle_{\Im\left(X_t(\omega)\right)-\Re\left(X_t(\omega)\right)+1}=\mathfrak{I}\left[\langle\tilde{L}\rangle_{X_t}\right]-\mathfrak{R}\left[\langle\tilde{L}\rangle_{X_t}\right] $$

(45)$$ \langle L\rangle_{\mathfrak{S}\left(X_t(\omega)\right)+\mathfrak{R}\left(X_t(\omega)\right)-1}=\mathfrak{I}\left[\langle\tilde{L}\rangle_{X_t}\right]+\mathfrak{R}\left[\langle\tilde{L}\rangle_{X_t}\right] $$

(46)$${\langle L\rangle_{1 - \frac{{AC}}{{DC}}}} ={-} \frac{{AC}}{{DC}}[{{\langle L\rangle_{AC}} - {\langle L\rangle_{CW}}} ]$$

(47)$${\langle L\rangle_{CW}} \equiv {\langle L\rangle_{DC}} \equiv \langle{\tilde{L}\rangle_{FD}}({\omega = 0} )$$

(48)$${\langle L\rangle_{AC}} = \mathrm{\Re }[{{\; }{{\langle \tilde{L}\rangle}_{FD}}} ]$$

The expression of ${\langle\tilde{L}\rangle_{FD}}$ used in this work (from which all the other global Jacobians are derived) is given in appendix C. The partial pathlengths (local Jacobians) can be derived from ${\langle\tilde{l}\rangle_{FD}}$ which is given by Eq. (14). For example, from Eq. (39) the opposite of the Jacobian of the data type $\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ]$ for a local change in the absorption coefficient in a voxel centered at ${{\mathbf r}_i}$ and volume $d{V_i}$ is:

(49)$$ -\frac{\partial\left(\phi-\Im\left[X_t(\omega)\right]\right)}{\partial \mu_{a i}}=\langle l\rangle_{\phi-\Im\left[X_t(\omega)\right], i}=\mathfrak{\Im}\left[\langle\tilde{l}\rangle_{F D, i}-\langle\tilde{l}\rangle_{X_t, i}\right] $$

where ${\langle\tilde{l}\rangle_{FD,i}}$ is given by Eq. (14) and ${\langle \tilde{l} \rangle_{{{\rm X}_t},i}}$:

(50)$${\langle \tilde{l}\rangle_{{{\rm X}_t},i}} = {{\rm X}_t}(\omega )[{{{\langle\tilde{l}\rangle}_{FD,i}} - {\langle l\rangle _{CW,i}}} ]$$

We remind that ${\langle l\rangle_{CW,i}}$ is a particular case of ${\langle\tilde{l}\rangle_{FD,i}}$ with $\omega = 0$. Equation (49) is valid for an arbitrary geometry and distribution of the optical properties of a medium. In the semi-infinite medium geometry with homogeneous optical properties the three Green’s functions of Eq. (14) are given in appendix C. More in detail: ${\mathrm{{\cal G}}_{FD,\mathrm{\Phi }}}({{{\mathbf r}_s},{{\mathbf r}_i},\omega } )$ of Eq. (14) is given by Eq. (61) with ${{\mathbf r}_i} = {\mathbf r}$. ${\mathrm{{\cal G}}_{FD}}\; ({{{\mathbf r}_s},{{\mathbf r}_d},\omega } )$ of Eq. (14) is given by Eq. (60). Last, ${\mathrm{{\cal G}}_{FD}}\; ({{{\mathbf r}_i},{{\mathbf r}_d},\omega } )$ of Eq. (14) is given formally by Eq. (60) provided that the positive source is at ${{\mathbf r}_i}$ instead of ${{\mathbf r}_s}$. Combining all the pieces we find the expression for ${\langle l\rangle_{\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ],i}}:$

(51)$$\begin{aligned} & {\langle l\rangle_{\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ],i}} = \\ & \begin{array}{l} {\Im }\left\{ {\frac{{\frac{1}{{4\pi D}}\left[ {\frac{{{e^{ - {{\tilde{\mu }}_{eff}}{r_1}}}}}{{{r_1}}} - \frac{{{e^{ - {{\tilde{\mu }}_{eff}}{r_2}}}}}{{{r_2}}}} \right]\left[ {{z_0}\left( {\frac{1}{{{r_{i1d}}}} + {{\tilde{\mu }}_{eff}}} \right)\frac{{{e^{ - {{\tilde{\mu }}_{eff}}\; {r_{i1d}}}}}}{{r_{i1d}^2}} + \left( {{z_0} + 2{z_b}} \right)\left( {\frac{1}{{{r_{i2d}}}} + {{\tilde{\mu }}_{eff}}} \right)\frac{{{e^{ - {{\tilde{\mu }}_{eff}}\; {r_{i2d}}}}}}{{r_{i2d}^2}}} \right]}}{{\left[ {{z_0}\left( {\frac{1}{{{r_{1d}}}} + {{\tilde{\mu }}_{eff}}} \right)\frac{{{e^{ - {{\tilde{\mu }}_{eff}}\; {r_{1d}}}}}}{{r_{1d}^2}} + \left( {{z_0} + 2{z_b}} \right)\left( {\frac{1}{{{r_{2d}}}} + {{\tilde{\mu }}_{eff}}} \right)\frac{{{e^{ - {{\tilde{\mu }}_{eff}}\; {r_{2d}}}}}}{{r_{2d}^2}}} \right]}}} \right.\left( {1 - {{{\rm X}}_t}\left( \omega \right)} \right) + \\ \left. {\frac{{\frac{1}{{4\pi D}}\left[ {\frac{{{e^{ - {\mu _{eff}}{r_1}}}}}{{{r_1}}} - \frac{{{e^{ - {\mu _{eff}}{r_2}}}}}{{{r_2}}}} \right]\left[ {{z_0}\left( {\frac{1}{{{r_{i1d}}}} + {\mu _{eff}}} \right)\frac{{{e^{ - {\mu _{eff}}\; {r_{i1d}}}}}}{{r_{i1d}^2}} + \left( {{z_0} + 2{z_b}} \right)\left( {\frac{1}{{{r_{i2d}}}} + {\mu _{eff}}} \right)\frac{{{e^{ - {\mu _{eff}}\; {r_{i2d}}}}}}{{r_{i2d}^2}}} \right]}}{{\left[ {{z_0}\left( {\frac{1}{{{r_{1d}}}} + {\mu _{eff}}} \right)\frac{{{e^{ - {\mu _{eff}}\; {r_{1d}}}}}}{{r_{1d}^2}} + \left( {{z_0} + 2{z_b}} \right)\left( {\frac{1}{{{r_{2d}}}} + {\mu _{eff}}} \right)\frac{{{e^{ - {\mu _{eff}}\; {r_{2d}}}}}}{{r_{2d}^2}}} \right]}}{X_t}\left( \omega \right)} \right\}d{V_i} \end{array} \end{aligned}$$

All the symbols of Eq. (51) are explained in appendix C and ${r_{i1d}}$, ${r_{i2d}}$ are the distances between the positive source (at the location of the defect ${{\mathbf r}_i}$) and the negative source (image source) with the detector, respectively. The coordinate of the positive and negative sources is ${{\mathbf r}_i} = ({{x_i},{y_i},{z_i}\; } )$ and $({{x_i},{y_i}, - {z_i} - 2{z_b}} )$, respectively, while the location of the detector is $({{x_d},0,\; 0\; } )$. More details are provided in appendix C. Equation (51) can be rewritten as:

(52)$${\langle l\rangle_{\phi - \mathrm{\Im }[{{{\rm X}_t}(\omega )} ],i}} = \mathrm{\Im }\{{{{\langle\tilde{l}\rangle}_{FD,i}}({1 - {{\rm X}_t}(\omega )} )+ {\langle l\rangle_{CW,i}}{{\rm X}_t}(\omega )} \}$$

From Eq. (52) it is easy to find the local Jacobians of the six data of interest:

(53)$${\langle l\rangle_{\phi + \mathrm{\Im }[{{{\rm X}_t}(\omega )} ],i}} = \mathrm{\Im }\{{{{\langle\tilde{l}\rangle}_{FD,i}}({1 + {{\rm X}_t}(\omega )} )- {\langle l\rangle_{CW,i}}{{\rm X}_t}(\omega )} \}$$

(54)$${\langle l\rangle_{1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]}} ={-} \mathrm{\Re }\{{{{\rm X}_t}(\omega )[{{{\langle\tilde{l}\rangle}_{FD,i}} - {\langle l\rangle_{CW,i}}} ]} \}$$

(55)$${\langle l\rangle_{1 - \mathrm{\Re }[{{{\rm X}_t}(\omega )} ]+ \phi }} = \mathrm{\Im }\{{{{\langle\tilde{l}\rangle}_{FD,i}}} \}- \mathrm{\Re }\{{{{\rm X}_t}(\omega )[{{{\langle\tilde{l}\rangle}_{FD,i}} - {\langle l\rangle_{CW,i}}} ]} \}$$

(56)$${\langle l\rangle_{\mathrm{\Im }({{{\rm X}_t}(\omega )} )- \mathrm{\Re }({{{\rm X}_t}(\omega )} )+ 1}} = \mathrm{\Im }\{{{{\rm X}_t}(\omega )[{{{\langle\tilde{l}\rangle}_{FD,i}} - {\langle l\rangle_{CW,i}}} ]} \}- \mathrm{\Re }\{{{{\rm X}_t}(\omega )[{{{\langle\tilde{l}\rangle}_{FD,i}} - {\langle l\rangle_{CW,i}}} ]} \}$$

(57)$${\langle l\rangle_{1 - \frac{{AC}}{{DC}}}} ={-} \frac{{AC}}{{DC}}[{{\langle l\rangle_{AC,i}} - {\langle l\rangle_{CW,i}}} ]={-} |{{{\rm X}_t}(\omega )} |\{{\mathrm{\Re }({{{\langle\tilde{l}\rangle}_{FD,i}}} )- {\langle l\rangle_{CW,i}}} \}$$

We stress again that the form of the local Jacobians provided by Eq. (52 - 57) (together with Eq. (14)) are valid in arbitrary conditions of the medium.

Appendix C

In the semi-infinite medium geometry and homogeneous optical properties, we choose a reference system with the origin where a pencil beam is impinging the medium (Fig. 2). The TD Green’s function for the fluence rate with extrapolated boundary conditions is given by:

(58)$${\mathrm{{\cal G}}_{TD,\mathrm{\Phi }}}({{{\mathbf r}_{\textrm{s},}}{\mathbf r},t} )= \frac{{v{e^{ - {\mu _a}vt}}}}{{{{({4\pi Dvt} )}^{\frac{3}{2}}}}}\left[ {{e^{ - \frac{{r_1^2}}{{4Dvt}}}} - {e^{ - \frac{{r_2^2}}{{4Dvt}}}}} \right]$$

where ${{\mathbf r}_\textrm{s}} = ({0,0,{z_0}} )$ is the location of the positive isotropic point source within the medium (and ${z_0} = \frac{1}{{\mu _s^{\prime}}}$), ${\mathbf r} = ({x,y,z} )$ is the field point inside the medium, $\; v$ the speed of light in the medium, ${\mu _a}$ is the absorption coefficient, $D = \frac{1}{{3\mu _s^{\prime}}}$ the diffusion factor ($\mu _s^{\prime}$ is the reduced scattering coefficient), ${z_b} = 2AD$ is the distance between real and extrapolated boundaries (where A is the reflection coefficient; [22]); also ${r_1}$, and ${r_2}$ are the distances between ${\mathbf r}$ and the positive ${{\mathbf r}_{1s}} = {{\mathbf r}_s}$ and the negative sources ${{\mathbf r}_{2s}} = ({0,0, - {z_0} - 2{z_b}} )$, respectively. By applying the Fick’s law to Eq. (58), we obtain the TD Green’s function ${\mathrm{{\cal G}}_{TD}}({\rho ,t} )$

(59)$${\mathrm{{\cal G}}_{TD}}({\rho ,t} )= \frac{1}{2}\frac{{{e^{ - {\mu _a}vt}}}}{{{{({4\pi Dv} )}^{\frac{3}{2}}}{t^{\frac{5}{2}}}}}\left[ {{z_0}{e^{ - \frac{{r_{1d}^2}}{{4Dvt}}}} + ({{z_0} + 2{z_b}} ){e^{ - \frac{{r_{2d}^2}}{{4Dvt}}}}} \right]$$

where ${r_{1d}}$, ${r_{2d}}$ are the distances between ${{\mathbf r}_{1s}}$ and ${{\mathbf r}_{2s}}$ and the detector ${{\mathbf r}_d} = ({{x_d},0,\; 0} )$, respectively, and $\rho = |{{x_d}} |$. The FD Green’s function ${\mathrm{{\cal G}}_{FD}}({\rho ,\omega } )$ for the output flux is given by:

(60)$${\mathrm{{\cal G}}_{FD}}({\rho ,\omega } )= \frac{1}{{4\pi }}\left[ {{z_0}\left( {\frac{1}{{{r_{1d}}}} + {{\tilde{\mu }}_{eff}}} \right)\frac{{{e^{ - {{\tilde{\mu }}_{eff}}\; {r_{1d}}}}}}{{r_{1d}^2}} + ({{z_0} + 2{z_b}} )\left( {\frac{1}{{{r_{2d}}}} + {{\tilde{\mu }}_{eff}}} \right)\frac{{{e^{ - {{\tilde{\mu }}_{eff}}\; {r_{2d}}}}}}{{r_{2d}^2}}} \right]$$

where ${\tilde{\mu }_{eff}} = \sqrt {\frac{{{\mu _a}}}{D} - i\frac{\omega }{{vD}}} $. Note that ${\mathrm{{\cal G}}_{FD}}({\rho ,\omega } )$ is ${\mathrm{{\cal G}}_{FD}}\; ({{{\mathbf r}_s},{{\mathbf r}_d},\omega } )$ of Eq. (14).

${\mathrm{{\cal G}}_{FD}}({\rho ,\omega } )$ is obtained by applying Fick’s law to the FD Green’s function for the fluence:

(61)$${\mathrm{{\cal G}}_{FD,\mathrm{\Phi }}}({{{\mathbf r}_s},{\mathbf r},\omega } )= \frac{1}{{4\pi D}}\left[ {\frac{{{e^{ - {{\tilde{\mu }}_{eff}}{r_1}}}}}{{{r_1}}} - \frac{{{e^{ - {{\tilde{\mu }}_{eff}}{r_2}}}}}{{{r_2}}}} \right]$$

where all the parameters are defined above. The normalized TD Green’s function, needed for the calculation of the moments is:

(62)$${\mathrm{\hat{{\cal G}}}_{TD}}({\rho ,t} )= \frac{{{\mathrm{{\cal G}}_{TD}}({\rho ,t} )}}{{{\mathrm{{\cal G}}_{CW}}(\rho )}} = \frac{{{\mathrm{{\cal G}}_{TD}}({\rho ,t} )}}{{{\mathrm{{\cal G}}_{FD}}({\rho ,\; \omega = 0} )}}$$

By applying the definition of moments (Eq. (1)) we obtain after some calculations:

(63)$$\langle t\rangle(\rho )= \frac{{\left[ {\frac{{{z_0}}}{{{r_1}}}{e^{ - {\mu_{eff}}{r_1}}} + \; \frac{{({{z_0} + 2{z_b}} )}}{{{r_2}}}{e^{ - {\mu_{eff}}{r_2}}}} \right]}}{{8\pi vD{\mathrm{{\cal G}}_{CW}}(\rho )}}$$

(64)$$\langle{t^2}\rangle(\rho )= \frac{{[{{z_0}{e^{ - {\mu_{eff}}{r_1}}} + ({{z_0} + 2{z_b}} ){e^{ - {\mu_{eff}}{r_2}}}} ]}}{{16\pi {{({vD} )}^2}{\mu _{eff}}{\mathrm{{\cal G}}_{CW}}(\rho )}}$$

(65)$$\langle{t^3}\rangle(\rho )= \frac{{\left[ {{z_0}\left( {{r_1} + \frac{1}{{{\mu_{eff}}}}} \right){e^{ - {\mu_{eff}}{r_1}}} + ({{z_0} + 2{z_b}} )\left( {{r_2} + \frac{1}{{{\mu_{eff}}}}} \right){e^{ - {\mu_{eff}}{r_2}}}} \right]}}{{32\pi {D^2}{v^3}{\mu _a}{\mathrm{{\cal G}}_{CW}}(\rho )}}$$

(66)$$\langle{t^4}\rangle(\rho )= \frac{{\left[ {{z_0}\left( {r_1^2 + \frac{{3{r_1}}}{{{\mu_{eff}}}} + \frac{3}{{\mu_{eff}^2}}} \right){e^{ - {\mu_{eff}}{r_1}}} + ({{z_0} + 2{z_b}} )\left( {r_2^2 + \frac{{3{r_2}}}{{{\mu_{eff}}}} + \frac{3}{{\mu_{eff}^2}}} \right){e^{ - {\mu_{eff}}{r_2}}}} \right]}}{{64\pi {D^{\frac{5}{2}}}{{({{\mu_a}} )}^{\frac{3}{2}}}{v^4}{\mathrm{{\cal G}}_{CW}}(\rho )}}$$

where ${\mu _{eff}} = {\tilde{\mu }_{eff}}({\omega = 0} )$. If we assume a constant speed of light in the medium the previous equations yield also the pathlength moments. Note that the above expressions are also valid for the complex pathlength in FD, in particular ${\langle\tilde{L}\rangle_{FD}} = v\tilde{t}$ (appendix B), where $\tilde{t}$ is given by Eq. (63) when we use ${\tilde{\mu }_{eff}}$ instead of ${\mu _{eff}}$ and ${\mathrm{{\cal G}}_{CW}}(\rho )$ is substituted by ${\mathrm{{\cal G}}_{FD}}({\rho ,\omega } )$.

Appendix D

In this appendix we provide the expression of the Jacobian of the moments $\langle{t^k}\rangle$. The general expression of these Jacobians for a photon emitted at ${{\mathbf r}_s}$ and detected at ${{\mathbf r}_d}$ is:

(67)$$ \frac{\partial\left\langle t^k\right\rangle\left(\mathbf{r}_s, \mathbf{r}_d\right)}{\partial \mu_{a i}}=\left\langle l_i\right\rangle_{C W}\left(\mathbf{r}_s, \mathbf{r}_d\right)\left\langle t^k\right\rangle\left(\mathbf{r}_s, \mathbf{r}_d\right)-\frac{1}{D C\left(\mathbf{r}_s, \mathbf{r}_d\right)} \int_0^{\infty} t^k \mathcal{G}_{T D}\left(\mathbf{r}_s, \mathbf{r}_d, t\right)\left\langle l_i\right\rangle\left(t, \mathbf{r}_s, \mathbf{r}_d\right) d t $$

where $ \left\langle l_i\right\rangle\left(t, \mathbf{r}_s, \mathbf{r}_d\right) $ is the time-resolved mean pathlength in a sub-volume “$ i $” of the diffusive medium:

(68)$$\langle{l_i}\rangle(t,{{{\mathbf r}_s},{{\mathbf r}_d}} )= \frac{{\mathop \smallint \nolimits_{{V_i}}^{} [{{\mathrm{{\cal G}}_{TD,\mathrm{\Phi }}}({{{\mathbf r}_s},{{\mathbf r}_i},\tau } )\otimes {\mathrm{{\cal G}}_{TD}}\; ({{{\mathbf r}_i},{{\mathbf r}_d},\tau } )} ](t )d{{\mathbf r}_i}}}{{{\mathrm{{\cal G}}_{TD}}\; ({{{\mathbf r}_s},{{\mathbf r}_d},t} )}}$$

where ${\mathrm{{\cal G}}_{TD,\mathrm{\Phi }}}({{{\mathbf r}_s},{{\mathbf r}_i},\tau } )$ is the TD Green’s functions for the fluence rate (Eq. (58)) for a photon emitted at ${{\mathbf r}_s}$ at $\tau $=0 and found at ${{\mathbf r}_i}$ at time $\tau $; ${\mathrm{{\cal G}}_{TD}}\; ({{{\mathbf r}_i},{{\mathbf r}_d},\tau } )$ is the TD Green’s functions for the exiting flux (Eq. (59)) for a photon emitted at ${{\mathbf r}_i}$ at time $\tau $=0 and detected at ${{\mathbf r}_d}$ at time $\tau $; ${\mathrm{{\cal G}}_{TD}}\; ({{{\mathbf r}_s},{{\mathbf r}_d},t} )$ is the TD Green’s functions for the exiting flux for a photon emitted at ${{\mathbf r}_s}$ at time $t$=0 and detected at ${{\mathbf r}_d}$ at time $\langle t \rangle$; “$ \otimes$” is the convolution operator. Equation (68) is valid for a finite region within a diffusive medium. For a voxel of volume $d{V_i}$ Eq. (67) simplifies into:

(69)$$\frac{{\partial \langle{t^k}\rangle({{{\mathbf r}_s},{{\mathbf r}_d}} )}}{{\partial {\mu _{ai}}}} = \langle{l_i}\rangle_{CW}\langle{t^k}\rangle({{{\mathbf r}_s},{{\mathbf r}_d}} )- \frac{{d{V_i}}}{{DC({{{\mathbf r}_s},{{\mathbf r}_d}} )}}\mathop \smallint \nolimits_0^\infty {t^k}[{{\mathrm{{\cal G}}_{TD,\mathrm{\Phi }}}({{{\mathbf r}_s},{{\mathbf r}_i},\tau } )\otimes {\mathrm{{\cal G}}_{TD}}\; ({{{\mathbf r}_i},{{\mathbf r}_d},\tau } )} ](t )dt\; $$

The Jacobian for uniform absorption changes in the medium is:

(70)$$ \frac{\partial\left\langle t^k\right\rangle}{\partial \mu_a}=\langle L\rangle_{C W}\left\langle t^k\right\rangle-v\left\langle t^{k+1}\right\rangle=v\left(\langle t\rangle\left\langle t^k\right\rangle-\left\langle t^{k+1}\right\rangle\right) $$

The calculation of the sensitivity functions for $\langle{t^k}\rangle$ was carried out by using a numerical integration of Eq. (69) with a time axis having a step of 1 ps and 12 ns long.

Funding

National Institutes of Health (R01 EB029414, R01 NS095334).

Acknowledgements

We thank Dr. Tapan Das for the collection of the data on the solid phantom.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Supplemental document

See Supplement 1 for supporting content.

References

1. A. Torricelli, D. Contini, A. Pifferi, M. Caffini, R. Re, L. Zucchelli, and L. Spinelli, “Time domain functional NIRS imaging for human brain mapping,” NeuroImage 85, 28–50 (2014). [CrossRef]

2. A. Liebert, H. Wabnitz, J. Steinbrink, H. Obrig, M. Moller, R. Macdonald, A. Villringer, and H. Rinneberg, “Time-resolved multidistance near-infrared spectroscopy of the adult head: intracerebral and extracerebral absorption changes from moments of distribution of times of flight of photons,” Appl. Opt. 43(15), 3037–3047 (2004). [CrossRef]

3. M. Kacprzak, A. Liebert, P. Sawosz, N. Zolek, and R. Maniewski, “Time-resolved optical imager for assessment of cerebral oxygenation,” J. Biomed. Opt. 12(3), 034019 (2007). [CrossRef]

4. J. Steinbrink, H. Wabnitz, H. Obrig, A. Villringer, and H. Rinneberg, “Determining changes in NIR absorption using a layered model of the human head,” Phys. Med. Biol. 46(3), 879–896 (2001). [CrossRef]

5. A. Liebert, H. Wabnitz, D. Grosenick, M. Möller, R. Macdonald, and H. Rinneberg, “Evaluation of optical properties of highly scattering media by moments of distributions of times of flight of photons,” Appl. Opt. 42(28), 5785–5792 (2003). [CrossRef]

6. D. Milej, A. Abdalmalak, D. Janusek, M. Diop, A. Liebert, and K. St. Lawrence, “Time-resolved subtraction method for measuring optical properties of turbid media,” Appl. Opt. 55(7), 1507–1513 (2016). [CrossRef]

7. P. Sawosz and A. Liebert, “Method to improve the depth sensitivity of diffuse reflectance measurements to absorption changes in optically turbid medium,” Biomed. Opt. Express 10(10), 5031–5041 (2019). [CrossRef]

8. D. Milej, A. Abdalmalak, P. McLachlan, M. Diop, A. Liebert, and K. St Lawrence, “Subtraction-based approach for enhancing the depth sensitivity of time-resolved NIRS,” Biomed. Opt. Express 7(11), 4514–4526 (2016). [CrossRef]

9. A. Liebert, H. Wabnitz, and C. Elster, “Determination of absorption changes from moments of distributions of times of flight of photons: optimization of measurement conditions for a two-layered tissue model,” J. Biomed. Opt. 17(5), 057005 (2012). [CrossRef]

10. H. Wabnitz, D. Contini, L. Spinelli, A. Torricelli, and A. Liebert, “Depth-selective data analysis for time-domain fNIRS: moments vs. time windows,” Biomed. Opt. Express 11(8), 4224–4243 (2020). [CrossRef]

11. S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol. 37(7), 1531–1560 (1992). [CrossRef]

12. H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Meth. Eng. 25(6), 711–732 (2009). [CrossRef]

13. V. Toronov, A. Webb, J. H. Choi, M. Wolf, L. Safonova, U. Wolf, and E. Gratton, “Study of local cerebral hemodynamics by frequency-domain near-infrared spectroscopy and correlation with simultaneously acquired functional magnetic resonance imaging,” Opt. Express 9(8), 417–427 (2001). [CrossRef]

14. G. A. Perkins, A. T. Eggebrecht, and H. Dehghani, “Quantitative evaluation of frequency domain measurements in high density diffuse optical tomography,” J. Biomed. Opt. 26(5), 056001 (2021). [CrossRef]

15. M. Doulgerakis, A. T. Eggebrecht, and H. Dehghani, “High-density functional diffuse optical tomography based on frequency-domain measurements improves image quality and spatial resolution,” Neurophotonics 6(03), 1 (2019). [CrossRef]

16. A. Sassaroli, G. Blaney, and S. Fantini, “Dual-slope method for enhanced depth sensitivity in diffuse optical spectroscopy,” J. Opt. Soc. Am. A 36(10), 1743–1761 (2019). [CrossRef]

17. S. Brandt, Data Analysis (Springer, 1998).

18. S. Fantini, D. Hueber, M. A. Franceschini, E. Gratton, W. Rosenfeld, P. G. Stubblefield, D. Maulik, and M. R. Stankovic, “Non-invasive optical monitoring of the newborn piglet brain using continuous-wave and frequency-domain spectroscopy,” Phys. Med. Biol. 44(6), 1543–1563 (1999). [CrossRef]

19. A. Sassaroli, F. Martelli, and S. Fantini, “Perturbation theory for the diffusion equation by use of the moments of the generalized temporal point spread function. I. Theory,” J. Opt. Soc. Am. A 23(9), 2105–2118 (2006). [CrossRef]

20. G. Blaney, R. Donaldson, S. Mushtak, H. Nguyen, L. Vignale, C. Fernandez, T. Pham, A. Sassaroli, and S. Fantini, “Dual-slope diffuse reflectance instrument for calibration-free broadband spectroscopy,” Appl. Sci. 11(4), 1757 (2021). [CrossRef]

21. H. Gu, H. Lu, F. Q. Ye, E. A. Stein, and Y. Yanga, “Noninvasive quantification of cerebral blood volume in humans during functional activation,” NeuroImage 30(2), 377–387 (2006). [CrossRef]

22. D. Contini, F. Martelli, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. 36(19), 4587–4599 (1997). [CrossRef]

Novel data types for frequency-domain diffuse optical spectroscopy and imaging of tissues: characterization of sensitivity and contrast-to-noise ratio for absorption perturbations

Abstract

1. Introduction

2. Theory

2.1 Definition of new data types

2.2 Calculation of the Jacobians of the data types

2.3 Definition of single-distance (SD), single-slope (SS), and dual-slope (DS) sensitivities

2.4 Contrast-to-noise ratio considerations

3. Results

3.1 Comparison of FD data types and their fourth order moment expansion

3.2 Sensitivity of higher order moments

3.3 Sensitivity of FD novel data types

3.4 Contrast-to-noise ratio of FD novel data types

4. Discussion

5. Conclusions and future studies

Appendix A

Appendix B

Appendix C

Appendix D

Funding

Acknowledgements

Disclosures

Data availability

Supplemental document

References

Supplementary Material (1)

Data availability

Cited By

Figures (14)

Equations (70)

Biomedical Optics Express