Prediction of cortical theta oscillations in humans for phase-locked visual stimulation

BACKGROUND
The timing of an event within an oscillatory phase is considered to be one of the key strategies used by the brain to code and process neural information. Whereas existing methods of studying this phenomenon are chiefly based on retrospective analysis of electroencephalography (EEG) data, we now present a method to study it prospectively. New Method: We present a system that allows for the delivery of visual stimuli at a specific phase of the cortical theta oscillation by fitting a sine to raw surface EEG data to estimate and predict the phase. One noteworthy feature of the method is that it can minimize potentially confounding effects of previous trials by using only a short sequence of past data.


RESULTS
In a trial with 10 human participants we achieved a significant phase locking with an inter-trial phase coherence of 0.39. We demonstrated successful phase locking on synthetic signals with a signal-to-noise ratio of less than -20 dB. Comparison with Existing Method(s): We compared the new method to an autoregressive method published in the literature and found the new method was superior in mean phase offset, circular standard deviation, and prediction latency.


CONCLUSIONS
By fitting sine waves to raw EEG traces, we locked visual stimuli to arbitrary phases within the theta oscillatory cycle of healthy humans.


Introduction
As early as 1929 rhythmic variations of the electric potential were recorded from the surface of the human brain in the absence of external stimulation (Berger, 1929). Such spontaneous oscillations are ubiquitous in brains of humans and animals (Basar and Güntekin, 2008), can be observed both in sleep and awake states, and are influenced by exand intrinsic factors (Karakas and Barry, 2017;Dang-Vu, 2012). Some researchers have argued that oscillations are fundamentally involved in neural processing of sensory, motor, and cognitive information (Buzsáki, 2006;Engel and Fries, 2001).
Oscillations are typically measured from electroencephalogram (EEG) data. Note that even though there is no uniform definition of the EEG frequency bands (Newson and Thiagarajan, 2019), and functional correlates differ in frequency depending on species (Jacobs, 2014), in this paper we refer to any oscillation in the range of 4-8 Hz as theta and 8-12 Hz as alpha. Many associations between oscillations and behavioral patterns have been described. For example, the power (i.e., amplitude squared) and frequency in the alpha and theta range have been linked to attention and perception (Ishii et al., 1999;Thut et al., 2006;Kastner et al., 1999). Oscillations in the theta range have also been linked to memory load, formation, and retrieval (Jensen and Lisman, 2005;Klimesch, 1999). A number of studies have shown that the perception of a visual stimulus varies with not only oscillatory frequency and power but also the actual phase within an oscillatory cycle ( Van-Rullen, 2016;Busch et al., 2009;Klimesch et al., 2007;Drewes and VanRullen, 2011;Mathewson et al., 2009;Busch and VanRullen, 2010;Stefanics et al., 2010;Schyns et al., 2011;Helfrich et al., 2018;Klimesch, 2012;Nokia et al., 2015;Fiebelkorn et al., 2013;Samaha and Postle, 2015).
Phase-specific effects reported in the published literature were mostly studied using post-hoc analysis of EEG signals that were recorded during perceptual tasks Harris et al., 2018;Hanslmayr et al., 2013;Fiebelkorn et al., 2013;Milton and Pleydell-Pearce, 2016;Mathewson et al., 2009;Samaha and Postle, 2015). Some researchers used entrained oscillations, i.e. oscillations driven by rhythmic external stimuli, to show phase-specific behavioral effects (Nakayama and Motoyoshi, 2019;Stefanics et al., 2010;Romei et al., 2012;Ronconi and Melcher, 2017). We now describe a method to systematically study phase-specific effects in healthy human subjects by using surface EEG recordings to predict stimulus timing.

Materials and methods
In this section, we first introduce our phase-estimation procedure named FitSine and describe the hardware setup used. We then explain the simulation and evaluation procedures. These simulations are designed to show the properties of our phase estimation and prediction algorithm, and compare them to existing methods. Finally, we introduce the experiment we used to test our method in a trial with healthy human participants.

