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

: In functional near-infrared spectroscopy (fNIRS) superficial hemodynamics can mask optical signals related to brain activity. We present a method to separate superficial and cerebral absorption changes based on the analysis of changes in moments of time-of-flight distributions and a two-layered model. The related sensitivity factors were calculated from individual optical properties. The method was validated on a two-layer liquid phantom. Absorption changes in the lower layer were retrieved with an accuracy better than 20%. The method was successfully applied to in vivo data and compared to the reconstruction of homogeneous absorption changes.


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][6][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][20][21] or the concept of timedependent mean partial pathlenghts (TMPPs) [22][23][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].

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 −1 and µ a = 0.1 cm −1 . 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 −1 and high absorption. The absorption of the upper layer µ a,1 was increased from 0.1 cm −1 to about 0.2 cm −1 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 −1 . Then the absorption of the lower layer µ a,2 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,2 was changed from 0.1 cm −1 to about 0.2 cm −1 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 −1 and the remaining steps about 0.01 cm −1 .
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 6 s −1 per TCSPC module and 0.9⋅10 6 s −1 per wavelength. Anatomical T 1weighted 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].

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 is proportional to the detected light intensity. The n-th normalized moment m n of a DTOF N(t) is defined as: The first moment m 1 is the mean time of flight. The n-th centralized normalized moment m n,C is defined as: 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): 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 1 , the variance V and the third centralized moment m 3,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 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 1 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 1 , 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 1 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. 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 −1 in both layers, r sd = 3 cm −1 , n in = 1.33, µ a,1 = 0.1 cm −1 , µ a,2 = 0.2 cm −1 , 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 .
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 1 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: 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]: 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.  Table 1). To our best knowledge values for V and m 3,C are reported for the first time. Unlike data for m 1 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.  We further used the values of m 1 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 1 and V. Optical properties derived in this way are presented in the lower panels of Fig. 3 and used in the next section.  Table 1.
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 1 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: 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 : 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 −1 to 0.35 cm −1 in steps of 0.01 cm −1 and from 4 cm −1 to 19.5 cm −1 in steps of 0.5 cm −1 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 1 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 1 and V.
With increasing scattering the position of the maximum sensitivity (black lines on color maps) of m 1 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 1 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 1 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.
In an anatomical model of the head z 1 reflects the thickness of the scalp, z 2 the distance to the gray matter and z 3 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: , , 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]: ( ) Here K is the covariance matrix of ΔM and in our case determined by the photon noise as: 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 1 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 timeconsuming 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: 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.

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 1 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,1 and µ a,2 .
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 1 variance V as a function of µ a,1 or µ a,2 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 −1 (nominal absorption change) which corresponds to a 2.5% to 5.6% relative change in µ a depending on the absolute value of µ a,1 or µ a,2 In addition, the fit helps to compensate for fluctuations in the courses of m 1 and V for µ a,1 between 0.15 cm −1 and 0.19 cm −1 , 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 −1 .
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,1 is changed Δµ a,ΔA produces the closest result to the true Δµ a,1 ; if µ a,2 is changed Δµ a,Δm1 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,1 and Δµ a,2 vary around the true values in a range of about ± 0.001 cm −1 . 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 −1 ( ± 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,1 and Δµ a,2 are shifted with respect to the reconstructed courses in Fig. 5 in different ways. If the thickness is too low and µ a,1 is increasing, both Δµ a,1 and Δµ a,2 are overestimated while a too high thickness leads to an underestimation. This behavior is reversed if µ a,2 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,1 and Δµ a,2 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,2 caused by a change of µ a,1 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 −1 and 10 cm −1 , respectively. We found minor effects on the reconstructed values of Δµ a,1 and larger deviations for Δµ a,2 . In the worst case only 30% of Δµ a,nom were retrieved for Δµ a,2 at µ a,2 ≈0.2 cm −1 . This example demonstrates that by taking into account adapted homogeneous BOP the RLAC method achieves a more reliable quantification of deep absorption changes.

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 1 = 7.4 mm and z 2 = 12.8 mm. The depth z 1 was obtained as the averaged distance between the scalp surface and the skull at the position of fNIRS sources and detectors. Analogously, z 2 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 1 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.
Concentration changes of oxygenated hemoglobin obtained using the RHAC method from all three moments ΔA, Δm 1 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 1 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 depthdependent 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.