Method for the discrimination of superficial and deep absorption variations by time domain fNIRS

: A method for the discrimination of superficial and deep absorption variations by time domain functional near infrared spectroscopy is presented. The method exploits the estimate of the photon time-dependent pathlength in different domains of the sampled medium and makes use of an approach based on time-gating of the photon distribution of time-of-flights. Validation of the method is performed in the two-layer geometry to focus on muscle and head applications. Numerical simulations varied the thickness of the upper layer, the interfiber distance, the shape of the instrument response function and the photon counts. Preliminary results from in vivo data are also shown.


Introduction
Noninvasive monitoring and imaging of cortical hemodynamics in the human brain by optical methods can complement the current neuroimaging modalities such as electroencephalography (EEG) and functional magnetic resonance imaging (fMRI). Optical methods like near infrared spectroscopy (NIRS) [1] and diffuse correlation spectroscopy (DCS) [2] provide insight on the vascular activities occurring in the brain cortex in resting conditions or elicited by somatosensory, motor or cognitive stimuli. These optical methods have also the advantages of being applicable at the bedside, thanks to their relative compactness, and for long term or repeated measurements, thanks to their complete noninvasiveness. Such techniques typically employ continuous wave (CW) light sources (e.g. laser or LED) to inject near infrared (NIR) light in the head by means of a probe (i.e. an optical fiber) gently positioned on the scalp and firmly hold by a cap or a biocompatible adhesive tape. While the light signal propagates in the head, thanks to the relative transparency of the tissue in the NIR spectral range, it is attenuated by the combined effect of scattering and absorption from tissue structures and chromophores. Oxygenated and deoxygenated hemoglobin (O 2 Hb and HHb, respectively), lipid and water are the main chromophores interacting with the NIR light. The light signal can be detected by a proper sensor typically consisting of a large core optical fiber (to optimize the collection of photons) and a photodiode or a photomultiplier coupled to suitable electronics. The sensor is placed at a distance of a couple of centimeters from the injection to maximize the sensitivity to photons that have travelled deep to the brain cortex (crossing the scalp, the skull and the cerebrospinal fluid) before being eventually detected. NIRS exploits the changes in the detected light signal occurring in the brain cortex as a consequence of the so called neuro-vascular coupling effect: a hemodynamic response originates locally in the cortex to supply energy to the working neurons [3]. Changes in O 2 Hb and HHb concentration in the cortex can be estimated by means of simple physical models like the Lambert-Beer law [4]. When studies are conducted to explore the overall brain functioning, in the literature the term functional NIRS (fNIRS) is often adopted.
A major concern with fNIRS is the possibility that systemic physiological changes (e.g. blood pressure or heart rate changes) might contaminate the detected signal [5,6]. These systemic phenomena typically occur in the superficial tissues (mainly the scalp) and can confound the readings of the instruments, especially when task related physiological responses happen [6,7]. A multi-channel configuration with sensors placed at different distances from the injection point is the solution that is currently being adopted in CW fNIRS to reduce the crosstalk between the superficial systemic signal and the deep cortical signal [8][9][10]. Conversely, in the time domain (TD) approach to fNIRS, picosecond pulsed lasers and fast photodetectors are used to acquire the photon distribution of time-of-flights in the head. A direct and intuitive link between the photons time-of-flight and the pathlength traveled within the probed tissue does exist: the longer the time-of-flight the larger the penetration depth [11,12]. Thus, by exploiting the information encoded in the photon time-of-flight, signals coming from different depths of the tissue can be discriminated. This holds for each acquisition channel (i.e. for each source-detector pair), without the need of implementing a multi-channel approach like in CW fNIRS.
However, the accurate quantification of the hemodynamic changes in the intra-and extracerebral compartment is still an open issue in TD fNIRS. Several approaches have been reported in the literature to tackle this problem like the use of moments of the photon distribution of time-of-flights [13,14] and the use of Mellin-Laplace transform [15,16]. Since these methods operate on the whole photon distribution of time-of-flight we have the feeling that the original peculiar characteristics of the TD approach of relating photons to depth is somehow unexploited. We have always been in favor of a simpler heuristic approach based on the selection of early and late time-gates to more directly address superficial and deep structure in the head [17]. A similar approach has been followed by other researchers [18]. While this approach proved effective to detect cortical changes in adult, it however lacks accuracy, mainly because it is not so trivial to precisely compute the time-dependent photon pathlength that is the time spent by photons in different structures of the head.
In this work we present an improved method for the discrimination between superficial and deep variations of the absorption properties of a multilayered tissue. Key point is a refined computation of the pathlength traveled by photons within each layer the tissue is composed of, by taking into account the non-idealities of the system set-up used to perform the measurements (as described by the instrument response function, IRF) and the heterogeneous structure of the human head.
In the following sections we will present the theoretical basis of the method, the validation by means of simulations in a two-layer model of the head, a preliminary application to in vivo data and a discussion of the results.

