2.1.

Digital Instrument Simulator

The digital instrument simulator is based on a realistic digital phantom of an exposed cortex computed with Monte Carlo simulations. The code is publicly available on GitHub. This simulator aims to model the light propagation in the exposed brain. It can be used to take into account the inhomogeneities of the optical properties of the exposed brain cortex for functional brain mapping studies and to reduce the quantification errors of biomarkers of brain functionality. The Monte Carlo framework is composed of several steps, see Fig. 1. First, a color image of the exposed cortex was taken during the neurosurgery, see Sec. 2.1.1. This image was segmented into three classes (i.e., grey matter and large and small blood vessels, see Sec. 2.1.2). A brain volume was then modeled, and a white Monte Carlo approach was used to estimate the partial path length, the position, and the angle of exiting packets of photons, see Sec. 2.1.3. Finally, the information of the exiting packets of photons was converted into images of diffuse reflectance and mean path length, see Sec. 2.1.4.

Fig. 1

Flowchart of digital instrument simulator. (a) Acquisition of a color image during a neurosurgical operation. (b) Three-class segmentation of the color image (grey matter and large and small blood vessels). (c) 3D tissue modeling. (d) White Monte Carlo simulations of light propagation in the tissue. (e) Reconstruction of diffuse reflectance and mean path length images from white Monte Carlo simulations. (f) Diffuse reflectance image at 500 nm. (g) Mean path length image at 500 nm.

2.1.2.

Image segmentation and brain volume modeling

The color image was segmented into three classes (i.e., grey matter and large and small blood vessels) using morphological operations, see Fig. 2.

• 1. The color image ($A$) was converted into grayscale ($C$) using the Python library OpenCV (v4.8.1).³³

• 2. The grayscale image ($C$) was converted into a binary image ($D$) using a Gaussian adaptive threshold with a block size of $5\times 5\mathrm{mm}$. The block size was automatically converted into an odd number of pixels using the resolution of the image ($\text{ceil}(5/0.073)=69\text{pixels}$). A block size of 5 mm corresponded to the optimal size for detecting the blood vessels in the human brain. A smaller or larger size would impact the detection of the large and small blood vessels, respectively. The adaptive threshold was only applied to the exposed cortex by applying the mask ($B$). The latter is created manually. Saturated pixels due to specular reflection were removed from the image ($D$) with a simple thresholding of image ($C$).

• 3. The small vessels were removed from the binary image ($D$) using a morphological opening with a disk of 0.5-mm diameter (7 pixels) as a structuring element. Using the resulting image, a morphological closing was computed to get the image ($F$), the large blood vessel mask.

• 4. A binary exclusive disjunction (xor) was computed between images ($D$) and ($F$) to get the small blood vessel mask ($E$).

• 5. The mask of grey matter ($G$) corresponded to the rest of the pixels.

Fig. 2

Segmentation of a color image (a) into masks of large blood vessels (f), small blood vessels (e), and grey matter (g). 1: Conversion of color (a) into grayscale image (c). 2: Gaussian adaptive thresholding to obtain the binary image (d). 3: Computation of the large blood vessel mask (f). 4: Computation of the small blood vessel mask (e). 5: Computation of the grey matter mask (g).

Once the image was segmented into three classes, we modeled the brain volume. First, the binary segmentation masks were encased in a larger image having the label of grey matter. This operation is performed to prevent a blood vessel label from appearing at the edges of the image. The objective is to avoid photon loss due to the boundary effects during the Monte Carlo simulations, see Sec. 2.1.3. The binary segmentation masks were expanded along the vertical direction ($z$ axis) on 2 cm to convert the binary images into binary volumes, see Figs. 3(a) and 3(b). Then, we modeled the 3D blood vasculature with morphological erosion, see Fig. 3(c). The thickness of the blood vessels was calculated on the basis of the vessel’s diameter. The structuring element used for the erosion was a disk of 0-pixel diameter for $z=0$ (in pixels) and was increased to 1 pixel while increasing $z$ axis. The binary volumes of the three classes were finally merged together with a final isotropic resolution of $73\mu \mathrm{m}$, see Fig. 1. With this model of the cerebral vasculature, the blood vessels were only modeled on the surface of the volume. Deep blood vessels (which are not visible in color imaging) were not included.

Fig. 3

Steps for modeling a brain volume with blood vessels. (a) Segmentation mask of the grey matter and the blood vessels. (b) Replication of the binary segmentation masks along the vertical direction ($z$ axis) on 2 cm. (c) Morphological erosion to create the blood vessel depth.

