Investigation of the quantification of hemoglobin and cytochrome-c-oxidase in the exposed cortex with near-infrared hyperspectral imaging: a simulation study

Abstract. Significance: We present a Monte Carlo (MC) computational framework that simulates near-infrared (NIR) hyperspectral imaging (HSI) aimed at assisting quantification of the in vivo hemodynamic and metabolic states of the exposed cerebral cortex in small animal experiments. This can be done by targeting the NIR spectral signatures of oxygenated (HbO2) and deoxygenated (HHb) hemoglobin for hemodynamics as well as the oxidative state of cytochrome-c-oxidase (oxCCO) for measuring tissue metabolism. Aim: The aim of this work is to investigate the performances of HSI for this specific application as well as to assess key factors for the future design and operation of a benchtop system. Approach: The MC framework, based on Mesh-based Monte Carlo (MMC), reproduces a section of the exposed cortex of a mouse from an in vivo image and replicates hyperspectral illumination and detection at multiple NIR wavelengths (up to 121). Results: The results demonstrate: (1) the fitness of the MC framework to correctly simulate hyperspectral data acquisition; (2) the capability of HSI to reconstruct spatial changes in the concentrations of HbO2, HHb, and oxCCO during a simulated hypoxic condition; (3) that eight optimally selected wavelengths between 780 and 900 nm provide minimal differences in the accuracy of the hyperspectral results, compared to the “gold standard” of 121 wavelengths; and (4) the possibility to mitigate partial pathlength effects in the reconstructed data and to enhance quantification of the hemodynamic and metabolic responses. Conclusions: The MC framework is proved to be a flexible and useful tool for simulating HSI also for different applications and targets.


