Precision Timing with α–β Oscillatory Coupling: Stopwatch or Motor Control?

Abstract Precise timing is crucial for many behaviors ranging from conversational speech to athletic performance. The precision of motor timing has been suggested to result from the strength of phase–amplitude coupling (PAC) between the phase of alpha oscillations (α, 8–12 Hz) and the power of beta activity (β, 14–30 Hz), herein referred to as α–β PAC. The amplitude of β oscillations has been proposed to code for temporally relevant information and the locking of β power to the phase of α oscillations to maintain timing precision. Motor timing precision has at least two sources of variability: variability of timekeeping mechanism and variability of motor control. It is ambiguous to which of these two factors α–β PAC should be ascribed: α–β PAC could index precision of stopwatch-like internal timekeeping mechanisms, or α–β PAC could index motor control precision. To disentangle these two hypotheses, we tested how oscillatory coupling at different stages of a time reproduction task related to temporal precision. Human participants encoded and subsequently reproduced a time interval while magnetoencephalography was recorded. The data show a robust α–β PAC during both the encoding and reproduction of a temporal interval, a pattern that cannot be predicted by motor control accounts. Specifically, we found that timing precision resulted from the trade-off between the strength of α–β PAC during the encoding and during the reproduction of intervals. These results support the hypothesis that α–β PAC codes for the precision of temporal representations in the human brain.

Recent investigations have implicated beta (β) activity to be involved in timekeeping processes (Kononowicz, Roger, & van Wassenhove, 2018;Wiener, Parikh, Krakow, & Coslett, 2018;Kulashekhar et al., 2016;Bartolo & Merchant, 2015;. During the production of temporal intervals, the power of β activity scales with the duration of produced intervals (Kononowicz, Roger, & van Wassenhove, 2018;Kononowicz, Sander, & van Rijn, 2015). Yet, the consideration of only one frequency band could be oversimplifying the understanding of timekeeping mechanisms. For example, coupling between neural oscillations may be important for the coding (Lisman & Jensen, 2013) and transmission (Fries, 2005) of information in the brain. One particular cross-frequency coupling regime, known as phase-amplitude coupling (PAC), reflects the interaction between the phase and amplitude of two neural signals. PAC may support cognitive (Roux & Uhlhaas, 2014;Tort, Komorowski, Eichenbaum, & Kopell, 2010;Tort, Komorowski, Manns, Kopell, & Eichenbaum, 2009) and motor (De Hemptinne et al., 2013) functions. In a previous study, we found that the power of β oscillations was regulated by the phase of alpha (α) oscillations during motor timing (Grabot et al., 2019) so that the strength of α-β PAC was commensurate with the precision of the internally generated timed action: The narrower the behavioral distributions, the stronger the α-β PAC. These findings suggest that stronger α-β PAC allows for a better precision of temporal representations.
In the study by Grabot et al. (2019), participants were instructed to produce a time interval as close as possible to a required target interval using two button presses: The first button press started the interval, and the second terminated it. However, in time production tasks, it is difficult to distinguish between motor, perceptual, and central cognitive processes, as all processes intervene in the production of the time interval (Keele, Pokorny, Corcos, & Ivry, 1985). To the contrary, in a reproduction task, each trial starts with the presentation of a target duration that the participant needs to subsequently reproduce via a motor action. The separation in two task stages, namely, duration encoding and duration reproduction, provides a handle to investigate timing processes with a limited motor component in the encoding stage. This premise was supported by Baudouin, Vanneste, Isingrini, and Pouthas (2006) who showed that time production tends to be primarily associated with spontaneous tempo tasks, involving tapping at the preferred rate as regularly as possible, whereas time reproduction tends to be associated with working memory measurements, providing distinct cognitive blueprints in both tasks. Here, we will focus on the difference between duration encoding and duration reproduction in a time reproduction task. Focusing on time reproduction as opposed to time production may enable us to differentiate motor from nonmotor task components as if α-β PAC is linked to the precision of motor control, we should only observe effects during duration reproduction, whereas if the precision of more general, higher cognitive processes that are timing related drives α-β PAC, effects should be observed during both encoding and duration reproduction.
To differentiate motor control and higher cognitive processes contributing to time encoding, we reanalyzed previously published data ) from a temporal reproduction task. In every trial, participants were randomly presented with 2-, 3-, or 4-sec intervals, which they had to reproduce. We set out to investigate if α-β PAC is present in encoding and reproduction stages and whether α-β PAC during encoding or reproduction is related to the precision of temporal durations.

