NeuroRegulation http://www.isnr.org

The combination of modern machine learning and traditional statistical methods allows the construction of individual regression models for predicting the blood oxygenation level dependent (BOLD) signal of a selected region-of-interest within the brain using EEG signal. Among the many different models for motor cortex, we chose the EEG Fingerprint one-electrode approach, based on rigid regression model with Stockwell EEG signal transformation, used before only for the amygdala. In this study we demonstrate the way of finding suitable model parameters for the cases of BOLD signal reconstruction for five individuals: three of them were healthy, and two were after a hemorrhagic stroke with varying degrees of damage according to Medical Research Council (MRC) Weakness Scale. The principal possibility of BOLD restoring using regressor model was demonstrated for all the cases considered above. The results of direct and indirect comparisons of BOLD signal reconstruction at the motor region for healthy participants and for patients who suffered from a stroke are presented.


Introduction
There are certain advantages of combining functional magnetic resonance imaging (fMRI) and electroencephalographic (EEG) modalities. The EEG data has an excellent submillisecond resolution which allows for an adequate representation of rich temporal dynamics of the neural ensembles' activity recorded from the scalp surface. The issues arising from the limited number of electrodes used to record 2D signals and facial skin and muscles' conductivity result in a limited spatial resolution of the signal and inadequate signal-to-noise ratio. This is also known as an ill-posed inverse problem of reconstruction of the source of EEG activity (Custo et al., 2014). FMRI allows to improve the spatial resolution to a submillimeter level, while temporal dynamics of the blood oxygenation level dependent (BOLD) signal is generally recorded with few seconds repetition time (Ogawa et al., 1990) . It should be noted that the BOLD signal does not directly represent a neural activity, but rather is a currently poorly understood combination of hemodynamic response and metabolic processes characterized by a relatively slow temporal dynamic (Logothetis, 2002).
In addition, the analysis of BOLD signals reveals spatial patterns that are correlated by time and may represent functional connectivity Smith et al., 2009) . Lastly, it is a known fact that the parameters of the hemodynamic response function (HRF) differ depending on a certain brain area and varies from subject to subject (Handwerker et al., 2012). Thus, perspectives of studying dynamics of brain metabolic response on the smaller time scale will result in a better understanding of various basic brain processes.
Traditional methods of EEG signal reconstruction using EEG band power-based regressors convoluted with the HRF do not fully reflect the whole picture of nuances and complexities of the EEG phenomena (Laufs et al., 2003;Murta et al., 2017;Yin et al., 2016). Due to limitations of the time scale and the number of available electrodes, as well as predetermined frequency ranges (e.g., alpha, beta, gamma, etc.), the resulting structure of the EEG becomes oversimplified (Marecek et al., 2016). The authors consider a model presented in the study by Meir-Hasson et al. (2016) to be one of the most promising models. A group of subjects underwent scanning during the so-called "emotional regulation," and the resulting data showed the efficacy of the authors' method of reconstruction of BOLD signal from the region of interest (ROI) based on EEG data recorded during the biofeedback (Keynan et al., , 2019. In the described approach, there were no limitations to HRF shape, and respective time delays were factored into the model. The authors used the Stockwell transformation (ST) of EEG signal followed by a down sampling of frequency in the time-frequency domain (Stockwell, 2007). The increase of the temporal resolution of reconstructed BOLD signal was discussed as a part of the model which created the possibility of better understanding its temporal dynamics. This method was claimed as a universal one that could be applied for various functional brain networks and, more importantly, allowed for the simulation of signal of intracerebral ROIs and not just cortical structures.
However, to the best of our knowledge, the model was tested only for the amygdala; therefore, its versatility and applicability to other brain structures and predictive power in that case require further investigation. Several methodical issues such as specific frequency ranges, bands, electrodes, or other model hyperparameters are subject to further research and adjustment. The aim of the current study based on the aforementioned Meir-Hasson model (2016) was to investigate possibility and accuracy of the sensorimotor (SMR) BOLD signal reconstruction that reflects hand motion from EEG recordings. Relative homogeneity and simplicity of motor action should lead to effective prediction with adequate determination coefficients. We expect that in a case of motor area activity the predicted impact of specific EEG bands such as mu-and SMR-bands will be the greatest one.

