2.1.

Geometry and Optical Properties of the Domain

The MC framework implements a methodology to produce a realistic tetrahedral-mesh heterogeneous domain of a section of the exposed cerebral cortex of a mouse (including pial vasculature and subpial brain tissue) from a 2-D grayscale image acquired in vivo using a conventional charge-coupled device. The workflow diagram describing this methodology is illustrated in Fig. 1.

Fig. 1

Workflow diagram of the methodology used in the Monte Carlo HSI framework to create a 3-D meshed domain of the exposed cortex: from an in vivo 2-D image (in grayscale), a binary mask is first created (in black and white) identifying the two media; then a 3-D mesh of the pial vasculature (in red) is generated, as well as a slab of subpial gray matter (in gray) encasing it; finally a 2-D source (in gold) and a 2-D detector (in green) are added to the final domain, with an additional mesh made of air (in cyan) filling the gap between the source and the cortex mesh.

The grayscale image of the exposed cortex, showing a $1.2\times 1.2\mathrm{mm}$ field of view (FOV) of the surface of the brain of a mouse and composed of $400\times 400\text{pixels}$, is first manually segmented to obtain a binary mask that differentiates between blood vessels and the surrounding brain tissue. A 3-D binary volume of the pial vasculature ($1.2\times 1.2\times 0.1\mathrm{mm}$) is then generated by expanding the mask along the vertical direction while symmetrically eroding the sections of the vessels from the central plane. This is done to replicate the curvature of the vascular geometry. The 3-D binary volume of the pial vasculature is then converted into a meshed volume using iso2mesh, constituting the first medium of the final domain. The pial vasculature volume is the encased in a $2.4\times 2.4\times 1\mathrm{mm}$ slab reproducing the surrounding mouse subpial gray matter. The extra layers added to the $1.2\times 1.2\mathrm{mm}$ FOV have the purpose of minimising boundary effects during the MC simulations.

Both media in the domain are defined by their geometry as well as by the associated optical properties (absorption coefficient ${\mu }_a$, scattering coefficient ${\mu }_s$, anisotropy $g$, and refractive index $n$). The medium that replicates the mouse subpial gray matter is considered to be made of water (${\mathrm{H}}_2\mathrm{O}$), lipid (fat), different concentrations of ${\mathrm{HbO}}_2$ and HHb (according to the fraction of blood and oxygen saturation level in the tissue), and different concentrations of the redox states of CCO, namely oxCCO and reduced CCO (redCCO). The medium reproducing both major and minor pial vessels (about 100 and $20\mu \mathrm{m}$ in diameter, respectively) includes water, fat, as well as ${\mathrm{HbO}}_2$ and HHb in different concentrations, according to the oxygen saturation value selected for the pial vasculature.

The composition and the optical properties of the two media are based on equations and reference data by Jacques.²⁰ Standard values, characteristic of general biological tissues, are assumed for the anisotropy and the refractive index of all the media of the domain, setting $g$ equal to 0.9 and $n$ equal to 1.365.²¹ The scattering coefficient ${\mu }_s(\lambda )$ is considered to be dependent only on the given wavelength $\lambda$ of the incident photon packet.²⁰ The absorption coefficient ${\mu }_a(\lambda )$ of each medium of the simulated domain is estimated as the sum of the single absorption coefficients, ${\mu }_{a,{\mathrm{H}}_2\mathrm{O}}(\lambda )$, ${\mu }_{a,\text{fat}}(\lambda )$, ${\mu }_{a,{\mathrm{HbO}}_2}(\lambda )$, ${\mu }_{a,\mathrm{HHb}}(\lambda )$, ${\mu }_{a,\mathrm{oxCCO}}(\lambda )$, and ${\mu }_{a,\mathrm{redCCO}}(\lambda )$, at the given wavelength $\lambda$, of the major chromophores composing the medium, i.e., water, fat, ${\mathrm{HbO}}_2$, HHb, oxCCO, and redCCO, respectively, and weighted accordingly to their content in it.²⁰ The data for ${\mu }_{a,\mathrm{H}2\mathrm{O}}(\lambda )$ and ${\mu }_{a,\text{fat}}(\lambda )$ in the NIR range are taken from Matcher et al.,²² for water, and van Veen et al.,²³ for fat (Table 4 in Appendix). The values of ${\mu }_{a,{\mathrm{HbO}}_2}(\lambda )$, ${\mu }_{a,\mathrm{HHb}}(\lambda )$, ${\mu }_{a,\mathrm{oxCCO}}(\lambda )$, and ${\mu }_{a,\mathrm{redCCO}}(\lambda )$ are calculated from the molar extinction coefficients ${\epsilon }_{{\mathrm{HbO}}_2}(\lambda )$, ${\epsilon }_{\mathrm{HHb}}(\lambda )$, ${\epsilon }_{\mathrm{oxCCO}}(\lambda )$, and ${\epsilon }_{\mathrm{redCCO}}(\lambda )$ of ${\mathrm{HbO}}_2$, HHb, oxCCO, and redCCO, respectively. In particular, for the oxCCO and redCCO contributions, this is done according to their selected concentrations [oxCCO] and [redCCO] in the given medium, whereas for the contributions of ${\mathrm{HbO}}_2$ and HHb, the average molar concentration of hemoglobin [Hb] in blood, the content $B$ of blood in the specific medium and the oxygen saturation $S$ are taken into account.²⁰ The molar extinction coefficients ${\epsilon }_{{\mathrm{HbO}}_2}(\lambda )$ and ${\epsilon }_{\mathrm{HHb}}(\lambda )$ of ${\mathrm{HbO}}_2$ and HHb are taken from Matcher et al.,²⁴ whereas the molar extinction coefficients ${\epsilon }_{\mathrm{oxCCO}}(\lambda )$ and ${\epsilon }_{\mathrm{redCCO}}(\lambda )$ of oxCCO and redCCO were measured by John Moody at the University of Plymouth in the bovine heart⁶ (Table 4 in Appendix).

