Within-cycle instantaneous frequency profiles report oscillatory waveform dynamics

Non-sinusoidal waveform is emerging as an important feature of neuronal oscillations. However, the role of single cycle shape dynamics in rapidly unfolding brain activity remains unclear. Here, we develop an analytical framework that isolates oscillatory signals from time-series using masked Empirical Mode Decomposition to quantify dynamical changes in the shape of individual cycles (along with amplitude, frequency and phase) using instantaneous frequency. We show how phase-alignment, a process of projecting cycles into a regularly sampled phase-grid space, makes it possible to compare cycles of different durations and shapes. ‘Normalised shapes’ can then be constructed with high temporal detail whilst accounting for differences in both duration and amplitude. We find that the instantaneous frequency tracks non-sinusoidal shapes in both simulated and real data. Notably, in local field potential recordings of mouse hippocampal CA1, we find that theta oscillations have a stereotyped slow-descending slope in the cycle-wise average, yet exhibiting high variability on a cycle-by-cycle basis. We show how Principal Components Analysis allows identification of motifs of theta cycle waveform that have distinct associations to cycle amplitude, cycle duration and animal movement speed. By allowing investigation into oscillation shape at high temporal resolution, this analytical framework will open new lines of enquiry into how neuronal oscillations support moment-by-moment information processing and integration in brain networks.