Equipment, Materials, and Subjects
The study was performed at the tomography center using an Ingenia magnetic resonance machine (Philips, Amsterdam, Netherlands) with a magnetic field induction of 3.0T equipped with iViewBold software for real-time fMRI results. A head coil with a slanting mirror based on dStream technology and a mobile pixels-compatible monitor (NordicNeuroLab, Bergen, Norway) were used. Reference anatomical images were obtained using the T1 Mapping 3D turbo field echo (TFE) method, TR (repetition time) / TE (echo time) = 7.5 / 3.7 ms with a field of view of 250 x 250 x 180 mm 3 and a reconstructed voxel size (3D image element) of 1 x 1 x 1 mm 3 . The reconstruction included the creation of three packs of slices oriented along the main orthogonal planes (Multiplanar Reformation Procedure [MPR], i.e., multiplanar reconstruction). The main working T2 weighted images (T2WI) were obtained by the secure shell (SSH) protocol, Standard Portable Intermediate Representation (SPIR) echo planar imaging (EPI) method, TR / TE = 2500 / 35 ms with a sagittal orientation of a packet of 25 slices 5 mm thick and a field of view of 220 x 220 mm 2 , with a resolution on a 2 x 2 mm 2 plane. The uniform distribution of the slices in the TR interval provided a constant recording rate of 100 ms per slice, and their nonstandard sagittal arrangement simplified the identification of target zones. The relatively large thickness of the slices made it possible to limit their number and achieve optimal sensitivity at a time of TE = 35 ms. BrainVision (Brain Products GmbH, Gilching, Germany) encephalograph and a BrainCap MR-compatible helmet for EEG BrainCap MR (EASYCAP GmbH, Wörthsee, Germany), 32-channel/128-channel Ag/AgCl electrodes including the reference one, using an extended 10-20 system, were synchronously monitored for electrical activity of the brain. An electrode for recording ECG was placed under the shoulder blade of the subject. Before placing the participants in the tomograph, the level of electrode impedance < 20 kOhm was achieved.
In total, five right-handed subjects took part in the experiment, see Table 1. Three of five subjects were healthy individuals (S1, S2, S3 in Table 1 experiment had a block design: each subject used his hand in response to the stimulus on the monitor screen (the task was to squeeze a ball for the entire duration of stimulus). During the session, time of stimulus presentation was equivalent to 25 s, followed by a rest period of the same duration, and the total number of scan recordings was 240 s; i.e., 600 s per session, in total 24 blocks per session. This paradigm was repeated for the left hand only. The data obtained in each experiment (first phase) was subsequently used to build the EEG-fMRI regression model. The second phase of the experiment was conducted outside of the MRI scanner. Same subjects participated, and the second session of the first phase was executed with the same stimuli along with the recording of EEG data using Neuron-Spectr-5 commercial electroencephalograph (Neurosoft, Ivanovo, Russia). In the experiment the MCScap textile cap (Medical Computer Systems, MCS, Moscow, Russia) with 30 removable passive Ag/AgCl electrodes was placed according to the 10-20 system. Montage was referenced to Cz electrode.

Signal Refinement EEG
Recording EEG signal in the MR scanner is affected by gradient exposure due to the fast magnetic field alternations. The amplitudes of physiological component of the signal decline on a logarithmic scale, while prominence of the gradient (MR) component does not fade as fast. Thus, during the signal processing phase it is necessary to closely monitor the process of EEG signal refinement from MR gradient, so that the physiological component of the spectrum does not degrade. We used an EEGLab Software Package for that purpose (Delorme & Makeig, 2004). During the refinement process, Bergen plugin (Moosmann et al., 2009) was used for subtraction of the gradient artifacts, as well as Functional Magnetic Resonance Imaging of the Brain (FMRIB) plugin was applied for pulse artifacts removal Niazy et al., 2005) using standard settings. For the signal free of magnetic artifact, our cleaning procedure was close to the Makoto's preprocessing pipeline steps (Swartz Center for Computational Neuroscience, n.d.). To remove occasional large amplitude noise/artifact subspace reconstruction (ASR), clean raw data plugin was chosen (Chang et al., 2018). For the constant fixed source noise/artifact signals independent component analyses (ICA) was used, where an automatic classification label for independent components (IClabel) EEGLab plugin was run to classify which of the independent components had to be subtracted (Pion-Tonachini et al., 2019). All algorithms used have their inner parameters as the number of principal component analysis (PCA) components (FMRIB plugin) or level of probability (IClabel plugin); therefore, the cleaning procedure was done always under visual inspection of signal and its spectral characteristic.

