1.

Introduction

Noninvasive functional brain mapping is an imaging technique used to localize the functional areas of the patient brain. This technique is used during brain tumor resection surgery to indicate to the neurosurgeon the cortical tissues which should not be removed without cognitive impairment. Functional magnetic resonance imaging (fMRI)¹ is widely used to localize patient functional areas. However, after patient craniotomy, a brain shift invalidates the relevance of neuronavigation to intraoperatively localize the functional areas of the patient brain.² In order to prevent any localization error, intraoperative MRI has been suggested but it complicates the surgery gesture which makes it rarely used. For these reasons, electrical brain stimulation (EBS)³ is the gold standard for the identification of brain functional areas. For instance, Roux et al.⁴ demonstrated that language fMRI data obtained with naming or verb generation tasks, before and after surgery, were imperfectly correlated with electrical brain mapping. The overall results of their study demonstrated that language fMRI could not be used to make critical surgical decisions in the absence of EBS.

In 1977, Jöbsis⁵ demonstrated that the blood and tissue oxygenation changes in the brain can be measured using near-infrared (NIR) light. Then Chance et al.⁶ demonstrated that optical imaging (visible and NIR light) can monitor the brain activity with the determination of concentration changes of oxygenated hemoglobin ($\mathrm{\Delta }[{\mathrm{HbO}}_2]$), deoxygenated hemoglobin ($\mathrm{\Delta }[\mathrm{Hb}]$), blood volume ($\mathrm{\Delta }[\mathrm{Hb}]+\mathrm{\Delta }[{\mathrm{HbO}}_2]$), and cytochrome oxidase. Since then, numerous works have demonstrated the ability of optical imaging to detect functional areas thanks to hemodynamics.^7–13 The motivation for this paper is derived from the necessity to analyze the hemodynamics in brain tissue following the neuronal activation, which are closely linked to the blood oxygenation level-dependent (BOLD) contrast used in functional MRI studies.

During neurosurgery, the craniotomy gives a direct access to the brain cortex. Intrinsic optical imaging^14–16 can be used intraoperatively to localize the patient hemodynamic activity in the cerebral cortex. The intrinsic signal refers to the cortical reflectance changes^14,17 due to hemodynamic response. A hyperspectral camera,^18,19 or a single-wavelength illumination in conjunction with a low-noise CCD camera,^16,20 can be used to acquire the intrinsic signal. The time course of this signal is characterized by the early hemodynamic responses in brain tissue related to neuronal activity (initial dip) followed by a larger response that corresponds to the BOLD signal in fMRI studies.¹⁷ This technique is a powerful tool to understand the cognitive functions at the neural circuit level¹⁵ and to define more precisely the hemodynamic response following a physiological stimulus.^21,22 In some studies, a spectroscopic analysis of the intrinsic signal is computed to assess the cortical hemoglobin concentration changes. New approaches consist in using an RGB^21,23 or hyperspectral¹⁸ camera with a continuous wave white light source^18,21,24 or pulsed narrow bandpass illumination sources.^23,25 These setups have the main advantages of being usable in real time^18,25 and directly in the operative room. Steimers et al.²³ analyzed an exposed rat cortex with an RGB-LED light source and an RGB camera. The results of their work indicate that semiquantitative functional maps (in arbirary units) can be processed with the modified Beer–Lambert law. But the optical mean optical path lengths were not taken into consideration. Bouchard et al.²⁵ developed an ultrafast device made up of two-pulsed LEDs and a monochromatic camera to assess in real time the hemoglobin concentration changes in blood vessels using the modified Beer–Lambert law. Inhomogeneities of the optical properties were not taken into consideration since an averaged path length provided by Monte-Carlo simulation was used in the model. Pichette et al.¹⁸ used a hyperspectral camera and continuous wave white light illumination to assess in real time the absorbance changes. In their study, the blood vessels and the gray matter pixels were automatically segmented by comparing measured reflectance spectra to blood and gray matter simulated reflectance spectra. This allowed the blood vessels to be masked to reduce their influence on absorbance change measurements. In these works, qualitative or quantitative brain maps are not compared to the patient hemodynamic response.²⁶ So it remains difficult for the neurosurgeon to efficiently localize the functional brain areas with an RGB camera.

The objective of the present work is to supply the methodological tools for the construction of quantitative functional brain maps based on the analysis of the cortical reflectance spectra acquired by a digital RGB camera. For each camera pixel, the intrinsic signal was converted into oxy- and deoxygenated hemoglobin concentration changes using the modified Beer–Lambert law. This law was computed with estimated mean optical path lengths calculated by Monte-Carlo simulations. These mean optical path lengths are chosen according to the local optical properties of the patient brain. Three Monte-Carlo models have been investigated for the study of the light propagation in cortical tissue such as gray matter, surface blood vessel, and buried blood vessel. We propose in our study to complement the intrinsic optical imaging approaches with statistical analyses inspired by the BOLD fMRI method. The Pearson correlation coefficient is calculated between the expected hemodynamic response²⁶ and each measured concentration changes time courses to accurately localize the functional areas of the patient brain. The brain areas identified by intraoperative electrical stimulation showed a good correlation with the quantitative brain maps processed with the proposed method. These results could help to get robust intraoperative brain area identification based on RGB imaging.

2.

Materials and Methods

A schematic overview of the computation of the functional quantitative maps can be seen in Fig. 1. Once the video was acquired, the following computational steps were applied. The first image was manually segmented into three classes (gray matter, surface blood vessel, and buried blood vessel). This segmentation step aims to associate each camera pixel to the appropriate optical mean path length used in the Beer–Lambert law; see Sec. 2.4. For each frame of the video, the brain repetitive motion was compensated, then data were corrected and filtered. The details of each step are given in Sec. 2.3. A functional model was applied to preprocessed data to compute functional quantitative maps; see Sec. 2.4.

Fig. 1

Overview of the method.

2.1.

Experimental Setup

The imaging system is composed of an RGB CMOS camera (BASLER acA2000-165uc) in conjunction with an Edmund Optics camera lens ($f=50\mathrm{mm}f/2-f/22$), a continuous wave white light source (OSRAM Classic 116-W 230-V light bulb) and a laptop (processor: Intel Core i5-7200U, $2.50\mathrm{GHz}\times 4$, RAM: 15.3 GiB); see Fig. 2. During data acquisition, the camera also acquired residual light since the operative room lights were on. Data were directly acquired by the laptop via an USB link. A C++ software acquired and processed the images using open source tools such as Qt (v5.9.4), openCV (v3.2.0),²⁷ FFTW (v3.3.7),²⁸ and pylon (BASLER library). 8 bits RGB images were acquired every 33 ms (the sampling rate is set to 30 frames per second) with a resolution which at best is 400 ppi (the minimum size of a square pixel is $64\times 64\mu \mathrm{m}$).

Fig. 2

Schematic of the imaging system.

2.3.

Preprocessing

2.3.1.

Manual image segmentation

The first image of the video sequence was segmented with a semiautomated procedure into three classes: gray matter, surface blood vessel, and buried blood vessel. Pixels were clustered into four clusters using the $K$-means algorithm from the C++ library OpenCV (v3.2.0).²⁷ The components of each cluster were manually sorted and attributed to the three classes. The cluster components that have not been selected corresponded to saturated areas (specular reflection) and were discarded. As the video is motion compensated, the segmentation of the first image is valid during the whole video sequence. The segmentation method is able to detect blood vessels with a diameter larger than $500\mu \mathrm{m}$. The blood vessels with a diameter smaller than $500\mu \mathrm{m}$ are included in the gray matter class. This segmentation step aims to associate each camera pixel to the appropriate optical mean path length used in the modified Beer–Lambert law; see Sec. 2.4.

