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 muscle34.–5 and brain studies.67.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,1011.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.1818.104.22.168.–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.2526.–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 experiments30 or the analysis of the broadening of short laser pulses during their propagation in the tissue under investigation.25,27,3132.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.3536.–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 photons27,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 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 counting43 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 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 varying the number of measured quantities and the source-detector separations.
The equations that will be derived in this section are related to the general case of a medium with compartments. The specific results shown in Sec. 4 refer to a two-layered medium, i.e., a superficial layer with thickness 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 whereas the scattering properties remain constant ().
The time-of-flight distributions are obtained as histograms of photon counts in time channels . We start from the following statistical moments of measured DTOFs, the 0th moment or total photon count
The following depth-resolved analysis is based on the changes in attenuation , mean time of flight and variance 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 in the individual compartments of the model and changes of moments of DTOFs is linear,Ref. 31. The vectors , and contain changes of attenuation, mean time of flight and variance of the DTOFs at source-detector separations , i.e., , , . The vector contains changes of in different compartments of the model, i.e., .
The sensitivity matrices MPP (mean partial pathlength), MTSF (mean time of flight) and VSF (variance of the DTOF) are 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 simulations31 or by applying a perturbation method based on the diffusion equation.37,45,46 For a measurement carried out at source-detector separations the system of equations given by Eqs. (5)–(7) consists of equations with unknowns, i.e., the changes of the absorption coefficient in the compartments of the model. Let , and denote the estimates for , and , i.e., the measured data. An estimate may then be determined by the weighted least-squares approach according to4710) is an covariance matrix of the corresponding measurements.
The covariance matrix of is then obtained as11) are the variances. The expression represents the standard deviation of the estimated absorption changes for the j-th compartment. Note that care needs to be taken in the calculations when is nearly rank-deficient.47
In the following we give expressions for the covariance matrix (for details see Appendix), assuming that the (dominant) photon noise follows a Poisson distribution.12)–(14) and (17) indicate that the respective moments are statistically dependent in case they are derived from the same DTOF, i.e., from the same distance . The covariances between and as well as between and turn out to be zero even for [Eqs. (15) and (16)]. However, they do not vanish in general as a property of the moments, but rather because of the specific property of the Poisson distribution (see Appendix).
For the following application we focused on the particular case of 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 as well as the covariance matrix depend on the background optical properties. We presumed a semi-infinite homogeneous turbid medium with absorption coefficient , reduced scattering coefficient and refractive index which are values typically taken for simple head models.25,31,48
The sensitivity matrix [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 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 of the upper layer. Monte Carlo calculations were carried out for a total of 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 matrix [see Eq. (10)] requires the knowledge of the moments , and as well as the total photon count [see Eqs. (12)–(17)]. For each source-detector separation the higher moments were derived from the DTOFs as obtained from the Monte-Carlo simulations described above. Regarding , we introduced experimentally realistic conditions. We assumed the total number of detected photons to be , 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 depends only on the covariance matrix 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 of the data is available (or can be calculated) as well as the sensitivity matrix 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.
This section presents the results for the standard deviations of the absorption changes in the two layers of the turbid medium. We start from a dataset , or related to , or , 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 and 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 were repeated with sensitivity matrices corresponding to various values of the thickness of the upper layer.
Figure 2 displays the results obtained for various combinations of moments when all four inter-optode distances were included. The calculated standard deviations and of the changes of the absorption coefficient of the lower and upper layer, respectively, are shown as a function of thickness 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, increases and decreases for increasing thickness . 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 [, or , see Fig. 2(a) and 2(c)], is definitely most suitable when analyzing the upper layer [Fig. 2(c)], whereas and become similarly advantageous when determining , in particular for larger . This finding is in line with the different depth dependence of the sensitivity factors of the various moments31 that also changes with the inter-optode distance .
As can be seen in Fig. 2(b) and 2(d), in all cases is reduced when including higher moments. In particular, we compared the standard deviations obtained for the following three cases, (i) only, measured at four source-detector separations, (ii) combined with , (iii) the full set , and at four source-detector separations is taken into account. As an example, for , a realistic thickness of the extracerebral layer in an adult, the inclusion of improves the precision of by a factor of 2.5. Additional inclusion of 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 . As a start, a single detector at a large distance () was chosen, and shorter source-detector separations were incorporated step by step. The data set containing changes in two moments of DTOFs, i.e., and 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 . The influence of the inclusion in the analysis of data acquired at shorter is particularly small for the deeper layer in case of large [see Fig. 3(a)]. This finding is obvious, since detection of moments at small distances is rather insensitive to deep absorption changes. Analysis of the analogous data set containing all three moments of DTOFs, i.e., , and , showed a similar pattern of changes in as a function of for the various combinations of source-detector separations.
In order to study the potential of a single-distance time-resolved measurement, the analysis of 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 ( and ) results in high , 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 changes31 which is a consequence of their normalization to by definition [see Eqs. (2) and (3)]. Comparison of the results for the combination , with the results for all moments , , shows that both options lead to almost identical at large thickness [see Fig. 4(a) and 4(b)]. In contrast, at small the values of converge for the combinations , and , , .
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 increases by 1 to 2 orders of magnitude with decreasing source-detector separation. Moreover, the range of between small (2 mm) and large (18 mm) thickness 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 . 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.
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 shown above were calculated assuming the same number () of photons acquired at all source-detector separations. Since all elements of the covariance matrix scale with , the results can be simply rescaled to other common values of . 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 and in most cases. With an assumed background absorption coefficient of , a signal-to-noise ratio of 1 would be attained with a relative absorption change between 0.1% and 1%.
Most of the results were presented as a function of the thickness of the upper layer that was assumed to be known. This is realistic if anatomical images are available. In principle, could also be estimated from the measurements. However, the dependence of the measurands on 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 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 and 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 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 which is the major measurand in brain imaging. However, a reconstruction based on and 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 in addition to and 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.
Here the derivation of the elements of the covariance matrix is described in more detail. We consider the sub-matrices of where and stand for the vectors , and which contain the estimates of changes in attenuation, mean time of flight and variance measured at source-detector separations. Since these measurements are statistically independent it follows that each sub-matrix is diagonal, i.e.,1)–(4)]. The total attenuation is where the number of input photons is assumed to be constant. 12)–(17)] can be derived.
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.”