1: Introduction
Frequency, phase and amplitude have long been reported as important features of neuronal 3 oscillations with behavioural and electrophysiological relevance. Furthermore, neuronal 4 oscillations show non-sinusoidal waveform shapes that span a wide range of spatial and 5 temporal scales . Whilst waveform shape is emerging as a fourth 6 relevant feature of neuronal oscillations, many theories of neuronal oscillations currently 7 assume sinusoidal waveforms. This might be due to the fact that characterising and 8 quantifying non-sinusoidal waveforms remains a substantial analytic challenge (Amzica and 9 Steriade, 1998; . To uncover the role of waveform dynamics in 10 rapidly unfolding brain activity, there is a growing need for novel analysis methods that are 11 able to characterise a wide range of waveform shape features at the single-cycle level. 12 13 Waveform shape related parameters, such as skewness or asymmetry, can be estimated from 14 higher order Fourier spectra such as the bispectrum or bicoherence (Bartz et  2014). Whilst this approach is tractable on single cycles, the extracted features must be 20 defined a priori and are limited to the resolution of the selected cycle control points such as 21 the extrema and zero-crossings. 22 23 The temporal dynamics in oscillatory frequency can be quantified for a given waveform by 24 its instantaneous frequency computed from the differential of the signal's instantaneous phase 25 (Boashash, 1992;Huang et al., 2009). Such instantaneous frequency estimates have been 26 used previously in electrophysiology to explore dynamics in oscillatory peak frequency at 27 high temporal resolution (Cohen, 2014;Liang et al., 2005;Nelli et al., 2017;Rudrauf et al., 28 2006). Crucially, any non-sinusoidal waveform features in an oscillation will lead to within-29 cycle instantaneous frequency modulations in which the frequency of an oscillation changes 30 from moment-to-moment within a single cycle (Huang et al., 1998). The degree of non-31 linearity of an oscillation is related to the total amount of within-cycle frequency modulation 32 (Huang et  We introduce a novel approach which creates smooth waveform shape profiles that describe 35 non-sinusoidal features in single cycles with high temporal detail. To this end, we first 36 operationalise waveform shape as the profile of instantaneous frequency across the cycle's 37 instantaneous phase. We then identify when and how an ongoing cycle deviates from a 38 sinusoidal waveform by identifying points in the cycle where instantaneous frequency departs 39 from a flat profile. For example, a cycle with a wide peak has a relatively low instantaneous 40 frequency around the peak, and a cycle with a fast-ascending edge will have a relatively high 41 instantaneous frequency between the trough and peak. In order to allow between cycle 42 comparisons, we also need to account for how different cycles of an oscillation will play out 43 at different speeds, leading to differences in extrema timing and overall duration. To 44 overcome these problems, we introduce the process of 'phase-alignment' which reregisters 45 the instantaneous frequency profiles onto a regularly sampled set of points in phase. To 46 obtain the instantaneous phase time course of each cycle, we use the Empirical Mode 47 Decomposition (EMD). EMD decomposes the time-series of interest into its oscillatory 48 modes (Intrinsic Mode Functions; IMFs) that retain the non-stationary and non-linear signal 49 features. 50 51 We outline and validate our novel approach in simulated data before applying it to theta-band 52 oscillations recorded in the local field potentials (LFPs) of the mouse hippocampal CA1 53 during active exploratory behaviour. The hippocampal theta rhythm has a characteristic non-54 sinusoidal waveform shape (Belluscio et al., 2012;Buzsáki et al., 1985;Siapas et al., 2005) 55 that is modulated by movement (Sheremet et al., 2016) and changes in sleep or drug states 56 (Buzsáki et al., 1985). Using EMD to identify the theta rhythm, we show that phase-aligned 57 instantaneous frequency is able to robustly characterise a continuous waveform shape profile 58 for hippocampal theta. The results confirm the stereotyped fast-ascending and slow-59 descending shape in the cycle-wise average. Additionally, it reveals but with high amounts of 60 shape variation across single cycles of theta, which we describe using a set of data-driven 61 shape 'motifs'. Finally, we find that cycle-level shape motifs have differential associations 62 with theta amplitude, theta cycle duration and mouse movement speed. The next stage is to segment the IMFs into their constituent cycles and identify which cycles 137 will be included in further analysis. The start and end of theta cycles are located by the 138 differentials greater than six in the phase. The start and end point of cycles in this paper is the 139 ascending zero-crossing as this occurs at the point where the phase time-course wraps. Once 140 identified, some cycles will be 'bad' in the sense that the oscillation captured by the IMF is 141 not well represented, e.g., because the rhythm is not present over that time period or it is 142 poorly estimated, and will be excluded from subsequent analyses. This is important for 143 instantaneous frequency analyses as the differentiation step (see section 2.2.1) can be very 144 noise sensitive. Included cycles are identified from the wrapped instantaneous phase time-145 course of the IMF containing the oscillation to be analysed. As the instantaneous phase 146 computation via the Hilbert transform returns a value for every sample regardless of whether 147 a prominent rhythm is present, only 'good' theta cycles are retained for further analysis. A 148 good cycle is defined as having a phase with a strictly positive differential (i.e., no phase 149 reversals) that starts with a value 0≤x≤π/24 and end within 2π-π/24≤x≤2π and 4 control 150 points (peak, trough, ascending edge and descending edge). 151 152

2.2.3: Control Point analysis 153 154
One approach to quantify oscillatory waveforms is to base analyses around parts of the cycle 155 which can be straightforwardly identified in all cycles (Belluscio et al., 2012;Cole and 156 Voytek, 2019). For example, the peak, trough, ascending zero-crossing and descending zero-157 crossings of a cycle can define a set of control points describing the relative timing of 158 inflections and mid-points in a cycle. In the present analyses, the extrema (peaks and troughs) 159 are detected by finding zero-crossings in the differential of an IMF. This initial extrema 160 estimate is limited by the sampling frequency of the dataset and is refined using parabolic 161 interpolation (Rato et al., 2008 We present an alternative approach to control points that ensures that entire waveform 169 profiles can be combined and/or compared across cycles despite cycle-by-cycle differences in 170 progression rate and overall duration. To compare waveforms across cycles that play out at 171 different speeds, we use phase-alignment to register cycles onto a common grid. Phase 172 alignment is performed on the instantaneous phase of a cycle and a measure of interest, such 173 as the instantaneous frequency, which is observed at the same time-intervals. A linear one-174 dimensional interpolation function is fitted between the instantaneous phase (as x values) and 175 the instantaneous frequency (as y values). The interpolation function is evaluated on a 176 template set of instantaneous phase values with a linear spacing between 0 and 2pi, if any 177 points in the template fall outside the fitted range, the interpolator returns an extrapolated 178 value. This interpolated version of instantaneous frequency is then directly comparable across 179 cycles as each point in the phase will occur at the same time. We compute phase-alignment 180 using a linear interpolation across 48 fixed points across the zero to 2pi phase range. 181 182 Once an instantaneous frequency profile has been phase aligned, we can visualise a 183 normalised waveform by projecting the frequency content back to a phase-time course. This 184 is achieved by re-normalising the instantaneous frequency from Hertz back to radians in 185 order to create a profile of successive phase differences. The phase time-course is then 186 reconstructed from the cumulative summation of these phase differences. Where ∅ 6 is the uniform phase-grid used in phase-alignment and 444 is the phase-aligned 199 instantaneous frequency. This is similar to the mean-vector approach to computing phase-200 amplitude-coupling (Canolty et al., 2006). The mean-vector of a sinusoidal cycle will be zero 201 whilst a non-sinusoidal cycle will return a complex value whose angle indicates which phase 202 has the highest instantaneous frequency and the magnitude indicates the extent of the 203 frequency modulation through the cycle. This method provides a straightforward summary 204 but is only sensitive to unimodal deviations from a flat instantaneous frequency profile. The component motif matrix defines the axes of variability in waveform shape across cycles. 217 The shapes captured along each of these axes are visualised by defining a set of PC scores 218 containing the maximum or minimum observed score for the PC in question and zeros for all 219 others. These scores can be projected back into the original data space to provide exemplar 220 instantaneous frequency profiles for both extremes of the PC axes. Finally, these exemplar IF 221 profiles can be projected back into the time-domain to generate a normalised waveform that 222 preserves the shape depicted in the exemplar IF profiles. 223 224

226
2.3.1: Schematic cycle generation 227 228 To illustrate the relationship between waveform shape, instantaneous phase and frequency, a 229 set of noise free oscillations were generated. First, a linearly progressing phase time-course is 230 generated and sinusoid is created by taking a sine-transform of this wrapped phase. Different 231 non-sinusoidal cycles are generated by modulating the unwrapped phase time-course by sine 232 and cosine waves at different phases and frequencies. The cycles with extrema and edge 233 asymmetry are generated by modulating the phase with a 1 Hz sine or cosine respectively. 234 A more realistic noisy simulation was used for the results in figures 3, 4 and 5. A simulated 243 oscillation at 12 Hz was generated using an autoregressive oscillator with the following 244 transfer function: 245 246 Where is the angular frequency of the oscillator (in rads/sec) and is the magnitude of the 249 roots of the polynomial (0 < r < 1). For this simulation, we computed for r=0.95 and 250 equivalent to 12 Hz and used its parameters to filter (a forwards and backwards filter using 251 scipy.signal.filtfilt) white noise. This generates a noisy sinusoidal oscillation which contains 252 random dynamics in the frequency and amplitude of each oscillatory cycle. Sixty seconds of 253 data were generated at 512 Hz. 254 255 This simulated oscillation was then modulated by one of two equations defined in equations 256 50.24 and 50.25 in section 50-6 in Volume 1 of Feynman's Lectures of Physics (Feynman et  257 al., 2011). The first equation defines a linear system that scales the signal by a constant 258 leaving the waveform shape unchanged. 259 260 Whilst the second equation defines a nonlinear system that includes a term inducing a change 263 in waveform shape as well as scaling the signal. 264 265 The nonlinearity in the second equation makes the peak of the oscillations shorter and widens 268 the trough. To ensure that the detected theta cycles are physiologically interpretable theta activity, we 341 identified cycles in each recording during times where the speed of movement of the mouse 342 was greater than 1 cm/second. As faster movement is associated with stronger theta 343 oscillations, this restriction increases the probability that our detected cycles represent 344 physiologically interpretable theta events. We additionally restricted analyses to cycles in 345 IMF-6 where cycle duration corresponded to 4-11 Hz frequency range (i.e., 312 and 113 346 samples respectively) and cycle amplitude was above the bottom 10% of the amplitude 347 distribution. Finally, cycles that failed the cycle inclusion checks outlined in section 2.2.2 348 were removed from analysis at this point ( Figure 1D). 349 350

2.4.3: Cycle comparisons 351 352
We computed the temporally aligned instantaneous frequency profile, phase-aligned 353 instantaneous frequency profile and normalised waveform for each included theta cycle. The 354 average waveform shape within each dataset was estimated from the averaged phase aligned 355 instantaneous frequency and a group average constructed from the mean of the six individual 356 runs. Variability in waveform shape across single cycles in the group data is summarised 357 using the instantaneous frequency mean vector (section 2.2.5) and visualised as a distribution 358 in the complex plane in which the x-axis represents asymmetry between ascending and 359 descending edge frequency and the y-axis represents asymmetry between peak and tough 360 frequency. For comparison, we also identified the control points from each cycle of the theta 361 IMF and constructed the peak-to-trough and ascending-to-descending duration ratios (Cole 362 and Voytek, 2019). 363 364

2.4.4: Waveform Motifs and Relation to behaviour 365 366
We next look to explore the waveform shapes which are present in the phase-aligned 367 instantaneous frequency values. We use PCA (section 2.2.6) to identify the data-driven set of 368 shape components that explain the most variance in the shape of theta cycles in this dataset. 369 The first four principal components explaining 95% of variance defined our four shape motifs 370 and were retained for further analysis. The reproducibility of the PCA is validated across 500 371 split-half iterations assessing the proportion of variance explained by each PC and the 372 correspondence between the component shape in the two halves (Supplemental section 8.3; 373 figure S3). 374 375 The relationship between the shape motifs and a set of three cycle covariates (i.e., cycle 376 amplitude, cycle duration and mouse movement speed) was quantified using a General Linear 377 Model (

414
We next use simulations to illustrate how instantaneous frequency analyses can be conducted 415 on a noisy signal with a dynamic 12 Hz oscillation modified by a non-linearity which widens 416 the trough of each cycle (see section 2.3.2). This oscillation was isolated from the noisy 417 background using mask sift ( Figure 3B) before the Hilbert transform was used to compute the 418 instantaneous phase time-course ( Figure 3C). It is evident that the phase time-courses do not 419 progress linearly through all cycles; these deviations from monotonic phase progression are 420 quantified in the instantaneous frequency time-course ( Figure 3D). Instantaneous frequency 421 sweeps within a single cycle reflect the non-sinusoidal shapes of the time-domain waveforms. 422 For this simulation, the instantaneous frequency tends to be higher during the first half of the 423 cycle and lower in the second half, reflecting the non-linearity that shortens the peak and 424 widens the trough of these oscillations. 425 426 The HHT of the simulated signal ( Figure 3E) retains the high time-frequency resolution of 427 the instantaneous frequency time-course allowing within cycle frequency dynamics to be 428 visible. In contrast, while a standard 5-cycle Morlet wavelet transform identifies similar 429 power dynamics, variability in frequency within single cycles are not resolved ( Figure 3F). A 430 further disadvantage is that the non-sinusoidal waveform shape of this simulation introduces 431 a 24Hz harmonic component into the wavelet transform. 432 shown in Figure 4A. The ratios of peak-to-trough duration and ascending-to-descending time 453 of these cycles suggests longer troughs and shorter peaks, whilst the ascending and 454 descending portions of the cycle are approximately equal in duration ( Figure 4B). 455 456 An alternative approach for comparing cycles is to align the instantaneous frequency profiles 457 to one of the control points. For example, we aligned the 500 cycles to the ascending zero-458 crossing and computed their time-locked average ( Figure 4C). The time-locked instantaneous 459 frequency profile of these cycles is not flat, reflecting the presence of non-sinusoidal shape in 460 this simulated signal. However, the precise type of non-sinusoidal shape is ambiguous from 461 this average, due to variability in the location of different waveform features within single 462 cycles. In this case, the instantaneous frequency is highest around 10 samples after the 463 ascending zero-crossing; however, this time-lag might correspond to different points in the 464 waveform in different cycles. In addition, variability in the duration of cycles means that, 465 after a certain point, different numbers of cycles contribute to the average, making the 466 estimate unstable. 467 468 Here we present an approach that overcomes these shortcomings. In brief, phase-alignment 469 removes this ambiguity by visualising the instantaneous frequency of a cycle across a fixed 470 grid of points along its phase (see section 2.2.3). For instance, an oscillatory peak is 471 normalised to occur at the same phase value irrespective of the cycle's duration or shape. 472 This corresponds always to exactly one quarter of the phase of each cycle, but not necessarily 473 to one quarter of the duration of each cycle. By aligning the instantaneous frequency to the 474 phase, we remove the temporal distortions caused by varying shapes and cycle durations, and 475 express the shape with the phase-aligned instantaneous frequency values. The phase-aligned 476 instantaneous frequency of the simulated cycles ( Figure 4D) now unambiguously shows the 477 increased frequency around the peak of the 12Hz oscillation and decreased frequency around 478 the trough. The average across the phase-aligned cycles is then a smooth representation of the 479 shape of the entire cycle. Simulated data illustrating how waveform shape can be quantified using control point ratios, 484 instantaneous frequency and phase-alignment. 485 A : The durations between successive control points for each simulated cycle. 486 B : Top: distributions of peak-to-trough and ascent-to-descent durations. Bottom: The peak-487 to-trough and ascent-to-descent ratios for each cycle. Comparisons between sets of cycles is straightforward once waveform shape has been 493 estimated from instantaneous frequency and normalised through phase-alignment. This is 494 illustrated by contrasting the shape of a noisy a 12-Hz oscillation modulated by either a linear 495 or non-linear system (section 2.3.2; Figure 5A). The non-linear system has a wide-trough 496 shape whilst the linear system has a sinusoidal waveform ( Figure 5B & C). The average 497 phase-aligned instantaneous frequency values for the linear system correspond to a flat line at 498 12 Hz throughout the cycle. In contrast the non-linear system has an increased instantaneous 499 frequency in the first half of the cycle and a decreased frequency in the second half ( Figure  500 5D). We can compare the average instantaneous frequency profiles from both systems 501 ( Figure 5E) and compute conventional t-statistics ( Figure 5F) to quantify any differences in 502 waveform shape. In this simulation, we find that the non-linear system creates oscillations 503 with increased frequency peaks and decreased frequency troughs. 504 505 Finally, the phase-aligned average frequency profiles can be projected back into a 506 'normalised waveform' to more intuitively visualise the type of non-sinusoidal distortions. 507 These normalised waveforms have a constant duration and an amplitude of one, but retain 508 any distortions in waveform shape quantified in the instantaneous frequency profile. For this 509 simulation, the normalised waveform reveals the pinched peak and the widened generated by 510 the non-linear system ( Figure 5G). 511

524
LFP data recorded from the mouse hippocampus were analysed to explore the utility of 525 phase-aligned instantaneous frequency as a measure of waveform shape. Figure 6A shows a 526 three-second LFP recording from the pyramidal layer of the mouse dorsal CA1 (black) 527 overlaid with the EMD-extracted theta IMF (red, Figure 6B shows all IMFs). In this case, the 528 theta oscillation was isolated into IMF 6 with minimal disruption to its amplitude or 529 waveform shape dynamics. Many of the theta cycles within this window have prominent non-530 sinusoidal waveform shapes, which are qualitatively visible in both the raw data trace 531 (Buzsáki et al., 1986(Buzsáki et al., , 1985 and the EMD-extracted theta IMF. Importantly, the oscillatory 532 waveform shape varies between successive cycles, though the amplitude and duration of the 533 theta cycles are relatively consistent. 534 535 The instantaneous phase ( Figure 6C) and instantaneous frequency ( Figure 6D) were 536 computed from the theta oscillation in IMF-6. As with the simulation analysis, any within-537 cycle dynamics in the instantaneous frequency naturally represent the waveform of each 538 cycle. This was summarised with the standard deviation of instantaneous frequency values 539 within each cycle ( Figure 6E). As an illustration, cycles 5, 9 and 11 have relatively sinusoidal 540 shapes with flat instantaneous frequency profiles and low frequency variability. In contrast, 541 cycle 13 is relatively non-sinusoidal with a dynamic instantaneous frequency profile and high 542 frequency variability. The HHT provides a time-frequency description with sufficient 543 resolution to depict these within cycle instantaneous frequency sweeps ( Figure 6F). In 544 contrast, a 5-cycle wavelet transform of the same data was not able to resolve these dynamics 545 ( Figure 6G). 546 547 Looking at individual cycles illustrates how instantaneous frequency can characterise 548 waveform shape (Figure 7). Frequency increases and decreases correspond to slowing down 549 and speeding up of the cycle as its waveform shows non-sinusoidal behaviour. It is evident 550 that there are many observed shape profiles. For instance, the cycles labelled as i and iii in 551 figure 7 had frequencies that dip during the centre of the cycle, indicating an elongated, low 552 frequency descending edge. Cycle v had slowest frequency around -pi/2 corresponding to a 553 wide peak. In contrast, cycle vi had relatively high frequency around -pi/2 and lower 554 frequency at +pi/2 leading to a short, pinched peak and an elongated trough. Overall, the 555 phase-aligned instantaneous frequency profiles and normalised waveforms provide a rich 556 description of oscillatory waveform, despite wide variability in cycle amplitude, duration and 557 shape. 558

578
We next computed waveform shape from around 1500 hippocampal theta cycles from a 579 single recording using three different methods: control-point ratios (see section 2.2.6), 580 control-point locking, and phase-alignment (see section 2.2.3). The durations between 581 specified control points were computed for each cycle ( Figure 8A) and the peak-to-trough 582 ratio and the ascent-to-descent ratio computed ( Figure 8B). The peak-to-trough ratios are 583 evenly distributed around zero, whereas there is a bias in the ascent-to-descent ratios 584 suggesting that the descending edge of theta is longer than the ascending edge. The 585 instantaneous frequency profiles locked to the ascending zero-crossing show a wide variety 586 of shapes with a group average tendency for frequency to start around 9 Hz and to decrease 587 through the duration of the cycle ( Figure 8C). As described above, this average effect is 588 challenging to interpret due to within-cycle variability in the timing of cycle features and 589 between cycle variability in total cycle duration. Our proposed phase-aligned instantaneous 590 frequency profiles ( Figure 8D) resolve these ambiguities. This shows that theta cycle 591 instantaneous frequency in this single recording starts around 9 Hz at the ascending zero-592 crossing, decreasing to around 8.1 Hz at the descending zero-crossing, before increasing 593 again to 9 Hz at the end of the cycle. This is consistent with a fast-ascending and slow 594 descending cycle shape revealed by the control point analysis and in previous literature. The 595 phase-aligned instantaneous frequency approach is able to show this effect as a continuous 596 shape profile for single cycles, which can be straightforwardly compared at the group level. 597 instantaneous frequency profile for each cycle. 606 3.5: Theta has a stereotyped asymmetric shape with wide variability over cycles 607 608 We next summarised the average waveform across theta cycles from six recordings taken 609 from three mice. The average phase-aligned instantaneous-frequency profile is computed for 610 each recording and for the whole dataset. The overall group-level average waveform had a 611 cosine-type profile centred around an average instantaneous frequency of ~8.6 Hz ( Figure  612 9A, average in black and individual recording sessions in grey). On average, the 613 instantaneous frequency peaked within the cycle around 9 Hz at the ascending zero-crossing 614 and drops to just below 8.4 Hz between the peak and descending zero-crossing. These results 615 are consistent with previous studies showing an asymmetry between the fast-rising and slow-616 decaying halves of a theta cycle (Belluscio et al., 2012;Buzsáki et al., 1986;Cole and 617 Voytek, 2019). All six recordings across three animals showed a shape with a maximum 618 frequency around the ascending zero-crossing and a minimum on the descending edge, 619 though there was some variability in whether the lowest frequency was closer to the peak or 620 trough. 621 622 To visualise the variability in waveform shape across cycles and recording sessions, we 623 performed a complementary analysis using the instantaneous frequency mean-vector to see 624 the distribution of single-cycle waveforms across a simplified 2-dimensional shape-space 625 ( Figure 9B; section 2.2.4). The distribution had a non-zero mean on the x-axis for all 626 recordings, indicating that the highest frequencies in a cycle are typically at the ascending 627 edge, consistent with the average in Figure 9A and with previous literature on the theta cycle 628 (Belluscio et al., 2012). Though the overall mean shift in the distribution of cycles is robust, 629 there is substantial cycle-to-cycle variability indicated by the width of the distribution. 630

642
To further describe the variability in waveform shape across cycles and characterise its 643 relation to movement speed, theta amplitude and theta cycle duration we identify a set of 644 waveform shape 'motifs' using PCA. The PC values define a set of shape motifs that are each 645 expressed to different degrees in the observed cycles (see section 2.2.5 and 2.4.4). The first 646 four components describing 96% of variance are retained for further analysis (Figure 10). 647 The waveform shape represented by each PC motif is summarised by the normalised 648 waveforms ( Figure 10A). These normalised waveforms are computed from the instantaneous 649 frequency profiles of the PC component-vectors ( Figure 10B) projected onto the extreme 650 ends of the PC score distribution ( Figure 10C). The shape of each individual cycle can then 651 be described by a set of four PC scores, relating to the amount of each component which it 652 contains. 653 654 PC-1 (63.31% of variance) describes a continuum of shape from a sharp peak and wide 655 trough through to a wide peak and sharp trough. This shape is similar to the y-axis in the 656 mean vector distribution in Figure 9B. In contrast, PC-2 (22.56% of variance) describes 657 shapes ranging between an elongated ascending edge and an elongated descending edge, 658 similar to the x-axis of Figure 9B. The remaining components describe more complex shapes 659 with relatively small contributions to the variance explained. PC-3 (6.86% of variance) 660 captures shapes with a left or right 'tilt' around their extrema and PC-4 (3.33% of variance) 661 describes shapes with a sharper or flatter curvature around the extrema. 662 663 The control-point-based ascending-to-descending ratio and peak-to-trough ratio are computed 664 for each cycle. For each PC, these values are partitioned into cycles with positive or negative 665 PC scores (relating to distinct ends of the shape continuum for that component) and their 666 distributions plotted in Figure 10D. The peak-to-trough ratios are clearly separated in the two 667 ends of PC-1 whilst the ascending-to-descending ratios are similar for cycles with a positive 668 or negative score in PC-1. This is consistent with the normalised waveforms summarising 669 PC-1 in Figure 10A. PC-2 also shows the expected separation of ascending-to-descending 670 ratios by PC score, whilst the peak-to-trough ratios are unchanged. Whilst PC-3 and PC-4 671 describe around 10% of shape variability, they are not characterised by the control point 672 analyses. Neither peak-to-trough ratios nor ascending-to-descending ratios are changed by PC 673 score for PC-3 or PC-4. These shape profiles are robustly identified by the phase-aligned 674 instantaneous frequency method but are not distinguished by these control point-based 675 metrics as the shape distortions in PC-3 and PC-4 occur between the four specified control 676 points. 677 678 A general linear model was used to quantify the relationship between the different shape 679 motifs and theta amplitude, theta duration and mouse movement speed. This regression is 680 computed separately for each PC and the resulting parameter estimates converted into t-681 statistics. PC-1 codes for changes in average instantaneous frequency across the cycle with a 682 small shape distortion around the descending edge. This PC has a strong relationship with 683 cycle duration and amplitude but no significant covariation with movement speed. Longer 684 and higher amplitude cycles tend to have more positive scores in PC-1 relating to wide peak 685 shapes. PC-2 has a strong relationship with duration and movement speed. Specifically, 686 cycles with elongated descending edges have longer cycle durations and are more likely to 687 occur during fast animal movement. PC-3 shows significant covariance with cycle amplitude. 688 High amplitude cycles tend to have shapes in which instantaneous frequency is relatively 689 high just before the peak or trough. Finally, PC-4 varies strongly with duration and weakly 690 with movement speed. Cycles with flatter curvatures around the extrema have longer cycle 691 durations and are less likely to occur during faster animal movement. 692 Figure 10 Non-sinusoidal waveforms are often visible by eye in raw LFP traces of electrophysiological 712 datasets, yet discovering and quantifying these non-sinusoidal and non-linear features present 713 substantial analytic challenges. We utilise within-cycle variability in instantaneous frequency 714 to describe distortions in waveform ). Further, we introduce phase-715 alignment as a solution to comparing full-resolution waveforms between cycles of different 716 durations. Taken together, we establish that the phase-aligned instantaneous frequency profile 717 of an oscillation provides a flexible framework for complete characterisation of oscillatory 718 waveform shape. We demonstrate the utility of this approach by applying it to simulated data 719 and LFP recordings of theta oscillations of behaving mice. 720 721 In real data, we observed that theta oscillations have, on average, a fast-ascending and slow 722 descending waveform, in line with previous reports (Belluscio et al., 2012;Buzsáki et al., 723 1986Buzsáki et al., 723 , 1985Cole and Voytek, 2019). Though this average shape is robust across many 724 cycles, recording sessions and animals; the shape of individual cycles is highly variable. We 725 characterise this variability using PCA to identify a range of shape components, or 'shape 726 motifs', which maximally explain the variability in the dataset. The first two PCs quantify the 727 relative durations of the peak and trough (PC-1) and the ascending to descending edge (PC-728 2). These PCs broadly map onto the features described by the peak-to-trough and ascending-729 to-descending control point ratios. We show that these theta shape PCs have distinct patterns 730 of covariation with movement speed, theta amplitude and theta cycle duration. Critically, we 731 show that though PC-2 describes less variability overall, it most clearly co-varies with 732 movement speed. 733 734 PC-3 and PC-4 capture more complex waveform shapes. We show that the curvature around 735 the extrema of the waveform shape (PC-4) is wider in theta cycles occurring during faster 736 animal movement. This shape is naturally described by instantaneous frequency but not 737 visible to standard ascending-to-descending and peak-to-trough control-point ratios. More 738 generally, if the waveform shape of interest is known a priori, it is possible to construct 739 specific control point-based measures so that the waveform shape can be identified. For 740 instance, waveform sharpness can be explored by looking at the differential between the 741 extrema and the samples 5ms before and after . However, in real data, we 742 may not know the waveform shape of interest a priori, implying that many separate metrics 743 may need to be computed for each cycle. In contrast, the phase-aligned instantaneous 744 frequency can quantify any waveform shape as a within-cycle instantaneous frequency sweep 745 without pre-specifying the features which may be of interest. 746 747 The present results demonstrate that single-cycle dynamics in oscillations can be 748 meaningfully estimated using phase-aligned instantaneous frequency, and that specific shape 749 motifs are differentially related to the wider electrophysiological (theta amplitude and 750 duration) and behavioural (movement speed) context. Future models of theta function may 751 consider these dynamics in waveform shape that deviate from a canonical sinusoidal theta 752 template. Given that many sub-processes occur preferentially at different parts of the theta 753 cycle (Klausberger and Somogyi, 2008), we hypothesise that shape distortion may indicate or 754 reflect a change in the underlying theta-phase-nested sub-processes. 755 756 The outlined approach requires that each cycle is smooth in both its waveform and phase 757 profiles, as any jumps or discontinuities will lead to noisy or even negative instantaneous 758 frequency estimates. If the cycle is smooth, we can characterise very large distortions in 759 waveform shape as within-cycle dynamics in instantaneous frequency. Finally, we assume 760 that the features being analysed are well described as oscillations. If the features are non-761 sinusoidal and non-oscillatory, such as spiking activity, then descriptions using the language 762 of frequency may not be appropriate. With these improvements and caveats in hand, this 763 approach is readily generalisable to other datasets and provides a flexible framework for 764 investigating waveform shape oscillating systems. 765 766 In conclusion, the full-cycle waveform of single cycles of hippocampal theta can be 767 quantified and explored with phase-aligned instantaneous frequency. We use this approach to 768 confirm the characteristic fast-ascending waveform of theta oscillations; and to additionally 769 reveal that this is highly variable on the single-cycle level. Moreover, we are able to link this 770 variability with behavioural and electrophysiological states, suggesting that waveform shape 771 is a relevant feature of neuronal oscillations alongside frequency, phase and amplitude. 772 Finally, whilst we have illustrated this approach with hippocampal theta oscillations, it is 773 likely that this methodology will readily generalise to neuronal oscillation in other brain 774 regions, frequency bands and contexts. 775