2.3.2.

Motion compensation

After craniotomy, the brain surface undergoes a repetitive motion due to the patient breath and cardiac pulsation. This motion as well as potential video camera motion prevents accurate video analysis. The motion compensation aims to ensure that each camera pixel corresponds to the same cortical area all along data acquisition. Sdika et al.^29,30 proposed a repetitive motion compensation algorithm. This algorithm is split into two parts. First, the basis of the repetitive motion is learned from few initial frames of the video (50 frames), then, each video frame is registered to the reference image (first image of the video). To ensure that motion compensation worked in normal operation, each registered image is compared to the reference image using the normalized cross covariance (NCC):

Eq. (1)

$$\mathrm{NCC}(t)=\frac{n\sum _{}^{}{I}_{0}I(t)-\sum _{}^{}{I}_0\sum _{}^{}I(t)}{\sqrt{[n\sum _{}^{}{I}_0^2-{(\sum _{}^{}{I}_0)}^2][n\sum _{}^{}I{(t)}^2-{(\sum _{}^{}I(t))}^2]}}.$$

$I(t)$ is the image registered at time $t$ and $I_0$ the reference image (first frame). $n$ denotes the number of pixels of the images. The more the NCC value tends toward 1, the more the images $I(t)$ and $I_0$ are similar. The NCC curve can be compared with six reference NCC curves that have been computed with six registered videos, validated by Sdika et al.’s algorithm.^29,30 A linear regression is applied on each reference NCC curve to get their slope and intercept. The reference videos do not have the same recorded duration or the same sampling rate than the acquired video. The linear regression aims to rebuild the reference NCC curves in the temporal domain of the acquired video.

Let ${\mathrm{NCC}}_{\mathrm{ref}}^i$ be the reference NCC curve $i$ [$i\in (\mathrm{1,6})$] and ${\mathrm{NCC}}_{\mathrm{mes}}^j$ the NCC curve of the registered video $j$ [$j\in (\mathrm{1,5})$]. The acceptable NCC variation range is defined as the area between the horizontal lines $y_{{\mathrm{ref}}_1}=1$ and $y_{{\mathrm{ref}}_2}=\min [{\mu }_{\mathrm{ref}}(t)-{\sigma }_{\mathrm{ref}}(t)]$, with ${\mu }_{\mathrm{ref}}(t)=\frac{\sum _{i=1}^6{\mathrm{NCC}}_{\mathrm{ref}}^i(t)}{6}$ and ${\sigma }_{\mathrm{ref}}(t)=\sqrt{\frac{\sum _{i=1}^6{[{\mathrm{NCC}}_{\mathrm{ref}}^i(t)]}^2}{6}-{[{\mu }_{\mathrm{ref}}(t)]}^2}$. The NCC dispersion range of the registered videos is defined as the area between the curves $y_1(t)={\mu }_{\mathrm{mes}}(t)-{\sigma }_{\mathrm{mes}}(t)$ and $y_2(t)={\mu }_{\mathrm{mes}}(t)+{\sigma }_{\mathrm{mes}}(t)$. The motion compensation of the five videos is validated if the NCC dispersion range of the registered videos is included in the acceptable NCC variation range.

2.3.3.

Data filtering

According to Chance et al.,⁶ low-frequency modulation of light absorption is linked to cortical activity. These frequencies of interest can be visualized in the expected hemodynamic response spectrum. The hemodynamic response of the brain in relation to the neural activities (cortical stimulus) can be expressed by the hemodynamic impulse response function (HIRF).²⁶ In our application, the expected hemodynamic response can be obtained by convolving HIRF with the window function $P$ representing the experimental paradigm; see Fig. 3. The expected hemodynamic response is the ideal representation of the hemodynamic response of the patient whose validity is discussed in Sec. 4. A low-pass filtering was implemented by multliplying the signal by a Blackman window (cut off frequency: 0.05 Hz; see Fig. 3) in the Fourier domain.

Fig. 3

(a) HIRF.²⁶ (b) The green curve represents the experimental paradigm $P$ and the black curve the expected hemodynamic response which is obtained by convolving HIRF with $P$. (c) The black curve represents the expected heamodynamic response spectrum and the red curve, the transfer function of a Blackman window (cut off frequency: 0.05 Hz) in the Fourier domain.

2.3.4.

Data correction

According to Oelschlägel et al.,²⁴ data preprocessing is mandatory to correct the slow drift of the collected intensity due to the cortical tissue desiccation during the video acquisition. The assumption is that the beginning and the end of the video corresponded to the same patient physiological state. So we considered that these intensity values had to be identical. Linear regressions were computed on the measured RGB time courses, then the calculated regression lines was subtracted to original data to compensate the slow drift of the collected intensity.

2.4.

Functional Model

The standard approach for the analysis of the reflected spectra is based on the modified Beer–Lambert law. The assessment of the concentration changes depends on the determination of the mean optical path length of the detected photons; see Eq. (3). For the purpose of measuring more precisely the hemoglobin concentration changes, three different cortical tissues have been modeled; see Fig. 4. Monte-Carlo simulations were processed to estimate their mean optical path length which were assigned to pixels belonging to the segmented classes; see Sec. 2.3.1.

Fig. 4

Representation of modeled cortical tissues. Volumes are made up of $50\times 50\times 50$ voxels with a $1\text{-}{\mathrm{mm}}^3$ resolution. Red voxels represent large blood vessels and gray voxels cortical tissues. A black arrow symbolizes a Monte-Carlo simulation for an emission of ${10}^6$ packets of photons at a given position. Model 1 represents a 2-mm diameter blood vessel on the surface of the cortical tissue, model 2 a 2-mm diameter blood vessel buried under 1 mm of gray matter, and model 3 a cortical tissue without a large blood vessel.

2.4.1.

Modified Beer–Lambert law

The modified Beer–Lambert law can be expressed as a matrix system:³¹

Eq. (2)

$$\left[\begin{array}{c}\mathrm{\Delta }{A}_R(t)\\ \mathrm{\Delta }{A}_G(t)\\ \mathrm{\Delta }{A}_B(t)\end{array}\right]=\left[\begin{array}{cc}{E}_{\mathrm{R},\mathrm{Hb}}& E_{{\mathrm{R},\mathrm{HbO}}_2}\\ E_{\mathrm{G},\mathrm{Hb}}& E_{{\mathrm{G},\mathrm{HbO}}_2}\\ E_{\mathrm{B},\mathrm{Hb}}& E_{{\mathrm{B},\mathrm{HbO}}_2}\end{array}\right]\times \left[\begin{array}{c}\mathrm{\Delta }[\mathrm{Hb}](t)\\ \mathrm{\Delta }[{\mathrm{HbO}}_2](t)\end{array}\right],$$

where

Eq. (3)

$$E_{i,n}=\int {ϵ}_n(\lambda )D_i(\lambda )S(\lambda )L(\lambda )\mathrm{d}\lambda .$$

$\mathrm{\Delta }{A}_i(t)$ is the absorbance changes measured at time $t$:

Eq. (4)

$$\mathrm{\Delta }{A}_i(t)={\log }_{10}\left(\frac{R_i^0}{R_i(t)}\right),$$