The meshed domain of a section of mouse brain cortex is then integrated with a wide-field planar source for hyperspectral illumination at numerous wavelengths. The 2-D source has dimensions equal to $0.6\times 0.6\mathrm{mm}$ and is centered on the slab. It is also parallel to the top surface of the meshed domain, at a distance from it equal to 0.5 mm. The photon packets at each given wavelength $\lambda$ are launched from the surface of the planar source and evenly distributed over a $1.2\times 1.2\mathrm{mm}$ central section of the top surface of the domain, with a beam divergence of 90 deg. The MMC package implements the wide-field illumination source by mesh retessellation of the entire domain, creating an additional meshed medium between the source and the main domain, having the same optical properties of air [${\mu }_a(\lambda )$ and ${\mu }_s(\lambda )$ equal to $0{\mathrm{mm}}^{-1}$, and $g$ and $n$ equal to 1].^12,13,19

Finally, the Monte Carlo HSI framework also takes into account the detection and recording of information regarding the simulated photon packets by placing a $1.2\times 1.2\mathrm{mm}$ 2-D detector at the top surface of the mouse cortex domain, coextensive with the illumination field from the source. The choice of locating the detector precisely on the surface area of the domain has the advantage of maximizing the solid angle between the reflected photons and the detector, and thus the geometric detection efficiency of the configuration. This is not fully realistic, as it neglects the fraction of light that would be loss due to the distance between imaged target and detector (as well as the presence of the focusing optics), although such loss would only minimally affect the signal-to-noise ratio (SNR) of the results. Nonetheless, with this configuration, the MC framework does not have to take into account any lens or objective for focusing and collection of light in the simulations.

2.2.

Data Processing and Analysis

Hyperspectral illumination and imaging of the meshed domain representing the exposed brain cortex are reproduced using the MC framework by simulating photon incidence, diffusion, and reflection in each medium at different wavelengths in the NIR range, from 780 to 900 nm. At each execution of the MC code routine, 30 million ($3\times {10}^7$) photon packets are launched from the planar source, for each simulated wavelength. This number was chosen after performing convergence analysis. The photons reaching the detector surface after interacting with the domain are then recorded, in particular, the information about their final positions on the detector, their weights when they reached the detector and the partial pathlengths each of them have travelled in each medium. The detector is then divided into $185\times 185\text{pixels}$ ($6.5\text{-}\mu \mathrm{m}$ pixel size) and the detected photons for each wavelength are binned in these pixels according to their final position. The spatial images at each wavelength are then reconstructed by adding up the weights of all the photons binned in each pixel, in order to create a detected intensity map. A similar approach is used to reconstruct spatial maps of the average total photon pathlengths at each wavelength: these are obtained by summing up the partial pathlengths travelled in each medium by all the binned detected photons in each pixel, weighted by their corresponding weights, and then dividing this sum for the sum of the weights of the detected photons binned in that pixel. These maps provide the spatial distribution of the pathlength that a photon, arriving at a certain pixel, has travelled on average in the domain during a single run of the MC framework and for each wavelength. The reconstructed images at each wavelength are then stacked up to form 3-D spatiospectral datasets, called hyperspectral cubes or hypercubes. The same is done for the reconstructed spatial maps of the average total photon pathlengths to create 3-D average total photon pathlength distribution hypercubes.

For the computational studies reported here, two different brain physiological conditions are simulated, according to the different compositions of each medium of the mouse cortex model: (1) a baseline condition, representing the normal resting state of the brain and (2) an acute hypoxic condition, where cerebral oxygenation and metabolism drop significantly. Therefore, for each condition, the absorption properties of the media constituting the meshed domain of the exposed cortex are determined from their compositions. The scattering properties are only dependent on the selected wavelengths and thus are assumed constant between the two conditions. Water and fat contents are also assumed constant for each medium in both the two conditions. Furthermore, a significant decrease in oxygen saturation, as well as an increase in the total concentration of hemoglobin [to simulate an increase in cerebral blood volume (CBV)], are simulated in the pial vessels and in the subpial gray matter to recreate the hemodynamic response of the exposed cortex during the hypoxic conditions, leading to an overall decrease in the concentration of ${\mathrm{HbO}}_2$ and an increase in the concentration of HHb in the whole domain. Similarly, a reduction in the concentration of oxCCO and an increment in the concentration of redCCO are also applied only to the subpial gray matter medium, as to imitate the metabolic response to the lack of oxygen supply in the cerebral cortex. The concentration changes are selected so that the total sum of [oxCCO] and [redCCO] in the entire domain remains constant between the two conditions.⁶

For each simulated condition, image hypercubes and average total photon pathlength hypercubes are reconstructed. Light attenuation changes $\mathrm{\Delta }{A}_{k,l}$ ($\lambda$) between the simulated baseline and hypoxia are then calculated for each pixel $k$, $l$ (for $k$, $l=1\dots 185$) and each wavelength $\lambda$ from the photon intensities $I_{k,l}$ ($\lambda$) of the image hypercubes as

From Eq. (1), hemodynamic and metabolic maps charting the relative changes in concentrations $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, $\mathrm{\Delta }[\mathrm{HHb}]$, and $\mathrm{\Delta }[\mathrm{oxCCO}]$ of ${\mathrm{HbO}}_2$, HHb, and oxCCO, respectively, between the two conditions are estimated: this is done by applying the modified Beer–Lambert’s law (MBLL) to the simulated light attenuation changes $\mathrm{\Delta }{A}_{k,l}$ ($\lambda$), pixel by pixel.^6,25 Therefore, for each pixel $k$, $l$, the following system of algebraic equations is set as

3.1.

Study on HSI Performances and Accuracy (Study 1)

