Separation of superficial and cerebral hemodynamics using a single distance time-domain NIRS measurement

Alexander Jelzow; Heidrun Wabnitz; Ilias Tachtsidis; Evgeniya Kirilina; Rüdiger Brühl; Rainer Macdonald

doi:10.1364/BOE.5.001465

1. Introduction

Functional near-infrared spectroscopy (fNIRS) is an optical technique for non-invasive studies of brain function. Due to its mobility and applicability at the bedside as well as its cost effectiveness, fNIRS became popular in neuroscientific research. The development of the method, its advantages and potential fields of application have been reviewed in [1,2].

However, due to the fundamental principles of light propagation in tissue, signals measured by fNIRS exhibit a higher sensitivity to absorption changes in superficial tissues compared to absorption changes in the deeper lying brain compartment [3,4]. Therefore light absorption changes resulting from spontaneous as well as task-evoked changes in the scalp perfusion [5–7] can mask cerebral signals and may cause false-positive results in functional studies [8]. Thus, a reliable separation of superficial and cerebral hemodynamic signals is required for better accuracy of fNIRS.

In the past, numerous methods for the separation of superficial and cerebral signals, mostly for continuous-wave (cw) fNIRS, have been suggested. Examples are signal correction based on recording at two [9,10] or multiple [7] source-detector separations, application of principle or independent component analysis [11,12], adaptive filtering [13,14], general linear modeling (GLM) [6], wavelet coherence analysis [15] as well as methods exploring temporal correlation between oxy- and deoxyhemoglobin concentration changes induced by neuronal activation [16,17]. The latter approaches rely on assumptions about dynamic properties of extra- and intracerebral time series data.

Alternately, the time-domain (td) fNIRS technique can be applied in order to improve depth selectivity of fNIRS measurements. Using information about the time of flight of photons through tissue [18] in combination with suitable models and algorithms td-NIRS has the potential to achieve depth selectivity even from a measurement at a single source-detector separation and independently of the properties of the time series. From distributions of time of flight of photons (DTOF) measured by td-NIRS, depth-dependent absorption changes can be determined using the time-resolved Beer-Lambert law [19–21] or the concept of time-dependent mean partial pathlenghts (TMPPs) [22–24]. These approaches differ with respect to the calculation of the model related parameters, e.g., diffusion approximation based models vs. Monte-Carlo simulations, and the handling of the instrumental response function (IRF). However, a broad IRF was shown to hamper the quantification of the reconstructed absorption changes [22,24]. In addition, the poor knowledge of the optical properties of the head can introduce systematic errors to model related parameters.

In this work we present a new approach to address these challenges. We applied a data analysis method for time-domain fNIRS signals based on moments of DTOFs which strongly simplifies the handling of the IRF [3,25]. Our approach to separate extracranial and cerebral absorption changes is based on the work by Liebert et al. [3] which was adapted by Kacprzak et al. [26,27]. In comparison to the previously described approaches we introduce the following important modifications. First, a method is introduced to correct for systematic errors which result during the calculation of moments from experimental data due to the influence of the IRF and the finite integration range. Second, we present an algorithm for the reconstruction of absorption changes based on an approximation of the human head by a layered geometry, a forward calculation of model related parameters and thickness parameters derived from an anatomical magnetic resonance imaging (MRI) measurement. Both, the model and the thickness parameters are derived from individual experimental data. The approach is validated on a two-layered phantom. In addition, the retrieved changes in absorption for the two-layered geometry are compared to those obtained by a model assuming homogeneous absorption changes and exploiting the specific intrinsic depth sensitivity of the various moments only.

In order to demonstrate the potential of the proposed approach we used it for the analysis of functional td-NIRS data obtained on the forehead during the execution of a continuous performance task. This task was demonstrated to induce strong task-evoked hemodynamic changes in the skin thereby degrading the sensitivity of fNIRS to cerebral signals [6]. Applying the new approach we were able to separate the superficial and cerebral components without the use of skin-fMRI [6] or additional recordings of systemic physiological signals [15].

2. Experimental

The time-domain NIRS instrumentation used for all measurements presented in this paper was described previously [28,29]. Briefly, the device was equipped with three diode laser heads emitting picoseconds light pulses at wavelengths of 687 nm, 797 nm and 826 nm (Sepia, Picoquant, Germany). Light was collected by four fiber bundles with a diameter of 4 mm and a numerical aperture (NA) of 0.54. Photon detection was performed using four fast photomultipliers (R7400U-02, Hamamatsu Photonics, Japan) each connected to a separate module for time-correlated single photon counting (TCSPC, Becker & Hickl SPC-134, Germany). The td-NIRS instrument was designed for functional in vivo measurements but can also be used for phantom experiments. All IRF measurements were performed as recommended in [30].

A two-layered liquid phantom was used to test the methods for the reconstruction of absorption changes described in the next section. The phantom consisted of a 1 cm thick upper layer and a 6 cm thick lower layer which can effectively be assumed as infinite. Both layers were separated by a thin (30 µm) slightly scattering Mylar foil. The separation between the source fiber and the detector fiber bundle was 3 cm.

A stock solution for both layers was prepared from Intralipid, ink and water. We used Intralipid and pre-diluted ink which were accurately characterized at the University of Florence (Italy) and used in multi-laboratory studies within the European project “nEUROPt” [31,32]. Based on the known scattering and absorption coefficients of the components, the stock solution (S0) for the phantoms was prepared to have optical properties µ_s´ = 10 cm⁻¹ and µ_a = 0.1 cm⁻¹. The characterization and preparation procedures employed are described in detail in [33].

In the first part of the phantom experiment, the absorption of the upper layer was increased while the optical properties of the lower layer remained unchanged. To add absorber without inducing a change in scattering, another solution (S1) was prepared from a mixture of dilute ink, water and intralipid, with µ_s´ = 10 cm⁻¹ and high absorption. The absorption of the upper layer µ_a,₁ was increased from 0.1 cm⁻¹ to about 0.2 cm⁻¹ in 11 steps. Solution S1 was filled in syringes and a certain amount of drops was added to the upper layer. In the second part of the experiment the liquid in the upper layer was removed and replaced with the stock solution S0 so that its absorption coefficient was reset to 0.1 cm⁻¹. Then the absorption of the lower layer µ_a,₂ was changed in the same way as it was done before for the upper layer. Another solution (S2) with a higher absorption but the same scattering as S0 was added to the second layer. The absorption µ_a,₂ was changed from 0.1 cm⁻¹ to about 0.2 cm⁻¹ in 11 steps while the optical properties of the upper layer remained unchanged. In both parts of the experiment, the first two steps of absorption changes were about 0.005 cm⁻¹ and the remaining steps about 0.01 cm⁻¹.

To determine the exact amount of liquid added in a certain step each syringe was weighed by a precision balance before and after adding solution S1 or S2 to the phantom. The mass difference, dilution factors and the known intrinsic optical properties of the ink solution and Intralipid were used to calculate the experimentally prepared absorption change. These changes were treated as the “true” values of µ_a.

For the phantom measurement only one detector of the td-NIRS instrument and the laser emitting at 797 nm were used. A series of 100 DTOFs with a collection time of 1 s per DTOF was recorded for each step. An IRF measurement was performed after each change of absorption to track potential drifts of the IRF in time.

