The shaky ground truth of real-time phase estimation.

Instantaneous phase of brain oscillations in electroencephalography (EEG) is a measure of brain state that is relevant to neuronal processing and modulates evoked responses. However, determining phase at the time of a stimulus with standard signal processing methods is not possible due to the stimulus artifact masking the future part of the signal. Here, we quantify the degree to which signal-to-noise ratio and instantaneous amplitude of the signal affect the variance of phase estimation error and the precision with which "ground truth" phase is even defined, using both the variance of equivalent estimators and realistic simulated EEG data with known synthetic phase. Necessary experimental conditions are specified in which pre-stimulus phase estimation is meaningfully possible based on instantaneous amplitude and signal-to-noise ratio of the oscillation of interest. An open source toolbox is made available for causal (using pre-stimulus signal only) phase estimation along with a EEG dataset consisting of recordings from 140 participants and a best practices workflow for algorithm optimization and benchmarking. As an illustration, post-hoc sorting of open-loop transcranial magnetic stimulation (TMS) trials according to pre-stimulus sensorimotor μ-rhythm phase is performed to demonstrate modulation of corticospinal excitability, as indexed by the amplitude of motor evoked potentials.

Introduction

Oscillatory activity at different frequencies is a prominent feature of EEG recordings of brain activity (Buzsaki and Draguhn, 2004). The functional role of brain oscillations is demonstrated in time–frequency analysis of evoked EEG activity, averaged over many trials, showing brain region specific changes in spectral power with different brain states, such as those related to visual attention, memory retention, and motor behaviour (Pfurtscheller and Lopes da Silva, 1999). Here, we focus on single-trial single time-point analysis, where oscillations can be characterized according to instantaneous amplitude (Freeman, 2004a) and instantaneous phase (Freeman, 2004b) using the Hilbert transform (Freeman, 2007) or, equivalently, Fourier or wavelet based analysis (Bruns, 2004). The functional relevance of investigating single-trial instantaneous phase is motivated by effects of cortical excitability (Bergmann et al., 2012; Massimini et al., 2003; Thies et al., 2018; Zrenner et al., 2018) and sensory threshold (Ai and Ro, 2014) and by its correlation with the rhythmic neuronal activity of different circuits (Haegens et al., 2011; Miller et al., 2012).

Standard signal processing methods require data before and after the time of interest to determine instantaneous phase and amplitude of a signal and are therefore “non-causal” (Fig. 1, bracket labeled “2”). In the case of real-time phase estimation, or in post-hoc estimation of phase at the time of a stimulus that affects the signal, such as a TMS pulse but also in EEG potentials evoked from sensory stimulation (Dugue et al., 2011), only data preceding the time of interest is available (Fig. 1, bracket labeled “3”). In the case of causal phase estimation, we aim to construct the best estimator of phase, approximating the measure that we would have obtained with a non-causal estimator from the whole signal (Fig. 1, bracket labeled “2”), if this were available, but using only data preceding the time point of interest.

Fig. 1

Measures of phase, with the phase of interest at the center of an epoch of data (dotted red line). (A) The oscillation of interest is not normally accessible as a “clean” signal, except in the case of simulated data with a synthetic sinusoids, in which case phase is known by definition (labeled 1, “true phase”). (B) 1/f noise, adds irrelevant information to the signal of interest reducing estimator precision. The relative spectral amplitude at the frequency of interest and the amplitude of the sinusoid determines the signal-to-noise ratio. (C) The recorded signal contains the oscillation of interest and noise. Phase of a signal not containing a stimulation artifact can be estimated offline from data before and after the timepoint of interest using band-pass filtering (labeled 2, “non-causal phase estimate”), or causally, from the segment of data preceding the timepoint of interest only with a real-time algorithm (labeled 3, “causal phase estimate”).

Except for the case of a synthetic signal, the only available benchmark for any given causal estimator is simply non-causally estimated phase. But how meaningful is this benchmark? Given that an EEG recording is the result of global brain activity, any target oscillation of interest is affected by other activity in the form of noise. In the limit, when there is no oscillation of a given frequency present (no spectral peak above background noise), the phase value determined by band-pass filtering and Hilbert transformation is meaningless and the filtering process may even introducing spurious results (Widmann et al., 2015). With regard to the theoretical minimum error variance of any estimator, this is expressed as the Cramér–Rao Lower Bound. If the properties of the signal and noise are known, it can be derived analytically that error variance is inversely proportional to signal-to-noise ratio (SNR) (Peleg and Porat, 1991).