Phase estimation
The goal of our FitSine algorithm is to estimate the instantaneous frequency and phase of a surface EEG signal in order to deliver a stimulus at that point in time where the oscillation is at a predetermined phase. We will refer to this concept as phase-locked stimulation. The algorithm is based on fitting a discrete set of 41 sinusoids with equally spaced frequencies in the theta range (4-8 Hz) to the raw input signal, and then selecting the best fitting sinusoid. We then extrapolate the best fitting sinusoid to predict the oscillatory phase at the earliest possible stimulus onset.
This approach allows to work with less than one oscillatory cycle of past data and compensate for latencies by prediction. Here we use the term latency to describe the time period between the acquisition of the signal and the presentation of the stimulus on the screen.
We use a general sinusoid in the form y(t) = Asin(2πft + φ) + C, and rewrite it for discrete times to consist of two sinusoids that are offset in phase by one-quarter cycle i.e. an in-phase and a quadrature component.
We define our base function (1) at a discrete time point t n using scalar n ∈ {1, ..., N}, N ∈ N as time index, and use scalar j ∈ {1, ..., J}, J ∈ N as index for a discrete frequency out of the vector of test frequencies f ∈ R J . We denote the amplitudes belonging to frequency f j as A j and B j , where A ∈ R J belongs to the cosine part, B ∈ R J belongs to the sine part, and the constant offset is C ∈ R J .
We define the cost for fitting the sinusoid to the measured data x = [x 1 , ..., x N ] ∈ R N as the sum of the squared error: (2).
Evaluation of the optimal fit for the given data requires minimization (5) of the cost function (2), which can be done by differentiation and subsequent solving for zero, leading to (6).
The phase φ can now be estimated using (7).
As can be seen in (6) the inversion and the trigonometric functions can be precomputed, leaving only two matrix multiplications to evaluate the fit for a single sinusoid. Thereby computation time during runtime can be reduced.

Benchmark simulations
To compare the performance of our FitSine algorithm with existing methods, we performed a benchmark simulation using Matlab 2018b and synthetic data.
We used an autoregressive (AR) model to predict the data similar to the one described by Chen et al. (2013) using intracranial EEG, and implemented it as described by Blackwood et al. (2018).
In short, the AR based method is based on first filtering the signal using a second order Butterworth IIR filter. Afterwards, a 20th-order AR model (AR20) is trained individually for each data sequence using Burg's lattice-based method, which is then used to predict future data points. Hilbert transformation is then applied to the predicted data points to estimate the phase.
For the benchmark simulation we use a linear chirp that changes from 4 Hz to 8 Hz over the course of 1000 s, and added various levels of white noise to reach the SNR levels specified below. For each algorithm we simulated 1000 predictions at randomly sampled locations within the chirp signal for each combination of SNR level (− 50 dB, − 40 dB, − 30 dB, − 25 dB − 20 dB, − 10 dB, 0 dB, 10 dB), input buffer length (100 ms and 300 ms), and prediction time (10 ms, 25 ms, and 34 ms).

Implementation
To test our FitSine algorithm for phase-specific stimulation we have implemented a setup consisting of a BRAIN VISION actiCHamp EEG amplifier equipped with one 32-channel acquisition module, a PC, and an LCD monitor. We used 3 photodiodes connected to the auxiliary (AUX) input of the EEG amplifier to acquire timing information, one mounted to the top right corner of the monitor to detect the true frame switch timing, and the other two fixed to LEDs connected to the local parallel port (LPT) of the PC. An overview of the system is shown in Fig. 1.
The visual stimuli are presented on an Asus ROG Swift PG279Q monitor running at 144 Hz with a custom color calibration profile. We used a desktop computer system running on Windows 10 Pro, equipped with an Intel I7-8700k running at 4.7 GHz, 32 GB RAM, and an NVIDIA GTX 1080Ti as the platform for all measurements.
We use the functionality of the Lab Streaming Layer library (liblsl) project (LabStreamingLayer, 2019) and the pre-compiled binaries to create a LabStreaminglayer (LSL) stream to transfer the data from the EEG amplifier into Matlab. For that we use the "ActiChamp Connector" provided by LSL community on the same PC to configure our EEG amplifier and stream out the data at a native sampling rate of 10 kHz.
The phase estimation is calculated every 2 ms, using our Matlab implementation of the algorithm described in subsection 2.1. For the phase estimation the most recent 100 ms (N = 1001) of past data is used, without any pre-processing.