Introduction
Hyperspectral imaging (HSI) is an emerging optical technique for biomedical applications that can be potentially used to quantitatively monitor in vivo changes in the metabolic and hemodynamic states of the brain, specifically on the exposed cortex. HSI provides extensive spectral information, in addition to spatial data, by acquiring images over a broad range of the light spectrum at numerous and contiguous wavelength bands. 1,2 Changes in the concentrations of relevant biomarkers, such as oxyhemoglobin (HbO 2 ) and deoxyhemoglobin (HHb), can be retrieved by measuring the intensity changes of multiple different wavelengths of reflected light after having interacted with the cerebral tissue. These light intensity changes originate from variations in the optical properties of brain tissue during physiological processes, e.g., changes *Address all correspondence to Luca Giannoni, E-mail: l.giannoni@ucl.ac.uk in brain oxygenation and perfusion. 3,4 HSI can also be used to target the changes of the redox state of cytochrome-c-oxidase (CCO), which is a chromophore directly involved in the production of adenosine triphosphate in the mitochondria. 5 CCO has a high specificity as a biomarker for monitoring brain metabolism, due to its high concentration in the cortical tissue. 5,6 Metabolic monitoring through CCO is primarily performed noninvasively via broadband near-infrared spectroscopy (bNIRS), which, similar to HSI, analyzes spectroscopically a large number of wavelengths (tens to hundreds) in the near-infrared (NIR) range between 780 and 900 nm. 6 This specific range is chosen due to the presence of a predominant broad peak in the absorption spectrum of the copper CuA redox center of CCO, which enables a better differentiation of the CCO signal from those of HbO 2 and HHb. 6,7 However, bNIRS is not a wide-field imaging technique and it is limited in terms of spatial resolution, due to the high-scattering properties of biological tissue in the NIR range. Furthermore, bNIRS only provides information about changes in metabolism and hemodynamics that are averaged over relatively large volumes of tissue (typically from 1 to 100 cm 3 ), which can include both blood vessels as well as surrounding extravascular tissue, since it is based on measuring NIR light diffusing through the scalp, the skull, and the gray and white matter. For this reason, looking directly at the exposed cerebral cortex using HSI in the NIR range could provide additional and more exhaustive information about brain metabolism and hemodynamics, in particular by spatially differentiating between the regions where the two processes are primarily located, i.e., the pial vasculature and the surrounding subpial cortical tissue, for the hemodynamic response and the metabolic response, respectively. The NIR hyperspectral approach targeting the exposed cortex can thus be used to obtain a deeper understanding of brain physiology during different conditions, such as hypoxic-ischemia or other similar abnormal alterations in brain oxygenation.
We present here a Monte Carlo (MC) framework that simulates NIR HSI of the hemodynamic and metabolic states of the exposed cortex. To our knowledge, no MC computational analysis has been published before to reproduce wide-field HSI simultaneously targeting the changes in HbO 2 , HHb, and the oxidative state of CCO (oxCCO). This approach is used to investigate the feasibility and performances of using HSI in the NIR range to quantitatively measure changes in concentrations of the abovementioned biomarkers, by simulating a realistic portion of mouse brain cortex (created from an in vivo image) during changes from cerebral normoxia to acute hypoxia. In particular, the computational analysis focuses on: (1) assessing the capacity of HSI to reconstruct spatial maps of metabolic and hemodynamic activity; (2) evaluating the accuracy of HSI in quantitatively estimating relative changes in the concentrations of HbO 2 , HHb and oxCCO; (3) investigating what is the optimal selection and number of wavelength bands to use for HSI to simultaneously image HbO 2 , HHb, and oxCCO; and (4) studying the effects, influence, and magnitude of cross talk and partial pathlength effects affecting the hemoglobin and oxCCO signals. We define cross talk as the erroneous measured change in the concentration of a chromophore that is induced by the genuine concentration change of another chromophore. 6 Conversely, we define the partial pathlength effect as the erroneous measured change in a chromophore concentration due to large variance in the photon pathlengths or to incorrect estimates of the latter. 6,8,9 Finally, an alternative hyperspectral illumination and detection configuration, as well as different data processing methods, are also explored and tested to find which could be the ideal HSI methodology to efficiently and reliably monitor hemodynamics and metabolism in the exposed cortex. This last aspect is significant in the context of designing and operating an HSI benchtop system that can experimentally achieve the same level of performances in vivo on small animal models, such as mice and rats. complex domain using a volumetric mesh with triangular surfaces. This modeling technique greatly improves the accuracy of the solutions when modeling objects with curved and complex boundaries, as well as providing an efficient way to sample the problem domain. Thanks to that and to the use of a fast-ray tracing algorithm using Plücker coordinates for rapidly calculating tetrahedron intersections, MMC is also able to efficiently speed up computational time and use less memory during the simulation. 10 MMC is coupled with a mesh-generation and processing toolbox called iso2mesh, 17,18 used to create a volumetric meshed domain that replicates the geometry and structure of cerebral tissue and vasculature from a two-dimensional (2-D) in vivo image of the exposed cortex. Finally, recent releases of the MMC package have implemented the capability to also simulate arbitrary wide-field sources and detectors over large surface areas using mesh retessellation algorithms with high computational efficiency. 12,19 This aspect is crucial for the simulation of HSI, due to the requirement of accurate and reliable representation of 2-D illumination and detection patterns that are characteristic of this optical imaging technique.

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.
The grayscale image of the exposed cortex, showing a 1.2 × 1.2 mm field of view (FOV) of the surface of the brain of a mouse and composed of 400 × 400 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 × 1.2 × 0.1 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 × 2.4 × 1 mm slab reproducing the surrounding mouse subpial gray matter. The extra layers added to the 1.2 × 1.2 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 μ a , scattering coefficient μ s , anisotropy g, and refractive index n). The medium that replicates the mouse subpial gray matter is considered to be made of water (H 2 O), lipid (fat), different concentrations of 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 μm in diameter, respectively) includes water, fat, as well as 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. 20 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. 21 The scattering coefficient μ s ðλÞ is considered to be dependent only on the given wavelength λ of the incident photon packet. 20 The absorption coefficient μ a ðλÞ of each medium of the simulated domain is estimated as the sum of the single absorption coefficients, μ a;H 2 O ðλÞ, μ a;fat ðλÞ, μ a;HbO 2 ðλÞ, μ a;HHb ðλÞ, μ a;oxCCO ðλÞ, and μ a;redCCO ðλÞ, at the given wavelength λ, of the major chromophores composing the medium, i.e., water, fat, HbO 2 , HHb, oxCCO, and redCCO, respectively, and weighted accordingly to their content in it. 20 The data for μ a;H2O ðλÞ and μ a;fat ðλÞ in the NIR range are taken from Matcher et al., 22 for water, and van Veen et al., 23 for fat ( Table 4 in Appendix). The values of μ a;HbO 2 ðλÞ, μ a;HHb ðλÞ, μ a;oxCCO ðλÞ, and μ a;redCCO ðλÞ are calculated from the molar extinction coefficients ε HbO 2 ðλÞ, ε HHb ðλÞ, ε oxCCO ðλÞ, and ε redCCO ðλÞ of 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 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. 20 The molar extinction coefficients ε HbO 2 ðλÞ and ε HHb ðλÞ of HbO 2 and HHb are taken from Matcher et al., 24 whereas the molar extinction coefficients ε oxCCO ðλÞ and ε redCCO ðλÞ of oxCCO and redCCO were measured by John Moody at the University of Plymouth in the bovine heart 6 ( 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 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.
dimensions equal to 0.6 × 0.6 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 λ are launched from the surface of the planar source and evenly distributed over a 1.2 × 1.2 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 [μ a ðλÞ and μ s ðλÞ equal to 0 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 × 1.2 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.

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 × 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 × 185 pixels (6.5-μ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 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. 6 For each simulated condition, image hypercubes and average total photon pathlength hypercubes are reconstructed. Light attenuation changes ΔA k;l (λ) between the simulated baseline and hypoxia are then calculated for each pixel k, l (for k, l ¼ 1: : : 185) and each wavelength λ from the photon intensities I k;l (λ) of the image hypercubes as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 6 1 5 ΔA k;l ðλÞ ¼ −log 10 I k;l;hypoxia ðλÞ I k;l;baseline ðλÞ : From Eq. (1), hemodynamic and metabolic maps charting the relative changes in concentrations Δ½HbO 2 , Δ½HHb, and Δ½oxCCO of 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 ΔA k;l (λ), pixel by pixel. 6,25 Therefore, for each pixel k, l, the following system of algebraic equations is set as   (2) where ε diffCCO ðλÞ are the oxidized-reduced difference molar extinction coefficients of CCO 6 ( Table 4 in Appendix), PL k;l (λ) are the values, in each pixel, of the mean between the average total photon pathlengths in the baseline and hypoxic conditions obtained from the corresponding hypercubes, whereas M is the total number of wavelengths selected for the specific simulation. The hemodynamic and metabolic maps are finally obtained by solving in all the pixels the corresponding systems of algebraic equations in Eq. (2) for the three unknowns Δ½HbO 2 , Δ½HHb, and Δ½oxCCO, using the Moore-Penrose pseudoinverses of the matrices of the molar extinction coefficients. 26,27 3

Computational Studies
In the first computational study (study 1), feasibility and performances of HSI are assessed by running the MC framework for the maximum allowable number of wavelengths (121) in the range 780 to 900 nm and simulating the two conditions previously described. The capability of HSI to reconstruct correct hemodynamic and metabolic maps is evaluated, in particular regarding image quality, as well as the accuracy in the quantification of the relative changes in the concentrations of HbO 2 , HHb, and oxCCO. In addition, corrections in the algorithm for the analysis of the simulated data are introduced and explored to check if the accuracy in the calculated estimates of Δ½HbO 2 , Δ½HHb, and Δ½oxCCO can be improved, as well as to reduce any cross talk or partial pathlength effects during data postprocessing. The second study (study 2) with the MC framework is aimed at understanding how the performances of HSI in monitoring hemodynamics and metabolism are influenced by the specific choice of the wavelengths. Different combinations and numbers of wavelengths in the NIR range are tested in order to find an optimal selection of the spectral bands for maximizing precision of the quantitative data.
Cross talk between hemoglobin and CCO and partial pathlength effects are the main targets for the third study (study 3): the MC framework is used to examine the magnitude of the errors introduced by these factors in the reconstructed maps of hemodynamics and metabolism and to verify the physiological origin of the optical signals that are measured from the simulated data. This is done by comparing the realistic scenario tested in study 1 with ideal and hypothetical scenarios, where one or more concentrations of the chromophores remain constant between the two conditions.
The fourth and final study (study 4) explores the implementation of localized hyperspectral illumination and detection on the simulated domain, as a way to identify the best configuration to efficiently apply HSI to the measurement of the hemodynamic and metabolic states of the exposed cortex.

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, 20 from which the absorption properties used in the simulations are obtained.
At the onset of the acute hypoxic condition, an oxygen saturation drop Δ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. 31 These two simulated phenomena correspond to a theoretical increase Δ½HHb in the concentration of HHb of þ1162.79 and þ43.60 μM, in the pial vasculature and in the subpial gray matter, respectively, as well as to a theoretical decrease Δ½HbO 2 in the concentration of HbO 2 equal to −465.12 and −17.44 μM, in the pial vasculature and in the subpial gray matter, respectively. For the metabolic response, it is assumed that the relative concentration change Δ½oxCCO of oxCCO in the subpial gray matter is equal to −3 μM (this is mirrored by an equivalent increase in [redCCO]). 32,33

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

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)  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 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 HbO 2 and HHb or from the influence of the variance of the photon pathlengths.

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 × 1.2 to 0.2 × 0.2 mm (using the same number of pixels, 185 × 185, in the reconstruction). This is obtained by decreasing the dimension of the source area from 0.6 × 0.6 to 0.1 × 0.1 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. Figure 3 depicts the two hemodynamic maps, for Δ½HbO 2 and Δ½HHb, and the metabolic map of Δ½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 Similarly, a minor hemodynamic response from 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 Δ½HbO 2 and Δ½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 μm in diameter) are partially resolved in the metabolic map.

Study 1
Evaluation of the accuracy in quantifying the correct relative changes in the concentrations of HbO 2 , HHb, and oxCCO is performed by calculating and analyzing the spatial averages of the concentration changes Δ½HbO 2 , Δ½HHb, and Δ½oxCCO in specific regions of interest (ROIs)  Fig. 3(a). The concentrations changes for each chromophore are spatially averaged across the pixels of each ROI. The values of the averages hΔ½HbO 2 i, hΔ½HHbi, and hΔ½oxCCOi in the two ROIs are reported in Table 2 and compared with the   Table 2 Comparison between the spatial average changes hΔ½HbO 2 i, hΔ½HHbi, and hΔ½oxCCOi in the concentrations of HbO 2 , 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 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 Δ½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 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 pial vasculature, across all the wavelengths and between both simulated conditions. (e) Map of the correction factors CF 00 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. 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 Δ½HbO 2 and Δ½HHb in the pial vessels, as well as the occurrence of the spurious measured changes Δ½oxCCO in the same vascular medium, resulting from applying MBLL from Eq. (2) and using the average total pathlength of the detected photons.
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 HbO 2 and HHb in the pial vasculature. Two maps of correction factors, CF 0 k;l and CF 00 k;l (for each pixel k, l), are produced: (1) CF 0 k;l are the means across all the selected M wavelengths (in this case M ¼ 121) of the ratios between the average total pathlengths PL k;l ðλÞ travelled by the detected photons and the average partial pathlengths PL k;l;vessel ðλÞ they travelled in the pial vasculature (for each wavelength λ). (2) CF 00 k;l are the means across all the selected M wavelengths of the ratios between the average total pathlengths PL k;l ðλÞ travelled by the detected photons and the average partial pathlengths PL k;l;gray ðλÞ they travelled in the subpial gray matter (for each wavelength λ). Thus The two sets of correction maps can be obtained using Eq. 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 Δ½HbO 2 Ã and Δ½HHbÃ of the changes in concentration of HbO 2 and HHb in the hemodynamic maps only: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 2 7 5 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 Δ½HbO 2 Ã, Δ½HHbÃ, and Δ½oxCCOÃ of the changes in concentration of HbO 2 , HHb. and oxCCO: In both sets of Eqs. (4) and (5), CF 0 k;l and CF 00 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 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 Δ½HbO 2 and Δ½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 hΔ½HbO 2 i, hΔ½HHbi, and hΔ½oxCCOi of the concentration changes of HbO 2 , HHb, and oxCCO in the same two ROIs used for the uncorrected maps. Table 2 shows that the quantified values of hΔ½HbO 2 i and hΔ½HHbi in the pial vasculature are now much closer to the simulated theoretical values (−474.7 μM for hΔ½HbO 2 i and 1192.2 μM for hΔ½HHbi), with estimation errors of about 2.06% and 2.87%, for 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 hΔ½HbO 2 i, hΔ½HHbi, and hΔ½oxCCOi 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 Δ½HbO 2 , Δ½HHb, and Δ½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 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 Δ½oxCCO before and after correction are almost overlapping, as well as for the values of Δ½HbO 2 and Δ½HHb in the pixels on the subpial gray matter, in both Figs. 5(b) and 5(c).

Study 2
In study 2, similar hemodynamic and metabolic maps for Δ½HbO 2 , Δ½HHb and Δ½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), 34 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 hΔ½HbO 2 i, hΔ½HHbi, and hΔ½oxCCOi of the relative changes in the concentrations of 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 hΔ½HbO 2 i, hΔ½HHbi, and hΔ½oxCCOi 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 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. 34 that were estimated for bNIRS, extending them also to HSI targeting brain metabolism.

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 Δ½HbO 2 , Δ½HHb, and Δ½oxCCO. No postprocessing correction was performed in this case.
In the absence of any replicated hemodynamic response in the simulations, the hemodynamic maps do not display any contrast provided by the changes Δ½HbO 2 and Δ½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 hΔ½HbO 2 i, hΔ½HHbi, and hΔ½oxCCOi of the relative changes in the concentrations of 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 hΔ½HbO 2 i and hΔ½HHbi 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 hΔ½oxCCOi in the ROI corresponding to the pial vasculature (−3.018 μM). The magnitude of these is still in the same order of the simulated relative change in concentration of oxCCO (−3 μM). Finally, a more accurate quantification of the relative change Δ½oxCCO in the metabolic map emerges from the calculation of the spatial average hΔ½oxCCOi in the ROI corresponding to the gray matter (−3.098 μM), which appears closer to the effective simulated change in oxCCO than the result from study 1 (−3.157 μM). The corresponding estimation error is about 3.27% (against 5.2% in study 1).
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 Δ½HbO 2 , Δ½HHb, and Δ½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 Δ½HbO 2 and Δ½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. Table 3 Spatial average changes hΔ½HbO 2 i, hΔ½HHbi, and hΔ½oxCCOi in the concentrations of HbO 2 , 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.  The analysis of the spatial averages hΔ½HbO 2 i, hΔ½HHbi, and hΔ½oxCCOi 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 μ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 (hΔ½oxCCOi equal to −0.083 μ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.

Study 4
The new HSI configuration tested in the fourth and final study with the MC framework, implementing and simulating a 0.2 × 0.2 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 × 185 pixels. Since the FOV is now smaller (0.2 × 0.2 mm, against the previous 1.2 × 1.2 mm FOV), the size of the pixels decreases from 6.5 to 1.08 μ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 hΔ½HbO 2 i, hΔ½HHbi, and hΔ½oxCCOi of the relative changes in the concentrations of HbO 2 , HHb, and oxCCO are calculated from the hemodynamic and metabolic maps for an ROI of 60 × 60 pixels concentric with the 0.2 × 0.2-mm FOV. The ROI corresponds to a square region of about 65 × 65 μ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 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 hΔ½HbO 2 i and hΔ½HHbi in the ROI for the new HSI configuration stand at −17.67 μM for hΔ½HbO 2 i and 44.64 μM for hΔ½HHbi, against −19.39 and 48.72 μM, respectively for study 1. The new values are much closer to the theoretical simulated changes in the concentration of HbO 2 and HHb in the subpial gray matter (−17.44 μM for HbO 2 and 43.60 μM for HHb). The quantification errors for hΔ½HbO 2 i and hΔ½HHbi 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 μM for hΔ½oxCCOi with the 0.2 × 0.2 mm FOV, compared to −3.156 μM with the 1.2 × 1.2 mm FOV) and thus closer to the simulated metabolic response (−3 μM). The estimation error of hΔ½oxCCOi 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).