Here, we investigate the practical limits of phase estimation methods in general and how phase estimation error variance is affected by spectral properties and single trial parameters. We propose a framework to evaluate the accuracy of causal (using data preceding the timepoint of interest only) phase estimation methods under different experimental conditions and provide an open source software toolkit for causal estimation of instantaneous phase, PHASTIMATE, implementing an approach described previously (Chen et al., 2013; Zrenner et al., 2018). Method parameters are optimized for the extraction of sensorimotor μ-rhythm using a genetic algorithm. The utility of the scripts is demonstrated by post-hoc sorting single trials of a dataset of motor evoked potentials (MEP) according to pre-stimulus EEG μ-rhythm phase.

Methods

Subjects and data

This study analyses data from two different experiments. All experiments were performed after obtaining written informed consent from the participants in accordance with The Code of Ethics of the World Medical Association (Declaration of Helsinki) and having obtained approval from the ethics committee of the University of Tübingen (application 716/2014BO2).

To investigate of phase estimation accuracy, resting-state EEG is analyzed from a total of 140 participants (53 male, 87 female, age 24 ​± ​5 years, all right-handed) that were screened for EEG–TMS experiments between April 2018 and May 2019. There were no exclusion criteria based on EEG. Four datasets were removed because no spectral peak could be determined between 8 and 14 ​Hz, four additional datasets were removed as a negative SNR was obtained after fitting and subtracting 1/f noise, indicating that the noise removal failed (see supplementary data for spectra of included and excluded subjects), yielding 132 included datasets.

For the post-hoc sorting of MEPs according to pre-stimulus EEG phase, data from a separate experiment is analyzed. In this experiment, 1150 TMS pulses were applied to the hand-knob area of the left motor cortex of one healthy volunteer while simultaneously recording EEG as well as evoked motor responses through electromyography (EMG) of the right abductor pollicis brevis muscle.

EEG data was recorded with 64 channel EEG caps (Easycap GmbH, Germany) and 24 bit EEG amplifiers (NeurOne Tesla, Bittium, Finland) in DC mode at a sampling rate of 5 ​kHz (resting-state EEG experiment) and 10 ​kHz (MEPs experiment), respectively. In the MEP experiment, EMG data was recorded using the bipolar input channels of the EEG amplifier, EMG data was not downsampled. EEG data was spatially filtered with a C3-centered Hjorth-style (Hjorth, 1975) Laplacian (using FC1, FC5, CP1, and CP5 as surrounding electrodes) to isolate sensorimotor μ-rhythm, and down-sampled to 1 ​kHz after application of a low-pass anti-aliasing filter. The first 250 ​s of resting-state EEG data was extracted and stored in a raw data file. In the post-hoc sorting of MEPs experiment, biphasic TMS pulses (160 μs duration) were administered using a figure-of-eight coil (PMD70-pCool coil with PowerMAG Research 100, MAG & More GmbH, Germany) oriented perpendicular to the left precentral gyrus with the second phase of the induced electric field in the posterior-anterior direction. The stimulation intensity was 115% of resting motor threshold (i.e., the stimulation intensity necessary to evoke MEPs exceeding 50 ​μV amplitude with 50% probability), the inter-stimulus interval was 2.1 ​s on average with a maximum jitter of ±0.1 ​s.

Summary of data analysis pipeline

Spatial filtering and extraction of a 250 ​s resting-state EEG segment yielded a one-dimensional data vector for each subject, which was concatenated to form a 140 ​× ​250,000 (subject ​× ​sample) raw data matrix (available for download). The subsequent analysis followed the following steps, as detailed below:

Spectral analysis was performed to determine the peak frequency and SNR of sensorimotor μ-rhythm for each data record.

Each data record was divided into 500 overlapping epochs 2 ​s long. Phase was determined at the center of each epoch non-causally (from the whole epoch data) using different band-pass filter designs followed by the Hilbert transform, and the uncertainty of this estimate was quantified.

Phase was then estimated based on a window of data preceding the center of each epoch only, using the PHASTIMATE method. For each subject, parameters of the estimation method were optimized using a genetic algorithm minimizing the circular deviation of the difference between the causal and the non-causal phase estimate (across the 500 epochs, within each subject).

The PHASTIMATE script was applied again with different algorithm parameter sets: parameters used previously (Zrenner et al., 2018), the average of the optimized parameters across subjects, and the parameters optimized individually for each subject.

Phase accuracy was assessed under different conditions, specifically SNR (as a subject-by-subject parameter) and fluctuations in instantaneous amplitude (as an epoch-by-epoch parameter).