Hardware in the loop simulation
We tested the system performance and delay compensation methods using a hardware-in-the-loop-simulation. For that we used a Matlab script, running in a separate instance of Matlab on the same PC, to acquire the LSL stream every 2 ms. We then replace the samples from the unconnected EEG electrodes with the synthetic test-signals, and stream it forward to the main application using a second LSL stream. Exchanging only the EEG channels with synthetic data, and keeping the timing information in the AUX channels allows testing the system performance and timing compensation methods without the need for additional hardware.
We used a 6 Hz sinus with different levels of added white noise to test the online implementation of the FitSine algorithm, and used the same simulation to compare the accuracy of the offline phase-extraction method (see subsection 2.5) against a ground truth signal.

Analysis
To evaluate the performance of the two phase-estimation algorithms we compared the predicted phase with an off-line phase estimate obtained by 4-8 Hz bandpass filtering and Hilbert transformation. The finite impulse response (FIR) bandpass filter was designed in Matlab using separate high-and low-pass filters with a window length 3 times larger than the inverse of the respective frequency. Offline zero-phase filtering was performed on the GPU using a forward backward approach.
We use the inter-trial-coherence (ITC) value, as defined in formula (8), to measure the clustering of stimulus onset times towards a specific oscillatory phase. An ITC value of 1 stands for a "perfect" phase lock (a monophasic distribution wherein each stimulus is delivered at that point in time where the oscillation is at the predetermined phase), and 0 stands for uniform data (no clustering towards any angle).
Where K stands for the total number of trials and φ for the phase.

Online phase-specific stimulation experiment
In order to evaluate the performance of the system we conducted an in-vivo experiment studying the influence of the stimulus timing within the theta oscillatory cycle (4-8 Hz) on contrast-sensitivity. The task of the experiment was to report the orientation of a Gabor patch with 4 possible orientations.
For each subject we arbitrarily selected two 180 ∘ opposite phase angles as target for the stimulation onset. We use the term target angle to refer to the intended stimulation phases and refer to these two angles as The stimulus consists of a Gabor patch with stripes in one out of 4 possible orientations (0 ∘ , 90 ∘ , − 45 ∘ , 45 ∘ ), using a spatial resolution of 9 pixels ( ≈ 14.45 arcmin) for the sinusoidal base, and a width of the modulating Gaussian of 200 pixels ( ≈ 320.45 arcmin). The stimulus is presented in the center of the monitor and its contrast adjusted using QUEST (Watson and Denis, 1983) to a correct detection-rate of ≈ 60%. Psychtoolbox (PTB) (Kleiner et al., 2007) is used for stimulus generation and presentation.
At the beginning of each trial, a neutral gray screen is shown for a duration consisting of a uniformly distributed random time-period between 200 ms and 500 ms, and the time it takes until the requirement for phase-locked stimulation is fulfilled. For this requirement, the phase at the majority of the three occipital electrodes, as determined with FitSine and extrapolated to the next possible stimulus time-point, has to lie within ± 10 ∘ of the target angle. As soon as this is fulfilled a Gaborpatch is shown for ≈ 48.6 ms, followed by a neutral gray screen. The participants then have to report the orientation of the Gabor-patch without time limit using a forced choice test with four alternatives. The response of the participant on the number block of the keyboard starts the next trial. The trial sequence used for the experiment is shown in Fig. 2.
Additionally we performed a second experiment, denoted Random, with a similar task without phase locking in each participant. For the Random experiment we adapted the duration of the initial gray screen to match the inter-trial intervals of the phase-specific experiment. For this, we prolonged the initial random duration to a value between 200 ms and 700 ms.
10 healthy volunteers (3 females, ages 24-47 years, mean 31, SD 7) participated in the study. They all provided written informed consent. All procedures were approved by the local ethics committee (Kantonale Ethikkomission Bern (KEK), Basec PB_2016-00250), in accordance with the principles of the Declaration of Helsinki (World Medical Association Declaration of Helsinki, 2013).
Experiments were carried out in a dimly lit room, with participants Fig. 1. Block diagram of the setup used for phase-selective stimulation. EEG and trigger-signals are recorded using a BRAIN VISION actiCHamp amplifier using 32 + 8 channels, and transferred to the PC via USB2 using the Lab-StreamingLayer (LSL) app for device control. Data transfer into Matlab is achieved using a LSL stream for portability. Matlab is used for all data processing and experiment control. Visual stimuli are generated in Matlab using Psychtoolbox. Control-signals on Parallel-port (LPT) are used for delay measurements.

