1. Introduction

Near-infrared spectroscopy (NIRS) offers the ability to determine absolute absorption and scattering coefficients of biological tissues at multiple wavelengths. The retrieved absorptions allows quantification of different chromophores’ concentrations within the tissue, mainly oxy-hemoglobin (HbO₂) and deoxy-hemoglobin (Hb). Robust assessment of cerebral blood volume (CBV) and saturation (SO₂), derived from measured hemoglobin concentrations, is necessary for reliable studies of health and disease [1–5].

Different approaches of NIRS are continues-wave (CW) [6–8], frequency-domain (FD) [9–12] and time-resolved (TR) or time-domain (TD) [13–17]. Continuous-wave methods only measure the amplitude of light intensity. They need spatially or spectrally resolved information in order to disentangle the absorption and scattering contributions. Frequency-domain approach measures both amplitude and phase change of the input light. In the time-domain approach, the tissue is illuminated with a pulse of light with the duration of tens of picoseconds which will be detected at certain distance from the source after propagating a complicated path inside the tissue. The result is the histogram of time-of-flights of photons reaching the detector, known as photon distribution of time-of-flight (DTOF) or temporal point spread function (TPSF).

For reconstruction problem, a model must be considered for adult human head. For many years most common head model was a simple model described as a semi-infinite homogenous medium. This model had promising results for infants [1,4,10,12] but in case of adult human head, leads to strong contamination of optical properties of brain by extracerebral layers. For FD approach, Franceschini et al. [18] showed large errors due to a superficial layer in both phantom experiments and simulations on a two-layer slab geometry. Dehaes et al. [19] also found large errors (up to 45%), when investigating the FD method on adult head using simulated data on realistic head geometries. Using TD analytical solution of diffusion equation in homogenous medium, Kienle et al. [16] found strong contamination for data obtained from a two-layer phantom. Comelli et al. [20] fitted the simulated time-resolved data with a homogeneous model and found that estimated optical properties are much closer to superficial layers, similar to experimental results from Ohmae et al. [21] comparing TD-NIRS data with positron emission tomography (PET). To overcome the limitations of homogeneous model, two-layered models were developed. Pucci et al. [22] used broad-band CW approach in a two-layer model on phantom data. Kienle et al. [16] derived FD and TD analytical solution to diffusion equation in a two-layer geometry and studied its performance with Monte Carlo simulations and phantom data. Hallacoglu et al. [23] also, investigated FD solutions for two-layer geometry on phantom data and Monte Carlo simulations. Those studies showed improved results over homogeneous model.

Other complex models have been proposed to describe adult human head more accurately, including finite difference modeling of light propagation in true head anatomy obtained from a subject’s magnetic resonance imaging (MRI) scan [24,25], or a generic atlas model registered to a subject’s head [25,26] and Monte Carlo simulations on realistic head models [26–28]. While Monte Carlo method gives the most accurate description of the forward problem, its biggest drawback is being computationally intensive. In recent years several improvements have been introduced. Pifferi et al. [29] proposed a fitting routine based on a precomputed library of TPSFs for a number of reduced scattering coefficients in a homogenous medium. For any μ_s′ the TPSF could be interpolated from the library and then scaled for any absorption coefficient (μ_a) using Beer-Lambert law. Alerstam et al. [30] developed white Monte Carlo method that enable post-simulation scaling of scattering coefficient which also is limited to simple geometries as homogeneous semi-infinite or slab. Besides theoretical improvements, the technological developments in computer hardware led to production of powerful GPUs which can be used to increase the calculation capacity and reducing the simulation time significantly by a few orders of magnitude [31]. Fang et al. [32] introduced MCX software which is capable of utilizing several GPUs for Monte Carlo simulations in complicated geometries. Using precomputed TPSF libraries with GPU-based Monte Carlo simulations, Selb et al. [26] proposed a fitting routine for two different geometries, one based on a generic atlas model and the other a two-layered slab with variable thickness of superficial layer. While both models showed improved results over homogeneous model, their performances compared with each other were different for different subjects and at different locations on the head.