Theory
To a good approximation the head can be modeled as a heterogeneous structure, composed by different tissues, such as scalp, skull, cerebrospinal fluid, gray and white matter. In order to describe small and independent absorption changes in each domain, a different formulation of the modified Lambert-Beer's law can be exploited; for small absorption changes, the timeresolved reflectance (TRR) curve can be casted as [11,13]: where L j (t) (also known as time-dependent mean partial pathlength, TMPP) is the pathlength traveled in the j th domain by photons detected at time t, while Δµ a,j represents the absorption variation happening in the same domain, with respect to the unperturbed TRR curve R 0 (t). We recall that the TMPPs L j (t) can be derived as [12]: In this way L j (t) depends only on the characteristics of the unperturbed TRR curve R 0 (t) .
In an ideal experiment R(t) and R 0 (t) are measured and Δµ a,j are derived from Eq. (1) once the estimates of L j (t) are available. In a realistic time-resolved set-up, however, we hardly have access to the TRR curves R(t) and R 0 (t), that are relative to a δ-like time response of the instrumentation. More easily, we have to face an experimental set-up with an instrument response function s(t) that has stated temporal characteristics. As a matter of fact, the measured TRR curves ( ) R t  and 0 ( ) R t  can be related to the ideal TRR curves R(t) and R 0 (t) by a convolution with the IRF s(t): In these conditions, we still assume to be valid a relationship between perturbed and unperturbed TRR curves with the same form as Eq. (1): where ( ) R t  and 0 ( ) R t  are TRR curves measured by the experimental set-up and ( ) temporal functions that can be calculated following the same procedure as in Eq. (2): We will now demonstrate that ( ) j L t  represent the TMPPs in this realistic case. The TMPPs in the realistic case can be calculated as a weighted mean of the TMPPs in the ideal case; if we observe that 0 ( ') ( ') ' dN s t R t t dt = − represents the number of photons detected at time t which entered the diffusive medium at a time between t' and t' + dt', the TMPP j spent in the j th domain by photons detected at time t can be calculated as: where we exploited Eq. (2) for expressing the TMPP in the ideal case and swapped the derivative and integral operations. A comparison between Eq. (6) and Eq. (5) concludes our demonstration. Finally, in order to make the implementation of the described method for calculating the absorption changes Δµ a,j more reliable, the measured TRR curves can be sliced into a few time-windows, into which the photons are summed up (see Fig. 1): as a matter of fact, relaying on the raw time resolution, that can be also of a few picoseconds, is somehow useless and also worsening due to the bettering of signal-to-noise ratio that can be gained with temporal binning. Accordingly, g R  , i.e. the photons detected in the time windows (t g , t g + 1 ), can be expressed as: where 0,g R  are the photons detected in the same time window in the unperturbed case, while , g j L  are the mean partial pathlengths (MPPs) relative to the g th gate and can be calculated as: Next, Eq. (7) can be casted into a system of linear equations: that can be inverted in order to estimate the absorption variations Δµ a,j , provided that the MPPs , g j L  , which depend only on properties of the unperturbed condition, can be calculated.