Analysis of fMRI Data
Data in international DICOM imaging standard format were converted to Neuroimaging Informatics Technology Initiative (NIFTI) format using the MRIConvert utility (University of Oregon, Lewis Center for Neuroimaging, n.d.). To obtain BOLD signal the statistical parametric mapping (SPM) fMRI data analysis were used. Statistical analysis was performed in the SPM package (UCL Queen Square Institute of Neurology, n.d.). Standard settings were applied unless otherwise stated. The following processing steps were also performed: correction of mutual frame positions (Reslice) to eliminate the consequences of subjects' head movements and time differences in obtaining individual slices (Slice Timing), coregistration of structural and mean functional images (Coregister), division into gray and white matter and areas of cerebrospinal fluid with simultaneous transformation to the standard MNI (Montreal Neurological Institute) coordinate space with an isotropic distribution of 2 mm (Segment, Normalize: Write) and smoothing (Smooth) of functional tomograms using a spatial filter with a Gaussian nucleus of half-height width 6 mm. At the individual level, tasks were modeled by conversion of the standard function of the hemodynamic response and the rectangular function corresponding to the work and rest blocks. The design also included additional regressors: head movement correction parameters (three axes of displacement and three planes of rotation) to eliminate corresponding noise. This stage also employed a high-frequency filter with a cutoff period Vol. 8 (3)

Model Description
Availability of the signals with high spatial resolution (fMRI) and high temporal resolution (EEG) made it possible to build the model for BOLD signal reconstruction. It was shown by Meir-Hasson et al. (2014) that combining these two signals recorded from the amygdala, defining a connection mathematically, and then computing a regression model that predicts BOLD signal based on EEG data was possible. The certain steps to build such a model are shown in Figure 1. In the case of BOLD signal, the transformation procedure is rather trivial; the BOLD signal is normalized with a subsequent increase of the sampling rate to 4 Hz.
Transformation of EEG signal should be performed keeping in mind temporal dynamics of the BOLD signal. Hemodynamic response value reaches its peak to 5-to 6-seconds point from the stimulus onset, the whole response is somewhat about 10 seconds, and, strictly speaking, this function is heterogeneous with respect to a certain brain area. Thus, for the further EEG signal processing a 12-s time frame was selected for one of the electrodes for EEG signal with 80 Hz sampling frequency. For a better identification of the relevant EEG signal characteristics for a subsequent regression analysis, the 12-s frame was transformed into a timefrequency region using Stockwell transformation (ST) procedures implemented in the MNE-Python package (Gramfort et al., 2013).
To reduce number of the signal characteristics, the frequency band was limited to 30 Hz, and then divided into three frequency unknown subbands (to be found). In the time region the sampling rate was reduced from 80 Hz to 4 Hz, while the BOLD signal had the same frequency. After all transformations (five steps), one time frame of EEG signal resulted into a total of 144 parameters (12 s x 4 Hz x 3 bands) that could be used as a vector of independent values for the ridge regression ( Figure  1); i.e., to one value of the BOLD signal. Moving the whole 12-s time frame by 0.25 s, other 144 parameters (regressors) of transformed EEG signal can be used with the next value of the BOLD signal. The BOLD signal in this case is a dependent variable.

Optimizing Model Hyperparameters
At the beginning of the construction of the model, it was decided to limit the frequency band of the EEG signal to 6-30 Hz and consider only three subband frequency intervals. Thus, the model is characterized by nine unknown hyperparameters including: width -ST transformation coefficient; alpharidge regression coefficient; electrode; and six frequency parameters that define three subfrequency intervals.

Building Individual Models
The subjects performed the motor task "squeeze the ball" with a simultaneous recording of the fMRI and EEG signals. During the first phase, after removal of tomographic and miscellaneous artifacts from the data, a regression model for each subject was created for two halves of experimental data. Prior to BOLD signal reconstruction with a commercial electroencephalograph, we estimated the accuracy of the model using the r 2 metric for all subjects.

Electrode Selection, Alpha and ST Width
The procedure for determining the required values of hyperparameters was as follows. At the first stage, a restriction was introduced on the number of subfrequency intervals, namely, one frequency interval of 9-25 Hz was selected related to SMR activity, including classical mu-rhythm and upper beta rhythm (Emmert et al., 2016). The values of the ST width parameter and the alpha parameter of the ridge regressor were sorted out for each electrode.
The model building process went through the crossvalidation procedure with k-fold = 10. Data set was split into k = 10 consecutive folds (without shuffling). Then each fold was used as a validation set once while the k − 1 remaining fold formed the training set. The overall model was optimized as the mean of five subjects. The following task was set: for each value of alpha, ST width and electrode to find such values of the coefficient regression models that gives the maximum values of the correlations of the summary model: ( 1, 2) where ( 1, 2) is the individual ridge regressor model.
As expected, the best electrode for constructing a model for left-hand movement was C4 contralateral electrode located closest to the motor region of the right brain hemisphere; 0.1 width ST; and alpha 20.000 (see Figure 2).

