*in vivo*data, but their performances have never been characterized on realistic head structures. Here we implement a flexible fitting routine of time-domain NIRS data using graphics processing unit based Monte Carlo simulations. We compare the results for two different geometries: a two-layer slab with variable thickness of the first layer and a template atlas head registered to the subject’s head surface. We characterize the performance of the Monte Carlo approaches for fitting the optical properties from simulated time-resolved data of the adult head. We show that both geometries provide better results than the commonly used homogeneous model, and we quantify the improvement in terms of accuracy, linearity, and cross-talk from extracerebral layers.

## 1.

## Introduction

Various continuous-wave (CW),^{1}2.3.4.^{–}^{5} frequency-domain (FD),^{6}7.8.9.10.11.12.13.^{–}^{14} and time-domain (TD)^{15}16.17.18.^{–}^{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}20.21.^{–}^{22}

Continuous-wave methods, such as broadband or hyperspectral approaches originally proposed more than 15 years ago,^{2} require spatially^{1} or spectrally^{2}3.4.^{–}^{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^{6} has been extensively validated with Monte Carlo simulations,^{14} 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.^{22} In adults, the TD-NIRS technology has been applied to monitor variations in brain oxygenation and blood volume during cardiopulmonary bypass surgery^{25} or to detect vasospasms following subarachnoid hemorrhage.^{21}

Most of the *in vivo* studies mentioned above have relied on a simple head model described as a homogeneous semi-infinite medium.^{1}^{,}^{6}7.8.9.^{–}^{10}^{,}^{16}^{,}^{19}20.21.^{–}^{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.^{12} have shown, with simulations and phantom measurements in a slab geometry, that when a superficial layer thicker than $\sim 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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.^{28} 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.^{17} 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).^{16} 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.^{29} 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.^{4} Kienle et al. derived an analytical solution to the diffusion equation in a two-layer geometry for TD and FD.^{17} 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.^{17} 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.^{30} Paralleling the TD approach, Hallacoglu et al. demonstrated the improvement afforded by the FD two-layer geometry^{17} with Monte Carlo and phantom data, and applied the technique to *in vivo* data obtained on the adult head.^{31} 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,^{34} 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.^{35} 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.^{36} 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,^{39} 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.^{40} 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.^{40}

Contrary to a previous study where we reported the performance of the FDMD approach at one specific location on the head,^{28} 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.^{40} 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.^{43}

The true optical properties of the head were defined as follows. The reduced scattering coefficient was set to $10\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\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\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ in steps of $0.05\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$, while the absorption in the extracerebral layer (combined scalp and skull) was set to 0.08, 0.1, or $0.2\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\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.^{37} 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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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.^{45} 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.

## 2.2.

### Fitting Procedure

The simulated data for each source location (four TPSFs) were fitted using a Levenberg-Marquardt algorithm for nonlinear iterative least squares minimization (function *lsqcurvefit* from MATLAB^{®}, The MathWorks Inc., Natick, Massachusetts). We applied the least square cost function to the square root of the data to accentuate the relative weight of later delays. The fitting parameters were the absorption coefficients of the medium (one or two parameters depending on the model geometry), the reduced scattering coefficient (one parameter), and one scaling factor for each of the four source-detector separations. As we will describe later, in the case of the two-layer slab, the thickness of the first layer was not directly fitted for. Instead, we looped through all fixed values of the thickness and estimated the best value based on the minimal fit residual (defined as the minimization cost function, i.e., the sum of squared differences applied to the square root of the data). We fitted the TPSFs at the four source-detector separations simultaneously. The fitting range extended from 50% of the peak on the rising edge to 0.01% of the peak on the tail of each TPSF. We compared different forward models of light propagation, namely the analytical solution of the diffusion equation for a semi-infinite homogeneous medium, and Monte Carlo simulations of radiative transport for (1) a two-layer slab geometry with varying thickness of the first layer and (2) a generic atlas head geometry.

## 2.2.1.

#### Analytical homogeneous model

We first fitted the data with the analytical solution of the diffusion equation for a semi-infinite homogeneous medium, as is most commonly employed in the literature. This serves as a reference, and the performance of the other fitting approaches will be quantitatively characterized relative to this model. We used the expression first derived by Patterson et al.^{15} with an extrapolated boundary condition as detailed by Kienle et al.^{23} The direct line source–detector distances (i.e., without taking into account the head curvature) were rounded off to the closest millimeter.

## 2.2.2.

#### Monte Carlo models

While numerical approaches are slower than analytical ones, GPU-based Monte Carlo methods have reduced the computation times by three orders of magnitude^{37}^{,}^{38} compared to traditional CPU-based numerical simulations, rendering it a viable option for routine use in optical property fitting. Because the history of each detected photon (i.e., its partial path length in all tissue types) is saved, the detected fluence can be recomputed at every iteration of the fitting process for a new set of absorption values^{45} without the need to launch a new Monte Carlo simulation. On the contrary, fitting for scattering coefficient values requires either that a new simulation be run at each iteration of the fitting process or that the result be obtained from a library of prerun simulations for all possible scattering combinations. For both geometries we implemented (two-layer slab and atlas head), we generated a library of results for different scattering values.

### Monte Carlo fit on two-layer slab model

The two-layer slab model provides a simple medium geometry to account for the nonhomogeneous structure of the head in the absence of further information. For this slab geometry, all simulations were run beforehand, providing a library of data that are uploaded during the fitting process. One advantage of this approach is that the modeled geometry is independent of the true structure of the subject, so that the library of Monte Carlo data is only created once and can be used for any future subject.

We defined a large slab of lateral dimensions 180 mm by 180 mm and thickness 100 mm in order to minimize boundary effects. The first 20 mm from the surface are split in 20 tissue types, each one a layer of 1 mm. One source is located at the center of the slab surface and 64 detectors of 0.5 mm radius are located every 1 mm from 10 to 40 mm, along the $x$ and $y$ axes in both directions. In the fitting process, the real source-detector distances were therefore rounded to the closest multiple of 1 mm. The signals detected on detectors located at the same distance from the source (four of them) are subsequently averaged to improve the signal-to-noise ratio (SNR). The reduced scattering coefficient is homogeneous across all layers and varied homogeneously between 6 and $16\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$, in steps of $2\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$. One MCX simulation with ${10}^{9}$ photons was run for each ${{\mu}_{s}}^{\prime}$. Note that even though all layers were characterized by the same absorption for the simulations, the absorption can be modified *a posteriori* in each layer since the path length traveled by all detected photons has been recorded in a history file. We therefore created a flexible library of resulting history files to model a two-layer slab. The thickness of the superficial layer and the absorption coefficients of the two layers, ${\mu}_{a1}$ and ${\mu}_{a2}$, respectively, can be modified in the postprocessing of a single Monte Carlo simulation result: we assign absorption ${\mu}_{a1}$ to the first ${\mathrm{N}}_{1}$ tissue types (layers) and absorption ${\mu}_{a2}$ to the remaining layers, thus creating a two-layer slab with a first layer of thickness ${N}_{1}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$. Using a homogeneous scattering coefficient over the whole slab allows faster simulations and fitting process. For a first layer thickness typical of real head anatomy, the effect of the second layer ${{\mu}_{s}}^{\prime}$ is negligible on the resulting TPSF. We will review this assumption in more detail in the Discussion section below.

Even though it is feasible, we did not fit for ${\mu}_{a}$ and ${{\mu}_{s}}^{\prime}$ simultaneously in the Levenberg-Marquardt routine. Instead, we fix the value of ${{\mu}_{s}}^{\prime}$, fit for ${\mu}_{a}$, record the resulting fit residual, and repeat the process for all ${{\mu}_{s}}^{\prime}$. We then select the ${{\mu}_{s}}^{\prime}$ corresponding to the lowest residual. For a fixed ${{\mu}_{s}}^{\prime}$, the fit for ${\mu}_{a}$ takes $\sim 8\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$ because the large history file needs only be loaded once. On the contrary, when fitting ${\mu}_{a}$ and ${{\mu}_{s}}^{\prime}$ simultaneously, a new history file needs to be uploaded at each iteration of the nonlinear fit, each loading taking up $\sim 4$ to 5 s. In the following results, we looped through six values of ${{\mu}_{s}}^{\prime}$ for the fits from 6 to $16\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ in steps of $2\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$, for a total fitting time $<1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{min}$. Note that we limited our search to six scattering values because our goal here was to perform a time-consuming systematic study, where $\sim 500$ source locations and 18 absorption combinations were fitted for, or a total of $\sim \mathrm{10,000}$ individual fits per subject. In the case of a real-life experiment, where only a few datasets are recorded and need to be fit, a more refined search for the optimal ${{\mu}_{s}}^{\prime}$ can easily be implemented and would result in a reasonable fitting time of a few minutes.

### Monte Carlo fit on human brain atlas model

The use of a template atlas head in place of the true subject’s anatomy has shown great promises for reconstructing functional NIRS activation data.^{40}^{,}^{44}^{,}^{46} It allows the incorporation of a realistic brain structure without the need for costly individual MRI scans. If the optode locations have been registered with respect to specific superficial landmarks, such as the 10–20 reference points, the atlas head can be registered onto the subject’s head using these superficial references, as was detailed in Singh et al.^{43} The optical reconstruction can then be performed on this registered atlas with a known internal structure. Here we investigate the benefits of translating this atlas approach to characterizing the baseline optical properties of the brain. We are particularly interested in investigating whether it is beneficial to employ a realistic generic head structure, albeit not the true one, instead of a simpler two-layer model in cases where the true structure of the subject’s head is unknown.

We used the high-resolution Colin27 digital brain phantom described by Collins et al.^{47} as the segmented atlas volume. This anatomical atlas was first registered to each subject’s head using an affine transformation of 33 superficial landmarks, from the 10–20 reference points of the atlas onto the corresponding 10–20 reference points of each subject, using the method described by Singh et al.^{43} In a real-case experiment, the coordinates of these landmarks can be recorded using a 3-D digitizer. This registration process also adjusts the probe on the registered atlas surface. Monte Carlo simulations were then performed on the registered atlas head using the transformed probe coordinates. Similar to the data simulations, we segmented the atlas head into two tissue types: intra- and extracerebral tissue.

As for the slab, we kept ${{\mu}_{s}}^{\prime}$ homogeneous throughout the head, assuming that the contribution of the brain’s scattering is negligible, and we ran simulations for each optode and each ${{\mu}_{s}}^{\prime}$ between 6 and $16\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ in steps of $2\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$. Unlike the slab geometry, the atlas procedure requires new simulations to be run for each new subject since the atlas structure first undergoes an affine transformation to match the subject’s superficial landmarks. Remember, however, that each simulation takes only $\sim 5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{min}$ per scattering coefficient value. Therefore, the creation of a subject-specific library takes only 1 h per source for 12 scattering values (or less if simulations are run on parallel processors).

## 2.3.

### Performance Metrics of the Data Fitting

The performance of each fitting approach was characterized by four metrics of the retrieved brain absorption coefficient: the median relative error, the linearity (slope and ${R}^{2}$), and the cross-talk from extracerebral absorption. Each performance metric was computed at individual optode locations, and we present the median, 10-, 25-, 75-, and 90-percentiles over all locations and all three subjects.

## 2.3.1.

#### Relative error on brain absorption

At each location, we considered the retrieved brain absorption for all 12 combinations of the true brain absorption varying between 0.05 and $0.3\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ in steps of $0.05\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ and extracerebral absorption of 0.1 or $0.2\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$. We then computed the median of the absolute value of the relative error $|{\mu}_{a,\text{Retrieved}}-{\mu}_{a,\text{true}}|/{\mu}_{a,\text{True}}$ over all 12 values. This first metric characterizes the accuracy we can expect on the absolute values of the retrieved brain absorption.

## 2.3.2.

#### Linearity

For extracerebral absorption fixed at $0.1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$, and brain absorption varying between 0.05 and $0.3\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ in steps of $0.05\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$, we performed a linear fit of the retrieved brain absorption coefficient and characterized the linearity by the slope and the degree of freedom adjusted ${R}^{2}$ of the linear fit. The linearity metrics provide a measure of how accurately we can detect changes in the brain absorption for constant extracerebral contamination. Even if the retrieved absolute value for brain absorption is inaccurate, it can nonetheless be useful to detect relative changes accurately in individual subjects, for instance, to look at variations over time during an intervention or over several days during treatment.

## 2.3.3.

#### Cross-talk

For each fixed value of the brain absorption between 0.05 and $0.3\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$, we computed the relative change in the retrieved brain absorption for extracerebral absorption varying by 20% from 0.1 to $0.08\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$. The cross-talk is expressed as the corresponding percent change in retrieved ${\mu}_{a}$ averaged over all brain absorption values. We excluded from this computation cases where the retrieved brain absorption reached the lower or upper fitting boundaries (0.02 or $0.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$), which would artificially lead to null cross-talk. This metric is important to characterize how variations in the skin absorption will contaminate the retrieved brain absorption.

## 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\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ and brain absorption varying between 0.05 and $0.3\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ by steps of $0.05\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\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.

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.^{48} 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.

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.

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.1.

### Benefits of Monte Carlo-Based Fitting Approach

Using the Monte Carlo fitting approach combined with a layered model of the head (slab or atlas) improved the reconstruction for all metrics compared with a homogeneous medium. The improvement in terms of relative error is relatively modest, leading to a reduction in the error of 30 to 25% on average for slab and atlas, respectively. This value varies within each subject depending on the location, as well as between subjects. Cross-talk and linearity show a more dramatic improvement. While accuracy in the absolute values is a requirement for comparison of brain parameters between subjects, cross-talk and linearity are also essential parameters when measurements are to be performed within one subject over time. In this case, higher errors on the absolute values can be tolerated, provided that relative changes can be accurately detected with minimal contamination by changes in extracerebral tissue. These results suggest that it is beneficial to use a Monte Carlo-based approach combined with a two-tissue-type model of the head when fitting adult head data.

The benefits of the layered geometries, as well as the choice of the best model, appear to be subject-dependent. For instance, for subject 1, the two-layer slab reduces the average error by $\sim 50\%$ (from 25% error with homogeneous down to 14 and 13%). In contrast, for subject 3, the use of the Monte Carlo fitting does not improve the accuracy compared with a homogeneous model (error increasing slightly from 14 to 15%). Note that both cross-talk and linearity are improved, nonetheless, for subject 3. We postulate that the main contribution to this variability is the difference in the extracerebral layer thickness between each subject and the atlas head. This would explain why the two-layer slab, even though being a cruder representation of the head geometry and, in particular, ignoring its surface curvature, can perform better than the atlas by allowing a varying thickness of extracerebral tissue. While we are not directly fitting for the thickness of the extracerebral layer, we can estimate it by selecting the fit of smaller residual over all thicknesses. A finer sampling of the first layer thickness and of scattering coefficient values could easily be incorporated in this procedure, but was not implemented here because it would be too time-consuming in this systematic study. In our ongoing work, we will investigate the option to combine elements of both models, for instance, by performing Monte Carlo on ellipsoid-shaped medium with an extracerebral layer of varying thickness or by scaling the internal structure of the atlas brain.

Note that the Monte Carlo approach proposed here in combination with a layered slab or atlas head geometry could similarly be applied to spatially and/or spectrally resolved CW-NIRS data^{1}2.3.4.^{–}^{5} and to FD multiple-distance NIRS.^{6}7.8.9.10.11.12.13.^{–}^{14} We focused this study on evaluating the method for TD-NIRS data, in comparison with the traditional homogeneous model. Comparing the performance of this approach for different NIRS modalities (CW, FD, and TD), incorporating realistic noise characteristic for each technology, will require further studies beyond the scope of the present manuscript.

## 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\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\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\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\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\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ [Fig. 5(e)] and almost unnoticeable over the displayed five orders of magnitude for a higher absorption of $0.2\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ [Fig. 5(f)]. Furthermore, note that the range of scattering values we studied (4 to $20\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\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.

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.^{36} Instead, we used a library approach as was previously employed by different groups.^{35} Sampling ${{\mu}_{s}}^{\prime}$’ from 6 to $16\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-1}$ by steps of $2\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\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.

