Quantifying tissue optical properties of human heads in vivo using continuous-wave near-infrared spectroscopy and subject-specific three-dimensional Monte Carlo models

Abstract

Significance: Quantifying subject-specific optical properties (OPs) including absorption and transport scattering coefficients of tissues in the human head could improve the modeling of photon propagation for the analysis of functional near-infrared spectroscopy (fNIRS) data and dosage quantification in therapeutic applications. Current methods employ diffuse approximation, which excludes a low-scattering cerebrospinal fluid compartment and causes errors.

Aim: This work aims to quantify OPs of the scalp, skull, and gray matter in vivo based on accurate Monte Carlo (MC) modeling.

Approach: Iterative curve fitting was applied to quantify tissue OPs from multidistance continuous-wave NIR reflectance spectra. An artificial neural network (ANN) was trained using MC-simulated reflectance values based on subject-specific voxel-based tissue models to replace MC simulations as the forward model in curve fitting. To efficiently generate sufficient data for training the ANN, the efficiency of MC simulations was greatly improved by white MC simulations, increasing the detectors’ acceptance angle, and building a lookup table for interpolation.

Results: The trained ANN was six orders of magnitude faster than the original MC simulations. OPs of the three tissue compartments were quantified from NIR reflectance spectra measured at the forehead of five healthy subjects and their uncertainties were estimated.

Conclusions: This work demonstrated an MC-based iterative curve fitting method to quantify subject-specific tissue OPs in-vivo, with all OPs except for scattering coefficients of scalp within the ranges reported in the literature, which could aid the modeling of photon propagation in human heads.

1. Introduction

Human brain activities have long been of great interest to many fields stretching from basic neuroscience, branches of medicine, and education to even brain-computer interface. Functional magnetic resonance imaging (fMRI) is the gold standard method of brain activity monitoring with great three-dimensional (3D) spatial resolution over the whole brain. However, the temporal resolution of fMRI is low, and the subject or patient is restricted to the supine position in an MRI scanner, which is expensive and not readily accessible. Functional near-infrared spectroscopy (fNIRS) employs low-cost and compact instruments to monitor brain activities with moderate temporal and spatial resolution. Since fNIRS measurements can be taken while the subject or patient is in any position or even moving, their applications are more versatile than fMRI.

In fNIRS, red to near-infrared light in at least two wavelengths is shone on the head and intensity changes of the diffuse reflection light are measured. The cerebral cortex consumes oxygen during neuronal activation, which induces increases in local blood volume and blood flow. The concentrations of oxyhemoglobin (HbO) and deoxyhemoglobin (Hb) in the active area change consequently. Quantification of hemoglobin concentration changes ([Formula omitted. See PDF]) is achieved based on the modified Beer–Lambert law (MBLL)¹

[Formula omitted. See PDF] where [Formula omitted. See PDF] is the logarithm of the ratio of diffuse reflectance intensities between two-time points, PL is the light pathlength in the cerebral cortex, and [Formula omitted. See PDF] is the molar extinction coefficients of the hemoglobins. In essence, the concentration changes can be calculated from intensity measurements if the light pathlength in the corresponding tissue compartment can be quantified. Many fNIRS studies have used a pathlength calculated using a fixed differential pathlength factor,² which is the ratio between the pathlength in tissue and the distance between the source and detector at the scalp surface. Subject-specific pathlengths can be obtained based on a homogeneous tissue model and optical properties (OPs) including the absorption coefficient ([Formula omitted. See PDF]) and the transport scattering coefficient ([Formula omitted. See PDF])³ measured by for example time-resolved (TR)⁴ or frequency-domain (FD)⁵ reflectance spectroscopy.

Although MBLL is widely adopted to achieve acceptable results for qualitatively assessing the general trend in hemodynamic changes, its ability to accurately quantify the concentration changes is debatable because the cerebral cortex, the target region of fNIRS, only accounts for a small portion of light attenuation in the head, which is known as the partial volume effect.^6–8 To separate the cerebral contribution from the typically stronger influence of the extra-cerebral tissue on the measured intensity changes, it is important to separately estimate partial pathlengths in the brain and the extra-cerebral tissue compartment,^9,10 which can be achieved by Monte Carlo (MC) simulations of light propagations in a realistic head model such as that segmented from MRI scans of the head.¹¹ However, OPs of major tissue types in the head measured from ex-vivo tissue samples show wide variations among studies,¹² which can be attributed to different sample preparation procedures¹³ and intersubject variability. Moreover, there were substantial intersubject variations in pathlengths estimated by finite-element calculations of light propagation using the same OPs on 40 subject-specific head models.¹⁴ Therefore, it is imperative to develop a noninvasive method for quantifying tissue OPs in the head of individual subjects to account for intersubject variabilities.

Precisely measuring head tissue OPs for individual subjects can also help other applications. For example, photobiomodulation is a technique that uses red or near-infrared light to stimulate or even heal tissue and has been applied to the human brain against dementia, Alzheimer’s disease, and Parkinson’s disease. The beneficial effect of photobiomodulation shows the biphasic dose response,¹⁵ which requires the dosage to be within a proper range. Research has been done to simulate the fraction of energy delivered to the human brain under transcranial stimulation,¹² and found that the dosage in the brain is greatly affected by the head tissue OPs.

NIRS is often used to noninvasively measure the head OPs. The TR and FD NIRS systems have been used to measure head OPs in many studies,^16–18 but the instruments are more complicated and expensive than continuous-wave (CW) NIRS systems. CW NIRS systems have also been used to quantify OPs of piglets’ brains¹⁹ and human brain.²⁰ In these studies, the diffusion approximation of the radiative transfer equation for a semi-infinite homogeneous model¹⁹ or a two-layered (i.e., extra-cerebral and cerebral) model^16–18,20 was used to simulate diffuse reflectance spectra. A low-scattering cerebrospinal fluid (CSF) layer was not included in the models due to the limitation of diffusion approximation to highly scattering media. Indeed, using the diffusion approximation solution to simulate light propagation in a realistic head model containing a low-scattering CSF compartment could lead to significant errors.²¹ However, excluding the low-scattering CSF layer could also significantly affect the simulation results of photon propagation in the head.^22,23 The MC method, on the other hand, is considered the gold standard and has been implemented to simulate light propagation in realistic head models.²⁴ It has not been used in iteratively curve fitting of NIRS data for quantifying OPs in the human head due to its relatively long time to simulate photon propagation in a deep, complex model like the human head.