Finally, the PHASTIMATE script is applied post-hoc to EEG–TMS–EMG data to reproduce a previously reported relation between the pre-stimulus μ-rhythm phase and the MEP amplitude. MATLAB code of implementing the above analysis steps is available for download, with dependencies for the signal processing toolbox and the optimization toolbox, tested in Matlab R2017b (The MathWorks, Inc., MA, USA), and generating all the figures included in this manuscript.

Spectral analysis

Spectral analysis was performed implementing the Welch method: Data was segmented into 50% overlapping epochs of 2 ​s duration, which were linearly detrended (de-mean), Hann-windowed, Fourier-transformed, and averaged, resulting in an amplitude spectrum with 0.5 ​Hz resolution. Peak frequency was determined in the range between 8 and 14 ​Hz, 1/f noise was estimated by fitting a straight line to the log–log spectrum at fixed frequencies outside of known oscillations (0.5–7 ​Hz and 35–65 ​Hz) in the least-squares sense. SNR was defined as the difference between peak spectral amplitude and fitted noise at that frequency, on the log scale, in units of dB.

Non-causal phase estimation

To assess the performance of phase estimation at the edge of a window of data, a “benchmark” phase value is required. To this end, overlapping 2 ​s epochs of data were generated and phase was then determined at the center of each epoch by applying a band-pass filter followed by the Hilbert transform. In order to assess the reliability of the phase estimate, a family of equivalent estimators was used (Sameni and Seraj, 2017). Specifically, we defined 15 band-pass-filter designs (Table 1), consisting of seven finite impulse response (FIR) filters and eight infinite impulse response (IIR) filters. Passbands were centered on the individual peak frequency resulting in a magnitude response after zero-phase (forward and backward) filtering as shown in Fig. 3, panel A. As each of these designs could be used to determine a benchmark phase value with equal justification, we take the circular mean of the different phase estimates. A high variance between the estimates would indicate that it is not possible to determine a meaningful phase estimate for a given epoch, even with non-causal methods.

Table 1

Summary of the filters applied to estimate phase with standard non-causal signal processing methods.

Number of filters	Filter type	Design parameters
4	FIR	Windowed sinc with order 2,3,4, and 5 times the peak frequency period
3	FIR	Least squares with order 3,4, and 5 times the peak frequency period
3	IIR	Butterworth with order 4, 8, and 12
3	IIR	Chebyshev Type I with order 4, 6, and 8
2	IIR	Elliptic filter with 20-dB and 40-dB attenuation

Fig. 3

Non-causal phase estimation error with synthetic EEG data. (A) Magnitude responses of a set of 15 band-pass filter designs after forward and backward filtering used for phase estimation. (B) Analysis of simulated EEG data at 11 different levels of SNR between 4 and 22 ​dB, showing both the circular deviation of the 15 estimators (blue), and also the corresponding median absolute phase error of each estimator (dashed grey lines) as well as the median phase error when taking the average of all estimators (red line). Note that the phase error can only be determined because the data is synthesized using a sinusoid of known phase embedded in simulated 1/f noise; a similar analysis is not possible with physiological EEG data.

Synthetic EEG data generation

To assess the efficacy of this “family of equivalent estimators” approach, EEG data was simulated using synthetic sinusoids with known phase with additive 1/f noise. Such realistic background EEG noise was generated using the simulated EEG data generator toolbox1 which implements the method of summing sinusoids of randomly varying frequency and phase and scaling the amplitude of the sinusoid at each frequency to match a physiological 1/f EEG power spectrum (Yeung et al., 2004). The amplitude of the added 10 ​Hz sinusoid with known randomized phase was scaled to achieve different SNR levels between 4 and 22 ​dB, 1000 epochs each 1.5 ​s long were generated for each SNR level for subsequent non-causal phase estimation using the “family of equivalent estimators” described above.

Causal phase estimation

The following is a brief description of the autoregressive forward prediction approach for phase estimation (Chen et al., 2013) that was applied in this study and is implemented by PHASTIMATE (Fig. 2). A window of data extracted is band-pass filtered (forward and backward, resulting in zero phase shift) with edges removed. The size of the data window is a free parameter of the algorithm, typically containing 3–8 cycles of the oscillation of interest (for stable oscillations without phase resetting, a longer window would be expected to perform better, for oscillations with variable phase or phase reset, a shorter window would be expected to perform better). The signal is then extended into the future using autoregressive parameter estimation (Yule–Walker method) and forward prediction to encompass the time point of interest. The Hilbert transform is then applied to the resulting data segment yielding the complex analytic signal, the angle of which corresponds to phase.