Frequency Sub-Bands Selection
At the second stage for the further optimization of previously found hyperparameters C4 electrode, 0.1 width ST, and alpha 20000 with three subfrequency intervals were used, giving six unknown variables.
To search for a set of hyperparameters that would produce the best results the Hyperopt optimization library was chosen (Bergstra et al., 2013), that is faster than a simple grid search.
The optimization procedure was applied simultaneously to five subjects. The number of iterations was 5000. The BOLD signal obtained directly from fMRI data by extraction of the voxel intensities described earlier is highlighted orange. The BOLD signal obtained as a result of reconstruction is shown in blue; the resulting model from the first half of the experimental data was applied to the EEG signal of the second half of the experimental data. This procedure for constructing a model with its subsequent verification with a real BOLD signal was applied for the model obtained from the data of the second half of the experiment.
Frequency bands for the best result were found to be, consequently, 8-12, 12-14, and 16-22 Hz. Thus, a complete set of nine hyperparameters was found, which made it possible to start building individual models. The result of applying individual model to the selected frequency band (8)(9)(10)(11)(12)(12)(13)(14)(16)(17)(18)(19)(20)(21)(22) is shown at Figure 3. The correlation values are shown in Table 2.  Note. On the left 1.a) the model was built according to the first half of experimental data, on the right 1.b) according to the second half of data. 2.a), the individuals' models with three selected frequency bands were built according to the first half of experimental data, 2.b) according to the second part of data.

Reconstruction of BOLD-Dependent EEG Signal Building Pseudo-BOLD Signal
The same subjects further underwent a second session outside the scanner with not MRcompatible electroencephalograph Neuron-Spectr-5 used for the signal recording and for implementing a previously created model.
Due to the absence of the BOLD signal, it was not possible to directly compare the predicted BOLD signal with the real one. An averaged BOLD signal curve based on the analysis of actual BOLD signals was introduced, the middle point of which matched the label of the end of hand motor task (Figure 4). Correlation coefficient between this pseudo-BOLD signal and the reconstructed BOLD signal was accepted as an estimated metric of likelihood. Note. By labels, i.e., the moment when block associated with movement starts (dash lines). Left: the real BOLD signal; right: the signal reconstruction based on labels.

Reconstruction
Previously constructed individual regression models were applied to the EEG signals recorded with a commercial electroencephalograph. The result is presented in Figure 5, where reconstructions of BOLD signals built according to individual models of subjects according to data obtained with a commercial electroencephalograph are shown.

Figure 5. Reconstruction of the BOLD Signal by the EEG pattern.
Note. Only reconstruction marked by asterisks from Table 3  Orange color shows the model BOLD signal reconstructed by event labels. They cannot be compared directly; however, the shape time points of the signal values that increase and decrease in five experiments appear to be very similar to the normal dynamics of the BOLD signal registered in the MR scanner.
More information on reconstruction is presented in Table 3. Left column indicates patients depending on I or II part of data. Columns M(I2), M(I1), M(S1), M(S2), M(S3) show individual regression models for patients S1, S2, S3, I1, I2, where letter M is reserved for model. Channel columns (CH) represent EEG recording channels used for BOLD reconstruction for patients S1, S2, S3, I1, I2. To avoid such problems in the future, we plan to record the coordinates of the electrodes.