A functional in vivo measurement on the forehead was used to test the reconstruction methods described in the next section. This in vivo study was already presented in [6,15]. Briefly, the td-NIRS instrument was used to measure hemodynamic changes on four locations on the forehead of 15 subjects (5 female/10 male, mean age (35 ± 7) y). The laser light was split into two source fibers. One source fiber and two detection channels were placed on each hemisphere. The source-detector separation was 3 cm. A time series of DTOFs was acquired with a collection time of 49 ms per DTOF and 20 DTOFs per second. The mean photon count rate was about 2.7⋅10⁶ s⁻¹ per TCSPC module and 0.9⋅10⁶ s⁻¹ per wavelength. Anatomical T₁-weighted MR-images with isotropic spatial resolution of 1 mm, providing high contrast between skin, skull and brain tissue were recorded for all subjects. The subjects performed a series of nine blocks of a continuous performance task (CPT) which induces brain activation in the bilateral Brodmann Area 10. Two types of tasks were used alternately and formed a single block: (a) a 34.15 s long control task (word-CPT) and (b) a 34.15 s long semantic task (sem-CPT) that is more challenging. Both tasks were followed by 31.15 s of rest. The sem-CPT was designed to be more challenging than the word-CPT. A detailed description of the tasks can be found in [6].

3. Methods

In the present work we focus on the td-NIRS data analysis based on moments of a DTOF. In the following, we briefly review some basics of the method, present a novel approach for the correction of systematic errors during the calculation of moments and then describe the reconstruction of absorption changes from changes of moments based on a homogeneous and layered structure of a turbid medium.

Moments of DTOFs

In time-domain NIRS the result of a measurement is a distribution of time-of-flight of photons (DTOF). The total photon count is the integral over the DTOF N(t):

N_{T} (N) = \int N (t) d t

N_T is proportional to the detected light intensity. The n-th normalized moment m_n of a DTOF N(t) is defined as:

m_{n} (N) = \frac{1}{N_{T}} \int t^{n} N (t) d t

The first moment m₁ is the mean time of flight. The n-th centralized normalized moment m_n,C is defined as:

m_{n, C} (N) = \frac{1}{N_{T}} \int {(t - m_{1})}^{n} N (t) d t

The second centralized moment is the variance of a DTOF m_2,C = V. It should be noted that during analysis of experimental DTOFs finite sum representations of Eqs. (1)-(3) are used rather than integrals.

A measured DTOF is a convolution of the medium response R(t) with the IRF I(t):

N (t) = R (t) \otimes I (t)

When considering moments of N(t), R(t), and I(t), the convolution-deconvolution problem according to Eq. (4) is reduced to a summation, i.e. the mean time of flight m₁, the variance V and the third centralized moment m₃_,C of N(t) are calculated as the sum of the corresponding moments of R(t) and I(t). Note that this is not true for higher order moments [34]. A direct consequence of this summation rule is that changes in moments of N(t) and R(t) are identical and do not depend on the IRF. Taking the variance as an example we can write ΔV(R) = ΔV(N) = V – V_b, where V_b denotes a reference value of the variance, i.e. an adequately selected baseline. In addition to the changes in higher order moments we define the change in attenuation A based on N_T: ΔA = – ln(N_T / N_T,b), where N_T,b denotes the total photon count of an appropriate baseline. Below we will refer to the quantities ΔN_T, ΔA, Δm₁ and ΔV as the changes in moments.

In the ideal case the integration range in Eqs. (1)-(3) is infinite. In practice, however, we calculate sums over a finite range of time corresponding to a truncation of the integrated function. As a consequence values of moments are systematically too small due to this truncation, the error of which in general increases with the moment’s order [25]. Previously, it was found that this systematic error can be neglected with respect to the lower integration limit but strongly depends on the upper integration limit [25].

The integration limits are typically fixed in a relative manner, i.e. as a percentage of the maximum photon count N_max of a DTOF (cf [25].). In particular, the relative upper integration limit is defined as L_U = N_j / N_max, where N_j is the photon count in the j-th channel located on the tail of the DTOF. The time t_j assigned to this channel is the upper integration limit. Typical values for L_U range from few percent (stronger truncation) to a few per mille (weak truncation). Furthermore, when absolute values of moments are corrected for the influence of the IRF typically the same value of L_U is used for N(t) and I(t). In practice this leads to additional systematic errors which can be substantial if afterpeaks are present in the IRF.

Here we present a procedure which helps to identify the systematic influence as a function of L_U and provides an approximate correction of the systematic errors for the moments m₁, V and m_3,C. The key idea is to approximate the unknown response of the medium by a similar model function S(t) the moments of which would be known exactly. Within the context of this paper we used the reflectance of the homogeneous semi-infinite medium [25,35]. This simple model is frequently used when describing light propagation in the adult human head.

A schematic of the procedure is shown in Fig. 1. In the beginning, measured data is used to calculate initial values of moments m₁ and V which are affected by the systematic error discussed above. These values can be converted to optical properties [25], and from these, the known source-detector separation r_sd and the refractive index n_in an analytical reflectance curve S(t) is calculated. Then S(t) is convolved with the experimental IRF. Using the resulting curve N’(t) and the IRF I(t), the moments of S(t) are calculated for the initial relative integration limits L_L and L_U applied to S(t) as well as I(t). These values of moments include systematic errors due to the finite integration range and the influence of the IRF. However, by comparison with the true values of moments of S(t) estimates of these systematic errors can be obtained in terms of ratios of retrieved (with finite L_L and L_U) and true values for each type of moments. Subsequently, the initial values of the moments retrieved from N(t) can be corrected by the corresponding ratios. The whole correction process can be repeated iteratively to further improve the estimate of the optical properties used to simulate S(t). However, for typical conditions of our experiments and L_U < 0.05 we found N(t) and N’(t) to become very similar already after a single cycle.

Fig. 1 Workflow of the correction procedure for the calculation of moments for a specific set of relative limits L_L and L_U. Symbols are explained in the text.

Download Full Size | PDF

The performance of the correction was first tested on simulations for various configurations of a heterogeneous two-layered medium. For this test an analytical solution of the diffusion equation for the time-resolved light propagation in a layered medium was used [36]. We used the implementation kindly provided by the authors. Here we present an example with the following parameters: µ_s´ = 10 cm⁻¹ in both layers, r_sd = 3 cm⁻¹, n_in = 1.33, µ_a_,1 = 0.1 cm⁻¹, µ_a_,2 = 0.2 cm⁻¹, thickness of the upper layer d = 0.7 cm. The simulated curve R(t) was convolved with an experimental IRF measured by the td-NIRS instrument. This convolved curve was treated as a measured DTOF N(t) and the correction procedure was applied as described above for multiple values of L_U and L_L = 0.03. Figure 2 shows the results of the correction. The left panel contains the experimentally obtained IRF as well as the simulated and convolved DTOFs. The panels on the right show the uncorrected and corrected values of three moments as a function of L_U.