2.1.3.

White Monte Carlo simulations

Once the brain volume was defined, we computed Monte Carlo simulations of light propagation in the domain using MCX³⁴ software, see Fig. 4.

Fig. 4

White Monte Carlo simulations. (a) Schematics for the Monte Carlo simulations with Ref. 34. (b) Flowchart of white Monte Carlo simulations.

A planar light source was placed on the top of the volume to homogeneously illuminate the brain volume. We set the Fresnel reflection as boundary conditions for the top and bottom faces of the volume (i.e., an exiting packet of photons is lost). Cyclic boundary conditions were applied for the sides of the volume to avoid photon loss due to the boundary effects. Each voxel of the volume included optical properties, see Table 1.

Table 1

Optical properties used in the white Monte Carlo simulations. λ denotes the wavelength in nanometers.

	Grey matter	Large blood vessels	Small blood vessels
Absorption coefficient ${\mu }_a$ (${\mathrm{mm}}^{-1}$)	0	0	0
Scattering coefficient ${\mu }_s$ (${\mathrm{mm}}^{-1}$)	$4.08{\left(\frac{\lambda }{500}\right)}^{-3.089}$35	$2.2{\left(\frac{\lambda }{500}\right)}^{-0.66}$35	$2.2{\left(\frac{\lambda }{500}\right)}^{-0.66}$35
Anisotropy coefficient $g$	0.8536	0.93537	0.93537
Refractive index $n$	1.3638	1.437	1.437

For these simulations, we implemented a white Monte Carlo approach, and the absorption coefficient was set to $0{\mathrm{mm}}^{-1}$. With this technique, the packets of photons were not affected by the absorption. The advantage of the white Monte Carlo approach is that the absorption can be considered “a posteriori,” using the microscopic Beer–Lambert law.³⁹ Thus, when changing the absorption parameters of the volume (i.e., to simulate a brain activation for example), the simulation does not need to be run again, speeding up the processing. In this case, only one simulation per wavelength is required, considering that the scattering properties of the medium do not change. A scattering change can be considered, but a new simulation with a new set of ${\mu }_s$ would need to be rerun. For each wavelength (from 400 to 1000 nm by steps of 10 nm), three outputs were stored for every packet of photons exiting the top face of the volume: the position of the exiting packet of photons $(x,y)$, its exiting angle, and the partial path length (the length that each photon has spent in each class of the domain).

2.2.

Identification of the Optimal Wavelength for Hemodynamic and Metabolic Monitoring

Using the digital instrument simulator (see Sec. 2.1), we developed an optimization procedure based on a genetic algorithm (differential evolution) to identify the best wavelength combination in the visible and near-infrared range to minimize the quantification error (root mean square error in $M$) in ${\mathrm{HbO}}_2$, Hb, oxCCO, oxCytb, and oxCytc compared with ground truth values, see Fig. 6.

2.2.1.

Model of cerebral activation

Data used for the optimization procedure were generated with the digital instrument simulator. A simple model was considered, which includes non-activated and activated grey matter, see Fig. 5(a). In Fig. 5(b), we modeled a cerebral activity in the activated grey matter by changing the concentration of chromophores relative to their nominal values. The modeled concentration changes that are the ground truth in the optimization procedure are indicated $\mathrm{\Delta }{C}_{\mathrm{GT}}$ in the rest of the paper. A simple temporal perturbation of ${\mu }_a$ was performed: rest (nominal values, see Table 2) and activity (nominal value $+\mathrm{\Delta }C$). The $\mathrm{\Delta }C$ values were $+5$ and $-3.75\mu \mathrm{M}$ for ${\mathrm{HbO}}_2$ and Hb,^9,48 respectively, and $+0.5\mu \mathrm{M}$ for oxCCO,⁴⁹ oxCytb, and oxCytc, mirrored by concentration changes of $-0.5\mu \mathrm{M}$ for redCCO, redCytb, and redCytc (the concentration of cytochromes does not change during the cerebral activity).

Fig. 5

Monte Carlo model used for the optimization procedure. (a) Volume of non-activated grey matter (green) and activated grey matter (blue). (b) Spatial profile of the modeled concentration changes $\mathrm{\Delta }{C}_{\mathrm{GT}}$ indicated by the dotted line in panel (a).

2.2.2.

Optimization procedure