## 4.3.

### Limitations of the Present Study

## 4.3.1.

#### Monte Carlo noise versus experimental noise

Using simulated data enables us to create realistic data for known optical properties of the head. One limitation of this approach is that we restricted the noise of the data to that inherent to the Monte Carlo simulations. The Monte Carlo simulations yielded photon counts of 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, respectively (average over all subjects and all locations). The corresponding noise level at the peak of each TPSF was approximately $0.3\%$, 1.5%, 4%, and 7% of the signal. It is typical for *in vivo* TD-NIRS experiments to report photon counts from ${10}^{5}$ to ${10}^{6}$ per TPSF for integration times varying from a few hundreds of milliseconds to a few seconds.^{16}^{,}^{52}^{,}^{53} Since both Monte Carlo noise and experimental shot noise vary as the square root of the signal, our simulated SNR at 1 and 2 cm is realistic of experimental conditions, while the Monte Carlo SNR at 3 and 4 cm is lower than what is typically obtained in *in vivo* conditions. However, the Monte Carlo SNR at long delays might be optimistic compared to experimental conditions where other instrumental sources of noise, such as background noise arising from dark counts and stray light, deteriorate the signal. In the inverse problem, we fit the TPSF down to 0.01% of the peak, thus assuming that experimental conditions permit the recording of four orders of magnitude without reaching noise level. The range of TPSF time delays we used (from 50% of the peak on the rising edge to 0.01% of the peak on the tail) is typical of what is presented in the literature.

In summary, the noise inherent to our Monte Carlo simulations is similar to experimental conditions, but with differences at large source-detector separations and long delay times. These differences are likely to influence the uncertainty of the retrieved absorption values, but we did not investigate the effects of different noise levels or of different temporal ranges for the fitting, which could be the subject of further work.

## 4.3.2.

#### Head structure simplification

By using true head structures obtained from segmented MRI scans, we were able to simulate data that are more realistic than those obtained from homogeneous or two-layer slab geometries, in particular, by taking into account the true head curvature and the spatially varying thickness of extracerebral tissue. However, our head model remains a simplification of the heterogeneity of real heads. In particular, we assumed homogeneous scattering over the different head tissues. The skin and the different bone layers were also combined into a single extracerebral tissue type. Finally, we did not investigate the deterioration of the fitting results when considering a low-absorbing and low-scattering CSF layer. Further studies will require the description of more tissue types with heterogeneous optical properties.

## 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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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.