Simulations
In this section we describe the numerical simulations we performed to validate the proposed method and to test its accuracy in measuring absorption variations. In details, the influence of various parameters, related both to the geometry and optical properties of the sample under investigation, and to the experimental set-up configurations adopted for the measurements, was studied. In this way the range of applicability of the method has been assessed.
We have chosen the two-layer geometry for the validation of the proposed method (see Fig. 2(a)). To a first approximation in fact this can be used to mimic both NIRS measurements on the muscle, where a superficial layer (skin and a lipid tissue) cover the deeper muscle tissue, and fNIRS measurements on the human head, where an extra-cerebral layer (scalp, skull and cerebrospinal fluid) protects the intra-cerebral layer (gray and white matter). Reference optical properties of the medium were set at 0.1 cm −1 for the absorption coefficient µ a and at 10 cm −1 for the reduced scattering coefficient µ s '. µ a was then changed respectively in the upper or in the lower layer ranging from 0.01 to 0.2 cm −1 in steps of 0.01 cm −1 . The TRR curves were calculated by exploiting a suitable home-made software which supplies the solution of the Diffusion Equation for the adopted geometry [19]. Poisson noise was added to the simulated curves to mimic real measurements, while 100 repetitions for each experiment were considered to achieve a significant statistics. The TRR curves were than divided into 10 temporal gates with same width, using as the starting and ending points those with one thousandth of the peak counts on the leading and trailing edge of the TRR curve, respectively.
Next, by means of custom-made Matlab routines (MathWorks Inc., Natick, MA), the MPP relative to each layer was derived from Eq. (2) and Eq. (5) by calculating numerically the derivatives therein. An example of resulting TMPP, i.e. the mean path spent in each layer for every temporal gate, is plotted in Fig. 2(b) for a time-resolved set-up with an ideal time response (δ-like IRF) and for one with a real IRF (named s 2 , its properties will be described in paragraph 3.6). A typical gate width produced by this method is about 300-400 ps, depending on the optical properties of the medium, the interfiber distance and the IRF of the device. (b) Simulated mean pathlength followed by photons of each temporal gate within both layers, considering an ideal δ-like IRF (dashed lines) and a real IRF (solid lines). Interfiber distance is set at 3 cm, while the thickness of the upper layer is set at 1.5 cm.

Influence of the number of gates
Firstly the influence of the number of temporal gates was investigated (5, 10, 20, 40 gates). An ideal δ-like IRF was considered, upper layer thickness was set at 1.5 cm, interfiber distance at 3 cm and the curve counts at 10 8 ph. The outcomes demonstrate a negligible dependence of absorption estimates on this parameter (results not shown). The chosen value of 10 gates allows to maintain a reasonable computational time and a good signal-to-noise ratio in every single time window.