The most straightforward method to determine $k$ optimal wavelengths for cerebral monitoring is to test all possible $k$-element combinations out of 601 wavelengths (from 400 to 1000 nm) and to identify the combination that produces the smallest estimation errors when compared with the ground truth $\mathrm{\Delta }{C}_{\mathrm{GT}}$. As the investigated spectral range is large, such an exhaustive search results in a large time complexity. Thus, we proposed to use a genetic algorithm⁵⁰ to determine the $k$ optimal wavelengths through an optimization procedure, see Fig. 6. We used the differential evolution method that optimizes a problem by iteratively trying to improve a candidate solution with regard to the ground truth. We used the function differential evolution from the Python library Scipy (v1.10.0).⁵¹

Fig. 6

Flowchart of the optimization method, including definitions of the variables used in the procedure.

The optimization procedure consists of finding the best $k$-element combinations out of 601 wavelengths that minimize a cost function. This function takes as input a wavelength array, with values that are limited to the simulated spectral range (between 400 and 1000 nm). The cost function returns the root mean square error RMSE (in $M$) computed between noisy concentration changes measured with the modified Beer–Lambert law and the ground truth $\mathrm{\Delta }{C}_{\mathrm{GT}}$

Eq. (6)

$$\mathrm{RMSE}=\sqrt{\sum _{i=1}^{N_{\text{Noise}}}\frac{(\mathrm{\Delta }{C}_{\mathrm{mes}}^{\text{Noise}}-\mathrm{\Delta }{C}_{\mathrm{GT}}{)}^2}{N_{\text{Noise}}}},$$

where $N_{\text{Noise}}=1000$ is the number of iterations for the noise addition. Noise was added to the simulated intensities to identify the $k$-wavelengths the most robust to the noise. The intensities were converted into concentration changes $\mathrm{\Delta }{C}_{\mathrm{mes}}^{\text{Noise}}$ measured with the modified Beer–Lambert law

In Eq. (7), ${ϵ}_n({\lambda }_1)$ is the molar extinction coefficient of the chromophore $n$ (in $M^{-1}.{\mathrm{mm}}^{-1}$), see Fig. S1 in the Supplementary Material. $\mathrm{\Delta }A({\lambda }_1{)}^{\text{Noise}}$ is the noisy attenuation change measured for ${\lambda }_1$, such as

Eq. (8)

$$\mathrm{\Delta }A({\lambda }_1{)}^{\text{Noise}}={\log }_{10}(\frac{{\varphi }_{\text{rest}}({\lambda }_1)+\gamma }{{\varphi }_{\text{activity}}({\lambda }_1)+\gamma }),$$

where $\gamma$ is a zero-mean Gaussian noise with a standard deviation that was the same for each wavelength: $\sigma =\frac{{\varphi }_{\text{rest}}}{\mathrm{SNR}}$. With ${\varphi }_{\text{rest}}$, the diffuse reflectance simulated during rest and $\mathrm{SNR}=400$ the signal-to-noise ratio of the instrument. This value corresponds to the experimental devices we used in a previous study,²² but this could be adapted for another experimental device. The concentration changes were obtained through matrix inversion in the least square sense.

3.1.

Digital Instrument Simulator

The digital instrument simulator aimed to simulate hyperspectral images of the exposed brain cortex for 601 wavelengths (between 400 and 1000 nm). In Fig. 7, we represented the images of diffuse reflectance and mean path length at 500 and 900 nm reconstructed at the surface of the tissue considering that large and small blood vessels are arteries. The images were reconstructed with a $5\times 5$ binning, which leads to a resolution of $365\mu \mathrm{m}$. The red, grey, and magenta points indicate the position of points for a large blood vessel, for grey matter, and for a small blood vessel, respectively.

Fig. 7

Images of diffuse reflectance and mean path length reconstructed with a $5\times 5$ binning at 500 and 900 nm. (a) Diffuse reflectance image simulated at 500 nm. (b) Mean path image simulated at 500 nm. (c) Diffuse reflectance image simulated at 900 nm. (d) Mean path image simulated at 900 nm. The large and small blood vessels are considered arteries. The red, grey, and magenta points indicate the position of points for a large blood vessel, for grey matter, and for a small blood vessel, respectively.

