Optimization of phase prediction for brain-state dependent stimulation: a grid-search approach

Most experiments with non-invasive brain stimulation (NIBS) techniques, such as transcranial magnetic stimulation (TMS), aim at demonstrating the methodology's validity, understanding different stimulation parameters by comparing responses in different populations, and eventually suggesting new clinical applications. Yet, these investigations are treating the brain as a black box and do not consider the role of underlying brain activation, or state, variability. The resulting literature shows relevant heterogeneity, both intra- and inter-subject, which has triggered a shift in concept towards protocols, where the brain can no longer be considered as a black-box [1–4], but its ongoing activity, state, and oscillatory properties must be taken into account as they might strongly impact the response to TMS at the cortico-spinal (e.g. motor-evoked potentials, MEPs [5–7]) and cortical level (cortical-evoked potentials) [8–10]. Several studies have begun to include the brain state to reduce heterogeneity, especially when looking at MEPs peak-to-peak amplitudes. The research has followed two lines of approaches: offline and online. In the former, acquired data is post-hoc sorted according to oscillations-based biomarkers such as the power spectral density (PSD) [11–16], the phase of the wave [17] or both [7, 18–23]. The second approach, which has been applied only more recently, uses the same biomarkers to choose when to trigger the TMS pulse in real-time [5, 6, 16, 24–28]. In this methodology, to integrate endogenous brain oscillations, the instantaneous brain state in terms of phase and power at the stimulation time must be extracted and predicted from the last acquired information. To estimate power at the stimulation time, the latest acquired data can be used given the relatively low frequency resolution of power [29]. However, to stimulate at a specific phase of an oscillatory period, future prediction is required as the last samples acquired will already be in the past.

Only a few research groups have so far applied such an approach within TMS experiments targeting the motor cortex, however, notwithstanding a brain-state-dependent protocol, initial publications reported controversial results [5, 6, 16, 30]. Although the protocol, intended here as the general experimental design, is an important factor for explaining outcome differences, we believe that the processing pipelines (i.e. signal preprocessing and phase prediction) to forecast the brain state are also important determinants. So far, different research groups have applied different pipelines in terms of filtering procedures and phase forecasting methods. The preprocessing role has already proven essential in functional magnetic resonance imaging analyses [31], and great attention is invested as well in the electroencephalography (EEG) community [32, 33]. For different types of EEG analyses, standardized pipelines have begun to be proposed [20, 34] and analyses studying the effects of choices at different steps of those pipelines are being published [32]. With the goal of creating a standard pipeline for online brain-state-dependent stimulation setups, we analyzed if a specific combination of preprocessing parameters and forecasting algorithms would be outperforming the others in terms of accuracy in phase prediction and computational cost.

To evaluate if a specific combination of preprocessing parameters was more accurate than others, we used a grid-search approach: for each parameter of interest (e.g. filter type, order, etc), a set of values were chosen based on previous knowledge (table 1); then, all possible combinations of values and parameters were created. Each combination was tested on the same dataset. All the filters used were applied in a zero-phase manner to reduce errors deriving from filters phase shifts.

The grid-search was applied on a synthetic dataset and the combinations leading to the most-accurate results were verified on another synthetic dataset. Subsequently, the same combinations were validated on real resting-state EEG data. A similar grid-search was also applied directly to the real dataset to analyze the importance of subject-specifics. The code was run on Python 3.8 and the full scripts can be found at www.github.com/clovbig/grid_search_phase.

2.1. Synthetic dataset

A synthetic dataset was used to have a real ground truth and was considered comparable to EEG [35]. We generated a sinusoid at 10 Hz, to focus on alpha oscillations, on which we added pink noise, i.e. proportional to 1/f. The raw signal for each time point t was given by:

$$\begin{equation*}{\text{raw signal}}\left( t \right) = k \cdot \sin \left( {10 \cdot 2\pi \cdot t} \right) + {\text{pink noise}}\left( t \right)\end{equation*}$$

where the amplitude k of the underlying signal is related to the signal-to-noise ratio (SNR).

Different trials from the same signal were used for verification. For this step, we also increased the complexity of the raw signal by adding sinusoids near the 10 Hz frequency.

2.2. Real dataset

2.3. Phase-forecasting algorithms

2.4. Performance evaluation

2.5. Statistics