where $R_i(t)$ is the reflectance intensity measured at time $t$ by the camera color channel $i$ (red, green, or blue), $R_i^0$ is the reference reflectance intensity measured by the camera color channel $i$ (average of the reflectance intensity over the duration of the first rest step of the experimental paradigm; see Sec. 2.2). ${\epsilon }_n$ is the extinction coefficient of the chromophore $n$ (in $L{\mathrm{mol}}^{-1}{\mathrm{cm}}^{-1}$). $\mathrm{\Delta }[\mathrm{Hb}]$ is the deoxygenated hemoglobin molar concentration changes (in $\mathrm{mol}{\mathrm{L}}^{-1}$) and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ the oxygenated hemoglobin molar concentration changes (in $\mathrm{Mol}{\mathrm{L}}^{-1}$). Our model takes into consideration the receiving spectrum of the RGB camera and the emission spectrum of the light source. The spectral sensitivity of the detector $i$ of the RGB camera is represented by $D_i(\lambda )$ and $S(\lambda )$ is the normalized intensity spectrum of the light source. $L(\lambda )$ is the wavelength-dependent mean optical path length of the photons traveled in tissue. Hemoglobin concentration changes are obtained by matrix inversion once the matrix $E$ has been calculated. However, the mean optical path length $L(\lambda )$ needs to be determined.

2.4.2.

Mean optical path length determination

The mean optical path length $L(\lambda )$ was obtained by Monte-Carlo simulations using MCX software.³² Three cortical tissues were modeled under a homogeneous white light illumination; see Fig. 4. Model 1 represents a 2-mm diameter blood vessel on the surface of the cortical tissue, model 2, a 2-mm diameter blood vessel buried under 1-mm of cortical tissue, and model 3, a cortical tissue without a large blood vessel. The cortical tissue is composed of gray matter perfused by capillaries.

Each voxel of the modeled tissues included the information of optical parameters (absorption, reduced scattering, anisotropy coefficients, and refractive index). A homogeneous illumination was achieved by scanning the emission point of the light source along the entire model surface. A white light illumination has been simulated by scanning the optical parameters along the entire illumination spectrum (from 400 to 1000 nm in steps of 10 nm). A total of $50\times 50\times 61\times 3$ simulations were computed. $50\times 50$ represents the number of simulations used to illuminate the modeled cortical tissues in a homogeneous way. 61 denotes the number of emitted wavelength and 3 the number of modeled cortical tissues. The optical parameters were taken from the literature and correspond to a nominal physiological condition. References 33–35 have been used for the blood vessel parameters, and Refs. 35–38 for the gray matter parameters.

The size of the modeled tissues has been chosen in accordance with the photon sensitivity profile³⁹ computed for the detector situated at the center of the top face of model 3. A detector collects all the photons reaching the surface in a square area of $1{\mathrm{mm}}^2$. Model 3 has been chosen because the photons traveling through the cortical tissue without a large blood vessel have a greater probability to have a long trajectory than the photons traveling through the cortical tissue with a large blood vessel. The photon sensitivity profile is roughly represented as a half sphere of radius of 12.5 mm. To avoid any photon loss and inexact results due to the boundary conditions (a simulation of the travel of a packet of photon stops when this packet of photon leaves the volume), the size of the models is set to $50\times 50\times 50$ voxels with a resolution of $1{\mathrm{mm}}^3$. This allows us to precisely compute the mean optical path length of the photons reaching the detectors located at the top face of the volume.

For the three models, the mean optical path length $L$ of all detected photons is calculated; see Fig. 5. For $\lambda <580\mathrm{nm}$, the mean path length of model 1 is on average 39% smaller than the one of model 3. The one of model 2 is on average 4% smaller than the one of model 3. For $\lambda >580\mathrm{nm}$, the mean path length of model 1 is on average 42% smaller than the one of model 3. The one of model 2 is on average 34% smaller than the one of model 3. When the blood vessel is buried under 1 mm of gray matter (model 2), a smaller proportion of photons is absorbed in model 2 than in model 1. However, this blood vessel still has a large impact on the photon propagation.

Fig. 5

The red, blue, and gray curves represent the computed wavelength dependent mean optical path length of models 1, 2, and 3, respectively (see Fig. 4).

2.4.3.

Quantitative functional brain map

Two quantitative functional brain maps were computed by selecting the $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ time courses which were correlated with the expected hemodynamic response. The Pearson correlation measurement aims to highlight the hemodynamic response associated with the experimental paradigm and ignore any other hemodynamic variations. Only positive correlations were considered, so the $\mathrm{\Delta }[\mathrm{Hb}]$ time courses were multiplied by $-1$. Two quantitative functional brain maps were defined as the Hb and ${\mathrm{HbO}}_2$ concentration changes averaged over the duration of the patient hand stimulation (see Sec. 2.2). The quantitative maps were only processed for pixels for which the Pearson correlation coefficient was higher than the Pearson correlation coefficient threshold (PCCT) value:

Eq. (5)