In Fig. 8, we represented the diffuse reflectance and mean path length spectra measured for the three points identified in Fig. 7. Solid and dotted lines indicate that large and small blood vessels are considered arteries and veins, respectively. For both cases, we can observe the effect of the blood vessels on the simulated quantities. For wavelengths lower than 600 nm, the diffuse reflectance was almost $0\mu {\mathrm{m}}^{-2}$ on the large blood vessels; this means that hardly any photons left the large blood vessel. Diffuse reflectance values increased when measuring small blood vessels, indicating a greater contribution from outgoing photons. The higher values can be found for grey matter regions. For wavelengths higher than 600 nm, diffuse reflectance values were almost the same for grey matter and small blood vessels, and lower values can be found at the level of the large blood vessels. The mean path length spectra measured for grey matter and the small blood vessels were almost the same for all wavelengths, but lower values were obtained for the large blood vessels. We can observe that the values of the mean path length are not impacted if all vessels are considered either veins or arteries. However, we can observe differences among the diffuse reflectance spectra measured at the level of the large blood vessels for wavelengths higher than 600 nm. There is a drop in intensity in the red and near infrared for veins compared with arteries. This corresponds to the fact that veins appear darker or bluer than arteries.

Fig. 8

Diffuse reflectance and mean path length spectra measured on grey matter, a large blood vessel, and a small blood vessel. (a) Image of the cerebral cortex segmented into three classes (grey matter and large and small blood vessels). The three black points indicate the position of points for a large blood vessel, for grey matter, and for a small blood vessel. (b) Diffuse reflectance spectra collected at the location of the three black points. (c) Mean path length spectra collected at the location of the three black points. Solid and dotted lines indicate that large and small blood vessels are considered arteries and veins, respectively.

3.2.

Identification of the Optimal Wavelength for Hemodynamic and Metabolic Monitoring

In Fig. 9, we represented the optimal combination of 2, 4, 6, 8, and 10 wavelengths for monitoring $C_{{\mathrm{HbO}}_2}$ and $C_{\mathrm{Hb}}$ changes in activated grey matter. These wavelengths are represented in Fig. 9(a) by dots. We also represented the configuration of two wavelengths used by Bouchard et al.,²⁹ the four wavelengths used by White et al.³⁰ and the broadband spectral range between 780 and 900 nm used by Bale et al.^14,17 and Arifler et al.³¹ The quantification errors related to these configurations are plotted in Fig. 9(b). The bars represent the RMSE averaged over 1000 noisy measurements, see Eq. (6), and the vertical line represents the standard deviation in the error. Finally, the concentration changes in $\mathrm{\Delta }{C}_{{\mathrm{HbO}}_2}$ and $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ averaged over the 1000 noisy measurements are plotted in Figs. 9(c) and 9(d), respectively.

Fig. 9

(a) Optimal wavelength for monitoring $C_{{\mathrm{HbO}}_2}$ and $C_{\mathrm{Hb}}$ changes in activated grey matter. (b) Quantification errors (in $\mu \mathrm{M}$). (c) and (d) Concentration changes in ${\mathrm{HbO}}_2$ and Hb, respectively. Several configurations are represented: two wavelengths proposed by Bouchard et al.,²⁹ four wavelengths proposed by White et al.,³⁰ a broadband spectra proposed by Bale et al.,^14,17 and the combination of 2, 4, 6, 8, and 10 wavelengths identified with our optimization procedure.

Contrary to the configurations proposed by Bouchard et al. and White et al., the wavelengths identified with the optimization procedure do not include the hemoglobin isobestic point at 530 nm. The wavelengths were located in the visible and near-infrared range. Wavelengths were identified before 500 nm, where ${\mathrm{HbO}}_2$ absorption predominates, between 590 and 730 nm where Hb absorption predominates and on either side of the isobestic point at 800 nm. The minimum errors in $\mathrm{\Delta }{C}_{{\mathrm{HbO}}_2}$ and $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ are obtained with the combination of 10 wavelengths calculated with our optimization procedure. Among the literary configurations, White’s spectral configuration minimizes errors in $\mathrm{\Delta }{C}_{{\mathrm{HbO}}_2}$ and $\mathrm{\Delta }{C}_{\mathrm{Hb}}$. Errors in $\mathrm{\Delta }{C}_{{\mathrm{HbO}}_2}$ and $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ obtained with the optimal group of 10 wavelengths were 29% and 23% lower than those obtained with White’s spectral configuration. The mean and standard deviation of the RMSE decrease slightly as more spectral bands are considered. For all wavelengths, the standard deviation and the mean RMSE are lower than those obtained with White’s spectral configuration.