Fig. 2

Autoregressive forward prediction method for phase estimation. (A) A window of data extracted from a Laplacian montage centered at the C3 electrode is band-pass filtered. (B) The edges are removed and autoregressive parameters are determined. (C) The signal is extended into the future to encompass the time point of interest. (D) Hilbert transform is applied to yield the complex analytic signal, and thereby instantaneous phase and amplitude at the time point of interest (Figure taken from Zrenner et al., 2018).

Algorithm parameter optimization

A genetic algorithm was used to determine individually optimal values of the four parameters of the PHASTIMATE algorithm listed in Table 2 that minimized the circular deviation of the difference between the causal PHASTIMATE estimate and the non-causal benchmark, determined as described above. The parameters were constrained to be integers, phase was estimated from 500 overlapping epochs, and a genetic algorithm population size of 100 competing parameter sets in each generation was used with the optimization bounds shown in Table 2 and with the additional constraint that the window size (in samples) needed to be at least three times the filter order, the sampling rate being 1 ​kHz.

Table 2

Bounds within which the genetic algorithm could optimize the parameters of the PHASTIMATE algorithm. The parameter values used previously in Zrenner et al. (2018), are shown for comparison. These parameters apply to a sample rate of 1 ​kHz.

	Parameter values in Zrenner et al. (2018)	Optimization bounds
Window size (samples)	500	400 .. 750
Filter order	128	100 .. 250
Number of samples removed at edge	64	30 .. 120
Order of the autoregressive model	30	5 .. 60

PHASTIMATE toolbox

PHASTIMATE is available as an open source toolbox and implements the steps laid out in this report as a best-practice approach for experimental studies investigating relationships between phase of an EEG signal and evoked responses. (1) As a first step, the accuracy at which phase can be determined with standard non-causal signal-processing methods is estimated from EEG data without stimulation artifacts using the family of equivalent estimators approach (Sameni and Seraj, 2017). The function phastimate_truephase.m calculates estimator variance and also generates the benchmark data required for Step 2. (2) In the second step, the performance of the predictive estimator of choice is then assessed from the same data. Algorithm parameters can be optimized with phastimate_optimize.m using a genetic algorithm to minimize phase error variance. (3) Steps 1 and 2 assess the accuracy at which phase can be determined in principle and in practice in a given dataset. The predictive algorithm implemented in the script phastimate.m is then used for trial sorting, optionally with the optimized parameter set determined in Step 2 and only considering epochs where instantaneous amplitude is high. The entire procedure is demonstrated in a main script (main_script.m) which runs through the analysis of the supplied sensorimotor dataset performing the steps described above.

Circular statistics

Circular statistics formulas were adapted from the CircStat toolbox (Berens, 2009), with circular variance used as a measure of circular spread that is bounded between 0 and 1 (1-R, where R is the magnitude of the resultant vector averaging all the normalized moments, bounded between 0 and 1, and also known as phase locking value). Estimation error deviation is reported as circular standard deviation () rather than precision (1/variance) or phase locking value to ease interpretation.

Post-hoc sorting of trials according to pre-stimulus phase

EEG data in the 1 ​s preceding each TMS pulse was extracted, spatially filtered, and down-sampled to 1 ​kHz. MEP amplitude was determined as the peak to peak of the EMG recording in the period between 25 and 45 ​ms post-stimulus. For the sake of simplicity, no trials were discarded from the 1150 available stimuli, neither based on EEG artifact criteria nor based on EMG criteria such as pre-innervation.

We analyzed whether the MEP amplitudes were related to the μ-rhythm phase at the time of the stimulation by applying circular-to-linear regression analysis (Cox, 2006; Cremers et al., 2018; Kempter et al., 2012), which is a sensitive method to find such relations (Zoefel et al., 2019). The fitted regression model had the form a + b cos(φ) ​+ ​c sin(φ), where a, b, and c are the model parameters and φ is the phase. Prior to the analysis, the MEP amplitudes were log-transformed to reduce the skewness of the distribution. The forward-prediction phase-estimation algorithm was used as implemented by PHASTIMATE on data down-sampled to 1 ​kHz with the parameters derived from the optimization of the resting state EEG data (window length of 719 ​ms, FIR filter with order 192, fixed pass-band of 8–13 ​Hz, edge 65 samples, autoregressive model order 25, data segment length for Hilbert transform 128 samples). A significance value was obtained from an F-statistic comparing regression model fit to a constant model. Power analysis was performed at a significance value of 0.05 for between 20 and 200 samples using 2000 repetitions of bootstrapping without replacement from the available 1150 trials.