This study aimed to establish a practical procedure to estimate the OPs of cerebral and extracerebral tissue compartments in-vivo using MC-based iterative curve fitting of CW NIRS data. Realistic human head models consisting of six tissue compartments were constructed from MRI scans of each subject. To overcome the challenge of huge time cost in running MC simulations, we trained an artificial neural network (ANN) to replace MC simulations as the forward model that calculated diffuse reflectance from given OPs.²⁵ However, due to the number and ranges of OPs required for this study a very large number of MC-simulated reflectance was needed to train the ANN. Therefore, the efficiency of generating the training data was enhanced by several techniques to prevent prohibitively long simulation time. First, white MC (WMC) was used to quickly calculate the reflectance for various [Formula omitted. See PDF] from only one simulation result obtained using each combination of [Formula omitted. See PDF]. Second, WMC simulations were run for detectors with a high acceptance angle to increase their efficiency, and a linear regression model was built to convert the simulated reflectance to that corresponding to the actual acceptance angle of the detectors used to collect in-vivo NIRS data. Finally, the number of [Formula omitted. See PDF] combinations in the training data was increased by spline interpolation of simulated reflectance obtained from 2808 [Formula omitted. See PDF] combinations covering the reported ranges.¹² The trained ANN was six orders of magnitude faster than the MC simulation, making iterative curve fitting practical.

The fitting process was tested on simulated target spectra without noise and with proper noise to analyze the theoretical performance of the proposed method. In-vivo reflectance spectra were measured from five healthy volunteers with a broadband CW NIRS system built in-house, and [Formula omitted. See PDF] and [Formula omitted. See PDF] of the scalp, skull, and gray matter were quantified. Thanks to the highly efficient ANN forward model, uncertainty in the extracted OPs due to fluctuations in the measured in-vivo spectra were estimated.

2. Methods

A flow chart of the proposed method is shown in Fig. 1. T1-weighted head MRI scan was performed on nine subjects to obtain the anatomical structure of the head and build a 3D model for each subject. One ANN forward model was trained for each subject and used in iterative curve fitting to quantity tissue OPs of the scalp, skull, and gray matter. Due to changes in the NIRS system, only five of the nine subjects underwent in-vivo CW NIRS measurements on the right forehead and their tissue OPs were quantified using the proposed method. The human study was approved by the institutional review board of National Taiwan University, and all participants provided written informed consent.

2.1. Head Model Preparation

To make the simulation closer to reality, we performed T1-weighted MRI head scans on nine subjects to get anatomical structure specific to each subject. The voxel size was [Formula omitted. See PDF]. The MRI data were segmented using SPM²⁶ into five compartments including scalp, skull, CSF, gray matter (GM), and white matter (WM). The region with relatively low MRI signal intensity in the front head was considered as sinus and manually segmented to be the sixth compartment. One slice of the segmented head model and its corresponding MRI scan slice are shown in Figs. 2(a) and 2(b), respectively. The surface of the scalp was transformed into a mesh and the Mesh2EEG²⁷ function was performed to find locations of the EEG 10-5 system²⁸ on the surface. To measure reflectance spectra from the right forehead the source was placed at the EEG Fp2 position and the locations of the source and six detectors are shown in Fig. 2(c). The source-to-detector separation (SDS) for the six detectors was 0.8, 1.5, 2.12, 3, 3.35, and 4.5 cm, respectively.

OPs needed by the MC method to simulate photon energy propagation in tissue include absorption coefficient ([Formula omitted. See PDF]), scattering coefficient ([Formula omitted. See PDF]), refractive index ([Formula omitted. See PDF]), and anisotropic factor ([Formula omitted. See PDF]). Ranges of the OPs for each tissue type are listed in Table 1. The refractive index was assumed constant throughout the head to simplify MC simulations. The refractive index of the medium outside the head was set to 1.457, which equaled to that of detector fibers to simulate the exit angle of photons properly. The space between the source and detector fibers in the NIRS system was a soft pad with [Formula omitted. See PDF] around 1.42, which was close to the value used in the simulations.

Table 1

Ranges of OPs for each compartment or medium.

Tissue type	μa (1/cm)	μs′ (1/cm)	n	g
Outer medium	0	0	1.457	—
Scalp	0.1 to 0.6	5 to 35	1.4	0.9
Skull	0.05 to 0.45	5 to 35	1.4	0.9
CSF	0.015 to 0.1	1 to 3.7	1.4	0.9
GM	0.05 to 0.5	5 to 35	1.4	0.9
WM	0.025 to 0.25	15 to 105	1.4	0.9
Sinus	0	0	1	—

In the diffuse regime, effects of [Formula omitted. See PDF] and [Formula omitted. See PDF] on photon propagation in tissue can be combined into the transport scattering coefficient ([Formula omitted. See PDF]), which is expressed as

[Formula omitted. See PDF]

The anisotropy factor of tissue was assumed to be constant since its effect was considered by [Formula omitted. See PDF]. Wavelength-dependent [Formula omitted. See PDF] were assumed to follow an inverse power-law function of wavelength [Formula omitted. See PDF]²⁹ and calculated with coefficients [Formula omitted. See PDF] and [Formula omitted. See PDF] as

[Formula omitted. See PDF]

The wavelength dependence of tissue [Formula omitted. See PDF] was assumed to be known based on absorption coefficients of major tissue chromophores in the wavelength range analyzed in this study. Exploiting the wavelength dependences of [Formula omitted. See PDF] and [Formula omitted. See PDF] could improve the efficiency and robustness of the inverse model applied on broadband spectra since only a few unknown parameters need to be determined instead of two OPs for each wavelength. The [Formula omitted. See PDF] of a tissue compartment was calculated as

[Formula omitted. See PDF] where [Formula omitted. See PDF] and [Formula omitted. See PDF] are the molar extinction coefficient of Hb and HbO, respectively, [Formula omitted. See PDF] is the tissue oxygen saturation, tHB is the total hemoglobin concentration, [Formula omitted. See PDF] is the absorption coefficient of 100% (v/v) of the [Formula omitted. See PDF]’th absorber other than hemoglobin, and [Formula omitted. See PDF] is the volume fraction of the absorber. The [Formula omitted. See PDF] and [Formula omitted. See PDF] spectra used in this study are shown in Figs. 3(a) and 3(b).

[Formula omitted. See PDF] and [Formula omitted. See PDF] spectra of the CSF were assumed to have fixed values adopted from the literature as shown in Figs. 3(c) and 3(d) since they are relatively small and barely affect the reflectance spectra. The [Formula omitted. See PDF] of WM was set to be half of the [Formula omitted. See PDF] of GM, and the [Formula omitted. See PDF] of WM was set to be 3 times the [Formula omitted. See PDF] of GM.³⁵ Since WM depths in the segmented models are between 16 and 20 mm, the approximation of WM OPs was not expected to influence the results of quantified tissue OPs. The tissue parameters used to calculate [Formula omitted. See PDF] and [Formula omitted. See PDF] of the five tissue compartments are summarized in Table 2.

Table 2

Tissue parameters used to calculate OPs of the five major tissue compartments.