In Fig. 10, we represented the optimal combination of 4, 6, 8, and 10 wavelengths for monitoring $C_{{\mathrm{HbO}}_2}$, $C_{\mathrm{Hb}}$, and $C_{\mathrm{oxCCO}}$ changes in activated grey matter. These wavelengths are represented in Fig. 10(a) by dots. We also represented the configuration of four wavelengths proposed by Arifler et al.³¹ and the broadband spectral range used by Bale et al.^14,17 The quantification errors are represented in Fig. 10(b) and the measurements of $\mathrm{\Delta }{C}_{{\mathrm{HbO}}_2}$, $\mathrm{\Delta }{C}_{\mathrm{Hb}}$, and $\mathrm{\Delta }{C}_{\mathrm{oxCCO}}$ are plotted in Figs. 10(c), 10(d), and 10(e), respectively.

Fig. 10

(a) Optimal wavelength for monitoring $C_{{\mathrm{HbO}}_2}$, $C_{\mathrm{Hb}}$, and $C_{\mathrm{oxCCO}}$ changes in activated grey matter. (b) Quantification errors (in $\mu \mathrm{M}$). (c), (d), and (e) Concentration changes in ${\mathrm{HbO}}_2$, Hb, and oxCCO, respectively. Several configurations are represented: four wavelengths proposed by Arifler et al.,³¹ a broadband spectra proposed by Bale et al.,^14,17 and the combination of 4, 6, 8, and 10 wavelengths identified with our optimization procedure.

The minimum errors in $\mathrm{\Delta }{C}_{{\mathrm{HbO}}_2}$, $\mathrm{\Delta }{C}_{\mathrm{Hb}}$, and $\mathrm{\Delta }{C}_{\mathrm{oxCCO}}$ are obtained with the combination of 10 wavelengths calculated with our optimization procedure. With this digital instrument simulator, we modeled an exposed cortex. The simulated light is not absorbed by the skin or the skull of the patient, and the visible light can be collected by the camera. Thus, the optimization procedure can take profit of the peaks of oxidized and reduced CCO extinction spectra in the visible and the near-infrared range to monitor the three chromophores. The errors in $\mathrm{\Delta }{C}_{{\mathrm{HbO}}_2}$, $\mathrm{\Delta }{C}_{\mathrm{Hb}}$, and $\mathrm{\Delta }{C}_{\mathrm{oxCCO}}$ obtained with this configuration were 23%, 48%, and 23% lower than those obtained with broadband spectra between 780 and 900 nm, respectively. For all chromophores, the standard deviation of the RMSE decreases as more spectral bands are considered. For $\mathrm{\Delta }{C}_{{\mathrm{HbO}}_2}$, the mean RMSE decreases slightly with the number of wavelengths, whereas the mean errors for $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ and $\mathrm{\Delta }{C}_{\mathrm{oxCCO}}$ remain constant. For all wavelengths, the mean RMSE is lower than that computed with Bale’s spectral configuration. For $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ and $\mathrm{\Delta }{C}_{\mathrm{oxCCO}}$, the standard deviation is less than or equal to that obtained with Bale’s configuration. However, for $\mathrm{\Delta }{C}_{{\mathrm{HbO}}_2}$, the standard deviation is higher than that obtained with Bale’s configuration.

In Fig. 11, we represented the optimal combination of 6, 8, and 10 wavelengths for monitoring $C_{{\mathrm{HbO}}_2}$, $C_{\mathrm{Hb}}$, $C_{\mathrm{oxCCO}}$, $C_{\mathrm{oxCytb}}$, and $C_{\mathrm{oxCytc}}$ changes in activated grey matter. These wavelengths are represented in Fig. 11(a) by dots. The quantification errors are represented in Fig. 11(b), and the concentration changes are plotted in Figs. 11(c), 11(d), 11(e), 11(f), and 11(g).

Fig. 11

(a) Optimal wavelength for monitoring $C_{{\mathrm{HbO}}_2}$, $C_{\mathrm{Hb}}$, $C_{\mathrm{oxCCO}}$, $C_{\mathrm{oxCytb}}$, and $C_{\mathrm{oxCytc}}$ changes in activated grey matter. (b) Quantification errors (in $\mu \mathrm{M}$). (c), (d), (e), (f), and (g) Concentration changes in ${\mathrm{HbO}}_2$, Hb, oxCCO, oxCytb, and oxCytc, respectively.