$${\mathrm{QMap}}_C(x,y)=\{\begin{array}{ll}\overline{C}(x,y),& \mathrm{if}{r}_C(x,y)\ge \mathrm{PCCT}\\ \text{nonprocessed},& \text{otherwise}.\end{array},$$

where

Eq. (6)

$$\overline{C}(x,y)=\frac{\sum _{n=N_1}^{N_2}\mathrm{\Delta }[C](n,x,y)}{N_2-N_1}.$$

${\mathrm{QMap}}_C$ is the chromophore $C$ quantitative functional brain map. $C$ represents either ${\mathrm{HbO}}_2$ or Hb. $(x,y)$ denotes an image pixel position. $n$ is the frame index, $N_1$ and $N_2$ are the patient hand stimulation starting and ending frame indexes. $r_C(x,y)$ denotes the Pearson correlation coefficient calculated between the chromophore $C$ concentration changes time course at the image pixel position $(x,y)$ and the expected hemodynamic response. $\bar{C}$ represents the $\mathrm{\Delta }[C]$ value averaged over the patient activity period. A $5\times 5$ median filter is applied to the quantitative functional maps. This induced a decrease of the resolution of the quantitative functional map from 400 to 80 ppi.

2.4.4.

Definition of activated cortical areas

In our study, an activated cortical area is a functional area that is associated with the patient hand stimulation; see Sec. 2.2. The $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ time courses of an activated cortical area should be highly correlated with the expected hemodynamic response, whereas the ones of a nonactivated cortical area should not. The $\overline{\mathrm{Hb}}$ and $\overline{{\mathrm{HbO}}_2}$ values [see Eq. (6)] of an activated cortical area should be significantly different from the values of a nonactivated cortical areas. Based on these assumptions, an activated cortical area can be statistically defined. Let ${\mathrm{Cort}}_{\mathrm{NF}}$ be a cortical area which has not been identified as functional by the EBS, and ${\mathrm{Cort}}_{\mathrm{F}}$ a functional area of interest identified by EBS. The ${\mathrm{Cort}}_{\mathrm{NF}}$ and ${\mathrm{Cort}}_{\mathrm{F}}$ areas are defined according to four quantitative distributions [$\overline{\mathrm{Hb}}$, $\overline{{\mathrm{HbO}}_2}$, $r_{\mathrm{Hb}}$, and $r_{{\mathrm{HbO}}_2}$; see Eqs. (5) and (6)]. Welch’s $T$ tests were computed to score the differences between the mean values of the ${\mathrm{Cort}}_{\mathrm{NF}}$ and ${\mathrm{Cort}}_{\mathrm{F}}$ areas. These means were denoted $⟨\overline{\mathrm{Hb}}⟩$, $⟨\overline{{\mathrm{HbO}}_2}⟩$, $⟨r_{\mathrm{Hb}}⟩$, and $⟨r_{{\mathrm{HbO}}_2}⟩$. If all four tests reject the null hypothesis at 1% significance level, the ${\mathrm{Cort}}_{\mathrm{F}}$ area is defined as an activated cortical area.

3.

Results

3.1.

Motion Compensation

In Fig. 6, the blue rectangular area corresponds to the acceptable NCC variation range defined as the area between the horizontal lines $y_{{\mathrm{ref}}_1}=1$ and $y_{{\mathrm{ref}}_2}=0.995$ (see Sec. 2.3.2). The red area corresponds to the NCC dispersion range of the five unregistered videos (see Sec. 2.3.2). The green area corresponds to the NCC dispersion range of the five registered videos. The NCC dispersion range of the five registered videos is included in the acceptable NCC variation range, whereas the NCC dispersion range of the five unregistered videos is not. This validates the motion compensation of the five videos.

Fig. 6

Validation of the motion compensation. The blue area corresponds to the acceptable NCC variation range. The red area corresponds to the NCC dispersion range of the five unregistered videos and the green area corresponds to the NCC dispersion range of the five registered videos.

3.2.

Quantitative Functional Maps

The Hb and HbO₂ quantitative functional maps of videos 1–5 are represented in Fig. 7. The colorbar represents the scale of variation of the ${\mathrm{QMap}}_{\mathrm{Hb}}$ and ${\mathrm{QMap}}_{{\mathrm{HbO}}_2}$ values in $\mu \mathrm{mol}{\mathrm{L}}^{-1}$; see Eq. (5). $M_{i-j}$ designates the motor area $i$ of the patient $j$ identified by EBS; see Sec. 2.2. $S_{i-j}$ designates the sensory areas $i$ of the patient $j$ identified by EBS; see Sec. 2.2. $C_j$ designates a cortical area of the patient $j$ which has not been identified as a functional area by the EBS. These areas are delimited with white circles of 10 pixels diameter for video 1 and 30 pixels diameter for videos 2 to 5. For each video, Hb and ${\mathrm{HbO}}_2$ quantitative functional maps are plotted according to four PCCT values (0, 0.3, 0.5, and 1). Only the $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ time courses whose Pearson correlation coefficient is higher than the PCCT value are considered in Eq. (5). The highest ${\mathrm{QMap}}_{\mathrm{Hb}}$ and ${\mathrm{QMap}}_{{\mathrm{HbO}}_2}$ values are situated at the level of the blood vessels surrounding the motor and sensory areas.

Fig. 7

Hb and ${\mathrm{HbO}}_2$ functional maps computed for the five videos. For videos 1, 2, 3, and 5, the stimulation of the cortex was achieved through a repetitive and alternative hand opening and closing at $\approx 1\mathrm{Hz}$ (movement performed by the patient: videos 1, 2, and 5; movement induced by an external person: video 3). For video 4, the stimulation of the cortex was achieved through a repetitive fingers and palm caresses at $\approx 1\hspace{0.17em}\mathrm{Hz}$ (the caresses were performed by an external person). The Hb and ${\mathrm{HbO}}_2$ functional maps are computed for different PCCT values. The colorbar represents the scale of variation of the ${\mathrm{QMap}}_{\mathrm{Hb}}$ and ${\mathrm{QMap}}_{{\mathrm{HbO}}_2}$ values in $\mu \mathrm{mol}{\mathrm{L}}^{-1}$ [see Eq. (5)]. $M_{i-j}$ designates the motor area $i$ of the patient $j$ identified by EBS. $S_{i-j}$ designates the sensory areas $i$ of the patient $j$ identified by EBS. $C_j$ designates a nonactivated area of the cortex of the patient $j$. A $5\times 5$ median filter is applied to the functional maps.

3.2.1.

Video 1

The stimulation of the cortex was achieved through a repetitive and alternative hand opening and closing at $\approx 1\mathrm{Hz}$ (the movement was performed by the patient). For $\mathrm{PCCT}=0$, all QMap values except the ones associated with the saturated pixels are processed. For $\mathrm{PCCT}=0.3$, the QMap values situated at the level of the $M_{1-1}$, $M_{2-1}$, $M_{3-1}$, $S_{1-1}$, and $S_{2-1}$ areas are computed. Some QMap values are also processed at the left side of the image. For $\mathrm{PCCT}=0.5$, the QMap values situated at the level of the $M_{1-1}$, $M_{2-1}$, $M_{3-1}$, $S_{1-1}$, and $S_{2-1}$ areas are computed. Some ${\mathrm{QMap}}_{{\mathrm{HbO}}_2}$ values are processed at the left side of the image.

We defined five points of interest within video 1 to explore in more details the hemodynamic time courses in the different brain areas. In Fig. 8(a), the point $M_1$ is situated at the center of the $M_{1-1}$ area. The point $S_1$ is situated at the center of the $S_{1-1}$ area. The point $A$ is situated on the blood vessel above the $M_{1-1}$ area. The point $B$ is localized on a blood vessel which is not linked with the functional areas. The point $C$ represents a point gray matter situated far away from the functional areas. The $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ time courses of these points of interest can be compared in Fig. 8(b).

Fig. 8

(a) Localization of points of interest. (b) Hb and ${\mathrm{HbO}}_2$ concentration changes time courses of the points of interest defined in (a).

From 0 to 35 s, the $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ values of the point $M_1$ oscillate above and below $0\mu \mathrm{mol}{\mathrm{L}}^{-1}$. From 35 to 56 s, the $\mathrm{\Delta }[\mathrm{Hb}]$ values decrease from 0 to $-1.7\mu \mathrm{mol}\mathrm{L}{.}^{-1}$ and the $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ values increase from 0 to $3\mu \mathrm{mol}{\mathrm{L}}^{-1}$. From 56 to 68 s, the $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ values revert to $0\mu \mathrm{mol}{\mathrm{L}}^{-1}$. From 68 s to the end of the acquisition, the $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ values oscillate again above and below $0\mu \mathrm{mol}{\mathrm{L}}^{-1}$. From 0 to 35 s, the $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ values of the point $S_1$ oscillate above and below $0\mu \mathrm{mol}{\mathrm{L}}^{-1}$. From 35 to 45 s, the $\mathrm{\Delta }[\mathrm{Hb}]$ values decrease from 0 to $-1.0\mu \mathrm{mol}{\mathrm{L}}^{-1}$ and the $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ values increase from 0 to $1.9\mu \mathrm{mol}{\mathrm{L}}^{-1}$. From 45 to 68 s, the $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ values revert to $0\mu \mathrm{mol}{\mathrm{L}}^{-1}$. From 68 s to the end of the acquisition, the $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ values oscillate again above and below $0\mu \mathrm{mol}\hspace{0.17em}{\mathrm{L}}^{-1}$. From 0 to 35 s, the $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ values of the point $A$ oscillate above and below $0\mu \mathrm{mol}{\mathrm{L}}^{-1}$. From 34 to 43 s, the $\mathrm{\Delta }[\mathrm{Hb}]$ values decrease from 0 to $-7.4\hspace{0.17em}\mu \mathrm{mol}{\mathrm{L}}^{-1}$ and the $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ values increase from 0 to $14.3\mu \mathrm{mol}{\mathrm{L}}^{-1}$. From 43 to 68 s, the $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ values revert to $0\mu \mathrm{mol}{\mathrm{L}}^{-1}$. From 68 s to the end of the acquisition, the $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ values oscillate again above and below $0\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ curves of the points $B$ and $C$ oscillate without any understandable correlation with the experimental paradigm. The $\mathrm{\Delta }[\mathrm{Hb}]$ variations of the point $B$ are $\approx 3.6$ times greater than the ones of the curves $C$. The $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ variations of the point $B$ are $\approx 5.3$ times greater than the ones of the curves $C$.

The $\overline{\mathrm{Hb}}$, $\overline{{\mathrm{HbO}}_2}$, $r_{\mathrm{Hb}}$, and $r_{{\mathrm{HbO}}_2}$ values [see Eq. (5)] of the curves plotted in Fig. 8(b) are represented in Table 2. The points $M_1$ and $A$ have the highest $\overline{\mathrm{Hb}}$ and $\overline{{\mathrm{HbO}}_2}$ values ($\overline{\mathrm{Hb}}\le -1.16\mu \mathrm{mol}{\mathrm{L}}^{-1}$ and $\overline{{\mathrm{HbO}}_2}\ge 1.82\mu \mathrm{mol}{\mathrm{L}}^{-1}$). The $\overline{\mathrm{Hb}}$ and $\overline{{\mathrm{HbO}}_2}$ values of the points $S_1$ and $B$ are in the same order of magnitude. The point $C$ has the lowest $\overline{\mathrm{Hb}}$ and $\overline{{\mathrm{HbO}}_2}$ values. The points $M_1$, $S_1$, and $A$ have higher $r_{\mathrm{Hb}}$ and $r_{{\mathrm{HbO}}_2}$ values ($r_{\mathrm{Hb}}\ge 0.76$ and $r_{{\mathrm{HbO}}_2}\ge 0.68$) than $B$ and $C$ points ($r_{\mathrm{Hb}}\le 0.04$ and $r_{{\mathrm{HbO}}_2}\le 0.06$).

Table 2

Hb‾, HbO2‾ values (in μmol L−1) and rHb, rHbO2 values of the curves plotted in Fig. 8(b).

	Hb¯(μmol  L−1)	rHb	HbO2¯(μmol  L−1)	rHbO2
$M_1$	$-1.16$	0.76	1.82	0.68
$S_1$	$-0.60$	0.84	1.22	0.83
$A$	$-4.08$	0.78	7.8	0.79
$B$	$-0.79$	0.02	1.07	0.06
$C$	$-0.07$	0.04	0.10	0.08

3.2.2.

Video 2

The stimulation of the cortex was achieved through a repetitive and alternative hand opening and closing at $\approx 1\mathrm{Hz}$ (the movement was performed by the patient). In Fig. 7, for $\mathrm{PCCT}=0$, all QMap values except the ones associated with the saturated pixels are processed. High QMap values are computed at the level of the blood vessels situated at the left and right side of the images. For $\mathrm{PCCT}=0.3$, the QMap values situated at the level of the $M_{1-2}$, $S_{1-2}$, $S_{2-2}$, and $S_{3-2}$ areas are computed. Some QMap values are also processed at the right side of the image. For $\mathrm{PCCT}=0.5$, the QMap values situated at the level of the $S_{1-2}$ and $S_{2-2}$ areas are computed. Some ${\mathrm{QMap}}_{{\mathrm{HbO}}_2}$ values are processed at the level of the $M_{1-2}$ area.

3.2.3.

Video 3

The stimulation of the cortex was achieved through a repetitive and alternative hand opening and closing at $\approx 1\mathrm{Hz}$ (the movement was induced by an external person). In Fig. 7, for $\mathrm{PCCT}=0$, all QMap values except the ones associated with the saturated pixels are processed. High QMap values are computed at the level of the blood vessels situated at the left and right side of the images. For $\mathrm{PCCT}=0.3$, the QMap values situated at the level of the $M_{1-2}$, $S_{1-2}$, and $S_{3-2}$ areas are computed. Some QMap values are also processed at the right side of the image. For $\mathrm{PCCT}=0.5$, the QMap values situated at the level of the $M_{1-2}$ and $S_{1-2}$ areas are computed.

3.2.4.

Video 4

The stimulation of the cortex was achieved through a repetitive fingers and palm caresses at $\approx 1\mathrm{Hz}$ (the movement was induced by an external person). In Fig. 7, for $\mathrm{PCCT}=0$, all QMap values except the ones associated with the saturated pixels are processed. High QMap values are computed at the level of the blood vessels situated at the left and right side of the images. For $\mathrm{PCCT}=0.3$, the QMap values situated at the level of the $S_{1-2}$ and $S_{2-2}$ areas are computed. The QMap values are also processed at the level of the blood vessel drawing the junction between the motor and the sensory areas. Some QMap values are processed at the center of the images. For $\mathrm{PCCT}=0.5$, some QMap values situated at the level of the $S_{1-2}$ area are computed.

3.2.5.

Video 5

The stimulation of the cortex was achieved through a repetitive and alternative hand opening and closing at $\approx 1\mathrm{Hz}$ (the movement was performed by the patient). In Fig. 7, for $\mathrm{PCCT}=0$, all QMap values except the ones associated with the saturated pixels are processed. For $\mathrm{PCCT}=0.3$, the QMap values situated at the level of the $M_{1-3}$ and $M_{2-3}$ areas are computed. The QMap values are also processed at the level of the blood vessels surrounding these areas. QMap values are processed at the center and at the top of the images. For $\mathrm{PCCT}=0.5$, some QMap values situated at the level of the $M_{1-3}$ area are computed.

3.3.

Statistical Analysis

The functional areas of the patient exposed cortex which have been identified by EBS are compared to a reference cortical area $C$. This area $C$ has not been identified as functional by the EBS. The objective is to compute the activation criterion introduced in Sec. 2.4.4 on each functional area to evaluate our method sensitivity. For each video, an ANOVA and Kruskal–Wallis $H$-test have been separately computed on the $\overline{\mathrm{Hb}}$, $\overline{{\mathrm{HbO}}_2}$, $r_{\mathrm{Hb}}$, and $r_{{\mathrm{HbO}}_2}$ distributions of the areas defined in Fig. 7. These tests reject the null hypothesis that two or more distributions have the same population mean. Each video is studied separately. For each distribution of each video, a Welch’s $T$-test is computed between the sample $X_{i-j}$ and $C_j$ ($X$ represents either the motor $M$ or sensory $S$ areas $i$ of the patient $j$ defined in Fig. 7) with the null hypothesis that the two samples have identical mean values. In Figs. 9 and 10, the distribution of the $\overline{\mathrm{Hb}}$ and $\overline{{\mathrm{HbO}}_2}$ values of the areas defined in Fig. 7 are represented. In Figs. 11 and 12, the distribution of the $r_{\mathrm{Hb}}$ and $r_{{\mathrm{HbO}}_2}$ values of the areas defined in Fig. 7 are represented. In these graphs, the black diamonds represent the mean values of the distributions and the half length of the blue lines the standard deviation values. The notation * indicates that the null hypothesis is rejected with the respect of the Bonferroni correction ($p$ value$<0.01/N$, where $N$ is the number of samples per video). An area $X_{i-j}$ is defined as activated if the null hypothesis is rejected for the $\overline{\mathrm{Hb}}$, $\overline{{\mathrm{HbO}}_2}$, $r_{\mathrm{Hb}}$, and $r_{{\mathrm{HbO}}_2}$ distributions; see Sec. 2.4.4.

Fig. 9

Distribution of the $\overline{\mathrm{Hb}}$ values of the areas defined in Fig. 7. The black diamonds represent the mean values of the distributions ($⟨\overline{\mathrm{Hb}}⟩$) and the half length of the blue lines the standard deviation values. Each video is studied separately. The notation $i^{*}$ indicates a $T$-test’s statistical significance for the comparison of the means of the distribution $i$ and $C_j$ ($j$ represents the patient id).

Fig. 10

Distribution of the $\overline{{\mathrm{HbO}}_2}$ values of the areas defined in Fig. 7. The black diamonds represent the mean values of the distributions ($⟨\overline{{\mathrm{HbO}}_2}⟩$) and the half length of the blue lines the standard deviation values. Each video is studied separately. The notation $i^{*}$ indicates a $T$-test’s statistical significance for the comparison of the means of the distribution $i$ and $C_j$ ($j$ represents the patient id).

Fig. 11

Distribution of the $r_{\mathrm{Hb}}$ values of the areas defined in Fig. 7. The black diamonds represent the mean values of the distributions ($⟨r_{\mathrm{Hb}}⟩$) and the half length of the blue lines the standard deviation values. Each video is studied separately. The notation $i^{*}$ indicates a $T$-test’s statistical significance for the comparison of the means of the distribution $i$ and $C_j$ ($j$ represents the patient id).

Fig. 12

Distribution of the $r_{{\mathrm{HbO}}_2}$ values of the areas defined in Fig. 7. The black diamonds represent the mean values of the distributions ($⟨r_{{\mathrm{HbO}}_2}⟩$) and the half length of the blue lines the standard deviation values. Each video is studied separately. The notation $i^{*}$ indicates a $T$-test’s statistical significance for the comparison of the means of the distribution $i$ and $C_j$ ($j$ represents the patient id).

In Fig. 9, for video 1, the $⟨\overline{\mathrm{Hb}}⟩$ values of the motor and sensory areas are smaller than or equal to $-0.46\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{\mathrm{Hb}}⟩$ value of the $C_1$ area is equal to $0\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{\mathrm{Hb}}⟩$ values of the motor and sensory areas are statistically different from the one of the $C_1$ area. For video 2, the $⟨\overline{\mathrm{Hb}}⟩$ values of the motor and sensory areas are smaller than or equal to $-0.71\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{\mathrm{Hb}}⟩$ value of the $C_2$ area is equal to $0.10\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{\mathrm{Hb}}⟩$ values of the motor and sensory areas are statistically different from the one of the $C_2$ area. For video 3, the $\overline{⟨\mathrm{Hb}⟩}$ values of the $M_{1-2}$, $S_{1-2}$, and $S_{3-2}$ areas are smaller than or equal to $-0.40\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{\mathrm{Hb}}⟩$ value of the $S_{2-2}$ area is equal to $0.11\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{\mathrm{Hb}}⟩$ value of the $C_2$ area is equal to $-0.02\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{\mathrm{Hb}}⟩$ values of the motor and sensory areas are statistically different from the one of the $C_2$ area. For video 4, the $⟨\overline{\mathrm{Hb}}⟩$ value of the $M_{1-2}$ is equal to $0.06\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $\overline{⟨\mathrm{Hb}⟩}$ values of the $S_{1-2}$, $S_{2-2}$, $S_{3-2}$, and $C_2$ areas are smaller than or equal to $-0.18\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{\mathrm{Hb}}⟩$ values of the $M_{1-2}$ and $S_{1-2}$ areas are statistically different from the one of the $C_2$ area, whereas the ones of the $S_{2-2}$ and $S_{3-2}$ areas are not. For video 5, the $\overline{⟨\mathrm{Hb}⟩}$ values of the motor areas are smaller than or equal to $-0.28\mu \mathrm{mol}{\mathrm{L}}^{-1}$). The $\overline{⟨\mathrm{Hb}⟩}$ values of the $C_3$ area is equal to $0.11\mu \mathrm{mol}{\mathrm{L}}^{-1}$). The $⟨\overline{\mathrm{Hb}}⟩$ values of the motor areas are statistically different from the one of the $C_3$ area.

