## 1.

## Introduction

Near infrared spectroscopy (NIRS) is a noninvasive optical technique to study hemodynamic changes and oxygenation in tissue. It utilizes differences in the absorption spectra of oxygenated and deoxygenated hemoglobin and the relatively good transparency of tissues for light in the near-infrared region.^{1}^{,}^{2} NIRS is based on illumination of the tissue and detection of diffusely remitted light in a spot located a few centimeters apart from the source position. The NIRS technology has been mostly used in muscle^{3}4.^{–}^{5} and brain studies.^{6}7.8.^{–}^{9} Regarding the assessment of brain hemodynamic changes and oxygenation in the brain, there is a growing field of clinical applications.^{4}^{,}^{6}^{,}^{7}^{,}^{10}11.12.13.14.15.16.^{–}^{17}

In the early stage of development, NIRS was used for evaluation of changes of oxy- and deoxyhemoglobin in the brain of neonates.^{18}19.20.21.22.^{–}^{23} With advances of light detection efficiency the technique was applied in adult humans, as well. Many authors reported evidence that NIRS signals do contain information originating from the cerebral cortex.^{9}^{,}^{10}^{,}^{24} A major obstacle in routine application of NIRS instruments in clinical practice is the contamination of the measured signals of changes of oxy- and deoxyhemoglobin in the brain by changes in extracerebral tissue.^{25}26.^{–}^{27}

Several methods allowing for compensation or removal of the influence of the extracerebral tissue contamination were proposed. They are based on multi-distance continuous wave measurements,^{28}^{,}^{29} frequency-domain experiments^{30} or the analysis of the broadening of short laser pulses during their propagation in the tissue under investigation.^{25}^{,}^{27}^{,}^{31}32.33.34.^{–}^{35}

The present paper is related to time-resolved near-infrared spectroscopy which is currently extensively studied as a tool that allows separation of extra- and intracerebral changes of oxy- and deoxyhemoglobin.^{35}36.^{–}^{37} The depth-selective determination of absorption changes in the adult human head can be achieved by analysis of ratios of time windows of the distributions of times of flight (DTOF) of photons^{27}^{,}^{38}^{,}^{39} or by evaluation of their statistical moments.^{31} Previously we showed the significant advantages of using moments for assessing depth-resolved changes of absorption.^{31}^{,}^{35}^{,}^{40} A particular feature of the analysis based on moments is that the elimination of the influence of the instrumental response function is straightforward.^{41} We found that depth-resolved profiles of absorption changes can be derived from multi-distance time-resolved measurements and evaluation of the moments: integral, mean time of flight and variance.^{31}^{,}^{42} The changes in the absorption coefficient of the tissue ${\mu}_{a}$ can be connected to functional activation of the cortex and corresponding changes in the concentrations of oxy- and deoxyhemoglobin. Absorption changes also occur during a bolus of an optical contrast agent.

For our previously published study on depth-resolved assessment of absorption changes the moments at several source-detector separations were required.^{31}^{,}^{42} However, multi-distance time-resolved measurements may not always be feasible. In particular, multi-site and multi-distance continuous wave measurements allow rather high lateral spatial resolution in optical brain imaging to be achieved.^{29} However, in the case of time-resolved light detection a multi-distance approach with recording at various distances by the same detector is rather difficult because of the limited dynamic range of the fast detectors. Furthermore, time-resolved measurements that are necessary to record DTOFs and their moments, in particular variance, most commonly rely on time-correlated single photon counting^{43} and are instrumentally demanding. This technology remains expensive and is not widely used in NIRS brain studies. The less sophisticated and commercially available frequency domain technique in principle provides access to mean time of flight via phase measurements.^{44}

The present paper addresses experimental design considerations, in particular the selection of source-detector separations used, as well as the available types of measurements corresponding to continuous wave, frequency or time domain methods, represented by the various subsets of moments. We compare the various approaches with respect to the standard deviation of the obtained absorption changes $\mathrm{\Delta}{\mu}_{a}$ in a two-layered tissue model due to photon noise. We assume the thickness of the upper layer to be known. The upper layer corresponds to extracerebral tissue (skin, skull) whereas the lower layer mimics the brain tissue. For this model, the number of measured quantities, i.e., the three moments recorded at one or more distances, exceeds the number of unknown quantities to be determined, i.e., the changes of the absorption coefficients in the two layers. The goal is to investigate whether the number of different measurements can be significantly reduced without substantially increasing the standard deviation of the estimated absorption changes in the two layers. To this aim we study the standard deviations of estimates of $\mathrm{\Delta}{\mu}_{a}$ varying the number of measured quantities and the source-detector separations.

## 2.

## Theory

The equations that will be derived in this section are related to the general case of a medium with $J$ compartments. The specific results shown in Sec. 4 refer to a two-layered medium, i.e., a superficial layer with thickness $d$ on top of a semi-infinite compartment (see Fig. 1). This model mimics the layered structure of the human head consisting of extracerebral tissue (skin, skull, cerebro-spinal fluid) and brain tissue. We assume that changes in DTOFs are mainly due to changes in the absorption coefficients of the tissues $\mathrm{\Delta}{\mu}_{a}$ whereas the scattering properties remain constant ($\mathrm{\Delta}{\mu}_{s}^{\prime}=0$).

The time-of-flight distributions are obtained as histograms of photon counts ${N}_{i}$ in time channels $i$. We start from the following statistical moments of measured DTOFs, the 0th moment or total photon count

the first moment ${m}_{1}$ or mean time of flight $T$ and the second centralized moment or variance where is the second moment of the DTOF. Here and in the following the sum is extended over a finite number of time channels of the DTOF from $i=1$ to $i={i}_{\mathrm{max}}$.The following depth-resolved analysis is based on the changes in attenuation $\mathrm{\Delta}A=-\mathrm{ln}({N}_{\mathrm{tot}}^{*}/{N}_{\mathrm{tot}})$, mean time of flight $\mathrm{\Delta}T={T}^{*}-T$ and variance $\mathrm{\Delta}V={V}^{*}-V$ where the symbols with star denote the quantity after an absorption change.