Fig. 2 Performance of the correction procedure. Left: Simulated DTOF R(t), experimental IRF I(t) and the result of the convolution of the two. Right: values of moments m₁, V and m_3,C as a function of the relative upper integration limit L_U. Moments were calculated without (blue) and with (orange) application of the correction procedure illustrated in Fig. 1 for multiple values of L_U. The horizontal cyan lines correspond to the true values of moments calculated directly from the simulated R(t) without clipping.

Download Full Size | PDF

The correction procedure results in more accurate values of moments for nearly any value of L_U. Indeed, an inspection of the course of the uncorrected moments as a function of L_U in Fig. 2 exhibits peaks and high variability. Without prior knowledge about the true value it is even difficult to decide which L_U value is acceptable. In contrast, this is much easier if the correction algorithm is applied. For L_U values below 1% the corrected values of moments are nearly independent of L_U. In the example shown here the effect of correction for m₁ is remarkable: the remaining deviation from the true value is less than 3 ps (0.4%). For the higher order moments the relative deviation is larger but still acceptable. It should be noted that even the peak deviation around L_U = 0.03 - caused by a shoulder in the IRF - is considerably reduced. Nevertheless such L_U ranges should be avoided, in particular when calculating m_3,C. The application of the correction procedure to other configurations, i.e. with other values of optical properties and layer thickness, of the two-layered medium produced similar results to those shown in Fig. 2.

In the past others admonished that moments of order two or higher must be used with great care because of the systematic errors [25,34]. Here we demonstrated that m_3,C can be calculated with an acceptable deviation of less than 10%. This is an important result for the application of the reconstruction of homogeneous absorption changes as discussed in the next section.

Reconstruction of homogenous absorption changes (RHAC)

We first model the head as semi-infinite medium with homogeneous absorption changes Δµ_a, i.e. equal absorption changes in all compartments of the medium. For small Δµ_a we can assume a linear relationship between Δµ_a and a change of a moment ΔM:

Δ M \approx S_{M} Δ μ_{a} = \frac{\partial M}{\partial μ_{a}} Δ μ_{a}

S_M is the sensitivity factor for the change of a moment ΔM due to the absorption change Δµ_a. For the model of homogeneous absorption changes a single moment is enough to calculate the change in absorption. The corresponding sensitivity factors for changes of attenuation, mean time of flight and variance were calculated as [3,22,37]:

S_{A} = c_{m} m_{1} \approx r_{s d} D_{P F}

S_{m 1} = - c_{m} V

S_{V} = - c_{m} m_{3, C}

Here c_m denotes the speed of light in the medium which is assumed to have an equal index of refraction in all compartments. Note that sensitivity factors of a moment basically depend on the moment of the next higher order. Also note that S_A is related to the differential pathlength factor D_PF which is widely used in continuous-wave fNIRS.

In time-domain fNIRS sensitivity factors can be calculated using Eqs. (6)-(8) and the values of moments obtained directly from the measurement without employing any additional model of light propagation. Subsequently using Eq. (5) the corresponding absorption change can be calculated from changes of individual moments. In general, absorption changes reconstructed from different moments will have the same value as long as Δµ_a is homogeneous. This, however, is usually not the case in fNIRS where absorption changes are typically more localized. Consequently, the absorption change calculated in this way will be underestimated – an observation which is known as the partial volume effect [38,39].

In the latter case of inhomogeneous absorption changes the use of Eqs. (5)-(8) will usually result in different values of Δµ_a derived from different moments. This behavior reflects differences in the depth selectivity of the different moments [3,40,41]. For example, changes in variance are comparably more sensitive to absorption changes in the brain than changes in attenuation. If the RHAC method is applied to functional in vivo data a trained operator can compare time traces of hemoglobin concentration changes obtained from different moments and qualitatively distinguish between a superficial or a cerebral origin of these signals. However, in [6] we showed that also variance based signals can be seriously disturbed by task-evoked hemoglobin concentration changes in the scalp. Despite these drawbacks the RHAC approach can be useful. It is based solely on data available from a typical td-fNIRS measurement and does not require the use of sophisticated models of light propagation in tissue and assumptions regarding geometry and optical properties.

The upper panels of Fig. 3 show distributions of values of those moments which are relevant for the calculation of sensitivity factors by Eqs. (6)-(8) as obtained in vivo on the forehead of 15 healthy adult subjects. The mean values of m₁ obtained here coincide within ± 2σ with data published in [42,43] and are always lower (cf. Table 1). To our best knowledge values for V and m₃_,C are reported for the first time. Unlike data for m₁ these values show only negligible spectral dependency. On the other hand, they exhibit a much wider variability. Mean values and the corresponding coefficients of variation are summarized in Table 1.

Fig. 3 Distribution of values of moments and homogeneous background optical properties obtained in vivo on the forehead of 15 adults and four optodes (60 samples). Upper panels: Distribution of values of the moments m₁, V and m₃_,C for r_sd = 3 cm obtained at three wavelengths. Blue dots - individual values, gray area - interpolated histograms, black and red lines - box plots (red horizontal bars - median position, wide black bars - 1st and 3rd quartile, short black bars - extreme values), white square - mean value. Lower panels: distribution of individual homogeneous background optical properties µ_a and µ_s´ calculated from m₁ and V at the three wavelengths. White squares mark the mean values. Mean values of moments are summarized in Table 1.

Download Full Size | PDF

Table 1. Mean values of moments m₁, V and m₃_,C as well as homogeneous optical properties µ_a and µ_s´ derived from m₁ and V. Corresponding coefficients of variation (CV) are given in parentheses (60 samples). For m₁ values from literature are given for comparison.

View Table

We further used the values of m₁ and V to calculate individual homogeneous background optical properties µ_a and µ_s´, according to the related equations derived for the homogeneous semi-infinite medium [25]. It should be noted that these equations rely on the solution of the diffusion equation for zero boundary conditions while the extrapolated boundary condition is generally accepted as being more adequate [45]. However, the major difference occurs in the amplitude of the DTOF which does not influence m₁ and V. Optical properties derived in this way are presented in the lower panels of Fig. 3 and used in the next section.

In reality the head is a heterogeneous medium. However, it is difficult to reconstruct optical properties of all compartments in vivo. Therefore, the approximation of the head by a homogeneous SIM model is a reasonable and frequently used compromise. We derive properties which are most likely not an exact match of those of any tissue type included. But they are the best “explanation” of the shape of the DTOF based on two experimentally obtained parameters, i.e. m₁ and V.

Reconstruction of absorption changes based on a layered geometry (RLAC)

The reconstruction of absorption changes presented in this Section is based on a layered structure of the turbid medium and aims to improve the quantification and the separation of superficial and cerebral hemodynamics in fNIRS compared to the homogeneous approach described in the previous section. The key idea is to model the head as a layered semi-infinite medium with individual (subject-specific) homogeneous background optical properties (BOP), i.e. the absorption coefficient µ_a, the reduced scattering coefficient µ_s´ and the refractive index. A model of light propagation is then applied to this geometry to calculate sensitivity factors for absorption changes in layers.

For small changes of absorption the change of a moment is given by the linear approximation:

Δ M \approx \sum_{j} S_{M, j} Δ μ_{a, j}

Here Δµ_a,j denotes the absorption change in the j-th layer and S_M,,j is the sensitivity factor of moment M with respect to Δµ_a,j:

S_{M, j} = \frac{\partial M}{\partial μ_{a, j}} \approx \frac{Δ M}{Δ μ_{a, j}}