In Fig. 10, for video 1, the $⟨\overline{{\mathrm{HbO}}_2}⟩$ values of the motor and sensory areas are higher than or equal to $0.65\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{{\mathrm{HbO}}_2}⟩$ value of the $C_1$ area is equal to $-0.23\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{{\mathrm{HbO}}_2}⟩$ values of the motor and sensory areas are statistically different from the one of the $C_1$ area. For video 2, the $⟨\overline{{\mathrm{HbO}}_2}⟩$ values of the motor and sensory areas are higher than or equal to $1.57\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{{\mathrm{HbO}}_2}⟩$ value of the $C_2$ area is equal to $-0.16\mu \mathrm{mol}{\mathrm{L}}^{-1}$). The $⟨\overline{{\mathrm{HbO}}_2}⟩$ values of the motor and sensory areas are statistically different from the one of the $C_2$ area. For video 3, the $⟨\overline{{\mathrm{HbO}}_2}⟩$ values of the $M_{1-2}$, $S_{1-2}$, and $S_{3-2}$ areas are higher than or equal to $0.85\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{{\mathrm{HbO}}_2}⟩$ value of the $S_{2-2}$ area is equal to $-0.43\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{{\mathrm{HbO}}_2}⟩$ value of the $C_2$ area is equal to $0\mu \mathrm{mol}{\mathrm{L}}^{-1}$). The $⟨\overline{{\mathrm{HbO}}_2}⟩$ values of the motor and sensory areas are statistically different from the one of the $C_2$ area. For video 4, the $\overline{⟨{\mathrm{HbO}}_2⟩}$ values of the $S_{1-2}$ and $S_{2-2}$ areas are higher than or equal to $0.59\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{{\mathrm{HbO}}_2}⟩$ value of the $M_{1-2}$ area is equal to $0.08\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{{\mathrm{HbO}}_2}⟩$ values of the $S_{3-2}$ and $C_2$ areas are higher than or equal to $0.15\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{{\mathrm{HbO}}_2}⟩$ values of the $M_{1-2}$, $S_{1-2}$, and $S_{2-2}$ areas are statistically different from the one of the $C_2$ area, whereas the one of the $S_{3-2}$ is not. For video 5, the $⟨\overline{{\mathrm{HbO}}_2}⟩$ values of the motor areas are higher than or equal to $0.57\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{{\mathrm{HbO}}_2}⟩$ value of the $C_3$ area is equal to $-0.30\mu \mathrm{mol}{\mathrm{L}}^{-1}$. The $⟨\overline{{\mathrm{HbO}}_2}⟩$ values of the motor areas are statistically different from the one of the $C_3$ area.