METHODS
We reanalyzed the data collected in the study of  to investigate the role of α-β PAC during the encoding and reproduction of 2-, 3-, and 4-sec time intervals.

Participants
Eighteen students enrolled at the Humboldt, Freie, or Technical University of Berlin (Germany), with no selfreported hearing or vision loss or neurological pathologies, took part in the experiment in exchange for monetary compensation. A written informed consent approved by the Ethical Committee Psychology of the University of Groningen was provided to each participant before the beginning of the experiment. The data of two participants were discarded from the analyses: One participant fell asleep during the experiment, and the second one exhibited unreliable and spurious patterns of PAC. The final sample was composed of 16 participants (all right-handed; eight men; between 19 and 29 years old).

Stimuli and Procedure
In a given trial, participants were presented with a to-bememorized time interval, called the encoding interval (EI; Figure 1A), which they subsequently were asked to reproduce (reproduction interval [RI]).
Each trial started with the presentation of a fixation cross "+" and a waiting time randomly sampled among 1, 1.5, 2, or 2.5 sec. The EI was then presented as an empty interval delimited by two auditory tones (5-msec duration, 1 kHz, ∼75 dB) separated by 2, 3, or 4 sec. The fixation cross "+" lasted for the entire time of the EI presentation and for the 1.5-or 2.5-sec waiting period after the EI. After this waiting period, participants could start the RI. The start of the RI was indicated by the same tone as the EI and with a visual change in the fixation cross from "+" to "x." Participants pressed the mouse button with their right hand when they thought the "x" was on the screen for the same amount of time as the "+" was. Pressing the response device button initiated the onset of a tone presentation, and the "x" disappeared from the screen. After an intertrial interval of 2.5, 3.5, or 4.5 sec, the next trial started. In all trials, the EI was followed by the RI. One hundred twenty trials were presented in 12 blocks of 10 trials for each timing condition. Each block contained at least three repetitions of each EI pseudorandomly presented. Each EI was thus presented 40 times.
Between blocks, participants received adaptive feedback on their performance indicating how many trials were correct. The range of correct feedback was initially set to 20% deviation of the target interval and then dynamically adjusted by decreasing (−2.5%) or increasing (+0.5%) the range after each correct or incorrect trial, respectively. Participants were instructed to reproduce durations as accurately as possible and to maximize their number of correct trials in each block. Before the start of the MEG recordings, participants were presented with five practice trials. Stimuli were presented using a PC running Presentation software (Neurobehavioral Systems).

MEG Recordings
In a dimly lit, standard magnetically shielded room (Ak3b, Vacuumschmelze) located at the PTB Berlin, each participant laid in supine position with eyes open. Visual stimuli were presented on a screen via a projector located outside the magnetically shielded room. Participants crossed their arms in front of their chest, as earlier work has shown this position to be most comfortable during longer recording sessions, and responded by clicking a button on a computer mouse, which was held in the right hand, located on their left upper body. Measurements were carried out with a Yokogawa MEG system ( Yokogawa Electric Corporation) using 125 axial gradiometers with three reference magnetometers used for offline data denoising and ambient noise removal.
The head position with respect to the sensor helmet was measured using coils attached to the scalp at anatomical landmarks (nasion and preauricular points). The locations of the coils were digitized with respect to three anatomical landmarks with a 3-D digitizer (Zebris).
The brief auditory bursts indicating the start and end of the intervals were presented binaurally via MEGcompatible tube earphones (Etymotic Research) at sound levels set for each participant individually, at approximately 75 dB.

Data Preprocessing
The data were preprocessed using FieldTrip (Oostenveld, Fries, Maris, & Schoffelen, 2011;Version 20130713) and a custom-written MATLAB (The MathWorks, Inc.) code. The MEG data were denoised by computing noise cancelation weights using three reference channels, by default incorporated in the Yokogawa system. First, PCA was computed on the reference channels and then projected on the gradiometer channels and subtracted from the data of interest. This computation was performed separately for each trial using FieldTrip's ft_denoise_pca function. A 50-Hz notch filter and a bandpass finite impulse response filter (0.01-130 Hz) were applied. Trials containing excessive ocular artifacts, movement artifacts, amplifier saturation, or Superconducting QUantum Interference Device (SQUID) artifacts were selected by visual inspection and excluded from further processing. Eye blinks, heartbeat, and muscle artifacts were corrected using independent component analysis.

