1.

Introduction

Various continuous-wave (CW),^1–5 frequency-domain (FD),^6–14 and time-domain (TD)^15–19 near-infrared spectroscopy (NIRS) approaches offer the ability to determine the absolute absorption and scattering coefficients of biological tissue. The retrieved optical absorption measured at multiple wavelengths allows quantification of different chromophores’ concentrations within the tissue. For brain imaging, the robust assessment of cerebral blood volume (CBV) and oxygenation, derived from the measure of hemoglobin concentrations in the brain, is essential for reliable cross-sectional and longitudinal studies of health, disease, and disease progression.^8,9,19–22

Continuous-wave methods, such as broadband or hyperspectral approaches originally proposed more than 15 years ago,² require spatially¹ or spectrally^2–5 resolved information in order to disentangle the contributions from tissue absorption and scattering. The frequency-domain multidistance (FDMD) approach based on a homogeneous model⁶ has been extensively validated with Monte Carlo simulations,¹⁴ phantoms,^6,12 and animal models.^11,13 It has been successfully applied to the monitoring of brain oxygenation and metabolism in healthy and brain-injured infants.^8,9,20 Time-resolved approaches are generally based on the nonlinear fit of temporal point spread functions (TPSFs). They have been validated with Monte Carlo simulations and phantoms,^15,23,24 and have been applied to monitor developmental cerebral changes in infants.²² In adults, the TD-NIRS technology has been applied to monitor variations in brain oxygenation and blood volume during cardiopulmonary bypass surgery²⁵ or to detect vasospasms following subarachnoid hemorrhage.²¹

Most of the in vivo studies mentioned above have relied on a simple head model described as a homogeneous semi-infinite medium.^{1,6–10,16,19–22,25,26} This model has shown promising results in piglets and infants,^{5,8,9,11,20,22,27} but it is widely recognized that, in the case of the adult head, its oversimplification causes strong contamination of the brain optical properties by those of the extracerebral tissue. For the FDMD approach, Franceschini et al.¹² have shown, with simulations and phantom measurements in a slab geometry, that when a superficial layer thicker than $\sim 1\mathrm{cm}$ is present, the error on the retrieved absorption of the second layer can exceed 50%. We have previously investigated with simulated data the performance of the FDMD method on realistic head geometries at different ages.²⁸ We showed that, while it provides accurate results in infants up to 1 year of age (10 to 15% error), its application to adult heads introduces large errors (20 to 45%). The FDMD method is therefore not directly translatable to adult head measurements. For a simple two-layer phantom geometry, Kienle et al. showed that using the TD analytical solution of the diffusion equation for a homogeneous medium induces strong contamination of the second layer optical properties by those of the first layer.¹⁷ Similarly, based on Monte Carlo simulations guided by in vivo measurements on the adult forehead, Comelli et al. showed that time-resolved data fitted with a homogeneous model return absorption and reduced scattering coefficients much closer to superficial layer values (scalp and skull) than to those of deeper layers (white and gray matter).¹⁶ This was further confirmed experimentally by Ohmae et al., who compared baseline absolute measurement of CBV obtained by positron emission tomography (PET) and by TD-NIRS.²⁹ While the two modalities showed good correlation, the PET measures of CBV were 50% higher than those assessed by TD-NIRS, which in turn were more similar to the PET measure of scalp blood volume.

These results highlight the importance of developing more realistic models of light propagation in the adult head. For this purpose, the implementation of layered models instead of the homogeneous geometry is becoming more and more common in combination with the different NIRS approaches. Pucci et al. combined the CW broadband approach with a two-layer model in phantom data and retrieved the chromophore concentrations in the second layer within a 10% error.⁴ Kienle et al. derived an analytical solution to the diffusion equation in a two-layer geometry for TD and FD.¹⁷ This TD approach has shown improved results over the homogeneous model, as demonstrated by Monte Carlo simulations on two-layer slabs, and layered phantom experiments.¹⁷ We have applied it to in vivo data on the adult head, where it reported brain optical properties distinct from that of the extracerebral layer.³⁰ Paralleling the TD approach, Hallacoglu et al. demonstrated the improvement afforded by the FD two-layer geometry¹⁷ with Monte Carlo and phantom data, and applied the technique to in vivo data obtained on the adult head.³¹ However, in the few in vivo studies using a layered model of the adult head,^30,31 there was no validation of the retrieved brain values by complementary modalities. More generally, the performances of the two-layer analytical methods have only been assessed through simulations or phantom measurements on simple two-layer slab geometries and not on realistic head structures. Since the two-layer slab model remains a crude approximation of the complex head structure, it could lead to significant errors when employed on real human measurements. The accuracy of the layer geometry on realistic human data still remains to be evaluated.