In order to calculate the sensitivity factors S_M,,j we used a forward model for time-resolved light propagation in a layered turbid cylinder [36]. The cylinder was subdivided into 25 plane layers with a thickness of 2 mm each and with the deepest layer of 10 mm thickness (60 mm in total). The radius of the cylinder was set to 60 mm and the source-detector separation to 30 mm. Altogether this configuration of the cylinder is effectively identical to a semi-infinite medium. The BOP of the layers were varied in the range from 0.04 cm⁻¹ to 0.35 cm⁻¹ in steps of 0.01 cm⁻¹ and from 4 cm⁻¹ to 19.5 cm⁻¹ in steps of 0.5 cm⁻¹ for µ_a and µ_s´, respectively. The refractive index was fixed to 1.4 for all layers. To obtain the sensitivity factor of the j-th layer the absorption coefficient in the corresponding layer was changed by ± 2% and a DTOF was simulated in each case. Using these DTOFs the changes in attenuation, mean time of flight and variance were calculated. The sensitivity factors were obtained according to Eq. (10). The results of this calculation procedure for the whole range of the background optical properties were used to create a lookup table (LUT) of sensitivity factors. This LUT is then used to obtain individual sensitivity factors based on individual BOP by interpolation for each layer, moment’s type and wavelength. The LUT approach is faster than on-demand calculation of sensitivity factors (SF) and accurate enough if the LUT is adequately sampled with respect to BOP.

The sensitivity factors for a change in µ_a of a 2 mm thick layer as a function of BOP are visualized in Fig. 4. As expected the depth of the peak sensitivity is increasing from A to m₁ and to V [3]. The results demonstrate that the sensitivity factors strongly depend on the BOP of the medium. This dependence is more pronounced for changes in the mean time of flight and variance than for attenuation changes. All depth sensitivity profiles shift towards medium surface with increasing µ_a. Furthermore, increasing absorption leads to a strong decrease in the overall sensitivity of m₁ and V.

Fig. 4 Sensitivity factors (SF) for absorption changes obtained from a simulation. Plots in the upper, middle and lower rows refer to SF for changes in attenuation A, mean time of flight m₁ and variance V, respectively. Left three columns: color maps of SF as a function of the layer number j and µ_s´. The units for the color bars are printed on the far-right. The plots refer to fixed values of µ_a, i.e. 0.05 cm⁻¹, 0.1 cm⁻¹ and 0.2 cm⁻¹ (left to right). Thin black lines show the position of the extrema of the (Z-dependent) sensitivity profiles. Right column: selected depth sensitivity profiles shown for µ_s´ = 10 cm⁻¹ (white dashed lines on the color maps). The colors of the lines (red, green, blue) refer to the µ_a values of the three columns on the left. Gray lines show the shift of the median (50%) and the 90% percentile.

Download Full Size | PDF

With increasing scattering the position of the maximum sensitivity (black lines on color maps) of m₁ and V shifts towards surface while it remains nearly constant for A. For low scattering the difference in depth selectivity between ΔA on the one side and Δm₁ and ΔV on the other is maximized. These results confirm the need to consider individual optical properties in order to obtain more reliable sensitivity factors for the m₁ and V.

In principle, using the LUT of the sensitivity factors it is possible to reconstruct absorption changes in a multi-layered medium. However, it is not possible to perform the reconstruction for all 25 layers from changes in only three moments because the set of linear equations to be solved would be highly underdetermined. Thus, in our model we assume that changes of absorption related to the functional activation can appear in two layers only, namely the superficial and cerebral compartments.

We further assume that there are no functional absorption changes in the intermediate skull compartment. Effectively, the medium is subdivided into three layers but changes in absorption are allowed in two layers only. Thus the two “active” layers are not contiguous but separated by a “passive” intermediate layer. These assumptions are physiologically reasonable and reduce the crosstalk in the reconstruction. It should be noted that Kacprzak et al. [27] used a similar approach but reconstructed absorption changes in all three layers and disregarded the intermediate layer in their further analysis.

The sensitivity factors S_M_,_S and S_M_,_B of the superficial and brain compartments are obtained as follows. First, the depth sensitivity profile S_M_,_j for 25 layers is obtained from the LUT. This profile is then interpolated to obtain the quasi continuous depth sensitivity profile S_M(Z) which is a function of the depth Z. The superficial SF S_M_,_S and the brain SF S_M_,_B are calculated as integrals over S_M(Z) in certain ranges, e.g. from 0 to z₁ for S_M_,_S and from z₂ to z₃ for S_M_,_B. Altogether, a change of the moment is then modeled as:

Δ M \approx Δ μ_{a, S} \int_{0}^{z_{1}} S_{M} (Z) d Z + Δ μ_{a, B} \int_{z_{2}}^{z_{3}} S_{M} (Z) d Z = Δ μ_{a, S} S_{M, S} + Δ μ_{a, B} S_{M, B}

In an anatomical model of the head z₁ reflects the thickness of the scalp, z₂ the distance to the gray matter and z₃ can be set to infinity because the depth profiles S_M(Z) are quickly approaching zero for typical BOP observed in vivo. In order to consider changes in all measured moments Eq. (11) can be written in matrix form:

Δ M = (\begin{array}{l} Δ A \\ Δ m_{1} \\ Δ V \end{array}) = [\begin{array}{l} \begin{matrix} S_{A, S} & S_{A, B} \end{matrix} \\ \begin{matrix} S_{m_{1}, S} & S_{m_{1}, B} \end{matrix} \\ \begin{matrix} S_{V, S} & S_{V, B} \end{matrix} \end{array}] \cdot (\begin{array}{l} Δ µ_{a, S} \\ Δ µ_{a, B} \end{array}) = S \cdot Δ µ_{a}

This set of linear equations can be solved in a least-squares sense. We follow the approach outlined by Liebert et al. [46] but adapt it to our specific conditions. Since the changes in moments are calculated from the same DTOF the covariance of the moments must be considered. The solution of Eq. (12) may be obtained using the weighted least-squares approach [47]:

Δ µ_{a} = {(S^{T} K^{- 1} S)}^{- 1} S^{T} K^{- 1} Δ M

Here K is the covariance matrix of ΔM and in our case determined by the photon noise as:

K = [\begin{matrix} \frac{1}{N_{T}} & 0 & 0 \\ 0 & \frac{V}{N_{T}} & \frac{m_{3, C}}{N_{T}} \\ 0 & \frac{m_{3, C}}{N_{T}} & \frac{m_{4, C} - V^{2}}{N_{T}} \end{matrix}]

Note that the values of moments in this matrix refer to the moments of the DTOF N(t) as measured, i.e. including the influence of the IRF.