Calculation and Statistical Assessment of PAC Using the Modulation Index
To prevent the influence of evoked responses in the PAC estimation, we focused on the time segment starting from 0.4 sec after the tone onset initiating the EI or the RI until the end of the EI duration (i.e., 2, 3, or 4 sec). In the RI analysis, we focused on the time segment from 0.4 sec after the stimulus onset up to 1.2 sec in the 2-sec condition. We chose 1.2 sec, the lower tail of the RI distribution, to ensure that the data were free of components associated with movement execution. For the same reason, we focused on the period extending up to 2 sec in the 3-and 4-sec RIs.
To calculate PAC, we split the data of each duration condition into temporal bins of the same length. This procedure ensured an equal number of time samples (given the variability in produced duration lengths) and prevented a possible erroneous assessment of the PAC calculation (Dupré la Tour et al., 2017). For this, we used bins of 0.8 sec long, starting 0.4 sec after the onset of the first stimulus in the EI and RI. The number of bins depended on the duration condition and ranged from one to four bins for the 2-to 4-sec trials, respectively. Across all conditions, Bin 1 ranged from 0.4 to 1.2 sec (for 2-, 3-, and 4-sec intervals); Bin 2, from 1.2 to 2 sec (for 2-, 3-, and 4-sec intervals); Bin 3, from 2.2 to 3 sec (for 2-and 3-sec intervals); and Bin 4, from 3.2 to 4 sec (for a 4-sec interval).
The strength of PAC was assessed using the modulation index (MI; Dupré la Tour et al., 2017;Tort et al., 2009): Raw data were bandpass filtered with a low-frequency bandwidth of 2 Hz and a high-frequency bandwidth of 20 Hz. The instantaneous amplitude of the high-frequency activity and the phase of the slow fluctuations were extracted using the Hilbert transform applied to epoched data. This analysis resulted in distributions of high-frequency amplitude as a function of the phase of the analyzed low-frequency activity. To assess whether the distributions diverged from uniformity, the Kullback-Leibler distance was computed and then normalized to provide an estimate of the strength of MI. The low-frequency component ranged from 3 to 15 Hz in steps of 0.5 Hz, and the high-frequency component ranged from 14 to 100 Hz in steps of 2 Hz. A comodulogram was computed for each MEG sensor, providing a full estimate of the MI over the full spectrum.
To assess the statistical significance of PAC at the individual level, the MI was compared to a surrogate distribution (n = 100) computed by randomly shifting the lowfrequency signal by a minimum of 0.3 sec, as was previously proposed by Tort et al. (2010). For each participant, five sensors with the highest z scores were selected and averaged. The resulting averages were plotted separately for the EI and RI ( Figure 2). The averaging of z scores indicated the level of significance so that a z score of 4 corresponded to a p value of .00001.
To assess whether PAC was present in the EI and RI, we compared the MI computed during the baseline (from −0.8 to 0 sec before the EI onset) to the MI computed during the EI and RI. Separate tests were performed to compare baseline activity against every bin and every condition to ensure that PAC was computed on an equal number of samples. The average outcomes of this comparison across bins and across durations are plotted in Figure 2A. Statistical assessments for the  Grabot et al. (2019). Slight differences in topographical maps between the two studies could be because of different MEG sensors (axial gradiometers in the current study and magnetometers in Grabot et al., 2019) and different positioning of participants in the experimental setting with respect to the scanner: horizontal in the current study and vertical in the study by Grabot et al. (2019). In addition, both studies had different trial numbers. There were around 600 trials in the study by Grabot et al., whereas in the current study, we had 40 trials per 2-, 3-, and 4-sec conditions. (C) Alpha (8-12 Hz) peak-locking plot. The individual time series were locked on the peak of the α oscillations (black line). β power (15-40 Hz) was computed in the realigned data, normalized by subtracting its mean, and displayed as heatmaps. For illustration, we picked the MEG sensor showing maximal α-β MI for one representative participant.
MI were performed using MNE-Python permutation_ cluster_1samp_test (n = 1000 permutations), which implements the nonparametric statistical testing (Maris & Oostenveld, 2007). The procedure uses a cluster analysis with a permutation test for calculating corrected p values. MI values arranged in a matrix were used as the input to the permutation test, and participants were used as observation samples (Figure 2A, red scale plots).
To assess the stability of PAC over the course of the EI and RI, we compared consecutive bins within the EI and RI. Similar to the analysis against baseline, we used a cluster permutation test to assess the whole range of computed MI values. We also employed a second approach that focused on the average MI values of α-β PAC and used a t test for repeated measures (ttest_rel as implemented in SciPy Python library).
α Peak-locking Procedure for a Graphical Representation of α-β PAC To visualize α-β PAC in the time domain, we aligned the peak of the α oscillation with the high-frequency β activity. Data were filtered with respect to the frequencies observed in the PAC analysis: As we observed a robust α-β PAC, we focused on illustrating how β power aligned with the phase of the α frequency. For this, we used a finite impulse response filter implemented in the pactools Python package (v0.1; Dupré la Tour et al., 2017). The α peaks were detected by localizing the maxima in the phase of the low-frequency signal. Time domain (α) and frequency domain (β) signals were then locked with respect to localized maxima in the phase of the α signal. The power of the β activity was assessed using a Hilbert transform to compute the instantaneous amplitude (the Hilbert analytical signal was multiplied with its conjugate, and the real part was extracted). The resulting plots are displayed in Figure 2C.