Results

Synthetic EEG data analysis

The simulated sinusoids in 1/f noise dataset with known “true” phase enabled the result of the family of equivalent band-pass filter based estimators to be related to phase estimation error even at very low signal to noise ratios. When the different estimators are in agreement (low variance), the median absolute phase error is low, which is the case when SNR is high; when the different equivalent estimators are not in agreement (high variance), the actual phase error of each individual estimator is high, which is the case when the frequency of interest has a low SNR (Fig. 3, panel B).

For this dataset, at a SNR of 6 ​dB with a circular deviation among estimators of 10°, this still corresponded to a median absolute error of the non-causal “benchmark” measure (circular mean phase of equivalent estimators) with respect to true (synthetic) phase of 22° (Fig. 3 panel B). Due to artifacts and noise characteristics in a given epoch of data affecting all estimators similarly, a high correspondence among estimators does not necessarily indicate accuracy – the estimators can all be ‘wrong’ in the same way.

Limits of phase determination

In accordance with the case of sinusoids with additive noise, where estimator precision is proportional to SNR (Peleg and Porat, 1991), in the dataset considered in this study, median circular deviation of non-causal phase estimates decreased from 17.5° to 6.5° across participants as SNR increased from 0 to 20 ​dB (Fig. 4, panel A) with a corresponding Pearson correlation coefficient of R2 ​= ​0.74 (p ​= ​10−41). Within participants, periods of low vs. high μ-rhythm amplitude increased estimator circular deviation by ~20–30° between the lowest and highest amplitude quartile, depending on SNR (Fig. 4, panel B).

Fig. 4

Non-causal phase estimator variance with real EEG data. (A) Median circular deviation of “non-causal phase” estimates across 500 epochs of data and the SNR for each of the 132 subjects. The line represents a linear fit to the data. (B) Median circular deviation for each participant in four subgroups of 125 epochs sorted according to the estimated μ-rhythm amplitude (the blue and purple data points correspond to the group of the highest and lowest μ-rhythm amplitude, respectively). The lines represent linear fits to the data of the subgroups.

Optimized real-time phase estimators

Instantaneous oscillatory phase of an EEG signal is not defined to arbitrary precision and the above non-causal analysis may serve as both an indication of the theoretical limit, as well as a benchmark against which causal estimators that only have data before the time point of interest can be optimized. These predictive algorithms can be employed in real-time applications or when a stimulus artifact renders the signal after the time of interest unusable.

The optimization of the PHASTIMATE parameters with a genetic algorithm yielded an average optimized window size of 719 samples, filter order 192, edge 65, and order of the autoregressive forward prediction of 25. The main difference between the parameters used previously (Zrenner et al., 2018) and the results of the current analysis is the longer window size and higher filter order (see Table 3). Note that we also had an individualized filter passband of ±1 ​Hz around the individual peak frequency as opposed to the fixed 8–13 ​Hz passband used by Zrenner et al. (2018).

Table 3

Optimized phase-estimation parameter values in 132 participants. The median (and the interquartile range) of the individual optimal parameter values for window size, filter order, number of samples removed at each edge, and the order of the autoregressive model used for forward prediction are given along with the parameter values used by Zrenner et al. (2018), the sampling rate being 1 ​kHz.

	Parameter values used by Zrenner et al. (2018)	Optimized paramater values
Window size (samples)	500	719 (693–740)
Filter order	128	192 (158–229)
Number of samples removed at edge	64	65 (53–81)
Order of the autoregressive model	30	25 (17–34)

The results of running PHASTIMATE with the two parameter sets of Table 3 are shown in Fig. 5, both with filters having a fixed pass-band of 8–13 ​Hz and using filters with a ±1-Hz pass-band around the individual peak frequency. Increasing the window length of the data under consideration from 500 to 719 ​ms (and changing the other parameters according to Table 3) resulted in a median reduction of phase estimation error deviation of ca. 5° (Fig. 5, panel B), the use of filters with individualized passbands resulted in a modest improvement only when the signal had a low SNR (Fig. 5, panel C). In comparison with the parameter values used previously, the individually optimized parameters yielded an average reduction of ca. 10–15° in error deviation.

Fig. 5