Tissue	Variable tissue parameters	Fixed tissue parameters
Scalp	[Formula omitted. See PDF], [Formula omitted. See PDF], [Formula omitted. See PDF], [Formula omitted. See PDF], [Formula omitted. See PDF]	[Formula omitted. See PDF], [Formula omitted. See PDF]
Skull	[Formula omitted. See PDF], [Formula omitted. See PDF], [Formula omitted. See PDF], [Formula omitted. See PDF]
CSF	None	Literature reported values
GM	tHB, [Formula omitted. See PDF], [Formula omitted. See PDF], [Formula omitted. See PDF]	[Formula omitted. See PDF]
WM	None	[Formula omitted. See PDF], [Formula omitted. See PDF]

2.2. Monte Carlo Simulation and Acceleration

An open-source graphics-processing units (GPU) accelerated MC simulation tool, Monte Carlo eXtreme (MCX),³⁶ was used to perform simulations of photon propagation in voxel-based head models. To accelerate the simulations, we performed WMC in which [Formula omitted. See PDF] of all tissue types were set to zero and the pathlength of each detected photon in each tissue compartment was recorded. After the simulation with one set of [Formula omitted. See PDF] was done, the reflectance of tissue models with the same set of [Formula omitted. See PDF] and any [Formula omitted. See PDF] combination could be quickly calculated by summing the attenuated weight of every detected photon according to microscopic BLL³⁷

[Formula omitted. See PDF] where [Formula omitted. See PDF] is the calculated reflectance, [Formula omitted. See PDF] is the number of detected photons, [Formula omitted. See PDF] is the absorption coefficient of tissue compartment [Formula omitted. See PDF], [Formula omitted. See PDF] is the pathlength of the [Formula omitted. See PDF]’th photon in tissue compartment [Formula omitted. See PDF] from the simulation, and [Formula omitted. See PDF] is the number of launched photons. WMC simulations were performed using 2808 [Formula omitted. See PDF] combinations whose values for each of the four tissue compartments are listed in Table 3. The number of sampled [Formula omitted. See PDF] values for each tissue compartment was determined empirically and roughly followed the sensitivity of the reflectance to each compartment’s [Formula omitted. See PDF]. Results of WMC simulations, i.e., [Formula omitted. See PDF] and [Formula omitted. See PDF] of detected photons, were recorded in a lookup table, which contained 2808 entries corresponding to the simulated [Formula omitted. See PDF] combinations.

Table 3

μs′ of each tissue compartment used to construct the lookup table.

Tissue	μs′ (1/cm)	Number of points
Scalp	5, 7.5, 10, 12.5, 15, 17.5, 20, 22.5, 25, 27.5, 30, 32.5, 35	13
Skull	5, 8.75, 12.5, 16.25, 20, 23.75, 27.5, 31.25, 35	9
CSF	1, 1.9, 2.8, 3.7	4
GM	5, 11, 17, 23, 29, 35	6

To further reduce the number of launched photons and hence the time needed by MC simulations while keeping the results converged, we set the numerical aperture (NA) of detector fibers in MC simulations to be 1.0, which collected significantly more photons than the [Formula omitted. See PDF] detector fibers actually used in the CW NIRS system described in Sec. 2.3. An NA-conversion linear regression model was developed to convert simulation results of reflectance collected by [Formula omitted. See PDF] fibers into those collected by [Formula omitted. See PDF] fibers. The resultant NA-conversion model took six inputs including logarithm of the reflectance, logarithm of the number of detected photons, [Formula omitted. See PDF], [Formula omitted. See PDF], [Formula omitted. See PDF], and [Formula omitted. See PDF]. The output of the NA-conversion model was the ratio between the reflectance collected by [Formula omitted. See PDF] fibers and that by [Formula omitted. See PDF] fibers. Each data point in the lookup table was simulated using MCX with [Formula omitted. See PDF] photons and a detector [Formula omitted. See PDF].

After MC simulations were done, we used the lookup table to generate data for training ANN forward models. We randomly chose 3000 [Formula omitted. See PDF] combinations within the ranges reported in Table 1 and calculated the corresponding reflectance using Eq. (5) for each [Formula omitted. See PDF] combination in the table. Then the reflectance simulated by [Formula omitted. See PDF] fibers was converted into reflectance detected by [Formula omitted. See PDF] fibers using the NA-conversion model. Since the error introduced by interpolation was rather small under these [Formula omitted. See PDF] sampling intervals, we could augment the amount of training data by interpolation between the simulated [Formula omitted. See PDF] values. Based on the 2808 [Formula omitted. See PDF] combinations in the lookup table, more reflectance values for other [Formula omitted. See PDF] combinations could be readily obtained through spline interpolation. We randomly chose 3000 [Formula omitted. See PDF] combinations and calculated their reflectance values, producing a total of [Formula omitted. See PDF] million OP combinations with reflectance values for training ANN forward models in a reasonable time.

The 17.4 million sets of OP/reflectance data for each subject were split into training, validation, and test datasets in fractions of 75%, 10%, and 15% respectively. Inputs to the ANN models were [Formula omitted. See PDF] and [Formula omitted. See PDF] of the four major tissue compartments in Table 3, and outputs were reflectance values received by the six detectors. We used a fully connected neural network with four hidden layers consisting of 850, 550, 300, and 150 neurons, respectively.

2.3. Broadband Near-Infrared Spectroscopy Instrument

We built a broadband CW NIRS system to measure in-vivo diffuse reflectance spectra at six SDSs from the forehead of five healthy volunteers. A schematic diagram of the system is shown in Fig. 4(a). We used a quartz tungsten halogen lamp (66997-250Q-R085, Newport Corporation, Irvine, California, United States) as the light source. Two filters were used to filter out wavelengths below 650 nm or above 1050 nm, and the remaining light was guided to the subjects by a bundle of 16 fibers whose NA was 0.39. The outer diameter of the active area of the source bundle was about 4.2 mm. The light remitted by tissue was collected by six detector fibers with a core diameter of 0.4 mm, NA of 0.12, and SDS of 0.8, 1.5, 2.12, 3, 3.35, and 4.5 cm, respectively. For convenient reference in the manuscript, the detectors were numbered as 1 to 6, respectively. A lens and an objective lens were used to image the end-face of the detector fibers to the entrance slit of an imaging spectrograph (HOLOSPEC-F/1.8- NIR, Oxford Instruments, Abingdon, Oxfordshire, United Kingdom), and spectra were measured by an electron-multiplying charge-coupled device (EMCCD) (Newton DU970P, Oxford Instruments) with a spectral resolution of about 4 nm. To make sure that the contact between the fibers and scalp was stable, we 3D-printed a curved holder and used springs as shown in Figs. 4(b) and 4(c) to press the fibers against the scalp with consistent pressure. The holder interfaced the scalp surface with a soft black pad made of polydimethylsiloxane (PDMS) to ensure stable and tight contact with the scalp surface. Ultrasound gel was smeared on the scalp near the measured site before the holder was mounted to prevent any air gap between the fibers and the scalp. The holder was placed on the subjects’ right forehead with the source fiber bundle pointing to the EEG Fp2 position.