Fig. 2.
Illustration of the task. Trial starts with a button-press on the keyboard. After a variable delay (200-500 ms uniformly distributed plus the time until the specific oscillatory phase is met) the stimulus (Gabor-patch with 4 possible orientations) is shown for 7 consecutive frames (48.6 ms). This is followed by a neutral gray background. Infinite time is given to the participant to assign one of four possible orientations to the trial using the number block of the keyboard. The keyboard input is used to trigger the start of the next trial.
seated ≈ 50 cm away from the monitor. EEG was recorded from Cz, C1, and C2 (using a 32-channel EASYCAP EEG cap arranged according to the extended international 10-20 system (Nuwer et al., 1998) using Fz as reference and FPZ as ground electrode. Impedance during the participant preparation was monitored and SuperVisc High Viscosity Electrolyte-Gel for active electrodes (EASYCAP, Herrsching, Germany) applied to achieve a targeted impedance below 10 kΩ.

Simulations
We compared our FitSine approach as introduced in section 2.1 with the AR based method. Computation-time for both algorithms was evaluated using 10 ′ 000 consecutive runs with an input sequence consisting of 1001 random samples. If not otherwise noted, we report the results as mean ± standard deviation (mean ± SD). We measured a computation time of (133.22 ± 474) μs for the FitSine algorithm with 41 equally spaced frequency bins between 4 and 8 Hz. We measured the following values for the single steps required by the AR-20 method for a 10 ms prediction window: Second order IIR bandpass filtering (68.80 ± 779) μ s; Training time for the AR-20 model using Burg's lattice-based method (5.56 ± 48) ms; Iterative prediction (791.27 ± 4108) μs; Phase estimation using the Hilbert filter (156.72 ± 1623) μs. This resulted in a total computation time of 6.57 ms for the AR-20 method.
Next, we compared the phase-locking performance of our FitSine approach to the AR-20 based method. According to the steps described under Benchmark simulation in section 2.2. The results of this simulation are visualized in Fig. 3, which illustrates that FitSine shows a better phase-locking performance with a higher ITC and a smaller mean deviation of the target, which we termed "phase offset". Moreover, we found a smaller circular standard deviation (circ std) for the FitSine as compared to the AR-20 based prediction.
To test the implementation of our FitSine algorithm, and the whole setup as shown in Fig. 1, we performed a hardware-in-the-loop simulation. The resulting ITC-values (summarized in Table 1), indicate a reliable prediction for SNR values down to − 20 dB. The ITC values shown under "ITC post", are computed using the offline phase-detection method as reference, as used later on with the EEG measurements. The values under "ITC Reference" are computed with the input signal of the simulation as reference. We found a mean offset of the target phase of 6.30 ∘ with a circular standard deviation of 19.83 ∘ for the 0 dB SNR. The signal replacement step increased the input-latency to a total of (8.32 ± 061) ms.

Experiment
We measured an input-latency during the experiment of (6.34 ± 63) ms between the issuing of a signal to the local parallel-port and its detection in Matlab.
For the phase-locking conditions we achieved a circular mean phase offset of 3.50 ∘ , and a circular standard deviation of 66.55 ∘ , over all trials.
To illustrate the phase-locking performance of the whole system we combined all phase-selective trials into one plot (Fig. 4) by subtracting the target phase. This figure incorporates both the offsets in the mean angle between the single participants, and the variation around it. Therefore, the ITC of all trials combined, as shown in Fig. 4, is smaller with 0.373 than the average of the ITCs from each participant (Table 2) with 0.394.
We checked the efficacy of the phase-locking system by testing the resulting phase distribution as detailed above for deviation from the uniform distribution using Rao's test for circular uniformity (R version 3.4.3, and the CircStats package version 0.2.6). Significant deviations from the uniform distribution were found for both phase-specific conditions, and no significant deviation for the random case [A: p = < 0.001, B: p = < 0.001, Random: p = > 0.10].
Then we computed ITC values to measure the precision of the system (summarized in Table 2). One-way ANOVA was conducted to compare the phase locking for the 3 conditions. Normality checks and Levene's test were carried out and the assumptions met using a significance level of 0.05. There is a significant difference in the mean ITC [F(2, 27) Fig. 3. Benchmark simulation using a 4-8 Hz linear 1000 s chirp signal to compare our Fit-Sine method against the AR based algorithm using input buffer lengths of 100 ms (solid lines), and 300 ms (dashed lines). Blue traces correspond to FitSine, green to AR based. Prediction windows of 10 ms, 25 ms, and 34 ms are shown using markers (none, circle, cross). Top: Inter-Trial-Coherence (ITC), Middle: Phase offset in degrees between the target angle and the mean angle at the predicted time point in the future, Bottom: Circular standard deviation of the phase at the predicted time point.

Table 1
Inter Trial Coherence (ITC) values for target angles at 0 ∘ (A) and 180 ∘ (B) using simulated oscillations with a range of signal to noise (SNR) ratios. Oscillations were simulated using 100 ms sequences of a 10 kHz sampled 6 Hz sine wave with added white Gaussian noise. = 40.37, p = 7.69 × 10 − 9 ] between the conditions. Post hoc comparisons using Tukey's test shows a significant difference between the phaselocked conditions (A and B) and random: A vs. random (p = 39.03 × 10 − 9 ) and B vs. random (p = 128.36 × 10 − 9 ). There is no significant difference in the ITC between the two phase-locked conditions (p = 0.87). The phase locking for all participants is shown in Fig. 5 using polar histograms. We used Bold lines to indicate the target angles, the used target angles and the resulting ITC values are listed in Table 2.