Assuming that the changes of the absorption coefficient are small, the relationship between $\mathrm{\Delta}{\mu}_{a}$ in the individual compartments of the model and changes of moments of DTOFs is linear,

See Ref. 31. The vectors $\mathbf{\Delta}\mathbf{A}$, $\mathbf{\Delta}\mathbf{T}$ and $\mathbf{\Delta}\mathbf{V}$ contain changes of attenuation, mean time of flight and variance of the DTOFs at $M$ source-detector separations $\mathbf{r}=({r}_{1},\dots ,{r}_{M})$, i.e., $\mathbf{\Delta}\mathbf{A}=\phantom{\rule{0ex}{0ex}}{(\mathrm{\Delta}{A}_{1},\dots ,\mathrm{\Delta}{A}_{M})}^{T}$, $\mathbf{\Delta}\mathbf{T}={(\mathrm{\Delta}{T}_{1},\dots ,\mathrm{\Delta}{T}_{M})}^{T}$, $\mathbf{\Delta}\mathbf{V}={(\mathrm{\Delta}{V}_{1},\phantom{\rule{0ex}{0ex}}\dots ,\mathrm{\Delta}{V}_{M})}^{T}$. The vector $\mathbf{\Delta}{\mathbf{\mu}}_{\mathbf{a}}$ contains changes of ${\mu}_{a}$ in $J$ different compartments of the model, i.e., $\mathbf{\Delta}{\mathbf{\mu}}_{\mathbf{a}}={(\mathrm{\Delta}{\mu}_{a1},\dots ,\phantom{\rule{0ex}{0ex}}\mathrm{\Delta}{\mu}_{aJ})}^{T}$.The sensitivity matrices **MPP** (mean partial pathlength), **MTSF** (mean time of flight) and **VSF** (variance of the DTOF) are $M\times J$ matrices that represent the sensitivity of the measured quantities to small changes of the absorption coefficients for the various source-detector separations and compartments of the medium. These sensitivities can be obtained by Monte Carlo simulations^{31} or by applying a perturbation method based on the diffusion equation.^{37}^{,}^{45}^{,}^{46} For a measurement carried out at $M$ source-detector separations the system of equations given by Eqs. (5)–(7) consists of $3\xb7M$ equations with $J$ unknowns, i.e., the changes of the absorption coefficient $\mathbf{\Delta}{\mathbf{\mu}}_{\mathbf{a}}$ in the $J$ compartments of the model. Let $\mathbf{\Delta}\hat{\mathbf{A}}$, $\mathrm{\Delta}\hat{\mathbf{T}}$ and $\mathrm{\Delta}\hat{\mathbf{V}}$ denote the estimates for $\mathbf{\Delta}\mathbf{A}$, $\mathrm{\Delta}\mathbf{T}$ and $\mathrm{\Delta}\mathbf{V}$, i.e., the measured data. An estimate $\mathbf{\Delta}{\hat{\mathbf{\mu}}}_{\mathbf{a}}$ may then be determined by the weighted least-squares approach according to^{47}

## (8)

$$\mathbf{\Delta}{\hat{\mathbf{\mu}}}_{\mathbf{a}}={({\mathbf{X}}^{T}{\mathbf{Z}}^{-1}\mathbf{X})}^{-1}{\mathbf{X}}^{T}{\mathbf{Z}}^{-1}\left(\begin{array}{c}\mathbf{\Delta}\hat{\mathbf{A}}\\ \mathbf{\Delta}\hat{\mathbf{T}}\\ \mathbf{\Delta}\hat{\mathbf{V}}\end{array}\right),$$## (9)

$$\mathbf{X}=\left(\begin{array}{c}\mathbf{MPP}\\ \mathbf{MTSF}\\ \mathbf{VSF}\end{array}\right)$$## (10)

$$\mathbf{Z}=\left(\begin{array}{ccc}\mathrm{cov}[\mathbf{\Delta}\hat{\mathbf{A}},\mathbf{\Delta}\hat{\mathbf{A}}]& \mathrm{cov}[\mathbf{\Delta}\hat{\mathbf{A}},\mathbf{\Delta}\hat{\mathbf{T}}]& \mathrm{cov}[\mathbf{\Delta}\hat{\mathbf{A}},\mathbf{\Delta}\hat{\mathbf{V}}]\\ \mathrm{cov}[\mathbf{\Delta}\hat{\mathbf{T}},\mathbf{\Delta}\hat{\mathbf{A}}]& \mathrm{cov}[\mathbf{\Delta}\hat{\mathbf{T}},\mathbf{\Delta}\hat{\mathbf{T}}]& \mathrm{cov}[\mathbf{\Delta}\hat{\mathbf{T}},\mathbf{\Delta}\hat{\mathbf{V}}]\\ \mathrm{cov}[\mathbf{\Delta}\hat{\mathbf{V}},\mathbf{\Delta}\hat{\mathbf{A}}]& \mathrm{cov}[\mathbf{\Delta}\hat{\mathbf{V}},\mathbf{\Delta}\hat{\mathbf{T}}]& \mathrm{cov}[\mathbf{\Delta}\hat{\mathbf{V}},\mathbf{\Delta}\hat{\mathbf{V}}]\end{array}\right).$$The covariance matrix of $\mathbf{\Delta}{\hat{\mathbf{\mu}}}_{\mathbf{a}}$ is then obtained as

## (11)

$$\mathrm{cov}(\mathbf{\Delta}{\hat{\mathbf{\mu}}}_{\mathbf{a}},\mathbf{\Delta}{\hat{\mathbf{\mu}}}_{\mathbf{a}})={({\mathbf{X}}^{T}{\mathbf{Z}}^{-1}\mathbf{X})}^{-1}.$$*j*-th compartment. Note that care needs to be taken in the calculations when ${\mathbf{Z}}^{-1/2}\mathbf{X}$ is nearly rank-deficient.

^{47}

In the following we give expressions for the covariance matrix $\mathbf{Z}$ (for details see Appendix), assuming that the (dominant) photon noise follows a Poisson distribution.

## (12)

