Depth-selective data analysis for time-domain fNIRS: moments vs. time windows

: Time-domain measurements facilitate the elimination of the inﬂuence of extracerebral, systemic eﬀects, a key problem in functional near-infrared spectroscopy (fNIRS) of the adult human brain. The analysis of measured time-of-ﬂight distributions of photons often relies on moments or time windows. However, a systematic and quantitative characterization of the performance of these measurands is still lacking. Based on perturbation simulations for small localized absorption changes, we compared spatial sensitivity proﬁles and depth selectivity for moments (integral, mean time of ﬂight and variance), photon counts in time windows and their ratios for diﬀerent time windows. The inﬂuence of the instrument response function (IRF) was investigated for all measurands and for various source-detector separations. Variance exhibits the highest depth selectivity among the moments. Ratios of photon counts in diﬀerent late time windows can achieve even higher selectivity. An advantage of moments is their robustness against the shape of the IRF and instrumental drifts.


Introduction
Functional near-infrared spectroscopy (fNIRS) of the brain aims at investigating the hemodynamic response of the cerebral cortex to functional activation [1]. This technique realized by numerous instrumental developments including an increasing number of commercial devices is finding more and more widespread applications in research and clinics [2][3][4]. With fNIRS, changes in oxyand deoxyhemoglobin concentrations are measured that are derived from wavelength-dependent changes of the absorption coefficient of the tissue, while the scattering properties are assumed to be unchanged. The fNIRS technology relies on the recording of attenuation changes of light propagating through the tissue. The source and detector optodes are attached to the scalp a few centimeters apart from each other in reflection geometry. By using multiple sources and detectors, lateral mapping, optical topography or tomography of brain activation in a larger area, even the whole head, are feasible [5][6][7].
Due to the comparably low absorption in the near-infrared wavelength range, a certain fraction of photons propagating from source to detector is able to reach the cerebral cortex that is typically more than 10 mm beneath the skin surface in adult humans [8,9]. However, fNIRS always faces the problem that the photons must travel through the extracerebral tissue compartments. Signals measured in reflectance geometry contain components related to extracerebral hemodynamic changes that occur mainly in the scalp. Thus, in many of the fNIRS applications, especially in adult humans, it is essential to eliminate this superficial contamination which is typically related to changes in systemic physiological parameters [10][11][12][13][14][15][16].
Time-domain fNIRS is an emerging technology that offers a favorite solution to this issue [17]. By acquiring the distribution of times of flight (DTOF) of photons propagating from source to detector, the information content of time-domain measurements is generally superior to continuous wave (CW) and frequency-domain measurements and allows one to achieve depth discrimination of absorption changes. To obtain DTOFs, time-resolved diffuse reflectance is recorded with short (picosecond) light pulses and time-resolved detection systems. Several research groups constructed prototypes of time-domain optical brain imagers based on time-correlated single photon counting (TCSPC) [18][19][20][21][22][23] or time-gated CCD cameras [24,25]. They were applied in human in-vivo studies, aiming at better identification of cerebral signals and elimination of the extracerebral contamination [26][27][28][29][30][31][32]. By means of fast-gated detectors, non-contact scanning time-domain functional imaging of the brain became feasible [33]. The progress in technological development has led to more compact and less expensive photonic components that will foster a more widespread use of time-domain diffuse optical methods including time-domain fNIRS [34]. The performance of a multitude of time-domain fNIRS instruments has been assessed and compared on the basis of two standardized protocols, "Basic instrumental performance protocol" [35] and "nEUROPt protocol" [36].
Several approaches were investigated to analyze measured DTOFs in order to obtain depthresolved changes in absorption of the tissue. Changes in the full temporal profile of the DTOF were analyzed using the concept of time-dependent mean partial pathlengths [37][38][39] or applying the microscopic Beer-Lambert law to the tail of the profile [40]. In most cases, however, the analysis starts from a reduced number of quantities characterizing the DTOF profile, following the pioneering studies of Arridge and Schweiger on various data types in time-domain optical tomography [41,42]. The motivation for the applications of such temporal filtering is to reduce the redundancy in the data since the temporal profiles exhibit no high-frequency features, and to substantially improve signal-to-noise ratio. Various research groups employed, in particular, integrated photon counts in various time windows [24,[43][44][45], moments [9,46,47], Mellin-Laplace moments [48][49][50] or Fourier components of the DTOF [51].
The measurands in the focus of this work have been successfully employed in laboratory and in-vivo experiments. Variance has been demonstrated to efficiently suppress the influence of superficial absorption changes in functional stimulation experiments [52] as well as during indocyanine green bolus tracking for perfusion assessment in animals [53] and in stroke patients [27,28]. The use of moments enabled the robust reconstruction of absorption changes in two-layer systems [9,20,46,47]. Analysis of early and late time windows has played an important role in time-domain optical mammography to detect lesions that differ from the surrounding tissue in scattering or absorption coefficients [54,55]. A combination of photon counts in late and early time windows enabled depth discrimination in time-domain fNIRS of the brain [24,44,56]. In non-contact optical brain imaging, deep and superficial hemodynamic changes were analyzed based on time windows of DTOFs [25,33]. Moments as well as photon counts in time windows served as examples of measurands to which the "nEUROPt protocol" was applied [36].
It turned out that some of the moments and time-window measurands provide fNIRS signals that are more brain specific, i.e. less prone to contamination by superficial absorption changes, than others. To assess this property of various measurands we introduced depth selectivity which we defined as the ratio of sensitivities of deep vs. superficial absorption changes [36]. A measurand that exhibits high depth selectivity can be suitable for empirically separating absorption changes in the brain from superficial changes, even without the application of advanced reconstruction algorithms that solve the inverse problem.
The assessment of depth selectivity as well as robustness against experimental artefacts is crucial for the optimization and comparability of the various analysis procedures in time-domain fNIRS. So far, the different measurands, either moments or time windows, have been used and investigated separately in numerous studies. However, a comparison of their features on a common basis, in a systematic and quantitative manner, is still lacking, even though some aspects of comparison were touched in our previous study [36]. Specifically, this previous study did not deal with spatial sensitivity profiles, nor did it consider combinations of photon counts in different time windows.
In the present work we characterize and compare the various approaches to time-domain fNIRS data analysis that are based on (1) moments (integral, mean time of flight and variance), (2) photon counts in individual time windows and (3) photon counts in different time windows. The major goal of our study is to identify those quantities which exhibit the highest depth selectivity. We employ diffusion-based perturbation simulations to calculate spatial distributions of the sensitivities of these measurands to small localized absorption changes. These sensitivity profiles are first compared qualitatively. The determination of depth selectivity enables a quantitative comparison of the different measurands. Furthermore, we study the influence of the instrument response function (IRF) on depth-selective measurements. We conclude with an overview of the relevant features of moments and time-window measurands.

Definition of moments and time windows of time-of-flight distributions
When time-resolved diffuse reflectance is measured on a turbid medium such as the human head, a DTOF N(t), i.e. photon counts vs. photon arrival time at the detector, is obtained. Data reduction is aimed at deriving a few characteristic quantities (here also referred to as measurands) M that still contain the major information but have a considerably improved signal-to-noise ratio compared to the photon counts in the individual time channels of the original DTOF. It is advantageously performed by calculating moments [41,42] or summing up photon counts within time windows [54,57,58]. These operations can be represented by where the weight function G(t) has the following forms: • G ≡ 1 for the total photon count N tot (zeroth moment) • G ≡ 1 for t ak ≤ t ≤ t bk and 0 elsewhere for the photon count N k summed up within a time window between limits t ak and t bk • G = t p /N tot for m p , the normalized moments of order p • G = (t − t ) p /N tot for the normalized central moments of order p In particular, the first moment m 1 (mean time of flight, unit: s) and the variance V (second central moment, unit: s 2 ) of a DTOF are and respectively. The photon count N k summed up within a time window between limits t ak and t bk is In addition to individual time windows we consider ratios of photon counts in different time windows as measurands. Figure 1(a,b) illustrates a DTOF N(t), the weight functions G(t) and the related integrands G(t)N(t) for the moments m 1 and V. A DTOF sectioned into time windows of equal width of 500 ps is shown in Fig. 1(c). This width is applied throughout the present work. This seemingly somewhat arbitrary choice is motivated by the compromise between catching the essential temporal features of the DTOF, achieving a good signal-to-noise ratio, specifically in the range of late photons, and avoiding excessive redundancy.  [59]. N(t) and G(t)N(t) in (a,b) were normalized to their respective maxima. Poisson noise shown in (a,b) was simulated for N tot = 10 6 and 5 ps width of the histogram time bins.
The complete DTOF enters the calculation of moments in any case, although with a timedependent weight. For m 1 , the weight increases linearly with time. For V, the weight grows heavily toward late photons according to a parabolic dependence while the influence of the central part of the DTOF is suppressed. There is a contribution of early photons as well, but toward earlier times, the DTOF reaches zero before the weight function rises dramatically. For late photons, however, the DTOF changes more slowly with time as compared to the leading edge. The weight function effectively uplifts the contribution to the integral at low signal level (photon counts) more and more with increasing time. Eventually the decay of the DTOF dominates and the integrand approaches zero as long as there is no offset (e.g. residual background level) or background-related noise. The third central moment and even higher moments are increasingly prone to these artefacts as well as to statistical noise. Therefore, they do not play a significant role in practical applications.
Since noise considerations are crucial for fNIRS applications we included in Fig. 1(a,b) a representation of counting noise, the dominating source of noise in TCSPC measurements, employing Poisson random numbers. For all measurands considered here, the statistical uncertainty (standard deviation) σ can be easily estimated for counting noise that follows a Poisson distribution with σ N = √ N. This relationship applies for the individual time bins of N(t) as well as for N k and N tot . Photon counts in different individual time bins are statistically uncorrelated, and the same is true for N k in different time windows k as long as they do not overlap. For the moments m 1 and V, simple analytical expressions result [47,60]: σ m1 = V/N tot The time-dependent weight function for moments also implies an increasing contribution of noise of the decreasing amount of photons in the tail of the distribution. This fact should be taken into account when calculating these moments for measured DTOFs.
From the definitions provided before, several relevant features of moments and photon counts in time windows of measured DTOFs can be derived immediately: (1) Due to normalization to N tot , the moments m 1 and V are independent of intensity and thus independent of drift and fluctuations in the laser power as well as in the coupling of the optodes to the skin. Moreover, variance is independent of the origin of the time axis. Thus, V is not influenced by any timing drift of the laser pulses as long as such drift is slow compared to the collection time. Consequently, variance is a rather stable signal, with its uncertainty being mainly determined by statistical (photon) noise. This feature turned out to be advantageous in a study on stroke patients [28].
(2) For both moments m 1 and V, deconvolution of the IRF is simply achieved by subtraction. Consequently, changes in the moments (∆m 1 , ∆V) do no longer depend on the IRF, in particular on laser pulse width and dispersion in fiber optics which may differ from experiment to experiment. We note in passing that deconvolution by simple subtraction is still valid for the third central moment but no longer for higher central moments [61].
(3) Photon counts in time windows are sensitive to power instabilities, to timing drifts of the laser pulses as well as to the shape of the IRF. However, ratios of photon counts in different time windows are independent of intensity fluctuations.

Perturbation simulations
In fNIRS, the sensitivity of a measurand to small changes in the absorption coefficients in different compartments of the head is of interest. To quantify and compare the sensitivities for the various measurands we employed a simplified model, a perturbation approach based on the diffusion theory for photon migration and Born approximation for which analytical solutions are available for certain geometries. The general formalism and the concept of "photon measurement density functions" were described by Arridge [62]. A small localized absorption change ∆µ a in a small volume dv at a position (x, y, z) in an optically turbid, otherwise homogeneous medium leads to a small change in the DTOF measured at the boundary of the medium. The solution for an absorption perturbation is constructed as follows (see, e.g. [63]). It consists of the temporal convolution of two Green's functions, (1) the time-dependent fluence rate at (x, y, z) inside the medium due to a Dirac delta light pulse injected at the source position, (2) time-dependent fluence rate at the detector position on the surface due to a source (temporal delta function) at (x, y, z), followed by the integration over the volume of the absorbing inclusion and the calculation of the perturbed diffuse transmittance or reflectance [64,65]. The spatial integral simplifies to a factor ∆µ a dv for a small inclusion with homogeneous ∆µ a . Analytical expressions for the perturbation of time-resolved reflectance for a slab geometry and extrapolated boundary conditions were provided in Appendix B of [65] as well as in [63]. This solution was adopted in the present study, in the limit of large slab thickness corresponding to a semi-infinite medium.
The simulations were performed for the following geometry: The point-like (isotropic) source and the point-like detector are located at (0, 0, z 0 ) and at (ρ, 0, 0), respectively, with z 0 = 1/µ s and ρ being the (fixed) source detector separation. The z axis is perpendicular to the surface and pointing into the medium. The position of the point-like absorber (x, y, z) is varied inside the medium in all three dimensions.
The simulations were carried out for two values of the source-detector separation ρ, i.e. 3 cm and 2 cm. The homogeneous background optical properties chosen were the typical set frequently adopted for the human head, µ s = 10 cm −1 , µ a = 0.1 cm −1 , refractive index of the medium n = 1.4. A small absorption perturbation with ∆µ a = 0.001 cm −1 (relative change of 1%) extended over a volume dv = δ 3 was assumed. The strength of the related absorption perturbation for a single grid point is ∆µ a δ 3 .
The perturbations ∆R(ρ, t; x, y, z) of time-resolved reflectance were simulated in M for a grid of points (x, y, z) with a spacing of δ = 0.5 mm in all directions and a temporal step width of 5 ps. The time range (5 ns) as well as the ranges of the spatial coordinates were chosen large enough to avoid any influence of clipping. The temporal profiles ∆R(ρ, t; x, y, z) were calculated, e.g., for -2 cm ≤ x ≤ 5 cm, 0 ≤ y ≤ 3 cm (data for y < 0 were derived based on symmetry) and 0 < z ≤ 3 cm, i.e. for nearly 590 000 grid points.
The unperturbed (baseline) time-resolved reflectance R 0 (ρ, t) was obtained from the solution of the diffusion equation for the homogeneous semi-infinite medium with extrapolated boundary conditions [59].
For each absorber location, the perturbed time-resolved reflectance was calculated as Finally, all measurands M were derived for the perturbed and unperturbed cases, M(x, y, z) and M 0 (for the sake of simplicity, we omitted the dependence on ρ here and in the following) from the reflectance profiles R(ρ, t; x, y, z) and R 0 (ρ, t), respectively. The related changes were calculated for each grid point according to

Sensitivities of the various measurands
For the moments, Eq.
The coordinates (x, y, z) were omitted, the subscript 0 denotes the unperturbed case. Accordingly, attenuation changes in time windows (indexed by k) were derived, For small changes, the relationship between attenuation and photon counts simplifies according to ∆A ≈ −∆N/N, for N tot as well as N k . The use of attenuation is preferable over photon counts since it allows one to express the contrast in all moments consistently as differences. Moreover, we investigated ratios of photon counts within a late (index l) and an early (index e) time window. As long as the changes are small, relative changes in this ratio can be expressed as changes in the difference of the related attenuations, The sensitivity factors (S M ) for all these quantities M were obtained as The symbol sign M stands for the sign that was chosen such that for simplicity of the following plots, sensitivities would be positive for a homogenenous increase in absorption. Such absorption change leads to positive ∆A and ∆A k (increase in attenuation due to decrease in N tot and N k ), but negative ∆m 1 and ∆V (decrease in mean time of flight and variance), therefore we set sign A = sign Ak = 1, sign m1 = sign V = −1. Regarding the correspondence with previously introduced terminology [46], the sensitivities S A , S m1 and S V , after integration over the volume considered (voxel or layer) are equal to MPP (mean partial pathlength), −MTSF (mean-time-of-flight sensitivity factor) and −VSF (variance sensitivity factor), and S Ak is closely related to the time-dependent mean partial pathlength (TMPP) introduced in [37]. The only difference between S Ak and TMPP, beside the presence of dv in the denominator of Eq. (9), is that TMPP is obtained from S Ak in the limit of infinitesimal width of time windows (corresponding to the width of individual time bins in TCSPC measurements, typically on the order of 10 ps). For a layered absorption change, i.e. when ∆µ a depends only on depth z, the sensitivities (superscript L) are obtained by integrating over x and y across the whole sensitivity profile,

Depth selectivity
Rather than the absolute sensitivity, it is often more relevant to know how selective a measurement is to an absorption change in a certain tissue compartment. Specifically, in fNIRS of the brain, the time course of the signal change is of primary interest while its absolute magnitude is less relevant and depends on many factors (e.g. overlap of the sampled volume with the activated region in lateral and depth directions, coupling between optodes and skin). If a measurement is sensitive to changes in a deep compartment only, the application of a complex reconstruction procedure is not necessary.
To characterize and compare different measurands with respect to the specificity of detecting absorption changes in the brain, we previously introduced the characteristic "depth selectivity" [36] and defined it as the ratio of sensitivities of the measurand M to absorption changes in the lower (indexed by 2) and the upper (indexed by 1) layers. To simulate depth selectivity based on the 3D sensitivity profiles S M (x, y, z), the sensitivities are integrated over the volume of the layers considered, before taking the ratio, The scattering properties are assumed to be the same in both layers. High (low) depth selectivity means dominating sensitivity to absorption changes in the lower (upper) layer. The major advantage of such a relative (dimensionless) characteristic is its suitability for the comparison of measurands with different dimensions (units) as in the case of moments, specifically m 1 and V, vs. photon counts in time windows.
In the context of performance assessment of instrumentation [36], the related test was performed on a specified two-layer phantom with fixed thickness of the upper layer. In the present work we apply "depth selectivity" in a more generalized manner.

Spatial sensitivity profiles
The capability of the various measurands M to achieve depth discrimination of absorption changes can be compared by plotting their 3D sensitivity profiles S M (x, y, z). Such spatial distributions are represented in Fig. 2 for attenuation, mean time of flight and variance. The color scales were adjusted to display all figures for each measurand with a reasonable dynamic range and to enable qualitative comparison for the different measurands. For the y-z cuts (second row), the maximum sensitivity was assigned the same (red) color for all columns (measurands M) while the maximum sensitivity is saturated (black) in some of the plots. S M = 0 is represented by the same light blue color in all plots. In each column, points (x, y, z) that are common to two or more plots have the same color. The most relevant observation is the well-known increase of the effective probing depth with increasing order of the moment when comparing the plots of the various measurands (columns) in Fig. 2. In particular, the depth of the maximum sensitivity at x = ρ/2 increases from 0.5 cm for A (left column) to 1.1 cm for m 1 (middle column) to 1.4 cm for V (right column). For A, the sensitivity is highest in regions in the vicinity of the source and detector positions (left column, first row). This is no longer the case for V (right column) for which the values are nearly equal within the banana-shaped region of high sensitivity (red) in the x-z cut. The behavior of m 1 (middle column) is in between that of A and V.
For m 1 as well as for V there is a region of negative sensitivity between source and detector near the surface of the medium. This means that the presence of absorbers in these regions increases m 1 and V while a homogeneous increase in absorption throughout the medium always reduces m 1 and V. The negative regions are surrounded by positive sensitivities, as can be seen, e.g. in the section at z 1 . At the transition between positive and negative regions, a point-like absorption perturbation is virtually "invisible" to that measurand. Such effect was first discussed for frequency-domain measurements [66]. An absorber near the surface preferentially eliminates photons with near-surface trajectories that have a shorter time of flight than the average, thus increasing m 1 and V. The appearance of a negative region can also be understood to be a consequence of normalization to N tot in the definition of m 1 and V: nominator and denominator exhibit different sensitivity profiles, with the denominator being more sensitive near the surface than the nominator.
Focusing on the aspect of lateral spatial resolution, the sensitivity maps plotted in Fig. 2 also exhibit different lateral extensions, as can be derived from the x-z and y-z sections. Regions with positive sensitivity clearly extend beyond the source and detector positions in x direction (i.e. for x < 0 and x > ρ) for m 1 and V, and the extension in y direction increases from A to m 1 to V. This behavior can be explained by the increasing contribution from photons with long times of flight and trajectory lengths.
By plotting a measure of depth selectivity, the sensitivities for the three moments can be related to each other quantitatively (see Fig. 3). The comparison of the magnitude of this spatially resolved selectivity is instructive. Specifically, from the y-z cuts it is obvious that, in relation to a homogeneous superficial change, m 1 and V have considerably higher maximum sensitivity values than A. For more details on depth selectivity refer to Section 3.2. Fig. 3. Spatial distributions of depth selectivity for moments A, m 1 and V, for depths z > 0.5 cm, obtained by normalizing the 3D sensitivities plotted in Fig. 2 to the total sensitivity within the upper layer (z ≤ 0.5 cm). The third row refers to the cut at z = z 2 = 1 cm. The color scales are the same in all panels, ranging from -0.48 to 1.2 times the maximum for S V in the y-z cut. Figure 4 illustrates sensitivity maps for time windows, in comparison with maps for moments. Other than the 3D representations in Fig. 2, the sensitivity maps in Fig. 4 are 2D maps related to the sum of sensitivities along the y direction in each x-z position, corresponding to the effect of a line absorber extended in y direction on the measured signals. This projection effectively boosts the visibility of the central regions of the sensitivity profiles relative to the narrow regions in the vicinity of source and detector. In the top row, the visibility was further enhanced by identifying the color bar maximum with 50% of maximum sensitivity in each plot.
The sensitivity profiles for time windows (Fig. 4, top row, panels 1-6) exhibit the well-known increases of penetration depth with time. For short times of flight, the sensitivity remains confined to shallow regions. With increasing time of flight, the probability of photons visiting deeper region increases successively. The sensitivity profile for attenuation A related to the total photon count in the DTOF (top row, panel 7) is equivalent to the sum of the (non-normalized) sensitivities for all consecutive time windows covering the complete DTOF. It is similar to the profiles for time windows in the central part of the DTOF, cf. A 2 and A 3 .  Fig. 1(c). Bottom row, panels 1-5: Attenuation difference A l -A e between late (l) and early (e) time windows. The selection of the early time window was varied (1 to 5) while the late window was always the last one (6 th ). Right part of the figure: Related plots for moments, A (top row, panel 7), m 1 and V (bottom row, panels 6 and 7). The plot in each panel is normalized to its maximum value, the relative color bar in each row is valid for all panels in that row.
Taking the ratio of photon counts in late and early time windows is known to suppress the influence of absorption changes in superficial layers of a turbid medium [24,56]. Figure 4 (bottom row, panel 1) shows the sensitivity profile for the ratio of photon counts (difference of attenuation) of the latest and earliest time windows (2.5 ns to 3 ns and 0 to 0.5 ns, respectively). In the same row, panels 2-5, the early time window was chosen to be the 2nd to the 5th while the late window was kept the same. Interestingly, the region of highest sensitivity moves deeper for later "early" time windows, together with progressively detaching from the surface. The attenuation differences for pairs of time windows (bottom row, panels 1-5) exhibit a similar behavior as mean time of flight and variance (panels 6 and 7), including the appearance of a region of negative sensitivity near the surface between source and detector. This behavior is a consequence of taking the difference of quantities with different spatial sensitivity profiles. By contrast, other features of time-window and moment measurands are distinctly different, in particular, the dependence on source-detector separation (see following section) and the influence of a finite IRF (see Section 3.3).

Depth-dependent sensitivity and depth selectivity
The sensitivities S M L (z) to a small absorption change in an infinitesimal layer at depth z are presented for all measurands in Fig. 5(a,b). They can be understood as the result of summing up 2D sensitivities like those illustrated in Fig. 4 additionally along the x direction. The first observation is that all z-dependent sensitivities are positive, i.e. for measurands that show regions of negative sensitivities (m 1 , V, A l -A e ) in Fig. 4 these cancel out with positive sensitivities at the same depth z. In other words, for any layer parallel to the surface of the semi-infinite medium a small homogeneous increase in absorption throughout this layer leads to a decrease in photon counts in total (N tot ) and in all time windows (N k ) individually, a decrease in the late-to-early ratio N l / N e , as well as a decrease in m 1 and V. It should be noted that the influence of a finite IRF can invalidate this rule for the time-window measurands, see Section 3.3.1.
Since the sensitivities for the three moments considered have different units they were compared in Fig. 5(a) after normalizing them to their maxima for ρ = 3 cm. They exhibit the well-known behavior: Attenuation A has the most superficial sensitivity while m 1 and V probe progressively deeper regions. Another important feature of m 1 and V is their low sensitivity to superficial regions which makes them suitable for depth-selective measurements. For moments, the sensitivity profiles clearly depend on the source-detector separation. The maximum sensitivity is smaller and located closer to the surface for ρ = 2 cm in comparison to ρ = 3 cm. Interestingly, the shape of the sensitivity profile of V at ρ = 2 cm is nearly identical to that of m 1 at ρ = 3 cm. The results shown here are consistent with (ρ-dependent) sensitivity factors for moments presented, e.g., in [46] based on Monte Carlo simulations for layered media. The fact that Monte Carlo simulations are more exact in the vicinity of source and detector has no noticeable impact on the profiles shown. Figure 5(b) compares the z-dependent sensitivities for attenuation in time windows and attenuation differences between time windows. They are displayed for a single ρ only since they are all virtually independent of ρ [38]. With increasing time (window number k), the maximum of the sensitivity moves to greater depth z and grows in magnitude. This behavior can be explained by the monotonous growth of TMPP in a layer at depth z with time [37]. At shallow depth, the sensitivity is almost independent of time (and k) while at greater depth (e.g. z = 2 cm) it is essentially zero for early photons and grows for greater times (k > 3), and in between it increases nearly linearly with k.
The sensitivities for the attenuation differences in a late and an early time window exhibit a substantially different behavior. The maximum of the sensitivity moves toward greater depth z with increasing k (time) of the early time window while keeping the late window at l = 6. The maximum magnitude of the sensitivity decreases for later "early" time windows, remaining always smaller than the sensitivity for the late window (A 6 ) alone. The low sensitivity to superficial regions implies a good depth selectivity.
By means of depth selectivity, all moments and time-window measurands are compared in Fig. 5(c,d) on the same dimensionless scale. The selectivity was obtained as the ratio of the sensitivity in a deep layer of 5 mm width and variable position in depth z L to that in a superficial layer extending from z = δ to z = 5 mm. This choice was made in view of fNIRS measurements on the head, with hemodynamic changes occurring in the cerebral cortex and in the skin. The plots for moments in Fig. 5(c) can be understood as a result of integrating the spatially resolved selectivity (see Fig. 3) over the deep layer. The highest selectivity is found for variance at ρ = 3 cm, with the deep layer ranging from 9 mm ≤ z < 14 mm. With a magnitude of 5.2, the values for V at ρ = 2 cm as well as for m 1 at ρ = 3 cm are exceeded by almost a factor of 2. The selectivity for A stays below 1, confirming that CW measurements at a single distance are not capable of discriminating between deep and superficial changes. Attenuation differences between time windows, with the early time window at a rather late position (k > 3), can achieve a depth selectivity superior to attenuation in individual time windows (see Fig. 5(d), maximum 8.5 for A 6 -A 5 ) as well as to moments, see Fig. 5(c). Figure 6 displays descriptors for the extension and maximum position of the z-dependent sensitivity S M L (z) of the various measurands to layered absorption changes. In case the sensitivity at small z is below 50% of the maximum, an FWHM (full width at half maximum) can be calculated, i.e. a measure of "depth resolution" (spatial resolution in z direction) can be obtained in a similar way as a lateral resolution in diffuse optical imaging [36]. It depends on the specific measurand and for moments on source-detector separation. For all measurands discussed here and for the optical properties and ρ values chosen, the FWHM in z stays between 1 cm and 1.4 cm for the moments and between 1.2 cm and 1.6 cm for the time-window measurands. It should be noted that these specific values are valid for the fixed set of optical properties assumed in our simulations. For moments, the qualitative behavior as discussed with Fig. 5(a) remains the same while, e.g., the maximum position moves deeper for lower absorption as well as for lower scattering [9]. For time windows, the results of the study by Del Bianco [38] are relevant: Time-dependent contrast and thus sensitivity for time windows is independent of ρ as well as of µ a . The influence of µ s can be described by an approximate scaling relationship with time, i.e. the z-dependent sensitivity profile approximately depends on t /µ s . Another limitation of our analysis is that it refers to homogeneous layer-wide absorption changes. Different results are obtained for localized absorption changes for which the full 3D sensitivity profile is relevant. Scans with small absorbers [67,68] can also reveal regions with negative sensitivity.

Time windows
An important feature of time-domain fNIRS measurements is the independence of time-dependent contrast as well as the related TMPP for small absorption changes from source-detector separation as well as from background absorption [38]. However, this independence may be substantially violated when measurements with an IRF of non-negligible width are considered. Moreover, depth-dependent contrast may be considerably impaired [36].
We illustrate the IRF dependence of sensitivities of the time-window measurands by using two examples of measured IRFs from [36]. IRF A was obtained with a laboratory setup based on a supercontinuum laser (SC500-6, Fianium Ltd, UK) with acousto-optic tunable filter, a hybrid detector module with GaAs cathode HPM-100-50 (Becker&Hickl GmbH, Germany) connected to a TCSPC module SPC-130 (Becker&Hickl GmbH, Germany) and a 2 m long step-index detection fiber (NA 0.39). IRF B is related to a multi-channel time-domain optical brain imager developed for clinical applications, with upgraded detector modules (see Refs. 30-32 in [36]). This system relied on a picosecond diode laser Sepia (PicoQuant GmbH, Germany), a compact photomultiplier module H7422-50 (Hamamatsu Photonics, Japan) with GaAs cathode, SPC-130 modules and a 1.5 m long detection fiber bundle (NA 0.54). For details of the instruments and their performance characterization by phantom-based test measurements see [36]. The shape of the IRF is mainly determined by laser pulse width, time resolution of the detector in conjunction with TCSPC and dispersion in (step-index) fibers or fiber bundles. All these factors differ for both cases and explain the largely different shape and FWHM (IRF A: 225 ps, IRF B: 580 ps) of the IRF profiles.
Both IRFs are illustrated in Fig. 7(a,b). Apart from the smaller width, IRF A has a clean shape with an approximately exponential tail while IRF B exhibits a pronounced afterpeak about 2.5 ns after the main peak, a specific feature of the GaAs photomultiplier module. The effect of convolution is shown for the (unperturbed) DTOFs for the two source-detector separations considered. For IRF A, the convolution merely acts as a minor broadening of the DTOF and delay of the tail. But for IRF B, the imprint of the afterpeak is visible in the convolved DTOFs, the wider the DTOF, the smoother the shoulder appears. Figure 7(a,b) also have the time windows of 500 ps width marked. Prior to defining them, a time origin (t = 0) was identified for each IRF as the barycenter in the interval between the points at half maximum, in accordance with the procedure applied previously [36]. The influence of the afterpeak is evident in the 6 th window from 2.5 ns to 3 ns (marked). Figure 8 and Fig. 9 report the influence of these IRFs on the sensitivities of time-window measurands to small absorption changes. Figure 8 illustrates the consequences for the spatial sensitivity profiles for attenuation in time windows (top row) and attenuation differences between late and early windows (bottom row). For ρ = 2 cm (Fig. 8(a)) the impact on A 6 is striking. The maximum sensitivity is located much shallower than for A 5 , and the profile rather resembles that for A 2 , dominating over the sensitivity to deeper regions that merely appears as a slight shadow. The reason is that the convolution with a relatively extended IRF can be regarded as a transfer of early photons to late times where the DTOF itself has already decayed to ∼1%. The resulting effect on the sensitivity related to the difference A 6 -A 5 is devastating, the sensitivity it is mostly deeply negative and not useful anymore. For ρ = 2 cm (Fig. 8(b)) the effect is weaker but still present, the negative region is more pronounced than without IRF influence (cf. Figure 4).   The z-dependent sensitivities to layered absorption changes can deviate considerably from the plots in Fig. 5(b). Examples for late time windows are presented in Fig. 9. Note that these sensitivities are dimensionless. To exemplify their magnitude, the effect of a homogeneous absorption change ∆µ a = 0.001 cm −1 as obtained by integration of S L A6 (z) and multiplication with ∆µ a , results in ∆A 6 = 0.057 for a delta-pulse like IRF.
The influence of IRF A remains very minor, the sensitivities for both A 6 and A 6 -A 5 are very similar for IRF A and IRF 0 and almost indistinguishable for both ρ. This is no longer the case for IRF B where the sensitivity for A 6 becomes ρ-dependent, drops and shifts toward smaller z. The x,y-integrated sensitivities for A 6 -A 5 are no longer always positive as without IRF influence. We conclude that wide IRFs with long tails or afterpeaks may prevent one from achieving reasonable depth selectivity for time-window measurands without prior deconvolution, which is usually not done because it is difficult on its own and may lead to additional artefacts.

Moments
As opposed to time-window measurands, differences of the moments considered here do not depend on the IRF. This statement is theoretically valid for an infinite integration range only, i.e. provided the integration limits encompass the complete DTOF. We illustrate the effect of incomplete integration on contrast for an absorption change in a deep layer as well as on photon noise, with the different IRFs introduced in Section 3.3.1. Poisson noise was calculated assuming a total of 10 6 counts per DTOF, by means of the expressions from [60] provided in Section 2.1. It should be noted that these equations are valid for arbitrary integration limits. Figure 10 represents the buildup of the signal change (contrast), the noise of the unperturbed signal and their ratio for the moments while increasing the upper integration limit t up . Figure 10(ac) detail the decrease in N tot , by plotting the relative change ∆N tot / N tot,0 ≈ -∆A, and the decrease in m 1 and V. The related evolution of noise and contrast-to-noise ratio are shown in Fig. 10(d-f) and Fig. 10(g-i), respectively.
The contrast is zero initially, i.e. for early photons, and builds up while the integration range covers an increasing fraction of the DTOF, see Fig. 10(a-c). The slow convergence that occurs already for IRF 0 can be understood as a consequence of the truncation error which we had studied previously for absolute moments [60]. Since the DTOFs become wider upon convolution with an IRF, the convergence of contrast to the true value occurs the later the wider the IRF is. This influence of the IRF on the contrast buildup is most pronounced for V (see Fig. 10(c)) which is evident from the illustration of the integrand for V in Fig. 1(b). But eventually, after the DTOF has decayed completely, the true value is approximated in all cases. However, for late times, the further gain in contrast remains minor while noise still increases substantially, in particular for V and IRF B due to its slowly decaying tail and pronounced afterpeak. Therefore, it is advisable to consider the contrast-to-noise ratio that may exhibit a maximum before it approaches a constant value at large t up where the photons have ceased which contribute to signal as well as to noise. In such case, it may be reasonable to limit the integration before noise starts to dominate, i.e. to find a trade-off between incomplete contrast and acceptable noise. As an alternative, a correction procedure was introduced to eliminate such truncation error [9].
It should be noted that Poisson noise is not the only source of uncertainty in the calculation of moments. Time-independent background signals (e.g. due to dark counts or ambient light), if not correctly subtracted, lead to an over-or underestimation of moments and prevent their convergence to a constant value at large t up . Background noise remaining after correct subtraction does not distort the signal but causes a continuing increase of noise and thus a decrease in contrast-to-noise ratio with t up even after complete decay of the DTOF. These effects should also be considered when selecting integration limits for experimental data.

Conclusions
The present work contributes to an optimization of time-domain fNIRS of the brain by a systematic and quantitative characterization of different methods to analyze measured timeof-flight distributions of photons. We considered moments, photon counts in individual time windows and ratios of photon counts in different time windows. First, 3D spatial maps of sensitivity as well as of depth selectivity were analyzed, based on perturbation simulations of changes of time-resolved diffuse reflectance due to localized absorption changes in a homogeneous semi-infinite medium. Then a quantitative comparison was performed for depth(z)-dependent sensitivity and depth selectivity for layered absorption changes. To study the robustness of the different measurands under various experimental conditions, we employed realistic IRFs in the simulations.
While we previously presented experimental results on contrast and depth selectivity for moments and photon counts in time windows in the context of performance assessment of instruments by dedicated phantom measurements in [36], ratios of photon counts in time windows are here included in the comparison for the first time. Moreover, a direct comparison of 3D spatial sensitivity profiles and depth-dependent sensitivity for moments and time-window measurands has not been reported before. Another novel aspect of the present study is the consideration of noise for changes in moments in the context of optimization of the integration interval in the presence of a non-ideal IRF.
We decided to use a simple model and straightforward perturbation simulations and to focus on layered absorption changes to illustrate the essential features of our approach to assess depth selectivity of analysis methods. Our methodology can be potentially extended to more complex cases, i.e. heterogeneous media or larger changes of the optical properties. We did not elaborate on the influence of the background optical properties which has been done by other authors for time-resolved reflectance [38] and moments [9].
The comparison of various features of the moments and time-window methods is summarized in Table 1. Each of these approaches has specific advantages. In conclusion, there is not "the optimum" quantity but rather a recommendation when to use time windows (for highest contrast) and moments (for wide IRF and when timing is unstable). Our methodology can be employed to study other measurands, as, in particular, Mellin-Laplace moments [48], modified time windows [45] or combinations of measurements at different sourcedetector separations [69,70]. The finding that ratios of photon counts in two late time windows provide superior depth selectivity could advantageously be applied in fast-gated detection schemes with small source-detector separation [71], in particular in non-contact fNIRS [33]. In that case, only the late part of DTOFs is recorded which impedes the use of moments. Conversely, late photons are available with considerably increased signal-to-noise ratio.
We add that, for a more comprehensive understanding, the findings on depth selectivity need to be complemented by an analysis of the influence of (photon) noise and the related contrast-to-noise ratio. This analysis is part of a separate study [72].

Disclosures
The authors declare no conflicts of interest.