The first study on the performances and accuracy of HSI in reconstructing quantitative hemodynamic and metabolic maps of brain activity from the exposed cortex domain is conducted using the maximum allowable number of wavelengths in the NIR range from 780 to 900 nm, consisting of 121 wavelengths at 1-nm sampling, for both the baseline and the hypoxic condition. The compositions of the two media for both the simulated conditions are reported in Table 1,²⁰ from which the absorption properties used in the simulations are obtained.

Table 1

Different compositions of each medium in the meshed domain of the mouse brain cortex, for both the two simulated conditions (baseline and hypoxia).

At the onset of the acute hypoxic condition, an oxygen saturation drop $\mathrm{\Delta }S$ of $-35\%$ is mimicked in the pial vasculature and in the subpial gray matter, compared to the baseline.^29,30 Simultaneously, an increase of +30% in the total concentration [Hb] of hemoglobin in the pial vasculature and in the subpial gray matter is also simulated, as to replicate an overall increase in CBV during hypoxia.³¹ These two simulated phenomena correspond to a theoretical increase $\mathrm{\Delta }[\mathrm{HHb}]$ in the concentration of HHb of $+1162.79$ and $+43.60\mu \mathrm{M}$, in the pial vasculature and in the subpial gray matter, respectively, as well as to a theoretical decrease $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ in the concentration of ${\mathrm{HbO}}_2$ equal to $-465.12$ and $-17.44\mu \mathrm{M}$, in the pial vasculature and in the subpial gray matter, respectively. For the metabolic response, it is assumed that the relative concentration change $\mathrm{\Delta }[\mathrm{oxCCO}]$ of oxCCO in the subpial gray matter is equal to $-3\mu \mathrm{M}$ (this is mirrored by an equivalent increase in [redCCO]).^32,33

3.2.

Study on Optimal Selection of Wavelengths (Study 2)

For the second study, focused on evaluating the influence of the number and selection of NIR wavelengths on the quality and accuracy of the HSI data, the previous simulations for the two conditions (baseline and hypoxia) are repeated by changing the designated wavelengths for the illumination. Specifically, the following combinations of wavelengths are tested: (1) an arbitrary number of wavelengths in the range 780 to 900 nm, consisting of 25 wavelengths at 5-nm sampling and (2) an optimal selection of 8 wavelengths (784, 800, 818, 835, 851, 868, 881, and 894 nm) that was estimated by Arifler et al.³⁴ to be an ideal minimum combination of spectral bands for bNIRS to differentiate between the signals of hemoglobin and CCO with $<2\%$ mean error, compared to the “gold standard” of 121 wavelengths. The results of the two runs of the MC framework at different wavelengths are then compared with those of study 1, performed at the maximum allowable number of 121 wavelengths. This is intended to demonstrate that changing the number of wavelengths does not significantly affect the results of the quantification of the spatial changes in the concentrations of ${\mathrm{HbO}}_2$, HHb, and oxCCO (as long as the selected wavelengths are uniformly sampled in the NIR interval between 780 and 900 nm). Moreover, the outcomes of study 2 also aim at validating that the optimal selection of eight wavelengths for bNIRS is enough to obtain accurate results also in HSI of the exposed cortex with minimal differences from the results with 121 wavelengths.

3.3.

Assessment and Mitigation of Cross Talk and Partial Pathlength Effects (Study 3)

Evaluation of the presence and severity of cross talk and partial pathlength effects on the simulated hyperspectral data is conducted in the third study. This is done by re-running twice the simulations performed in study 1 while changing the optical properties of the hypoxic condition both times. Specifically: (1) first, simulations with the MC framework are run with only the metabolic response occurring (only the concentrations of redCCO and oxCCO change by $\pm 3\mu \mathrm{M}$, respectively), with no hemodynamic response (the saturation drop $\mathrm{\Delta }S$ and the increase in [Hb] are equal to zero in the whole domain, thus $[{\mathrm{HbO}}_2]$ and [HHb] do not change between the two conditions). (2) Second, the MC simulations are repeated this time with only the hemodynamic response occurring (the concentrations of ${\mathrm{HbO}}_2$ and HHb change according to the drop $\mathrm{\Delta }S$ in oxygen saturation equal to −35% and the increase in [Hb] equal to $+30\%$), whereas [oxCCO] and [redCCO] remain constant between the two conditions (no metabolic response is simulated). For both sets of simulations, the optimal combination of eight wavelengths tested in study 2 (784, 800, 818, 835, 851, 868, 881, and 894 nm) is selected, for each condition.

The new data from both runs of the MC framework are then compared with the results of study 1 to assess the influence of cross talk and partial pathlength effects in the reconstructed data, as well as to provide an indication of their potential sources. Simulating only the cerebral metabolic response during hypoxia, with no changes in the oxygenation of the tissues, though physiologically unrealistic and improbable, permits to isolate the single optical signature of CCO from those of hemoglobin, thus ideally limiting the occurrence of contamination effects from ${\mathrm{HbO}}_2$ and HHb to the minimum. Similarly, the simulation with only the brain hemodynamic response occurring should minimize any cross talk from CCO and related partial pathlength effects. Moreover, this approach can validate the simulated data in the realistic scenario from study 1 by demonstrating that the estimated changes in [oxCCO] are effectively obtained from true changes in the optical properties of the cerebral subpial tissue containing CCO between the two conditions, instead of arising from cross talk signals caused by changes in the concentrations of ${\mathrm{HbO}}_2$ and HHb or from the influence of the variance of the photon pathlengths.

3.4.

Alternative HSI Configuration (Study 4)