The approach for the reconstruction of absorption changes in a layered medium described above shares some similarities with approaches published before. Steinbrink et al. [22,23] used a layered structure for absorption changes, the concept of mean partial pathlengths, the analysis of the full DTOF and calculated sensitivity factors using a Monte-Carlo simulation. This powerful approach, however, employed only a single set of homogeneous background optical properties and required either a narrow IRF or a precise deconvolution. Liebert et al. [3,46] circumvent the difficulties with the IRF and the deconvolution by using moments of DTOFs, as mentioned above. However, they still used a Monte-Carlo simulation and, in addition, a multi-distance time-domain measurement. Finally, Kacprzak et al. [26,27] employed a td-NIRS measurement at a single source-detector separation, moments of DTOFs and calculated sensitivity factors using a perturbation approach. As mentioned by the authors this method yields significant errors in the vicinity of light sources and detectors. Furthermore, the authors used a strict subdivision of the medium into upper, intermediate and lower layers with a fixed thickness of 9 mm, 6 mm and 15 mm, respectively, and a reconstruction of absorption changes in all three layers. In our opinion this approach is prone to cross talk in the reconstruction because of the strong similarity of the sensitivity profiles of m₁ and V [48] and the resulting ill-posed inverse problem.

In our approach we also use a td-fNIRS measurement at a single source-detector separation, moments of DTOF and a layered structure of the medium. Sensitivity factors are calculated using a forward simulation of light propagation in tissue which is less time-consuming than Monte-Carlo calculations. The use of a LUT further accelerates the calculation of sensitivity factors and allows considering individual background optical properties. The subdivision of the layered medium according to the thickness of scalp and skull obtained from individual MR images is more adequate than fixed compartments. The integration of the interpolated sensitivity profiles in Eq. (11) avoids round up errors and reduces systematic errors due to thickness inaccuracies in the reconstruction [24].

Calculation of hemoglobin concentration changes

The absorption changes obtained at each wavelength and for each compartment using either the RHAC or the RLAC approach and the concentration changes of oxy- and deoxyhemoglobin, ΔC_HbO and ΔC_HbR, respectively, are related by:

Δ μ_{a} (λ) = (\begin{matrix} Δ μ_{a} (λ_{1}) \\ ⋮ \\ Δ μ_{a} (λ_{i}) \end{matrix}) = \ln (10) \cdot ε \cdot (\begin{array}{l} Δ C_{H b O} \\ Δ C_{H b R} \end{array}) = \ln (10) \cdot ε \cdot Δ C

Here ε denotes the matrix of specific molar absorbance coefficients. Throughout this work we used the hemoglobin absorption spectra measured by M. Cope [49,50]. Equations (14) and (15) were solved in a least squares sense using the lscov and mrdivide routines provided by MATLAB (Mathworks Inc.), respectively.

4. Validation of methods on a two-layer phantom experiment

The phantom measurements described in Section 2 were analyzed as follows. 100 DTOFs recorded for each absorption step were summed up. From these binned DTOFs and the corresponding IRFs moments were calculated. The results are shown in the upper panels of Fig. 5. In general, with increasing absorption in either layer the change in attenuation increases while m₁ and V decrease. We also compared the measured moments to the theoretically predicted values (not shown). These were calculated using the forward solution of light propagation in layered media [36] and the configuration of the two layers. We found a good agreement for the absolute values of all moments and their dependence on µ_a,₁ and µ_a,₂.

Fig. 5 Results of the measurement on the two-layered phantom. Upper panels: measured courses of changes in attenuation ΔA, mean time of flight m₁ and variance V while µ_a is changing in the upper (left) and lower (right) layer. Note the different color scale for the different measurands. Gray lines refer to cubic fits. Lower panels: changes in the absorption retrieved using reconstruction approaches based on a model of homogeneous (RHAC) and layered (RLAC) absorption changes. Gray lines refer to the true values of Δµ_a.

Download Full Size | PDF

In order to apply the RHAC and RLAC procedures to the data obtained from the phantom experiment small absorption changes are needed in order not to violate the linearization in Eqs. (9)-(10). Therefore, first the courses of the attenuation A, mean time of flight m₁ variance V as a function of µ_a,₁ or µ_a,₂ were fitted by a polynomial of third order. Second, using the data from the cubic fit it was possible to interpolate the course of moments and obtain small changes of moments corresponding to small changes of µ_a. The changes of moments for this procedure were calculated as a difference of adjacent points. The increment for interpolation was set to Δµ_a,_nom = 0.0056 cm⁻¹ (nominal absorption change) which corresponds to a 2.5% to 5.6% relative change in µ_a depending on the absolute value of µ_a,₁ or µ_a,₂ In addition, the fit helps to compensate for fluctuations in the courses of m₁ and V for µ_a,1 between 0.15 cm⁻¹ and 0.19 cm⁻¹, most likely due to insufficient mixing of the liquid in the upper layer.

Both reconstruction methods, i.e. RHAC and RLAC, were applied in such a way that the sensitivity factors were derived individually for each step in absorption. Thus, in general, reconstruction was applied to a state of the phantom where BOP are heterogeneous. Only the first point (leftmost point in the graphs of Fig. 5) describes a homogeneous situation with µ_a,1 = µ_a,2 = 0.1 cm⁻¹.

Reconstruction results achieved by both methods are shown in the lower panels of Fig. 5. The RHAC procedure yields underestimated Δµ_a values (lines with symbols) for all three moments and both layers. If µ_a,₁ is changed Δµ_a,_Δ_A produces the closest result to the true Δµ_a,₁; if µ_a,₂ is changed Δµ_a,_Δ_m₁ and Δµ_a,_Δ_V perform better than Δµ_a,_Δ_A but the retrieval still remains poor. This behavior to some extent reflects differences in the depth sensitivity of moments and the partial volume effect.

The RLAC method yields better results and it can separate absorption changes of both layers (red and blue lines in the lower panels of Fig. 5). In both cases the retrieved Δµ_a,₁ and Δµ_a,₂ vary around the true values in a range of about ± 0.001 cm⁻¹. The dependence of the retrieved values on the absolute absorption coefficients is moderate. Some modulations can be seen. The origin of these deviations is not completely clarified. Potentially, they result from the interplay of the fitting of the moment’s data, the calculation of sensitivity factors and the crosstalk during the reconstruction. Taking these deviations into account the accuracy of the layered reconstruction can be estimated to be better than ± 0.001 cm⁻¹ ( ± 18% if normalized to the nominal change).

The data obtained on the phantom can also be used to study the influence of mismatching parameters on the results of the reconstruction by the RLAC method. First, the RLAC method was applied assuming a wrong thickness of the upper layer. We found that the reconstructed Δµ_a,₁ and Δµ_a,₂ are shifted with respect to the reconstructed courses in Fig. 5 in different ways. If the thickness is too low and µ_a,₁ is increasing, both Δµ_a,₁ and Δµ_a,₂ are overestimated while a too high thickness leads to an underestimation. This behavior is reversed if µ_a,₂ is increased. The shift was approximately linear in the thickness range from 8 mm to 12 mm. On average and compared to the nominal absorption change Δµ_a,_nom the reconstructed Δµ_a,₁ and Δµ_a,₂ were off by 4% and 12% per mm of thickness change, respectively. The observed shift leads to a crosstalk between reconstructed absorption changes in both layers. In the in vivo case a non-zero Δµ_a,₂ caused by a change of µ_a,₁ might be wrongly interpreted as false positive brain activation.

