Time-resolved phase-amplitude coupling in neural oscillations

Cross-frequency coupling between neural oscillations is a phenomenon observed across spatial scales in a wide range of preparations, including human non-invasive electrophysiology. Although the functional role and mechanisms involved are not entirely understood, the concept of interdependent neural oscillations drives an active field of research to comprehend the ubiquitous polyrhythmic activity of the brain, beyond empirical observations. Phase-amplitude coupling, a particular form of cross-frequency coupling between bursts of high-frequency oscillations and the phase of lower frequency rhythms, has recently received considerable attention. However, the measurement methods have relatively poor sensitivity and require long segments of experimental data. This obliterates the resolution of fast changes in coupling related to behavior, and more generally, to the non-stationary dynamics of brain electrophysiology. We propose a new measure of phase-amplitude coupling that can resolve up to one, optimally two, cycles of the underlying slow frequency component. The method also provides a measure of the coupling strength, for augmented insight into the mechanisms involved. We demonstrate the technique with synthesized data and compare its performances with existing methods. We also show that the method reveals rapid changes in coupling parameters in data from the entorhinal cortex of a free-behaving rat. The time-resolved changes revealed are compatible with behavior and complement observed modulations of oscillatory power. We anticipate that this new measure of dynamic phase-amplitude coupling will contribute to accelerate research into the dynamics of inter-dependent oscillatory components related to brain functions and dysfunctions.


Introduction
Classical studies of the role of neural oscillations in brain functions and behavior have reported on oscillatory rhythms within distinct, bandlimited frequency ranges (Goldman et al., 2002;Klimesch, 1998;Tallon-Baudry and Bertrand, 1999) (see (Cohen, 2008) for a review). The oscillations that compose brain rhythms are known to be interdependent across frequency bands (Buzs aki and Draguhn, 2004). This form of interaction known as cross-frequency coupling is a phenomenon readily observed in electrophysiology at multiple scales and with a range of experimental techniques (Buzsaki, 2006). Phase-amplitude coupling (PAC) is one of the best-studied subtypes of cross-frequency coupling (Canolty and Knight, 2010), accounting for the more frequent occurrence of higher-frequency bursts at preferred phases of underlying low-frequency cycles. PAC has received tremendous attention recently, with several studies revealing modulations of such coupling depending on task and resting-state conditions in health and disease (Tort et al., , 2009Axmacher et al., 2010;Fell and Axmacher, 2011;Schutter and Knyazev, 2012;Florin and Baillet, 2015;De Hemptinne et al., 2015).
The physiological relevance for PAC is in the assumption that slow oscillations mark the cycles of relative net excitability of neural ensembles (Von Stein and Sarnthein, 2000;Fries, 2005;Haider et al., 2006;Lisman and Buzs aki, 2008), which in turn pace the occurrence of neural spiking and that of faster post-synaptic activity, marked by high-frequency and often broadband bursts (Canolty and Knight, 2010).
One can further hypothesize that such coupling is transient by nature, reflecting the elusive dynamics of polyrhythmic brain activity. Some early task-related evidence of this assumption was demonstrated in humans by Tort et al. (2008). Ideally, measures of PAC need to provide the best possible frequency estimates of the oscillatory components related in phase and amplitude through this form of coupling. We would also need to assess the strength (intensity) of such coupling, to evaluate how it might be affected by behavior or physiopathological mechanisms. Finally and ideally, these measures would need to be accessible at the best possible temporal resolution, to detect and track such PAC changes at the natural "speed" of brain activity. This latter aspect has proven to be methodologically challenging, essentially because of the relatively poor signal-to-noise ratio affecting the higher-frequency portion of electrophysiological brain signalsin particular with non-invasive measures such as electro (EEG) and magnetoencephalography (MEG).
Event-related phase-amplitude coupling (ERPAC) was recently introduced in an attempt to bridge that methodological gap (Voytek et al., 2013). The method indeed provides a PAC measure with high temporal resolution, however under the constraint that expected PAC changes are time-locked to a stimulus or event of interest across repeated trials. This aspect restricts the sensitivity and the range of applications of ERPAC to event-related experimental designs and analyses.
More recently, Dvorak and Fenton (2014) proposed to estimate PAC over two complementary global and local time scales. The global time scaledefined over 10 s or more, typicallyis to identify the frequencies of the most coupled pair of oscillatory components: a slower oscillation (frequency for phase, f P ), which phase modulates the amplitude of faster bursts expressed at frequency f A (frequency for amplitude). In turn, a local time scaledefined over the f A cyclesis for detecting time variations in coupling strength. One identified issue with the approach though is that it assumes stationarity in coupled frequency pairs (a.k.a. modes) (f A , f P ) over possibly long periods of observations: this is an unlikely eventuality in neurophysiology.
The strength of PAC is another measure of interest to obtain dynamically (Tort et al., , 2009). However, existing methods to assess PAC strength are also challenged by poor temporal resolution (see (Tort et al., 2010) for a review).
We propose a new approach and practical method to address these issues in the widest range of experimental settings (task-related and spontaneous, ongoing electrophysiological data).