In the fourth and final study, the Monte Carlo HSI framework is used to explore the implementation of a more localized and selective hyperspectral illumination and detection configuration, designed to improve the accuracy of the quantification of the hemodynamic and metabolic responses in the subpial gray matter, as well as to further mitigate cross talk effects with hemoglobin and partial pathlength effects. In particular, this configuration consists in reducing the illumination area and the FOV of the 2-D detector from $1.2\times 1.2$ to $0.2\times 0.2\mathrm{mm}$ (using the same number of pixels, $185\times 185$, in the reconstruction). This is obtained by decreasing the dimension of the source area from $0.6\times 0.6$ to $0.1\times 0.1\mathrm{mm}$ and then moving its center to align it to the new detector FOV, as depicted in Fig. 2(a). Such configuration enables to selectively illuminate only a portion of the domain outside the vasculature [Fig. 2(b)], which contains only subpial gray matter, as well as to collect only information from photon packets arriving in the same region. Simulations with the MC framework are run again using the same optical properties used in the study 1 (from Table 1) and with the optimal combinations of eight wavelengths (784, 800, 818, 835, 851, 868, 881, and 894 nm) used in both study 2 and study 3.

Fig. 2

(a) New meshed domain implementing a 2-D source (in gold) and detector (in green) producing a localized $0.2\times 0.2\mathrm{mm}$ illumination and FOV. (b) Position of the localized $0.2\times 0.2\mathrm{mm}$ FOV of the detector (in green) on the simulated domain, compared to the $1.2\times 1.2\mathrm{mm}$ illumination field and detection FOV used in the previous studies.

4.1.

Study 1

Figure 3 depicts the two hemodynamic maps, for $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HHb}]$, and the metabolic map of $\mathrm{\Delta }[\mathrm{oxCCO}]$ tracking the relative changes in concentration of the three targeted chromophores during the acute hypoxic condition that was simulated in study 1, using 121 wavelengths between 780 and 900 nm. The hemodynamic maps of the relative changes in concentration of ${\mathrm{HbO}}_2$ [Fig. 3(b)] and HHb [Fig. 3(c)] present high image quality and spatial resolution, compared to the actual depiction of the FOV of the simulated domain [Fig. 3(a)]. The large vascular hemodynamic response related to both chromophores is accurately localized within the boundaries of the pial vasculature, resolving both major (about $100\mu \mathrm{m}$ in diameter) and minor vessels (about $20\mu \mathrm{m}$ in diameter), as well as showing a decrease in the concentration of ${\mathrm{HbO}}_2$ and an increase in the concentration of HHb, as theoretically expected. Similarly, a minor hemodynamic response from ${\mathrm{HbO}}_2$ and HHb is also reconstructed in the surrounding tissue that is consistent with the simulate changes in oxygen saturation and blood volume in the subpial gray matter. However, from the hemodynamic maps, a large underestimation in the quantification of both $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HHb}]$ in the pial vasculature clearly emerges. The metabolic map of the relative changes in concentration of oxCCO [Fig. 3(d)] shows a poorer image quality than the hemodynamic maps, due to lower SNR in the processed data for CCO and the presence of spurious measured changes in concentration of CCO in the pial vasculature. These factors make difficult to fully localize the metabolic response and to differentiate between pial vasculature and surrounding tissue with high spatial resolution. Only the major pial vessels (about $100\mu \mathrm{m}$ in diameter) are partially resolved in the metabolic map.

Fig. 3

(a) Picture of the $1.2\times 1.2\mathrm{mm}$ FOV of the detector on the simulated domain, showing the position of two $65\times 65\mu \mathrm{m}$ ROIs used in the data analysis, one including only pial vasculature (blue square) and the other only subpial gray matter (black square). (b) Hemodynamic map charting the relative changes $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ in the concentration of ${\mathrm{HbO}}_2$ between baseline and hypoxia. (c) Hemodynamic map showing the relative changes $\mathrm{\Delta }[\mathrm{HHb}]$ in the concentration of HHb between baseline and hypoxia. (d) Metabolic map showing the relative changes $\mathrm{\Delta }[\mathrm{oxCCO}]$ in the concentration of oxCCO between baseline and hypoxia. (e) New hemodynamic map of the relative changes $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ in the concentration of ${\mathrm{HbO}}_2$ between baseline and hypoxia, after postprocessing correction. (f) New hemodynamic map of the relative changes $\mathrm{\Delta }[\mathrm{HHb}]$ in the concentration of HHb between baseline and hypoxia, after postprocessing correction. (g) New metabolic map of the relative changes $\mathrm{\Delta }[\mathrm{oxCCO}]$ in the concentration of oxCCO between baseline and hypoxia, after postprocessing correction.