We used the NIRS system to measure reflectance spectra from five healthy male subjects aging between 24 and 47. We repeated the procedure of applying the holder/fibers and acquiring spectra for at least three times on each subject to make sure that the contact between the fibers and scalp was stable and to estimate the variations in measured reflectance intensities. Coefficients of variation (CVs) of the repeatedly measured reflectance intensities at each SDS were averaged over all wavelengths and subjects, and results are listed in Table 4.

Table 4

Coefficients of variation (CV) of repeated in-vivo reflectance spectra measurements on each subject. The numbers shown are the average over all wavelengths and subjects.

SDS (cm)	0.8	1.5	2.12	3	3.35	4.5
CV (%)	3.0	4.2	5.1	5.2	5.4	12.1

2.4. Iterative Curve Fitting

Measured in-vivo reflectance spectra were calibrated using six homogeneous phantoms made of PDMS, [Formula omitted. See PDF] and India ink. After [Formula omitted. See PDF] and [Formula omitted. See PDF] of the phantoms were measured by an integrating sphere and the inverse adding-doubling method, we simulated reflectance spectra of the phantoms with CUDAMCML³⁸ and performed calibrating linear regression³⁹ of the simulated reflectance against the measured reflectance of the phantoms. The resultant calibrating regression equation was applied to calibrate in-vivo reflectance spectra to the absolute reflectance, which were compared with MC simulated reflectance. The average of multiple measured spectra of the same subject after calibration was taken as the target spectra.

We used iterative curve fitting to find values for the variable tissue parameters listed in Table 2 that minimized the root-mean-square (RMS) percent error between simulated spectra and calibrated in-vivo or simulated test spectra. The spectral error was calculated as

[Formula omitted. See PDF] where [Formula omitted. See PDF] is the number of detectors, [Formula omitted. See PDF] is the number of wavelengths considered, [Formula omitted. See PDF] is the simulated reflectance in the current iteration, and [Formula omitted. See PDF] is the measured reflectance after calibration or simulated reflectance of the test spectra. We chose 22 wavelengths nonevenly distributed between 700 and 880 nm to perform the fitting. Initial values of the tissue parameters to be determined were chosen from a pool of 5000 sets of randomly selected parameters. Each set of measured spectra to be fitted was first compared with presimulated spectra corresponding to the 5000 sets of parameters, from which 20 sets of parameters resulting in the smallest spectral error to the target spectra were chosen as initial values. The target spectra were fitted 20 times to address the local minimum problem. In each iteration, the temporarily guessed tissue parameters were converted to OPs, which in turn were sent to the ANN model to get simulated spectra. The iterative curve fitting function “fmincon” with a sequential quadratic programming algorithm provided by MATLAB^® (MathWorks, Inc., Natick, Massachusetts, United States) was used to minimize the spectral error within 2000 iterations. During the fitting process, the tissue parameters were nonlinearly constrained to make sure that the corresponding OPs did not exceed the ranges shown in Table 1.

In our previous preliminary study,⁴⁰ the proposed method was validated to quantify OPs of a three-layered tissue phantom. Details are described in the Supplementary Material, and errors for extracted OPs of all layers were below 15%.

To assess the accuracy of the iterative curve fitting method using realistic head models, we generated 15 sets of test spectra using randomly picked tissue parameters and each of the nine subjects’ ANN model. We added to the test spectra 14 sets of random noise whose CV was similar to those estimated from in-vivo measurements to simulate the fluctuations of in-vivo spectra due to various sources including noise of the EMCCD readings, instability and variation of the probe-tissue coupling, and physiological changes of hemoglobin concentrations in the probed region. Therefore, there were 210 sets of test spectra with noise per subject, and 1890 sets of test spectra in total. Each set of the test spectra was fitted 20 times using different initial values as described above.

2.5. Multiple Solution

The fitting result with the smallest spectral error among the 20 times of fittings was typically chosen as the final result. However, since there was noise in the measured reflectance spectra, when multiple fitting results showed spectral errors to be within the noise level of each other, we could not verify which one was the best solution. The noise level of the NIRS system was estimated to be 2% by measuring static tissue phantoms presenting reflectance intensities similar to in-vivo forehead. When the difference between spectral errors of multiple fitting results was [Formula omitted. See PDF], we considered those fitting results as potential multiple solutions. To reduce the number of multiple solutions and increase the chance of getting a unique solution of OPs, we employed spectra measured from detectors that were not used in the iterative curve fitting as tie breakers (see Table 9 for detector combinations chosen for fitting in-vivo spectra). An example of choosing the best solution is shown in Fig. 5. Assume only spectra of detectors one to four were used in the fitting. We compared spectral errors of the 20 fittings and chose the results whose spectral errors were within 2% of the smallest error. Then we compared the spectral error of detectors not included in the fitting, i.e., detectors 5 and 6, and again chose the results with errors within 2% of the lowest error. If two or more solutions were chosen, they were considered multiple solutions.

3. Results

3.1. Performance of the ANN Model

Using a larger detector fiber NA to run MC simulations increases the number of photons collected by the fiber and hence the simulation efficiency. By using the NA-conversion model described in Sec. 2.2, we could reduce the launched photon number by about 10 times while keeping the noise at a similar level as shown in the first and third rows in Table 5.

Table 5

Comparison of CVs of MC-simulated reflectance between NA = 1.0 with NA-conversion and NA = 0.12.

Photon number	NA of fiber	Conversion to low NA	SDS (cm)
			0.8	1.5	2.12	3	3.35	4.5
[Formula omitted. See PDF]	1.0	Yes	0.42%	0.80%	1.42%	3.60%	4.15%	8.21%
[Formula omitted. See PDF]	0.12	No	0.87%	2.60%	2.51%	7.21%	8.42%	25.80%
[Formula omitted. See PDF]	0.12	No	0.45%	0.87%	0.86%	1.79%	2.64%	8.30%
[Formula omitted. See PDF]	0.12	No	0.10%	0.19%	0.19%	0.40%	0.59%	1.86%

The WMC simulations of 2808 [Formula omitted. See PDF] combinations in the lookup table took about 24 h using seven NVIDIA GeForce RTX 2080. The generation of the 17.4 million training data from lookup table took 3 h. The training of ANN was done in 2 h. To assess the combined error introduced by the lookup table, NA-conversion model and trained ANN forward model, we simulated three spectra, each consisting of 35 wavelengths, for two subjects using MC simulations with [Formula omitted. See PDF] photons and [Formula omitted. See PDF] detector fibers to serve as ground-truth reflectance values. The same OPs were input into the trained ANNs and resultant reflectance values were compared with the MC-simulated ground-truth values. The RMS spectral error for each detector SDS is listed in Table 6. We found that the error between the ANN output and the ground-truth MC simulation was comparable to the noise of MC simulations with about [Formula omitted. See PDF] photons (see Table 5), while the speed of ANN was six orders of magnitude faster than MC simulations. For example, it took about 50 min to run MCX simulations with [Formula omitted. See PDF] photons each to generate reflectance spectra consisting of 22 wavelengths, while using the ANN only took 2.6 ms.