$$\mathrm{cov}{(\mathbf{\Delta}\hat{\mathbf{A}},\mathbf{\Delta}\hat{\mathbf{A}})}_{k,l}=\frac{2}{{N}_{\mathrm{tot},k}}{\delta}_{k,l}$$## (13)

$$\mathrm{cov}{(\mathbf{\Delta}\hat{\mathbf{T}},\mathbf{\Delta}\hat{\mathbf{T}})}_{k,l}=2\frac{{V}_{k}}{{N}_{\mathrm{tot},k}}{\delta}_{k,l}$$## (14)

$$\mathrm{cov}{(\mathbf{\Delta}\hat{\mathbf{V}},\mathbf{\Delta}\hat{\mathbf{V}})}_{k,l}=2\frac{{m}_{4,k}^{c}-{V}_{k}^{2}}{{N}_{\mathrm{tot},k}}{\delta}_{k,l}$$## (15)

$$\mathrm{cov}{(\mathbf{\Delta}\hat{\mathbf{A}},\mathbf{\Delta}\hat{\mathbf{T}})}_{k,l}=\mathrm{cov}{(\mathbf{\Delta}\hat{\mathbf{T}},\mathbf{\Delta}\hat{\mathbf{A}})}_{k,l}=0$$## (16)

$$\mathrm{cov}{(\mathbf{\Delta}\hat{\mathbf{A}},\mathbf{\Delta}\hat{\mathbf{V}})}_{k,l}=\mathrm{cov}{(\mathbf{\Delta}\hat{\mathbf{V}},\mathbf{\Delta}\hat{\mathbf{A}})}_{k,l}=0$$## (17)

$$\mathrm{cov}{(\mathbf{\Delta}\hat{\mathbf{T}},\mathbf{\Delta}\hat{\mathbf{V}})}_{k,l}=\mathrm{cov}{(\mathbf{\Delta}\hat{\mathbf{V}},\mathbf{\Delta}\hat{\mathbf{T}})}_{k,l}=2\frac{{m}_{3,k}^{c}}{{N}_{\mathrm{tot},k}}{\delta}_{k,l},$$## 3.

## Implementation

For the following application we focused on the particular case of $J=2$ layers, i.e., we assumed that the tissue consists of an upper and a lower layer which represent extra- and intracerebral tissue compartments. The sensitivity matrix $\mathbf{X}$ as well as the covariance matrix $\mathbf{Z}$ depend on the background optical properties. We presumed a semi-infinite homogeneous turbid medium with absorption coefficient ${\mu}_{a}=0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, reduced scattering coefficient ${\mu}_{s}^{\prime}=1{\text{\hspace{0.17em}}\text{\hspace{0.17em}}mm}^{-1}\text{\hspace{0.17em}}$ and refractive index $n=1.4$ which are values typically taken for simple head models.^{25}^{,}^{31}^{,}^{48}

The sensitivity matrix $\mathbf{X}$ [see Eq. (9)] was determined by Monte-Carlo simulations. Details of the Monte Carlo code for a layered turbid medium used for these calculations were described elsewhere.^{25}^{,}^{31} For the calculation of sensitivity factors the medium was divided into nine layers each 2 mm thick on top of a semi-infinite bottom layer. The sensitivity matrices **MPP**, **MTSF**, and **VSF** were calculated for all layers of the model according to algorithms given in.^{31} In these calculations the DTOFs were sampled in ${i}_{\mathrm{max}}=150$ time channels, each 25 ps wide. For the two-layered model we summed up the values of **MPP**, **MTSF**, and **VSF** obtained for an upper and a lower group of layers of the 10-layered structure. This approach allows us to perform the analysis for various values of the thickness $d$ of the upper layer. Monte Carlo calculations were carried out for a total of ${5\xb710}^{8}$ photon packages detected at source–detector separations between 5 and 62 mm. The DTOFs, as well as the sensitivities, were calculated simultaneously for a set of 19 concentric, consecutive, ring-shaped detectors, each 3 mm wide. The source-detector separations mentioned below in the description of results refer to the inner diameter of the detector rings.

The calculation of the $\mathbf{Z}$ matrix [see Eq. (10)] requires the knowledge of the moments ${V}_{k}$, ${m}_{3,k}^{c}$ and ${m}_{4,k}^{c}$ as well as the total photon count ${N}_{\mathrm{tot},k}$ [see Eqs. (12)–(17)]. For each source-detector separation ${r}_{k}$ the higher moments were derived from the DTOFs as obtained from the Monte-Carlo simulations described above. Regarding ${N}_{\mathrm{tot}}$, we introduced experimentally realistic conditions. We assumed the total number of detected photons to be ${N}_{\mathrm{tot},k}={N}_{\mathrm{tot}}={10}^{6}$, independently of the inter-optode distance, corresponding to a measurement at a count rate of 1 MHz with an acquisition time of 1 s. These conditions are typical for measurements based on state-of-the-art time-correlated single photon counting technology that is capable of recording at count rates as high as a few MHz.^{43}

The covariance matrix in Eq. (11) for the least-squares estimate $\mathbf{\Delta}{\hat{\mathbf{\mu}}}_{\mathbf{a}}$ depends only on the covariance matrix $\mathbf{Z}$ of the data Eq. (10) and on the design of the experiment. It does not, however, depend on the data, as the problem is linear. We can therefore calculate the covariance matrix of the least-squares estimate without carrying out the experiment, provided that the covariance matrix $\mathbf{Z}$ of the data is available (or can be calculated) as well as the sensitivity matrix $\mathbf{X}$ in Eq. (9). Both can be (approximately) calculated using the Monte Carlo simulations as described above. In doing so, we were able to investigate the effect of the selection of source-detector separations and of the number of measured quantities on the accuracy of the obtained estimate of the absorption change. Subsequently we present the obtained standard deviations of these estimates in dependence on the selected source-detector separations and the chosen number of measured quantities.

## 4.

## Results