Evaluation of the accuracy in quantifying the correct relative changes in the concentrations of ${\mathrm{HbO}}_2$, HHb, and oxCCO is performed by calculating and analyzing the spatial averages of the concentration changes $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, $\mathrm{\Delta }[\mathrm{HHb}]$, and $\mathrm{\Delta }[\mathrm{oxCCO}]$ in specific regions of interest (ROIs) in the hemodynamic and metabolic maps. Two ROIs of $10\times 10\text{pixels}$, both corresponding to square regions of $65\times 65\mu \mathrm{m}$ in the FOV, are selected: (1) one including only pial vasculature and (2) one including only subpial gray matter. The position and size of the two ROIs on the FOV of the detector are shown in Fig. 3(a). The concentrations changes for each chromophore are spatially averaged across the pixels of each ROI. The values of the averages $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$, $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$, and $⟨\mathrm{\Delta }[\mathrm{oxCCO}]⟩$ in the two ROIs are reported in Table 2 and compared with the corresponding theoretical values. The values of the average concentration changes $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$ and $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$ in the ROI associated with the pial vessels reproduce the trend of the expected temporal hemodynamic response to the simulated lack of oxygenation to brain tissue, although they are significantly lower ($-32.89\mu \mathrm{M}$ for ${\mathrm{HbO}}_2$ and $82.78\mu \mathrm{M}$ for HHb) than the theoretical simulated changes ($-465.12\mu \mathrm{M}$ for ${\mathrm{HbO}}_2$ and $1162.79\mu \mathrm{M}$ for HHb). The corresponding quantification error is equal to about 92.9%. In addition, an erroneous decrease in $\mathrm{\Delta }[\mathrm{oxCCO}]$ is also estimated for the same ROI in the pial vasculature ($-2.929\mu \mathrm{M}$), which is in the same order of magnitude of the actual simulated change in the concentration of oxCCO ($-3\mu \mathrm{M}$). Finally, the quantification of the relative changes $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, $\mathrm{\Delta }[\mathrm{HHb}]$, and $\mathrm{\Delta }[\mathrm{oxCCO}]$ in the concentration of ${\mathrm{HbO}}_2$, HHb, and oxCCO in the subpial gray matter shows smaller estimation errors of about 10.4%, 11.3%, and 5.2%, respectively, as demonstrated by the value of the average concentration changes $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$, $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$, and $⟨\mathrm{\Delta }[\mathrm{oxCCO}]⟩$ in the central ROI ($-19.32\mu \mathrm{M}$ for ${\mathrm{HbO}}_2$, $48.54\mu \mathrm{M}$ for HHb, and $-3.156\mu \mathrm{M}$ for oxCCO), which are all close to the theoretical changes in the simulated chromophores ($-17.44\mu \mathrm{M}$ for ${\mathrm{HbO}}_2$, $43.60\mu \mathrm{M}$ for HHb, and $-3\mu \mathrm{M}$ for oxCCO).

Table 2

Comparison between the spatial average changes ⟨Δ[HbO2]⟩, ⟨Δ[HHb]⟩, and ⟨Δ[oxCCO]⟩ in the concentrations of HbO2, HHb, and oxCCO, respectively, and the corresponding theoretical values in two different ROIs of 10×10  pixels on the reconstructed hemodynamic and metabolic maps: (1) for the results of study 1, at 121 NIR wavelengths, both before (A) and after correction (B); and (2) for the results of study 2, at 25 NIR wavelengths (C) and at 8 optimal NIR wavelengths (D), both after correction.

Both the large underestimation in the changes of concentration of vascular ${\mathrm{HbO}}_2$ and HHb in the hemodynamic maps, as well as the occurrence of spurious measured changes in the concentration of oxCCO in the pial vasculature, could be connected to partial pathlength effects. The latter should not be confused with cross talk, because the erroneous measured values of $\mathrm{\Delta }[\mathrm{oxCCO}]$ are not induced by a genuine change in the concentration of this chromophore (since it is not present in the ROI), but they are due to the significant difference between the partial pathlengths of the detected photons that travelled in the pial vasculature and those that travelled in the surrounding subpial brain tissue. This is further investigated and validated by the results of study 3.

Figure 4 shows examples, at 835 nm, of the average total pathlength maps of the detected photons across the entire domain [Fig. 4(a)], as well as the average partial pathlength maps of the detected photons in the pial vasculature [Fig. 4(b)] and in the subpial gray matter [Fig. 4(c)], respectively. The maps compare the fractions of the average total pathlengths travelled by the detected photons in each of the two media, during the baseline condition. It can be seen that the partial pathlengths of the detected photons in the pial vasculature are considerably shorter than the partial pathlengths the same photons travelled in the subpial gray matter. The latter also account for more than 97% of the average total pathlength. Moreover, the comparison between the average partial pathlength maps reveals that the majority of photons that were detected in pixels located on the pial vasculature have effectively travelled mostly in the subpial gray matter. This could explain both the significant underestimation of $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HHb}]$ in the pial vessels, as well as the occurrence of the spurious measured changes $\mathrm{\Delta }[\mathrm{oxCCO}]$ in the same vascular medium, resulting from applying MBLL from Eq. (2) and using the average total pathlength of the detected photons.

Fig. 4

(a) Average total photon pathlength map at 835 nm. (b) Average partial photon pathlength map in the pial vasculature at 835 nm. (c) Average partial photon pathlength map in the subpial gray matter at 835. (d) Map of the correction factors ${\mathrm{CF}}_{k,l}^{\prime }$ obtained from the mean ratios between the average total pathlengths of the detected photons and the average partial pathlengths of the same photons in the pial vasculature, across all the wavelengths and between both simulated conditions. (e) Map of the correction factors ${\mathrm{CF}}_{k,l}^{″}$ obtained from the mean ratios between the average total pathlengths of the detected photons and the average partial pathlengths of the same photons in the subpial gray matter, across all the wavelengths and between both simulated conditions.

A postprocessing correction of the hemodynamic and metabolic maps using the information about the average partial photon pathlengths is here proposed, to primarily improve the quantification of the changes of concentrations of ${\mathrm{HbO}}_2$ and HHb in the pial vasculature. Two maps of correction factors, ${\mathrm{CF}}_{k,l}^{\prime }$ and ${\mathrm{CF}}_{k,l}^{″}$ (for each pixel $k$, $l$), are produced: (1) ${\mathrm{CF}}_{k,l}^{\prime }$ are the means across all the selected $M$ wavelengths (in this case $M=121$) of the ratios between the average total pathlengths ${\mathrm{PL}}_{k,l}(\lambda )$ travelled by the detected photons and the average partial pathlengths ${\mathrm{PL}}_{k,l,\text{vessel}}(\lambda )$ they travelled in the pial vasculature (for each wavelength $\lambda$). (2) ${\mathrm{CF}}_{k,l}^{″}$ are the means across all the selected $M$ wavelengths of the ratios between the average total pathlengths ${\mathrm{PL}}_{k,l}(\lambda )$ travelled by the detected photons and the average partial pathlengths ${\mathrm{PL}}_{k,l,\text{gray}}(\lambda )$ they travelled in the subpial gray matter (for each wavelength $\lambda$). Thus