Discussion
Preliminary studies with the MC framework proved the suitability of HSI as an optical imaging modality for spatially and quantitatively monitoring the hemodynamic and metabolic response of the exposed cortex to hypoxia: the hemodynamic response was correctly localized in the pial vasculature with high spatial resolution, whereas changes in the concentrations of HbO 2 , HHb, and oxCCO were accurately estimated in the subpial gray matter. The results obtained with the MC framework, regarding the monitoring of hemoglobin oxygenation and blood perfusion on the exposed cortex during hypoxia, are comparable and consistent with previous in vivo HSI studies using primarily visible and NIR light. [35][36][37][38] We speculate that both the underestimation in the quantification of the changes in the concentrations of hemoglobin in the pial vasculature, as well as the occurrence of spurious measured changes in the concentration of oxCCO in the same region (where the MC framework did not simulate such concentration change), are only caused by the large differences in the partial pathlengths of the detected photons between the two media. These differences are not taken into account when applying MBLL since this only consider the total pathlength of the photons. Relevant insights on this phenomenon are highlighted by the findings from the first part of the third study: the hemodynamic maps obtained from simulating only metabolic response during hypoxia demonstrated that negligible relative changes in the concentrations of HbO 2 and HHb occur in both the pial vessels and the surrounding subpial gray matter in the absence of hemodynamic response. Spurious signals from oxCCO still appear in the pial vasculature, in the same order of magnitude of the relative changes in concentrations of oxCCO due to actual metabolism. Nevertheless, quantification of Δ½oxCCO in the central subpial tissue was still accurate and closer to the actual change in oxCCO than the results of the study 1. This suggests that partial pathlength effects do not affect significantly the quantification of the metabolic response in the same region and the changes in CCO do not arise as a cross talk from the hemoglobin signals. Thus the measured data obtained for Δ½oxCCO in the subpial gray matter are primarily connected to the optical signature of CCO, proving the efficacy of HSI to retrieve metabolic signal in the exposed cortex. This conclusion is furtherly supported by the results from the second part of third study, where oppositely only the hemodynamic response in the domain was simulated, showing no changes in Δ½oxCCO in both the hemodynamic and metabolic maps, as expected.
We then proposed a postprocessing, spatially selective correction taking into account the differences in the partial pathlengths of the detected photons, which enhanced image contrast in the hemodynamic maps and the accuracy of the quantification of the hemodynamic response in the pial vasculature with an estimation error of <3%.
No major differences were found in the outcomes of the second study using different numbers and combinations of wavelengths for hyperspectral illumination, compared to the results of the first study using the maximum allowable number of 121 wavelengths. Thus we showed that reducing the number of simulated wavelengths (as long as they are evenly sampled across in the selected NIR range) down to an optimal combination of only eight does not significantly affect the quality of the hyperspectral data, nor provide significant differences in the accuracy of the quantification of both the hemodynamic and metabolic responses. This is consistent with the results by Arifler et al. 34 on wavelength optimization for simultaneous monitoring of hemoglobin and CCO via bNIRS.
The findings in study 2 can be advantageous for designing an experimental benchtop HSI system for monitoring hemodynamic and metabolism in the exposed cortex of small animals since reducing the number of necessary wavelengths needed to obtain accurate data decreases complexity and cost of the instrumentation, as well as computational burden to process a smaller volume of hyperspectral data.
Finally, the results of the fourth study provided a preliminary proof of concept for the use of an alternative HSI imaging approach based on localized and selective hyperspectral illumination and detection, to increase the accuracy in the quantification of the hemodynamic and metabolic response in the subpial gray matter. This approach improved the accuracy of the quantification of relative changes in the concentrations of HbO 2 , HHb, and oxCCO in that region (estimation errors of <2.5% for hemoglobin and 4% for CCO) without the need of postprocessing correction.
The new hyperspectral illumination and detection approach, as well as the hyperspectral processing algorithms here reported, can be implemented in a benchtop HSI system and validated under controlled experimental conditions, e.g., using blood and yeast liquid phantoms. 39 Moreover, the tested HSI configuration could be further explored and developed in the future, e.g., by spatially scanning larger FOVs including both vasculature and gray matter or by applying modulated illumination techniques similar to those used in spatial frequency-domain imaging and structured illumination imaging. 40,41 The findings of the four studies reported here can be translated into an experimental setting and could improve the performances of any benchtop NIR HSI system that targets the relative changes of concentration of HbO 2 , HHb, and oxCCO on the exposed cortex, whereas the MC framework can be easily coupled with the instrumentation to aid hyperspectral data acquisition and reconstruction for in vivo applications.
Further studies and developments of the MC HSI framework can be explored in the future, which can consist of: (1) refining the simulated domain to include also subpial microvasculature; (2) considering potential differences in the scattering properties between pial vasculature and subpial gray matter; and (3) simulating and investigating additional cerebral physiological conditions besides hypoxia, such as hypercapnia, hyperemia, and other abnormal brain hemodynamic and metabolic responses.