In Fig. 11, for video 1, the $⟨r_{\mathrm{Hb}}⟩$ values of the motor and sensory areas are higher than or equal to 0.49. The $⟨r_{\mathrm{Hb}}⟩$ value of the $C_1$ area is equal to 0.15. The $⟨r_{\mathrm{Hb}}⟩$ values of the motor and sensory areas are statistically different from the one of the $C_1$ area. For video 2, the $⟨r_{\mathrm{Hb}}⟩$ values of the motor and sensory areas are higher than or equal to 0.41. The $⟨r_{\mathrm{Hb}}⟩$ value of the $C_2$ area is equal to 0.09. The $⟨r_{\mathrm{Hb}}⟩$ values of the motor and sensory areas are statistically different from the one of the $C_2$ area. For video 3, the $⟨r_{\mathrm{Hb}}⟩$ values of the $M_{1-2}$, $S_{1-2}$, and $S_{3-2}$ areas are higher than or equal to 0.48. The $⟨r_{\mathrm{Hb}}⟩$ values of the $S_{2-2}$ and $C_2$ areas are smaller than or equal to 0.09. The $⟨r_{\mathrm{Hb}}⟩$ values of the $M_{1-2}$, $S_{1-2}$, and $S_{3-2}$ areas are statistically different from the one of the $C_2$ area, whereas the one of the $S_{2-2}$ area is not. For video 4, the $⟨r_{\mathrm{Hb}}⟩$ value of the $S_{1-2}$ area is equal to 0.50. The $⟨r_{\mathrm{Hb}}⟩$ values of the $M_{1-2}$, $S_{2-2}$, $S_{3-2}$, and $C_2$ areas are smaller than or equal to 0.23. The $⟨r_{\mathrm{Hb}}⟩$ values of the $S_{1-2}$ and $S_{2-2}$ areas are statistically different from the one of the $C_2$ area, whereas the ones of the $M_{1-2}$ and $S_{3-2}$ areas are not. For video 5, the $⟨r_{\mathrm{Hb}}⟩$ values of the motor areas are higher than or equal to 0.30. The $⟨r_{\mathrm{Hb}}⟩$ value of the $C_3$ area is equal to 0.02. The $⟨r_{\mathrm{Hb}}⟩$ values of the motor areas are statistically different from the one of the $C_3$ area.