Discussion
By using a sine wave fitting procedure to raw EEG signals, we were able to lock visual stimulus presentation to any phase angle within the theta cycle. Our system makes no constraints on the target phase such as approaches based on peak or level detectors. This allows prospective testing of hypotheses derived from correlative evidence concerning arbitrary theta phases.
This was achieved using efficient computing combined with commercially available equipment consisting of an EEG-amplifier, a monitor, and a computer. Such a closed-loop visual stimulation system in the theta band using surface electrodes in humans is to our knowledge unprecedented. A great number of post-hoc analysis of stimulus timing in rodents, monkeys, and humans demonstrated that the location of a stimulus within an oscillatory cycle critically determines synaptic plasticity (Law and Stan Leung, 2018;Teyler et al., 2005;Hyman et al., 2003;Christian et al., 1997), memory formation and memory retrieval (Klimesch, 2012;Jensen and Mazaheri, 2010;Kaiser and Lutzenberger, 2005;Osipova et al., 2008;Lisman and Idiart, 1995;Tseng et al., 2018;Reinhart and Nguyen, 2019;Solomon et al., 2017).
Our phase prediction method allows for phase-locked visual stimulation in the theta range. The phase difference between the target and stimulus angle is illustrated in Fig. 4 where a clear clustering of the stimuli around 0 • is visible. The clear clustering around the target angle is also visible in Fig. 5 for the single participants. The averaged ITC value for the phase-locked condition in each subject is bigger than the biggest ITC value for random stimulation, indicating a successful phase locking on the subject level.
We consider the mean phase offset of 3.5 • as small compared to the much larger variability of the phase offset in individual trials (standard deviation = 70 • ). This relatively high accuracy is possible since we do not induce systematic errors with non-zero-phase filters. The frequency dependent phase shift of such filters cannot be compensated as the frequency varies between participants and also during the recording. On a trial by trial basis the target offsets were variable with a standard deviation of 70 • . In Fig. 3 we show that the low SNR is the primary cause of this variability.
The mean angle close to 0 of the AR method should only be interpreted in combination with the circular standard deviation, as the averaging over the frequency sweep leads to a cancellation of the systematic offsets between 45 • at 4 Hz and − 45 • at 8 Hz introduced by the second order IIR filter. The ITC value is a pure precision metric and does not report accuracy metrics for which we use the mean phase offset.
We limited the usage of past data to a short sequence of 100 ms. The goal of this is to prevent event related potentials and artefacts from previous trials from influencing the phase prediction when using short inter-trial intervals. In contrast to the iterative prediction procedure of the AR method the computation time for the phase estimation with FitSine is independent of the prediction window. Comparing the total time of all required computation steps for 10 ms prediction, our approach of using sine fitting to raw EEG signal in combination with precomputation resulted in a 50 times smaller computation time than the AR based method (133 μs vs. 6.57 ms).
We arbitrarily selected a frequency resolution of 0.1 Hz, leading to 41 discrete frequencies for the evaluation of the fit. The quantification error we introduce with this 0.1 Hz step size leads to an angular error of maximally 0.18 • , which is insignificant in comparison to other sources of errors in the system. We refrained from using less frequencies as the gain in computational time would be negligible compared to the input delay which is in the millisecond range. As the FitSine algorithms only requires matrix multiplications, it can be easily and efficiently implemented on a wide range of hardware and software platforms.
More sophisticated models such as the AR-20 method may model the dynamics of a signal with more complexity than our FitSine procedure and could thus achieve a better prediction. However, in our case the AR-20 cannot be applied on unfiltered raw signals. In the design of the filter required to pre-process the signal for the AR method compromises between filter-steepness, group-delay, and non-linear phase-transfer properties have to be weighed against each other for online processing (Widmann et al., 2015). For example, the second-order IIR filter used in the Benchmark simulation introduces a phase shift of ∓ 45 • between 4 and 8 Hz. Digital filters with less degrading properties are possible, but are of higher order and their associated group delay largely outweighs the better prediction performance of AR models, as the introduced delay has to be compensated for by increasing the prediction interval.
The fast computation in combination with the measured input delay of our setup of less than 7 ms allowed us to use a basic representation of  the signal for the required prediction window. With the AR model we were not able to achieve a whole system delay below the 10 ms used in the simulation due to a significantly higher computation time. In Fig. 3 we showed the influence of different prediction times on the phaselocking performance.
The 10 ms prediction with the AR method would therefore not be possible using our setup. Due to the closed loop, our setup allows to adapt to variable latencies in the system, as shown in the hardware-inthe-loop simulation.
A general limitation in the analysis of EEG-signals is the lack of a ground truth signal. The theta band we used in this experiment has a low SNR, and the values we used for evaluation might be corrupted by noise. Tests with simulated signals showed that our used post-hoc comparison signal is valid at least up to a SNR of − 20 dB. Table 1 shows that we were not able to reliably reconstruct the signals below a SNR of − 20 dB, neither with the FitSine nor the FIR-filter. Better performance in the synthetic test case could be achieved by choosing longer sequences of past data, which leads to an improvement of the SNR due to the uncorrelated nature of the noise.
EEG recordings are noisy and analysis usually involves different filtering and signal processing steps. Filtering introduces artefacts which vary with window size, and filter-type (Vanrullen, 2011;Widmann and Schröger, 2012). Such artefacts are particularly problematic if the same filters are used in the online prediction and in the post-hoc analysis to determine the performance of the algorithm. Given that we have used sine wave fitting and extrapolation on the raw signal for the prediction, and then used a classical FIR filter in the post-hoc analysis procedure, we believe that artefacts and bias we may have introduced in the prediction should not propagate through the post-hoc analysis part.
Our FitSine approach can also be applied to frequency bands other than theta. In an unpublished and unreviewed pilot study with 9 participants we used the same protocol in the alpha range, and measured an average ITC of (0.57 ± 0.16). We attribute the significantly better phase-locking performance to the higher SNR in the alpha band for surface EEG signals (see Fig. 3 for the dependency between SNR and the ITC).
An analysis of the impact of a phase-specific visual stimulus on the corresponding psychophysical response will be the scope of future publication.
In conclusion we show a relatively simple way of providing visual stimuli at a specific oscillatory phase of a cortical theta cycle.

Data availability
Code will be provided upon request to interested parties. EEG data will not be shared at this stage, as further studies on this dataset are ongoing.