Influence of the thickness of the upper layer
In this section we have studied the effect of the thickness of the upper layer in the estimation of the absorption variations calculated by the proposed method.
The distance between the skin and the underlying brain cortex, considered as the thickness of the upper layer, presents high variability between subjects. Okamoto et al.
[20] estimate the depth from head surface to brain being 16 ± 2.3 mm in the regions close to C3 position of the 10/20 EEG International System; however Okada et al. [21] studied the effect of the superficial layer thickness on the NIRS sensitivity by considering the head as a multilayered medium, using data of previous papers [22][23][24], and reported a wide range of values (skincortex distance spanning from 7.25 to 20 mm). We have therefore considered the following values for the thickness of the upper layer: 0.5, 1, 1.5 and 2 cm.  Figure 3 shows the results obtained when an ideal δ-like IRF is considered. An elevate photon number was used in the simulations (10 8 ph), while the source-detector distance was set at 3 cm. The graphs are disposed as follows: the estimated values of the absorption coefficients are reported in the left column, calculated for both layers with the new proposed method; the corresponding relative errors between estimated values and nominal values are depicted in the right column. We note that, due to the ideal conditions we adopted for the IRF in this case, the error reported in the right column can be considered as the systematic error affecting the absorption retrieval caused by the proposed method. In the top row the graphs show the outcomes obtained by varying the absorption in the upper layer while keeping unchanged the absorption in the lower layer; in the bottom row the results are relative to a change of absorption in the lower layer, while maintaining constant the absorption of the upper one. The square-shape markers indicate µ a calculated for the upper layer, while the diamond-shape series represents µ a of the lower layer. In the graphs are also reported the error bars for each retrieved value of µ a , calculated over the 100 repetitions of the measurement. This kind of variability affecting µ a represents the stochastic error. Finally, the reference values are represented with a black dashed line for the superficial medium and with a black solid line for the deeper one. The color code is for identifying different thicknesses of the upper layer.
The absorption coefficients relative to the superficial layer are always well assessed, except a weak underestimation for thin thicknesses (0.5 cm) when the absorption increases in the underlying medium. Conversely, µ a of the deeper medium presents always a systematic underestimation error, both for positive or negative µ a changes happening in the upper or lower layer. As expected for the deep layer, the larger the thickness of the upper layer the higher the systematic error of retrieved µ a values. Also the stochastic errors affecting µ a values increase with the increasing of the thickness of the upper layer, but they are very small in any case, due to the high photon number considered.
Finally, from Fig. 3 it emerges that the retrieved absorption values are affected by some systematic errors also in a set-up with an ideal temporal response: as a matter of fact, they are due to the approximation assumed in Eq. (1) for the description of the time-resolved curves, where the pathlengths spent by photons in each domain of the sample are assumed small and independent from the absorption changes. Then, these errors are somehow intrinsic of the method presented in this paper, and Fig. 3 can be seen as representative the limit of validity of Eq. (1). The effect of varying the source-detector interfiber distance ρ was examined for an ideal δlike IRF in this section. We focused on interfiber distances of 1, 2 and 3 cm, since they are the usual values adopted in real measurements. Curve counts were set at 10 8 ph and the upper layer thickness at 1.5 cm. The obtained outcomes showed that the mean absorption assessment is basically insensitive to a change of ρ (see Fig. 4). However an increase of results variability in the lower layer can be noticed for the smallest interfiber distance, mainly due to the fact that this short separation is more sensitive to superficial variations than deeper ones.

Influence of the bulk scattering
The influence of modifications in the scattering properties of the whole medium was investigated in this section, by assuming µ s ' values equal to 7.5, 10 and 12.5 cm −1 . An ideal δlike IRF was considered, upper layer thickness was set at 1.5 cm, interfiber distance at 3 cm and the curve counts at 10 8 ph. All the simulations demonstrate a negligible dependence of absorption estimates on the scattering coefficient of the examined probe (results not shown). In this section numerical simulations were performed in order to determine the effect of the curve counts on the absorption retrieved values. In real applications the count rate is much less than 10 8 ph/s and consequently the photon number per curve much less than 10 8 assumed so far. Then, TRR curves were simulated by setting the photon number at 10 5 , 10 6 , 10 7 and 10 8 ph. Here the thickness of the superficial layer was set at 1.5 cm, while the interfiber distance at 3 cm.

Influence of the TRR curve counts
We note that the decrease of the photon number in the TRR curves does not affect appreciably the systematic error of the retrieved µ a values (data not shown). On the contrary, by an inspection of Fig. 5, we can see how a decreasing of the photon number worsens their coefficient of variation (CV), defined as the ratio between the standard deviation of the data distribution and its mean.

