## 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.”