Discussion
The new tendencies in the neurological research originate from the desire to preserve the benefits of the fMRI-EEG data without using MR scanner. This research task is labeled as a problem of the BOLDdependent EEG. Some research studies involving simultaneous records of fMRI and EEG signals utilize EEG spectral features that may be useful in predicting specific BOLD dynamics (de Munck et al., 2009;Goense & Logothetis, 2008;Goldman et al., 2002;Kilner et al., 2005;Mantini et al., 2007;Murta et al., 2015;Rosa et al., 2010;Wan et al., 2006). Such innovation necessity is motivated by the immobility and the cost of bimodal platform. Current studies attempt to retain its advantages outside the scanner and produce a "shadow" copy of fMRI using the BOLD-dependent EEG. The key idea is to use the obtained data on fMRI-EEG interdependencies with a modern commercial electroencephalograph.
Our study of SMR cortex BOLD signal prediction using EEG signal was based on the EEG-channelspecific BOLD prediction technique described in the articles by Meir-Hasson et al. (2014. Unlike the original study, the motor region of the cortex responsible for the real movements was chosen. To the best of our knowledge, there were not published any similar studies targeted to any cortical region of the brain except of amygdala. Despite the fact that this approach suits to different brain structures, it is necessary to adjust unique model hyperparameters; namely, certain electrodes, division of the EEG spectrum into frequency bands, providing the best reconstruction of the BOLD signal, ST width parameter, and, finally, rigid regressor alpha parameter. Indeed, the results of our study confirm the effectiveness of this approach and demonstrate that BOLD signal reconstructions by EEG have statistically significant levels of correlation coefficients with "ground truth" real and pseudo-BOLD signals. This means that the technique allows to avoid or significantly limit the usage of MR scanner without dramatic loss of quality of the reconstructed BOLD signal and to maintain its spatial precision in certain tasks.
In this study the individual EEG-fMRI models were built for each participant. A total of five subjects took part in the experiment. Three of them were healthy, one had a mild motor impairment, and one was with a severe poststroke motor impairment.
This article describes the procedure for obtaining hyperparameters needed for building a regression model for reconstruction of BOLD-dependent signal by EEG pattern that is valid for all five subjects. A direct comparison of the BOLD reconstruction with a real statistically significant correlation coefficient was carried out (9 of 10 comparisons have correlation coefficients more than 0.59, including a patient with severe paresis). An indirect comparison of data obtained on a commercial electroencephalograph showed correlation coefficients greater than 0.5 for all subjects. For patient I2, according to indirect comparison, the best electrode was found to be not C4, but CP2. A possible reason is different sizes of the cap used at the tomographic center where data in bimodal EEG + fMRI mode were obtained, and the cap used at the clinical electroencephalograph. To avoid this problem, it is planned to record coordinates of the electrodes in the future. The next step is moving away from individual models to a small number of universal ones (or even to one model fitting all) and thus limiting the use of the fMRI. According to data from Table 3, we may assume that this is possible. At this stage, the number of subjects was not sufficient to provide a statistically significant answer. Differences in individual physiological characteristics of participants make the construction of such models a complicated task. There are already proposed techniques such as hierarchical cluster analysis to overcome these difficulties (Meir-Hasson et al., 2014Wei et al., 2018).
As there are no methods of direct comparison of the reconstructed BOLD signal with a real one when only one electroencephalograph is used, it is planned to include objective monitoring tools (EMG sensors) in our future experiments. However, some open questions remain: the applicability of the model for different temporal designs, the relationship of hand's pressing force and ball manipulation, and model's adequacy in case of motor imaging.
The long-term goal is a usage of the BOLDdependent EEG in the rehabilitation of people who suffer from a stroke that resulted in damage of the motor area. There is a need to improve the existing ineffective paradigm of neurorehabilitation which results in up to 44% of disability after stroke (Katan & Luft, 2018). We hope that the technologies for reconstructing the BOLD using one or only a few EEG electrodes will potentially make it possible to Vol. 8 (3)

Conclusion
Demonstration of the possibility of restoring the BOLD-dependent signal, not only in healthy subjects but also in the patients with stroke, shows that this approach can be applicable to people with damage to the motor regions of the brain. The question of how static the obtained models are and whether they can be applied at the entire stage of rehabilitation of poststroke patients remains open, as well as the question of building universal models.
Our findings demonstrated that the model proposed by Meir-Hasson et al. (2016) and the hypothesis of its universality with respect to different brain structures are generally correct and may be used at least for the SMR cortex, although some additional tuning of frequencies and subband selection procedure, as well as different electrode should be used in this case.

Author Acknowledgement
Sponsored by Russian Science Foundation grant N 21-15-00209.

Compliance with Ethical Standards
Conflict of Interest. None of the authors have potential conflicts of interest to be disclosed.
Ethical Approval. All procedures in this research study were performed in accordance with the ethical standards compliant with the 1964 Helsinki declaration. Signed informed consent approved by the Federal Research Center of Fundamental and Translational Medicine Ethics Committee was obtained from all individuals who participated in the study.