Second, we applied the RLAC method using exemplary non-adapted homogeneous BOP. To this end µ_a and µ_s´ were fixed to 0.1 cm⁻¹ and 10 cm⁻¹, respectively. We found minor effects on the reconstructed values of Δµ_a,₁ and larger deviations for Δµ_a,₂. In the worst case only 30% of Δµ_a,_nom were retrieved for Δµ_a,₂ at µ_a,₂ ≈0.2 cm⁻¹. This example demonstrates that by taking into account adapted homogeneous BOP the RLAC method achieves a more reliable quantification of deep absorption changes.

5. Application to in vivo data

Both reconstruction methods were applied to in vivo functional activation data obtained on the forehead of adult subjects in a combined fNIRS/fMRI study [6]. There we found that the CPT induces strong task-evoked hemodynamic changes localized in the veins of the forehead skin. These changes were time-locked to the task and their time courses differed from the cerebral functional signals.

Exemplary results of the application to data from one subject and two detectors are shown in Fig. 6. The individual thickness parameters as required in Eq. (11) were derived from segmented anatomical MR images as z₁ = 7.4 mm and z₂ = 12.8 mm. The depth z₁ was obtained as the averaged distance between the scalp surface and the skull at the position of fNIRS sources and detectors. Analogously, z₂ was the average distance from the scalp surface to the gray matter. The time-series of moments’ changes were first processed by a low-pass filter (cutoff frequency 0.3 Hz) and a high-pass filter (cutoff frequency 0.004 Hz). Absorption changes were calculated using the RHAC and RLAC approaches and then converted to hemoglobin concentration changes using Eq. (15). Finally, the data of the nine repetitions of the task was block averaged. In order to identify brain activation we qualitatively compared the time course of block averages with the group averaged BOLD signal shown in Fig. 2(b) in [6]. It exhibits a sharp increase within the first five seconds after the onset of stimulation followed by a plateau phase and a slow decrease after the end of stimulation. Moreover, concentration changes are expected to be positive for HbO and negative for HbR.

Fig. 6 Application of reconstruction methods to in vivo data from a single subject and two channels. The periods of the first task (T1, word-CPT), rest (R) and the second task (T2, sem-CPT) are separated by dashed vertical lines. Left: concentration changes in oxy- (red) and deoxyhemoglobin (blue) obtained from changes in ΔA, Δm₁ and ΔV and using the RHAC method. Right: concentration changes obtained from the same in vivo data as before but using the RLAC. Results are reported for the upper and lower layers. Shadowed areas around the lines depict the standard error of mean obtained from the nine repetitions of the tasks.

Download Full Size | PDF

Concentration changes of oxygenated hemoglobin obtained using the RHAC method from all three moments ΔA, Δm₁ and ΔV show an increase during the performance of the tasks while deoxyhemoglobin exhibits a decrease. However, the temporal behavior of the related time series is different. Concentration changes of HbO obtained from ΔA exhibit a pronounced monotonic increase during the task period and a decrease during rest. Visually this course forms a shape similar to a triangle. The course of HbR also shows such a triangular shape but with a reverse sign. In contrast to the data from ΔA, time courses of ΔC_HbR calculated from Δm₁ and ΔV rather have a shape as it is expected for typical brain activation. This finding reflects the higher depth sensitivity of the higher order moments. Time courses reconstructed by the RHAC method are always a superposition of signals originating from the scalp and the brain, with relative weights depending on the order of the moment.

Using the RLAC method a separation of the signal contributions can be achieved. The right part of Fig. 6 shows such a separation. The concentration changes in the upper layer exhibit a pronounced triangular shape which was found before in the ΔA signals as well as in a skin fMRI measurement [6]. This shape is observed for both hemoglobin species but is more prominent in ΔC_HbO. The ΔC_HbR signal obtained for the lower layer and channel 3 shows clearly identifiable brain activation which is slightly stronger for the sem-CPT task (T2) compared to T1. This observation is compatible with the fact that the sem-CPT requires a higher effort and is expected to lead to stronger brain activation than word-CPT. The ΔC_HbO signals in the upper layer can be attributed to task-evoked changes in forehead skin perfusion [6], with a magnitude also depending on the complexity of the task. The ΔC_HbO signals in the lower layer show additional fast signal fluctuations. These are most likely induced by photon noise. It should be noted that in this example the time courses of the lower-layer signals are rather similar to those of the variance-based signals. This is not surprising if the depth-dependent sensitivity profile of variance matches the anatomical conditions to essentially select the cerebral region. However, often the variance-based signals still contain a superficial contribution, as, e.g., in the case of the group-averaged HbO signals shown in Fig. 4 of [6]. This is an indication that the variance sensitivity profile that depends on the optical properties does not perfectly suppress the superficial influence. In this case the individual moments alone (RHAC method) cannot accomplish the retrieval of the cerebral response.

The amplitudes of the hemoglobin concentration changes derived by the RLAC method are larger than those obtained by RHAC. This is a direct consequence of the partial volume effect because RHAC always underestimates the local absorption changes and thus the hemoglobin concentration changes.

The RLAC approach has a few intrinsic limitations. First, it approximates the head as a semi-infinite medium with homogeneous optical properties and considers absorption changes in plane layers. However, the optical properties are derived individually for each subject. Thus we regard the RLAC approach as a good compromise with respect to feasibility and robustness. Second, we use thickness parameters derived from an anatomical MR scan. If such data is not available these parameters might be reasonably approximated using mean thickness parameters obtained from literature. In this case, however, the real thickness may be different and thus a systematic crosstalk between the layers is expected to appear as discussed in Section 4. Alternatively, the required thicknesses can be measured by transcranial ultrasound instrumentation available in many clinical facilities. In this case the whole approach would still be applicable at the bedside.

Conclusion

We presented a method for separation of extracranial and cerebral hemodynamics in td-fNIRS. Our method (RLAC) employs the reconstruction of absorption changes, and consequently hemoglobin concentration changes, in two layers of a turbid medium using a time-domain measurement at a single-source detector separation. RLAC is based on moments of DTOFs and accounts, to some extent, for the individual variability of optical properties by approximating the head with a homogeneous semi-infinite medium with individual optical properties. RLAC reconstructs absorption changes for each point in a time series independently.

We further presented a method to correct values of moments for systematic errors due to a limited integration range and the resulting influence of the IRF. Using this correction method, even centralized moments of third order can be calculated from measured DTOFs with a deviation of less than 10% and can reliably be used within reconstruction procedures.

The RLAC method was validated on experimental data obtained on a two-layered liquid phantom and applied to functional in vivo data. It was shown that a separation of superficial and deeper lying absorption changes can be achieved using only td-NIRS data measured at a single source-detector separation and additional information from structural MRI. In addition, RLAC results were compared to results obtained by a method reconstructing homogeneous absorption changes (RHAC). RHAC always underestimates absorption changes and thus also the hemoglobin concentration changes.

The RLAC approach has the potential to improve the inter-subject comparability of results obtained in fNIRS studies. The quantification of cerebral hemoglobin concentration changes is enhanced by taking into account individual parameters and by separating interfering superficial contamination from cerebral signals without additional skin-fMRI measurements and recordings of systemic physiological signals.

Acknowledgments