Influence of the IRF temporal shape
Till now we considered a time-resolved set-up with an ideal temporal response. In this section, instead, we studied the effect of a not-ideal temporal response of the instrumental setup on the absorption retrieved values. To this purpose, numerical simulations were performed by adopting two realistic IRFs, with different temporal shapes (see Fig. 6). The IRF features depend on the instrumental set-up: in particular, the delivered power of the pulsed lasers and the arrangement of the optical coupling between sample and detector are crucial to this extent. Here, two characteristics of the temporal shape of the IRFs were considered: the full width at half maximum (FWHM) and the time decaying constant τ, calculated by fitting the tail of the IRFs with an exponential function. The first realistic IRF s 1 (t) we consider is characteristic of an instrumental set-up employing a hybrid photomultiplier [25]. This IRF presents a FWHM of 260 ps and its decaying constant τ is about 130 ps. The second realistic IRF s 2 (t) is relative to an instrument having 4-channel photomultiplier tubes as detectors, with a delivered power of 1 mW [26]: in this case the FWHM is 510 ps and τ is estimated being about 350 ps. In the simulations the TRR curve counts were set at 10 8 ph, while the upper layer thickness at 1.5 cm and interfiber distance ρ at 3 cm. Finally, for the sake of comparison, also the ideal IRF was considered in the simulations.  From Fig. 7, one can notice that the IRF s 1 shows a behavior very similar to the ideal δlike IRF. On the contrary, the IRF s 2 presents a larger systematic error in the estimation of the absorption values. In order to understand which of the two parameters characterizing the IRF, that is FWHM and the decaying time constant τ, is responsible for the worsening observed for IRF s 2 , further simulations have been performed.
In order to probe the effect of FWHM, rectangular-shaped IRFs with different widths were considered. In particular, rectangular IRFs with width set at 250 and 500 ps were simulated and compared to data computed by convolving with a δ-like IRF (null width). The change in width of the rectangular IRFs does not affect at all the retrieve values of absorption (results are not shown).
Then, a decreasing exponential function with different time constant τ was used to probe the influence of this parameter. Unlike FWHM, the change of τ strongly affects the results (see Fig. 8), in a way very similar to that previously noticed for the realistic IRF s 2 . In Fig. 8 results are reported for the following values of τ: 100, 200 and 400 ps, together with the ideal IRF (null τ). One can notice that if τ values are less than about 200 ps, the time-resolved setup has almost an ideal behavior: as a matter of fact, the estimated absorption coefficients have the same systematic error observed for a δ-like IRF. This threshold is mainly dependent on the tail slope, and thus on the analyzed absorption properties. Finally, we can now state that the worsening in the systematic errors observed for IRF s 2 is due to its high value of the decay time constant.

Dependence on a priori information
In the characterizations previously presented we assumed that the estimation of the different studied parameters was not affected by error. However, in a real in vivo measurement, the data are usually deteriorated by inaccuracies due both to a measurement error and to a wrong assessment of parameters values.
Hence, further numerical simulations have been performed on a bilayered medium in order to address the issue of how the accuracy of the method could worsen with a wrong estimate of different parameters, such as the interfiber distance, the baseline optical properties and the geometry of the probed medium.
As before, reference optical properties of the medium were set at 0.1 cm −1 for µ a and at 10 cm −1 for µ s '. The thickness of the upper layer was set at 1.5 cm, while the interfiber distance at 3 cm. The MPP relative to each layer was then derived from Eq. (2) and Eq. (5) for this configuration.

Uncertainty in the thickness estimation
As discussed in paragraph 3.2, the correct estimation of the geometry of probed media during in vivo fNIRS experiments is still an open issue. In fact, without the anatomical knowledge of the probed biological tissues, given for example by MRI data, the thickness of the different layers a tissue can be composed of is usually estimated by literature values, often presenting a wide range of variability. An investigation about this aspect is therefore of interest in the interpretation of in vivo results.
An uncertainty of about 30% was considered with respect to the reference value (1.5 cm): the influence of thicknesses of the upper layer of 1 and 2 cm has been evaluated.
Results in Fig. 9 show that, as expected, the choice of the geometry severely influences the accuracy of the outcomes. From Fig. 9, we can assess the lack in accuracy due to an underestimation or an overestimation of the thickness value of the upper layer, when the simulated data thickness is set at 2 cm or 1 cm, respectively, while the MPPs were computed for the reference thickness set at 1.5 cm.