In Fig. 12, for video 1, the $⟨r_{{\mathrm{HbO}}_2}⟩$ values of the motor and sensory areas are higher than or equal to 0.38. The $⟨r_{{\mathrm{HbO}}_2}⟩$ value of the $C_1$ area is equal to 0.12. The $⟨r_{{\mathrm{HbO}}_2}⟩$ values of the motor and sensory areas are statistically different from the one of the $C_1$ area. For video 2, the $⟨r_{{\mathrm{HbO}}_2}⟩$ values of the motor and sensory areas are higher than or equal to 0.38. The $⟨r_{{\mathrm{HbO}}_2}⟩$ value of the $C_2$ area is equal to 0.11. The $⟨r_{{\mathrm{HbO}}_2}⟩$ values of the motor and sensory areas are statistically different from the one of the $C_2$ area. For video 3, the $⟨r_{{\mathrm{HbO}}_2}⟩$ values of the $M_{1-2}$, $S_{1-2}$, and $S_{3-2}$ areas are higher than or equal to 0.46. The $⟨r_{{\mathrm{HbO}}_2}⟩$ values of the $S_{2-2}$ and $C_2$ areas are smaller than or equal to 0.13. The $⟨r_{{\mathrm{HbO}}_2}⟩$ values of the motor and sensory areas are statistically different from the one of the $C_2$ area. For video 4, the $⟨r_{{\mathrm{HbO}}_2}⟩$ value of the $S_{1-2}$ area is equal to 0.41. The $⟨r_{{\mathrm{HbO}}_2}⟩$ values of the $M_{1-2}$, $S_{2-2}$, $S_{3-2}$, and $C_2$ areas are smaller than or equal to 0.23. The $⟨r_{{\mathrm{HbO}}_2}⟩$ values of the $S_{1-2}$ and $S_{2-2}$ areas are statistically different from the one of the $C_2$ area, whereas the ones of the $M_{1-2}$ and $S_{3-2}$ area are not. For video 5, the $⟨r_{{\mathrm{HbO}}_2}⟩$ values of the motor areas are higher than or equal to 0.25. The $⟨r_{{\mathrm{HbO}}_2}⟩$ value of the $C_3$ area is equal to 0.06. The $⟨r_{{\mathrm{HbO}}_2}⟩$ values of the motor areas are statistically different from the one of the area $C_3$.

According to the criteria defined in Sec. 2.4.4, the functional areas represented in Fig. 7 are designated as activated or nonactivated cortical areas; see Table 3.

Table 3

Definition of the functional areas identified by EBS (represented in Fig. 7) as activated or nonactivated.

Video ID	Functional areas identified by EBS	Statistical significance				Cortical activation measurement
		⟨Hb‾⟩	⟨HbO2‾⟩	⟨rHb⟩	⟨rHbO2⟩
Video 1	$M_{1-1}$	✓	✓	✓	✓	✓
	$M_{2-1}$	✓	✓	✓	✓	✓
	$M_{3-1}$	✓	✓	✓	✓	✓
	$S_{1-1}$	✓	✓	✓	✓	✓
	$S_{2-1}$	✓	✓	✓	✓	✓
Video 2	$M_{1-2}$	✓	✓	✓	✓	✓
	$S_{1-2}$	✓	✓	✓	✓	✓
	$S_{2-2}$	✓	✓	✓	✓	✓
	$S_{3-2}$	✓	✓	✓	✓	✓
Video 3	$M_{1-2}$	✓	✓	✓	✓	✓
	$S_{1-2}$	✓	✓	✓	✓	✓
	$S_{2-2}$	✓	✓	✗	✓	✗
	$S_{3-2}$	✓	✓	✓	✓	✓
Video 4	$M_{1-2}$	✓	✓	✗	✗	✗
	$S_{1-2}$	✓	✓	✓	✓	✓
	$S_{2-2}$	✗	✓	✓	✓	✗
	$S_{3-2}$	✗	✗	✗	✗	✗
Video 5	$M_{1-3}$	✓	✓	✓	✓	✓
	$M_{2-3}$	✓	✓	✓	✓	✓

4.

Discussion

In videos 1, 2, and 5 (see Table 1), the statistical results of Table 3 indicate that the motor and sensory areas identified by EBS correspond to the highlighted cortical areas of Fig. 7. This demonstrates the ability of our functional model to identify in a robust manner these functional areas. In videos 3 and 4, the motor and sensory cortex was explored in a more subtle way. In video 3, the hand movement was induced by an external person. All motor and sensory areas are defined as activated cortical areas except the $S_{2-2}$ area. In video 4, the stimulation of the cortex was achieved through a repetitive fingers and palm caresses. Only the $S_{1-2}$ area is defined as an activated cortical area. These results seem to indicate that the $S_{1-2}$ area could be linked to somatosensory functions, whereas the $S_{2-2}$ and $S_{3-2}$ areas to sensorimotor functions. This indicates that our model seems to be rather sensitive to subtle differences in the physiological stimuli of the hand.

These results should be taken with caution because our functional model needs to be improved. In particular, a physiological a priori on the hemoglobin concentration change values could be incorporated in the model. However, there is still no consensus in the literature on such threshold values linked with a physiological stimulus. In video 3, when the $\overline{\mathrm{Hb}}$ and $\overline{{\mathrm{HbO}}_2}$ distributions of the $S_{2-2}$ area are statistically significant (see Table 3), the $⟨\overline{\mathrm{Hb}}⟩$ and $⟨\overline{{\mathrm{HbO}}_2}⟩$ values of this area differ from the area $C_2$ in the sense of a deactivation of the cortical area, which is difficult to interpret (the $⟨\overline{\mathrm{Hb}}⟩$ value of the $S_{2-2}$ area is higher than the one of the $C_2$ area and the $⟨\overline{{\mathrm{HbO}}_2}⟩$ value of the $S_{2-2}$ area is smaller than the one of the $C_2$ area). The $S_{2-2}$ area is defined as nonactive because its $⟨r_{\mathrm{Hb}}⟩$ value is not statistically different from the one of the $C_2$ area. This result shows that our functional model is able to eliminate this nonphysiologic hemodynamic event by comparing the similarity of the hemoglobin changes time courses to the expected hemodynamic response. This confirms that our model is more robust than models with a nonphysiological a priori. The same phenomenon can be observed for the $M_{1-2}$ area of video 4.

Our functional model has to be used with the EBS. Indeed, it requires a reference cortical area which is defined as a nonfunctional area by the EBS. The main limits are that the results of the statistical analysis depend of the choice of the reference and the definition of the cortical activation is not obtained for each camera pixel but for a group of pixels. A more robust definition of cortical areas will be explored in future works, notably with the implementation of a SPM⁴⁰ like analysis.