Principles
Formally, let us define PAC as the modulation of the amplitude A fA of an oscillatory component of frequency f A by the phase ϕ fP of a slower rhythm of frequency f P , with f P < f A . Our approach to derive a timeresolved measure of PAC (tPAC) was inspired by Cohen (2008) and Canolty et al. (2006), which methods we combined and optimized to achieve the best possible temporal resolution and sensitivity to coupling strength. The methodological steps are detailed below and summarized in Fig. 1.
The principle of the tPAC procedure is that it searches for the f P oscillation with strongest PAC coupling with f A bursts, over time windows that slide along the electrophysiological data signal xðtÞ. The user is first required to define a spectral range of interest for f P and f Afor instance, ½f Pmin ; f Pmax (e.g., [2,12] Hz) for the f P range. The range for f A is then subdivided into twenty centre frequencies either linearly or logarithmically, determined by the user. The other parameter that should be defined by the user is the length of sliding time window. This window should be long enough to cover at least one full cycle of f Pmin .
For each f A centre frequency tested, xðtÞ is bandpass filtered around f A with a zero-phase-shift, even-order, non-causal finite impulse response (FIR) filter. The bandwidth of each filter around each f A centre frequency is defined as the maximum between i) the difference between consecutive f A centre frequencies, and ii) the highest tested frequency for f P . Thus, the filter bandwidth spans the interval between consecutive f A frequency candidates and is inclusive of the range of interest for possible f P oscillations. Note that to minimize the filter edge effects and increase frequency resolution, bandpass filtering needs to be performed on the full-length signal before extracting its components in the sliding time window.
The amplitude envelope A fA ðtÞ of the bandpass filtered signal around the f A centre frequency, x fA ðtÞ is then extracted using the Hilbert transform before a sliding time window is applied to assess dynamic changes in PAC.
If in the current time-window, the amplitude of f A oscillations is coupled to the phase of a slower rhythm oscillating at frequency f P , then the power spectrum density (PSD) P A of the windowed A fA ðtÞ is expected to feature a peak at f P . Further, to avoid registering spurious manifestations of PAC, a peak around f P shall also be expected in P x , the PSD of the original data signal xðtÞ limited to the same time window, as recommended by Aru et al. (2015). P x and P A are estimated using the Discrete Fourier Transform (DFT) magnitude coefficients. The best candidate for the frequency for phase (f Ã P ) corresponds to the frequency of the highest peak of P A coinciding with a peak of P x . A tolerance threshold for the correspondence between peaks of P A and P x is set to the maximum between 1.5 the power spectrum's resolution (i.e. the inverse of the time window length) and 1.5 Hz. If no actual peak or peak coincidence is found, the coupling intensity is set to zero and an arbitrary valueoutside the range of interestis assigned to f P . For robustness purposes, the peaks in P x with an amplitude below 10% of that of the highest peak (2) the envelope of the resulting bandpass filtered signal is extracted using standard Hilbert transform; A short temporal window slides on the resulting envelope: (3) the power spectrum of the windowed envelope (PA) is estimated and its peaks are extracted; (4) the power spectrum of the original signal in the same time window is estimated, and is used for finding the highest peak in PA that co-occurs with a peak in Px. This determines the dominant frequency for phase f Ã P ; (5) windowed xðtÞ (with a signal buffer on both sides) bandpass filtered around the selected frequency for phase f Ã P , (6) and its analytic phase is obtained via Hilbert transform; (7) for each time point, the amplitude of the fast-oscillation envelope and the instantaneous phase of the slow oscillation are reported using a polar vector; (8) the Euclidean norm of the summed vectors averaged over an integer multiply of fP cycles is a measure of coupling strength within these time and frequency intervals. This value is further normalized with respect to the magnitude of x fA ðtÞ, to minimize the influence of signal magnitude in the measurement of coupling strength. These steps are repeated for all predefined fA frequency bins, and for all sliding time windows. in the interval of interest for f P are ignored. The reason for adding this step was that, as typical for the purpose of optimizing the FFT's computational time, the number of time points was set to the closest power of 2 (equal or larger than) of the data length, with zero-padding. This latter corresponds to signal windowing with an unknit boxcar in the time domain. Consequently, this produces spurious secondary peaks in the spectrum of the DFT magnitude coefficients. To avoid the detection of these side lobes as true signal peaks in P x , we propose the aforementioned threshold applying procedure for the robust detection of the dominant slow oscillation in P x . Furthermore, to prevent missing the peaks on the border of f P range of interest, in peak detection from P x and P A , one extra point from both sides of the f P interval were included. Note that decreasing the window length would potentially improve the temporal resolution of PAC estimation, but is detrimental to the detection of f Ã P , as it imposes a lower resolution to the f P search.
Next, to extract the phase time-series of slow f Ã P oscillation, xðtÞ in the time-window of interest is narrowband filtered around f Ã P (zero phaseshift, even-order, non-causal finite impulse response (FIR) filter, bandwidth ¼ 3 Hz) to yield x fP ðtÞ, which analytic phase ϕ fP ðtÞ is extracted using the Hilbert transform. The FIR filters used in tPAC analysis have transient responses that should be discarded; a data buffer is a tool that is used for this purpose: 2 s of extra data from both sides are included analysis. This data buffer is neither used in the detection of f P nor in the estimation of coupling strength. In the first and last few windows where extra data are not available, the signal was zero-padded.
To estimate the strength of coupling and the preferred phase of PAC occurrence along the f P cycle, the modulation phase and amplitude at each time point are pooled into the complex signal zðtÞ ¼ A fA ðtÞ:e ϕ f P ðtÞ (Canolty et al., 2006). A normalized estimation of coupling strength is obtained by computing the Euclidean norm of z. z is the vector sum of zðtÞ across all time samples (in the current sliding time window), divided by the average amplitude A fA ðtÞ in the same time window (Eq. (1)). Only full cycles of f P oscillations are considered in the estimation of z. The remaining part of the signal towards the end is cropped to avoid biasing the vector sum towards phase samples that would de facto repeat more often than others. To achieve this aim, we identified the phase angle at the start of each sliding time window, and used as effective signal length, that of the maximum number of cycles that would fit within the limits of the time window. In other words, while tracking the rotation of zðtÞ vectors in the polar plane, we excluded the time vectors that came after the last occurrence of the phase measured at the onset of the time window. This was to prevent non-uniform sampling of f P phase angles. The preferred coupling phase (ϕ tPAC ) was then estimated as the phase of the normalized vector z.
(1) t P denotes the data length used in averaging (i.e. equal or shorter than sliding window length, Δ); and k is the number of f P 's full cycles included in the sliding window. The signal time window is then moved forward with an overlap that can be adjusted by the user. The procedurefrom detecting the best f P to estimating coupling intensityis repeated for all successive time windows. This iterative process is subsequently repeated for all f A centre frequency candidates. The resulting measures are pooled in a summary array, which can be conveniently represented as the map, as illustrated in Fig. 1.
In a real experimental context, the actual value of the driving frequency for phase (f Ã P ) is unknown a priori. As previously mentioned, the minimum tested value for f P is the most important factor for defining the temporal resolution of the PAC measure, as it imposes a minimum duration to the window length applied to the data. Indeed, at least one full cycle of the tested f P is required to assess the possible interdependence of the f P phase cycle with the amplitude of f A bursts. Indeed, the presence of a dominant slow oscillation at the frequency for phase is necessary for genuine PAC coupling (Aru et al., 2015). Hence, peaks in the power spectrum of the data help identify potential frequency candidates for f P . A relatively wide f P range can help capture multiple, coexisting modes of coupling. If necessary and if the detected f Ã P s were notably faster than the minimum f P of interest, a second pass of tPAC can be run with a narrower range of f P of interest, around the previously identified f Ã P . This would also provide an opportunity to increase the temporal resolution of the tPAC analysis, with shorter time windows.