Table 6

RMS percent error between MC-simulated and ANN-generated spectra average across two subjects and three spectra.

SDS (cm)	0.8	1.5	2.12	3	3.35	4.5
Error	2.00%	4.74%	5.02%	3.29%	3.03%	7.89%

Instead of training ANN forward models, the OP/reflectance data generated by MC simulations described in Sec. 2.2 could also be used to build an eight-dimensional lookup table. We tested using linear interpolation of the lookup table as the forward model, and found that spectral errors of the lookup table-based forward model were significantly larger than those of the ANN models for SDS equal to or longer than 3 cm, while the speed was forty times slower than ANN.

3.2. Test Spectra Fitting Results

Fitting one set of test spectra with 20 different initial values was done in 4 min. Relative errors in each of the six OPs estimated from iterative curve fitting of 1890 test spectra, 22 wavelengths per spectrum, were compiled to form a histogram for one OP as shown in Fig. 6. Multiple solutions as defined in Sec. 2.5, if existed, were included as separate data points. The histogram of errors could be taken as an estimation of the probability density function of the fitting errors. We defined the range of errors that covered 68% or 95% of the fitting results as confidence intervals of each estimated OP. The errors in OPs corresponding to 68% confidence intervals for each detector combination and OP are listed in Table 7.

Table 7

Ranges of errors (%) in the fitted OPs corresponding to 68% confidence intervals.

Detector combinations	μa,scalp	μs,scalp′	μa,skull	μs,skull′	μa,GM	μs,GM′
346	−19∼21	−9∼12	−27∼42	−15∼19	−30∼25	−28∼118
1234	−14∼20	−7∼7	−29∼41	−14∼18	−28∼25	−27∼111
2345	−17∼21	−9∼10	−26∼38	−14∼18	−29∼25	−28∼99
12345	−15∼20	−7∼8	−27∼40	−14∼17	−27∼23	−28∼89
123456	−15∼17	−7∼7	−19∼42	−15∼12	−28∼21	−27∼78

To compare errors between quantified OPs using subject-specific models and a general model, we used Colin 27 average brain model⁴¹ to train an ANN. The 1890 test spectra were fitted using the Colin 27 ANN and OPs’ errors were calculated. Table 8 shows a comparison of ranges of errors in the fitted OPs between each subject’s own ANN model and the Colin 27 ANN model. The ranges of errors are significantly larger for the Colin 27 model. We also fitted in-vivo spectra with the Colin 27 ANN model, and fitted OPs were different from those fitted with each subject’s own ANN even considering the OPs’ error ranges. The results suggest that the subject-specific model is needed for more accurate quantification of tissue OPs.

Table 8

Comparison between fitting test spectra with subject-specific models and the general Colin 27 model. Ranges of errors (%) in the fitted OPs correspond to 68% confidence intervals. Spectra of detectors 1 to 6 were used.

Fitting model	μa,scalp	μs,scalp′	μa,skull	μs,skull′	μa,GM	μs,GM′
Each subject’s own model	−15∼17	−7∼7	−19∼42	−15∼12	−28∼21	−27∼78
Colin 27	−32∼92	−34∼13	−1∼162	−42∼70	−44∼57	−73∼42

3.3. In-Vivo Fitting Results

When fitting in-vivo spectra, we found that in some subjects not all detectors showed a small (< 15%) spectral error. This might be caused by mismatches between actual human heads and the 3D tissue models built. The model assumes homogeneous OPs within each compartment, while in reality, there could be some inhomogeneities such as blood vessels under a detector or multiple detectors, making it hard to achieve small spectral errors for all detectors simultaneously. Table 9 shows spectral errors for each subject when spectra from different detector combinations were fitted. We chose the combination that included the most detectors and achieved a spectral error below 15% as the final fitting result. The chosen detector combination for each subject is marked in bold in Table 9, and the corresponding modeled and measured spectra are shown in Figs. S1–S5.

Table 9

Average spectral errors of each subject’s fitting results using different detector combinations. Only errors contributed by the fitting detectors are counted. The detector combination chosen for each subject is shown in bold.

Detector combination	Subject
	1	2	3	4	5
346	4.9%	5.1%	5.1%	25.1%	22.0%
1234	7.1%	26.3%	13.2%	9.5%	7.2%
2345	11.9%	37.5%	14.8%	7.3%	18.4%
12345	11.9%	35.5%	19%	10.2%	17.9%
123456	11.3%	32.8%	19.3%	29.3%	29.2%

Fitting in-vivo spectra by the procedure described in Secs. 2.4 and 2.5 resulted in only one solution for all subjects except for subject 1. The two solutions of subject 1’s fitting results showed similar scalp OPs. In Fig. 7, we compare the skull OPs of the two solutions to literature-reported values. Solution 1 had a slightly smaller [Formula omitted. See PDF] compared with values reported in the literature, while solution 2 had a slightly smaller [Formula omitted. See PDF] compared with values reported in the literature. Since the [Formula omitted. See PDF] from both solutions were close to literature-reported values, we considered both of them as final results for subject 1. As a summary, OPs of the five subjects are compared with the literature-reported values in Fig. 8.

4. Discussion

4.1. Fitting Results

We validated the process of choosing the fitting result with the smallest spectral error out of the 20 fittings based on results of fitting the 1890 test spectra sets. For each test spectra set the 20 fitting results were ranked according to the spectral error. To compare errors of estimated OPs between the 20 fitting results we calculated an RMS percentage error of OPs estimated by each fitting over the 22 fitted wavelengths and six OPs, and normalized the RMS errors of the 20 fitting results to the range between 0 and 1. Normalizing the OP errors enabled combining results of the 1890 test spectra sets since errors in OPs varied substantially across the test data. Afterward, the normalized OP errors of the fitting result with the same spectral error ranking from each of the 1890 test spectra sets were compiled together. Figure 9 shows a boxplot of normalized OP errors from results of all the 1890 test spectra sets to compare the 20 fitting results ranked according to spectral errors. It can be seen that the normalized OP errors are smaller in fitting results with smaller spectral errors. The result suggests that choosing the fitting result with the smallest spectral error leads to overall smaller errors in the six estimated OPs.