This section presents the results for the standard deviations $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,j})$ of the absorption changes in the two layers of the turbid medium. We start from a dataset $\mathbf{\Delta}\hat{\mathbf{A}}$, $\mathrm{\Delta}\hat{\mathbf{T}}$ or $\mathrm{\Delta}\hat{\mathbf{V}}$ related to $\mathrm{\Delta}A$, $\mathrm{\Delta}T$ or $\mathrm{\Delta}V$, each obtained at four source-detector separations (8, 17, 26, and 35 mm). In such case the system of equations [Eqs. (5)–(7)] contains 12 equations to evaluate the two unknowns $\mathrm{\Delta}{\mu}_{a,1}$ and $\mathrm{\Delta}{\mu}_{a,2}$ and is highly over-determined. Hence the number of input quantities can be reduced. In particular, several scenarios were studied in which a subset of measurements was selected from the full data set. All calculations of $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,j})$ were repeated with sensitivity matrices corresponding to various values of the thickness $d$ of the upper layer.

Figure 2 displays the results obtained for various combinations of moments when all four inter-optode distances $r$ were included. The calculated standard deviations $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,2})$ and $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,1})$ of the changes of the absorption coefficient of the lower and upper layer, respectively, are shown as a function of thickness $d$ of the upper layer. The results for the lower layer (index 2) are always shown first since this case is of major importance in brain imaging. Note that all standard deviations are given in absolute units. Here and within each subsequent figure, the scales of all panels are the same.

In general, $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,2})$ increases and $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,1})$ decreases for increasing thickness $d$. This observation is connected to a “partial volume” effect, i.e., the relative magnitude of the contribution of the absorption change in the respective layer to the measured signal change. Comparing the results for single moments [$\mathrm{\Delta}A$, $\mathrm{\Delta}T$ or $\mathrm{\Delta}V$, see Fig. 2(a) and 2(c)], $\mathrm{\Delta}A$ is definitely most suitable when analyzing the upper layer [Fig. 2(c)], whereas $\mathrm{\Delta}T$ and $\mathrm{\Delta}V$ become similarly advantageous when determining $\mathrm{\Delta}{\mu}_{a,2}$, in particular for larger $d$. This finding is in line with the different depth dependence of the sensitivity factors of the various moments^{31} that also changes with the inter-optode distance $r$.

As can be seen in Fig. 2(b) and 2(d), in all cases $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,j})$ is reduced when including higher moments. In particular, we compared the standard deviations obtained for the following three cases, (i) $\mathrm{\Delta}A$ only, measured at four source-detector separations, (ii) $\mathrm{\Delta}A$ combined with $\mathrm{\Delta}T$, (iii) the full set $\mathrm{\Delta}A$, $\mathrm{\Delta}T$ and $\mathrm{\Delta}V$ at four source-detector separations is taken into account. As an example, for $d=12\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$, a realistic thickness of the extracerebral layer in an adult, the inclusion of $\mathrm{\Delta}T$ improves the precision of $\mathrm{\Delta}{\mu}_{a,2}$ by a factor of 2.5. Additional inclusion of $\mathrm{\Delta}V$ yields in another improvement by about 4%.

Next, we analyzed the influence of the number of source-detector separations at which the moments were measured, on $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,j})$. As a start, a single detector at a large distance ($r=35\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$) was chosen, and shorter source-detector separations were incorporated step by step. The data set containing changes in two moments of DTOFs, i.e., $\mathrm{\Delta}A$ and $\mathrm{\Delta}T$ was considered for this analysis. The results presented in Fig. 3 confirm that the increase of the number of source-detector separations always leads to a reduction of $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,j})$. The influence of the inclusion in the analysis of data acquired at shorter $r$ is particularly small for the deeper layer in case of large $d$ [see Fig. 3(a)]. This finding is obvious, since detection of moments at small distances $r$ is rather insensitive to deep absorption changes. Analysis of the analogous data set containing all three moments of DTOFs, i.e., $\mathrm{\Delta}A$, $\mathrm{\Delta}T$ and $\mathrm{\Delta}T$, showed a similar pattern of changes in $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,j})$ as a function of $d$ for the various combinations of source-detector separations.

In order to study the potential of a single-distance time-resolved measurement, the analysis of $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,j})$ was performed for a single source-detector separation and various combinations of moments (see Fig. 4). In general, the best option is always to take all three moments into account. But any combination of two of them is sufficient to retrieve the two unknown absorption changes. The exclusive use of changes of mean time of flight and variance ($\mathrm{\Delta}T$ and $\mathrm{\Delta}V$) results in high $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,j})$, in particular for the upper layer [see Fig. 4(b)]. This effect can be explained by the insensitivity of these two moments to superficial absorption changes^{31} which is a consequence of their normalization to ${N}_{\mathrm{tot}}$ by definition [see Eqs. (2) and (3)]. Comparison of the results for the combination $\mathrm{\Delta}A$, $\mathrm{\Delta}V$ with the results for all moments $\mathrm{\Delta}A$, $\mathrm{\Delta}T$, $\mathrm{\Delta}V$ shows that both options lead to almost identical $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,j})$ at large thickness $d$ [see Fig. 4(a) and 4(b)]. In contrast, at small $d$ the values of $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,j})$ converge for the combinations $\mathrm{\Delta}A$, $\mathrm{\Delta}T$ and $\mathrm{\Delta}A$, $\mathrm{\Delta}T$, $\mathrm{\Delta}V$.

Since the investigation of cortical activation is the major goal of time-domain brain imaging, the retrieval of the absorption changes in the lower layer from measurements at a single distance is of particular interest. Therefore, the analysis illustrated in Fig. 4(a) was extended to other source-detector separations (see Fig. 5). As expected, the absorption change in the lower layer is advantageously determined at large source-detector separations. Comparing Fig. 5(a)–5(d), it can be seen that $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,2})$ increases by 1 to 2 orders of magnitude with decreasing source-detector separation. Moreover, the range of $\sigma (\mathrm{\Delta}{\hat{\mu}}_{a,2})$ between small (2 mm) and large (18 mm) thickness $d$ is much smaller (one order of magnitude) for larger inter-optode distance (35 mm) than for short source-detector separation (up to two orders of magnitude for 8 mm). That means the thicker the upper layer, the more important it is to increase $r$. Nevertheless, the relative order of curves obtained for different combinations of moments, as well as the improvement when including more input quantities, are generally preserved for different source-detector separations.