Uncertainty in the interfiber distance estimation
An uncertainty in the interfiber distance estimate was considered as possibly affecting a real experiment. However, being this parameter more easily verifiable than previous one, a smaller variability range was considered. Source-detector distances were thus simulated at 2.7 or 3.3 cm (over-or underestimation of 10% from the reference value). The results are not influenced at all by an uncertainty on this parameter (outcomes not shown).

Uncertainty in the baseline absorption properties estimation
In a real measurement, a wrong assessment of the baseline optical properties of a medium could be related to different factors: the non-idealities of the experimental setup, an inaccurate estimate of the geometry, the intrinsic error committed by the physical model and so on. A study of the performance of the model in relation to this aspect is thus of practical interest when considering an in vivo investigation.
We firstly tested the influence of an inaccuracy in the estimate of the absorption coefficient of the superficial layer. A set of synthetic data was simulated with the baseline µ a of the superficial layer modified of ± 25% from the reference value, while keeping the baseline µ a of the lower layer at the reference value (0.1 cm −1 ); the three series displayed in A corresponding data set was simulated by modifying the baseline µ a of the deep layer of ± 25% from the reference value (baseline µ a at 0.075 and 0.125 cm −1 ) while maintaining at the reference value the µ a of the upper layer. Results are depicted in Fig. 11.
From both experiments one can notice that the computation of the µ a variations in the superficial layer are not affected at all by an uncertainty in the baseline µ a values. On the other hand, this latter strongly influences the estimated µ a changes in the deeper layer. Fig. 10. Calculated µ a variations (left column) and the relative error with respect to nominal absolute value (right column) as a function of µ a variation happening in the superficial layer (top row) or in the lower layer (bottom row). Upper µ a is marked by squares, lower µ a by diamonds. Each data series refers to a different value of the baseline µ a of the upper layer of the simulated data. Fig. 11. Calculated µ a variations (left column) and the relative error with respect to the nominal absolute value (right column) as a function of µ a variation happening in the superficial layer (top row) or in the lower layer (bottom row). Upper layer absorption is marked by squares, lower layer absorption by diamonds. Each data series refers to a different value of the baseline µ a of the lower layer of the simulated data.

Uncertainty in the scattering properties estimation
Simulations were also performed to study the influence of a wrong estimate of the reduced scattering coefficient of the probed medium.
Firstly data series have been simulated by modifying µ s ' of the upper layer of ± 25% from the reference value (µ s ' set at 7.5 and 12.5 cm −1 ) while keeping unvaried µ s ' of the lower layer at 10 cm −1 . Results are depicted in Fig. 12.
Moreover, synthetic data have been conversely created by maintaining the reduced scattering coefficient of the upper layer at 10 cm −1 while varying µ s ' of the lower layer (7.5 and 12.5 cm −1 ). Results are shown in Fig. 13, showing that the inaccuracy in the estimate of µ s ' of the upper layer introduces an error in the results of both layers, while the uncertainty in the lower µ s ' seems to weakly affect them.

In vivo applications
The method described so far has been applied on in vivo data, acquired by a multi-channel dual-wavelength TD fNIRS instrument [26], to preliminary validate its performance on such a kind of data. One right handed subject underwent the experiment, consisting of 10 repetitions of 20 s of baseline, 20 s of finger tapping with the right hand and 40 s of recovery, for a total duration of the experiment of 800 s. A pair of source and detector fibers (2 cm interfiber distance) was placed on the contralateral motor area over the hot spot position previously identified by Transcranial Magnetic Stimulation. The photon count rate was set to about 4•10 5 ph/s per wavelength. In order to have a reliable signal-to-noise ratio the acquisition time was set to 1 s and a moving average of 5 s was applied: in this way we processed TRR curves with about 2•10 6 total counts.
Firstly, the unperturbed absorption at each wavelength was calculated during the baseline period of each task repetition by a least square fitting of the TRR curves with a proper theoretical model, assuming the head as a homogeneous semi-infinite medium. The absorption variations at each wavelength with respect to each corresponding baseline were, then, computed applying the new method. In this case the head was modeled as a two-layer semi-infinite medium with a thickness of the upper layer set at 1.5 cm, chosen by visual inspection of the anatomical MRI map of the subject.
The concentration of oxygenated (O 2 Hb) and deoxygenated hemoglobin (HHb) was calculated by the Lambert Beer's law assuming that these are the only chromophores contributing to the absorption. In Fig. 14