In summary, using a two-layer slab model, an atlas head model or a realistic model of subject’s head based on MRI scans will improve the results compared with homogenous model. However, these methods have their own drawbacks. Although a true model based on MRI scans can describe the forward problem most accurately, obtaining one is time consuming, costly and not applicable to all situations. A two-layer slab model can take into account the layered structure of the head and can sample different thicknesses of superficial layer but its flat geometry is still far away from true shape of the head with all its curvatures. On the other hand, an atlas registered to a subject’s head is capable of taking into account the shape of the subject’s head to some extent, however, it is unable to consider varying thickness of superficial layers in different subjects.

Here we propose our novel head model. Since the brain occupies the intracranial space, its shape roughly resembles the cranium and surface of the head. We can create a layered-head model solely based on the surface geometry of a subject’s head to estimate the anatomical structure of the brain and extra-cerebral layers without having a priori information on structural details inside the head. This model is potentially capable of taking into account both, the thickness of extra-cerebral layers and subject’s head shape, while former was absent in a generic atlas model and latter was not feasible in a slab model.

In this study we adopted the TD approach with the so-called null-distance source-detector (S-D) separation. TD is likely a candidate for future diffuse optical imaging devices because of its intrinsic advantages as well as theoretical and technological advancements in recent years [15,33]. The null-distance S-D method was proposed and tested experimentally which allowed source detector configurations with only few millimeters apart [34,35]. This opened the possibility for devices with high-density of optodes. It is also important since it increases the lateral resolution of the imaging system.

Therefore, the goals of the present study are: (1) to propose a new subject specific layered-head model to overcome the shortcomings of previous models and (2) to assess its performance on realistic TD-NIRS data obtained from a realistic head model and compare it with two-layered slab model with variable thickness presented by Selb et al. [26]. Realistic TD-NIRS data simulated using MCX software at several locations all over a subject’s head. We fitted this data with our proposed layered-head model, and compared the results to a two-layer slab model with variable thickness and characterized the performance of this method in terms of relative error.

2. Materials and methods

2.1 Monte Carlo simulations

In this study, we used Monte Carlo method to simulate 2 × 10⁸ photons using MCX software [32] for each source. The output of each simulation is a history file for each detector which can be used to calculate TPSF for each detector. This history file stores the partial path lengths (PPLs) of each photon. PPLs are the lengths that a photon travels inside different regions of the medium, for example different parts of a head like scalp, skull or brain. Having PPLs enables us to calculate new TPSFs for different absorption coefficients of each region by post-simulation scaling using Beer-Lambert law without the need for new simulations. However, changing the scattering coefficient requires a new simulation. The time of flight of each photon reaching the detector is calculated by,

To perform a light propagation simulation, first, it is necessary to define the geometry and optical properties of the medium. Then, the locations of the sources and the detectors must be defined along with their properties. Both will be discussed in details in the following sections.

2.2 Simulation of actual data set

2.2.1 True head model