## 5.

## Discussion and Conclusions

The paper is devoted to the determination of absorption changes in several tissue compartments that are derived from various moments of DTOFs measured at a single or multiple distances. The reconstruction of the absorption changes that are assumed to be small is performed by solving the corresponding system of linear equations. The calculation of the standard deviations of the results is done by propagating the standard deviations of the measurements that are not necessarily statistically independent. In particular, various moments derived from the same measured DTOF may be correlated, and their nonvanishing covariances have to be included. The calculation of the covariance matrix of the absorption changes according to Eq. (11) accounts for such correlations. Our analysis may serve as an example for the application of this approach to similar problems, e.g., the reconstruction of absorption changes based on photon counts in time windows of measured DTOFs.

Necessary pieces of information for the analysis performed in this paper are the sensitivities that relate the measured changes in moments to the underlying absorption changes, and covariance matrices of the measurements. It is important to note that due to the linearity of the problem the knowledge of the magnitude of the absorption changes in the tissue compartments is not required for the calculation of their standard deviations. Thus the standard deviations calculated should apply to any absorption changes, as long as the assumption of linearity is not violated.

The sensitivities depend on the geometry, in particular on the thickness of the upper layer and the source-detector separation, as well as on the background optical properties of the medium. The sensitivity factors used in the analysis were obtained from Monte Carlo simulations (see Ref. 31). Expressions to calculate the covariance matrix of the measurements, i.e., attenuation, mean time of flight and variance of the DTOF, were derived based on the assumption of Poisson-distributed photon noise of the DTOF. They contain moments of the DTOF up to the 4th order. All standard deviations of $\mathrm{\Delta}{\mu}_{a}$ shown above were calculated assuming the same number (${10}^{6}$) of photons acquired at all source-detector separations. Since all elements of the covariance matrix $\mathit{Z}$ scale with $1/{N}_{\mathrm{tot}}$, the results can be simply rescaled to other common values of ${N}_{\mathrm{tot}}$. Individual integral photon counts for different source-detector separations can also be taken into account easily.

As to the absolute magnitude of the standard deviations obtained in the analysis presented above, the optimal values range between ${10}^{-5}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and ${10}^{-4}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ in most cases. With an assumed background absorption coefficient of $0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, a signal-to-noise ratio of 1 would be attained with a relative absorption change $\mathrm{\Delta}{\mu}_{a}/{\mu}_{a}$ between 0.1% and 1%.

Most of the results were presented as a function of the thickness of the upper layer $d$ that was assumed to be known. This is realistic if anatomical images are available. In principle, $d$ could also be estimated from the measurements. However, the dependence of the measurands on $d$ is nonlinear. In such a case the linear approach outlined above would not apply. In particular, the standard deviations of the estimates would depend on the estimates themselves.

Our analysis has been restricted to the consideration of standard deviations due to photon noise. The major factor that limits the accuracy of the quantification of absorption changes is the insufficient knowledge of the optical properties of the various tissues of the head that are required to calculate the sensitivities. In addition, systematic deviations can originate from the necessary limitation of the integration range in the calculation of moments, especially of higher order, from measured DTOFs.

For the practical application of the approach presented here it is important to discuss the influence of an instrument response function (IRF) of finite width, which has so far been neglected. The results of the measurement, i.e., the absorption changes themselves, do not depend on the IRF since the changes in mean time of flight as well as in variance are virtually independent of the IRF.^{41} On the contrary, the elements of the covariance matrix $Z$ of the input measurands contain (absolute) moments of the DTOFs [see Eqs. (12)–(14) and (17)] that definitely depend on the IRF. Consequently, width and shape of the IRF will influence the evaluated standard deviations of absorption changes. In this sense, the values for $\sigma (\mathrm{\Delta}{\mu}_{a,1})$ and $\sigma (\mathrm{\Delta}{\mu}_{a,2})$ shown in this paper represent lower bounds. However, an in-depth quantitative study on the effect of realistic IRFs would go beyond the scope of this paper. When analyzing real measurements, the moments needed to calculate the covariance matrix $\mathbf{Z}$ can be easily obtained from the measured DTOFs.

The analysis of standard deviations of absorption changes as presented above may facilitate the design of experiments in functional brain imaging. The consideration of various subsets of moments may lead to conclusions regarding the choice of continuous-wave (cw), frequency-domain or time-domain methods. However, different sources of noise present in these techniques should be carefully considered. Moreover, the choice of source-detector separations can be optimized, and the effect of photon count rates can be studied. On the other hand, the quality of a given set of measured data can be assessed and combinations of moments can be identified for which the standard deviation of the absorption changes is lowest.

In general, inclusion of additional moments and source-detector separations leads to better estimates of the absorption changes in the two layers. Mean time of flight and variance of the DTOFs are important to estimate $\mathrm{\Delta}{\mu}_{a,2}$ which is the major measurand in brain imaging. However, a reconstruction based on $\mathrm{\Delta}T$ and $\mathrm{\Delta}V$ alone is not the optimum choice (see Fig. 4). This combination would have the advantage to be independent of amplitude and movement artifacts. But the depth sensitivity profiles of both measurands are similar, and they are statistically correlated. Another observation is that the inclusion of $\mathrm{\Delta}V$ in addition to $\mathrm{\Delta}A$ and $\mathrm{\Delta}T$ did not improve the analysis considerably [see Figs. 2(b), 2(d), 4 and 5]. Thus inclusion of variance seems to be of minor relevance. It should be noted, however, that the variance of the DTOF has several advantages that are of practical relevance for brain imaging but were not subject of the present analysis.^{37}^{,}^{40}^{,}^{49} These are primarily its depth selectivity, i.e., high sensitivity to deep and low sensitivity to superficial absorption changes, and its robustness against amplitude and timing drifts.

## Appendices

### Appendix

