Abstract
Brain rhythms have been proposed to facilitate brain function, with an especially important role attributed to the phase of lowfrequency rhythms. Understanding the role of phase in neural function requires interventions that perturb neural activity at a target phase, necessitating estimation of phase in realtime. Current methods for realtime phase estimation rely on bandpass filtering, which assumes narrowband signals and couples the signal and noise in the phase estimate, adding noise to the phase and impairing detections of relationships between phase and behavior. To address this, we propose a state space phase estimator for realtime tracking of phase. By tracking the analytic signal as a latent state, this framework avoids the requirement of bandpass filtering, separately models the signal and the noise, accounts for rhythmic confounds, and provides credible intervals for the phase estimate. We demonstrate in simulations that the state space phase estimator outperforms current stateoftheart realtime methods in the contexts of common confounds such as broadband rhythms, phase resets, and cooccurring rhythms. Finally, we show applications of this approach to in vivo data. The method is available as a readytouse plugin for the Open Ephys acquisition system, making it widely available for use in experiments.
Introduction
Rhythms, consistently periodic voltage fluctuations, are an ubiquitous phenomenon observed in brain electrophysiology across scales and species. Many studies have described a relationship between rhythms and behavior (Buzsáki, 2006; Haegens and Zion Golumbic, 2018; VanRullen, 2016). A prominent feature of rhythms proposed as relevant for neural processing is the phase. Phase has been proposed to coordinate neural spiking locally and potentially, through coherent networks, even globally (Fries, 2015; Maris et al., 2016). For example, the phase of slow rhythms (3–15 Hz) has shown relationships to dynamics in perception (Busch et al., 2009; Gaillard et al., 2020; Gregoriou et al., 2009; Helfrich et al., 2018). Phaseamplitude synchrony is common (Canolty and Knight, 2010; Hyafil et al., 2015), and has been suggested to change in neurodegenerative disease (Hemptinne et al., 2013). Crossregional, lowfrequency phase synchrony is consistently observed as a correlate of topdown executive control (Widge et al., 2019). Most work examining the importance of phase for neural dynamics and behavior has been correlative in nature. To better understand the functional relevance of the phase of a rhythm requires an ability to monitor and perturb phase relationships in realtime.
Several methods exist to estimate phase in realtime (Blackwood et al., 2018; Chen et al., 2013; Mansouri et al., 2017; Rodriguez Rivero and Ditterich, 2021; Rutishauser et al., 2013; Siegle and Wilson, 2014; Zrenner et al., 2020). At their core, these methods rely on bandpass filtering, a mechanism to forward predict data (using autoregressive [AR] models or linear extrapolation of phase) and the Hilbert transform to estimate the phase. Brain rhythms are often nonsinusoidal (Cole and Voytek, 2017) and broadband (Buzsaki, 2004; Roopun et al., 2008), making development of an accurate bandpass filter difficult. While the contemporary modeling approach using filters accurately estimates the phase across a range of contexts, several limitations exist that limit the phase estimate accuracy (Matsuda and Komaki, 2017; Siegle and Wilson, 2014): (1) By depending on bandpass filters, existing realtime phase estimators are susceptible to nonsinusoidal distortions of the waveform, and inappropriate filter choices may miss the center frequency or total bandwidth of the rhythm. (2) Phase resets, moments when the phase slips because of stimulus presentation or spontaneous dynamics, cannot be tracked using filters, which control the maximum possible instantaneous frequency. (3) By filtering the observed data directly, these approaches do not model signal and noise processes separately. Phase estimates thus represent signal and noise phase together. (4) Many existing approaches depend on buffered processing (e.g., Fourier or related transforms), meaning that current phase estimates are delayed in ways that may not support realtime intervention. Further, bufferbased processing is susceptible to edge effects that are not of concern in traditional, offline analysis. (5) Contemporary approaches do not define a measure of confidence in the phase estimate. Despite these limitations, realtime phase simulations using current techniques have produced interesting results (Cagnan et al., 2017; Desideri et al., 2019; Hyman et al., 2003; Kundu et al., 2014; Schaworonkow et al., 2018; Siegle and Wilson, 2014; Zrenner et al., 2018), motivating the value of generating better methods to estimate phase causally.
To address the limitations of existing approaches, we implement here a dynamical systems approach to realtime phase estimation. Rather than filtering the data to identify the rhythm of interest, we focus on tracking the analytic signal (a complex valued entity, representing the amplitude and phase of a rhythm Kramer and Eden, 2016) as a latent process in the observed time series using a model proposed by Matsuda and Komaki, 2017. In doing so, we remove the requirement of bandpass filtering and instead model the signal as a state space harmonic oscillator (or combination of such oscillators) driven by noise. This modeling approach allows: (1) a separate estimate of observation noise, (2) a procedure to model and reduce the confounding activity of other rhythms, and (3) a principled method for assessing uncertainty in the realtime phase estimates. After first fitting parameters of the state space harmonic oscillator model acausally, the method then tracks the phase in realtime assuming spectral stationarity. We call this approach the state space phase estimate (SSPE).
To benchmark the SSPE technique against the current stateoftheart methods of realtime phase estimation, we design several simulations and analyze case studies of in vivo data. We limit our comparison to the realtime methods proposed in Zrenner et al., 2020, and Blackwood et al., 2018, methods which provide an estimate of phase on receipt of each new sample of data and have readily available reference code. We additionally apply a standard method for nonrealtime phase estimation (an acausal finite impulse response [FIR] filter followed by a Hilbert transform; Lepage et al., 2013). We note that these and other existing realtime phase estimation techniques (Chen et al., 2013; Mansouri et al., 2017; Rodriguez Rivero and Ditterich, 2021; Rutishauser et al., 2013; Siegle and Wilson, 2014) apply a bandpass filter and forecasting technique, and, by that measure, the acausal FIR approach by utilizing more information suggests the room for improvement available for these estimators. We test the ability of these different methods to track the phase of rhythms with different spectral profiles – including narrowband rhythms, broadband rhythms, multiple rhythms, and rhythms with phase resets, and in different signaltonoise settings. We show that the SSPE method outperforms the existing realtime phase estimation methods in almost all cases. Finally, we illustrate the application of the SSPE method to in vivo local field potential (LFP) and electroencephalogram (EEG) recordings, demonstrate the utility of SSPE for predicting behavior, and introduce a plugin of the SSPE method for Open Ephys (Siegle et al., 2017), making the method widely available for use in realtime experiments. Our work demonstrates that the SSPE method improves the ability to track phase accurately, in realtime, across a diverse set of contexts encountered in data.
Materials and methods
State space model framework
Request a detailed protocolTo estimate phase in realtime, we utilize a datadriven model (Matsuda and Komaki, 2017) that operates based on principles from dynamical systems and Markov models. This type of model (called a state space model) separates what is observed (called the observation equation) from what we wish to estimate (called the state equation). The state equation attempts to capture the underlying, unobserved dynamics of the system, while the observation equation transforms the state into the observed signal. We generate an optimal prediction for the state using a Kalman filter that compares the initial prediction of the observation with the actual observed value.
Each state variable represents an oscillator with a fixed frequency that undergoes small random perturbations. As these perturbations accumulate, the instantaneous phase and frequency of the oscillator become more uncertain, since the oscillator is not directly observed. We assume that the observed data are noisy versions of these oscillators. The state space model framework combines information from two models: (1) a state model that describes the oscillator dynamics and (2) an observation model that describes how the oscillator generates the data. The methods we develop use both of these sources of information to optimally estimate the instantaneous frequency and phase of each oscillator at every moment in time.
We implement a linear state equation that models the real and imaginary part of the complex valued analytic signal of rhythms present in the neural data. For each rhythm modeled, the state equation dynamics rotate the analytic signal at that rhythm’s frequency. The observation equation sums the real part of the analytic signals of every modeled rhythm to predict the observed neural activity. We treat as realvalued the real and imaginary parts of the complex valued analytic signal for ease of analysis and interpretation. This state space model structure allows us to derive the amplitude and phase of each modeled rhythm, and estimate the analytic signal in the temporal domain, thus avoiding complications associated with windowing (i.e., the need to choose large enough window sizes to support sufficient frequency resolution and the associated unavoidable temporal delay in phase estimates) that occur in related frequency domain approaches (Kim et al., 2018).
The state space model is linear in the state evolution (dynamics of state transitions) and observation equations. Further, the covariance for the state (which determines the bandwidth for the rhythms) and noise for the observation is Gaussian (white noise). These assumptions (linearity and Gaussianity) imply that we can use the Kalman filter to correct the state estimate based on the true observation (i.e., to filter the state) (Kalman, 1960). This type of modeling approach has been successfully applied in many contexts to track latent structure from observed data while accounting for state and observation noise effects (Brockwell et al., 2004; Eden et al., 2018; Galka et al., 2004; Yousefi et al., 2019).
To illustrate this approach, we first consider the case of tracking a single, 6 Hz rhythm. The state equation for this scenario is:
where $x}_{t$ has both a real and imaginary part. At time instant t−1, the modeled rhythm is in state $x}_{t1$, that is, the rhythm is at some point along its trajectory. The state equation provides a prediction for the next point along this trajectory ($x}_{t$) by rotating the current state ($x}_{t1$) and adding noise ($u}_{t$). The rotation matrix $O(\omega )$ rotates the real and imaginary parts of the state by the amount $\omega$ appropriate for a 6 Hz rhythm:
where Fs is the sampling frequency in Hz. After rotation, the state is scaled (or rather, damped) by a factor $a$, and then Gaussian noise (with mean and variance $Q$) is added. We estimate the observation at time t ($y}_{t$) as the real part of the state at time t with additive noise:
where $\nu$ is Gaussian noise with mean 0 and variance $\sigma}^{2$. Note that this model for the state is akin to a damped harmonic oscillator driven by noise, a model that has been shown to be relevant and useful for gamma rhythms in the visual cortex (Burns et al., 2010; Spyropoulos et al., 2019) and also for EEG (Franaszczuk and Blinowska, 1985).
To model N rhythms, we define X as the state variable, a 2N × T matrix over total time T. An individual time point with state ${X}_{t}$ (2N × 1 vector) is rotated according to each rhythm’s (${x}_{j}^{t}$ , a 2 × 1 complex vector for rhythm j) central frequency (${\omega}_{j}=2\pi {f}_{j}\u2206t\in \left(0,\pi \right)$ , where ${f}_{j}$ is frequency in Hz and $\u2206t$ is sampling interval) and damped by a constant (${a}_{j}^{}$). For each rhythm ${x}_{j}^{t}$ we define ${\theta}_{j}^{t}$ as a scalar quantity tracking the phase, and Θ as an N × T matrix of phase values for all rhythms over all time. For the N rhythms, O is the rotation matrix, a 2N × 2N and block diagonal matrix, defined by the central frequency of each rhythm. A is the scaling matrix, a diagonal 2N × 2N matrix, with scaling factor ${a}_{j}\in \left(\mathrm{0,1}\right)$ to ensure stability, the same for each rhythm, and different across rhythms. The driving noise for the state U is Gaussian, a 2N × T matrix with covariance structure Q, a 2N × 2N diagonal matrix with the same variance for the real and imaginary parts of each rhythm. The state is collapsed into the observation $Y$(a T × 1 vector) through M (a 1 × 2N vector) which summates the real parts of the analytic signals of all rhythms. For any time point, the observation ${y}_{t}$ is a scalar. The noise for the observation ${v}_{t}$ is Gaussian with covariance. $\sigma}_{R}^{2$ The model equations are:
We note that the state (Q) and observation ($\sigma}_{R}^{2$) covariance matrices are estimated from the data, as described below. While we term ${{x}_{j}}^{t}$ as tracking the analytic signal, this nomenclature differs from the signal processing literature. In signal processing, the analytic signal is defined as a signal that has zero power at negative frequencies. However, given the stochastic frequency modulation permitted in the state vector ${{x}_{j}}^{t}$ , we cannot guarantee that ${{x}_{j}}^{t}$ has zero power at negative frequencies. Nevertheless, since for sufficiently complex signals there are (potentially) several ways one could represent the dynamics of the amplitude and phase (Boashash, 1991; Rosenblum et al., 1997), we use the term analytic signal here for ease of understanding our intent with this model and comparing it to classical Hilbertbased methods.
Realtime phase estimation
Request a detailed protocolWe perform realtime estimation of phase as follows. First, we use an existing interval of data to acausally fit the parameters of the state space model ${a}_{j}$ , ${\omega}_{j}$ , ${Q}_{j}$ , $\sigma}_{R}^{2$ using an expectationmaximization (EM) algorithm as proposed by Soulat et al., 2019; for details, see the supplementary information in Soulat et al., 2019. Under the EM approach, optimization of parameters follows a twostep process. First, initial values for parameters are selected, usually from prior knowledge such as an examination of the power spectrum of an initial data sample. This allows the algorithm to estimate the expectations for the state and the observation. Second, using the state and observation estimates, an analytic solution exists (Shumway and Stoffer, 1982) for the parameters that maximizes the likelihood. We repeat the expectation and maximization procedures until the parameter estimates do not change between two iterations beyond a threshold. The rate of convergence of the EM algorithm (speed of computation) depends on the initial parameter estimates and the signaltonoise ratios (SNR) of the rhythms present in the signal.
In simulations, unless otherwise specified, we use the first 2 s of data to fit the model parameters. With these model parameters estimated and fixed, we then apply a Kalman filter to predict and update the state estimates, and estimate the phase and amplitude for each oscillator, each representing a different rhythm, for every sample (Figure 1). We note that the Kalman filter is an optimal filter, that is, in the presence of Gaussian noise, the Kalman filter achieves minimum error and is an unbiased estimate of the mean phase. However, if the noise is not Gaussian, the Kalman filter remains the optimal linear filter for phase estimation under this model. Given the stochastic frequency modulation that is possible under the model, the Kalman filter is robust to small shifts in the central frequency of the rhythm.
We now define the stages of Kalman filtering. The state is initialized with zeros and the state covariance to a diagonal matrix with 0.001 along the diagonal. To predict the future state and state estimation error at time t prior to filtering (xtt−1 and Ptt−1) given the filtered state and estimation error ($x}_{t1}^{t1$ and $P}_{t1}^{t1$) from the previous sample (t−1), we compute:
We then estimate the Kalman gain for the current sample (t):
Finally, we filter the current sample using the Kalman gain to compute filtered state ($x}_{t}^{t$) and state estimation error ($P}_{t}^{t$) at time t:
The SSPE method also tracks the state estimation error ($P}_{t}^{t$), which allows estimation of confidence in the mean phase estimates. To do so, we sample 10,000 samples from the posterior distribution of each state j, and from this estimate the width of the 95% credible intervals for the phase. The credible intervals are derived from the posterior distribution of the state, which is Gaussian; however, following their transformation into angles, the intervals are not constrained to be symmetric. Note that, alternatively, one could attempt to define a closed form for the phase distribution to estimate confidence bounds, as we estimate the posterior distribution for the analytic signal under the Kalman filter. However, this strategy requires an expectation of asymptotic convergence and the assumption of independence of amplitude and phase (Withers and Nadarajah, 2013).
ARbased forecasting to estimate phase
Request a detailed protocolWhile current techniques depend upon filtering to isolate a narrowband rhythm from which to estimate phase, different methods exist to forecast the phase. One category of methods forecasts the phase as a linear extrapolation from the current moment while assuming a singular central frequency, which is estimated from the data (Rodriguez Rivero and Ditterich, 2021; Rutishauser et al., 2013). Other methods estimate phase by using an AR model. After estimating the AR parameters (either asynchronously or synchronously), these methods forecast data to limit filter edge effects, then use a filter and Hilbert transform to estimate phase (see Blackwood et al., 2018, and Zrenner et al., 2020). We compare the SSPE to this latter set of techniques as these are more flexible to variability in rhythms than methods that assume a constant central frequency when estimating phase. We implement the version of the Zrenner et al., 2020, algorithm accessible at https://github.com/bnplab/phastimate: and term this method as Zrenner. Please refer to Zrenner et al., 2020, for algorithm details. The Zrenner algorithm uses an FIR filter (applied forward and backward) and the Hilbert transform as the core phase estimation technique. The AR model used to forecast the data and allow realtime phase estimation (by expanding the window for the FIR filter) is fitted to the 30th order and efficiently updated on every new sample. Compared to the default settings, we set the algorithm to track a 6 Hz rhythm by increasing the window size for the FIR filter to 750 ms, the filter order to 192, and the frequency band to 4–8 Hz. We utilize the implementation of the Blackwood et al., 2018, algorithm as provided by the authors. In Blackwood et al., 2018, the authors develop a forward bandpass filter – ARbased (fifthorder AR model whose parameters are fit every second) forecasting and then the Hilbert transform – to compute the phase. Here, we instead apply a Hilbert transformer (of order 18) that estimates the phase by directly computing the Hilbert transform while filtering the data as implemented in Blackwood, 2019. We do so because this implementation is more computationally efficient, functionally equivalent to the original method, and the code is made widely available. We verified in sample data that application of the original and modified Blackwood methods produced consistent phase estimates. We refer to this method as Blackwood.
Noncausal phase estimation
Request a detailed protocolExisting studies compare causal, realtime phase estimation approaches to an acausal estimate of phase computed by applying an (acausal) FIR filter and the Hilbert transform (Blackwood et al., 2018; Chen et al., 2013; Zrenner et al., 2020). To maintain consistency with the existing literature, we do so here and apply a leastsquares linearphase FIR filter of order 750, with bandpass 4–8 Hz. After forward and backward filtering of the data (using MATLAB’s filtfilt function), we apply the Hilbert transform to compute the analytic signal (using MATLAB’s hilbert function). We estimate phase using the fourquadrant arctan function implemented by MATLAB’s angle function and estimate the amplitude envelope using MATLAB’s abs function. We use the angle function throughout the analysis to estimate the phase from a complex valued analytic signal. This acausal estimate of phase has been proposed to serve as a lower bound in phase estimation error for methods that apply a filterbased approach to realtime phase estimation (Zrenner et al., 2020).
Defining error as circular standard deviation
Request a detailed protocolIdentifying a common error metric across alternative algorithms to causally estimate phase is crucial for accurate comparison. We follow Zrenner et al., 2020, and use the circular standard deviation of the difference between the estimated phase and true phase:
where n is number of samples, θ is true phase, $\hat{\theta}$ is the phase estimate, and   indicates the absolute value. We transform the final result from radians to degrees. The circular standard deviation (a transformed version of the phase locking value) captures the variance in the circular error for the phase estimate and ignores any potential bias (nonzero mean shift).
Simulated data
Narrowband to broadband rhythms
Request a detailed protocolOur initial simulation tests the capability of different realtime phase estimators to track different frequency band profiles. We shift the bandwidth of the target rhythm from narrow (sharply peaked in the spectral domain at 6 Hz) to wide (broadly peaked in the spectral domain; peak bandwidth 4–8 Hz). We simulate all rhythms for 10 s at a sampling rate of 1000 Hz, with a central frequency of 6 Hz. For each simulated rhythm, we apply each phase estimation method, as described in Results.
The first type of rhythm we simulate is a pure sinusoid with added white noise:
where T is time, and N(0,1) indicates a standard normal distribution. The second type of rhythm we simulate is a pure sinusoid with added pink noise, that is, the power spectral density (PSD) of the noise decreases with the reciprocal of frequency (f):
The third type of rhythm we simulate is filtered pink noise with added pink noise. To create a broadband spectral peak, we apply the same FIR filter utilized in the acausal phase estimation analysis (see Materials and methods: Noncausal phase estimation) to pink noise ($PSD$~ $1/{f}^{1.5}$), normalizing by the standard deviation and amplifying the resulting signal 10fold. The resulting signal consists of a broad spectral peak centered at 6 Hz. To this broadband signal we again add independent pink noise (${\u03f5}_{1}$) as in Equation 13.
For the final rhythm we simulate a signal ${y}_{t}$ from the SSPE model (see Materials and methods: State space model framework) which is nonsinusoidal and broadband. The state in the state space model is an AR model of order 2 (West, 1997) (equivalently, a damped harmonic oscillator driven by noise), which has been proposed to capture features of gamma rhythms in visual cortex and rhythms in EEG (Franaszczuk and Blinowska, 1985; Xing et al., 2012; Burns et al., 2010; Spyropoulos et al., 2019). The addition of observation noise under the state space model can be analogized to measurement noise. Thus, the state space model provides a simple possible model for brain rhythms observed in electrophysiological data. For the model, we fix one latent state (j = 1) with frequency 6 Hz (${\omega}_{j}=6$), and set a (scaling factor) = 0.99, variance Q = 10, and observation noise variance ($\sigma}_{R}^{2$) to 1. For this rhythm, phase was estimated from the state directly as described in Materials and methods: State space model framework. We show example instances of all four simulated rhythm types in Figure 2A.
Two simultaneous rhythms with nearby frequencies
Request a detailed protocolIn many physiologic use cases (Tort et al., 2010; VanRullen, 2016), we may wish to track one rhythm while ignoring another at a nearby frequency. This may occur because of low spatial resolution (e.g., from mixing of signals due to volume conduction; Nunez and Srinivasan, 2006) or because there exist overlapping neural populations expressing different rhythms (Kocsis et al., 1999). To simulate this case, we consider two sinusoids. We define the target rhythm to track as a 6 Hz sinusoid with consistent amplitude. We define the confounding rhythm as a sinusoid with varying frequency and amplitude. We define the observation (signal plus noise) as:
where $A\in \left[\mathrm{0.2,2}\right]$ and $F\in \left[\mathrm{1,11}\right]$ . The confounding rhythm frequency assumes a range of integer values from 1 to 11 Hz. We note that replacing the constant phase shift ($\pi /4$) with a randomly selected (uniform $\left[\pi ,\pi \right]$) phase shift produces qualitatively consistent results. The confounding rhythm amplitude ranges from 0.2 to 2 times the amplitude of the target rhythm, in increments of 0.2.
Phase reset
Request a detailed protocolOngoing slow rhythms in the brain occasionally undergo a sudden shift in phase. This can occur both in response to a stimulus or as a result of ongoing spontaneous neural dynamics (Voloh and Womelsdorf, 2016). To simulate sudden changes in a rhythm’s phase, we implement a pure sinusoid (6 Hz) with additive pink noise and impose phase slips at four distinct time points, separated by at least 1 s. We fix the extent of phase slip, relative to the previous phase, at $\pi /2$ . However, the phase at which the slip occurs differs across the four time points (see example in Figure 4A). We simulate the signal in five intervals, as follows:
Signaltonoise ratio
Request a detailed protocolTo investigate the impact of SNR on phase estimation, we again simulate a filtered pink noise signal (centered at 6 Hz) with added pink noise (the third type of rhythm simulated in Materials and methods: Narrowband to broadband rhythms). To modulate the SNR, we first compute:
where ${\sigma}_{signal}$ (${\sigma}_{noise}$) is the standard deviation of the filtered pink noise signal (added pink noise). We then divide the filtered pink noise signal by SNR_{0}, multiply the filtered pink noise signal by the chosen SNR (1, 2.5, 5, 7.5, 10), and add this scaled signal to the pink noise.
In vivo data
Request a detailed protocolWe consider two datasets to demonstrate how the SSPE method performs in vivo. The first dataset consists of LFP activity collected from the infralimbic cortex in rats. LFP were acquired continuously at 30 kHz (Open Ephys, Cambridge, MA). An adaptor connected the recording head stage (RHD 2132, Intan Technologies LLC, Los Angeles, CA) to two Millmax malemale connectors (eight channels each, Part number: ED90267ND, DigiKey Electronics, Thief River Falls, MN). For additional details on data collection, please refer to Lo et al., 2020. We analyze here 250 s of LFP data from one channel. Before realtime phase estimation, we downsample the data to 1000 Hz (with a low pass filter at 100 Hz to prevent aliasing) for computational speed. To model these data, we fit three oscillators with the SSPE method, with the goal of tracking a target theta rhythm (defined as 4–11 Hz). Motivated by an initial spectral analysis of these data, we choose one oscillator to capture lowfrequency delta band activity (${\omega}_{j}=1$ Hz, based on spectral peak), one to capture the theta band rhythm of interest (${\omega}_{j}=7$ Hz), and a final oscillator to track highfrequency activity (${\omega}_{j}=40$ Hz). We estimate these model parameters (so center frequencies may change) using the first 10 s of LFP data, and then fix these model parameters to estimate future phase values (up to 250 s). In doing so, we assume that the same oscillator models in the SSPE method remain appropriate for the entire duration of the LFP recording.
The second dataset was derived from the publicly available data in Zrenner et al., 2020. The data consist of spatially filtered and downsampled scalp EEG recordings from electrode C3 (over central areas), which records a strong somatosensory mu rhythm (8–13 Hz). For details on data collection and preprocessing, see Zrenner et al., 2020. We analyze 250 s of data from a single participant for demonstration purposes. For these data, before the realtime phase estimation step, we first apply a moving average causal filter (using 3 s of past data for each new data point) to remove slow drifts. We again track three oscillators in the SSPE method, with the goal of estimating phase from the mu rhythm: one in the delta band (${\omega}_{j}=2$ Hz), one centered at the mu rhythm (${\omega}_{j}=10$ Hz), and one in the beta band (${\omega}_{j}=22$ Hz). We note that these frequency values initialize the model and are then tracked in time. In analyzing the mu rhythm, we ensured that phase was estimated from an oscillator whose frequency was in the mu rhythm range, that is, 8–12 Hz.
For the in vivo data, the true phase is unknown. To estimate phase error, we define the true phase as follows: first, we select a 10 s segment centered on the 1 s interval of the phase estimate. We then fit the state space model (Equations 14) acausally to this 10 s segment and estimate the phase in the 1 s interval at the segment’s center. We note that here, the phase estimate is acausal because we use forward and backward smoothing of the state estimate when estimating the analytic signal. We also note that, in some segments, the frequency band of interest may not be identified through EM due to low signal strength. When this occurred, we omit error in the phase estimates, as these intervals lack the frequency of interest.
Realtime implementation in TORTE
Request a detailed protocolTo make SSPE more broadly available, we implemented SSPE into the realtime TORTE (Toolkit for Oscillatory Realtime Tracking and Estimation) system (https://github.com/tnelab/phasecalculator, 10.5281/zenodo.2633295). TORTE is built upon the Open Ephys platform (https://openephys.org/) and can be used with a broad range of electrophysiologic acquisition devices. Phase outputs and triggers from the system can be used to drive presentation of psychophysical stimuli or delivery of optical/electrical/magnetic stimulation. We directly integrated the MATLABbased SSPE algorithm described above into C++ for improved efficiency. TORTE uses a bufferbased approach to receive data, based on the buffered approach of the underlying Open Ephys system. New buffers are repeatedly received by TORTE containing the most recent neural data. These buffers are subject to system latency causing a delay between activity in the brain being received by the toolkit. Both the buffer size and latency individually vary from microseconds to a few milliseconds depending on the acquisition system used. The phase of each sample within the buffer is estimated using SSPE. A detector continuously monitors the estimated phase and can trigger output logic at any desired phase. We note that the observation noise parameter (Q) is defaulted to 50, but should be updated by the user to the square root of the amplitude. Two modifications were needed to integrate SSPE into TORTE. For computational efficiency, data were downsampled before application of the SSPE algorithm, and then the resulting phase estimates were upsampled back to the native sampling rate. The Open Ephys system can provide buffer sizes that are not exactly divisible into downsampled buffers. When this occurs, the phase estimates at the end of the upsampled buffer cannot be generated from the downsampled SSPE result. To account for this, an additional predicted state sample (with credible intervals) is computed by forward filtering using Equations 07; 08. This extra value is inserted at the end of the downsampled buffer. This provides an extra data point in the SSPE result that can be used to generate the remainder of samples for the upsampled phase output buffer. The second modification results from the realtime data acquisition by TORTE. In this case, low pass filtering must be causal, not the acausal (i.e., forwardbackward filtering) used in the MATLAB implementation. To account for the consistent phase shift introduced in a realtime system, TORTE includes an online learning algorithm that adjusts target phase thresholding to improve event timings in closed loop applications. We test the realtime implementation on an AMD FX8350 CPU with 16 GB RAM.
Code availability
Request a detailed protocolThe code to perform this analysis is available for reuse and further development at MATLAB (https://github.com/EdenKramerLab/SSPEpaper, copy archived at swh:1:rev:6d80fe9c5f610d0dfffe23b5eef2012f780ba621, Mark, 2021b) and Open Ephys (https://github.com/tnelab/phasecalculator, copy archived at swh:1:rev:e11421f23403399ca6b2d85132a31e2c1d5b1397, Mark, 2021a).
Results
We first examine in simulations the viability of different causal phase estimation algorithms, and how these existing methods compare to the proposed method – the SSPE. In these simulations, we vary the nature of the underlying rhythm of interest – either the spectral bandwidth or the noise/confounding signal present. Through these simulations, we examine the assumptions under which different methods perform best at estimating the phase in realtime. We then illustrate the application of the SSPE method to in vivo data in two case studies: a rodent LFP and a human EEG.
SSPE accurately tracks the phase of a single oscillation
To test the accuracy of causal estimates of phase, we simulate rhythms with different spectral profiles in four different ways. For the first two rhythms, we simulate a 6 Hz sinusoid and add either white noise (constant power spectrum, ‘sine in white noise’, Figure 2A.i) or pink noise (power spectrum decreases with the reciprocal of frequency, ‘sine in pink noise’, Figure 2A.ii). For the next simulation scenario, to generate a broadband spectral peak, we filter pink noise into the band of interest using an FIR filter (4–8 Hz, ‘filtered pink noise’, see Materials and methods), amplify this signal, and then add additional (unfiltered) pink noise (Figure 2A.iii). For the last simulation, to generate the signal of interest we sample data from the model underlying the SSPE method. The state space harmonic oscillator model (‘state space model’) consists of a damped harmonic oscillator driven by noise (the signal) with additive white noise (Figure 2A.iv). We note that, for the first two simulation approaches in which the signal is a pure sinusoid, the true phase is known. For the ‘filtered pink noise’ simulation, we estimate the phase by applying the Hilbert transform to the signal before adding noise (see Materials and methods). For the ‘state space model’ simulation, we determine phase directly from the model (see Materials and methods). Thus, in all simulation scenarios, we have an independent measure of the true phase.
We repeat each simulation 1000 times with different noise instantiations and calculate error as the circular variance of the difference between the estimated phase and true phase (see Materials and methods). We find that for rhythms with narrowband spectral peaks (sine in white noise and sine in pink noise, Figure 2B and C), all estimators perform well, and those estimators that presume a narrowband signal (i.e., the filterbased approaches) perform best. However, for rhythms with broad spectral peaks (filtered pink noise and state space model, Figure 2B and C), the SSPE method outperforms the other causal approaches. Moreover, the SSPE method performs as well as (filtered pink noise), or better than (state space model), the acausal FIR method. In what follows, we further explore the implications of this result and the impact on assessment of phase estimator performance. We conclude that the SSPE method performs well in all four simulation scenarios and outperforms existing causal phase estimation methods in two instances with broad spectral peaks.
SSPE accurately tracks the phase of two simultaneous oscillations
To replicate a situation that has been documented in EEG (e.g., multiple alpha rhythms coexisting) (VanRullen, 2016) and in LFP (e.g., multiple theta/alpha rhythms in hippocampus) (Tort et al., 2010), we simulate two rhythms with similar frequencies. In this situation, phase estimation requires an approach to identify the presence of two rhythms and target the rhythm of interest. This simulation examines the impact of choosing a specific, previously specified, frequency response to isolate the rhythm of interest, as required by the two comparative causal phase estimation methods implemented here. In some cases, a finely tuned bandpass filter may separate the two rhythms, although such a filter requires sufficient prior knowledge and may be difficult to construct, particularly if the center frequencies are not perfectly stationary. The SSPE method provides an alternative approach to identify the two (or potentially more) rhythms without filtering and then allows tracking of the appropriate target rhythm. Further, while we focus on the phase of the target rhythm here, we note that the SSPE method provides estimates of the amplitude and phase of all modeled rhythms.
To illustrate this scenario, we simulate two slow rhythms at different frequencies. We set the target rhythm as a 6 Hz sinusoid with fixed amplitude, and the confounding rhythm as a 5 Hz sinusoid with 1.5 times the 6 Hz amplitude (Figure 3A). The observed signal consists of the combined target and confounding rhythms, plus white noise, simulated for 15 s. We apply each phase estimation method to the observed signal, and estimate the phase of the target rhythm after fitting parameters with 5 s of initial data. Visual inspection for this example (Figure 3B) suggests that only the SSPE method accurately tracks the target phase.
Repeating this analysis with the confounding amplitude and frequency varied across a range of values (frequencies 1–11 Hz and amplitudes 20–200% of the target rhythm; see Materials and methods), we find consistent results. For all methods, the error increases as the frequency of the confounding rhythm approaches the frequency of the target rhythm (6 Hz, Figure 3C). However, the increased error is restricted to a narrower frequency interval for the SSPE method. For the existing causal phase estimators, the error likely reflects choices in the filtering procedure (i.e., type of filter, features of passband and stopband, window size, see Sameni and Seraj, 2017) and modeling procedure (i.e., AR order, number of samples forecasted; Chen et al., 2013; Zrenner et al., 2020). Similarly, for the acausal FIR estimator, the error likely reflects the filter passband of 4–8 Hz, with increasing error as the confounding rhythm nears the 6 Hz target rhythm. For the existing causal and acausal methods, error also increases as the confounding oscillation amplitude increases (Figure 3C).
Previously, we showed that all estimation methods performed well for a signal consisting of a single sinusoid in white noise (see Figure 2, sine in white noise). Here, we show that for a signal consisting of two sinusoids in white noise, the SSPE method outperforms the existing methods. Unlike these existing methods, the SSPE method estimates parameters for both the target and confounding rhythms, and tracks each oscillation independently. By doing so, the SSPE method accounts for the confounding effect of the nearby oscillation and allows more accurate tracking of the target rhythm. We conclude that the SSPE method outperforms these existing methods in the case of two simultaneous rhythms with nearby frequencies.
Phase estimation following phase reset
Phase resets are defined as consecutive samples of data that show a sudden shift in the oscillatory phase, proposed to represent the presence of synchronization of a neural mass in response to a stimulus, for example (Lakatos et al., 2009). Phase resets are common and well documented in electrophysiologic recordings (Fiebelkorn et al., 2011; Makeig et al., 2004). Recent work delivering electrical stimulation based on phase of the theta rhythm has been forced to avoid segments with phase resets to reduce error in realtime phase estimation (Gordon et al., 2021). Further, accurate tracking of phase resets would provide better evidence for the relevance of phase resetting in generating event related potentials (Sauseng et al., 2007) – an analysis that could also be performed acausally. We examine the performance of each phase estimation method in the presence of simulated phase resets.
To do so, we simulate phase reset as a shift in the phase of a 6 Hz sinusoid by 90 degrees at four points in a 10 s signal (example in Figure 4A). We then compute the error of each phase estimation method in the 167 ms (one period of a 6 Hz rhythm) following the phase reset. Visual inspection of an example phase reset (Figure 4B) suggests that only the SSPE method accurately tracks the phase in the interval following the reset. We note that we fit SSPE parameters using an initial 2 s interval containing no phase resets (Equation 15). Repeating this analysis for 1000 simulated signals (4000 total phase resets), we find consistent results; the error of the SSPE method is significantly lower than the error of the other methods (p ~ 0 in all cases using twosided ttest, see Table 1). These results demonstrate that the existing causal and acausal methods implemented here are unable to estimate the phase accurately immediately after a phase reset. To determine how quickly each method returned to accurate phase estimation, we first estimated the prestimulus average error for each method, prior to the first phase reset, over a 500 ms interval. We then determined the time between each subsequent phase reset and the moment when the error reduced to 1.5 times the prestimulus error. Note that the level of prestimulus error differs across methods, so, for example, we determine when the SSPE converges to a lower error than the Blackwood method. We find that each method eventually converged to accurate phase estimation after a phase reset, with the SSPE method converging most rapidly (Table 1). The differences between methods likely result from the windowing used by each method when estimating the phase; smaller windows will lead to less impact of the phase reset on the error (as seen with the Blackwood approach converging more quickly than the Zrenner method), with the SSPE needing a window of a single sample.
We note that these existing methods all require a bandpass filtering step, which limits the fastest rate at which the phase can change (instantaneous frequency Boashash, 1991). The SSPE method, however, does not require the phase change within a limited range of instantaneous frequencies and therefore can track rapid changes in phase. We conclude that the SSPE method is able to most rapidly track the true phase after a phase reset and thus provides the most accurate estimate.
Phase estimation with different SNR
To demonstrate the impact of noise on phase estimation, we simulate signals with different SNR. To generate a signal we filter pink noise into the band of interest (see Materials and methods: Signaltonoise ratio), and then add additional (unfiltered) pink noise. The ratio of the standard deviations of the signal and noise defines the SNR; we modulate the standard deviation of the signal to achieve the desired SNR. Finally, we normalize the full observation so that amplitudes across simulations are consistent. In these simulations the SNR varies from 1 (signal and noise have the same total power) to 10 (signal has 100 times the total power of the noise); see example traces of the observation under each SNR in Figure 5A.i–iv, and see example spectra under each SNR in Figure 5—figure supplement 1. We repeat each simulation 1000 times with different noise instantiations for each SNR.
We compute two metrics to characterize the impact of signaltonoise on phase estimation. First, we consider the credible intervals estimated by the SSPE, here measured as the average credible interval width throughout a simulation. As expected, the credible intervals track with the SNR (Figure 5B.i); signals with higher SNR correspond to smaller credible intervals. Computing the average amplitude (as derived from the acausal FIR, see Methods: NonCausal Phase Estimation) throughout each simulation, we find no overlap between the distributions of amplitudes at SNR = 1 and SNR = 2.5 (Figure 5.ii), while the distributions of credible intervals overlap. This suggests that, at low SNR, an arbitrary amplitude threshold (e.g., amplitude >1) would eliminate all data at SNR = 1 and preserve all data at SNR = 2.5. Alternatively, using the credible interval width allows identification of transient periods of confidence (e.g., credible interval width <40 degrees) at both SNR values. We also note that the credible interval width – in units of degrees – provides a more intuitive threshold than the amplitude, whose intrinsic meaning at any particular value is only understood relatively. These results suggest that, in this simulation study, amplitude serves as a less accurate criterion for confidence in phase compared to the credible interval width.
Second, we consider the error in phase estimates using four different phase estimation methods (Figure 5C). We find that the error for the SSPE method reduces more quickly than for the other methods, and at an SNR of 1, the SSPE method is comparable to or better than the contemporary methods of realtime phase estimation. At SNR = 10, the error in the SSPE phase estimation is comparable to the acausal FIR, while the alternative realtime algorithms continue to perform with error similar to an SNR of 1. These results suggest that, for the two existing realtime phase estimation methods, broadband rhythms are difficult to track even with high SNR, consistent with Figure 2C.
Example in vivo application: rodent LFP
Having confirmed the SSPE method performs well in simulation, we now illustrate its performance in example in vivo recordings. We first apply the SSPE method to an example LFP recording from rodent infralimbic cortex (see Materials and methods), which contains frequent intervals of theta band (4–8 Hz) activity. Unlike the simulated data considered above, in the in vivo signal, the target rhythm changes in time (e.g., time intervals exist with, and without, the theta rhythm). Accurate application of the SSPE method requires accurate estimation of the parameters defining the oscillators in the SSPE model. In the simulated data, the stationary target rhythm allowed consistent SSPE model parameter estimation from (any) single interval in time; in those simulations, we estimated the SSPE model parameters from a single, early time interval (see Materials and methods) and then used those (fixed) model parameters for future causal phase estimates.
For the in vivo data, in which the target rhythm may change in time, model estimates from an initial time may be inappropriate – and produce inaccurate phase estimates – at a later time. However, comparing the estimated phase to the true phase (see Materials and methods), we find that the phase error tends to remain consistent in time (blue curve in Figure 6A; circular standard deviation 95% interval = [33.11, 75.39]; for error of phase estimates computed using alternative methods, see Figure 6—figure supplement 1). We find no evidence of a linear trend in error with time (linear fit, slope = 0.01, p = 0.53); that is, we find no evidence that the error increases in time. Instead, visual inspection suggests that the error increases when the spectral content of the LFP changes (e.g., near 120 s in Figure 6B). We note that the error bounds are within the bounds of the range expected for realtime phase estimates in Zrenner et al., 2020. Finally, we note that visual inspection of the PSD (Figure 6C) demonstrates that the theta band power appears more closely matched by a broadband peak, rather than a narrowband peak.
An alternative approach is to reestimate the SSPE model parameters over time. By allowing the parameters of the model oscillators to change in time, this approach may improve the accuracy of the phase estimates. To examine this in the in vivo data, we refit the SSPE model parameters in each 10 s interval of the LFP data. More specifically, we refit model parameters in every 10 s interval of the 250 s recording, incrementing the 10 s interval by 1 s, to create 240 models. We then use each of these models to causally estimate the phase for all future times. For example, we use the SSPE model with parameters estimated from the interval 40–50 s to then estimate the phase for data after 50 s. We compute the error in the phase estimates for all 240 models to assess whether refitting the SSPE model parameters impacts future phase estimates. We find no evidence that the distribution of phase errors changes for the refit models compared to the original fixed parameter model for the immediate 30 s that followed (p = 0.21, Wilcoxon rank sum test) or for the entire time period analyzed (red interval in Figure 6A, circular standard deviation 95% interval = [30.68, 88.7] degrees; p = 0.35 using Wilcoxon rank sum test). The original fixed parameter SSPE model (fit using the first 10 s of data; blue curve in Figure 6A) performs as well as the models refit with future data. We conclude that, under the assumption that parameters for the frequency band of interest are stable, we are able to estimate the SSPE model parameters whenever there exists a consistent rhythm of interest and then apply this model to estimate phase at future times (here on the scale of hundreds of seconds). In other words, for these data, any 10 s interval with sufficiently strong theta rhythm can be used to fit parameters before the model is applied causally. This result has an important practical implication; the SSPE model parameters can be fit using an initial segment of data, and then applied to derive causal estimates of future phase values. In other words, the model has robust estimation properties; parameters estimated early in the recording remain appropriate for future observations.
In addition, these results allow a novel insight into the rhythmic activity of the data, beyond those available using traditional methods. Typical analysis of neural rhythms either presumes the presence of a rhythm (for subsequent filtering) or operates in the frequency domain using the PSD and statistical testing to identify significant rhythms. Here, through direct fitting of the state space model to the time series data, we find consistent model parameters that accurately track the observed theta rhythm. This consistency is maintained through variations in the rhythm’s amplitude. These observations suggest that the theta rhythm in infralimbic LFP maintains a consistent underlying rhythmic structure over the entire duration of the observation.
Confidence in the phase estimate
Past approaches to causally estimate phase have not attempted to quantify confidence in the phase estimate beyond utilizing the amplitude as a surrogate (Zrenner et al., 2020). However, using the amplitude assumes a monotonic linear relationship between the phase confidence bounds and the amplitude, which may not be true (e.g., for an AR(2) process with complex roots, Spyropoulos et al., 2019). Unlike existing methods, the SSPE method provides realtime estimates of the credible interval estimates for the phase (resulting from estimating the parameters of the SSPE model and tracking the phase using the Kalman filter; see Materials and methods). To illustrate this, we show the phase and credible interval estimates of the SSPE method for an example 3 s interval of the LFP data (Figure 7; for comparison to phase estimates computed using alternative methods, see Figure 7—figure supplement 1). We note that, during a prominent theta rhythm, the credible intervals tightly track the mean phase estimate; at other times – when the theta rhythm is less obvious – the credible intervals expand. We can restrict error in the phase estimate by examining only samples with small CI width. When thresholding at 10 degrees of credible interval width, we find that the error decreases from (mean) 46.9 (s.e. 0.89) degrees to 26.91 (s.e. 1.07) degrees, while still retaining 27% of the data. This is within the current stateoftheart in phase estimation with minimum of 20–40 degrees error as a function of SNR (when applying a sufficiently high amplitude threshold, Zrenner et al., 2020).
Previous studies have proposed a relationship between the theta amplitude envelope, estimated using the traditional approach of narrowband filtering followed by a Hilbert transform, and confidence in the phase (Zrenner et al., 2020); when the theta signal is strong (i.e., the amplitude is large), the theta phase estimate improves. To examine this relationship here, we plot the theta amplitude envelope versus the phase credible intervals derived from the SSPE method (Figure 7C.i). Visual inspection reveals that, when the amplitude is large, small decreases in amplitude coincide with increased credible intervals, as expected. However, if we choose to approximate certainty in the phase estimate through thresholding the amplitude, we find that fixed amplitude thresholds produce wide ranges of credible intervals (Figure 7C.ii). This is especially apparent at lower amplitude thresholds; choosing a threshold of 65% results in credible intervals that range from 3.8 (2.5th percentile) to 96.8 (97.5th percentile) degrees. Therefore, a fixed amplitude threshold does not guarantee an accurate phase estimate. Instead, we propose that inspection of the credible intervals better guides certainty in the phase estimate. These credible intervals are directly derived from the modeling procedure, as opposed to surrogate measures of confidence in the phase estimate, such as amplitude. In addition, credible intervals are directly interpretable unlike amplitude thresholds, which depend on the measurement scale (e.g., a credible interval threshold of 10 degrees) versus an amplitude threshold of 10 mV. We conclude that amplitude is linked to confidence in the phase, however, a fixed amplitude threshold is not directly interpretable and encompasses a wide range of credible interval values.
Example in vivo application: human EEG
To illustrate application of the SSPE method in another data modality, we analyze an example EEG recording of the mu rhythm (8–13 Hz) from Zrenner et al., 2020. To do so, we apply the SSPE method to track three oscillators (1, 10, 22 Hz) in 10 s intervals (see Materials and methods),and examine the phase estimates of the 8–12 Hz oscillator. Consistent with the rodent LFP data, we find that refitting of the SSPE model parameters is not required. Instead, the SSPE model parameters estimated from an initial 10 s segment of the EEG produce phase estimation errors consistent with SSPE models refit at future times (Figure 8A). Like the theta rhythm in the LFP rodent data, the mu rhythm in the human EEG data maintains a consistent spectral profile in time (Figure 8B) and further, it appears to be narrowband (Figure 8C). Therefore, estimating SSPE model parameters early in the data accurately captures the mu rhythm tracked later in the data, and there is no need to perform the (computationally expensive) refitting of the model.
The phase estimates and credible intervals for the mu rhythm reveal how this rhythm waxes and wanes in time (example in Figure 8D and E; for comparison to phase estimates computed using the alternative methods, see Figure 8—figure supplements 1 and 2). When a strong mu rhythm emerges, tight credible intervals surround the mean phase estimate tracking the rhythm of interest; as the mu rhythm wanes, the credible intervals expand. When we restrict analysis to samples with credible interval width less than 25 degrees, the error drops from 40.5 (s.e. 1.1) to 11.2 (s.e. 0.46) degrees while still retaining 16.4% of the data. As with the LFP, we are able to assess certainty in our phase estimates using the credible intervals. Finally, we note that computing the state estimate including both the 10 and 22 Hz components of the model better approximates the nonsinusoidal waveform of the mu rhythm (Figure 8—figure supplement 3).
Prestimulus mu rhythm phase predicts MEP amplitude
To illustrate the utility of realtime phase estimation, we again analyze an example EEG recording from Zrenner et al., 2020, in which the authors stimulate a subject at random intervals using transcranial magnetic stimulation to elicit a motor evoked potential (MEP) from the right abductor pollicis brevis muscle. We analyze the MEP amplitude to assess the influence of prestimulus mu rhythm phase on the response (replicating the analysis performed by Zrenner et al., 2020). We estimate mu rhythm phase using the SSPE method with model parameters derived from a different subject. We find that the prestimulus phase of the mu rhythm, at the instant immediately prior to stimulation, predicts (using circular linear regression) 13% of the logtransformed MEP amplitude (Figure 9A). Furthermore, thresholding the trials using credible interval width, we find improved prediction of the MEP amplitude. When using trials with CI width less than 50 degrees, the mu rhythm phase explains 30% of the variance in the MEP amplitude (Figure 9B), and with CI width less than 25 degrees, 48% of the variance (Figure 9C). These results show that modulation of the MEP amplitude depends on the prestimulus phase of the mu rhythm, and suggest that accurate estimation of the mu rhythm (e.g., using the SSPE method with narrow CI width) may enhance the evoked motor response.
Speed and accuracy of TORTE realtime implementation
To assess the validity of the integration of SSPE into TORTE (Figure 10A), we generated a 5 Hz sine wave with no additive noise and estimated the phase using the MATLAB and C++ implementations. The parameter estimates for the fitted model were computed on the first 10 s of the data. The computed parameters were nearly identical between the two implementations and the time to estimate the parameters (which, as above, would need to be done only once per experiment) ranged from 2 to 100 s. The phase estimates of the two implementations have a mean circular difference of 2.06 degrees (Figure 10B) with a circular standard deviation of 0.65 degree; this variation in phase estimates likely arises from differences in filtering in the TORTE realtime (causal filtering on individual buffers) implementation and the offline MATLAB realtime (acausal filtering across all data) implementation. Further, given that it is a near constant shift in phase, the online learning algorithm in TORTE will correct for the biased phase estimates. We also examined the speed of evaluating different sized buffers with SSPE within TORTE. Across five buffer sizes from 420 to 1700 samples (encompassing 9–37 ms of data, see Figure 10C), estimating phase with SSPE took less than 0.35 ms in the worst case. We note that the system latency of TORTE is in the range of a few milliseconds (driven mainly by the USBbased implementation of the current Open Ephys system), whereas the buffers were evaluated on the scale of microseconds. That is, the realtime bottleneck lies in signal amplification and digitization, not this processing. In conclusion, we present a realtime implementation of SSPE that could be easily deployed in various experimental setups.
Discussion
We have introduced a realtime phase estimator derived from a time series state space model. This approach – the SSPE method – addresses limitations of existing methods (namely, the requirements associated with bandpass filtering in phase estimation) and provides a realtime estimate of confidence in the phase estimate. We demonstrated the benefits of the SSPE method in simulations with varying breadth of the peak rhythm frequency and SNR, a confounding rhythm with a nearby frequency, and phase resets. In two case studies of in vivo data, we showed that reestimation of model parameters is unnecessary and that the SSPE method provides clear estimates of confidence in phase estimates. We propose that the SSPE serves as a useful tool for realtime phase applications, for example, realtime phaselocked stimulation or realtime phaseamplitude synchrony estimation.
Narrowband versus broadband signals
While brain rhythms tend to be broadband (Buzsaki, 2004; Roopun et al., 2008), the most common methods in neuroscience to track a rhythm’s phase require appropriate choices of filter parameters (e.g., selection of bandwidth) and generally expect a narrowband oscillation. However, realtime methods that require a narrowband oscillation struggle to capture the phase accurately for broadband rhythms (e.g., Figure 2). Moreover, when the broadband rhythm has nonnegligible power across all frequencies (in which case the phase cycle is more unstable), the standard nonrealtime approach (acausal FIR and Hilbert transform) tracks the phase with limited accuracy. For forecasting needed in realtime analyses, implementing the same assumption (of a broadband peak) in the model may improve phase estimation accuracy. In the SSPE method, the assumption of a narrowband oscillation is not required, as the covariance allows flexibility in modeling peaks of different bandwidths.
Noise in the signal
By applying a bandpass filter directly on the observed signal, existing realtime phase estimators do not disambiguate the signal and the noise, thus folding the noise into the phase estimate. The impact of structured noise (e.g., a competing rhythm or sudden phase shift) is particularly visible in Figure 3 and Figure 4. To address this, the SSPE approach explicitly separates the signal and noise, and models two separate components of noise: state noise and observation noise. The SSPE state consists of harmonic oscillators driven by (white) noise. The resulting spectrum consists of broadband peaks at the oscillator frequencies upon a brown aperiodic component, similar to spectra observed in neural field activity. Confounding rhythms (a particular type of noise) can be explicitly modeled and accounted for in the SSPE approach. We consider the observation noise as additive white noise (flat spectrum) due to unknown noise in the measurement process. This model structure is more flexible than methods reliant on bandpass filters, allowing different expectations for the noise (e.g., modeling a confounding rhythm as a pure harmonic oscillator) or adding additional model elements to account for power at higher frequencies (e.g., adding moving average terms; Matsuda and Komaki, 2017).
Detecting a rhythm for phase estimation
To estimate phase for a target rhythm, there should be evidence that the rhythm exists in the data. A few contemporary realtime phase estimators tackle this problem by using an amplitude threshold to decide when a rhythm is present (Rodriguez Rivero and Ditterich, 2021; Rutishauser et al., 2013). An amplitude criterion depends on accurate knowledge of the distribution of amplitudes possible for the data, relative to the background noise. The SSPE method, in contrast, directly represents the presence/absence of a rhythm in the credible intervals for the phase, as demonstrated for the in vivo data. The result is directly interpretable and permits a more intuitive decision metric (e.g., a threshold depending on the credible interval width) than an arbitrary amplitude criterion. Further, while we expect a central frequency for a rhythm estimated from the data, the Gaussian assumptions of the harmonic oscillator driven by noise allow for the instantaneous frequency to be stochastic. This permits greater flexibility in tracking a rhythm’s central frequency compared to prior phase estimators and supports analysis of broadband rhythms with drifting central peaks, and abrupt changes in phase of a fixed rhythm (e.g., Figure 4). However, the SSPE method does require identification of target and confounding rhythms for accurate estimation. If details about the rhythms (central frequencies and bandwidth) are unknown, or the target rhythm is known to exist only for short intervals in a long recording, then alternative methods that allow realtime visualization and userdirected intervention in estimation may be pursued (Rodriguez Rivero and Ditterich, 2021; Rutishauser et al., 2013).
Computational speed
For realtime phase applications, speed of estimation is critical. We focus here on algorithmic speed as opposed to hardware speed, which is a function of the computing device. Algorithmic speed can be measured at two different points: parameter estimation and application. In the SSPE method and other realtime phase estimators (Chen et al., 2013; Zrenner et al., 2020), the parameter estimation stage is the slower aspect of phase estimation. Past work has required 10 s of minutes for parameter estimation (Chen et al., 2013). Here, for the SSPE method using a window size of 10 s and a sampling rate of 1 kHz, we find that parameter estimation time rarely exceeds a minute or two. In the application stage, we find that speed is orders of magnitude faster than the sampling rate (Figure 8), considerably improving on the algorithmic delays present in existing algorithms (Rodriguez Rivero and Ditterich, 2021; Rutishauser et al., 2013).
Limitations and future directions
Like all realtime phase estimators, the proposed method possesses specific limitations and opportunities for future improvement. First, the SSPE method treats the generative model for any rhythm as a harmonic oscillator driven by noise. A pure sinusoid is outside this model class, and as such the SSPE method performs slightly worse than filterbased algorithms at tracking the phase of a pure sine. Nonetheless, the error remains small for the SSPE tracking a pure sinusoid in noise (Figure 2). We note that a simple alteration in the SSPE model (removing the state covariance Q) would allow more accurate tracking of a sinusoid, but not the broadband rhythms more commonly observed in vivo. Second, slow drifts (<0.1 Hz) in the baseline of the data can impair the ability of the SSPE method to track the phase. Slow drifts could be modeled as an additional confounding oscillator, however, estimating its presence would require long segments of time (possibly exceeding 40 s). Alternatively, the SSPE model could be updated to include a model of the slow drift dynamics (perhaps as a piecewise linear term), or a moving average mean trend could be causally removed as we demonstrated with the in vivo EEG analysis. Third, in applying the SSPE method, we assume that the model we estimate for a rhythm remains appropriate for the extent of the experimental duration. While we showed in two in vivo case studies that this assumption may be reasonable, in other cases this assumption may fail. Rhythms, observed over time, may be better represented by models with changing parameters. Indeed, nonstationarity is an important issue to consider when tracking brain rhythms. The SSPE method is robust (by virtue of the model structure and application of the Kalman filter) to small changes in the central frequency or the bandwidth of a rhythm (e.g., Figure 2). However, nonstationary rhythms require new algorithms to be developed (such as the PLSO method in Song et al., 2020). An extension of the SSPE method that could potentially track changing central frequencies is to apply an extended Kalman filter (Schiff, 2012) that simultaneously estimates the frequencies of interest while filtering the state. Alternatively, the SSPE method could also be extended to implement a switching model that utilizes multiple sets of parameters and switches between parameter sets as necessary, or by refitting the SSPE model as time evolves. Finally, we note that data sampled at a high rate can reduce computational efficiency when fitting parameters for SSPE. To address this, a quasiNewton optimization – instead of EM, as implemented here – may improve performance, as suggested in Matsuda and Komaki, 2017. Determining the optimization framework best suited for practical use remains an important direction for future work.
Conclusion
We introduced a new algorithm – the SSPE method – to track phase in realtime that addresses several limitations of current methods. We demonstrated in simulation that SSPE performs better than current stateofthe art techniques for realtime phase estimation and comparably, or better, than a commonly used acausal FIR approach. We showed in in vivo data that we can estimate SSPE model parameters and apply them for hundreds of seconds without reestimation. Finally, we showed how credible intervals from the SSPE method track confidence in the realtime phase estimate, and may provide an intuitive, userdetermined decision metric for stimulation at a particular phase. We have provided reference implementations suitable for immediate incorporation into phaseaware in vivo experiments. Future work using the SSPE method in practice may help provide critical evidence on the potential functional importance of phase in neural dynamics.
Data availability
All data generated or analyzed during this study, or were used to create the figures are included in the supporting files, or are available through already public data archives (https://gin.gnode.org/bnplab/phastimate, https://doi.org/10.6084/m9.figshare.14374355).
References

ConferenceAnnual International Conference of the IEEE Engineering in Medicine and Biology Society. IEEE Engineering in Medicine and Biology Society. Annual International ConferenceContinuous Phase Estimation for PhaseLocked Neural Stimulation Using an Autoregressive Model for Signal Prediction. pp. 4736–4739.https://doi.org/10.1109/EMBC.2018.8513232

Estimating and interpreting the instantaneous frequency of a signalpart 1: FundamentalsProceedings of the IEEE 80:520–538.https://doi.org/10.1109/5.135376

Recursive bayesian decoding of motor cortical signals by particle filteringJournal of Neurophysiology 91:1899–1907.https://doi.org/10.1152/jn.00438.2003

The Phase of Ongoing EEG Oscillations Predicts Visual PerceptionJournal of Neuroscience 29:7869–7876.https://doi.org/10.1523/JNEUROSCI.011309.2009

Neuronal Oscillations in Cortical NetworksScience 304:1926–1929.https://doi.org/10.1126/science.1099745

BookRhythms of the BrainOxford, UK: Oxford University Press.https://doi.org/10.1093/acprof:oso/9780195301069.001.0001

The functional role of crossfrequency couplingTrends in Cognitive Sciences 14:506–515.https://doi.org/10.1016/j.tics.2010.09.001

RealTime Brain Oscillation Detection and PhaseLocked Stimulation Using Autoregressive Spectral Estimation and TimeSeries Forward PredictionIEEE Transactions on Biomedical Engineering 60:753–762.https://doi.org/10.1109/TBME.2011.2109715

Brain Oscillations and the Importance of Waveform ShapeTrends in Cognitive Sciences 21:137–149.https://doi.org/10.1016/j.tics.2016.12.008

BookCharacterizing Complex, MultiScale Neural Phenomena Using StateSpace ModelsIn: Chen Z, Sarma SV, editors. Dynamic Neuroscience. Springer International Publishing. pp. 29–52.https://doi.org/10.1007/9783319719764_2

Linear model of brain electrical activity—EEG as a superposition of damped oscillatory modesBiological Cybernetics 53:19–25.https://doi.org/10.1007/BF00355687

Prefrontal attentional saccades explore space rhythmicallyNature Communications 11:925.https://doi.org/10.1038/s41467020146497

Prefrontal ThetaPhase Synchronized Brain Stimulation With RealTime EEGTriggered TMSFrontiers in Human Neuroscience 15:691821.https://doi.org/10.3389/fnhum.2021.691821

Rhythmic facilitation of sensory processing: A critical reviewNeuroscience & Biobehavioral Reviews 86:150–165.https://doi.org/10.1016/j.neubiorev.2017.12.002

Neural CrossFrequency Coupling: Connecting Architectures, Mechanisms, and FunctionsTrends in Neurosciences 38:725–740.https://doi.org/10.1016/j.tins.2015.09.001

A New Approach to Linear Filtering and Prediction ProblemsJournal of Basic Engineering 82:35–45.https://doi.org/10.1115/1.3662552

Interdependence of multiple theta generators in the hippocampus: A partial coherence analysisThe Journal of Neuroscience 19:6200–6212.https://doi.org/10.1523/JNEUROSCI.191406200.1999

BookCase Studies in Neural Data Analysis: A Guide for the Practicing NeuroscientistCambridge, MA: MIT Press.

Prestimulation phase predicts the TMSevoked responseJournal of Neurophysiology 112:1885–1893.https://doi.org/10.1152/jn.00390.2013

Some sampling properties of common phase estimatorsNeural Computation 25:901–921.https://doi.org/10.1162/NECO_a_00422

Paired Electrical Pulse Trains for Controlling Connectivity in EmotionRelated Brain CircuitryIEEE Transactions on Neural Systems and Rehabilitation Engineering 28:2721–2730.https://doi.org/10.1109/TNSRE.2020.3030714

Mining eventrelated brain dynamicsTrends in Cognitive Sciences 8:204–210.https://doi.org/10.1016/j.tics.2004.03.008

Diverse Phase Relations among Neuronal Rhythms and Their Potential FunctionTrends in Neurosciences 39:86–99.https://doi.org/10.1016/j.tins.2015.12.004

SoftwareAnalytic signal plugin, version swh:1:rev:e11421f23403399ca6b2d85132a31e2c1d5b1397Software Heritage.

Time Series Decomposition into Oscillation Components and Phase EstimationNeural Computation 29:332–367.https://doi.org/10.1162/NECO_a_00916

BookElectric Fields of the Brain: The Neurophysics of EEGOxford, UK: Oxford University Press.

A userfriendly algorithm for adaptive closedloop phaselocked stimulationJournal of Neuroscience Methods 347:108965.https://doi.org/10.1016/j.jneumeth.2020.108965

Temporal interactions between cortical rhythmsFrontiers in Neuroscience 2:034.https://doi.org/10.3389/neuro.01.034.2008

Phase synchronization in driven and coupled chaotic oscillatorsIEEE Transactions on Circuits and Systems I 44:874–881.https://doi.org/10.1109/81.633876

A method for closedloop presentation of sensory stimuli conditional on the internal brainstate of awake animalsJournal of Neuroscience Methods 215:139–155.https://doi.org/10.1016/j.jneumeth.2013.02.020

A robust statistical framework for instantaneous electroencephalogram phase and frequency estimation and analysisPhysiological Measurement 38:2141–2163.https://doi.org/10.1088/13616579/aa93a1

BookNeural Control Engineering: The Emerging Intersection between Control Theory and NeuroscienceCambridge, MA: MIT Press.

An Approach to Time Series Smoothing and Forecasting Using the Em AlgorithmJournal of Time Series Analysis 3:253–264.https://doi.org/10.1111/j.14679892.1982.tb00349.x

Open Ephys: An opensource, pluginbased platform for multichannel electrophysiologyJournal of Neural Engineering 14:045003.https://doi.org/10.1088/17412552/aa5eea

Cortical Networks Produce Three Distinct 7–12 Hz Rhythms during Single Sensory Responses in the Awake RatThe Journal of Neuroscience 30:4315–4324.https://doi.org/10.1523/JNEUROSCI.605109.2010

Perceptual CyclesTrends in Cognitive Sciences 20:723–735.https://doi.org/10.1016/j.tics.2016.07.006

A Role of PhaseResetting in Coordinating Large Scale Neural Networks During Attention and GoalDirected BehaviorFrontiers in Systems Neuroscience 10:18.https://doi.org/10.3389/fnsys.2016.00018

The distribution of the amplitude and phase of the mean of a sample of complex random variablesJournal of Multivariate Analysis 113:128–152.https://doi.org/10.1016/j.jmva.2012.05.017

Stochastic Generation of GammaBand Activity in Primary Visual Cortex of Awake and Anesthetized MonkeysJournal of Neuroscience 32:13873–13880.https://doi.org/10.1523/JNEUROSCI.564411.2012

Decoding Hidden Cognitive States From Behavior and Physiology Using a Bayesian ApproachNeural Computation 31:1751–1788.https://doi.org/10.1162/neco_a_01196
Decision letter

Frances K SkinnerReviewing Editor; Krembil Research Institute, University Health Network, Canada

Laura L ColginSenior Editor; University of Texas at Austin, United States

Milad LankaranyReviewer; The Hospital for Sick Children, University of Toronto, Canada
Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Acceptance summary:
This paper will be of interest to theoretical and experimental neuroscientists. It proposes a novel approach based on state space modeling to estimate the phase of neural signals like EEG and LFP in realtime. It is expected to help drive the improvement and proliferation of phaserelated concepts in neurobiology.
Decision letter after peer review:
Thank you for submitting your article "A State Space Modeling Approach to RealTime Phase Estimation" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Laura Colgin as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Milad Lankarany (Reviewer #1).
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential revisions:
1) Consider the impact of noise in the phase estimation, include specifications of the Kalman filter and its robustness, and consider the performance of the estimated phase relative to other methods.
2) Demonstrate the advantage of using this method for exploring the link between rhythms and behaviour and for phaselocked stimulation, and separate applications that require real time phase estimation from those that can be performed offline using acausal methods.
3) Provide some help to readers from nontechnical background by linking the method back to applications and invivo data more.
4) Expand details regarding computing credible intervals, slow drift considerations, and physiological motivation of the simulations.
Further details are provided below in the 'Recommendations for the authors'.
Reviewer #1 (Recommendations for the authors):
I believe the proposed method in this paper is very novel and timely for the neuroscience community. The paper is well organized and written very clearly. I have some comments that might improve the quality of this paper.
1. Impact of Noise in the phase estimation. Although there is a section in the Discussion related to Noise in Signal, the authors did not address the impact of Signal to Noise Ratio (SNR) on the performance of their algorithms. I think this will add value to your algorithm, specifically because the Kalman Filter can be much more reliable than other methods in this context.
2. Specifications of the Kalman Filter and its robustness. specifications of the Kalman filter (KF) should be described in the Method Section. For example, the choice of state and observation covariance matrices (Q and R) should be clarified. In this regard, the authors can discuss the robustness of the KF in phase estimation.
3. For invivo data, how the performance of the estimated phase might be compared to other methods (used in the simulation study)? It can be helpful to visually compare the estimated phase calculated by SSPE and other methods which were used in the simulation study.
Reviewer #2 (Recommendations for the authors):
I think the manuscript can be improved by
1) clearly demonstrating the advantage of using this method for exploring the link between rhythms and behaviour and for phaselocked stimulation. The authors should separate applications that require real time phase estimation from those that can be performed offline using acausal methods.
2) Building on my previous point:
– In figure 2, proposed method clearly outperforms others for signals generated using the "state space model". How representative is this for invivo recordings?
– In figure 4, it is demonstrated that other methods cannot track phasereset, however, is this offset sustained or is it a transient slip in phase estimation for one cycle. In other words, if we zoom out in time, do the other methods converge back to true phase? And if so, when would it be critical to track phasereset for one cycle?
– The performance of SSPE should be demonstrated with respect to other methods for invivo data (figures 5 and 7) so that the reader can evaluate the utility of this method.
3) The paper is well written but comes across as very technical. Readers from nontechnical backgrounds may struggle with following some of the concepts. This could be improved by linking the method back to applications and invivo data more, and also by reviewing the language used in the methods segment.
Reviewer #3 (Recommendations for the authors):
1) Credible intervals
An interesting aspect of the manuscript is the possibility of computing credible intervals, making it possible to establish bounds on time points when the phase can be estimated with sufficient accuracy. In previous approaches, this is usually done via an amplitude threshold, with the assumption that a high amplitude will also result in high phase estimation accuracy. In Figure 6 the authors therefore study the relationship between amplitude and credible intervals. It is tough to interpret this evidence; the argument is that for an amplitude threshold using the 65th percentile the range of credible intervals is from the minimum to maximum possible. But looking at the major bulk of probability mass, it's centered around intervals with an acceptable looking range? Could more details regarding computation of confidence intervals be provided, e.g. are they symmetric? Figure 6C is a scatter plot with lots of points that are correlated, maybe it would be good to visualize this in a different way reducing the number of points per cycle? Can the credible intervals be related to error in degrees? They look pretty negligible in the left part of the figure (given that the realtime procedure also introduces a lag, Figure 8B),is the plotted xrange the relevant one? How does SNR/amplitude relate to phase accuracy to credible intervals? If the spread of credible intervals can be largely explained by amplitude (aka high correlation between credible interval width and amplitude), then maybe this is not such a big factor?
2) Slow drifts
A practical reason why existent algorithms rely on buffering is the frequent existence of slow drifts in the data, for which detrending needs to be applied. For the EEG data, the following methodological detail of the current approach is given "Before realtime phase estimation, we highpass filter the data at 0.5 Hz to remove slow drifts." This procedure is not possible when estimating the phase in realtime. The authors claim the benefit of their approach is the absence of needing buffers, but in the light of this I am not sure about whether this is the case. In the discussion, the authors discuss the possibility of modeling the slow drifts via an additional rhythm in the state space, but state that long segments of data are needed for this. In my view, this is a crucial aspect of any algorithm targeting realtime phase and more information and analysis resolving slow drifts would greatly improve the actual applicability of the present approach.
3) Data generation process and simulations
I find some of the simulations not sufficiently motivated from a physiological perspective.
– Two simultaneous rhythms with nearby frequencies: In this scenario the authors simulate the signal as a summation of two sinusoids with a frequency offset with a fixed phase shift. This will result in a signal with amplitude modulation with the envelope given by the difference of the frequencies, resulting in the appearance of oscillatory bursts. I would argue that for physiological signals that present like this, it is not of interest to track a single frequency in this signal, as e.g. the trough phase of this compound signal corresponds to a specific physiological state. I would also suggest showing the amplitude for this case in Figure 3. Possibly the simulation could be modified to be better physiologically motivated, as they are certainly cases where there are neighboring rhythms which present on the same electrode through volume conduction, but it is unclear (but would be of interest) to me if the proposed method could be helpful in this case.
– From my perspective phase is only welldefined if there are clearcut oscillatory bursts. In that sense, the signal given as an example in Figure 2A iv does not have clearcut oscillatory states, with no clear spectral peak and the true phase (blue curve) changing rapidly. The authors show that the algorithm performs well here as this is exactly according to their datagenerating process assumptions (e.g. Figure 2C rightmost), whereas for the other algorithms errors accumulate especially in fast phase changing periods. In my view, this simulated situation conceptually does not capture processes that are relevant to investigating neural oscillations in the brain, which would follow a conceptual model more similar to: multiplicative 1/f contribution * oscillator + additive 1/f contribution, so that presumably corresponds to the presented filtered pink noise case here. Can the authors clarify the motivation for this simulated case?
– The authors speak of broadband oscillations, for instance relating to the murhythm (l. 554). In my view this is a narrowband oscillation with a clear peak judging from the power spectrum, with the broadness of the peak influenced by the presence of strong amplitude envelope modulations. Can authors clarify their criteria for defining an oscillation as broadband, in my view this is nonstandard terminology.
– One limitation of bandpass filtered approaches that the authors mention is nonsinusoidality (l. 48). In my view, this is manifested in differences in waveform shape, e.g. the arcshape of the murhythm in Figure 7D. It does not seem that problems in phase estimation for this type of signals (e.g. nonuniform phasevelocity) is addressed by the method as the estimated state does not capture the arcshape of the rhythm. Can the authors clarify in what way they contribute to resolve this limitation?
4) I would welcome additional details in the text regarding following methodological procedures, which should yield a more selfcontained manuscript:
– While bandpass filtering and Hilberttransform are quite common operations for electrophysiologists, a state space approach may be unfamiliar to the typical reader. Possibly it would be helpful to have clearer paragraphs describing the framework (aka rework l. 98ff), I think the attempt here was to first explain the procedure without formulas, which appear further down, I am not quite sure that this succeeded in providing clarity. Maybe it would be helpful to first make an example with only 1 rhythm, so that dimensionality of the involved entities is not as complex as for tracking N rhythms. Maybe it would also be helpful to extend Figure 1 with a conceptual diagram. How is the number and frequency of oscillators chosen for empirical data?
– E.g. regarding the parameter estimation the reader is referred to the supplementary material of Soulat et al. 2019. Can an intuition about this procedure be given? I find this interesting in the context that a quite large range for expected computation time is given (2100 s, l. 573). Can the authors comment on factors influencing this? Maybe a more detailed description of the procedure will be helpful for understanding this.
– Currently descriptions of the comparison algorithms are lacking, e.g. what is the general procedure and difference between Zrenner and Blackwood algorithm, what are the crucial parameters here? The authors give codeaccessibility as the main reason for selecting those two specific instances. Because recently a lot of different methods (mostly variants involving bandpass filtering) were proposed for estimating the phase in realtime, maybe those 2 selected ones could be better discussed in the context of others in the literature.
– typo, l.55 extant > existant
– the dotted lines for phase estimates may be not optimal visually (e.g. Figure 6B), maybe shading of credible intervals would look nicer?
– Can the authors provide more details on how the power spectra were estimated? The simulation ones look like Welch's method and the data ones like multitaper were used.
– If I understood correctly, the employed twosided ttest is calculated using lots of correlated data points (using all time points in an interval), resulting in reduced degree of freedom. I think the evidence that the SSPE approach works better here is actually strong, but the ttest looks slightly weird because of wrong degrees of freedom, distracting from the result. Maybe the information from table 1 can just be plotted in a distributional way, more like Figure 2C?
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "A State Space Modeling Approach to RealTime Phase Estimation" for further consideration by eLife. Your revised article has been evaluated by Laura Colgin (Senior Editor) and a Reviewing Editor.
The manuscript has been considerably improved and this interesting paper provides inspiration for the field, bringing state space models and phase estimation for neural recordings together. However, there are a few remaining points that are mainly related to presentation results that could be improved. It will not be sent back to the reviewers for final evaluation.
Reviewer #3 (Recommendations for the authors):
I thank the authors for extensive replies to my comments, which have clarified many aspects and made the method description more accessible. I like the new waveform shape figure showing the simultaneous tracking of the base + harmonic frequency. I find the SSPE method interesting from a principal perspective, but am not yet sure in what way the performance of the algorithm is superior to other algorithms in physiologically relevant cases. My remaining comments:
1) Plausibility of simulations R3.Q4
– The state space model simulations (Figure 2A.iv) generating model are not a common data generation model for the rhythms the authors actually investigate (theta and mu). The newly added citations are for the γrhythm (and a quite obscure EEG paper). > I would say this is not a generally accepted data generation model for mu and theta. The fact that the SSPE performs well on data generated from the exact state space model does not seem to entirely relevant for the chosen application, with the simulation lacking stability in cycle periods that is present for mu and theta. Possibly, a note can be added in the Discussion section in the narrowband vs broadband signal section, reflecting on broadband phase would be beneficial.
– The filtered pink noise seems to result in a broader spectral peak (green line in Figure 2B.ii) than is present in the empirical data. As the bandpass filter parameters of the alternative algorithms are possibly adjusted for more narrowband signals, not adjusting the bandwidth before e.g., doing a Hilbert transform may yield an unfair comparison. Can you comment on that?
– My suggested approach using multiplicative 1/f corresponds to the second of the simulations provided by the authors in R3.A4. To clarify this, I attach some python lines at the end to show a simulated signal, if there is interest (parameters to be adjusted). The phase estimation errors in the provided plot in the response seem really low, so maybe the oscillatory SNR is too high for this simulated case (also the amplitude modulation is not really prominent?) I find it confusing that the SSPE method does not yield any benefit here. I am not requesting specifically for this type simulation to be included, but I would like to understand generally the conditions under which SSPE yields better results than the alternative algorithms.
From my understanding, SSPE should yield benefits for phase estimation in intermediate SNR ranges, but at the moment the authors have added more evidence that all algorithms show very similar behavior at high SNR (Figure Supplement I for Figure 7 and 8) and the simulation type provided in Figure 5 has some issues with phase stability as it is filtered pink noise (see comments above) and the modelled SNR range (there is not much change in phase estimation errors which I would have expected to vary with increasing SNR).
R3.Q3:
I can acknowledge the scenario of the two simultaneous rhythms with nearby frequencies which the authors attribute to two different rhythms instead of seeing as amplitude modulated bursts. For this case, maybe the modeling scenario would be more general, if the rhythms had varying phase lag (random for each trial) instead of a fixed phase lag (pi/4), which seems like the case in which the rhythms are most easily distinguishable.
2) The relationship between credible intervals and amplitude measures. R3.Q1.
– The authors argue for a dissociation between credible intervals and amplitude values, with evidence provided by e.g., Figure 7. Given the nature of the data, with oscillations occurring in bursts and not as a continuous oscillation, I am wary of phase estimates for time points where there is no oscillation (e.g., for the theta data, oscillation probability is heavily dependent on animal speed and there will be periods in the data without an oscillation). My main point here is that credible intervals seem to have a correspondence to amplitudes in a range where an oscillation is clearly present (high amplitude) and phase is not defined when an oscillation is absent. Sure, the acausal FIR will return a value for every time point, but it's not necessarily meaningful in a physiological sense. Taking this as the ground truth is problematic for the credible interval and amplitude measures dissociation. For instance, is there a discernible oscillation for the data in Figure 7 for an amplitude criterion > 65th percentile? Is a wide or narrow credible interval around a value that is not defined meaningful?
Figure 7C is a scatter plot with lots of points that are correlated, maybe it would be good to visualize this in a different way reducing the number of points per cycle
possibly by subsampling or a 2D histogram. In the current presentation, the plot basically just shows the range of the data, amplifying lowoccurrence data points.
This suggests that, at low SNR, an arbitrary amplitude threshold (e.g., amplitude > 1) would eliminate all data at SNR=1, and preserve all data at SNR=2.5.
Typically (e.g., in the Zrenner studies), amplitude is defined as a relative criterion, e.g., with percentile threshold or mean amplitude + 2.5 SD using some training interval, so I am not sure about this argument.
– We find that fixed amplitude thresholds produce wide ranges of credible intervals (Figure 7C.ii). This is especially apparent at lower amplitude thresholds; choosing a threshold of 65% results in credible intervals that range from 0.042.29 (minimum) to 174.75 (maximum) degrees.
Can authors list the 95% range here, rather than min/max?
– Figure 5: Can a spectrogram can also be plotted in this figure?
I would expect the height of the spectral peak over the 1/fcontribution to vary here.
3) Code availability
I wanted to check the used bandpass filter parameters (since it's not described for the Blackwood algorithm) and noticed that the linked repository does not feature the code to recreate the plots in the manuscript, possibly it would be helpful for other researchers to provide that, also for comparing to other algorithms in the future.
https://doi.org/10.7554/eLife.68803.sa1Author response
Essential revisions:
1) Consider the impact of noise in the phase estimation, include specifications of the Kalman filter and its robustness, and consider the performance of the estimated phase relative to other methods.
We address these suggestions in detail below.
2) Demonstrate the advantage of using this method for exploring the link between rhythms and behaviour and for phaselocked stimulation, and separate applications that require real time phase estimation from those that can be performed offline using acausal methods.
We address these suggestions in our response in R2.A2.
3) Provide some help to readers from nontechnical background by linking the method back to applications and invivo data more.
Please see our responses in R3.A7 for a detailed response to this suggestion.
4) Expand details regarding computing credible intervals, slow drift considerations, and physiological motivation of the simulations.
Please see our responses in R2.A3, R3.A1, R3.A2, R3.A3 and R3.A4 where we address these points.
Further details are provided below in the 'Recommendations for the authors'.
Reviewer #1 (Recommendations for the authors):
I believe the proposed method in this paper is very novel and timely for the neuroscience community. The paper is well organized and written very clearly.
Thank you for your kind comments, they helped us further describe how performance links to signaltonoise ratio and show the performance of SSPE relative to other methods.
I have some comments that might improve the quality of this paper.
1. Impact of Noise in the phase estimation. Although there is a section in the Discussion related to Noise in Signal, the authors did not address the impact of Signal to Noise Ratio (SNR) on the performance of their algorithms. I think this will add value to your algorithm, specifically because the Kalman Filter can be much more reliable than other methods in this context.
R1.A1. As recommended by the Reviewer, we have updated the manuscript to include a new simulation study characterizing the impact of signal to noise on real time phase estimation. We show that the SSPE method outperforms the other real time methods, and that the credible intervals (CI) vary although signal amplitude is fixed (in response to R3.Q1). We have updated the Results and Methods as follows:
Results, after Section Phase estimation following phase reset
“Phase estimation with different signaltonoise ratios: To demonstrate the impact of noise on phase estimation we simulate signals with different signaltonoise ratios. To generate a signal we filter pink noise into the band of interest (FIR filter, 48 Hz, “Filtered Pink Noise”, see Methods), and then add additional (unfiltered) pink noise. The ratio of the standard deviations of the signal and noise defines the signaltonoise ratio; we modulate the standard deviation of the signal to achieve the desired SNR. Finally, we normalize the full observation so that amplitudes across simulations are consistent. In these simulations the SNR varies from 1 (signal and noise have the same total power) to 10 (signal has 100 times the total power of the noise); see example traces of the observation under each SNR in Figure 5A.i – iv. We repeat each simulation 1000 times with different noise instantiations for each signal to noise ratio.
We compute two metrics to characterize the impact of signaltonoise on phase estimation. First, we consider the credible intervals estimated by the SSPE, here measured as the average credible interval width throughout a simulation. As expected, the credible intervals track with the SNR (Figure 5.B.ii); signals with higher SNR correspond to smaller credible intervals. Computing the average amplitude (as derived from the acausal FIR, see Methods: “NonCausal Phase Estimation”) throughout each simulation, we find no overlap between the distributions of amplitudes at SNR=1 and SNR=2.5 (5B.ii), while the distributions of credible intervals overlap. This suggests that, at low SNR, an arbitrary amplitude threshold (e.g., amplitude > 1) would eliminate all data at SNR=1, and preserve all data at SNR=2.5. Alternatively, using the credible interval width allows identification of transient periods of confidence (e.g., credible interval width < 40 degrees) at both SNR values. We also note that the CI width – in units of degrees – provides a more intuitive threshold than the amplitude, whose intrinsic meaning at any particular value is only understood relatively. These results suggest that, in this simulation study, amplitude serves as a less accurate criterion for confidence in phase compared to the CI width.
Second, we consider the error in phase estimates using four different phase estimation methods (Figure 5C). We find that the error for the SSPE method reduces more quickly than for the other methods, and at a SNR of 1, the SSPE method is comparable to or better than the contemporary methods of realtime phase estimation. At SNR = 10, the error in the SSPE phase estimation is comparable to the acausal FIR, while the alternative realtime algorithms continue to perform with error similar to an SNR of 1. These results suggest that, for the two existing realtime phase estimation methods, broadband rhythms are difficult to track even with high SNR, consistent with Figure 2C.”
Methods, Simulated data, following the subsection Phase Reset;
Signal to Noise Ratio
“To investigate the impact of signaltonoise ratio (SNR) on phase estimation, we again simulate a filtered pink noise signal (centered at 6 Hz) with added pink noise (the third type of rhythm simulated in Methods: Narrow Band to Broad Band Rhythms). To modulate the SNR, we first compute,
where $\sigma}_{\text{signal}$ ($\sigma}_{\text{noise}$) is the standard deviation of the filtered pink noise signal (added pink noise). We then divide the filtered pink noise signal by $\text{SN}{R}_{0}$, multiply the filtered pink noise signal by the chosen SNR (1, 2.5, 5, 7.5, 10), and add this scaled signal to the pink noise.”
2. Specifications of the Kalman Filter and its robustness. specifications of the Kalman filter (KF) should be described in the Method Section. For example, the choice of state and observation covariance matrices (Q and R) should be clarified. In this regard, the authors can discuss the robustness of the KF in phase estimation.
R1.A2. As recommended, we have updated the Methods to clarify the choice of the state and observation covariance matrices (Q and R), as well as the starting values for the state and state covariance as follow:
Methods, following equation (5):
“We note that the state (Q) and observation (R) covariance matrices are estimated from the data, as described below.”
Methods, Section Realtime Phase Estimation:
”We now define the stages of Kalman filtering. The state is initialized with zeros and the state covariance to a diagonal matrix with 0.001 along the diagonal. To predict the future state and state estimation error …”
We also now discuss the robustness of the Kalman filter in phase estimation:
Methods, Section Realtime Phase Estimation:
“ … With these model parameters estimated and fixed, we then apply a Kalman filter to predict and update the state estimates, and estimate the phase and amplitude for each oscillator, each representing a different rhythm, for every sample (Figure 1). We note that the Kalman filter is an optimal filter, i.e., in the presence of Gaussian noise, the Kalman filter achieves minimum error and is an unbiased estimate of the mean phase. However, if the noise is not Gaussian, the Kalman filter remains the optimal linear filter for phase estimation under this model. Given the stochastic frequency modulation that is possible under the model, the Kalman filter is robust to small shifts in the central frequency of the rhythm.”
3. For invivo data, how the performance of the estimated phase might be compared to other methods (used in the simulation study)? It can be helpful to visually compare the estimated phase calculated by SSPE and other methods which were used in the simulation study.
R1.A3. To address this, we have updated the Results to include the estimated phase computed with the other realtime phase estimation methods. We show that phase estimates converge when the rhythm has a high signaltonoise ratio (SNR), but diverge at low SNR for both the theta and the mu rhythm. Additionally, we show error estimates for both the LFP theta and the EEG mu rhythm data, demonstrating that the SSPE performs comparably to the acausal FIR on the data. We note that the error estimate critically depends on the generative model assumed for the in vivo signal (the SSPE method and existing methods assume different underlying models for the rhythm). We have updated the Results in the revised manuscript as follows,
Results, Section Example in vivo application: rodent LFP:
…”However, comparing the estimated phase to the true phase (see Methods), we find that the phase error tends to remain consistent in time (blue curve in Figure 6A; circular standard deviation 95 percent interval = [33.11, 75.39]; for error of phase estimates computed using alternative methods, see Figure Supplement 6I). “…
“To illustrate this, we show the phase and credible interval estimates of the SSPE method for an example 3 s interval of the LFP data (Figure 7; for comparison to phase estimates computed using the existing methods, see Figure Supplement 7I). We note that, …”
Reviewer #2 (Recommendations for the authors):
I think the manuscript can be improved by
1) clearly demonstrating the advantage of using this method for exploring the link between rhythms and behaviour and for phaselocked stimulation. The authors should separate applications that require real time phase estimation from those that can be performed offline using acausal methods.
R2.A2. To address this question, we again utilize the dataset shared by Zrenner and colleagues in Zrenner et al. 2020. There, the authors stimulate using transcranial magnetic stimulation (TMS) at random phases of the mu rhythm and use a posthoc causal analysis to examine whether there is modulation of the motor evoked potentials’ (MEP) amplitudes as a function of the phase of the mu rhythm. They show that using their realtime phase estimation method, MEP amplitude is selectively higher closer to the trough of the mu rhythm. We test whether the SSPE method can replicate these findings. To do so, we apply the SSPE method, using parameters estimated on a different subject, to the prestimulus EEG data to estimate the phase at the moment of TMS stimulation. We find that phase as estimated from the SSPE method is a strong predictor of the MEP amplitude (when fitting a circular linear regression model), with the peak close to the mu rhythm trough as predicted.
We add a new figure and results to the Results section of the manuscript, as follows:
Results: Prestimulus MuRhythm Phase Predicts MEP amplitude
“To illustrate the utility of realtime phase estimation, we again analyze an example EEG recording from (Zrenner et al., 2020) in which the authors stimulate at random intervals using transcranial magnetic stimulation (TMS) a subject using to elicit a motor evoked potential (MEP) from the right abductor pollicis brevis muscle. We analyze the MEP amplitude to assess the influence of prestimulus mu rhythm phase on the response (replicating the analysis performed by Zrenner et al. 2020). We estimate mu rhythm phase using the SSPE method with model parameters derived from a different subject. We find that the prestimulus phase of the mu rhythm, at the instant immediately prior to stimulation, predicts (using circular linear regression) 13 percent of the logtransformed MEP amplitude (Figure 9A). Furthermore, thresholding the trials using credible interval width, we find improved prediction of the MEP amplitude. When using trials with credible interval width less than 50 degrees, the mu rhythm phase explains 30% of the variance in the MEP amplitude (Figure 9B), and with credible interval width less than 25 degrees, 48% of the variance in the MEP amplitude (Figure 9C). These results show that modulation of the MEP amplitude depends on the prestimulus phase of the mu rhythm, and suggest that accurate estimation of the mu rhythm (e.g., using the SSPE method with narrow credible interval width) may enhance the evoked motor response.”
Further, we have ensured that in the manuscript, we clarify whether realtime analysis is required:
Phase Estimation Following Phase Reset: “… Further, accurate tracking of phase resets would provide better evidence for the relevance of phase resetting in generating event related potentials (Sauseng et al., 2007) – an analysis that could also be performed acausally.”
2) Building on my previous point:
– In figure 2, proposed method clearly outperforms others for signals generated using the "state space model". How representative is this for invivo recordings?
R2.A3. To clarify this point, we have updated the manuscript to motivate the state space model as representative of features in in vivo recordings as follows:
Methods: Narrow Band to Broadband Rhythms: “For the final rhythm we simulate a signal $y}_{t$ from the SSPE model (see Methods: State Space Model Framework) which is nonsinusoidal and broadband. … The state in the state space model is an autoregressive model of order 2 (West 1997a) (equivalently, a damped harmonic oscillator driven by noise), which has been proposed to capture features of γ rhythms in visual cortex and rhythms in EEG (Franaszczuk and Blinowska 1985, Xing et al., 2012, Burns et al. 2010, Spyropoulos et al. 2020). The addition of observation noise under the state space model can be analogized to measurement noise. Thus, the state space model provides a simple possible model for brain rhythms observed in electrophysiological data.”…
New references:
West, M. (1997a) Time series decomposition. Biometrika, 84, 489494
Xing, D., Shen, Y., Burns, S., Yeh, C. I., Shapley, R., and Li, W. (2012). Stochastic generation of γband activity in primary visual cortex of awake and anesthetized monkeys. Journal of Neuroscience, 32(40), 1387313880a.
Burns, S. P., Xing, D., Shelley, M. J., and Shapley, R. M. (2010). Searching for autocoherence in the cortical network with a timefrequency analysis of the local field potential. Journal of Neuroscience, 30(11), 40334047.
Franaszczuk, P. J., and Blinowska, K. J. (1985). Linear model of brain electrical activity—EEG as a superposition of damped oscillatory modes. Biological cybernetics, 53(1), 1925.
Spyropoulos, G., Dowdall, J. R., Schölvinck, M. L., Bosman, C. A., Lima, B., Peter, A., … and Fries, P. (2020). Spontaneous variability in γ dynamics described by a linear harmonic oscillator driven by noise. bioRxiv, 793729.
– In figure 4, it is demonstrated that other methods cannot track phasereset, however, is this offset sustained or is it a transient slip in phase estimation for one cycle. In other words, if we zoom out in time, do the other methods converge back to true phase? And if so, when would it be critical to track phasereset for one cycle?
R2.A4. Indeed, it is the case that the methods all return to accurately track the rhythm after some time. This simulation suggests that, when stimuli appear for a human or nonhuman animal during normal behavior (when phase resets might occur), the SSPE method continues to perform well. In fact, recent work attempting to perform phase stimulation based on the theta rhythm has been forced to avoid segments with phase resets to reduce error in realtime phase estimation (Gordon et al. 2021). While bandpass filter approaches struggle with phase estimation in the presence of phase resets as a function of the window size used, the SSPE method (by virtue of the Kalman Filter) continues to perform optimally.
To address this important question, we include additional analysis and explanation in the revised manuscript as follows:
Results: Phase estimation following phase reset: “…Phase resets are common and well documented in electrophysiological recordings (Fiebelkorn et al., 2011; Makeig et al., 2004). Recent work delivering electrical stimulation based on phase of the theta rhythm has been forced to avoid segments with phase resets to reduce error in realtime phase estimation (Gordon et al., 2021). Further, accurate tracking of phase resets …
…“These results demonstrate that the existing causal and acausal methods implemented here are unable to estimate the phase accurately immediately after a phase reset. To determine how quickly each method returned to accurate phase estimation, we first estimated the prestimulus average error for each method, prior to the first phase reset, over a 500 ms interval. We then determined the time between each subsequent phase reset and the moment when the error reduced to 1.5 times the prestimulus error. Note that the level of prestimulus error differs across methods, so, for example, we determine when the SSPE converges to a lower error than the Blackwood method (see top row of Table 1). We find that each method eventually converged to accurate phase estimation after a phase reset, with the SSPE method converging most rapidly (Table 1). The differences between methods likely result from the windowing used by each method when estimating the phase; smaller windows will lead to less impact of the phase reset on the error (as seen with the Blackwood approach converging more quickly than the Zrenner method), with the SSPE needing a window of a single sample.”
– The performance of SSPE should be demonstrated with respect to other methods for invivo data (figures 5 and 7) so that the reader can evaluate the utility of this method.
R2.A5. As recommended, we now include this comparison in the revised manuscript. Please see R1.A3 for a detailed response (including new figures) to this comment.
3) The paper is well written but comes across as very technical. Readers from nontechnical backgrounds may struggle with following some of the concepts. This could be improved by linking the method back to applications and invivo data more, and also by reviewing the language used in the methods segment.
R2.A6. To address this comment, we have updated the manuscript to include additional intuitive descriptions (including a new schematic diagram) in the Methods; please see R3.A7 for details. In addition, we now include a new in vivo application; please see R2.A2 for details.
Reviewer #3 (Recommendations for the authors):
1) Credible intervals
An interesting aspect of the manuscript is the possibility of computing credible intervals, making it possible to establish bounds on time points when the phase can be estimated with sufficient accuracy. In previous approaches, this is usually done via an amplitude threshold, with the assumption that a high amplitude will also result in high phase estimation accuracy. In Figure 6 the authors therefore study the relationship between amplitude and credible intervals. It is tough to interpret this evidence; the argument is that for an amplitude threshold using the 65th percentile the range of credible intervals is from the minimum to maximum possible. But looking at the major bulk of probability mass, it's centered around intervals with an acceptable looking range? Could more details regarding computation of confidence intervals be provided, e.g. are they symmetric? Figure 6C is a scatter plot with lots of points that are correlated, maybe it would be good to visualize this in a different way reducing the number of points per cycle? Can the credible intervals be related to error in degrees? They look pretty negligible in the left part of the figure (given that the realtime procedure also introduces a lag, Figure 8B),is the plotted xrange the relevant one? How does SNR/amplitude relate to phase accuracy to credible intervals? If the spread of credible intervals can be largely explained by amplitude (aka high correlation between credible interval width and amplitude), then maybe this is not such a big factor?
R3.A1. We address these important points in the updated manuscript as follows:
Updates to original Figure 6 (now Figure 7)
To help improve the interpretation of original Figure 6 (now Figure 7), we have updated the xaxis to show values in degrees. Additionally, we have edited the Figure caption to reflect the autocorrelation of points across time:
Figure 7 … (C.i): Phase credible intervals versus theta rhythm amplitude for each time point (gray dot). On the xaxis we plot an example cycle of a rhythm with credible intervals, see xaxis of (C.ii) for numerical values of credible intervals. Note that the phase and amplitude estimates are autocorrelated in time so that each gray dot is not independent of every other dot.
Could more details regarding computation of confidence intervals be provided, e.g. are they symmetric?
To address this, we have updated the Methods as follows:
Methods, after Equation (10): “The SSPE method also tracks the state estimation error $({P}_{t}^{t})$, which allows estimation of confidence in the mean phase estimates. To do so, we sample 10,000 samples from the posterior distribution of each state j, and from this distribution estimate the width of the 95% credible interval for the phase. The credible intervals are derived from the posterior distribution of the state, which is Gaussian; however, following their transformation into angles, the intervals are not constrained to be symmetric. Note that, alternatively, one could attempt to define a closed form for the phase distribution to estimate confidence bounds …”
Can the credible intervals be related to error in degrees?
Yes, the credible intervals are expressed in units of radians. We have updated Figure 6 to units of degrees.
They look pretty negligible in the left part of the figure (given that the realtime procedure also introduces a lag, Figure 8B),is the plotted xrange the relevant one?
The realtime procedure in C++ introduces a small lag due to buffering implemented in OpenEphys. This lag is not present if the user does a nonbuffered analysis as provided in the MATLAB implementation.
How does SNR/amplitude relate to phase accuracy to credible intervals?
To address this question, we now include in the revised manuscript a new simulation comparing SNR, credible interval width, and phase accuracy. Please see R1.A1 for details.
If the spread of credible intervals can be largely explained by amplitude (aka high correlation between credible interval width and amplitude), then maybe this is not such a big factor?
To address this question, we have updated the Results as follows:
Results: Confidence in the Phase Estimate: “… Therefore, a fixed amplitude threshold does not guarantee an accurate phase estimate. Instead, we propose that inspection of the credible intervals better guides certainty in the phase estimate. These credible intervals are directly derived from the modeling procedure, as opposed to surrogate measures of confidence in the phase estimate, such as amplitude. In addition, credible intervals are directly interpretable unlike amplitude thresholds, which depend on the measurement scale (e.g., a credible interval threshold of 10 degrees versus an amplitude threshold of 10 millivolts). We conclude that amplitude is linked to confidence in the phase, however, a fixed amplitude threshold is not directly interpretable and encompasses a wide range of credible interval values.”
2) Slow drifts
A practical reason why existent algorithms rely on buffering is the frequent existence of slow drifts in the data, for which detrending needs to be applied. For the EEG data, the following methodological detail of the current approach is given "Before realtime phase estimation, we highpass filter the data at 0.5 Hz to remove slow drifts." This procedure is not possible when estimating the phase in realtime. The authors claim the benefit of their approach is the absence of needing buffers, but in the light of this I am not sure about whether this is the case. In the discussion, the authors discuss the possibility of modeling the slow drifts via an additional rhythm in the state space, but state that long segments of data are needed for this. In my view, this is a crucial aspect of any algorithm targeting realtime phase and more information and analysis resolving slow drifts would greatly improve the actual applicability of the present approach.
R3.A2. To address this important issue, we have updated the procedure to apply a realtime method that avoids acausal filtering. We have updated the Methods to describe this new approach:
Methods: in vivo Data: “…We analyze 250 s of data from a single participant for demonstration purposes. For these data, before each realtime phase estimation step, we first apply a moving average causal filter (using 3 s of past data for each new data point) to remove slow drifts. We track three oscillators in the SSPE method, …”
Discussion: ”Alternatively, the SSPE model could be updated to include a model of the slow drift dynamics (perhaps as a piecewise linear term), or a moving average mean trend could be causally removed as we demonstrated with the in vivo EEG analysis…”
Using this new realtime approach, we find results consistent with the original manuscript. We include the updated Figure 8 (previously Figure 7) and revised caption here:
3) Data generation process and simulations
I find some of the simulations not sufficiently motivated from a physiological perspective.
– Two simultaneous rhythms with nearby frequencies: In this scenario the authors simulate the signal as a summation of two sinusoids with a frequency offset with a fixed phase shift. This will result in a signal with amplitude modulation with the envelope given by the difference of the frequencies, resulting in the appearance of oscillatory bursts. I would argue that for physiological signals that present like this, it is not of interest to track a single frequency in this signal, as e.g. the trough phase of this compound signal corresponds to a specific physiological state. I would also suggest showing the amplitude for this case in Figure 3. Possibly the simulation could be modified to be better physiologically motivated, as they are certainly cases where there are neighboring rhythms which present on the same electrode through volume conduction, but it is unclear (but would be of interest) to me if the proposed method could be helpful in this case.
R3.A3. We appreciate the Reviewer’s concerns regarding this simulation. Indeed, the Reviewer’s interpretation is correct: the summed sinusoids can be interpreted as a single neural population oscillating at the average of the two frequencies, with amplitude modulation. With that interpretation, a different analysis approach – perhaps focused on crossfrequency coupling – might be preferred, However, we note that multiple interpretations exist for such a situation, i.e., we believe the inverse solution here is degenerate making our prior expectation important. In many experimental situations, limited spatial resolution prevents distinguishing whether an observed signal represents a single neural population or rhythms from distinct neural populations. In addition, we expect situations in which multiple rhythms exist at nearby frequencies (such as multiple α in EEG over sensorimotor cortex or multiple theta in the hippocampus). Consistent with the latter examples, in this simulation, we assume the prior expectation that the observed signal results from two rhythms.
We have added text to clarify this assumption, and added the observed signal to the revised Figure as suggested:
Methods: Two Simultaneous Rhythms With Nearby Frequencies: “… In many physiologic use cases (e.g., in EEG (VanRullen, 2016), and in LFP (Tort et al., 2010)), we may wish to track one rhythm while ignoring another at a nearby frequency. This may occur because of low spatial resolution (e.g., from mixing of signals due to volume conduction (Nunez and Srinivasan, 2006)) or because there exist overlapping neural populations expressing different rhythms (Kocsis et al., 1999). To simulate this case, we consider two sinusoids. …"
– From my perspective phase is only welldefined if there are clearcut oscillatory bursts. In that sense, the signal given as an example in Figure 2A iv does not have clearcut oscillatory states, with no clear spectral peak and the true phase (blue curve) changing rapidly. The authors show that the algorithm performs well here as this is exactly according to their datagenerating process assumptions (e.g. Figure 2C rightmost), whereas for the other algorithms errors accumulate especially in fast phase changing periods. In my view, this simulated situation conceptually does not capture processes that are relevant to investigating neural oscillations in the brain, which would follow a conceptual model more similar to: multiplicative 1/f contribution * oscillator + additive 1/f contribution, so that presumably corresponds to the presented filtered pink noise case here. Can the authors clarify the motivation for this simulated case?
R3.A4. In the revised manuscript, we now better motivate the statespace model in terms of neural oscillations in the brain. Please see R2.A3 above.
Additionally, we simulated the model suggested by the Reviewer in two ways. In the first, we implemented a model where (1/f) noise was simulated and directly multiplied by a cosine oscillator at 6 Hz, sample by sample. We plot an example realization of this process in Author response image 1:
Sign changes in the pink noise result in discontinuous jumps in the signal. Given this, we attempted a second simulation where the multiplicative (1/f) component was understood to mean that the amplitude had a power spectrum that followed a (1/f) trend, while the phase followed a 6 Hz rotation, with added (1/f) noise. We plot an example realization of this process in Author response image 2:
Testing the error in phase estimation for this signal with added pink noise, we find that the acausal FIR, Blackwood and SSPE methods perform comparably, while the Zrenner approach performs slightly worse than the others (see Author response image 3).
We include this result here for the Reviewer, and if recommended, we will update the manuscript to include it.
– The authors speak of broadband oscillations, for instance relating to the murhythm (l. 554). In my view this is a narrowband oscillation with a clear peak judging from the power spectrum, with the broadness of the peak influenced by the presence of strong amplitude envelope modulations. Can authors clarify their criteria for defining an oscillation as broadband, in my view this is nonstandard terminology.
R3.A5. As suggested by the Reviewer, we now refer to the mu rhythm as narrowband:
Results: Example invivo application: human EEG: “… Like the theta rhythm in the LFP rodent data, the mu rhythm in the human EEG data maintains a consistent spectral profile in time (Figure 8B) and further, it appears to be narrowband (Figure 8C).”
Figure 8, caption: “(C) Average power spectral density (PSD) of (B) over all times in (A). The mean PSD (dark blue) and 95% confidence intervals (light blue) reveal a peak in the mu rhythm near 11 Hz. (D) Example human EEG data (red, solid) with a consistent peak in the mu rhythm (813 Hz). …”
– One limitation of bandpass filtered approaches that the authors mention is nonsinusoidality (l. 48). In my view, this is manifested in differences in waveform shape, e.g. the arcshape of the murhythm in Figure 7D. It does not seem that problems in phase estimation for this type of signals (e.g. nonuniform phasevelocity) is addressed by the method as the estimated state does not capture the arcshape of the rhythm. Can the authors clarify in what way they contribute to resolve this limitation?
R3.A6. We thank the Reviewer for this nice observation. For the human EEG, the SSPE method models both the 10 Hz and the 22 Hz (a potential harmonic) rhythms that together capture the nonsinusoidality of the murhythm. In Figure 7D of the original manuscript, we only showed the 10 Hz component of the model. When we include the summation over the (realvalued portions of the) state estimates of both rhythms, the model estimates capture the arcshape of the mu rhythm (see Figure 8—figure supplement 1). We also show the phase for the two oscillators, demonstrating the capability of the SSPE approach to estimate simultaneously the phase of the 22 Hz harmonic and the 10 Hz rhythm.
We have updated the manuscript to include this result as a Supplementary Figure as follows:
Results: Example invivo application: human EEG: “As with the LFP, we are able to assess certainty in our phase estimates using the credible intervals. In Figure Supplement 8III, we show the state estimate including both the 10 Hz and 22 Hz components of the model and demonstrating that it recreates the nonsinusoidal waveform of the mu rhythm.…”
4) I would welcome additional details in the text regarding following methodological procedures, which should yield a more selfcontained manuscript:
– While bandpass filtering and Hilberttransform are quite common operations for electrophysiologists, a state space approach may be unfamiliar to the typical reader. Possibly it would be helpful to have clearer paragraphs describing the framework (aka rework l. 98ff), I think the attempt here was to first explain the procedure without formulas, which appear further down, I am not quite sure that this succeeded in providing clarity. Maybe it would be helpful to first make an example with only 1 rhythm, so that dimensionality of the involved entities is not as complex as for tracking N rhythms. Maybe it would also be helpful to extend Figure 1 with a conceptual diagram. How is the number and frequency of oscillators chosen for empirical data?
R3.A7. As recommended, we have update Figure 1 to include a conceptual diagram to help improve understanding of the SSPE method:
Further, we have edited the beginning of the Methods section to improve accessibility, and now include an example with only 1 rhythm, as recommended:
Methods: State Space Model Framework:
“To estimate phase in realtime, we utilize a datadriven model that operates based on principles from dynamical systems and Markov models. This type of model (called a statespace model) separates what is observed (called the observation equation) from what we wish to estimate (called the state equation). The state equation attempts to capture the underlying, unobserved dynamics of the system, while the observation equation transforms the state into the observed signal. We generate an optimal prediction for the state using a Kalman filter that compares the initial prediction of the observation with the actual observed value. […] Note that this model for the state is akin to a damped harmonic oscillator driven by noise, a model that has been shown to be relevant and useful for γ rhythms in the visual cortex (Burns et al., 2010; Spyropoulos et al., 2020) and also for EEG (Franaszczuk and Blinowska, 1985).
To model N rhythms, we define …”
How is the number and frequency of oscillators chosen for empirical data?
As we discuss in the Methods section, we chose the number and frequency of oscillators for the empirical data based on a preliminary analysis of the power spectral density. We have reproduced the relevant text below:
In vivo data: To model these data, we fit three oscillators with the SSPE method, with the goal of tracking a target theta rhythm (defined as 4 – 11 Hz). Motivated by an initial spectral analysis of these data, we choose one oscillator to capture low frequency δ band activity (${\omega}_{j}=1$ Hz, based on spectral peak), one to capture the theta band rhythm of interest (${\omega}_{j}=7$ Hz), and a final oscillator to track high frequency activity (${\omega}_{j}=40$ Hz).
– E.g. regarding the parameter estimation the reader is referred to the supplementary material of Soulat et al. 2019. Can an intuition about this procedure be given? I find this interesting in the context that a quite large range for expected computation time is given (2100 s, l. 573). Can the authors comment on factors influencing this? Maybe a more detailed description of the procedure will be helpful for understanding this.
R3.A8. To address these questions, we have added text to the Methods section to help further explain the EM procedure and offer potential reasons for the range of times to convergence:
Methods: Realtime Phase Estimation:
“We perform realtime estimation of phase as follows. First, we use an existing interval of data to acausally fit the parameters of the state space model $a}_{j$, $\omega}_{j$, $Q}_{j$, $\sigma}_{R}^{2$ using an expectationmaximization (EM) algorithm as proposed by (Soulat et al., 2019); for details see the Supplementary Information in (Soulat et al., 2019). Under the EM approach, optimization of parameters follows a two step process. First, initial values for parameters are selected, usually from prior knowledge such as an examination of the power spectrum of an initial data sample. This allows the algorithm to estimate the expectations for the state and the observation. Second, using the state and observation estimates, an analytic solution exists (Shumway and Stoffer, 1982) for the parameters that maximizes the likelihood. We repeat the expectation and maximization procedures until the parameter estimates do not change between two iterations beyond a threshold. The rate of convergence of the EM algorithm (speed of computation) depends on the initial parameter estimates and the signaltonoise ratios of the rhythms present in the signal.”
– Currently descriptions of the comparison algorithms are lacking, e.g. what is the general procedure and difference between Zrenner and Blackwood algorithm, what are the crucial parameters here? The authors give codeaccessibility as the main reason for selecting those two specific instances. Because recently a lot of different methods (mostly variants involving bandpass filtering) were proposed for estimating the phase in realtime, maybe those 2 selected ones could be better discussed in the context of others in the literature.
R3.A9. Thank you for this suggestion. We have updated the manuscript to place the two algorithms better in the context of others in the literature as follows:
Methods: ARbased forecasting to Estimate Phase: “While current techniques depend upon filtering to isolate a narrow band rhythm from which to estimate phase, different methods exist to forecast the phase. One category of methods forecasts the phase as a linear extrapolation from the current moment while assuming a singular central frequency, which is estimated from the data (Rivero and Ditterich, 2021; Rutishauser et al., 2013). Other methods estimate phase by using an autoregressive (AR) model. After estimating the AR parameters (either asynchronously or synchronously), these methods forecast data to limit filter edge effects, then utilize a filter and Hilbert transform to estimate phase (Blackwood et al., 2018 and Zrenner et al., 2020). We compare the SSPE to this latter set of techniques as these are more flexible to variability in rhythms than methods that assume a constant central frequency when estimating phase. We implement the version of the Zrenner et al. (2020) algorithm accessible at …”
We have also added information about the critical hyperparameters of the Zrenner and Blackwood approaches:
Methods: ARbased forecasting to Estimate Phase: “ … We implement the version of the Zrenner et al. (2020) algorithm accessible at https://github.com/bnplab/phastimate and term this method as Zrenner. Please refer to Zrenner et al. (2020) for algorithm details. The Zrenner algorithm uses an FIR filter (applied forward and backward) and the Hilbert transform as the core phase estimation technique. The AR model used to forecast the data and allow realtime phase estimation (by expanding the window for the FIR filter) is fitted to the 30th order and efficiently updated on every new sample. We set the algorithm to track a 6 Hz rhythm by increasing the window size for the FIR filter to 750 ms, the filter order to 192, and the frequency band to 4 to 8 Hz. We utilize the implementation of the Blackwood et al. (2018) algorithm as provided by the authors. In Blackwood et al. (2018), the authors develop a forward bandpass filter – ARbased (5th order AR model whose parameters are fit every second) forecasting and then the Hilbert transform – to compute the phase. Here, we instead apply a Hilbert transformer (of order 18) that estimates the phase by directly computing the Hilbert transform while filtering the data as implemented in (Blackwood, 2019). We do so because this implementation is more computationally efficient, functionally equivalent to the original method and the code is made widely available. We verified in sample data that application of the original and modified Blackwood methods produced consistent phase estimates. We refer to this method as Blackwood.”
– typo, l.55 extant > existant
R3.A10. Thank you, we have edited the manuscript, replacing “extant” with “existing”.
Introduction L55: “Many existing approaches depend on buffered processing (e.g., Fourier or related transforms), …”
– the dotted lines for phase estimates may be not optimal visually (e.g. Figure 6B), maybe shading of credible intervals would look nicer?
R3.A11. Thank you for the suggestion. We have updated the figures to use shading in the revised manuscript. We show the relevant components of the revised figures here:
– Can the authors provide more details on how the power spectra were estimated? The simulation ones look like Welch's method and the data ones like multitaper were used.
R3.A12. As recommended, we have edited the caption of Figure 2 to include these details:
“… (B) For each scenario, example spectra of (B.i) the signal, and (B.ii) the observation (i.e., signal plus noise). Spectra were estimated for 10 s segments using the function “pmtm” in MATLAB, to compute a multitaper estimate with frequency resolution 1 Hz and 9 tapers.”
For the in vivo data, the Figure captions include details of the power spectrum calculation, reproduced here:
Figure 6 “… (B) Corresponding spectrogram of the rodent LFP data (multitaper method with window size 10 s, window overlap 9 s, frequency resolution 2 Hz, and 19 tapers). (C) Average power spectral density (PSD) of (B) over all times with theta band oscillator identified in (A). The mean PSD (dark blue) and 95% confidence intervals (light blue) reveal a broadband peak in the theta band.
Figure 8 “… (B) Corresponding spectrogram for the EEG data (multitaper method with window size 10 s, window overlap 9s, frequency resolution 1 Hz, and 19 tapers). (C) Average power spectral density (PSD) of (B) over all times in (A). The mean PSD (dark blue) and 95% confidence intervals (light blue) reveal a broadband peak in the murhythm near 11 Hz. …”
– If I understood correctly, the employed twosided ttest is calculated using lots of correlated data points (using all time points in an interval), resulting in reduced degree of freedom. I think the evidence that the SSPE approach works better here is actually strong, but the ttest looks slightly weird because of wrong degrees of freedom, distracting from the result. Maybe the information from table 1 can just be plotted in a distributional way, more like Figure 2C?
R3.A13. The statistical test compares the error estimates across all simulation iterations (N=1000), wherein each iteration has an independent source of noise. Error within each iteration is calculated using all time points within one cycle of the phase reset. We note that the error distributions for each method are approximately normally distributed (see Author response image 4 and 5), and we find that plotting these error estimates as distributions is rather uninformative (Author response image 4) given the tight bounds of the error distributions, and the large differences in means. We therefore think a table is a better approach for conveying this result.
We also examined the error to confirm approximate Gaussianity (Author response image 5):
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
Reviewer #3 (Recommendations for the authors):
I thank the authors for extensive replies to my comments, which have clarified many aspects and made the method description more accessible. I like the new waveform shape figure showing the simultaneous tracking of the base + harmonic frequency. I find the SSPE method interesting from a principal perspective, but am not yet sure in what way the performance of the algorithm is superior to other algorithms in physiologically relevant cases. My remaining comments:
1) Plausibility of simulations R3.Q4
– The state space model simulations (Figure 2A.iv) generating model are not a common data generation model for the rhythms the authors actually investigate (theta and mu). The newly added citations are for the γrhythm (and a quite obscure EEG paper). > I would say this is not a generally accepted data generation model for mu and theta. The fact that the SSPE performs well on data generated from the exact state space model does not seem to entirely relevant for the chosen application, with the simulation lacking stability in cycle periods that is present for mu and theta. Possibly, a note can be added in the Discussion section in the narrowband vs broadband signal section, reflecting on broadband phase would be beneficial.
We agree with the Reviewer that (i) the state space model has not been shown to reproduce theta and mu rhythms, and (ii) it is not surprising that the SSPE method performs well on data generated from the same state space model. We propose that the results shown here provide initial evidence that the state space model can produce theta and mu rhythms, consistent with the in vivo data. We also propose that modeling rhythms as damped harmonic oscillators driven by noise is both simple and physically reasonable. Finally, we note while the parameters used for the state space model instantiation in Figure 2 do not produce phase stability, the model can sustain such phase stability, as demonstrated by the intervals of transient phase stability the Reviewer identifies for the mu and theta rhythms (Figure 7,8).
As suggested, we have updated the Discussion to note that the phase cycles of the state space model and filtered pink noise are more unstable:
“Moreover, when the broadband rhythm has nonnegligible power across all frequencies (in which case the phase cycle is more unstable), the standard nonrealtime approach (acausal FIR and Hilbert transform) tracks the phase with limited accuracy.”
– The filtered pink noise seems to result in a broader spectral peak (green line in Figure 2B.ii) than is present in the empirical data. As the bandpass filter parameters of the alternative algorithms are possibly adjusted for more narrowband signals, not adjusting the bandwidth before e.g., doing a Hilbert transform may yield an unfair comparison. Can you comment on that?
To address this, we now emphasize in the Discussion the importance of filter bandwidth parameters,
“Narrowband vs Broadband Signals
While brain rhythms tend to be broadband (Buzsaki, 2004; Roopun et al., 2008), the most common methods in neuroscience to track a rhythm’s phase require appropriate choices of filter parameters (e.g., selection of bandwidth) and generally expect a narrowband oscillation. …”
We also note that the SNR in Figure 2B.ii is high, increasing the apparent bandwidth. Increasing the background pink noise (and making the simulated signal more consistent with the in vivo data), would decrease the apparent bandwidth in Figure 2B.ii (see Figure 5 —figure supplement 1).
Finally, we now make public all the code for the simulations and for plotting the results of the simulations so parameters across approaches can be compared. It is available here: https://github.com/EdenKramerLab/SSPEpaper
– My suggested approach using multiplicative 1/f corresponds to the second of the simulations provided by the authors in R3.A4. To clarify this, I attach some python lines at the end to show a simulated signal, if there is interest (parameters to be adjusted). The phase estimation errors in the provided plot in the response seem really low, so maybe the oscillatory SNR is too high for this simulated case (also the amplitude modulation is not really prominent?) I find it confusing that the SSPE method does not yield any benefit here. I am not requesting specifically for this type simulation to be included, but I would like to understand generally the conditions under which SSPE yields better results than the alternative algorithms.
Thank you again for recommending this interesting model. The SSPE method does not yield any benefit here because the simulated model of phase is a sinusoid. Thus, methods that assume a narrowband oscillation do better, even when the rhythm waxes and wanes via amplitude modulation; after an accurate estimate of phase when the SNR is high, the true and estimated phases progress linearly, even when the SNR is low (i.e., no rhythm is visible). This is most obvious for the Blackwood and acausal FIR approaches, and less so for the Zrenner method which forecasts more data prior to phase estimation.
We provide here results from the new simulation, where we set the SNR = 0.5 (code also available at https://github.com/EdenKramerLab/SSPEpaper/blob/main/simulations/sim_multOneOverF_amp.m):
We conclude – based on the existing simulations in Figure 2, and the new simulation suggested here – that these examples illustrate a more complicated version of the “Sines in Pink Noise” simulation shown in Figure 2; when the true signal is a sinusoid, those estimators that presume a narrow band signal (i.e., the filter based approaches) perform best.
From my understanding, SSPE should yield benefits for phase estimation in intermediate SNR ranges, but at the moment the authors have added more evidence that all algorithms show very similar behavior at high SNR (Figure Supplement I for Figure 7 and 8) and the simulation type provided in Figure 5 has some issues with phase stability as it is filtered pink noise (see comments above) and the modelled SNR range (there is not much change in phase estimation errors which I would have expected to vary with increasing SNR).
We propose that the primary benefits of the SSPE method are due to the state model, more so than the SNR. If the state model is inaccurate, then a phase estimation method will often perform poorly across a broad range of SNR. This is illustrated in Figure 5; estimates from the Blackwood and Zrenner methods do not improve with increasing SNR because the narrowband state model is inaccurate (for this simulated broadband rhythm). For the in vivo data examples (Figure 7,8 Supp I), the phase estimates agree when the rhythm of interest is clearly present. At these times, different state models (sinusoids or noise driven oscillators) produce qualitatively consistent phase estimates.
R3.Q3:
I can acknowledge the scenario of the two simultaneous rhythms with nearby frequencies which the authors attribute to two different rhythms instead of seeing as amplitude modulated bursts. For this case, maybe the modeling scenario would be more general, if the rhythms had varying phase lag (random for each trial) instead of a fixed phase lag (pi/4), which seems like the case in which the rhythms are most easily distinguishable.
The Reviewer raises a good point. To address this, we implemented this recommendation and repeated the simulations in Figure 3 such that the phase lag between the two rhythms varied randomly (uniformly from 0 to 2*pi) for each trial. We find qualitatively consistent results, Author response image 8:
We now state this results in the revised manuscript as follows,
Following Equation (14): “… where [0.2,2] and [1,11]. We note that replacing the constant phase shift ($\pi /4$) with a randomly selected (uniform $[\pi ,\text{}\pi ]$) phase shift produces qualitatively consistent results. The confounding rhythm frequency assumes a range of integer values from 1 to 11 Hz. …”
2) The relationship between credible intervals and amplitude measures. R3.Q1.
– The authors argue for a dissociation between credible intervals and amplitude values, with evidence provided by e.g., Figure 7. Given the nature of the data, with oscillations occurring in bursts and not as a continuous oscillation, I am wary of phase estimates for time points where there is no oscillation (e.g., for the theta data, oscillation probability is heavily dependent on animal speed and there will be periods in the data without an oscillation). My main point here is that credible intervals seem to have a correspondence to amplitudes in a range where an oscillation is clearly present (high amplitude) and phase is not defined when an oscillation is absent. Sure, the acausal FIR will return a value for every time point, but it's not necessarily meaningful in a physiological sense. Taking this as the ground truth is problematic for the credible interval and amplitude measures dissociation. For instance, is there a discernible oscillation for the data in Figure 7 for an amplitude criterion > 65th percentile? Is a wide or narrow credible interval around a value that is not defined meaningful?
To address this, we provide here additional examples of the LFP data shown in Figure 7. These examples illustrate cases in which the amplitude is low (below the 65th percentile), but the credible interval is narrow. At these times, visual inspection suggests a theta rhythm exists, despite the low amplitude. These examples illustrate how an amplitude threshold may fail to capture transient bursts of theta. Conversely, when a theta burst exceeds an amplitude threshold, existing methods may be overly confident in the stability of the phase cycle; these examples show that, even when the amplitude is high, the credible intervals can be large.
Figure 7C is a scatter plot with lots of points that are correlated, maybe it would be good to visualize this in a different way reducing the number of points per cycle
possibly by subsampling or a 2D histogram. In the current presentation, the plot basically just shows the range of the data, amplifying lowoccurrence data points.
As recommended by the Reviewer, we have updated Figure 7C to include a 2D histogram replacing the scatterplot:
This suggests that, at low SNR, an arbitrary amplitude threshold (e.g., amplitude > 1) would eliminate all data at SNR=1, and preserve all data at SNR=2.5.
Typically (e.g., in the Zrenner studies), amplitude is defined as a relative criterion, e.g., with percentile threshold or mean amplitude + 2.5 SD using some training interval, so I am not sure about this argument.
The Reviewer identifies a nice solution to avoid the choice of an arbitrary amplitude threshold eliminating or including all data. However, we propose that this criterion still lacks the direct interpretability of a credible interval threshold (e.g., 30 degrees, which is easily interpreted and may include data across amplitudes).
– We find that fixed amplitude thresholds produce wide ranges of credible intervals (Figure 7C.ii). This is especially apparent at lower amplitude thresholds; choosing a threshold of 65% results in credible intervals that range from 0.042.29 (minimum) to 174.75 (maximum) degrees.
Can authors list the 95% range here, rather than min/max?
Thank you for the suggestion, we have updated the test to reflect the 95% range:
“This is especially apparent at lower amplitude thresholds; choosing a threshold of 65% results in credible intervals that range from 3.8 (2.5th percentile) to 96.8 (97.5th percentile) degrees.”
– Figure 5: Can a spectrogram can also be plotted in this figure?
I would expect the height of the spectral peak over the 1/fcontribution to vary here.
To address this, we now include a supplement to Figure 5 that includes the spectra for each simulated SNR.
We have updated the manuscript as follows:
“… In these simulations the SNR varies from 1 (signal and noise have the same total power) to 10 (signal has 100 times the total power of the noise); see example traces of the observation under each SNR in Figure 5A.i – iv, and see example spectra under each SNR in Figure 5 —figure supplement 1. We repeat each simulation 1000 times with different noise instantiations for each signal to noise ratio.”
3) Code availability
I wanted to check the used bandpass filter parameters (since it's not described for the Blackwood algorithm) and noticed that the linked repository does not feature the code to recreate the plots in the manuscript, possibly it would be helpful for other researchers to provide that, also for comparing to other algorithms in the future.
Thank you for raising this important point. We have created an independent repository with code to run the simulations and to recreate the plots in the manuscript. Please see, https://github.com/EdenKramerLab/SSPEpaper.
https://doi.org/10.7554/eLife.68803.sa2Article and author information
Author details
Funding
National Institutes of Health (R01 EB026938)
 Anirudh Wodeyar
 Alik S Widge
 Uri T Eden
 Mark A Kramer
National Institutes of Health (R01 MH119384)
 Alik S Widge
National Institutes of Health (R01 MH123634)
 Alik S Widge
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Senior Editor
 Laura L Colgin, University of Texas at Austin, United States
Reviewing Editor
 Frances K Skinner, Krembil Research Institute, University Health Network, Canada
Reviewer
 Milad Lankarany, The Hospital for Sick Children, University of Toronto, Canada
Publication history
 Preprint posted: March 25, 2021 (view preprint)
 Received: March 26, 2021
 Accepted: September 24, 2021
 Accepted Manuscript published: September 27, 2021 (version 1)
 Version of Record published: October 22, 2021 (version 2)
Copyright
© 2021, Wodeyar et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 830
 Page views

 114
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.