Sparse estimation of time-resolved PAC
The parameters of interest in PAC estimation span time, frequency for phase, and frequency for amplitude. Coupling strength can in principle be estimated at all points in this 3-D subspace of unknowns, however at a very high computational cost that is certainly not justified and necessary for all applications (N fP Â N fA Â N time windows PAC estimations, per data time series). We can also safely assume redundancy in the range of unknowns, for instance by considering that oscillatory bursts at frequency f A are PAC modulated by one unique slow f P oscillation within each sliding time window.
In practical terms, the method assumes that each f A candidate is coupled to a single f P frequency. Therefore, the resulting comodulogram i.e. a 2-D map indicating the strength of coupling between different oscillations, with f P values along the x-axis, and f A values along the y-axis over a sliding time window does not contain more than one peak in each f A row across all possible f P tested frequencies. Therefore, instead of testing several f P frequencies for each f A , tPAC detects the single most coupled f P oscillation coupled to an f A centre frequency candidate, and estimates the coupling strength for ðf P ; f A Þ. This yields a sparse representation of the comodulogram.
For visualization purposes, the sparse comodulogram array can be projected on 2-dimensional hyperplanes, by averaging over the third dimension. The possible produces of such projections are: tvs:f A , tvs:f P , or f P vs:f A representations.