Here the derivation of the elements of the covariance matrix $\mathbf{Z}$ is described in more detail. We consider the $M\times M$ sub-matrices $\mathrm{cov}(\mathrm{\Delta}\hat{\mathbf{Q}},\mathrm{\Delta}\hat{\mathbf{R}})$ of $\mathbf{Z}$ where $\mathrm{\Delta}\hat{\mathbf{Q}}$ and $\mathbf{\Delta}\hat{\mathbf{R}}$ stand for the vectors $\mathbf{\Delta}\hat{\mathbf{A}}$, $\mathbf{\Delta}\hat{\mathbf{T}}$ and $\mathbf{\Delta}\hat{\mathbf{V}}$ which contain the estimates of changes in attenuation, mean time of flight and variance measured at $M$ source-detector separations. Since these $M$ measurements are statistically independent it follows that each sub-matrix is diagonal, i.e.,

## (18)

$$\mathrm{cov}(\mathrm{\Delta}{\hat{Q}}_{k},\mathrm{\Delta}{\hat{R}}_{l})=\mathrm{cov}(\mathrm{\Delta}{\hat{Q}}_{k},\mathrm{\Delta}{\hat{R}}_{l}){\delta}_{k,l}.$$## (19)

$$\mathrm{cov}(\mathrm{\Delta}{\widehat{Q}}_{k},\mathrm{\Delta}{\widehat{R}}_{k})=\mathrm{cov}[({\hat{Q}}_{k}^{*}-{\widehat{Q}}_{k}),({\hat{R}}_{k}^{*}-{\widehat{R}}_{k})]\phantom{\rule{0ex}{0ex}}=\mathrm{cov}({\hat{Q}}_{k}^{*},{\hat{R}}_{k}^{*})+\mathrm{cov}({\widehat{Q}}_{k},{\widehat{R}}_{k}).$$## (20)

$$\mathrm{cov}(\mathrm{\Delta}{\hat{Q}}_{k},\mathrm{\Delta}{\hat{R}}_{k})\approx 2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{cov}({\hat{Q}}_{k},{\hat{R}}_{k}).$$## (21)

$$\mathrm{cov}(\hat{T},\hat{T})=\mathrm{var}(\hat{T})=\sum _{i}{\left(\frac{\partial T}{\partial {N}_{i}}\right)}^{2}\mathrm{var}({N}_{i})\phantom{\rule{0ex}{0ex}}=\sum _{i}{[\frac{1}{{N}_{\mathrm{tot}}}\frac{\partial (\sum {t}_{i}{N}_{i})}{\partial {N}_{i}}\phantom{\rule{0ex}{0ex}}-\frac{1}{{{N}_{\mathrm{tot}}}^{2}}(\sum {t}_{i}{N}_{i}\left)\frac{\partial {N}_{\mathrm{tot}}}{\partial {N}_{i}}\right]}^{2}\mathrm{var}({N}_{i})\phantom{\rule{0ex}{0ex}}=\sum _{i}{(\frac{1}{{N}_{\mathrm{tot}}}{t}_{i}-\frac{T}{{N}_{\mathrm{tot}}}\xb71)}^{2}{N}_{i}\phantom{\rule{0ex}{0ex}}=\frac{1}{{N}_{\mathrm{tot}}^{2}}\sum _{i}({t}_{i}^{2}{N}_{i}-2T{t}_{i}{N}_{i}+{T}^{2}{N}_{i})\phantom{\rule{0ex}{0ex}}=\frac{1}{{N}_{\mathrm{tot}}}({m}_{2}-2{T}^{2}+{T}^{2})=\frac{V}{{N}_{\mathrm{tot}}}$$## (22)

$$\mathrm{cov}(\hat{A},\hat{T})=\sum _{i}\frac{\partial A}{\partial {N}_{i}}\frac{\partial T}{\partial {N}_{i}}\mathrm{var}({N}_{i})\phantom{\rule{0ex}{0ex}}=\sum _{i}\left(\frac{1}{{N}_{\mathrm{tot}}}\frac{\partial {N}_{\mathrm{tot}}}{\partial {N}_{i}}\right)[\frac{1}{{N}_{\mathrm{tot}}}\frac{\partial (\sum {t}_{i}{N}_{i})}{\partial {N}_{i}}\phantom{\rule{0ex}{0ex}}-\frac{1}{{{N}_{\mathrm{tot}}}^{2}}(\sum {t}_{i}{N}_{i}\left)\frac{\partial {N}_{\mathrm{tot}}}{\partial {N}_{i}}\right]\mathrm{var}({N}_{i})\phantom{\rule{0ex}{0ex}}=\frac{1}{{N}_{\mathrm{tot}}}\sum _{i}(\frac{1}{{N}_{\mathrm{tot}}}{t}_{i}-\frac{T}{{N}_{\mathrm{tot}}}){N}_{i}\phantom{\rule{0ex}{0ex}}=\frac{1}{{N}_{\mathrm{tot}}}(\frac{1}{{N}_{\mathrm{tot}}}\sum _{i}{t}_{i}{N}_{i}-\frac{T}{{N}_{\mathrm{tot}}}\sum _{i}{N}_{i})=0$$## Acknowledgments

The research leading to these results has received funding from the European Community’s Seventh Framework Programme (No. FP7/2007-2013) under Grant Agreement FP7-HEALTH-F5-2008-201076. The study is also partly financed by the Polish–German project: “Hypo- and Hyperperfusion during Subacute Stroke. Tracking Perfusion Dynamics in Stroke Patients with Optical Imaging.”

## References

S. Wrayet al., “Characterization of the near infrared absorption spectra of cytochrome aa3 and haemoglobin for the noninvasive monitoring of cerebral oxygenation,” Biochim. Biophys. Acta 933(1), 184–192 (1988).BBACAQ0006-3002http://dx.doi.org/10.1016/0005-2728(88)90069-2Google Scholar

F. F. Jobsis, “Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science 198(4323), 1264–1267 (1977).SCIEAS0036-8075http://dx.doi.org/10.1126/science.929199Google Scholar