The two sets of correction maps can be obtained using Eq. (3) for both the baseline and the hypoxic condition, respectively. The two final sets of correction maps [Fig. 4(d) for ${\mathrm{CF}}_{k,l}^{\prime }$ and Fig. 4(e) for ${\mathrm{CF}}_{k,l}^{″}$] to apply to the hemodynamic and metabolic maps are calculated from the values of the mean between the corresponding correction factors of the two conditions (for each pixel $k$, $l$), respectively.

The postprocessing correction of the hemodynamic and metabolic maps via the two correction maps in Figs. 4(d) and 4(e) is performed selectively using the segmented binary map of the FOV (utilised during the mesh domain creation, as shown in Fig. 1) as a guide. Thus for pixels $k$, $l$ corresponding to the pial vasculature medium in the binary mask, the following correction is applied to obtain the corrected values $\mathrm{\Delta }[{\mathrm{HbO}}_2]*$ and $\mathrm{\Delta }[\mathrm{HHb}]*$ of the changes in concentration of ${\mathrm{HbO}}_2$ and HHb in the hemodynamic maps only:

The correction in Eq. (4) is not applied to the metabolic map since theoretically no CCO in the pial vasculature is simulated. Contrarily, for pixels $k$, $l$ corresponding to the subpial gray matter medium in the binary mask, this other correction is applied to both the hemodynamic and the metabolic maps to obtain the corrected values $\mathrm{\Delta }[{\mathrm{HbO}}_2]*$, $\mathrm{\Delta }[\mathrm{HHb}]*$, and $\mathrm{\Delta }[\mathrm{oxCCO}]*$ of the changes in concentration of ${\mathrm{HbO}}_2$, HHb. and oxCCO:

In both sets of Eqs. (4) and (5), ${\mathrm{CF}}_{k,l}^{\prime }$ and ${\mathrm{CF}}_{k,l}^{″}$ correspond to the two sets of correction factors (for each pixel $k$, $l$) calculated from Eq. (3). This selective correction permits one to weight the hyperspectral data by taking into account the large differences in the average partial photon pathlengths between the two media of the domain.

The new hemodynamic and metabolic maps resulting from the postprocessing correction with Eqs. (3)–(5) are depicted in Figs. 3(e)–3(g) for ${\mathrm{HbO}}_2$, HHb, and oxCCO, respectively. A large enhancement in image contrast, as well as in the accuracy in the localization of the hemodynamic response in the pial vasculature, is evident from the new hemodynamic maps, due to the significant improvement in the quantification of $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HHb}]$ in the pial vessels. Contrarily, the correction in the subpial gray matter produces minimal effects in the hemodynamic and metabolic maps. This is due to the similarity between the average total pathlength and the average partial pathlength of the detected photons in the subpial gray matter, as visible by comparing Figs. 4(a) and 4(c).

The efficacy of the postprocessing correction is further corroborated by the values of the spatial averages $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$, $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$, and $⟨\mathrm{\Delta }[\mathrm{oxCCO}]⟩$ of the concentration changes of ${\mathrm{HbO}}_2$, HHb, and oxCCO in the same two ROIs used for the uncorrected maps. Table 2 shows that the quantified values of $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$ and $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$ in the pial vasculature are now much closer to the simulated theoretical values ($-474.7\mu \mathrm{M}$ for $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$ and $1192.2\mu \mathrm{M}$ for $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$), with estimation errors of about 2.06% and 2.87%, for ${\mathrm{HbO}}_2$ and HHb, respectively. However, the analysis in the ROI localized on the subpial gray matter demonstrates negligible differences ($<1\%$) in the estimates of $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$, $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$, and $⟨\mathrm{\Delta }[\mathrm{oxCCO}]⟩$ between the corrected and the uncorrected maps, for all the three chromophores. This suggests that the postprocessing correction is only necessary for the pial vasculature in the hemodynamic maps.

A comparison between the cross-section views of the theoretical values of $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, $\mathrm{\Delta }[\mathrm{HHb}]$, and $\mathrm{\Delta }[\mathrm{oxCCO}]$ and the corresponding reconstructed values, both before correction and after postprocessing correction, is provided in Fig. 5. This analysis on a line of pixels offers additional insight on the partial pathlength effects in all the three maps: in the metabolic map, a large variance characterises the spurious estimated changes in the concentration of oxCCO in the pial vasculature. Figure 5 also further highlights how the postprocessing correction produces: (1) a considerable improvement in spatial localization of the hemodynamic response; (2) a substantial enhancement in the accuracy of the quantification of the relative changes in concentrations of ${\mathrm{HbO}}_2$ and HHb; and (3) insignificant differences in the quantification of the relative changes in concentrations of all the three chromophores in the subpial gray matter, in both the hemodynamic and metabolic maps. This last aspect can be clearly seen in Fig. 5(d), where the values of $\mathrm{\Delta }[\mathrm{oxCCO}]$ before and after correction are almost overlapping, as well as for the values of $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HHb}]$ in the pixels on the subpial gray matter, in both Figs. 5(b) and 5(c).

Fig. 5

(a) Position on the FOV of the detector of the line of pixel (in blue) used in the data analysis. (b) Relative changes $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ in the concentration of ${\mathrm{HbO}}_2$ along the line of pixels. (c) Relative changes $\mathrm{\Delta }[\mathrm{HHb}]$ in the concentration of HHb along the line of pixels. (d) Relative changes $\mathrm{\Delta }[\mathrm{oxCCO}]$ in the concentration of oxCCO along the line of pixels. Values are depicted both before and after the correction.

4.2.

Study 2