The $\overline{\mathrm{Hb}}$, $\overline{{\mathrm{HbO}}_2}$, $r_{\mathrm{Hb}}$, and $r_{{\mathrm{HbO}}_2}$ noise levels in our measurements were calculated by taking the standard deviation of the mean values of $10\times 5$ cortical areas which have not been identified as functional areas by the EBS (10 measurements were realized in each of the five videos). The $\overline{\mathrm{Hb}}$ and $\overline{{\mathrm{HbO}}_2}$ noise levels are equal to $0.26\mu \mathrm{mol}{\mathrm{L}}^{-1}$, and the $r_{\mathrm{Hb}}$ and $r_{{\mathrm{HbO}}_2}$ noise levels are equal to 0.07. In video 4, the $S_{2-2}$ area is not defined as an activated area only because the $⟨\overline{\mathrm{Hb}}⟩$ value is not statistically different from the one of the $C_2$ area. The mean values of the other distributions are statistically different from the ones of the $C_2$ areas and the sign of the differences corresponds to the sense of an activation of the cortical area. When comparing the $⟨\overline{\mathrm{Hb}}⟩$, $⟨\overline{{\mathrm{HbO}}_2}⟩$, $⟨r_{\mathrm{Hb}}⟩$, and $⟨r_{{\mathrm{HbO}}_2}⟩$ values of the $S_{2-2}$ and $C_2$ areas, the differences between these values are approximately equal to the noise levels. Thus by increasing the signal-to-noise ratio, the activated cortical areas could be determined in a more robust way.

The five videos have been studied separately, so no general PCCT value was proposed to localize the activated cortical areas in the five videos. The $\overline{\mathrm{Hb}}$, $\overline{{\mathrm{HbO}}_2}$, $r_{\mathrm{Hb}}$, and $r_{{\mathrm{HbO}}_2}$ distributions of the functional areas identified by EBS have different mean and standard deviation values. For instance, the $r_{{\mathrm{HbO}}_2}$ values of the $M_{1-2}$ area in video 2 range between $0.42\pm 0.25$ and the ones of the same motor area range between $0.53\pm 0.22$ in video 3. The $r_{{\mathrm{HbO}}_2}$ values of video 3 are higher than the ones of video 2. This may be due to a more constant hand movement in video 3 than in video 2. Indeed, in video 2, the patient has moved himself his hand, whereas the movement was induced by an external person in video 3. In our study, the stimulation step of the experimental paradigm (see Sec. 2.2) is not exactly constant during 30 s (patient sleepness and irregular movement). This implies that the experimental paradigm function plotted in Fig. 3 does not exactly fit a window function. Since an ideal expected hemodynamic response (convolution of the HIRF²⁶ with a window function) is compared to the $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ time courses, the $r_{\mathrm{Hb}}$ and $r_{{\mathrm{HbO}}_2}$ distributions may differ depending on the quality of the stimulus. Furthermore, habituation is known to interfere with the linearity of the hemodynamic response during a 30-s plateau. This would further modify the actual stimulus from the ideal window function.

The model of the HIRF may be improved. The positive BOLD signal used in MRI studies occurs in superficial cortical layers (0 to 1 mm). This positive BOLD signal is often accompanied with a negative BOLD signal occurring in deeper cortical layers (1 to 2 mm).⁴¹ The HIRF may be a conjunction of these two BOLD signals. In our study, we consider that the theoretical Hb impulse response function correspond to the opposite of the ${\mathrm{HbO}}_2$ impulse response function. This assumption may introduce some errors since it has been shown that these HIRFs do not correspond to an opposite operation.⁴² It also has been shown that the HIRF is different in the gray matter and in the arterioles.⁴² To improve our functional model, several HIRFs may be determined for each camera pixel depending on its biological attribution (gray matter, surface blood vessel, and buried blood vessel). For this purpose, the segmentation step (see Sec. 2.3.1) could be used.

The highest $\overline{\mathrm{Hb}}$ and $\overline{{\mathrm{HbO}}_2}$ values are localized at the level of the blood vessels (see Table 2). Indeed, the partial volume effect has less impact on the blood vessels than on the gray matter. So the underestimation of the concentration changes is lower at the level of the blood vessel than at the level of the gray matter. It is also interesting to notice that the concentration changes of the blood vessels which are not directly linked to the functional areas are higher than those of gray matter pixels associated with the sensory and motor areas; see Fig. 8. Similar observations can be found in the literature.^18,25 Pichette et al.¹⁸ proposed to mask the blood vessels in order to render more clearly the smaller hemodynamic changes in gray matter. The computation of the quantitative maps in conjunction with the measurement of the Pearson correlation coefficient between the concentration changes time courses with the expected hemodynamic response allows us to visualize the activated cortical areas without masking the blood vessels. With the results of this study, it is complicated to understand the role of the blood vessels. Indeed, it is impossible to get the direction of the blood flow and thus conclude that a blood vessel is associated with a functional area. A complementary imaging modality such as speckle imaging⁴³ could help to interpret the role of these blood vessels.

In Sec. 2.4.2, three simple models for the light transport of photons in cortical tissue are defined for the analysis of reflectance spectra. These models are applied to each camera pixel depending on its biological attribution (gray matter, surface blood vessel, and buried blood vessel). This segmentation allows us to approach the optical properties inhomogeneities of a cortical tissue. If the mean path length of the gray matter model (see the curve of model 3 in Fig. 5) was used for the computation of the $\mathrm{\Delta }[\mathrm{Hb}]$ and $\mathrm{\Delta }[{\mathrm{HbO}}_2]$ time courses of the surface blood vessel point $A$, this would imply a 40% underestimation of the $\overline{\mathrm{Hb}}$ value and 35% underestimation of the $\overline{{\mathrm{HbO}}_2}$ value. In our model, the distance of a gray matter pixel to the nearest blood vessel is not taken into consideration. We consider that the mean path length associated with this pixel is model 3 curve in Fig. 5. This is an approximation since the blood vessel still has an impact on the photon propagation. To solve this problem, a pixel-wise determination of the optical mean path length can be implemented. This could be realized by fitting theoretical reflectance spectra (computed for a great number of Monte-Carlo models) with measured reflectance spectra. This comparison appears not to be an easy task since only three spectral values can be obtained which are integrated over the broad wavelength ranges covered by each detector. A hyperspectral camera seems to be a more appropriate solution to answer this problem with the acquisition of high-resolution reflectance spectra. A similar method has been implemented by Pichette et al.¹⁸ using a hyperspectral imaging device and a two class image segmentation but does not take account of the distance between cortical tissues and blood vessel.

The quantitative maps are displayed after the end the data acquisition. Unlike Bouchard et al.’s²⁵ work, our method is not processed in real time. But it would be possible to develop a real-time processing using a multithreading software architecture and a code optimization using GPU (graphics processing units).

5.

Conclusion

A quantitative and noninvasive method for imaging the oxygenated hemoglobin concentration changes ($\mathrm{\Delta }[{\mathrm{HbO}}_2]$) and deoxygenated hemoglobin concentration changes ($\mathrm{\Delta }[\mathrm{Hb}]$) of a human exposed cortex using a digital RGB camera was demonstrated in the present report. The results showed that the method can be used by the neurosurgeon to localize the motor and sensory areas before a tumor removal surgery. The brain areas defined as activated cortical areas by our functional model are well correlated with the functional areas localized by EBS. Our model is based on quantitative measurements and physiological and statistical comparisons. Each step of the model can be improved, such as the noise level reduction, the consideration of a physiological a priori in the statistical comparisons, or the precision of the quantitative measurements. Indeed, the estimation of three different path lengths associated with three different cortical tissues was used to process the Beer–Lambert law. The Monte-Carlo models used in our study remain an approximation of the optical properties inhomogeneities of the patient brain. To increase the accuracy of our method, a detailed tissue structure that best matches the cortical pattern of the acquired images can be implemented in the Monte-Carlo simulations. However, this would imply the determination of a three-dimensional cortical tissue from an RGB image and would require considerably more computational time and a large amount of simulated data. These directions will be explored in future works. The current work demonstrates that an RGB camera combined with a quantitative modeling of brain hemodynamics biomarkers could evaluate in a robust way the functional areas during neurosurgery. This strengthens the relevance of using a classical RGB camera for functional intraoperative brain imaging.