B. Chanceet al., “Time-resolved spectroscopy of hemoglobin and myoglobin in resting and ischemic muscle,” Anal. Biochem. 174(2), 698–707 (1988).ANBCA20003-2697http://dx.doi.org/10.1016/0003-2697(88)90076-0Google Scholar

T. Hamaokaet al., “Near-infrared spectroscopy/imaging for monitoring muscle oxygenation and oxidative metabolism in healthy and diseased humans,” J. Biomed. Opt. 12(6), 062105 (2007).JBOPFO1083-3668http://dx.doi.org/10.1117/1.2805437Google Scholar

A. KienleT. Glanzmann, “In vivo determination of the optical properties of muscle with time-resolved reflectance using a layered model,” Phys. Med. Biol. 44(11), 2689–2702 (1999).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/44/11/301Google Scholar

M. WolfM. FerrariV. Quaresima, “Progress of near-infrared spectroscopy and topography for brain and muscle clinical applications,” J. Biomed. Opt. 12(6), 062104 (2007).JBOPFO1083-3668http://dx.doi.org/10.1117/1.2804899Google Scholar

E. M. Hillman, “Optical brain imaging in vivo: techniques and applications from animal to man,” J. Biomed. Opt. 12(5), 051402 (2007).JBOPFO1083-3668http://dx.doi.org/10.1117/1.2789693Google Scholar

H. ObrigA. Villringer, “Beyond the visible—imaging the human brain with light,” J. Cereb. Blood Flow Metab. 23(1), 1–18 (2003).JCBMDN0271-678Xhttp://dx.doi.org/10.1097/00004647-200301000-00001Google Scholar

A. VillringerB. Chance, “Non-invasive optical spectroscopy and imaging of human brain function,” Trends Neurosci. 20(10), 435–442 (1997).TNSCDR0166-2236http://dx.doi.org/10.1016/S0166-2236(97)01132-6Google Scholar

A. VillringerJ. SteinbrinkH. Obrig, “Editorial comment—cerebral near-infrared spectroscopy: how far away from a routine diagnostic tool?,” Stroke 35(1), 70–72 (2004).SJCCA70039-2499http://dx.doi.org/10.1161/01.STR.0000110122.57772.C3Google Scholar

C. Zweifelet al., “Continuous assessment of cerebral autoregulation with near-infrared spectroscopy in adults after subarachnoid hemorrhage,” Stroke (2010).SJCCA70039-2499Google Scholar

M. J. Herrmannet al., “Reduced prefrontal oxygenation in Alzheimer disease during verbal fluency tasks,” Am. J. Geriatr. Psych. 16(2), 125–135 (2008).http://dx.doi.org/10.1097/JGP.0b013e3180cc1fbcGoogle Scholar

C. Hocket al., “Decrease in parietal cerebral hemoglobin oxygenation during performance of a verbal fluency task in patients with Alzheimer’s disease monitored by means of near-infrared spectroscopy (NIRS)—correlation with simultaneous rCBF-PET measurements,” Brain Res. 755(2), 293–303 (1997).BRREAP1385-299Xhttp://dx.doi.org/10.1016/S0006-8993(97)00122-4Google Scholar

H. E. A. van Beek Arendaet al., “Cortical oxygen supply during postural hypotension is further decreased in Alzheimer’s disease, but unrelated to cholinesterase-inhibitor use,” J. Alzheimers Dis. 21(2), 519–526 (2010). JADIF91387-2877http://dx.doi.org/10.3233/JAD-2010-100288Google Scholar

R. CarandangD. W. Krieger, “Near infrared spectroscopy: finding utility in malignant hemispheric stroke,” Neurocrit. Care 6(3), 161–164 (2007).NCEACB1541-6933http://dx.doi.org/10.1007/s12028-007-0023-yGoogle Scholar

D. R. Hargroveset al., “Near-infrared spectroscopy in stroke: from research to clinical practice,” Stroke 35(11), 2430; author reply 2430–2431 (2004).SJCCA70039-2499http://dx.doi.org/10.1161/01.STR.0000144656.77330.18Google Scholar

C. W. Pennekampet al., “The value of near-infrared spectroscopy measured cerebral oximetry during carotid endarterectomy in perioperative stroke prevention: a review,” Eur. J. Vasc. Endovasc. Surg. 38(5), 539–545 (2009).1078-5884http://dx.doi.org/10.1016/j.ejvs.2009.07.008Google Scholar

P. RolfeY. A. WickramasingheM. Thorniley, “The potential of near infra-red spectroscopy for detection of fetal cerebral hypoxia,” Eur. J. Obstet. Gyn. Reprod. Biol. 42(Suppl. S), S24–S28 (1991).EOGRAL0301-2115Google Scholar

P. M. O’BrienP. M. DoyleP. Rolfe, “Near infrared spectroscopy in fetal monitoring,” Br. J. Hosp. Med. 49(7), 483–487 (1993).BJHMAB0007-1064Google Scholar

D. T. Delpyet al., “Cerebral monitoring in newborn infants by magnetic resonance and near infrared spectroscopy,” Scand. J. Clin. Lab. Invest. 47(Suppl. 188), 9–17 (1987).Google Scholar

M. Ferrariet al., “Cerebral blood volume and hemoglobin oxygen saturation monitoring in neonatal brain by near IR spectroscopy,” Adv. Exp. Med. Biol. 200, 203–211 (1986).AEMBAP0065-2598http://dx.doi.org/10.1007/978-1-4684-5188-7Google Scholar

M. CopeD. T. Delpy, “System for long-term measurement of cerebral blood and tissue oxygenation on newborn infants by near infra-red transillumination,” Med. Biol. Eng. Comput. 26(3), 289–294 (1988).MBECDY0140-0118http://dx.doi.org/10.1007/BF02447083Google Scholar

J. S. Wyattet al., “Quantification of cerebral oxygenation and haemodynamics in sick newborn infants by near infrared spectrophotometry,” Lancet 328(8515), 1063–1066 (1986).LANCAO0140-6736http://dx.doi.org/10.1016/S0140-6736(86)90467-8Google Scholar

G. LitscherG. Schwarz, Transcranial Cerebral Oximetry, Pabst-Science Publisher, Lengerich (1997).Google Scholar