In study 2, similar hemodynamic and metabolic maps for $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, $\mathrm{\Delta }[\mathrm{HHb}]$ and $\mathrm{\Delta }[\mathrm{oxCCO}]$ are reproduced: (1) first using an arbitrary number of 25 NIR wavelengths between 780 and 900 nm at 5-nm sampling and then (2) using an optimal selection of eight NIR wavelengths (784, 800, 818, 835, 851, 868, 881, and 894 nm),³⁴ for the same two simulated brain conditions (baseline and hypoxia). The same postprocessing correction of the hemodynamic and metabolic maps from study 1 is also applied, using Eqs. (3)–(5). Calculation of the spatial averages $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$, $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$, and $⟨\mathrm{\Delta }[\mathrm{oxCCO}]⟩$ of the relative changes in the concentrations of ${\mathrm{HbO}}_2$, HHb, and oxCCO is also performed in the same two ROIs used in study 1, for both the two sets of new maps at 25 and 8 wavelengths. All these values are shown in Table 2: they differ marginally from the corresponding ones in study 1, relative to the data at 121 wavelengths. The differences in the corresponding estimates of $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$, $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$, and $⟨\mathrm{\Delta }[\mathrm{oxCCO}]⟩$ in both ROIs, between the three combinations of selected wavelengths, varies from 0% to a maximum of 2.1%. In particular, accuracy in quantifying the relative changes in the concentration of ${\mathrm{HbO}}_2$, HHb, and oxCCO do not seem to be significantly affected by reducing the spectral information from the maximum allowable number of 121 wavelengths in the selected NIR range down to an optimal combination of only 8. Therefore, these results corroborate the findings of Arifler et al.³⁴ that were estimated for bNIRS, extending them also to HSI targeting brain metabolism.

4.3.

Study 3

The results of the first run of the MC framework in study 3, where only the metabolic response is simulated during the hypoxic condition with no changes in hemoglobin in the domain, provide an insight on the influence of cross talk and partial pathlength effects on the reconstruction of the hyperspectral data. The set of hemodynamic and metabolic maps calculated from the hypercubes simulated in this scenario, using the optimal combination of eight wavelengths (784, 800, 818, 835, 851, 868, 881, and 894 nm) from study 2, are shown in the top row of Fig. 6, for $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, $\mathrm{\Delta }[\mathrm{HHb}]$, and $\mathrm{\Delta }[\mathrm{oxCCO}]$. No postprocessing correction was performed in this case.

Fig. 6

Top row: Hemodynamic and metabolic maps of the relative changes (a) $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, (b) $\mathrm{\Delta }[\mathrm{HHb}]$, and (c) $\mathrm{\Delta }[\mathrm{oxCCO}]$ in the concentration of ${\mathrm{HbO}}_2$, HHb, and oxCCO between baseline and a hypoxic condition, during which only the metabolic response occurs (no simulated changes in hemoglobin). Bottom row: hemodynamic and metabolic maps of the relative changes (d) $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, (e) $\mathrm{\Delta }[\mathrm{HHb}]$, and (f) $\mathrm{\Delta }[\mathrm{oxCCO}]$ in the concentration of ${\mathrm{HbO}}_2$, HHb, and oxCCO between baseline and a hypoxic condition, during which only the hemodynamic response is simulated (no changes in CCO occur).

In the absence of any replicated hemodynamic response in the simulations, the hemodynamic maps do not display any contrast provided by the changes $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HHb}]$ in hemoglobin, neither in the pial vasculature nor in the surrounding tissue, as expected. Therefore, no cross talk effect from CCO is found. The image quality of the metabolic map is higher compared to the results from study 1, due to the lower influence of the optical signatures of hemoglobin in the data. However, nonzero relative changes in the concentration of oxCCO are still estimated in the vessels, even though the pial vasculature does not contain any CCO.

Further validation to these deductions is obtained by looking at the spatial averages $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$, $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$, and $⟨\mathrm{\Delta }[\mathrm{oxCCO}]⟩$ of the relative changes in the concentrations of ${\mathrm{HbO}}_2$, HHb, and oxCCO, again calculated for the same two ROIs analyzed in the previous studies. As reported in Table 3, the values of the averages $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$ and $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$ in both ROIs are close to zero. Larger and still non-negligible spurious measurements are estimated for oxCCO from the analysis of the spatial averages $⟨\mathrm{\Delta }[\mathrm{oxCCO}]⟩$ in the ROI corresponding to the pial vasculature ($-3.018\mu \mathrm{M}$). The magnitude of these is still in the same order of the simulated relative change in concentration of oxCCO ($-3\mu \mathrm{M}$). Finally, a more accurate quantification of the relative change $\mathrm{\Delta }[\mathrm{oxCCO}]$ in the metabolic map emerges from the calculation of the spatial average $⟨\mathrm{\Delta }[\mathrm{oxCCO}]⟩$ in the ROI corresponding to the gray matter ($-3.098\mu \mathrm{M}$), which appears closer to the effective simulated change in oxCCO than the result from study 1 ($-3.157\mu \mathrm{M}$). The corresponding estimation error is about 3.27% (against 5.2% in study 1).

Table 3

Spatial average changes ⟨Δ[HbO2]⟩, ⟨Δ[HHb]⟩, and ⟨Δ[oxCCO]⟩ in the concentrations of HbO2, HHb, and oxCCO, respectively, in two different ROIs of the reconstructed maps in study 3 (3), obtained for: (A) the simulations with only metabolic response during hypoxia (no correction applied) and (B) the simulations with only hemodynamic response during hypoxia (no correction applied). Both datasets were simulated using eight optimal NIR wavelengths between 780 and 900 nm.

Additional insight on the phenomenon of partial pathlength effect is inferred from the results of the second run of the MC framework in this study: contrarily to the previous run, only the hemodynamic response during the hypoxic condition is simulated in this case, in both the pial vasculature and the subpial gray matter, whereas no change in CCO occurs in the subpial gray matter. Again, the simulations are run for the same optimal combination of eight wavelengths and no postprocessing correction is applied to the three maps.