Precision of the causal phase-estimation algorithm. (A) The median circular deviation of the phase estimates obtained with different autoregressive forward-prediction algorithm parameter sets and of the corresponding phase estimates across 500 epochs of data and the SNR for each of the 132 subjects. The lines represents linear fits to the data. The data correspond to a fixed filter with a passband of 8–13 ​Hz (solid line, markers visible) and an individualized filter with a passband of ±1 ​Hz around the individual peak frequency (dotted line, markers not shown). (B) Histogram showing the phase estimation accuracy improvement when changing from the 500 ​ms window to the 719 ​ms window parameter set (see Table 3) for both fixed 8–13 ​Hz and individualized passband filters. (C) Histogram showing the phase estimation accuracy improvement when changing from fixed 8–13 ​Hz to individualized passband filters for both the 500 ​ms window and the 719 ​ms window parameter set.

As was the case for the spread of non-causal phase estimates (Fig. 4), the estimation error of the predictive phase estimation method implemented in PHASTIMATE was strongly affected by instantaneous amplitude (Fig. 6): The quartile of epochs with the lowest μ-rhythm amplitude had a median error deviation of ~70–100° depending on SNR, whereas the quartile of epochs with the highest μ-rhythm amplitude had a median error deviation of only ~20–50°.

Fig. 6

The median circular deviation of the phase estimates obtained with group optimized algorithm parameters (window length: 719 ​ms, filter order: 192, pass-band: individual peak frequency ±1 ​Hz, edge: 64 samples, autoregressive model order: 30, at a sample rate of 1 ​kHz) for four subgroups of 125 epochs sorted according to the estimated μ-rhythm amplitude and the SNR for each of the 132 subjects. The lines represent linear fits to the data of the subgroups. The plots on the right visualize the overall error distribution including all trials in all subjects within the respective amplitude subgroup.

MEP size is significantly larger when TMS is applied at the negative peak of the oscillation vs. at the positive peak (Zrenner et al., 2018) in a real-time closed-loop experiment. Here, the same question was addressed through post-hoc sorting of trials, where TMS pulses were applied open-loop, and therefore having a random μ-rhythm phase at the time of the stimulus. Circular-to-linear regression analysis between the forward-predicted phase as estimated with PHASTIMATE and the log-transformed MEP amplitudes showed a highly significant correlation with the sinusoidal model having a significantly better fit than a constant model (p ​= ​10−32), and with the largest evoked responses coinciding approximately to the trough of C3-centered μ-rhythm phase (Fig. 7, panel A). A power analysis was performed (Fig. 7, panel B) showing that 75 trials would be the minimum required in this subject to demonstrate a phase effect (power 80%, alpha 0.05).

Fig. 7

Results of the post-hoc analysis. (A) MEP amplitude as a function of estimated phase of sensorimotor μ-rhythm as extracted by a C3 Hjorth-style Laplacian at the time of stimulation (note the logarithmic amplitude axis). The solid sinusoid shows the result of a circular regression fit. The dotted sinusoid shows the pre-stimulus phase corresponding to the horizontal axis. The overlaid boxplots show the median and interquartile range of the data corresponding to 36°-wide bins. Below, the average pre-stimulus EEG signal in a 150 ​ms time window preceding the stimulus is shown for each phase bin, with the post-stimulus artifact as a grey rectangle. (B) Power analysis at a significance value of 0.05 for between 20 and 200 samples using 2000 repetitions of bootstrapping without replacement from the available 1150 trials. The dashed line indicates that approx. 75 trials are required in this subject for an 80% likelihood of showing a significant phase effect.

Discussion

Relevance

An estimate of phase of a given oscillation, extracted by surface EEG, at any particular point of time, is just a scalar in the range from −π to +π. Nevertheless, under the right conditions, this very crude measure of instantaneous brain state predicts evoked responses (Fig. 7). In this study, we considered the limits within which such a scalar state value may be a meaningful measure of some aspect of brain state. Especially when multiple spatially distributed measures of oscillatory phase serve as the basis to derive more complex metrics such as cortical waves (Alexander et al., 2013, 2015) or connectivity state (Stefanou et al., 2018), it is important to understand the precision at which a phase measure can be theoretically and practically determined given various conditions.

An estimate of oscillatory phase in EEG is also just an estimate of a signal parameter in noise. A perhaps banal observation is that a signal parameter can only be estimated if the signal is present and distinguishable from the noise within which it is embedded. However, as the analysis of the accuracy of non-causal phase estimates in synthetic EEG data (Fig. 3, panel B) and the analysis of the variance of equivalent estimators in real EEG data (Fig. 4, panel B) show, for datasets with low SNR and during periods of low amplitude, even a non-causal phase estimate cannot be accurately determined, let alone be meaningfully approximated with causal methods that use only data preceding the timepoint of interest.