To investigate sources of errors in quantified OPs, we also fitted the same 135 sets of simulated test spectra without added noise. Ranges of errors in OPs (shown in Table S2) were not significantly lower than those with added spectral noise (Table 7). Since the test spectra were generated using the same ANN as used in the iterative curve fitting, ideally a perfect fit should be achievable. However, the average spectral error of test spectra fitting results was 1.78% rather than 0%. This is probably caused by noise of MC simulations, which the trained ANN models also inherit. During iterative curve fitting the optimization algorithm needs to calculate partial derivatives of the reflectance with respect to individual OPs to find the global minimum. As shown in Fig. 10, the sensitivity of CW NIRS measurements to deep tissue OPs is very low, which means that the reflectance hardly changes given a small perturbation in one OP. When the noise in MC-simulated or ANN output reflectance is comparable to or even larger than the reflectance changes due to the small perturbation in an OP, the ability of iterative curve fitting to find the global minimum could be hindered. See Fig. S6 for an example of trends in the ANN output reflectance at the six SDSs when only one of the six OPs varies and the other OPs remain constant. It can be seen that the cases with low sensitivity (all SDSs in the case of [Formula omitted. See PDF] and short SDSs in cases of [Formula omitted. See PDF] and [Formula omitted. See PDF]) show rugged or not monotonic trends in the reflectance with varying OP values. If the sensitivity of measured NIRS data to deep tissue OPs could be improved, e.g., using TR or FD NIRS systems, the influence of MC noise and ANN output errors could be reduced, and the errors of the fitted OPs may be improved.

The comparison of the fitted OPs and the literature reported values is shown in Fig. 8. All OPs except for [Formula omitted. See PDF] are within the ranges reported in the literature. Although [Formula omitted. See PDF] is much lower than the literature reported values, [Formula omitted. See PDF] is not overestimated. Therefore, the underestimation of [Formula omitted. See PDF] is not due to cross-talk between [Formula omitted. See PDF] and [Formula omitted. See PDF]. Considering that the SDSs used in this study ranged from 0.8 to 4.5 cm and the spatial resolution of the 3D model was relatively low with a 0.93-mm edge length, the scalp was not segmented into finer layers such as epidermis, dermis, and subcutaneous tissue. The scattering phase functions used in MC simulations were set to the conventional Henyey–Greenstein phase function with a constant anisotropy factor. Quantification of [Formula omitted. See PDF] could be improved by detecting spectra at shorter SDS and using more realistic structural and optical parameters for skin.

It is worth mentioning that the quantified [Formula omitted. See PDF], [Formula omitted. See PDF], [Formula omitted. See PDF], and [Formula omitted. See PDF] show high variabilities among subjects even considering the uncertainties shown in Fig. 8. The confidence interval for [Formula omitted. See PDF] is too large, which prevents reliable quantification of [Formula omitted. See PDF]. The high intersubject variability underpins the importance of quantifying the tissue OPs of each subject to more accurately calculate partial pathlengths for fNIRS applications or determine the fraction of photon energy delivered to the brain for therapeutic applications.

4.2. Feasibility

The process of building a 3D head model is done manually and costs no more than 1 h. Simulations of 2808 [Formula omitted. See PDF] combinations in the lookup table for generating ANN training data cost about 24 h using 7 NVIDIA GeForce RTX 2080. The measurement of NIRS spectra can be done in 1 h and can be done at the same time while the simulation is running. The generation of training data from the lookup table can be done in 3 h. The training of ANN can be done in 2 h. The fitting for one target spectra can be done in 4 min. The whole quantifying process can be done within 2 days after taking the MRI scan.

4.3. Future Work

We use iterative curve fitting in this study, whose performance is affected by the smoothness of the ANN output. Random noise in MC simulation results could be reduced using filtering techniques.⁴⁸ Using other optimization algorithms, e.g., genetic algorithm, which do not rely on the gradient of the loss function might help reduce the number of multiple solutions.

The uncertainties of [Formula omitted. See PDF] and [Formula omitted. See PDF] are relatively large compared with the other OPs. To improve the sensitivity of the GM OPs, one may consider using TR or FD NIRS systems, as in Refs. 1617.–18. With the help of each subject’s own model and the precise MC simulations, the accuracy of quantified OPs could be improved.

The edge length of one voxel in our voxel-based model is 0.93 mm, which is relatively large compared with the 0.8 cm SDS. Tran et al.⁵¹ proposed a pipeline for generating mesh-based models from MRI images with smooth surfaces for each layer. With a mesh-based MC simulation tool, they showed that the travel distance of the photons in GM is overestimated compared with the voxel-based model. This may affect the simulated reflectance spectra. Switching to the mesh-based model may help improve the match between the real head and the model.

5. Conclusion

This study aimed to use a multidistance CW NIRS system to quantify the OPs of the scalp, skull and GM in the human head in-vivo. MC simulations were performed to generate reflectance spectra given OPs and a voxel-based tissue model of each subject’s head to improve the consistency between the subject and the model. Lookup tables and white MC method were employed to efficiently generate numerous data for training ANNs that replaced MC simulations in iterative curve fitting of measured reflectance spectra. The average error between ANN outputs and corresponding MC simulations was under 10%. Iterative curve fitting was used to optimize the solution. The curve fitting could be done in several minutes with the help of the ANNs. We performed fittings on thousands of simulated test spectra to estimate the confidence interval of each fitted OP. The OPs of the scalp, skull, and gray matter at the right forehead were quantified for five subjects. The scalp [Formula omitted. See PDF] was underestimated while the other OPs were within the ranges reported in the literature. Differences in estimated OPs between different subjects are larger than the confidence intervals, which suggests the importance of quantifying the OPs for each subject.

Disclosures

Authors have no conflicts of interest to declare.

Acknowledgments

We thank Ministry of Science and Technology of Taiwan for financial support (Grant No. 108-2221-E-002-075-MY3).

Tzu-Chia Kao received his BS degree in electrical engineering and his MS degree in biomedical electronics and bioinformatics, both from National Taiwan University, Taipei, Taiwan, in 2018 and 2021, respectively. During 2018–2021, he was a graduate research assistant with Biomedical Spectrum and Optical Imaging Laboratory, National Taiwan University. His research interests include Monte Carlo simulations of the photon-tissue interaction and diffuse reflectance spectroscopy.

Kung-Bin Sung received MS and PhD degrees in biomedical engineering from the University of Texas at Austin in 1999 and 2003, respectively. He worked as a research scientist at Intel Corporation from 2003 to 2006. He has been a faculty member of the Graduate Institute of Biomedical Electronics and Bioinformatics at National Taiwan University since 2006. His research focuses on developing noninvasive optical spectroscopy and microscopy techniques for early diagnosis of precancers and monitoring physiological parameters.

References

1. L. Kocsis, P. Herman and A. Eke, “The modified Beer–Lambert law revisited,” Phys. Med. Biol., 51 (5), N91 –N98 (2006). https://doi.org/10.1088/0031-9155/51/5/N02 PHMBA7 0031-9155