The bottom row of Fig. 6 illustrates the hemodynamic and metabolic maps of $\mathrm{\Delta }[{\mathrm{HbO}}_2]$, $\mathrm{\Delta }[\mathrm{HHb}]$, and $\mathrm{\Delta }[\mathrm{oxCCO}]$ obtained from the simulations with only changes in hemoglobin. As expected, the hemodynamic maps, when only changes in hemoglobin occur, again measure and localize the simulated hemodynamic response in both the pial vasculature and the subpial gray matter, similar to the results obtained for study 1 with 121 wavelengths (before postprocessing correction), as seen in Figs. 3(a) and 3(b). The changes $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ and $\mathrm{\Delta }[\mathrm{HHb}]$ in the pial vessels are still greatly underestimated, suggesting that the underestimation is not affected by the presence of the metabolic response of CCO and thus is not generated by cross talk. Furthermore, no contrast between pial vasculature and subpial gray matter appears in the metabolic map.

The analysis of the spatial averages $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$, $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$, and $⟨\mathrm{\Delta }[\mathrm{oxCCO}]⟩$ in the two selected ROIs (Table 3) clearly demonstrates the extent of the relative changes in the concentrations of oxCCO still present in the metabolic map is very minimal ($-0.086\mu \mathrm{M}$ on average in the ROI including only the subpial gray matter), as well as any occurrence of spurious measured changes in the concentration of oxCCO in the pial vessels ($⟨\mathrm{\Delta }[\mathrm{oxCCO}]⟩$ equal to $-0.083\mu \mathrm{M}$ in the pial vasculature), compared to the results in study 1 (Table 2).

The findings in study 3 further validate the assumption that the spurious measured signals in a region of the maps are not affected by the presence of the concentration change of another chromophore (cross talk between the chromophores) but purely arise from partial pathlength effects.

4.4.

Study 4

The new HSI configuration tested in the fourth and final study with the MC framework, implementing and simulating a $0.2\times 0.2\mathrm{mm}$ illumination field and detection FOV, explores the possibility to improve accuracy in quantifying brain hemodynamic and metabolic response in the subpial gray matter during the hypoxic condition, compared to the earlier results obtained in study 1 and study 2, without the need of postprocessing correction. The MC framework is run using the optimal combination of eight wavelengths (784, 800, 818, 835, 851, 868, 881, and 894 nm). The reconstruction is performed in the same way as for the previous studies, providing hemodynamic and metabolic maps composed of $185\times 185\text{pixels}$. Since the FOV is now smaller ($0.2\times 0.2\mathrm{mm}$, against the previous $1.2\times 1.2\mathrm{mm}$ FOV), the size of the pixels decreases from 6.5 to $1.08\mu \mathrm{m}$. No postprocessing correction is applied to the reconstructed hyperspectral data obtained with the new simulated HSI configuration.

The reconstructed maps for the localized illumination and imaging are not reported: this is because they do not show any significant spatial contrast since the new configuration involves the FOV being located entirely over a homogeneous area of subpial gray matter, with no features to be differentiated. Spatial averages $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$, $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$, and $⟨\mathrm{\Delta }[\mathrm{oxCCO}]⟩$ of the relative changes in the concentrations of ${\mathrm{HbO}}_2$, HHb, and oxCCO are calculated from the hemodynamic and metabolic maps for an ROI of $60\times 60\text{pixels}$ concentric with the $0.2\times 0.2\text{-}\mathrm{mm}$ FOV. The ROI corresponds to a square region of about $65\times 65\mu \mathrm{m}$ of subpial gray matter. This is done to conduct the spatial average analysis on exactly the same portion of subpial gray matter that was targeted in all the previous studies.

An improvement in the quantification of the concentrations of both ${\mathrm{HbO}}_2$ and HHb in the subpial gray matter is achieved with the new configuration, without postprocessing correction, compared to the corresponding values obtained in study 1 for the same ROI (Table 2). The spatial averages $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$ and $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$ in the ROI for the new HSI configuration stand at $-17.67\mu \mathrm{M}$ for $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$ and $44.64\mu \mathrm{M}$ for $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$, against $-19.39$ and $48.72\mu \mathrm{M}$, respectively for study 1. The new values are much closer to the theoretical simulated changes in the concentration of ${\mathrm{HbO}}_2$ and HHb in the subpial gray matter ($-17.44\mu \mathrm{M}$ for ${\mathrm{HbO}}_2$ and $43.60\mu \mathrm{M}$ for HHb). The quantification errors for $⟨\mathrm{\Delta }[{\mathrm{HbO}}_2]⟩$ and $⟨\mathrm{\Delta }[\mathrm{HHb}]⟩$ with the new configuration decrease to about 0.75% and 2.11%, respectively, against 10.4% and 11.3% for study 1. This is due to targeting a smaller volume of cerebral tissue, thus reducing the influence of scattering on the estimated average photon pathlengths, as well as to avoiding the illumination of the pial vasculature, which significantly reduces the possibility that a photon may have travelled through that region.

Finally, the quantification of the relative changes in the concentration of oxCCO achieved with the alternative hyperspectral configuration is also more accurate than the one obtained in study 1 ($-3.121\mu \mathrm{M}$ for $⟨\mathrm{\Delta }[\mathrm{oxCCO}]⟩$ with the $0.2\times 0.2\mathrm{mm}$ FOV, compared to $-3.156\mu \mathrm{M}$ with the $1.2\times 1.2\mathrm{mm}$ FOV) and thus closer to the simulated metabolic response ($-3\mu \mathrm{M}$). The estimation error of $⟨\mathrm{\Delta }[\mathrm{oxCCO}]⟩$ with the new HSI configuration is about 4.03% (compared to 5.2% with the larger FOV used in study 1). This is the most accurate quantification of the metabolic response from oxCCO obtained among all the reported studies (excluding the unrealistic cases of study 3).