Conclusion
A MC framework simulating NIR HSI quantitative monitoring of the hemodynamic and metabolic states of the exposed cortex has been here described and tested for a realistic meshed domain, generated from in vivo data and replicating mouse cerebral pial vasculature and subpial gray matter. We demonstrated its efficacy for modeling hyperspectral illumination and data acquisition, using up to 121 wavelengths in the NIR range between 780 and 900 nm, as well as for reproducing measurements of the relative changes in the concentrations of HbO 2 , HHb, and oxCCO in the form of hemodynamic and metabolic maps. The MC framework can be also flexibly tuned for different applications, as well as different numbers and ranges of wavelengths, making it a powerful and reusable tool for simulating HSI, due to its ability to reproduce complex meshed domains of various types of tissue from real data. Table 4 in this section reports the absorption coefficients μ a;H 2 O ðλÞ and μ a;fat ðλÞ of water and fat, respectively, and the molar extinction coefficients ε HbO 2 ðλÞ, ε HHb ðλÞ, ε oxCCO ðλÞ, and ε redCCO ðλÞ of HbO 2 , HHb, oxCCO, and redCCO, respectively, which were used to calculate the absorption properties of the simulated domain, with corresponding references or sources. The oxidizedreduced difference molar extinction coefficients ε diffCCO ðλÞ of CCO used in Eq. (2) are also reported in Table 4. Table 4 Values in the NIR range between 780 and 900 nm of: (1) the absorption coefficients μ a;H 2 O and μ a;fat of water and fat, respectively; (2) the molar extinction coefficients ε HbO 2 , ε HHb , ε oxCCO , and ε redCCO of HbO 2 , HHb, oxCCO, and redCCO, respectively; and (3) the oxidizedreduced difference molar extinction coefficients ε diffCCO of CCO. References and sources of these values are also provided.

Appendix
Matcher et al. 22 van Veen et al. 23 Matcher et al. 24 Matcher research interest is to develop near-infrared spectroscopy instruments and methodologies for biomedical applications.
Ilias Tachtsidis received his PhD from the UCL, UK, in 2005. He is a senior member of the Biomedical Optics Research Laboratory at UCL, a senior Wellcome Trust fellow, an associate professor in biomedical engineering, and a leader of the Multimodal Spectroscopy Group. He works on the development and application of NIRS techniques to monitor the function of the brain both in health and disease, including adults and neonatal brain injury.