Localization of PAC Using a Jackknife Procedure
A specific PAC regime is described by its phase and frequency parameters. To determine the location of the PAC peaks, we used a jackknife procedure, which is well suited to discover small shifts in the parameters of interest. Previously, this technique was used to estimate small differences in the latencies of evoked brain components (Ulrich & Miller, 2001;Miller, Patterson, & Ulrich, 1998). Here, we used this procedure to tease apart differences in the PAC peak parameters between conditions. Per individual, we identified the five sensors that exhibited, after averaging, the largest α-β PAC. Each "jackknife sample" was obtained by removing a single sample and averaging the remaining samples (n − 1; Ulrich & Miller, 2001;Miller et al., 1998). The coordinates of amplitude frequency and phase frequency of the MI were obtained on the basis of jackknife samples. F values, obtained by running the ANOVA (AnovaRM as implemented in statsmodels Python library, v0.9) on jackknife-based samples, were corrected by dividing obtained F values by (n − 1) 2 .

Regression Analyses between Behavioral Precision and PAC
The sensor level analyses were performed using a linear regression and a model comparison to check whether factors other than PAC were needed to account for timing behavior. All statistical analyses were performed using R v3.3.2 statistical programming language (R Core Team, 2000).
Similar to the within-participant approach, MEG data were randomly split into four trial-based bins to provide a maximal number of within-participant samples for the regression analyses while preserving a minimal number of samples for the computation of the coefficients of variation (CVs; ∼10 samples per bin; van Belle, 2008). We again selected five sensors with the highest MI for α-β PAC. To maximize PAC contribution to the selected sensors, the sensors were selected per participant, condition (i.e., 2, 3, and 4 sec), and task (EI or RI) stage. This ensured that each sample entered in the regression resulted from the most prominent PAC contributions.
The inverse CV (invCV: mean duration production divided by the standard deviation in a set of approximately 10 trials) and the comodulograms were computed for each data bin. Each outcome served as a single observation in the regression model. We choose invCV as opposed to CV for its convenient positive association with precision, which facilitates interpretation of regression results in relation to α-β PAC. Thus, invCV captures the width of a distribution of time reproduction in a straightforward manner.
As we split the data per block and per individual, we used linear mixed-effects models (e.g., Gelman & Hill, 2007;Pinheiro & Bates, 2000) to account for sample dependency (multiple per-participant observations in the data). Linear mixed-effects models are regression models that model the data by taking into consideration multiple levels. Participants and trial-based bins were random effects in the model and were allowed to vary in their intercept. p Values were calculated based on a Type 3 ANOVA with Satterthwaite approximation of degrees of freedom, using the lmerTest package in R (Kuznetsova, Brockhoff, & Christensen, 2017). The mixed-effects model approach was combined with model comparisons that allow for selecting the best fitting model in a systematic manner.
Whenever a Shapiro-Wilk normality test indicated a deviation from normality, we transformed the data using the Lambert W function. The Lambert W function provides an explicit inverse transformation, which removes heavy tails from the observed data (Goerg, 2011(Goerg, , 2015. The data are transformed into a heavytailed form using log-likelihood decomposition; subsequently, the heavy-tailed form is transformed back into a Gaussian distribution. All transformations were performed using the LambertW R package.

α-β PAC Is Robustly Observed during Encoding and Reproduction of Temporal Intervals
In a previous work, we showed a robust α-β PAC in a temporal production task and found that incorporating EEG data, in addition to MEG recordings, did not qualitatively change the outcomes (Grabot et al., 2019). Thus, for simplicity, we will focus here on the MEG signals.
To quantify the presence of PAC during the encoding and reproduction of temporal intervals, we computed PAC over a broad spectrum of frequencies and collapsed the data across reproduced durations, participants, and sensors. The MI (Tort et al., 2009) quantifies PAC by estimating a deviation of the amplitude distribution with respect to a certain phase from the uniform distribution. We quantified the MI separately for the EI and RI to produce comodulograms. This allowed us to identify the peaks in the shuffled distribution of the MI. This analysis revealed significant peaks of α-β PAC (6-12 Hz for the carrier and 15-40 Hz for the modulated frequency) in both the EI and RI (see Figure 2A for z scores and t values). Whereas previous work identified an 8-to 12-Hz low-frequency carrier (Grabot et al., 2019), the significant lower PAC frequencies showed a lower center of gravity at 7-8 Hz (Figure 2A). Therefore, all analyses used an extended lower PAC frequency range (6-12 Hz). This difference could be because of a supine or seated position as body position significantly affects oscillatory brain activity (Spironelli, Busenello, & Angrilli, 2016).
We inspected the α-β PAC coupling regime by aligning β power to α oscillations in the time domain. An example of realigned time series for a single participant is provided in Figure 2B: For this participant, β power was highest at the trough of the α oscillation and lowest at the peak of the α oscillation. The comparison of each α-β PAC within each bin against the baseline interval converged with the permutation analysis of the MI and visual inspection of peak-locked data. Each bin during encoding and during reproduction showed a significant increase in the strength of α-β PAC as compared to the baseline (Figure 2A, right panels for EI and RI). The observation that α-β PAC during the EI significantly increased as compared to the baseline corroborated the prior working hypothesis that α-β PAC may play a role in timing (Grabot et al., 2019): Considering that α-β PAC was present when participants were not required to perform motor actions, α-β PAC may go beyond precision in motor execution and motor preparation. A significant increase of α-β PAC was also observed during reproduction.
As stronger couplings in the EI or RI could indicate predominant support of α-β PAC in one of the two stages of the task, we then contrasted the MI in the EI and RI using a cluster permutation test. We found no evidence of a significant α-β PAC difference between the EI and RI ( p > .1). Permutation tests were accompanied by Bayesian t test focused on α-β PAC (temporal Bin 1 in EI was tested against temporal Bin 1 in RI, etc.). None of the tested comparisons showed a Bayes factor exceeding 1, suggesting that α-β PAC in the EI and RI could equally contribute to the behavioral outcomes.

The Strength of α-β PAC Is Constant during Encoding and Reproduction of Temporal Intervals
Climbing neuronal activity measured as a gradual change of amplitude from onset to offset (e.g., contingent negative variation; Walter, Cooper, Aldridge, McCallum, & Winter, 1964) has been associated with dynamic processes related to duration perception (Boehm, van Maanen, Forstmann, & van Rijn, 2014;Kononowicz & van Rijn, 2014;Wiener et al., 2012;Macar et al., 1999). Along the same lines, we sought to investigate whether α-β PAC underwent gradual changes within the EI or the RI or whether it was stable within the timed interval. Specifically, if α-β PAC leverages the precision of duration (Grabot et al., 2019), α-β PAC should be stable within the RI; alternatively, if α-β PAC supports dynamical features of on-line timing-such as the gradual integration of temporal information-we should observe a covariation of α-β PAC with time, within a timed interval (EI or RI). To test this, we assessed the magnitude of α-β PAC separately during the EI and RI (Figure 3) using two approaches.
In the first approach, we focused on α-β PAC using collapsed values of five sensors with maximal MI per condition. We then performed a series of t test comparisons between consecutive temporal bins, within the EI, within the RI, and for all duration conditions, also using a Bayesian t test. None of the comparisons showed significant results (all ps > .1, all Bayes factors < 1). This series of comparison was not adjusted for multiple comparisons. In the second approach, we collapsed the data across all sensors and performed the permutation tests on a broad MI frequency spectrum, comparing consecutive bins in the EI and RI. Again, none of the comparisons showed significant results (all ps > .1).
Hence, across all durations, neither approach supported the notion of dynamical changes in the strength of α-β PAC during the encoding period or during the reproduction period. This suggested a relatively stable α-β PAC within the timed intervals, in line with the precision hypothesis, that is, the notion that α-β PAC leverages the precision of endogenous timing goal.
The Characteristics of α-β PAC Could Not Differentiate Durations, Encoding, or Reproduction Considering that three durations were tested, we questioned whether the spectral characteristics of α-β PAC varied as a function of duration, during encoding and reproduction. According to the "precision hypothesis" (Grabot et al., 2019), one would not expect significant shifts of the PAC to be involved in the timing of different durations.
To explore the characteristics of α-β PAC as a function of durations, we used the data depicted in Figure 4: For each time bin, participant, EI, RI, and duration, we quantified the coordinate of the maximum MI using a jackknife method (Figure 4). The obtained maximal amplitude of the high-frequency MI modulated by the Figure 4. Stable peak of α-β PAC during trials and across reproduced durations. The location of the maximal amplitude frequency did not vary as a function of the reproduced time interval or task stage. The trend for the phase frequency to decrease during the RI did not reach significance. The individual data points were smoothed using 2-D kernel density estimate as implemented in Seaborn Python package. maximal phase frequencies was entered into a 2 (Stage: EI, RI) × 3 (Durations: 2, 3, and 4 sec) repeatedmeasures ANOVA. No significant (main or interaction) effects for amplitude frequency and phase frequency were found (all ps > .1). Although Figure 4 suggests a possible small effect of low-frequency oscillation of PAC between the EI and RI, a statistical test did not confirm this observation, F(1, 16) = 1.7, p = .212. Hence, we Figure 5. The difference of α-β coupling strength between encoding and reproduction informs on the precision of internal timing. (A) α-β PAC difference between Bin 1 of the EI and Bin 1 of the RI was significantly correlated with the precision of temporal reproduction (invCV). The model outcomes are also plotted in B. The topographical plots at the bottom of A illustrate the probability of selecting a given sensor for the statistical analysis. Namely, a sensor that was picked for statistical analysis was labeled as 1; and 0, otherwise. Mean values of 0 and 1 vector for each sensor are plotted. (B) The fitted values of the precision are plotted against the precision values predicted by the mixed model, informed by the difference of α-β PAC strength. (C) The top and bottom rows illustrate two samples of the effect plotted in A and B. Comodulograms in the EI and RI (left and right) for two subsets of trials (rows) with high reproduction precision (top) and with low reproduction precision (bottom). The subset with a higher behavioral precision is shown on the top row and displays larger α-β PAC in the EI than in the RI. TP = time production. found no evidence in our data for a shift of PAC across the EI, RI, or duration.

α-β PAC and Timing Precision
According to the precision hypothesis, α-β PAC may reflect the precision with which an endogenous timing goal may be maintained (Grabot et al., 2019). In our current reproduction paradigm, the encoding was separate from the subsequent reproduction, yielding two possible stages at which timing precision may play a role. We thus reasoned that the α-β PAC in both the EI and RI may separately reflect the precision during encoding and reproduction, respectively. Alternatively, the overall behavioral precision may depend on the relative contributions of the precision in sensory encoding and motor reproduction.
To test this hypothesis, we used the coupling strength of the observed α-β PAC in the EI and RI as a function of behavioral precision measured by invCV (see Methods: Regression analyses between behavioral precision and PAC). Given the stability of α-β coupling strength within and across the different duration trials, we focused our regression analyses on the first and last bins when assessing the precision in the EI and on the first bin only when assessing the precision in the RI. We chose the first and last bins under the assumption that the beginning and end of the EI and the end of RI should represent the same cognitive steps irrespective of the target duration. Each observation was calculated using the average of five sensors with maximal α-β MI for each participant and each condition. A regression model using α-β PAC to predict invCV showed no significant effects in the EI alone (all ps > .05) or RI alone ( p > .1).
We then tested whether timing precision depended on the relative α-β PAC in the EI and RI: If the target duration was encoded with high precision during the EI (high α-β coupling strength), the likelihood of a precise reproduction should increase, yielding a larger PAC in the RI than EI; conversely, if the precision during encoding was poorly achieved (low α-β coupling strength), precise reproduction could only be achieved with higher α-β coupling strength.
Alternatively, similarity between the EI and RI could play a role in behavioral precision. Relative α-β PAC in the EI and RI can be considered from the perspective of similarity. Performance in a match/nonmatch type of task may depend on similarity between the activity during the encoding and probe trial (analogically the EI and RI in the current paradigm). Recent work showed that neural representations do not remain stable between encoding and probe, drifting toward less accurate states (Samrani, Bäckman, & Persson, 2019). In our paradigm, the similarity could be quantified as a relative α-β PAC in the EI and RI. This account suggests that, in time reproduction tasks, most precise trials would show equal strength of α-β PAC between the EI and RI. The α-β PAC ratio imbalance between the EI and RI should be associated with a decrease of behavioral precision. The similarity account predicts a U-shape association between the difference in α-β PAC during the EI and RI and behavioral precision.
To test these predictions, we used the difference between the strength in α-β PAC, indexed by the MI (MI α-β ), in the EI and RI [MI α-β (EI) − MI α-β (RI)] as a predictor for invCV. The regression model, using the difference in α-β MI between the initial bin in the EI and RI, showed that the strength of α-β oscillatory coupling significantly predicted the behavioral invCV, F(172) = 2.8, p = .006 ( Figure 5A and 5B). To ensure that the effect was stable over time, we tested the difference of α-β coupling strength between the last bin in the EI and the first bin in the RI: The regression model showed that the strength of α-β oscillatory coupling significantly predicted the behavioral invCV, F(179) = 2.5, p = .014 (Supplementary Figure 1A and 1B). The analysis, based on Akaike information criterion (AIC; Wagenmakers & Farrell, 2004), showed that both models containing the difference of α-β PAC strength between the EI and RI as predictors were justified as compared to the model including only the intercept (ΔAIC = −5.5, p = .006; ΔAIC = 4.0, p = .014, respectively).
In both models containing the difference between the EI and RI, including the α-β PAC strength as a predictor in addition to their interaction with duration was not warranted (ΔAIC = −1.7, p > .1; ΔAIC = −1.6, p > .1; the first and last bins, respectively), indicating that the association between precision and the EI-RI α-β PAC strength did not depend on the duration being reproduced.
Previous studies investigating PAC have indicated that the MI estimation could be confounded by the estimation of phase and power frequencies (e.g., Aru et al., 2015). To assess whether the association between PAC strength and behavioral invCV was exclusively driven by oscillatory coupling and not confounded by the power in α or β bands, we tested whether the inclusion of α or β power was justified in the model predicting invCV. Similar to α-β PAC, the α and β power were computed as a difference between the EI and RI. The inclusion of α and β power was not justified in the model predicting invCV for the first bin (ΔAIC = −0.7, p > .1; ΔAIC = −1.8, p > .1, respectively) and the last bin (ΔAIC = −1.1, p > .1; ΔAIC = −1.8, p > .1, respectively).

DISCUSSION
We investigated the role of PAC in the encoding and reproduction of time intervals. Our results provide evidence that the strength of α-β PAC leverages the precision of timing, replicating, and largely extending previous observations in a temporal production task (Grabot et al., 2019). Here, we specifically show an implication of α-β PAC during the encoding of a duration and during the reproduction of the encoded temporal interval. We also show the stability of this effect across three different durations, ranging from 2 to 4 sec. Crucially, we found that the difference in α-β PAC between the encoding and reproduction of temporal intervals was associated with the precision of reproduced intervals ( Figure 5C). This result suggested that behavioral precision not only is a function of α-β PAC strength in the reproduction stage but also incorporates the precision of evidence during the encoding of the time interval. With this study, we thus provide additional evidence that the strength of α-β PAC indexes the precision with which temporal information is maintained in the brain.
α Oscillations and PAC have been implicated in the maintenance of task-relevant information in working memory (Roux & Uhlhaas, 2014). In line with this view and with recent work (Grabot et al., 2019), α-β PAC was interpreted as reflecting the active maintenance (α) in the timing network of a neural code for duration (β), thereby controlling the precision of timing information. The precision hypothesis predicted an association between α-β PAC and behavioral precision to occur not only in the RI but also in any period during which the maintenance of temporal precision would be needed. The current results support this hypothesis as behavioral precision was best predicted by the α-β PAC difference between the EI and RI. More specifically, behavioral precision was a function of the strength of PAC during the encoding stage relative to the strength of PAC in the reproduction stage: Precise timing was supported by a high α-β PAC in the EI relative to a low α-β PAC in the RI. This seems to indicate that, if the temporal goal was encoded with a limited degree of precision in the EI, its maintenance in the RI would not lead to precise temporal reproduction: To the contrary, stronger PAC in the EI did support higher precision of reproduced durations. This imbalance suggests that α-β PAC primarily supports the encoding of a temporal target as opposed to the sole precision of motor commands during temporal production. In line with the precision hypothesis, α-β PAC appears to index the maintenance of higher-level temporal representations distinct from the control of motor commands.

α-β PAC Does Not Support the Motor Control Hypothesis
Although recent studies have reported β power to play a role in time estimation, the functional role of β power originally described in motor tasks (Pfurtscheller, Stancak, & Neuper, 1996; for a review, see Kilavik, Zaepffel, Brovelli, MacKay, & Riehle, 2013) was associated with the control of motor commands (Swann et al., 2011). Here, we considered the hypothesis that α-β PAC could support motor control, as opposed to precision of endogenous timing goals (Grabot et al., 2019).
We report several sources of evidence against the motor control hypothesis. First, the motor control hypothesis predicts that α-β PAC should only be observed in the RI. Considering there is still a possibility that the presentation of tones could elicit covert motor responses (Schnitzler, Salenius, Salmelin, Jousmäki, & Hari, 1997), a slightly more liberal version of the motor control hypothesis would suggest that α-β PAC may be stronger in the RI than EI. Despite plausible covert motor responses in the EI, motor control process should be more vivid in the RI, in which participants do execute motor action. We found no evidence in our data for a difference of PAC between the EI and RI. Second, according to the motor control hypothesis, the association between α-β PAC and behavioral precision would only be observed in the RI. Finally, if α-β PAC reflects motor accuracy, it seems improbable that high α-β PAC at encoding and low α-β PAC at retrieval should yield high timing precision.
Neither predictions of the motor control hypothesis were verified in our study: We found clear α-β PAC in the EI that did not differ from α-β PAC in the RI, and the difference in α-β PAC between the encoding and reproduction of temporal intervals was associated with the precision of reproduced intervals.
Taken together, these results strengthen the suggestion that α-β PAC in motor timing is not only a function of motor control precision. Instead, the obtained results are in line with the hypothesis that α-β coupling supports timing-related processes beyond motor control.

α-β PAC in Timing
To better understand the specificity of representations coded by α-β coupling, future studies could investigate α-β coupling in different timing contexts, for example, in explicit and implicit timing, in which the key distinction would be whether task instructions require participants to provide an explicit estimate of duration or not (Herbst & Obleser, 2017;Coull & Nobre, 2008). In current and previous work (Grabot et al., 2019), α-β PAC was found during explicit timing tasks, but such coupling may extend to implicit timing (e.g., foreperiod paradigm; Praamstra, 2010) or temporal expectation (Breska & Deouell, 2014;Nobre, Correa, & Coull, 2007). As previously suggested, the precision hypothesis could reflect a generic property of PAC strength indexing the temporal precision of information maintenance in brain networks (Grabot et al., 2019). There is emerging evidence that delta-beta PAC (θ-β PAC) is sensitive to the occurrence, or absence, of a predicted external event in time (Cravo, Rohenkohl, Wyart, & Nobre, 2011). θ-β PAC has been found to increase in response to predictable cues in temporal orienting (Mento, Astle, & Scerif, 2018). Together, these results suggest that PAC is involved in anticipatory processes in implicit timing. Future studies should explore the specific functionality of coupling regimes jointly in explicit and implicit timing and further disentangle its role in sensoryversus motor-based timing.
Cerebral Origins of α-β PAC Grabot et al. (2019) found large variations in the spatial position of sensors involved in α-β PAC. On average, the largest α-β PAC values were found in the motor and sensorimotor regions as well as in the frontal regions. The precision hypothesis was not limited to any particular brain region. Here, the main goal was to compare the strength of α-β PAC in the EI and RI. Therefore, we choose to compare the selected sensors that responded strongly in both the EI and RI on a per-individual basis, meaning that, within a participant, a set of selected sensors could differ. Although less commonly used, that approach guaranteed that sensors representing the largest PAC were compared between the EI and RI, in line with the objective.
The sensor selection approach suggests that α-β PAC was considered as a local property of brain activity. However, α oscillations are typically seen as a global property of brain activity (Nunez, Wingeier, & Silberstein, 2001). To reconcile the global and local views on the cerebral basis of α-β PAC, Grabot et al. (2019) showed that an individual α peak frequency correlated with the individual's α peak frequency observed in α-β PAC. We suggest that the phase of α oscillations may provide an optimal window for the reactivation of β activity. α-β PAC could originate as the interplay between local and global brain networks.

Conclusion
In this work, we report that α-β PAC is present during the encoding and reproduction of time intervals and that timing precision results from the balance between the strength of α-β PAC during the encoding and reproduction of intervals. We therefore suggest that α-β PAC reflects the precision with which representations are maintained in the brain.