## 2 Materials and Methods
### Data Processing  
All our custom signal processing methods were implemented in MATLAB (Mathworks, Natick, MA, USA, Version R2020b). At the pre-processing stage the data from all spectrometers were resampled in both time and spectral domains to the time step of 0.5 s and wavelength step of 1 nm between 650 and 1024 nm, filtered in the time domain using a band-pass filter with the window of 0.01-0.3 Hz, and smoothed in the spectral domain with the median filter of the 20 nm width. To analyze the chromophores, we used two custom data processing approaches. Our first algorithm was based on the analytical solution to the diffusion equation [30,31] for the semi-infinite homogeneous medium. This allows for measuring the bulk absolute concentrations of tissue HbO\({}_{2}\) and HHb, and changes in rCCO concentration without accounting for the layered tissue structure. The baseline concentrations of HbO\({}_{2}\) and HHb were calculated by performing a non-linear least square fitting [30,31] of the measured absorbance spectrum at each moment of time by the analytical solution to the diffusion equation, in which the optical absorption coefficient \(\mu_{a}(\lambda)\) was modeled as  
\[\mu_{a}(\lambda)=[\text{Hb}]\varepsilon(\lambda)_{\text{Hb}}+[\text{HbO}_{2}] \varepsilon(\lambda)_{\text{HbO}_{2}}+\eta(\lambda) \tag{1}\]  
According to Beer-Lambert law, the reduced scattering coefficient \(\mu_{s}{{}^{\prime}}(\lambda)\) as a function of wavelength \(\lambda\) was modeled using the power law as described in [30]. In Equation (1) and below the square brackets denote concentrations measured in micromoles (\(\upmu\)M), \(\eta(\lambda)\) is the absorption by water assuming 80% water concentration. The temporal changes in the hemoglobin concentrations HbO\({}_{2}\), Hb, and rCOO were resolved using a two-step data-fitting algorithm [30,31] (also based on the same analytical solution to the diffusion equation) by relating changes in HbO\({}_{2}\), Hb, and rCOO to the changes in the optical absorbance as  
\[\Delta\mu_{a}(\lambda,t)=\Delta[\text{Hb}](t)\varepsilon(\lambda)_{\text{Hb} }+\Delta[\text{HbO}_{2}](t)\varepsilon(\lambda)_{\text{HbO}_{2}}+\Delta[\text {CCO}](t)\varepsilon(\lambda)_{\text{ccO}}, \tag{2}\]  
where \(\varepsilon(\lambda)_{x}\) were the spectra of the extinction coefficients of HbO\({}_{2}\), Hb, and rCCO (see Figure 2) [30,31]. The data-fitting was performed in two steps: first \(\Delta\)[HbO\({}_{2}\)] and \(\Delta\)[HHb] were calculated assuming \(\Delta\)[CCO] = 0, and after that \(\Delta\)[CCO] was calculated and retained non-zero only if the addition of \(\Delta\)[CCO] resulted in an improvement in the fit quality of at least 20% in terms of the norm of residuals. Such a stepwise regression worked as an adaptive filter reducing the noise in the time-course of rCCO without suppressing changes recognized from significant spectral distortions corresponding to \(\varepsilon(\lambda)_{\text{cco}}\). Additional details on recovering the absolute values and changes of chromophore concentrations from the hNIRS data can be found in [32; 33]. As the fraction of oxygenated hemoglobin relative to the total hemoglobin in the blood, the cerebral tissue saturation of oxygen was calculated as:  
\[t\text{SO}_{2}=\frac{[\text{HbO}_{2}]}{[\text{HbO}_{2}]+[\text{Hb}]}(\%). \tag{3}\]  
In order to assess differences between the cerebral and scalp responses to BH we used another common NIRS model--the modified Lambert-Beer law for a two-layer medium [32; 33].