The minimum errors to resolve all chromophores are obtained with the combination of 8 and 10 wavelengths. It is interesting to note that the use of a larger number of wavelengths (10) does not significantly reduce quantification errors. With the use of 6, 8, and 10 optimal wavelengths, the total quantification errors were 29, 28, and $28\mu \mathrm{M}$, respectively. The errors in $\mathrm{\Delta }{C}_{{\mathrm{HbO}}_2}$ and $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ were higher than those obtained in Fig. 10. However, with the consideration of cytochromes $b$ and $c$ in the modified Beer–Lambert law, the errors in $\mathrm{\Delta }{C}_{\mathrm{oxCCO}}$ were 19% lower than those obtained in Fig. 10. The noise in the measurements of $\mathrm{\Delta }{C}_{\mathrm{oxCytb}}$ and $\mathrm{\Delta }{C}_{\mathrm{oxCytc}}$ is very important, making it difficult to differentiate between activated and non-activated grey matter.

As proposed in our previous study,²² we can also identify different combinations of wavelengths to independently quantify the concentration changes in ${\mathrm{HbO}}_2$, Hb, oxCCO, oxCytb, and oxCytc, see Fig. 12. For instance, if we want to identify the optimal wavelength to resolve $C_{\mathrm{oxCytb}}$ changes only, the optimization procedure can be run to minimize $\mathrm{\Delta }{C}_{\mathrm{oxCytb}}$ values by only considering three chromophores in the modified Beer–Lambert law [${\mathrm{HbO}}_2$, Hb, and oxCytb, see Eq. (7)].

Fig. 12

Quantification errors (a) and concentration changes obtained with the optimal combination of 4, 6, 8, and 10 wavelengths for a separate quantification. The concentration changes in $C_{{\mathrm{HbO}}_2}$, $C_{\mathrm{Hb}}$, $C_{\mathrm{oxCCO}}$, $C_{\mathrm{oxCytb}}$, and $C_{\mathrm{oxCytc}}$ are plotted in panels (b), (c), (d), (e), and (f), respectively.

With this approach, we represented in Fig. 12(a) the quantification errors obtained with the optimal wavelength combination to separately resolve ${\mathrm{HbO}}_2$, Hb, oxCCO, oxCytb, and oxCytc changes. We also plotted the concentration changes in ${\mathrm{HbO}}_2$, Hb, oxCCO, oxCytb, and oxCytc in Figs. 12(b), 12(c), 12(d), 12(e), and 12(f). In Fig. S2 in the Supplementary Material, we represented the optimal combinations of 4, 6, 8, and 10 wavelengths to separately resolve ${\mathrm{HbO}}_2$, Hb, oxCCO, oxCytb, and oxCytc changes. The separate optimization procedure leads to different spectral configurations for each chromophore.

For all chromophores (expect oxCytb), the mean RMSE decreases with the number of wavelengths. We obtained smaller quantification errors when we applied the optimization procedure on each chromophore separately than when we tried to resolve all the chromophores at the same time, see Table 3.

Table 3

Average of the quantification errors (RMSE in μMol) obtained with two Beer–Lambert systems proposed in the literature and optimal combination of six wavelengths.

4.1.

Digital Instrument Simulator

The digital instrument allows the modeling of intensity maps collected by a camera sensor as well as the estimation of the mean path length of traveled photons through the tissue. This framework could be used to improve clinical and preclinical optical devices for brain mapping applications. In these studies, the modified Beer–Lambert law is used to monitor chromophore changes in an animal^29,30 or a human^7,9,11,52 brain. However, measurements are highly subject to quantification errors if incorrect path lengths are used to resolve chromophore changes. This could be the case in a lot of applications where the inhomogeneities of the optical properties are not taken into consideration. Indeed, a single path length is considered for the whole cerebral cortex, which is usually estimated with a homogeneous volume of grey matter with Monte Carlo simulations²⁹ or using the analytical solution to the diffusion approximation of the radiative transfer equation in a semi-infinite geometry.^30,52 In our study, we proposed to take into account the inhomogeneities of the optical properties with a pixel-wise estimation of the mean path length. As we can see on Figs. 7 and 8, the blood vessels have a great impact on the path length. For large blood vessels having a 2-mm diameter, the path length is almost null in the visible spectra (below 600 nm) and is lower than the path length estimated for grey matter in the red and near-infrared range. However, the path length measured at the level of the small blood vessels ($\approx 0.2\mathrm{mm}$) is almost the same as that measured at the level of the grey matter. This means that the detected packets of photons are mostly propagated in the grey matter under these small blood vessels. This result is consistent with the results presented by Giannoni et al.²⁴ and indicates that large blood vessels have a large impact on the path length. We also observed that there were no significant differences between the mean path length measured on arteries and veins. Indeed, when comparing the absorption coefficients of the artery and the vein ($\overline{{\mu }_{a,\text{Artery}}-{\mu }_{a,\text{Vein}}}=-15{\mathrm{cm}}^{-1}$), we do not observe as great a difference as that between the artery and the grey matter ($\overline{{\mu }_{a,\text{Artery}}-{\mu }_{a,\text{Grey matter}}}=172{\mathrm{cm}}^{-1}$) or the vein and the grey matter ($\overline{{\mu }_{a,\text{Artery}}-{\mu }_{a,\text{Grey matter}}}=188{\mathrm{cm}}^{-1}$, with $\overline{{\mu }_a}$ the average absorption coefficient over the wavelength range 400 to 1000 nm).