A healthy adult head model (Subject #4) selected from Brainweb library [36] which is created based on anatomical MRI scans. It is segmented into ten different tissues. We combined seven tissues together to form a single superficial layer including skin, skull, muscle, fat, marrow, dura and connective tissue. Combining three remaining tissues, gray matter, white matter and vessels formed the brain region. Now we have a head model with two distinct regions, superficial and brain as depicted in Fig. 1 with light and dark regions respectively. This model has anatomical features similar to a real human head such as shape of the head, brain structure especially the complex cortex structure and different thicknesses of the superficial layer at different locations on the head.

 

Fig. 1 True head model consists of a superficial layer and brain region used for generating realistic TD-NIRS data. (a) locations of optodes on the head. (b) and (c) 45 selected optodes for this study. Blue circles show the location of sources. The number next to each circle is the ID number of that optode. Yellow circles show the location of nearest neighbouring detectors of each source. These detectors were used to record the TPSFs.

Download Full Size | PDF

2.2.2 Optode positioning

We considered a hemisphere large enough so it can be put on the head. 1666 points were placed on the hemisphere in an equidistance formation. Then, each point were projected in the radial direction onto the head as demonstrated in Fig. 1(a) similar to the method used by Strangman et al. [37]. New optode coordinates were calculated and stored. These coordinates will be used later for simulations. In practice, this hemisphere can be thought of as a helmet with mentioned points being the locations of the housings that hold the optodes – such as optical fibers. The final position of the fibers can be achieved with a mechanism that allows the optical fibers to slide through the housing and get fixed after touching the subject’s head. The new locations are no longer necessarily equidistant, so we found and stored the Euclidean distance between each optode and all the other optodes. The stored distances have two purposes; first, it can be used to define the neighborhoods of each optode and second, it will be used for placing the detectors at correct distance in layered-slab model which will be discussed later. Each optode can function as a source while all other optodes act as detector. Finally, after projection, some optodes may end up somewhere unsuitable like on the ears or close to the eyes and should be removed from study. In current study we selected 45 optodes – from 1666 optodes, and defined them as sources to sample different regions of the head as depicted Figs. 1(b) and 1(c). TPSFs recorded by surrounding detectors. Here we only used detectors in nearest neighborhood. In this configuration, most of the optodes have 6 surrounding detectors while a few of them have 5 and those on the perimeter have 4. Maximum and minimum S-D separations are 8.9 mm and 4.7 mm. The diameter of the detectors is set to 3 mm.

2.2.3 Optical properties

TPSFs are calculated for different combinations of superficial absorption coefficient, ${\mu }_{a,super}^{true}$, and brain absorption coefficient, ${\mu }_{a,brain}^{true}$. ${\mu }_{a,super}^{true}$ had 4 different values; 0.01, 0.015, 0.02 and 0.025 mm⁻¹. Values of ${\mu }_{a,brain}^{true}$ changed from 0.005 to 0.035 mm⁻¹ in 0.005 mm⁻¹ steps (7 values). It amounts to 4 × 7 = 28 different combinations. Scattering coefficients ${\mu }_s^{true}$ were set to 10 mm⁻¹ with anisotropy, g = 0.9. We considered the same ${\mu }_s^{true}$ for both layers since the scattering coefficient of the brain has small effect in results as discussed in [26]. The refractive index, n, set to 1.35 for both layers. These values covers the typical range of optical properties for different head tissues in NIR region from ~700-900 nm [38,39].

2.3 Reconstruction

True TPSFs were fitted using the TPSFs of corresponding S-D channels for both reconstruction models namely, layered-slab model (see Fig. 2(a)) and layered-head model (see Fig. 2(b)). We considered three different scattering coefficient for reconstruction; ${\mu }_s^{model}$ = 8, 10 and 12 mm⁻¹. As said before, each new ${\mu }_s^{model}$ requires a new simulation. In both models the thickness of the superficial layer can vary from 1 to 20 mm. However, we only changed the thickness from 10 to 20 mm assuming that there is at least a 10 mm superficial layer above the brain for an adult human. It reduced the number of fitting iterations and made calculations faster.

 

Fig. 2 Reconstruction models. (a) Layered-slab model; there are twenty sublayers on top. Large region below is considered as brain tissue. Blue circle is source and red circles show four of detectors in nearest neighborhood. (b) Layered-head model; the surface geometry is the same as true head model. The contour is the superficial layer divided into twenty sublayers. Large middle region considered as brain. In both models, the thickness of superficial layer can be changed from 1 to 20 mm.

Download Full Size | PDF

2.3.1 Reconstruction model 1: layered-slab model

We created our layered-slab model with similar procedure as described in [26]. First we considered a 10x10x8 cm (L × W × H) slab. Then, 2 cm of the slab from top divided into twenty sublayers each with a thickness of 1 mm as can be seen in Fig. 2(a). Thus, by manipulating the optical properties of sublayers we can change the thickness of the superficial layer from 1 to 20 mm. The source is positioned at the top center of the slab. We placed detectors around the source using the S-D Euclidean distances as discussed above. Figure. 2(a) shows four of such detectors.

2.3.2 Reconstruction model 2: layered-head model

To create the layered-head model, first, we found the boundary voxels of the true-head model. It is as if we have a priori information regarding the subject’s head surface. Using this surface, we created sublayers as we did with the slab. Unlike the slab geometry, it was not a simple surface and required a more complicated algorithm, but essentially we wanted to come up with the same thing which was a contour of sublayers with 1 mm thickness (see Fig. 2(b)). Since the surface of layered-head model has exactly the same shape as true model, then the locations of the sources are also the same as true model.

2.4 Fitting procedure

TPSFs were fitted by Levenberg-Marquardt method which finds the least square error between true TPSF and the ones from reconstruction models. TPSFs have 20ps-wide temporal bins and fitting range was from 1000 to 3000 picoseconds. Since the S-D separation is small and there is a huge burst of early photons, we applied this delay so the initial peak passes which is also common in practice to suppress the initial photons [34,35]. Moreover, late photons are more sensitive to deeper region of the tissue like brain and it has been shown that using a delay will improve the results [40–42]. In this time range the TPSF amplitude has a drop of about four orders of magnitude corresponding to 40dB drop in signal strength.

For each combination of ${\mu }_{a,superf}^{true}$ and ${\mu }_{a,brain}^{true}$, first, we fixed ${\mu }_s^{model}$ and change the thickness of superficial layer, d, from 10 to 20 mm. For each thickness, we fit for ${\mu }_{a,super}^{model}$ and ${\mu }_{a,brain}^{model}$ and stored the residuals. Then we changed ${\mu }_s^{model}$ and repeated this routine until we obtained ${\mu }_{a,super}^{model}$ and ${\mu }_{a,brain}^{model}$ and corresponding residuals for all the cases. Then we compared the residuals and the values corresponding to the smallest residual were reported as ${\mu }_{a,super}^{estim}$, ${\mu }_{a,brain}^{estim}$, ${\mu }_s^{estim}$ and $d^{estim}$. In this way we obtained the results for one S-D channel. Repeating this procedure for all desired S-D channels gives us the complete data set. Finally, we represents the ${\mu }_{a,brain}^{estim}$ at the location of each source by averaging over all detectors in the nearest neighborhood of that source.

3. Results

There are 45 sources with each source having 6 neighboring detectors. For each source-detector channel, there are 28 combinations (4 ${\mu }_{a,superf}^{true}$ × 7 ${\mu }_{a,brain}^{true}$) which amounts to ~7560 different cases. We performed the fitting procedure and obtained ${\mu }_{a,brain}^{estim}$ for all the cases. After averaging over 6 detectors the number of cases reduced to 1260.

Figures 3(a) and 3(b) show the estimated brain absorption, ${\mu }_{a,brain}^{estim}$, for layered-slab and layered-head models respectively at the location of source number 665. Being closer to the solid black line means more accurate estimation. For this case layered-slab model showed better performance. Figures 3(c) and 3(d) show the results for source number 1270 where layered-head model performed better. Different colors correspond to different ${\mu }_{a,superf}^{true}$. The observed underestimation is a common phenomenon which occurs when ${\mu }_{a,brain}^{true}>{\mu }_{a,superf}^{true}$.

 

Fig. 3 Estimated brain absorptions for (a) layered-slab model and (b) layered-head model at source number 665, an example where layered-slab had better performance. (c) Estimation of brain absorption for layered-slab model and (d) layered-head model at source number 1270 where layered-head model showed better performance. Each color corresponds to a different superficial layer absorption.

Download Full Size | PDF

Above examples suggest that the performance of reconstruction models depends on the location on the head. Figure 4(a) demonstrates the boxplots of relative error for individual sources for both layered-slab model (green boxes) and layered-head model (blue boxes). Relative error is calculated by,

 

Fig. 4 (a) Boxplot of relative error for 45 individual sources. Blue and green boxes correspond to layered-slab and layered-head models respectively. Black dots show the median values. Relative error at source locations on the head for (b) layered-slab model and (c) layered-head model.

Download Full Size | PDF

Thus far, we compared layered-slab and layered-head model at different locations on the head. To assess the overall performance, we consider relative errors of all cases, so each boxplot in Fig. 5(a) contains 1260 values. Relative error was 11.5% for layered-slab and 9.15% for layered-head. The layered-head performs better than layered-slab with reducing the error by 2.35 or 20%. The upper tails of the boxplots are equal and less than 46%. As we said before, when ${\mu }_{a,brain}^{true}>{\mu }_{a,superf}^{true}$ a considerable underestimation occurs which increases the relative error, so it is useful to find relative error for different values of ${\mu }_{a,superf}^{true}$ as demonstrated in Fig. 5(b) for both layered-slab model (green boxes) and layered-head model (blue boxes). Here each box contain 45(sources) × 7(${\mu }_{a,brain}^{true}$) = 315 values. It can be seen that relative error reduces as ${\mu }_{a,superf}^{true}$ increases, besides, layered-head models shows an improvement in all cases as presented in Table 1.

 

Fig. 5 (a) Total relative errors in ${\text{μ}}_{\text{a},\text{brain}}^{\text{estim}}$. (b) Relative errors for different ${\text{μ}}_{\text{a},\text{superf}}^{\text{true}}$. Blue and green boxes correspond to layered-slab and layered-head models respectively and red lines indicate the medians. Layered-head model has smaller medians in all cases.

Download Full Size | PDF

Table 1. Relative errors for different absorption coefficient of superficial layer

View Table

4. Discussion

4.1 Benefits of proposed source-detector configuration

4.1.1 Null-source detector separation

In CW method the depth sensitivity is related to the S-D separation [43], thus for probing deeper layers the distance between source and detector needed to be large in the order of few centimeters and the drawback is that the lateral resolution will become lower. In contrast, in time-domain approach, it is the time of flight of photons which corresponds to the depth sensitivity [43,44]. Because of that, detectors can be placed very close to the source. In our study the S-D separations are just a few millimeters. It enables us to produce a functional activity map with higher resolution.

4.1.2 Averaging over nearest neighbors

For most cases of ${\mu }_{a,brain}^{true}>{\mu }_{a,superf}^{true}$, majority of the detectors exhibit an underestimation while few detectors had more accurate estimation, thus averaging over all detectors may reduce the underestimation effect producing more accurate results. However, if a large error is present in one of the detectors, using the average could lead to large error in estimation. In the current study we used data from all detectors and did not removed any outliers, but doing so may improve the results further. We also compared the results using the median values instead of mean values for a few sources, but the results degraded so we chose the mean values.

4.2 Benefits of the layered-head model

We proposed our layered-head model to take into account shape of the head. At first glance, head surface may seem to be flat enough to approximate it with a layered slab for S-D separation of few millimeters which is the case in this study. However, to see the significance of the head geometry it is useful to explore the S-D sensitivity maps [45]. Figures 6(a) and 6(b) show the sensitivity maps of the superficial layer for true model and layered-head model respectively. Different rows show S-D channels at different locations on the head. Similarly, Figs. 6(c) and 6(d) show sensitivity maps inside the brain region for true and layered-head model respectively for the same S-D channels. Here we set the superficial layer thickness of the layered-head model similar to the true model. The superficial layer thickness defined as the radial distance between the source and surface of the brain. Regarding the superficial layer, it can be seen that layered-head model has a very similar sensitivity map with true model showing the effects of the curvatures. It was expected since this model was based on the surface geometry of the subject’s head. Inside the brain region the difference will become greater due to the complexities of the brain structure. In spite of that, layered-head model still resembles the shape of the brain to some extent. It is a feature that is absent in the flat geometry of the slab model. Having sensitivities more similar to true model for both superficial and brain regions is probably the reason that layered-head model had better performance than the slab model.

 

Fig. 6 Sensitivity maps of superficial layer for S-D channels at three different locations on the head for (a) true model and (b) layered-head model. Sensitivity maps of brain region for the same S-D channels for (c) true model and (d) layered-head model.

Download Full Size | PDF

4.3 Underestimation

4.3.1 Nature of the underestimation

To investigate the nature of the underestimation presents in our results, we generated a second set of data using our true head model and fitted the detector data using the first set. This is as if we know the exact structural details of the subject’s head for the inverse problem. The underestimation was still present in the results. Since reconstruction model was exactly the same as true model, the error and underestimation must be due to the noise. For ${\mu }_{a,superf}^{true}$ = 0.01, 0.015, 0.02, 0.025 mm⁻¹, the relative error in ${\mu }_{a,brain}^{estim}$ was 5.1%, 2.7%, 1.6% and 1.2% respectively. These results showed an increased error when ${\mu }_{a,brain}^{true}>{\mu }_{a,superf}^{true}$ similar to what we saw in our two reconstruction models. These values are the minimum reachable error using this fitting method and for the present noise level in this study.

4.3.2 Correction factor

The underestimation is a recurring pattern in the results and highly depends on the absorption coefficient of the superficial layer. Since fitting procedure, estimates both ${\mu }_{a,superf}^{estim}$ and ${\mu }_{a,brain}^{estim}$ simultaneously, it may be possible to find an underestimation correction factor based on the difference between ${\mu }_{a,superf}^{estim}$ and ${\mu }_{a,brain}^{estim}$ to improve the results. To do this, layered-head model should be used for both forward and inverse problem. This can be investigated further in future studies.

4.4 Fitting procedure

4.4.1 Time interval

In our fitting routine, we considered TPSFs in time interval 1000-3000 picoseconds. With this choice the initial burst of photons has passed. In this interval light attenuates by four orders which sets the experimental condition for recording the signal. Reducing the initial delay will lead to a greater attenuation by a few orders requiring a recording device with much higher dynamic range. To see the effect of time delay, we also fitted the data in time interval 0-3000 picoseconds. The overall relative error decreased from ~9% to ~8% for layered-head model. In spite this improvement, using the time delay could still be preferable because of technical issues in experimental measurements.

4.4.2 Nearest neighboring detectors

Here TPSFs from each detector fitted separately. Instead, it is possible to fit TPSFs of all detectors simultaneously, in another word finding the ${\mu }_{a,brain}^{estim}$ that minimizes the total residual for all TPSFs. This method can be investigated to see if it can improve the results.

Additionally, it is possible to include more detectors in the fitting procedure. But, by doing so we must use detectors with greater S-D separations which leads to lower lateral resolution because the volume that probed by the light will become larger. Thus, we only considered the first neighborhood. However, in case of using only the first neighborhood, S-D separations are almost the same for all detectors, so technically it is not a multi-distance time-domain (MDTD) approach. Farther detectors – such as second nearest neighborhood, can be used for more complex MDTD analysis.

4.5 Limitations

4.5.1 Surface geometry of the subject’s head

The layered-head model we introduced in this paper does not require a priori information on subject’s head internal structure. However, it still relies on the knowledge of shape of the head’s surface and it must be obtained by some other means. In contrast, layered-slab model does not require any of that which is an advantage.

4.5.2 True head geometry

For this study we used an anatomical MRI scan and divided it into two regions namely superficial layer and brain by combining different tissues. We used this simpler model to evaluate the proposed null-distance S-D configuration and also investigate the improvement from layered-head model compared with layered-slab model. This model is still an approximation and does not include all the complexities of a human’s head anatomy, for example transparent CSF layer, vessels and hair follicles [46]. Moreover, the systemic absorption change in scalp may be greater than skull. These differences are very likely to affect the results and needs more investigation. One suggestion is to use a model with more than two layers although it is computationally more intensive and time consuming.

5. Conclusion

In this paper we proposed a null-distance S-D configuration for TD-NIRS measurements and implemented it on two reconstruction models; first, a layered-slab model and second a novel layered-head model which is based on the surface shape of the subject’s head. Our results showed total relative errors 11.5% for layered-slab and 9.15% for layered-head model. We also found that at 32 out of 45 (70 percent) of the locations layered-head model had less error than the slab model with a maximum error of 24% for layered-slab model and 15% for layered-head model. Additionally, we found that the performance highly depends on the absorption coefficient of superficial layer; however, layered-head model showed improved results for all superficial absorptions. Moreover, we observed an underestimation behavior which can be studied more extensively to find a possible compensation method to improve the results. Finally, our study was based on a null-distance configuration which is important since it allows for placing optodes with higher density and also increases the lateral resolution of obtained absorption maps.