2. M. Ferrari and V. Quaresima, “A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application,” Neuroimage, 63 (2), 921 –935 (2012). https://doi.org/10.1016/j.neuroimage.2012.03.049 NEIMEF 1053-8119

3. S. Fantini et al., “Non-invasive optical monitoring of the newborn piglet brain using continuous-wave and frequency-domain spectroscopy,” Phys. Med. Biol., 44 (6), 1543 –1563 (1999). https://doi.org/10.1088/0031-9155/44/6/308 PHMBA7 0031-9155

4. D. T. Delpy et al., “Estimation of optical pathlength through tissue from direct time of flight measurement,” Phys. Med. Biol., 33 (12), 1433 –1442 (1988). https://doi.org/10.1088/0031-9155/33/12/008 PHMBA7 0031-9155

5. A. Duncan et al., “Optical pathlength measurements on adult head, calf and forearm and the head of the newborn infant using phase resolved optical spectroscopy,” Phys. Med. Biol., 40 (2), 295 –304 (1995). https://doi.org/10.1088/0031-9155/40/2/007 PHMBA7 0031-9155

6. M. Hiraoka et al., “A Monte Carlo investigation of optical pathlength in inhomogeneous tissue and its application to near-infrared spectroscopy,” Phys. Med. Biol., 38 (12), 1859 –1876 (1993). https://doi.org/10.1088/0031-9155/38/12/011 PHMBA7 0031-9155

7. K. Uludag et al., “Cross talk in the Lambert–Beer calculation for near-infrared wavelengths estimated by Monte Carlo simulations,” J. Biomed. Opt., 7 (1), 51 –59 (2002). https://doi.org/10.1117/1.1427048 JBOPFO 1083-3668

8. G. Strangman, M. A. Franceschini and D. A. Boas, “Factors affecting the accuracy of near-infrared spectroscopy concentration calculations for focal changes in oxygenation parameters,” Neuroimage, 18 (4), 865 –879 (2003). https://doi.org/10.1016/S1053-8119(03)00021-1 NEIMEF 1053-8119

9. W. B. Baker et al., “Pressure modulation algorithm to separate cerebral hemodynamic signals from extracerebral artifacts,” Neurophotonics, 2 (3), 035004 (2015). https://doi.org/10.1117/1.NPh.2.3.035004

10. F. Fabbri et al., “Optical measurements of absorption changes in two-layered diffusive media,” Phys. Med. Biol., 49 (7), 1183 –1201 (2004). https://doi.org/10.1088/0031-9155/49/7/007 PHMBA7 0031-9155

11. A. Custo et al., “Effective scattering coefficient of the cerebral spinal fluid in adult head models for diffuse optical imaging,” Appl. Opt., 45 (19), 4747 –4755 (2006). https://doi.org/10.1364/AO.45.004747 APOPAI 0003-6935

12. L.-D. Huang et al., “Simulation study on the optimization of photon energy delivered to the prefrontal cortex in low-level-light therapy using red to near-infrared light,” IEEE J. Sel. Top. Quantum Electron., 27 1 –10 (2021). https://doi.org/10.1109/JSTQE.2021.3051671 IJSQEN 1077-260X

13. H. Arslan and B. Pehlivanoz, “Effect of purification, dehydration, and coagulation processes on the optical parameters of biological tissues,” Chin. Opt. Lett., 19 (1), 011701 (2021). https://doi.org/10.3788/COL202119.011701 CJOEE3 1671-7694

14. K. Nakamura et al., “Estimation of partial optical path length in the brain in subject-specific head models for near-infrared spectroscopy,” Opt. Rev., 23 (2), 316 –322 (2016). https://doi.org/10.1007/s10043-016-0179-9 1340-6000

15. M. R. Hamblin, “Shining light on the head: photobiomodulation for brain disorders,” BBA Clin., 6 113 –124 (2016). https://doi.org/10.1016/j.bbacli.2016.09.002

16. J. Choi et al., “Noninvasive determination of the optical properties of adult brain: near-infrared spectroscopy approach,” J. Biomed. Opt., 9 (1), 221 –229 (2004). https://doi.org/10.1117/1.1628242 JBOPFO 1083-3668

17. A. Farina et al., “In-vivo multilaboratory investigation of the optical properties of the human head,” Biomed. Opt. Express, 6 (7), 2609 –2623 (2015). https://doi.org/10.1364/BOE.6.002609 BOEICL 2156-7085

18. B. Hallacoglu, A. Sassaroli and S. Fantini, “Optical characterization of two-layered turbid media for non-invasive, absolute oximetry in cerebral and extracerebral tissue,” PLOS ONE, 8 (5), e64095 (2013). https://doi.org/10.1371/journal.pone.0064095 POLNCL 1932-6203

19. H. Z. Yeganeh et al., “Broadband continuous-wave technique to measure baseline values and changes in the tissue chromophore concentrations,” Biomed. Opt. Express, 3 (11), 2761 –2770 (2012). https://doi.org/10.1364/BOE.3.002761 BOEICL 2156-7085

20. J. D. Veesa and H. Dehghani, “Hyper-spectral recovery of cerebral and extra-cerebral tissue properties using continuous wave near-infrared spectroscopic data,” Appl. Sci., 9 (14), 2836 (2019). https://doi.org/10.3390/app9142836

21. Y. Oki, H. Kawaguchi and E. Okada, “Validation of practical diffusion approximation for virtual near infrared spectroscopy using a digital head phantom,” Opt. Rev., 16 153 –159 (2009). https://doi.org/10.1007/s10043-009-0026-3 1340-6000

22. C. Mansouri et al., “Depth sensitivity analysis of functional near-infrared spectroscopy measurement using three-dimensional Monte Carlo modelling-based magnetic resonance imaging,” Lasers Med. Sci., 25 (3), 431 –438 (2010). https://doi.org/10.1007/s10103-010-0754-4

23. E. Okada and D. T. Delpy, “Near-infrared light propagation in an adult head model. I. Modeling of low-level scattering in the cerebrospinal fluid layer,” Appl. Opt., 42 (16), 2906 –2914 (2003). https://doi.org/10.1364/AO.42.002906 APOPAI 0003-6935

24. D. A. Boas et al., “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express, 10 (3), 159 –170 (2002). https://doi.org/10.1364/OE.10.000159 OPEXFF 1094-4087

25. S.-Y. Tsui et al., “Modelling spatially-resolved diffuse reflectance spectra of a multi-layered skin model by artificial neural networks trained with Monte Carlo simulations,” Biomed. Opt. Express, 9 (4), 1531 –1544 (2018). https://doi.org/10.1364/BOE.9.001531 BOEICL 2156-7085

26. W. Penny et al., Statistical Parametric Mapping: The Analysis of Functional Brain Images, , (2007).