Discussion
We have described and validated a method for the discrimination of superficial and deep absorption variation by TD fNIRS. The method exploits the estimate of the TMPP in different domains of the sampled medium and makes use of an approach based on time-gating of the photon distribution of time-of-flights. The two-layer geometry was chosen for validation to focus on muscle and head applications. A quite wide range of variation for the absorption properties of the diffusive medium was used in the simulations, with consequently a quite large systematic error in the estimated parameters. In order to assess more precisely the range of applicability of the method, we observe that from previous in vivo fNIRS measurements the relative variation (contrast) of signal intensity between resting state and task period (subject performing a motor task, source and detector fibers centered on the contralateral motor region) is hardly larger than 10%. The simulations performed to validate the proposed method revealed that a contrast below 10% is reached in the range between 0.08 and 0.13 cm −1 (0.095 and 0.11 cm −1 ) for the lower (upper) layer absorption values, when the thickness of the upper layer is assumed to be 1.5 cm. Therefore, we will have in mind these absorption intervals, which can be considered of interest on the physiological viewpoint, when we will critically evaluate the outcomes in this study in next statements.
The analysis of the influence of different thicknesses of the upper layer revealed that the absorption properties of the superficial stratus are computed with a very weak error (less than 1%). Conversely those ones related to the deep layer are systematically underestimated and present, as expected, an error increasing with the thickness of the upper layer. However the error always remains below 5% in the physiological absorption interval previously defined. The interfiber distance and the bulk scattering seem not to influence the accuracy of the results as expected from the theory [12,27]. A larger interfiber distance (e.g. 3 cm) can be used to reduce the number of photons that travels in the superficial layers. This choice however determines a strong reduction in the SNR and in the contrast, therefore shorter interfiber distances (<1 cm) should be considered, provided a method to suppress early (superficial) photon is used [28]. Regarding the influence of the counts in the TRR curves, as expected, the lower the counts the higher is the stochastic error committed by the proposed method. However, the absorption assessment in the superficial layer seems not much influenced by the curve counts, with a CV always below 1%. Conversely, the drop in counts corresponds to a worsening in the precision assessment of µ a in the lower layer (especially for elevated thicknesses of the superficial layer). However if the thickness of the upper layer is not larger than 1.5 cm, in the physiological range discussed above the CV does not exceed the 6% for TRR curves with at least 10 5 counts. As for the influence of the IRF temporal shape, we can state that the crucial parameter is the decay time constant τ of the IRF, which has to be less than 200 ps to approach ideal conditions. On the contrary, the FWHM of the IRF seems not to have effect on the absorption assessment. These considerations can be useful when a TD instrumentation has to be set-up: a right tradeoff between detected signal and temporal broadening of the IRF, which both depend on the choice of the detectors and of the relay optical system, must be done by taking into account also this issue.
Further, we investigated the lack of accuracy of the method introduced in case of a wrong assessment of the geometrical and optical reference parameters used for the MPP calculation: this aspect belongs to the framework of the not yet resolved issue of the knowledge of priors in reconstruction problems. The influence of a wrong estimate of the upper layer thickness was firstly studied. Results show that a strong variation is noticed in the calculus of µ a changes for both layers. An underestimate of the upper thickness corresponds to the overestimate of upper layer µ a variations and to the underestimate of the deeper µ a variations (and vice versa). The knowledge of the anatomy of the subject's head, as well as the use of a more complex geometrical model (such as a three-layered model) will certainly improve the accuracy of the presented approach for in vivo applications. Another experiment was focused on the assessment of the effect of the uncertainty in the source-detector distance estimate. As expected, the penetration probability in the time resolved approach does not depend on the interfiber distance. Thus results are not affected by a wrong evaluation of this factor. As for the tissue optical properties, numerical simulations revealed a strong increase of the computation error of the µ a variations of the lower layer in presence of an uncertainty in the baseline absorption properties of both strati. Equivalent results are produced by the imbalance of absorption: an overestimate of the upper stratus brings to results comparable to those obtained by the underestimate of the lower one. The investigation about the effects of an uncertainty in the diffusion properties estimation revealed that an error committed in the evaluation of the upper µ s ' modifies the results in both layers, while the same error committed for the lower µ s ' faintly affects them. These outcomes suggest that some efforts must be spent during an in vivo experiment in order to minimize the error on these parameters. In a real experiment, an accurate estimate of the optical properties could be achieved for example by fitting the TRR curves acquired during a rest period by a theoretical model for bilayered geometry. The robustness of the estimate could be improved also by a preliminary measure of the tissue of interest by exploiting a multidistance setup.
Finally, the preliminary results obtained for in vivo data confirm the performance of the method in recording a functional activation. We note that the estimated hemoglobin concentration changes occurring during the motor task are about of 10 μM, in line with some results in the literature [29]. We notice that this value is almost 10 times larger than the results obtained by the Authors with the previous method [30, 31] and other results reported in the literature [for a review see 32]. The authors think that hemoglobin variations of 10 μM rather than 1μM are closer to what actually happens in vivo because there are evidences from phantom measurements that the values for hemoglobin changes during a motor task reported in literature were somehow underestimated [33]. However, we have to consider that the accuracy of the proposed method is affected by the choice of the geometrical model, as revealed by simulations discussed in paragraph 3.7.1. In fact, as one can expect, the same variation in the optical intensity registered on the head surface can be produced by a certain hemoglobin concentration change in a superficial layer or by a larger hemodynamic effect in a deeper layer. To this extent, time resolution of the measurement can help to distinguish with respect to a CW approach. Furthermore, the assumption that the superficial variations (actually happening on the scalp) come from a volume with a thickness equal to the skincortex distance, leads to an unavoidable error in modeling the structure of the probed tissues. Finally we observe that an advantage of the proposed method with respect to the previous one [17], beside the possibility to include the effect of the IRF and an overall better accuracy, is that there is no more the need to manually identify an early time-gate and a late time-gate in the TRR curve. This operation was in fact somehow dependent on the expertise of the operator. The whole TRR is now included into the analysis allowing for a more robust analysis.

Conclusions
A method for the computation of absorption variations in a heterogeneous medium, that can be interesting in TD fNIRS, is proposed. It is based on a precise computation of the mean photon pathlengths within the different domains of the medium, by taking into account the non-ideal temporal resolution of the instrumentation used for the measurements (that is its IRF). The performance of the new method has been validated by numerical simulations considering a bilayer semi-infinite medium. Both the influence of geometrical and optical parameters of the medium, and the temporal features of the instrumental set-up have been tested, in order to assess the range of applicability of the proposed method. Also a preliminary application of the method to in vivo TD fNIRS measurements has been successfully performed. The results obtained show a potential accuracy of about 5% in the retrieval of the absorption properties of the heterogeneous medium in the geometry considered. The outcomes of this study can be exploited as guidelines for the development of new devices intended for TD fNIRS. Future perspectives of this work are the use of more complex geometrical models, e.g. a three-layer model, to better mimic the actual head's structure and improving in this way the accuracy of the method when applied on in vivo data, and a comparison with the other methods used for discrimination of deep absorption changes.