J. Steinbrinket al., “Determining changes in NIR absorption using a layered model of the human head,” Phys. Med. Biol. 46(3), 879–896 (2001).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/46/3/320Google Scholar

Y. Hoshiet al., “Reevaluation of near-infrared light propagation in the adult human head: implications for functional near-infrared spectroscopy,” J. Biomed. Opt. 10(6), art. no.-064032 (2005).JBOPFO1083-3668http://dx.doi.org/10.1117/1.2142325Google Scholar

J. Selbet al., “Improved sensitivity to cerebral hemodynamics during brain activation with a time-gated optical system: analytical model and experimental validation,” J. Biomed. Opt. 10(1), 11013 (2005).JBOPFO1083-3668http://dx.doi.org/10.1117/1.1852553Google Scholar

V. Quaresimaet al., “Noninvasive measurement of cerebral hemoglobin oxygen saturation using two near infrared spectroscopy approaches,” J. Biomed. Opt. 5(2), 201–205 (2000).JBOPFO1083-3668http://dx.doi.org/10.1117/1.429987Google Scholar

B. W. Zeffet al., “Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography,” Proc. Natl. Acad. Sci. 104(29), 12169–12174 (2007).0370-0046http://dx.doi.org/10.1073/pnas.0611266104Google Scholar

J. Steinbrinket al., “Relevance of depth resolution for cerebral blood flow monitoring by near-infrared spectroscopic bolus tracking during cardiopulmonary bypass,” J. Thorac. Cardiovasc. Surg. 132(5), 1172–1178 (2006).0022-5223http://dx.doi.org/10.1016/j.jtcvs.2006.05.065Google Scholar

A. Liebertet al., “Time-resolved multidistance near-infrared spectroscopy of the adult head: intracerebral and extracerebral absorption changes from moments of distribution of times of flight of photons,” Appl. Opt. 43(15), 3037–3047 (2004).APOPAI0003-6935http://dx.doi.org/10.1364/AO.43.003037Google Scholar

D. Continiet al., “Multi-channel time-resolved system for functional near infrared spectroscopy,” Opt. Express 14(12), 5418–5432 (2006).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.14.005418Google Scholar

J. C. HebdenT. Austin, “Optical tomography of the neonatal brain,” Eur. Radiol. 17(11), 2926–2933 (2007).EURAE31432-1084http://dx.doi.org/10.1007/s00330-007-0659-1Google Scholar

J. SelbD. K. JosephD. A. Boas, “Time-gated optical system for depth-resolved functional brain imaging,” J. Biomed. Opt. 11(4), 044008 (2006).JBOPFO1083-3668http://dx.doi.org/10.1117/1.2337320Google Scholar

H. Wabnitzet al., “Time-resolved near-infrared spectroscopy and imaging of the adult human brain,” Adv. Exp. Med. Biol. 662 143–148 (2010).AEMBAP0065-2598http://dx.doi.org/10.1007/978-1-4419-1241-1Google Scholar

D. Continiet al., “Design and characterization of a two-wavelength multichannel time-resolved system for optical topography,” in Biomedical Optics 2006, Technical Digest, Optical Society of America, Washington, DC (2006).Google Scholar

M. Kacprzaket al., “Time-resolved optical imager for assessment of cerebral oxygenation,” J. Biomed. Opt. 12(3), 034019 (2007).JBOPFO1083-3668http://dx.doi.org/10.1117/1.2743964Google Scholar

J. Selbet al., “Discrimination between superficial and cerebral signals during functional brain imaging with a time-gated system,” in European Conferences on Biomedical Optics (ECBO), Munich, Germany (2005).Google Scholar

D. Continiet al., “Novel method for depth-resolved brain functional imaging by time-domain NIRS,” Proc. SPIE 6629, 662908 (2007).http://dx.doi.org/10.1117/12.728104Google Scholar

A. Liebertet al., “Bed-side assessment of cerebral perfusion in stroke patients based on optical monitoring of a dye bolus by time-resolved diffuse reflectance,” Neuroimage 24(2), 426–435 (2005).NEIMEF1053-8119http://dx.doi.org/10.1016/j.neuroimage.2004.08.046Google Scholar

A. Liebertet al., “Evaluation of optical properties of highly scattering media by moments of distributions of times of flight of photons,” Appl. Opt. 42(28), 5785–5792 (2003).APOPAI0003-6935http://dx.doi.org/10.1364/AO.42.005785Google Scholar

A. Liebertet al., “Intra- and extracerebral changes of hemoglobin concentrations by analysis of moments of distributions of times of flight of photons,” in Proc. SPIE 5138, 126–130 (2003).http://dx.doi.org/10.1117/12.500578Google Scholar

W. Becker, Advanced Time-Correlated Single Photon Counting Techniques, Springer, Berlin, Heidelberg (2005).Google Scholar

B. Chanceet al., “Phase measurement of light absorption and scatter in human tissue,” Rev. Sci. Instrum. 69(10), 3457–3481 (1998).RSINAK0034-6748http://dx.doi.org/10.1063/1.1149123Google Scholar

S. Carraresiet al., “Accuracy of a perturbation model to predict the effect of scattering and absorbing inhomogeneities on photon migration,” Appl. Opt. 40(25), 4622–4632 (2001).APOPAI0003-6935http://dx.doi.org/10.1364/AO.40.004622Google Scholar

H. Wabnitzet al., “A time-domain NIR brain imager applied in functional stimulation experiments,” Proc. SPIE 5859, 58590H (2005).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.632837Google Scholar

C. GolubC. van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore and London (1996).Google Scholar

A. Liebertet al., “Monte Carlo algorithm for efficient simulation of time-resolved fluorescence in layered turbid media,” Opt. Express 16(17), 13188–13202 (2008).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.16.013188Google Scholar

O. Steinkellneret al., “Optical bedside monitoring of cerebral perfusion: technological and methodological advances applied in a study on acute ischemic stroke,” J. Biomed. Opt. 15(6), 061708 (2010).JBOPFO1083-3668http://dx.doi.org/10.1117/1.3505009Google Scholar