27. P. Giacometti, K. L. Perdue and S. G. Diamond, “Algorithm to find high density EEG scalp coordinates and analysis of their correspondence to structural and functional regions of the brain,” J. Neurosci. Methods, 229 84 –96 (2014). https://doi.org/10.1016/j.jneumeth.2014.04.020 JNMEDT 0165-0270

28. V. Jurcak, D. Tsuzuki and I. Dan, “10/20, 10/10, and 10/5 systems revisited: their validity as relative head-surface-based positioning systems,” Neuroimage, 34 (4), 1600 –1611 (2007). https://doi.org/10.1016/j.neuroimage.2006.09.024 NEIMEF 1053-8119

29. J. Mourant and I. Bigio, Elastic- Scattering Spectroscopy and Diffuse Reflectance, , (2003).

30. S. K. Sekar et al., “Diffuse optical characterization of collagen absorption from 500 to 1700 nm,” J. Biomed. Opt., 22 (1), 015006 (2017). https://doi.org/10.1117/1.JBO.22.1.015006 JBOPFO 1083-3668

31. L. Kou, D. Labrie and P. Chylek, “Refractive indices of water and ice in the 0.65- to 2.5-μm spectral range,” Appl. Opt., 32 (19), 3531 –3540 (1993). https://doi.org/10.1364/AO.32.003531 APOPAI 0003-6935

32. S. L. Jacques and D. J. McAuliffe, “The melanosome: threshold temperature for explosive vaporization and internal absorption coefficient during pulsed laser irradiation,” Photochem. Photobiol., 53 (6), 769 –775 (1991). https://doi.org/10.1111/j.1751-1097.1991.tb09891.x PHCBAP 0031-8655

33. R. Nachabé et al., “Effect of bile absorption coefficients on the estimation of liver tissue optical properties and related implications in discriminating healthy and tumorous samples,” Biomed. Opt. Express, 2 (3), 600 –614 (2011). https://doi.org/10.1364/BOE.2.000600 BOEICL 2156-7085

34. N. Okui and E. Okada, “Wavelength dependence of crosstalk in dual-wavelength measurement of oxy- and deoxy-hemoglobin,” J. Biomed. Opt., 10 (1), 011015 (2005). https://doi.org/10.1117/1.1846076 JBOPFO 1083-3668

35. S. C. Gebhart, W. C. Lin and A. Mahadevan-Jansen, “In vitro determination of normal and neoplastic human brain tissue optical properties using inverse adding-doubling,” Phys. Med. Biol., 51 (8), 2011 –2027 (2006). https://doi.org/10.1088/0031-9155/51/8/004 PHMBA7 0031-9155

36. Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express, 17 (22), 20178 –20190 (2009). https://doi.org/10.1364/OE.17.020178 OPEXFF 1094-4087

37. A. Sassaroli and F. Martelli, “Equivalence of four Monte Carlo methods for photon migration in turbid media,” J. Opt. Soc. Am. A, 29 (10), 2110 –2117 (2012). https://doi.org/10.1364/JOSAA.29.002110 JOAOD6 0740-3232

38. E. Alerstam, T. Svensson and S. Andersson-Engels, “Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration,” J. Biomed. Opt., 13 (6), 060504 (2008). https://doi.org/10.1117/1.3041496 JBOPFO 1083-3668

39. K. B. Sung et al., “Accurate extraction of optical properties and top layer thickness of two-layered mucosal tissue phantoms from spatially resolved reflectance spectra,” J. Biomed. Opt., 19 (7), 077002 (2014). https://doi.org/10.1117/1.JBO.19.7.077002 JBOPFO 1083-3668

40. T.-X. Lin and K.-B. Sung, “Quantifying optical properties of multi-layered human head model from in-vivo spatially resolved near-infrared spectra (Conference Presentation),” Proc. SPIE, 10876 108760D (2019). https://doi.org/10.1117/12.2510963 PSISDG 0277-786X

41. C. J. Holmes et al., “Enhancement of MR images using registration for signal averaging,” J. Comput. Assist. Tomogr., 22 (2), 324 –333 (1998). https://doi.org/10.1097/00004728-199803000-00032 JCATD5 0363-8715

42. A. Bashkatov et al., “Optical properties of human cranial bone in the spectral range from 800 to 2000 nm,” Proc. SPIE, 6163 616310 (2006). https://doi.org/10.1117/12.697305 PSISDG 0277-786X

43. F. Bevilacqua et al., “In vivo local determination of tissue optical properties: applications to human brain,” Appl. Opt., 38 (22), 4939 –4950 (1999). https://doi.org/10.1364/AO.38.004939 APOPAI 0003-6935

44. M. Firbank et al., “Measurement of the optical properties of the skull in the wavelength range 650–950 nm,” Phys. Med. Biol., 38 (4), 503 –510 (1993). https://doi.org/10.1088/0031-9155/38/4/002 PHMBA7 0031-9155

45. C. R. Simpson et al., “Near-infrared optical properties of ex vivo human skin and subcutaneous tissues measured using the Monte Carlo inversion technique,” Phys. Med. Biol., 43 (9), 2465 –2478 (1998). https://doi.org/10.1088/0031-9155/43/9/003 PHMBA7 0031-9155

46. A. N. Bashkatov et al., “Optical properties of human skin, subcutaneous and mucous tissues in the wavelength range from 400 to 2000 nm,” J. Phys. D: Appl. Phys., 38 (15), 2543 –2555 (2005). https://doi.org/10.1088/0022-3727/38/15/004 JPAPBE 0022-3727

47. A. Bashkatov, E. Genina and V. Tuchin, “Optical properties of skin, subcutaneous, and muscle tissues: a review,” J. Innov. Opt. Health Sci., 4 (2011). https://doi.org/10.1142/S1793545811001319

48. P. van der Zee, M. Essenpreis and D. Delpy, “Optical properties of brain tissue,” Proc. SPIE, 1888 454 –465 (1993). https://doi.org/10.1117/12.154665 PSISDG 0277-786X

49. A. N. Yaroslavsky et al., “Optical properties of selected native and coagulated human brain tissues in vitro in the visible and near infrared spectral range,” Phys. Med. Biol., 47 (12), 2059 –2073 (2002). https://doi.org/10.1088/0031-9155/47/12/305 PHMBA7 0031-9155

50. E. K. Chan et al., “Effects of compression on soft tissue optical properties,” IEEE J. Sel. Top. Quantum Electron., 2 (4), 943 –950 (1996). https://doi.org/10.1109/2944.577320 IJSQEN 1077-260X

51. A. P. Tran, S. Yan and Q. Fang, “Improving model-based functional near-infrared spectroscopy analysis using mesh-based anatomical and light-transport models,” Neurophotonics, 7 (1), 015008 (2020). https://doi.org/10.1117/1.NPh.7.1.015008