We thank Dr. A. Liemert and Dr. A. Kienle for the software implementation of light propagation in layered turbid media. The authors are thankful to Dr. A. Heine for developing the task used in the functional in vivo measurement. We thank Dr. M. Niessing for running the segmentation procedure of the anatomical MR images. The research leading to these results has received funding from the European Community's Seventh Framework Programme [FP7/2007-2013] under grant agreement n° FP7-HEALTH-F5-2008-201076. Ilias Tachtsidis is supported by the Wellcome Trust (088429/Z/09/Z).

References and links

1. M. Ferrari and V. Quaresima, “A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application,” Neuroimage 63(2), 921–935 (2012). [CrossRef] [PubMed]

2. F. Scholkmann, S. Kleiser, A. J. Metz, R. Zimmermann, J. Mata Pavia, U. Wolf, and M. Wolf, “A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology,” Neuroimage 85(Pt 1), 6–27 (2014). [CrossRef] [PubMed]

3. A. Liebert, H. Wabnitz, J. Steinbrink, H. Obrig, M. Möller, 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] [PubMed]

4. P. G. Al-Rawi, P. Smielewski, and P. J. Kirkpatrick, “Evaluation of a Near-Infrared Spectrometer (NIRO 300) for the Detection of Intracranial Oxygenation Changes in the Adult Head,” Stroke 32(11), 2492–2500 (2001). [CrossRef] [PubMed]

5. T. Takahashi, Y. Takikawa, R. Kawagoe, S. Shibuya, T. Iwano, and S. Kitazawa, “Influence of skin blood flow on near-infrared spectroscopy signals measured on the forehead during a verbal fluency task,” Neuroimage 57(3), 991–1002 (2011). [CrossRef] [PubMed]

6. E. Kirilina, A. Jelzow, A. Heine, M. Niessing, H. Wabnitz, R. Brühl, B. Ittermann, A. M. Jacobs, and I. Tachtsidis, “The physiological origin of task-evoked systemic artefacts in functional near infrared spectroscopy,” Neuroimage 61(1), 70–81 (2012). [CrossRef] [PubMed]

7. N. M. Gregg, B. R. White, B. W. Zeff, A. J. Berger, and J. P. Culver, “Brain specificity of diffuse optical imaging: improvements from superficial signal regression and tomography,” in Front. Neuroenergetics (2010).

8. I. Tachtsidis, T. S. Leung, A. Chopra, P. H. Koh, C. B. Reid, and C. E. Elwell, “False positives in functional near-infrared topography,” Adv. Exp. Med. Biol. 645, 307–314 (2009). [CrossRef] [PubMed]

9. R. B. Saager and A. J. Berger, “Direct characterization and removal of interfering absorption trends in two-layer turbid media,” J. Opt. Soc. Am. A 22(9), 1874–1882 (2005). [CrossRef] [PubMed]

10. R. B. Saager, N. L. Telleri, and A. J. Berger, “Two-detector Corrected Near Infrared Spectroscopy (C-NIRS) detects hemodynamic activation responses more robustly than single-detector NIRS,” Neuroimage 55(4), 1679–1685 (2011). [CrossRef] [PubMed]

11. J. Virtanen, T. Noponen, and P. Meriläinen, “Comparison of principal and independent component analysis in removing extracerebral interference from near-infrared spectroscopy signals,” J. Biomed. Opt. 14(5), 054032 (2009). [CrossRef] [PubMed]

12. T. Funane, H. Atsumori, T. Katura, A. N. Obata, H. Sato, Y. Tanikawa, E. Okada, and M. Kiguchi, “Quantitative evaluation of deep and shallow tissue layers’ contribution to fNIRS signal using multi-distance optodes and independent component analysis,” Neuroimage 85(Pt 1), 150–165 (2014). [CrossRef] [PubMed]

13. Q. Zhang, E. N. Brown, and G. E. Strangman, “Adaptive filtering for global interference cancellation and real-time recovery of evoked brain activity: a Monte Carlo simulation study,” J. Biomed. Opt. 12(4), 044014 (2007). [CrossRef] [PubMed]

14. Q. Zhang, G. E. Strangman, and G. Ganis, “Adaptive filtering to reduce global interference in non-invasive NIRS measures of brain activation: How well and when does it work?” Neuroimage 45(3), 788–794 (2009). [CrossRef] [PubMed]

15. E. Kirilina, N. Yu, A. Jelzow, H. Wabnitz, A. M. Jacobs, and I. Tachtsidis, “Identifying and quantifying main components of physiological noise in functional near infrared spectroscopy on the prefrontal cortex,” Front Hum Neurosci 7, 864 (2013). [PubMed]

16. T. Yamada, S. Umeyama, and K. Matsuda, “Separation of fNIRS Signals into Functional and Systemic Components Based on Differences in Hemodynamic Modalities,” PLoS ONE 7(11), e50271 (2012). [CrossRef] [PubMed]

17. X. Cui, S. Bray, and A. L. Reiss, “Functional Near Infrared Spectroscopy (NIRS) signal improvement based on negative correlation between oxygenated and deoxygenated hemoglobin dynamics,” Neuroimage 49(4), 3039–3046 (2010). [CrossRef] [PubMed]

18. 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(Pt 1), 28–50 (2014). [PubMed]

19. Y. Nomura, O. Hazeki, and M. Tamura, “Exponential Attenuation of Light Along Nonlinear Path Through the Biological Model,” in Oxygen Transport to Tissue XI, K. Rakusan, G. P. Biro, T. K. Goldstick, and Z. Turek, eds., Advances in Experimental Medicine and Biology No. 248 (Springer US, 1989), pp. 77–80.

20. Y. Nomura and M. Tamura, “Picosecond Time of Flight Measurement of Living Tissue: Time Resolved Beer-Lambert Law,” in Oxygen Transport to Tissue XIII, T. K. Goldstick, M. McCabe, and D. J. Maguire, eds., Advances in Experimental Medicine and Biology No. 316 (Springer US, 1992), pp. 131–136.

21. Y. Nomura, O. Hazeki, and M. Tamura, “Relationship between time-resolved and non-time-resolved Beer-Lambert law in turbid media,” Phys. Med. Biol. 42(6), 1009–1022 (1997). [CrossRef] [PubMed]

22. J. Steinbrink, “Near-infrared-spectroscopy on the adult human head with picosecond resolution,” PhD Thesis, FU Berlin (2000).

23. 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] [PubMed]

24. L. Zucchelli, D. Contini, R. Re, A. Torricelli, and L. Spinelli, “Method for the discrimination of superficial and deep absorption variations by time domain fNIRS,” Biomed. Opt. Express 4(12), 2893–2910 (2013). [CrossRef] [PubMed]

25. 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] [PubMed]

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

27. M. Kacprzak, A. Liebert, W. Staszkiewicz, A. Gabrusiewicz, P. Sawosz, G. Madycki, and R. Maniewski, “Application of a time-resolved optical brain imager for monitoring cerebral oxygenation during carotid surgery,” J. Biomed. Opt. 17(1), 016002 (2012). [CrossRef] [PubMed]

28. H. Wabnitz, M. Moeller, A. Liebert, A. Walter, R. Erdmann, O. Raitza, C. Drenckhahn, J. P. Dreier, H. Obrig, J. Steinbrink, and R. Macdonald, “A time-domain NIR brain imager applied in functional stimulation experiments,” Proc. SPIE 5859, 58590H (2005). [CrossRef]