The segmentation step in the digital instrument simulator could also be improved, see Fig. 1 and Sec. 2.1.2. Segmentation is based on morphological operations calculated in a mask delimiting the surgical window. The latter was created manually. This manual operation has been chosen because automatic thresholding, using Otsu’s method⁵³ for example, does not allow the cerebral cortex to be identified repeatedly. In some patients, bleeding around the edges of the surgical window or bloody compresses prevents the cerebral cortex from being correctly identified. To overcome this limitation, segmentation methods based on deep learning (U-Net)⁵⁴ could be used for segmenting the blood vessels and delineate the contour of the surgical window. For the moment, this approach has not been implemented because we do not have enough data to train the neural network. In this study, we modeled a perfect optical device. The real spectral characteristics of the light source and the sensor could be integrated into the digital simulator by normalizing the estimated quantities with spectral sensitivity curves of the light sources and the camera, as proposed in our previous study.²² We modeled a real sensor in terms of resolution and definition, but this could be improved with the conversion of the diffuse reflectance into sensor intensities by considering parameters of the optics systems⁵⁵ (lens magnification, lens transmission, quantum efficiency). We can also integrate the effect of the lenses on the image creation such as geometric distortion, vignetting, or chromatic aberrations. A lens design model can be used for this purpose.⁵⁶ Using this digital instrument simulator, we modeled a simple cerebral activation by changing the concentrations of the chromophores in a portion of activated grey matter (see Fig. 5). However, we do not model the complex physiologic events related to the neuro-vascular coupling. To answer this limitation, we plan to model the spatio-temporal concentration changes in arteriole, capillary, and venous compartments related to blood volume, flow velocity, and oxygen consumption with the implementation of the dynamic model proposed by Fantini.²³

Although the strength of this digital instrument is to use the white Monte Carlo approach, this method has also a limitation, which is the size of the generated files. Indeed, for one wavelength, 5 GB of data is generated, which leads to 305 GB for the 61 wavelengths modeled in our study. To reduce the memory size, the absorption coefficients need to be fixed in the modeled tissue to directly estimate the radiative quantities. However, with this approach, we cannot reconstruct the temporal changes of the absorption coefficient, so a simulation has to be performed per wavelength and per temporal index, which increases drastically the computation time.

4.2.

Identification of the Optimal Wavelength for Hemodynamic and Metabolic Monitoring

Using the data obtained with the digital instrument simulator, we proposed an optimization procedure based on the genetic algorithm to identify the best wavelength combinations in the visible and near-infrared range to quantify changes in ${\mathrm{HbO}}_2$, Hb, oxCCO, oxCytb, and oxCytc. We added noise to the simulated quantities to identify the wavelength that is the most robust to noise.

For hemodynamic monitoring, see Fig. 9, we can see that the optimal groups of 2, 4, 6, 8, or 10 wavelengths aim to reduce the quantification errors in $\mathrm{\Delta }{C}_{{\mathrm{HbO}}_2}$ and $\mathrm{\Delta }{C}_{\mathrm{Hb}}$ compared with that obtained with the configurations proposed in the literature. Contrary to the spectral configurations proposed by Bouchard et al.²⁹ and White et al.,³⁰ the wavelengths identified with the optimization procedure do not include the hemoglobin isobestic point at 530 nm. The wavelengths were located in the visible and near-infrared range, where ${\mathrm{HbO}}_2$ or Hb absorption predominates. This could be explained by the addition of noise in the optimization procedure, which makes it difficult to interpret the attenuation changes measured at the isobestic points. Moreover, the optimization procedure identifies the best wavelengths but does not take into account the bandwidth of the light source spectra in hyperspectral devices based on spectral-scanning technology.²⁸ Taking into account a spectral bandwidth could modify the value of the central wavelengths by a few nanometers but would not introduce changes greater than the spectral bandwidth.