Finally, a few more complex models of the adult head have been proposed, including analytical solutions of the diffusion equation or of the radiative transfer equation for multiple layer slabs,^32,33 finite difference modeling of light propagation in the true head anatomy of a subject obtained from a magnetic resonance imaging (MRI) scan,³⁴ and Monte Carlo approaches^35,36 in different geometries. Monte Carlo methods provide the most accurate description of the forward problem, but because of long computation times, they have rarely been implemented in nonlinear fitting routines that require a new simulation for each iteration of the fitting process. Pifferi et al. proposed a fitting routine based on Monte Carlo simulations in a homogeneous medium.³⁵ They used a library of TPSFs precomputed for a number of scattering coefficient values, which could then be interpolated for any ${{\mu }_s}^{\prime }$ value and scaled using the Beer-Lambert law for any ${\mu }_a$ value. The method provided more accurate fitting results than the analytical solution of the diffusion equation but was limited to homogeneous semi-infinite or slab media. Truly ${{\mu }_s}^{\prime }$-scalable or white Monte Carlo methods have been developed that enable postsimulation scaling of the scattering coefficient, but they are currently limited to very simple geometries, such as a homogeneous semi-infinite or slab medium.³⁶ While Monte Carlo approaches theoretically offer flexibility in terms of head geometry, they were, in practice, limited to simple homogenous geometries for which a library of Monte Carlo results could be computed ahead of time, in order to avoid the long computation time required for Monte Carlo simulations. However, with the advent of fast Monte Carlo approaches based on graphical processing unit (GPU) computation,^37,38 the use of numerical approaches to model the forward problem of light propagation has become a viable alternative to analytical solutions.

In summary, the CW-, FD-, or TD-NIRS estimations of the adult brain optical properties based on a homogeneous model of the head are known to introduce significant contamination from extracerebral layers. More complex models have been proposed and occasionally applied to in vivo data, but their performances were never characterized on realistic head structures. Furthermore, the early implementation of Monte Carlo approaches has been restricted to simple geometries and, therefore, has not taken full advantage of the flexibility this modeling offers in terms of geometry.

Therefore, the goals of the present study are twofold: (1) to implement a flexible Monte Carlo-based fitting routine of TD-NIRS data to retrieve the brain optical properties and (2) to characterize its performance on realistic time-resolved adult head data and compare them to the homogeneous analytical solution of the diffusion equation. Multidistance TD data were generated with a GPU-based Monte Carlo code at different locations over the whole head on three subjects, including the frontal, parietotemporal, and occipital cortices. The simulated data were then fitted with different models of light propagation and head geometries. Specifically, we compared the homogeneous analytical solution of the diffusion equation and Monte Carlo simulations on two types of head geometries: a two-layer slab where the thickness of the first layer can be fixed or estimated and a generic atlas head registered to the subject’s scalp using superficial landmarks. We characterized the performance of each fitting approach in terms of relative error on the retrieved brain absorption, linearity, and cross-talk from extracerebral tissue. The present work is a continuation of our preliminary study,³⁹ extended here to more head locations and more subjects, and with a thorough characterization of the different fitting approaches through new metrics.

2.

Methods

2.1.

Simulated Data