29. H. Wabnitz, M. Moeller, A. Liebert, H. Obrig, J. Steinbrink, and R. Macdonald, “Time-resolved near-infrared spectroscopy and imaging of the adult human brain,” Adv. Exp. Med. Biol. 662, 143–148 (2010). [CrossRef] [PubMed]

30. A. Liebert, H. Wabnitz, D. Grosenick, and R. Macdonald, “Fiber dispersion in time domain measurements compromising the accuracy of determination of optical properties of strongly scattering media,” J. Biomed. Opt. 8(3), 512–516 (2003). [CrossRef] [PubMed]

31. L. Spinelli, M. Botwicz, N. Zolek, M. Kacprzak, D. Milej, A. Liebert, U. Weigel, T. Durduran, F. Foschum, A. Kienle, F. Baribeau, S. Leclair, J.-P. Bouchard, I. Noiseux, P. Gallant, O. Mermut, A. Pifferi, A. Torricelli, R. Cubeddu, H.-C. Ho, M. Mazurenka, H. Wabnitz, K. Klauenberg, O. Bodnar, C. Elster, M. Bénazech-Lavoué, Y. Bérubé-Lauzière, F. Lesage, P. D. Ninni, F. Martelli, and G. Zaccanti, “Inter-Laboratory Comparison of Optical Properties Performed on Intralipid and India Ink,” in Biomedical Optics and 3-D Imaging (Optical Society of America, 2012), p. BW1A.6.

32. H. Wabnitz, A. Jelzow, M. Mazurenka, O. Steinkellner, R. Macdonald, A. Pifferi, A. Torricelli, D. Contini, L. M. G. Zucchelli, L. Spinelli, R. Cubeddu, D. Milej, N. Zolek, M. Kacprzak, P. Sawosz, A. Liebert, S. Magazov, J. C. Hebden, F. Martelli, P. Di Ninni, and G. Zaccanti, “Performance assessment of time-domain optical brain imagers: a multi-laboratory study,” Proc. SPIE 8583, 85830L (2013). [CrossRef]

33. F. Martelli and G. Zaccanti, “Calibration of scattering and absorption properties of a liquid diffusive medium at NIR wavelengths. CW method,” Opt. Express 15(2), 486–500 (2007). [CrossRef] [PubMed]

34. C. A. Laury-Micoulaut, “The n-th centered moment of a multiple convolution and its applications to an intercloud gas model,” Astron. Astrophys. 51, 343–346 (1976).

35. M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. 28(12), 2331–2336 (1989). [CrossRef] [PubMed]

36. A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. II. Layered case,” Opt. Express 18(9), 9266–9279 (2010). [CrossRef] [PubMed]

37. 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] [PubMed]

38. M. Hiraoka, M. Firbank, M. Essenpreis, M. Cope, S. R. Arridge, P. van der Zee, and D. T. Delpy, “A Monte Carlo investigation of optical pathlength in inhomogeneous tissue and its application to near-infrared spectroscopy,” Phys. Med. Biol. 38(12), 1859–1876 (1993). [CrossRef] [PubMed]

39. D. A. Boas, T. Gaudette, G. Strangman, X. Cheng, J. J. A. Marota, and J. B. Mandeville, “The Accuracy of Near Infrared Spectroscopy and Imaging during Focal Changes in Cerebral Hemodynamics,” Neuroimage 13(1), 76–90 (2001). [CrossRef] [PubMed]

40. H. Wabnitz, M. Moeller, A. Liebert, A. Walter, R. Macdonald, H. Obrig, J. Steinbrink, R. Erdmann, and O. Raitza, “A Time-Domain NIR Brain Imager Applied in Functional Stimulation Experiments,” in Photon Migration and Diffuse-Light Imaging II, K. and C. Licha, ed., Proc. SPIE (Optical Society of America, 2005), Vol. 5859, p. WA5.

41. H. Wabnitz, A. Liebert, D. Contini, L. Spinelli, and A. Torricelli, “Depth Selectivity in Time-Domain Optical Brain Imaging Based on Time Windows and Moments of Time-of-Flight Distributions,” in Biomedical Optics, OSA Technical Digest (CD) (Optical Society of America, 2008), p. BMD9.

42. A. Duncan, J. H. Meek, M. Clemence, C. E. Elwell, L. Tyszczuk, M. Cope, and D. T. Delpy, “Optical pathlength measurements on adult head, calf and forearm and the head of the newborn infant using phase resolved optical spectroscopy,” Phys. Med. Biol. 40(2), 295–304 (1995). [CrossRef] [PubMed]

43. M. Essenpreis, C. E. Elwell, M. Cope, P. van der Zee, S. R. Arridge, and D. T. Delpy, “Spectral dependence of temporal point spread functions in human tissues,” Appl. Opt. 32(4), 418–425 (1993). [CrossRef] [PubMed]

44. M. Ferrari, Q. Wei, R. A. De Blasi, V. Quaresima, and G. Zaccanti, “Variability of human brain and muscle optical pathlength in different experimental conditions,” Proc. SPIE 1888, 466–472 (1993). [CrossRef]

45. A. Kienle and M. S. Patterson, “Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,” J. Opt. Soc. Am. A 14(1), 246–254 (1997). [CrossRef] [PubMed]

46. 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] [PubMed]

47. G. H. Golub and C. F. Van Loan, Matrix Computations, 3rd ed, Johns Hopkins Studies in the Mathematical Sciences (Johns Hopkins University Press, 1996).

48. A. Jelzow, “In vivo quantification of absorption changes in the human brain by time-domain diffuse near-infrared spectroscopy,” PhD Thesis, TU Berlin (2013).

49. M. Cope, “The application of near infrared spectroscopy to non invasive monitoring of cerebral oxygenation in the newborn infant,” PhD Thesis, University College London (1991).

50. S. J. Matcher, M. Cope, and D. T. Delpy, “Use of the water absorption spectrum to quantify tissue chromophore concentration changes in near-infrared spectroscopy,” Phys. Med. Biol. 39(1), 177–196 (1994). [CrossRef] [PubMed]

λ / nm	m₁ / ps		V / (10⁻³ ns²)	m₃_,C /	µ_a / cm⁻¹	µ_s´ / cm⁻¹
	this study	Lit. [Ref.]		(10⁻³ ns³)
689	858 (5.8%)	911 (17%) [42]	126 (20%)	54 (32%)	0.14 (16%)	9.5 (9.7%)
797	794 (6.1%)	885 (7.2%) [43]	121 (19%)	58 (29%)	0.13 (16%)	8.0 (9.0%)
		910 (7.7%) [44]
828	759 (6.4%)	820 (17%) [42]	113 (20%))	57 (30%)	0.14 (16%)	7.5 (9.8%)

Separation of superficial and cerebral hemodynamics using a single distance time-domain NIRS measurement

Abstract

1. Introduction

2. Experimental

3. Methods

Moments of DTOFs

Reconstruction of homogenous absorption changes (RHAC)

Reconstruction of absorption changes based on a layered geometry (RLAC)

Calculation of hemoglobin concentration changes

4. Validation of methods on a two-layer phantom experiment

5. Application to in vivo data

Conclusion

Acknowledgments

References and links

Cited By

Figures (6)

Tables (1)

Equations (15)

Biomedical Optics Express