What does this mean for the design of closed-loop EEG studies? One possibility is to include only subjects where the chosen spatial EEG filter (in case of the sensorimotor μ-rhythm, typically a C3-centered Hjorth-style Laplacian) yields a sufficient SNR to enable phase to be targeted accurately. Alternatively, various approaches exist to calculate optimized individual spatial filters based on anatomy and a dipole of interest (e.g., linearly constrained minimum variance beamforming) but with mixed success (Madsen et al., 2019), based on behavioural tasks such as motor imagery or fist clenching (e.g., common spatial patterns), or based on spectral signal properties (e.g., spatial–spectral decomposition (Nikulin et al., 2011; Schaworonkow et al., 2018), which may enable a signal of interest to be extracted at sufficient SNR for phase detection in a larger proportion of subjects.

Implications

However, given that the accuracy with which phase can be determined even in principle varies more strongly within subjects than between subjects (Fig. 4, panel B), selecting the right moment to stimulate may be more important than selecting the right participant. Since oscillatory power is task-dependent, another promising approach may be to target phase during a task that amplifies the oscillation of interest (e.g., the relaxation of a clenched fist for the sensorimotor μ-rhythm). The power of an oscillation of interest may be increased also through pharmacological intervention (e.g., the GABA reuptake inhibitor tiagabine increases oscillations in the alpha and beta band (Darmani et al., 2019).

An important constraint that arises from the decrease of phase targeting accuracy during periods of low amplitude is that phase and amplitude cannot be investigated independently. A comparison of phase-dependent evoked responses at high amplitude with responses at low amplitude would be confounded by the phase estimation accuracy being significantly lower at low amplitude. When targeting a specific phase in real-time, accuracy is improved by including an instantaneous amplitude threshold in the trigger condition, but this will also lengthen the inter-trigger-interval which will be influenced by slow fluctuations in spectral amplitude thus providing indirect amplitude-based feedback to the participant.

The purpose of this report is to provide a framework to verify that phase can be meaningfully targeted within a given experimental condition. We make the PHASTIMATE scripts available as an open-source toolbox to determine the limits of phase estimation and to facilitate the design of experiments that optimize conditions for phase targeting. In terms of best practice, we recommend the procedure laid out in this manuscript as implemented in the example script: (1) Spectral analysis enables the quality of the signal to be checked and the SNR to be estimated, including the reporting of single subject spectra. (2) Variance of phase estimates should be checked and reported (Sameni and Seraj, 2017) to establish the reliability of a non-causal estimator. (3) The performance of a predictive (causal, pre-stimulus data only) phase estimator can then be tested and optimized using EEG data without stimulation artifacts. Care should be taken that the properties of the EEG signal used to optimize the algorithm matches the properties of the signal used for real-time triggering or post-hoc sorting (e.g., if an amplitude threshold is used, amplitude should be matched in both datasets).

If the above tests yield satisfactory error variances, this provides confidence that the selected experimental conditions (frequency of interest, spatial filter, behavioural task) allow for meaningful phase targeting. Phase effects at the single-trial level can then be investigated using a post-hoc sorting approach as in the demonstration code generating Fig. 7. In terms of quantitative targets, based on the results of the synthetic EEG data with the parameters used in this study, an SNR of ≥9 ​dB would indicate a median error of 15° in the non-causal phase estimate, which would typically be acceptable in a benchmark for the causal estimator.

One further possible complication is the presence of multiple neighboring spectral peaks within the frequency range of interest. This indicates that a superposition of multiple oscillators from different sources (for example, a 10 ​Hz peak originating in visual cortex, and a 12 ​Hz peak originating in sensory cortex) is present in the signal, which results in a compound phase value that depends on the relative amplitudes and is not related to either signal. Here, it may be possible to design a EEG montage (spatial filter) that is sensitive to the cortical region of interest while attenuating specific regions generating an interfering oscillation (Hauk and Stenroos, 2014), or to perform a behavioural intervention (e.g. opening or closing of eyes).

Limitations

Although we believe the specific conceptual framework and the PHASTIMATE code to be generally applicable to targeting brain oscillations with EEG, the specific parameter optimizations presented here do not necessarily generalize beyond the sensorimotor μ-rhythm extracted with a C3-centered Hjorth-style Laplacian. Longer time windows are not expected to be generally better, especially if the oscillation of interest is not stable. It would in any case be advisable to rerun the optimization step with the specific signals of interest.

We also did not rigorously compare the effect of different filter designs on PHASTIMATE performance and only report results for simple FIR filters generated using a windowed sinc design method. Elsewhere, elliptic IIR filters have been reported to perform well, but they did not improve accuracy in initial tests using PHASTIMATE. It is nevertheless possible that IIR filter designs exist that outperform the FIR filter used in our study. The bound for the maximum window length was set to 750 ​ms; for many participants the genetic algorithm resulted in a window length at the maximum bound. We nevertheless decided against increasing this further since a longer time window would constrain the minimum inter trigger interval.

Furthermore, our spectral analysis method is comparatively simple, using the Welch method to generate the periodogram and fitting 1/f noise using a straight line in the log–log scale at frequencies where no physiological oscillations are expected. Due to the 2 ​s windows, our peak frequency estimate has a 0.5 ​Hz resolution. Using more advanced spectral analysis methods (such as the multi-taper method) and parametric methods for fitting 1/f noise such as the algorithm implemented by the “fitting of one-over-f” (FOOOF) toolbox (Haller et al., 2018) may result in better peak frequency and SNR estimates but would likely not change the qualitative results of this report.

In terms of the expected benefit of the genetic algorithm-based individual optimization, it should be noted that we are using the same dataset for calibration and testing, which will likely overestimate the benefit of the individual algorithm calibration. However, since we argue that fluctuating signal properties (instantaneous amplitude, SNR, slow drifts, presence of artifacts) have a far more important effect on the accuracy of the causal estimate as compared to a tweaking of the algorithm parameters, we believe that this analysis, which represents a “hypothetical best case” is acceptable.

Finally, we did not consider the effect of eye blinks, muscle, or movement artifacts in EEG beyond excluding a small number of subjects with spectra containing noise to a degree where no μ-rhythm peak could be found or 1/f noise fitting failed resulting in a negative SNR. Epochs with artifacts would result in falsely high amplitude of the signal of interest, and yet a low accuracy. We therefore expect automated artifact rejection methods to improve phase accuracy.

Outlook

Whereas our focus was the calibration of causal phase estimators using real EEG data, a further in-depth analysis of realistic synthetic EEG data with known phase but different confounders would clearly be merited. It would also be interesting to analyze a simple simulated case, where the Cramér–Rao Bound can be analytically derived from sinusoids in 1/f noise, following a similar approach to the derivations considering Gaussian noise (Peleg and Porat, 1991; Sameni and Seraj, 2017), which would provide a theoretical upper limit to the performance of any estimator. On the other hand, realistic head modeling would also enable scenarios to be studied where multiple sources oscillating at similar frequencies (e.g. occipital and somatosensory 8–13 ​Hz oscillations) contribute to the recorded signal, enabling different spatial filter configurations to be compared with respect to their ability to differentiate the two sources.

Finally, a further expected benefit of the PHASTIMATE toolbox and the large sensorimotor rhythm dataset is to facilitate the development of improved real-time phase estimation algorithms by providing a common benchmark. A similar class of algorithms using the Fourier transform (Mansouri et al., 2017; Zrenner et al., 2015) or wavelets (Madsen et al., 2018) that are mathematically analogous to the method implemented here might be more efficient. A distinct class of approaches that makes use of additional temporal information would use prior expectations (e.g., peak frequency, shape of the oscillation) or the information from previous time steps (e.g., by implementing a Kalman filter) to improve upon the algorithm tested here. Finally, the spatial dimension of the oscillation is also expected to confer additional information based on phase differences (Alexander et al., 2006) which should increase the performance of a local estimate, but a spatial approach would be considered out of scope of the current work and indeed, if the spatial dimension is a relevant measure in a given experiment, this would obviate the need to determine spatially localized phase.

Data availability

All data used in this report is available for download: https://gin.g-node.org/bnplab.

The PHASTIMATE toolbox and the analysis scripts used to generate the figures in this report are available for download: https://github.com/bnplab/phastimate.

Disclosures

CZ is coordinator of and partially funded through an EXIST Transfer of Research grant by the German (grant 03EFJBW169s). The goal of this grant is the commercialization of a real-time EEG analysis device through a spin-off start-up to enable therapeutic brain-oscillation synchronized stimulation.

CRediT authorship contribution statement

Christoph Zrenner: Conceptualization, Methodology, Investigation, Formal analysis, Writing - original draft. Dragana Galevska: Investigation, Project administration. Jaakko O. Nieminen: Investigation, Writing - review & editing, Funding acquisition. David Baur: Investigation. Maria-Ioanna Stefanou: Investigation. Ulf Ziemann: Conceptualization, Funding acquisition, Writing - review & editing.