Realistic TD-NIRS datasets were generated using GPU-based Monte Carlo simulations on adult heads whose structures were obtained from MRI anatomical scans. Three subjects were selected from a library of 32 head volumes, previously used by Cooper et al. to validate the use of atlas-based image reconstructions of functional data.⁴⁰ We chose three subjects which resulted in poor, average, and good performance in the functional NIRS study. The head volumes were segmented using FreeSurfer ( http://surfer.nmr.mgh.harvard.edu)^41,42 into four tissue types: extracerebral tissue (skin and skull), cerebrospinal fluid (CSF), gray matter, and white matter. CSF, gray, and white matters were subsequently combined into a single brain tissue type. The details of the MRI scan preprocessing steps and of the segmentation procedure have been previously described.⁴⁰

Contrary to a previous study where we reported the performance of the FDMD approach at one specific location on the head,²⁸ in the present work, we show a systematic characterization of the TD-NIRS fitting methods at various locations over the whole head. We defined three two-dimensional (2-D) arrays of optodes with 1 cm spacing, to be placed over the frontal ($6\times 20$ optodes), left parietotemporal ($12\times 14$), and occipital ($8\times 20$) regions. For each array, four additional dummy optodes were defined to anchor the probe on a specific location on the surface of the head. The probe was then wrapped onto the scalp using an iterative, spring-relaxation algorithm that has been previously described by Cooper et al.⁴⁰ The algorithm returns the three-dimensional (3-D) coordinates of all optodes on the surface of the head. The resulting probes wrapped on the head surface of subject 1 are displayed in Fig. 1. Note that in a real-life experiment, these coordinates could be obtained from 3-D digitization of the optode locations with respect to superficial landmarks on the scalp, such as 10–20 reference points.⁴³

Fig. 1

Three-dimensional view of the segmented head of subject 1, with the virtual arrays of optodes wrapped onto the surface of the head, over the frontal, left, and occipital regions. The middle figure also presents an example probe geometry at one optode location, consisting of one source and four detectors in a row at $\sim 10$, 20, 30, and 40 mm.

The true optical properties of the head were defined as follows. The reduced scattering coefficient was set to $10{\mathrm{cm}}^{-1}$ in all tissue types. The CSF, gray, and white matters were combined into a single brain tissue type. The brain absorption coefficient was varied between 0.05 and $0.3{\mathrm{cm}}^{-1}$ in steps of $0.05{\mathrm{cm}}^{-1}$, while the absorption in the extracerebral layer (combined scalp and skull) was set to 0.08, 0.1, or $0.2{\mathrm{cm}}^{-1}$. These absorption properties were chosen to cover the broad range of values reported in the literature over the 690- to 850-nm spectral range.^{8,16,22,30,44} For each subject, and all three arrays of optodes, we simulated time-resolved NIRS data with the GPU-based code Monte Carlo eXtreme (MCX) developed by Fang et al.³⁷ Each optode of a probe was set successively as a source, while all other optodes from the same probe acted as detectors with 1 mm radius. For each source, we launched ${10}^9$ photons, which took $<5\min$ on an NVIDIA C2050 Tesla GPU. A total of 448 MCX simulations (one per optode location) were run in parallel on a 16 GPU node cluster, for a total duration of $<4\mathrm{h}$ per subject. For each source, the MCX simulation returns a history file, which contains a list of all detected photons and their individual partial path lengths in each tissue type.⁴⁵ The resulting detected time-resolved fluence for any absorption coefficient combination in the different tissue types can be computed from this history file by applying the Beer-Lambert law, without launching a new simulation. Note, however, that a set of different scattering properties would require a different simulation.

To simulate a typical multidistance measurement setup, we considered the TPSFs from one source and a line of four detectors at 1, 2, 3, and 4 cm (see example on the middle panel of Fig. 1). The resulting TPSFs were computed by integrating the detected photons received over 50-ps-wide temporal bins. Note that for each source, we recorded the photons detected on all other optodes, but only four detectors were considered in the fitting process. It is possible that different geometries would be more advantageous and improve the performance of the different fitting procedures, but we did not investigate this aspect. For ${10}^9$ photons launched at each source, the resulting total numbers of detected photons were approximately $8\times {10}^5$, $1\times {10}^5$, $3\times {10}^4$, and $1\times {10}^4$ at 1, 2, 3, and 4 cm separations (average over all subjects and all locations), with a corresponding noise level at the peaks of the TPSF of 0.3%, 1.5%, 4%, and 7%, respectively. We used the MCX-generated TPSFs as the virtual data without additional noise or convolution by an instrumental response function in order to estimate the performance of the different fitting procedures in a best-case scenario.

3.

Results

Figure 2 shows the fitting results for one source location on the left probe of subject 2, for the extracerebral absorption set to $0.1{\mathrm{cm}}^{-1}$ and brain absorption varying between 0.05 and $0.3{\mathrm{cm}}^{-1}$ by steps of $0.05{\mathrm{cm}}^{-1}$. The homogeneous model yields underestimated values for brain absorption, consistent with high contamination by the lower extracerebral absorption. The Monte Carlo fit based on the atlas geometry returns better results. The absorption of the brain is slightly overestimated, but shows good linearity with respect to the true values. Note that this is not always the case and that at different locations the atlas fit can retrieve under- or overestimated brain absorption. In the case of the two-layer slab fit, the accuracy of the retrieved brain absorption depends on the assumed thickness of the first layer. The 6-mm layer results in highly underestimated absorption (similar to the homogenous case), and at the other extreme, the 14-mm layer yields highly overestimated brain absorption. Interestingly, selecting the thickness that results in the smallest residual for each fit (diamonds) yields good results, suggesting that fitting for the extracerebral layer thickness is possible.

Fig. 2

Example of fitting results for subject 2 at one probe location on the left temporal region. The absorption coefficient in the extracerebral layer was set at $0.1{\mathrm{cm}}^{-1}$, while the brain absorption varied between 0.05 and $0.3{\mathrm{cm}}^{-1}$. For the Monte Carlo based on the layered slab geometry, the results of the fit for each thickness of the extracerebral layer are presented. The diamonds show the retrieved absorption for the thickness corresponding to the minimal residual.

Figure 3 presents, for all three regions of one subject, the 2-D maps of the median error on the retrieved extracerebral and brain absorptions when applying the homogeneous model. The magnitude of the error varies smoothly over the whole brain, with regions of higher brain inaccuracy corresponding to regions of better accuracy on the extracerebral absorption. This is most probably due to the regionally varying thickness of extracerebral tissue resulting in varying contamination.⁴⁸ Notice, for instance, how the error on the brain absorption is maximal (median 45%) on the lower part of the frontal probe, where the brain is further away from the scalp.

Fig. 3

Two-dimensional maps of the median error on the retrieved extracerebral and brain absorptions, when applying the analytical solution of the diffusion equation for a homogeneous semi-infinite medium on the simulated data for subject 1. Each pixel represents the median error computed at the corresponding optode location over all combinations of ${\mu }_{a,\text{Extra-cerebral}}$ of 0.1 or 0.2 and ${\mu }_{a,\text{Brain}}$ varying between 0.05 and $0.3{\mathrm{cm}}^{-1}$ by steps of $0.05{\mathrm{cm}}^{-1}$.

The overall performance of the different fitting approaches is summarized in the box-and-whiskers plots of Fig. 4. The thick horizontal line in each box represents the median, and the box limits present the 25- and 75-percentiles of each metric over all locations and all three subjects. The whiskers extend from the 10- to the 90-percentiles of the data. The smaller box plots on the right of each plot show the performance for the individual subjects. In the case of the Monte Carlo slab model, we show the results for the thickness corresponding to the lowest fit residual, i.e., when we estimate the extracerebral layer thickness based on the goodness of individual fits (as represented by diamonds in Fig. 2). The analytical homogeneous model yields the worst performance as characterized by all metrics: the median error on ${\mu }_{a,\text{Brain}}$ is $\sim 20\%$; linearity is poor with a median slope of 0.4 and low $R^2$ of 0.92; and cross-talk is $\sim 8\%$ for a 20% variation in the extracerebral absorption.

Fig. 4

Box-and-whisker plot of all performance metrics computed over all locations and all subjects: relative error on brain absorption; cross-talk from extracerebral layer; and linearity (slope and $R^2$) of brain absorption, as detailed in Sec. 2.3. For the two-layer slab geometry, we chose the thickness corresponding to the minimal residual for each fit. For each metric, the smaller graphs on the right show the individual subject results.

The Monte Carlo approach on layered models improves all metrics, both for the two-layer slab and for the atlas head geometries. Accuracy is only moderately improved, from $\sim 20\%$ error with the homogeneous model down to $\sim 14$ and 15% for the slab and the atlas model, respectively. However, the linearity is more dramatically improved with the layered geometries, both in terms of slope (from 0.41 to 0.77 for the slab and 0.87 for the atlas) and $R^2$ of the linear fit (from 0.92 to $\sim 0.97$ for both the slab and the atlas). Finally, the cross-talk from extracerebral layer is decreased from 7.8% with the homogeneous model down to 5.5 and 4.7% for the slab and atlas geometries, respectively.

The performance of the Monte Carlo approach on layered models for both geometries (two-layer slab and atlas) varies with subject (as can be seen on the smaller box plots of Fig. 4), and within each subject with location (as illustrated by large whiskers for each subject). For instance, in subject 1, the use of the atlas is very advantageous (relative error and cross-talk almost divided by two), while for subject 3, the homogeneous model results in good accuracy (14% error on average), which is not improved by any of the Monte Carlo fitting. Note, however, that linearity and cross-talk are improved in all cases. While the linearity of the brain absorption is strongly increased with the Monte Carlo methods (slope increasing from 0.4 to 0.8 and 0.9), it presents a relatively large interquartile range from 0.6 to 1.4, reflecting large variations between locations. Similarly, we speculate that these variations are due to the local head geometry, and more specifically to the thickness of the extracerebral layer.

While we note relatively high variability between subjects, as well as spatial variability over the head within each subject, we did not observe systematic differences between the three probe locations we studied (frontal, left, occipital, individual results not shown). The forehead is a region of particular interest for the characterization of baseline brain optical properties because of the underlying prefrontal cortex involved in diverse cognitive processes. Furthermore, it is unencumbered by hair, which facilitates measurements. We studied the results of the different fitting approaches over a region of interest on the forehead, limited to the four upper rows of the frontal optodes and leaving out the three outer columns on each side of the probe (where more curvature is likely to deteriorate the results). The performances of all fitting approaches over this specific region were slightly better than average, with medians of the relative error of 13.6% (homogeneous), 11.8% (slab), and 14.2% (atlas). The corresponding values for cross-talk were 5.5% (homogeneous), 4.0% (slab), and 3.3% (atlas). The linearity was characterized by slopes of 0.55 (homogeneous), 0.85 (slab), and 1.40 (atlas) and $R^2$ of fit of 0.95 (homogeneous), 0.98 (slab), and 0.97 (atlas).

4.

Discussion

4.2.

Assumption of Homogeneous Scattering

For this study, we used a homogeneous scattering over the whole medium in the reconstruction process (slab or atlas head). We justify this assumption by the low sensitivity of TD-NIRS measurements to the scattering coefficient of the bottom layer of a two-layer medium for thicknesses typical of the adult brain (6 to 14 mm).

As a simple illustration, Fig. 5 shows the evolution of the simulated TPSFs (Monte Carlo) for a source-detector separation of 2 cm on a two-layer medium with a 1-cm-thick first layer, for ${\mu }_{a1}$ [Fig. 5(a)] or ${\mu }_{a2}$ [Fig. 5(d)] varying between 0.05 and $0.40{\mathrm{cm}}^{-1}$, and ${{\mu }_{s1}}^{\prime }$ [Figs. 5(b) and 5(c)] or ${{\mu }_{s2}}^{\prime }$ [Figs. 5(e) and 5(f)] varying between 4 and $20{\mathrm{cm}}^{-1}$. The effect of ${\mu }_{a1}$ can be seen from early time delays [Fig. 5(a)], while ${\mu }_{a2}$ affects the slope of the TPSF at later delays [Fig. 5(d)]. Increasing ${{\mu }_{s1}}^{\prime }$ shifts and broadens the TPSF [Figs. 5(b) and 5(c)], whereas the effect of ${{\mu }_{s2}}^{\prime }$ is very small for a low absorption of $0.1{\mathrm{cm}}^{-1}$ [Fig. 5(e)] and almost unnoticeable over the displayed five orders of magnitude for a higher absorption of $0.2{\mathrm{cm}}^{-1}$ [Fig. 5(f)]. Furthermore, note that the range of scattering values we studied (4 to $20{\mathrm{cm}}^{-1}$) is broader than that typically reported in the literature. Similarly, we (data not shown) and others^31,49 have observed that while TD or FD two-layer models enable reasonable recovery of ${\mu }_{a1}$, ${\mu }_{a2}$, and ${{\mu }_{s1}}^{\prime }$ on layered media, in contrast, ${{\mu }_{s2}}^{\prime }$ can only be estimated with high uncertainty due to the low sensitivity of the TPSF to this parameter.

Fig. 5

Monte Carlo simulated temporal point spread functions (TPSFs) for a 20 mm source-detector separation on a slab with a 10-mm-thick first layer. The absorption and reduced scattering coefficients of each layer are varied independently to visualize their effect on the resulting TPSF. Upper row shows the effect of the optical properties of the superficial layer: (a) varying absorption coefficient, (b) varying scattering coefficient for ${\mu }_a=0.1{\mathrm{cm}}^{-1}$, (c) varying scattering coefficient for ${\mu }_a=0.2{\mathrm{cm}}^{-1}$. Lower row shows the effect of the optical properties of the deep layer: (d) varying absorption coefficient, (e) varying scattering coefficient for ${\mu }_a=0.1{\mathrm{cm}}^{-1}$, (f) varying scattering coefficient for ${\mu }_a=0.2{\mathrm{cm}}^{-1}$. The default values are ${\mu }_a=0.1{\mathrm{cm}}^{-1}$ [(a), (b), (d), and (e)] or ${\mu }_a=0.2{\mathrm{cm}}^{-1}$ [(c) and (f)], and ${{\mu }_s}^{\prime }=10{\mathrm{cm}}^{-1}$. Absorption ${\mu }_a$ is varied from 0.05 to $0.40{\mathrm{cm}}^{-1}$ in steps of $0.05{\mathrm{cm}}^{-1}$, while reduced scattering ${{\mu }_s}^{\prime }$ is varied from 4 to $20{\mathrm{cm}}^{-1}$ in steps of $2{\mathrm{cm}}^{-1}$.

Using a homogeneous scattering coefficient over the whole medium allows faster simulations and fitting. Recall that each set of scattering coefficients requires a new simulation. Scalable Monte Carlo methods have been proposed, where the effect of a change in scattering on the detected fluence can be computed with a single Monte Carlo simulation.^50,51 This approach is, however, limited to small changes in scattering and, in practice, is implemented to model sensitivity to scattering changes but not its absolute value. Fully scalable Monte Carlo, also known as white Monte Carlo, uses postprocessing of temporal and spatial binning, but they currently require a homogeneous infinite or semi-infinite model.³⁶ Instead, we used a library approach as was previously employed by different groups.³⁵ Sampling ${{\mu }_s}^{\prime }$’ from 6 to $16{\mathrm{cm}}^{-1}$ by steps of $2{\mathrm{cm}}^{-1}$ required six simulations. Instead, if scattering is different between the two layers, the numbers of simulations would increase to 36 for a fixed layer thickness, or 180 combinations for the thickness of the first layer varying between 6 and 14 mm in 2-mm increments.

5.

Conclusion

We implemented a Monte Carlo-based baseline optical property fit of TD-NIRS data using either a two-layer slab geometry with an extracerebral layer of varying thickness or an atlas head registered to the subject’s surface. We generated a library of Monte Carlo data for the slab geometry, which can be used for any subject. For the atlas approach, new Monte Carlo sets need to be created for each subject, which can be done in $<1\mathrm{h}$ per source with GPU-based computation. The fitting procedure itself takes only a few minutes. The new approach was tested on extensive datasets of simulated measurements realistic of the adult human head, as opposed to previous studies relying on simulations or phantom data in layered slabs. We found that both Monte Carlo–based approaches offered improved performance compared to the homogeneous solution in terms of accuracy (25% error reduction), linearity (slope of 0.85 instead of 0.4), and cross-talk (40% reduction). The best option (slab or atlas) seems to be subject-dependent, suggesting the possibility of further improvement based on a combined approach.