Statistical tests for tPAC
Statistical significance of tPAC parameter estimates is verified with a non-parametric resampling techniques. We followed the recommendations by Aru et al. (2015), suggesting to generate surrogate datasets using a block-resampling approach. A benefit is that phase distortion is minimized, which reduces the rate of false-positive detections. The envelope time-series of the f A oscillation in each time-window were first split into five blocks. These blocks were then randomly permuted to yield a surrogate dataset that realizes the assumption of absence of PAC beyond chance levels. This approach were acknowledged to produce relatively conservative assessments of statistical significance of PAC measures (Aru et al., 2015).
In the experimental work reported below, n ¼ 500 surrogate trials were produced per time-window and f A centre frequency candidate. The coupling strengths within the 95 th percentile of the surrogate distribution were considered statistically significant. Further, when assessing significance of PAC values reported in the time vs. f A , and to compensate for multiple-comparison testing, the null distributions were that of the maximum statistic across all time windows and tested f A frequency candidates (Pantazis et al., 2005).

Experimental data
The comparison of tPAC against other methods was conducted with controlled, synthesized data. We then used electrophysiological recordings in animal preparations to relate observed tPAC effects to behavior. All electrophysiological data processing was performed using Brainstorm (Tadel et al., 2011), using default parameters, which distribution also contains the open Matlab source code of tPAC. We also share the Matlab scripts for reproducing and the synthesise of more test data (http://neuroimage.usc.edu/brainstorm/Tutorials/TutPac).

Synthesized data
Data simulations were used to provide ground-truth signals, with controlled PAC parameters. It consisted of times series obtained from the sum of two sinusoids, with additive noise. The amplitude of the highfrequency component was modulated by the phase of the low frequency component according to Eq. (2).

> > > > <
> > > > : K fP and K fA are fixed scalars that determine the peak amplitude of f P and f A sinusoidal components, respectively. ϕ c indicates the phase of x fP ðtÞ where the magnitude of x fA ðtÞ bursts is maximum. χ 2 ½0; 1 determines the fraction of x fA ðtÞ that is not modulated by x fP ðtÞ and therefore controls for coupling strength, defined as 1 À χ. Fig. 2 shows two examples of synthesized signals, with different coupling strengths.
The additive noise εðtÞ was generated from two components: random samples drawn from a power law, to mimic background brain activity, and white Gaussian realizations for simulating instrumental noise. The power of the Gaussian component was set to half that of the powerlaw samples.
We provide more details in the Appendix for advanced simulations where the duty cycle of x fP ðtÞ is manipulated, to test whether asymmetrical cycles affect tPAC and other PAC measures (Belluscio et al., 2012).

Electrophysiological data
Local field potential (LFP) recordings from multiple mesio-temporal brain regions of Long-Evan rats were collected using eight-shank multisite silicon probes (200 μm inter-shank distance). The freely-available data was published elsewhere (Mizuseki et al., 2009) and shared on the Collaborative Research in Computational Neuroscience website (http://crcns.org) by the Buzsaki laboratory. Recordings were carried out while the animals ran on a linear track (length: 250 cm), digitally sampled at 20 kHz (16-bit resolution, RC Electronics) and bandpass filtered (1 Hz-5 kHz). The processed LFPs were downsampled to 1250 Hz (anti-aliasing filter applied). Two LEDs were connected to the animal's head. The position of the rat was extracted from the video file recorded during free behavior on the linear track, with a sampling rate of 39.06 Hz, and time aligned with the LFP recordings.

Time-resolved estimation of PAC properties: frequency for phase, amplitude and coupling strength
The results of phase-amplitude coupling analysis of a synthesized data set using tPAC are illustrated in Fig. 3. This experiment covered two aspects of PAC analysis: the detection of coupling modes (Fig. 3a-c) and the estimation of coupling strength (Fig. 3d).
In the first part, the input signal was a synthesized data trace, generated using the model explained in section 2.4.1 (Fig. 3a). Multiple modes of coupling were present: during the first half of the signal (10 s duration), the phase of slow oscillation at f Ã P1 ¼ 9 Hz was coupled to the amplitude of a faster rhythm at f Ã A1 ¼ 115 Hz. In the second half, the first coupling mode was terminated and two other modes appeared simulta- Hz, and f Ã A3 ¼ 87 Hz, respectively. The signal-to-noise ratio was set to 6 dB, and the preferred coupling phase in the three modes are 270, 0, and 180 , respectively. The coupling strength was kept constance and identical in all modes.
The bands of interest for f P and f A in tPAC analysis were defined as [3,15] Hz and [20,200] Hz, respectively. The duration of the sliding time window was set to 0.75 s, hence 2.25 cycles of minimum f P of interest. Coupling strength was estimated for 20 centre frequencies distributed linearly along the f A range of interest.
In Fig. 3b, the top panel depicts the initial outcome of tPAC algorithm, which shows time variations of coupling strength with f A . This map reveals the time course of all three coupling modes expressed in the data. The time variations of coupling strength for each f P estimated from the sparsely estimated coupling strengths (see Sec. 2.2). The resulting f P vs. time map is shown in the bottom panel of Fig. 3b. The comodulogram was also reconstructed (Fig. 3c), as a typical, non time-resolved appreciation of PAC in signals.
To further illustrate tPAC's ability to estimate phase-amplitude coupling strength, a synthesized data trace with continuously timeevolving coupling strength was generated: f Ã P ¼ 4Hz, f Ã A ¼ 73Hz, SNR ¼ 5 dB, duty cycle ¼ 0.35, sampling rate ¼ 1000 Hz, 270 s signal duration. tPAC was applied to the data with f P and f A frequency ranges of interest, [2,15] Hz and [50, 140] Hz, respectively. The length of the sliding window was set to 0.53 s, which includes only 2 cycles of f Ã P , and at least one cycle of the minimum f P of interest (i.e. 2 Hz as explained above). The normalized time-resolved tPAC strength is shown in Fig. 3d. The correlation between the estimated time-variation of tPAC coupling and ground truth was 0.95 (p < 0:001). Increasing the sliding window length from 2 cycles of f Ã P to 4 and 8 cycles led to similar results (correlation coefficient of 0.97 and 0.99, respectively).
The dataset used consisted of 500 trials of short asymmetric (duty cycle ¼ 0:35) synthesized time series, with f Ã P ¼ 4Hz, f Ã A ¼ 73Hz, SNR ¼ 5 dB, and sampling rate of 1000 Hz. The ranges of interest for f P and f A were set to [2,15] Hz and [50, 140] Hz, respectively. The signal duration was 1.07 cycles of minimum f P of interest (f Pmin )i.e. 0.53 s ¼ 2.15 cycles of f Ã P . The signal length was a non-integer multiple of f P cycles, to make the situation more challenging and realistic, when, with physiological data, the signal length is not expected to be a multiple of f P cycles.
All tested methods (MVL, PLV, KL-MI, and APSD) were implemented in house, from the details provided in the respective publications. The methods parameters were identical to those used for tPAC. Zero-phase lagged FIR filters with 3-Hz and 15-Hz bandwidths were used for filtering in the f P and f A ranges of interest, respectively. Eighteen centre frequencies for f A were selected linearly in the [50, 140] Hz range. For the methods that scan on f P , fourteen centre frequency candidates were linearly distributed in the [2, 15] Hz range. The duration of filter transients (containing 99 percent of the energy from the filter's transient effect) for all f A filters was less than 75 ms (29-74 ms depending on the f A centre frequencynote that f A filtering was performed on the full-length input signal, not on the shorter, sliding time windows). The transient effects for f P filters were less than 1.03 s (289-1029 ms depending on f P centre frequencya data buffer of 2 s was added prior to filtering). Since MVL was originally introduced using a wavelet filter bank (Canolty et al., 2006), and wavelet filter banks are common in the PAC literature (Lakatos et al., 2005(Lakatos et al., , 2008Demiralp et al., 2007;Cohen et al., 2009;van der Meij et al., 2012), we also used an MVL implementation based on wavelet filtering. We used the actual code shared by Canolty et al. (2006) available in Brainstorm. Results from MVL-wavelet were obtained using a Morlet wavelet filter bank with 150 frequency candidates for f A and f P , that were distributed logarithmically within [2,150] Hz.
In all tested methods, the frequency pairs (f P , f A ) indicating the highest coupling strength was compared to the actual values used to generate the synthesized data. The relative errors in f A and f P were estimated (Eq. (4)). For each trial, the relative errors in f P and f A detection were averaged together, as a summarizing measure of accuracy. Fig. 4 shows the performances of the tested methods in terms of identification of the coupled (f P , f A ) pair, and quantification of coupling strength. Fig. 4a shows the relative errors in coupling detection for all algorithms. On these challenging short data segments, tPAC and APSD produced the smallest ( < 5%) relative errors.
We performed a systematic examination of the sensitivity of data length on PAC estimation accuracy, Fig. 4b. The simulation procedure described above was repeated with signal lengths increasing from 2.15 to 10.15 f P cycles. We found that for all tested methods, performances in coupling detection improved with increasing signal length, although at different rates. Overall, for signal durations longer than two f P cycles, tPAC and APSD produced the best estimates of f P and f A . As explained in Sec. 2.1, at least one full cycle of the minimum tested f P is required for PAC analysis. In this experiment, since the slowest tested f P was 2 Hz (a 500-ms cycle) with data where the unknown true f Ã P was 4 Hz (a 250-ms cycle), the shortest time window duration was 2-f Ã P cycles. Hence, we furthered the investigations concerning the temporal resolution of tPAC and other PAC estimators, by conducting another set of experiments where the smallest f P of interest (3.5 Hz) was brought closer to f Ã P ¼ 4 Hz (f P range: [3.5, 15] Hz; see Fig. 4b-right panel). These experiments confirmed that the observation of too few f Ã P cycles affects the accuracy of PAC measures. They also indicate that the performances of methods such as APSD that rely on the detection of f P peaks in the signal spectrum are negatively affected when the minimum tested f P is too close to f Ã P . Although tPAC also relies on such peak detection, we have implemented several procedures that augment the robustness of the method, in that respect, as demonstrated in Fig. 4b. For example, if f Ã P is close to the minimum f P of interest, APSD requires window lengths that cover at least 3 f Ã P -cycles, for best performances. Under the same conditions, coupling detection with MVL-wavelet became less stable with increasing timewindow lengths. MVL, PLV, and KL were less sensitive to the proximity of f Ã P with minimum tested f P , because these methods scan systematically and linearly the f P range of interest for f Ã P and do not rely on spectral estimations. Interestingly, this new frequency range for f P (narrower band, f Ã P closer to the bands lower bound) improved the performances of these algorithms slightly, for short signal duration of 2 f P cycles.
We tested how signal-to-noise ratio (SNR) affected PAC estimation performances. The same simulation procedure as shown Fig. 4a's was repeated, with SNR increased from À5 to þ7 dB (negative SNR values correspond to larger noise power compared to the signal's). Our results demonstrate that the performances of all tested methods go through the  (a)). (d) estimation of coupling strength (signal synthesis was identical as for panel (a)) with 3 tested coupling strengths: (1 À χ ¼ 0:2; 0:55 and 0:9), shown with horizontal dashed lines. Average estimates of coupling strength and corresponding SEM over 500 trials. (e) relative error in estimated coupling strength averaged at three tested coupling strengths. All methods showed significantly different performances (Kruskal-Wallis rank test, Dunn's test as post-hoc, ps < 0:001). See text for details on synthesized signals. same three behaviors depending on SNR: i) a flat response segment for SNR levels below À3 dB; ii) a dramatic improvement in performances for SNR between À3 and 3 dB; and iii) an asymptotic plateau in performances for SNR above 3 dB (3 dB means that signal's power was twice of the noise's power).
MVL-wavelet was the most robust method against noise, however with relatively large errors in estimating coupling parameters ( > 20%), even in high SNR conditions. In positive SNR conditions, APSD and tPAC had the best performances compared to MVL, PLV, and KL (Fig. 4c).
The ability of PAC methods to recover the coupling strength was also tested quantitatively (see Fig. 4d-e). Again, the simulation scheme of Fig. 4a was reproduced, with three different coupling strengths ranging from low (1 À χ ¼ 0:2), to intermediate (1 À χ ¼ 0:55), then high (1 À χ ¼ 0:9). Fig. 4d shows the coupling strengths recovered from the tested methods at the actual ðf P ; f A Þ mode. Bars indicate the SEM across trials. All PAC measures were normalized to their maximum possible value, as obtained in noiseless and maximum coupling strength (1 À χ ¼ 1) conditions.
The relative errors in coupling strength, on all trials with three coupling strength values, are depicted in Fig. 4e. The data demonstrate that MLV-wavelet and tPAC manifested the best ability to recover coupling strength quantitatively, while PLV was the most challenged method in that respect. Although tPAC uses an algorithm similar to MVL for coupling strength estimation, several optimized schemes in tPAC (namely the normalization to A fA power, window lengths reduced to multiple cycles of f P , as described in section 2.1) resulted in significantly higher performances for tPAC (error: 14:88±0:27%) in coupling strength estimation, compared to MVL-wavelet (error: 19:33±0:43%) and MVL (error: 27:41±0:49%), with short input signals (All methods showed significantly different performances: Kruskal-Wallis rank test, Dunn's test as post-hoc, ps < 0:001).
Taken together, tPAC demonstrated the best consolidated performances in estimating the frequency modes expressing PAC (f P and f A ), and in recovering coupling strength. APSD showed high performances in detection of coupled frequency pair (f P , f A ). However the APSD estimation of coupling strength was poor (error > 40%). MVL performed well in estimating the coupling strength, but remained poorly accurate in recovering the coupled frequency pair, especially on short data lengths (error > 15%). In MVL, using a wavelet filter bank produced mixed effects. More accurate coupling strength estimation, and less sensitivity to noise were the positive effects of using the wavelet filter bank. However, the wavelet version of MVL produced higher coupling detection error for short input ( < f Ã P -cycle), not stable results with f Ã P very close to the minimum f P of interest.
These performances obtained with the shortest lengths of data input indicate that using tPAC, all phase-amplitude coupling parameters can be estimated in a time-resolved fashion, using short (length of two f P cycles minimum) time windows. The next section details how we extended testing to electrophysiological data related to behavior.

Electrophysiological data
The time-resolved tPAC coupling strength in the third layer of entorhinal cortex (EC3) of a rat freely moving on a linear track was obtained, after setting the frequency of ranges of interest to [2,15] Hz for f P , and [35,215] Hz for f A . Coupling strength was calculated for 20 centre frequencies distributed linearly in the f A range of interest. Data analysis was performed at two different temporal scales: one coarse assessment using 10-s overlapping windows; and a finer-grained investigation using sliding windows of 2.5-s duration. In both cases, the sliding windows were moved along the data time series with a 50% of overlap.
We show the tPAC coupling strengths both as comodulograms (f P vs. f A ), and time-resolved maps (Figs. 5a and 6a, respectively -all maps were smoothed using 2-D interpolation, with details provided in the supplementary material). The overall strongest PAC mode of the coarse analysis was found between the phase of theta (7.5 Hz) and the amplitude of fast gamma ([90À130] Hz) oscillations (Fig. 5a). We verified the tPAC findings with the comodulogram obtained from MVL -wavelet (Canolty et al., 2006) (Fig. 5b): the results were consistent in terms of the principal ðf P ; f A Þ mode detected. Further, since tPAC detects the most coupled f P to each time-window and f A sub-band, only the dominant f P modes were represented in the tPAC comodulograms. This, by construction, makes the resulting coupling maps look less smeared than the ones obtained by the other tested methodssee section 2.1 for more details). Also akin to (Canolty et al., 2006) and for visualization purposes, the spectrogram and the epoched LFP signal were averaged time locked to the troughs of the 7.5-Hz f P cycle (Fig. 5c spectrogram is depicted relative to the averaged power across time). This confirmed that the power of faster oscillations was indeed modulated by the phase of the slow rhythm, as detected by tPAC.
The dynamics of cross-frequency coupling related to behavior were investigated further, using time-resolved tPAC maps. Our results show that burst-like variations in coupling strength occurred in the recordings ( Fig. 6a; displayed values are above 95th percentile of surrogate data distributioncorrected for multiple comparisons). The animal's position on the track (showed with gray dashed lines in Fig. 6a) indicated that such bursts of higher coupling occurred when the animal was running. We confirmed these observations by computing the PSD during respective periods of rest and running. The power spectra showed peaks in the theta band ([5, 10] Hz) in both cases, although at distinct frequencies (6 Hz and 7.5 Hz, respectively. Fig. 6c). The time-frequency decomposition of the LFP signal also showed increase in fast-oscillation power during movement (Fig. 6b). However, such increases spanned a wide frequency range ([80, 350] Hz), and were not specific of the f A band we found was coupled to the phase of theta oscillations.
The PSD data further revealed that in addition to the dominant peaks in theta range, there was at least another smaller peak in the delta band ( < 4 Hz) in both behavioral conditions (Fig. 6c). Previous studies reported on such slower rhythms in similar behavioral conditions, as possibly generated by the midline thalamic nucleus (Zhang et al., 2012).
All data segments consisting of at least 10 consecutive seconds of running followed by 10 s of resting at one end of the track were extracted. They were temporally aligned, with time 0 ms corresponding to when the animal stopped at the end of the track.
To detect possible event-related changes of cross-frequency coupling, tPAC was computed on 9 epochs of 20 s, collected on 8 channels in EC3, with a sliding window length of 2.5 s and 50% overlap, using the same ranges of interest for f P and f A as previous analysis. Our decision to set the minimum f P of interest at 2 Hz was guided by the available slow peaks in the power spectrum of the data, and behavioral aspects we wished to relate to dynamic PAC changes, namely the episodes of animal running vs. resting. These latter unfold over several seconds (about 10s each typically), hence the selection of a relatively slow minimum f P , and therefore long tPAC window of 2.5s to sample these events in time. We also chose a larger time window because of its potential for revealing multiple, co-existing PAC modes, at slower f P . Fig. 6d illustrates the time-resolved fluctuations in coupling strength along the epoch duration. The figure also shows the respective comodulograms of the running and resting periods. Our findings confirmed that the dominant coupling mode during running was between theta (7.5 Hz) and fast gamma ([90,130] Hz), and that the strongest coupling was observed in the running behavior. Two other coupling modes were revealed: between delta and low gamma ([40, 60] Hz)lower, but significant strengthand between lower theta and midrange gamma ([60,100] Hz) only during rest. These two modes were not clearly distinguishable from the non time resolved PAC analysis. The frequency of theta oscillations that drove fast gamma during running was close to the dominant global theta mode (7.5 Hz). The theta coupled with midgamma during rest was lower in frequency (about 6 Hz). This observation was consistent with the dominant peak of theta observed in the power spectrum of the full data. Our results also indicate that delta oscillations in EC3 were coupled to low-gamma rhythms, specifically during the resting behavior. The frequency range of the three gamma components identified were consistent with the three modes of gamma in hippocampus reported previously by Buzs aki and Schomburg (2015).
These results illustrate the insight enabled by time-resolved PAC measures with tPAC. Further investigations would be required to actually comprehend the functional relevance of the distinct, behaviorallydependent coupling modes observed.

Discussion
We introduced tPAC, a method to measure time-resolved parameters of phase-amplitude coupling in electrophysiological time series. tPAC performs comparatively better than existing methods for estimating at once the frequencies of coupled oscillations, the strength of such coupling, all using the shortest input signal length, with relative immunity to noise. One benefit of tPAC is therefore to resolve PAC measures in time. This property alleviates the constraint of working previously with long data epochs, which was not compatible with the expected non stationarity of electrophysiological signals, for instance, in relation to behavior. Importantly, non stationarity may induce spurious spectral correlation between Fourier coefficients, and therefore yields artifactual cross-frequency coupling (Aru et al., 2015). This adds to the significance of time-resolved PAC methods.
tPAC was compared to selected published methods, with an emphasis on data length. The results using ground-truth synthesized data showed that tPAC is sensitive and accurate with at minimum a signal duration of two f P cycles. Our data also showed that MVL-wavelet and KL-MI PAC methods performed relatively poorly on shorter input signalsespecially when minimum f P of interest is very close to f Ã P . This finding is consistent with claims by Tort et al. (2010), indicating that a minimum of 200 f P cycles (about 30 s) are required for optimal performance of these methods. The necessity for long input signals for MVL was also discussed by Penny et al. (2008).
One practical solution to work against poor temporal resolutions was previously to concatenate short recording epochs. Although this indeed increases data length, the discontinuities induced by concatenation can lead to spurious high-frequency components, which may affect coupling estimation. For example, when epochs of identical lengths are concatenated, discontinuities induce an artifactual PAC regime that relates sharp variations (broadband high-frequency spectral contents) to the epoch edges, with an implicit oscillatory cycle corresponding to the epoch duration. Consequently, spurious PAC may occur .
Further, the estimation of coupling strength is another asset of the tPAC method: most published methods do not properly recover this parameter (Tort et al., 2010). Hence, tPAC provides higher temporal resolution to the analysis of phase-amplitude coupling of ongoing data (e.g., resting-state, epilepsy, sleep recordings, etc.), and time-resolved sensitivity to coupling strength. Note that in the context of event-related studies, with enough trial repetitions (above 150, typically), ERPAC (Voytek et al. (2013)) can provide the fastest possible temporal resolution, reaching up to the data's sampling rate.

Conclusion
We propose tPAC as a new time-resolved method for estimating crossfrequency phase-amplitude coupling in electrophysiological signals. tPAC combines the highest temporal resolution, the capacity of estimating coupling strength, and the lowest sensitivity to noise conditions, even for shorter data lengths. These properties are key to reveal transient coupling variations related to behavior, from continuous and eventrelated data collected with a range of electrophysiological techniques that are relevant to the neuroimaging community, including EEG and MEG sensor and source data series.