The examination of preprocessing hyperparameters for different forecasting algorithms through a grid-search revealed that many combinations can lead to accurate results, with a preference for padding the signal before filtering and then cut the original window to avoid edge effects. Regarding the original window duration, longer ones (i.e. 0.7, 0.75 s) seemed to provide better results. Algorithms-wise, PF-PT [40] and SF-AR [5] worked best, whereas the SF-FFT [37] performed poorly, no matter the preprocessing. Supplementary tables 2(A) and (B) reports the ten best combinations for PF-PT and SF-AR and supplementary figure 2 the effect of each parameter of the grid-search while keeping the other constants. For most forecasting methods, the accuracy increased with the SNR (Spearman's corr = −0.56, p < 0.0001); above an SNR of 0.5 results were not significantly different (figure 1(a)). The verification on synthetic signal showed that addition of near-by sinusoids did not jeopardize the results. Computational cost wise, the SF-AR was the most demanding, taking up to 100 ms, whereas the other approaches were in the units of few milliseconds (figure 1(b)).

Zoom In Zoom Out Reset image size

3.1. Validation on real dataset

The ten-best combinations per algorithm (see supplementary tables 2(A) and (B) for examples) were first validated on healthy resting-state EEG. Phase errors increased for all methods compared to the synthetic dataset and the best approach for the latter was not necessarily the most performant on real data. For subsequent validation, we chose the new best combination and run it on both cohorts. Results are shown in figure 2; the preprocessing parameters chosen are reported in table 2. Due to its poor performance, SF-FFT was not considered for the following analyses.

Zoom In Zoom Out Reset image size

Table 2. Combination of preprocessing performing at best for each forecasting algorithm after verification on synthetic data and initial validation on healthy resting state EEG. Abbreviations used: PF-PT = point-forecasting peak-trough, PF-FFT = point-forecasting using FFT, SF-AR = signal-forecasting with autoregressive model, SF-FFT = signal-forecasting using FFT; RS = stop-band ripple, RP = pass-band ripple (these parameters are necessary only for Chebyshev and elliptic filters).

Algorithm	Filter name	Filter order	Bandwidth (Hz)	Edge cut (%)	Padding	RS	RP	FS	Window (s)
PF-PT [40]	Bessel	2	1	0.3	1	na	na	5000	0.75
PF-FFT	Butter	4	2	0.1	1	na	na	5000	0.75
SF-AR [5]	Butter	3	2	0.1	1	na	na	5000	0.75
SF-FFT [37]	Bessel	2	2	0	1	na	na	5000	0.75

The LMM on the healthy data looking at categorical variables only, reported a significant main effect of channel, algorithm, and eye condition, but of no interaction (supplementary table 3(A)). Regarding the former, Fz was the worst performant (supplementary table 3(C), p < 0.001 for all pair-wise comparisons), whereas other channels did not differ significantly. When looking at differences due to forecasting algorithms, PF-FFT was statistically worse than the other approaches (p < 0.001 for all pair-wise comparisons); moreover, there was a trend of PF-PT being more accurate than SF-AR (EMM = 3.3° or 0.06 radians, p = 0.041, supplementary table 3(D)). When subjects have their eyes closed, the accuracy also seems to increase, but with a negligible effect size (from supplementary table 3(B) we can observe that the EMM is of around 2.3° or 0.04 radians, p = 0.037). Because according to eye condition and channel, a variable that is expected to change is the PSD in alpha (e.g. higher alpha PSD is expected in occipital areas during eyes closed), and because PSD is an important factor when filtering data, we have added to the equation the relative PSD in alpha (i.e. the area under the curve of the power spectra in alpha over the area under the curve in the 1–50 Hz band) and the bandwidth size. Indeed, the relative power in the frequency of interest is a significant factor that may lower the error in phase prediction, as reported by the results of the LMM (p < 0.001, supplementary table 4). In general, it seems that the inverse relation of PSD and phase error holds true independently of the categorical factors studied (i.e. channel, eye condition, algorithm, and bandwidth used for computing the PSD) as can be extrapolated for supplementary figure 3. This point becomes even more evident when taking as factor the absolute alpha PSD divided into quartiles per subject. By taking as an example PF-PT on young healthy subjects, we observe that if only very high PSD trials are considered (i.e. only PSD above the 75th percentile for the different subjects), the interquartile difference is decreased to less than 45°, compared to almost 90° when all trials are considered (supplementary figure 4(b)). Nonetheless, some outliers are still present, probably due to the subject, channel, and eye condition variabilities. Similar effects were seen in the stroke dataset (supplementary tables 5–9) and the two populations had overall similar accuracies when the same approach was applied in the same conditions (figure 2).

Zoom In Zoom Out Reset image size

3.2. Grid-search on real data

3.3. Comparison with original algorithms' preprocessing

Protocol differences, such as stimulation parameters and other experimental design choices, are often held responsible for heterogeneity in brain-state-dependent experiments' results. The present work highlights the relevance of signal processing, often arbitrarily decided by the researcher. Indeed, different choices in EEG filtering can lead to importantly diverse results in phase-forecasting accuracy (supplementary figure 5). Two main points should be emphasized: (a) outcome comparisons with different processing pipelines must be analyzed with care as applied functions on the raw signal can change its look at a very early stage of the data analysis; (b) no matter the forecasting method chosen, optimization of preprocessing is critical and can significantly improve precision (figure 3(a)).

Our grid-search approach for finding a suited preprocessing combined background knowledge and empirical data. All the combinations were valid, but tests on synthetic and real data rejected some while keeping others. In particular, this work highlighted important steps that should be followed in the EEG preprocessing: padding the signal before filtering is probably the most important step to perform at the beginning to avoid most of the edge effects that will be introduced by the filter itself. Regarding the latter, infinite impulse response filters with low-order are to be preferred; bandwidth-wise narrower bands were suggested by the grid-search, although wider ones may still be tested on real data. In terms of sampling frequency, this work suggests one between 1 and 5 kHz. Finally, algorithm-wise FFT-based ones (both point and signal forecasting) had the lowest performance, which may be due to the intrinsic stationarity of FFT, from which minimal discrepancies from the ground-truth can lead to large errors when moving away from the initial window. Rejection of tested combinations was based on accuracy, but also on computation cost given that time is a main constrain in online set-ups. Except for the SF-AR, forecasting approaches took less than 5 ms, a good value that accounts for only a 25th of cycle in the alpha band (and maximum a 6th of cycle in the highest beta band). The present analyses did not aim at optimizing the speed of processing, but rather looked at the cost of the different pipelines when the same level of code optimization is used. Regarding phase-forecasting precision, results from this grid-search showed that newly found preprocessing combinations significantly increased the accuracy of some of the published forecasting methods in both the healthy and the stroke cohort data. Outcomes from the latter were indeed very positive: to our knowledge this is the first time that phase-forecasting algorithms have been tested on such a large dataset of stroke patients in different time points after the ictal accident. The achieved accuracy, comparable to that of healthy subjects, implies that these algorithms can be applied even in situations with impaired brain functions, such as after a stroke. Moreover, we have observed that personalizing a grid-search for single subjects does not significantly improve the performance. Suggesting that the optimized method found works similarly on different subjects. Altogether, this provides an important step forward using state-dependent TMS-approaches as interventional tools to e.g. enhance the effects of neurorehabilitation.

In addition to highlight some important aspects and effects of EEG preprocessing, the goal of this study was to find optimized parameters for phase forecasting. The synthetic grid-search returned similar filters for different prediction algorithms that were also confirmed and strengthen by the individual grid-search on real data. Another relevant outcome of the analysis is the importance of the PSD, which can be read on two levels: firstly, for a filter and subsequently for a forecasting algorithm to perform well, a relatively strong signal (and SNR) is required (supplementary figures 2 and 3); secondly, the PSD strength can help defining if a neural oscillation is actually present. In this study PSD was only analyzed as a factor on accuracy as our aim was to determine the performance of the respective approaches with unbiased/unselected real data. However, given the strong effect PSD has on accuracy, we believe that in future steps towards online experiments, PSD may be used to define if a trial is eligible, as suggested by [42].

Finally, we acknowledge that other parameters and forecasting approaches (e.g. machine-learning based [43]) could have been added to the grid-search, but were beyond the scope of the present study. The grid-search approach as an optimization tool is one of several others available. We believe it is a good starting point due to its simplicity and when most of the parameters studied are categorical, rather than continuous, and are based on theoretical knowledge.

As future perspective, we believe that original and optimized approaches should be applied on the same TMS-EEG-EMG dataset and evaluate the actual effect of pre-processing and thus phase labeling on outcomes such as the MEP amplitude.

The present grid-search results underline the key role of EEG preprocessing steps and demonstrate that this methodology can be a tool for optimizing signal processing pipelines crucially important for brain-state-dependent or closed-loop NIBS approaches. Although the analyses focused on phase forecasting for brain-state-dependent stimulation experiments, the experience from it can be extended to other EEG analytical approaches. Finally, brain-state-dependent approaches were successfully tested on resting-state EEG data from patients after a stroke, paving the way to being able to determine critical factors relevant for brain-state-dependent novel interventional approaches even in patients with brain lesions.

We acknowledge access to the facilities and expertise of the Neuromodulation facilities of the Human Neuroscience Platform of the Fondation Campus Biotech Geneva.

The scripts used for these analyses are publicly available and can be found at www.github.com/clovbig/grid_search_phase. The synthetic dataset can be regenerated from the scripts; while the real dataset used for this manuscript can be provided upon reasonable request.

The research was partially funded by the Wyss Center for Bio and Neuroengineering (AVANCER WCP030; Genève, Switzerland) to FCH; the Strategic Focus Area 'Personalized Health and Related Technologies (PHRT, #2017-205)' of the ETH Domain (CH) to FCH; the Defitech Foundation (Morges, CH) to FCH.

C B conceptualized the method, wrote the scripts for the grid-search and different EEG preprocessing and performed the statistical analysis. A C M and T M acquired the data. F C H and T M supervised the work. All authors read and the approved the final manuscript.

The authors declare that they have no competing interests.

All subjects from whom the data was acquired signed a written informed consent and the study was approved by the ethics committee of the Canton of Vaud, Switzerland (No. 2018-01355).