The optimization procedure took profit of the peaks of oxidized and reduced cytochrome extinction spectra in the visible and the near-infrared range to monitor hemodynamic and metabolic changes. For hemodynamics and oxCCO monitoring, the optimal wavelengths aim to reduce the quantification errors compared with that obtained with broadband spectra between 780 and 900 nm. In this study, we also proposed to quantify the changes of the cytochrome $b$ and $c$. We can see in Fig. 11 that the addition of the cytochromes $b$ and $c$ in the modified Beer–Lambert system helps to better resolve oxCCO changes compared with the model based on three chromophores (${\mathrm{HbO}}_2$, Hb, and oxCCO), see Fig. 10. Indeed, we can see that the quantification errors in oxCCO are lower when all cytochromes are considered, see Table 3. This result is interesting as it may help to obtain robust devices to monitor brain metabolism, which could be a hyperspectral device for intraoperative brain mapping or even a NIRS device for bedside monitoring. The monitoring of the oxidation state of CCO could help to obtain a more direct biomarker of neuronal activity, to detect brain injuries,^16,57 and to better understand the neurovascular coupling.⁵⁸

The optimization procedure does not help to resolve changes in cytochrome $b$ and $c$ in a significant way. Indeed, we can see in Figs. 11 and 12 that oxCytb and oxCytc measurements are really noisy. A two-sample $T$-test tells us that changes measured at the level of the activated grey matter are not significantly different from those measured on non-activated grey matter ($p_{\text{values}}$ are mainly higher than 0.1). This is maybe due to the noise addition, coupled with the fact that the attenuation changes due to the oxidized cytochrome $b$ and $c$ are rather low compared with those of hemoglobin and oxCCO.

The optimization procedure helps to obtain smaller quantification errors compared with the literature configurations. For hemodynamic monitoring, increasing the number of wavelengths from 2 to 10 helps to reduce linearly the quantification error from 8% to 26%, see Fig. 9. However for hemodynamic and metabolic monitoring, we do not observe a significant decrease in the quantification error compared with Bale’s spectral configuration if using 6, 8, or 10 optimal wavelengths (quantification errors of 32%, 33%, and 33%, respectively, see Fig. 10). Thus, the number of wavelengths of a hyperspectral camera could be reduced to six to resolve hemodynamic and CCO changes. It could help to optimize several critical parameters for intraoperative functional brain mapping applications such as the signal-to-noise ratio. However, the availability of light sources or spectral filters at the identified wavelengths may limit the development of the hyperspectral camera. To address this limitation, we plan to restrict the optimization procedure to a number of wavelengths selected from the available light sources or spectral filters. For each wavelength, the spectral bandwidth can also be incorporated into the procedure. As proposed by Leadley et al,³² we will also investigate in a future study the influence of extinction coefficient uncertainties on chromophore quantification.

We also showed that a separate quantification of hemodynamics and oxCCO leads to a better estimation of the chromophore changes, see Table 3. The optimization procedure helps to identify the best wavelength combination in the range 400 to 1000 nm, which reduces the quantification errors in ${\mathrm{HbO}}_2$, Hb, and oxCCO by 47%, 57%, and 57% compared with the gold standard of 121 wavelengths between 780 and 900 nm. The configuration of six wavelengths is 468, 482, 610, 753, 814, and 907 nm to monitor ${\mathrm{HbO}}_2$ and Hb changes and 608, 650, 661, 672, 820, and 865 nm to monitor oxCCO changes. The separate quantification method should be taken with caution because it distorts the link between chromophores. Moreover, this could lead to incorrect interpretations between the chromophore changes and the physiological status of the patient. However, this method could be interesting for some clinical applications. For example, the method can be used to identify a spectral configuration for precise monitoring of Hb changes with the idea to define a blood oxygen level dependent-like contrast.⁵⁸ In intraoperative brain mapping studies, this approach has been proposed by several research groups. The authors used a single illumination at 605⁵⁹ or 610 nm⁶⁰ combined with a monochrome camera to assess Hb absorption. For these two wavelengths, the authors considered that the signal changes were mainly due to Hb because the absorption related to ${\mathrm{HbO}}_2$ is negligible compared with that of Hb (the ratio between Hb molar extinction coefficient to that of ${\mathrm{HbO}}_2$ is $\approx 6.26$). Although the contribution of ${\mathrm{HbO}}_2$ is small compared with that of Hb, the signal measured with this type of device cannot accurately measure variations in Hb concentration. Our method solves this problem.