Complex interplay between spectral harmonicity and different types of cross-frequency couplings in nonlinear oscillators and biologically plausible neural network models

Damián Dellavale, Osvaldo Matías Velarde, Germán Mato, and Eugenio Urdapilleta

Phys. Rev. E 102, 062401 – Published 1 December 2020

Abstract

Cross-frequency coupling (CFC) refers to the nonlinear interaction between oscillations in different frequency bands, and it is a rather ubiquitous phenomenon that has been observed in a variety of physical and biophysical systems. In particular, the coupling between the phase of slow oscillations and the amplitude of fast oscillations, referred as phase-amplitude coupling (PAC), has been intensively explored in the brain activity recorded from animals and humans. However, the interpretation of these CFC patterns remains challenging since harmonic spectral correlations characterizing nonsinusoidal oscillatory dynamics can act as a confounding factor. Specialized signal processing techniques are proposed to address the complex interplay between spectral harmonicity and different types of CFC, not restricted only to PAC. For this, we provide an in-depth characterization of the time locked index (TLI) as a tool aimed to efficiently quantify the harmonic content of noisy time series. It is shown that the proposed TLI measure is more robust and outperforms traditional phase coherence metrics (e.g., phase locking value, pairwise phase consistency) in several aspects. We found that a nonlinear oscillator under the effect of additive noise can produce spurious CFC with low spectral harmonic content. On the other hand, two coupled oscillatory dynamics with independent fundamental frequencies can produce true CFC with high spectral harmonic content via a rectification mechanism or other post-interaction nonlinear processing mechanisms. These results reveal a complex interplay between CFC and harmonicity emerging in the dynamics of biologically plausible neural network models and more generic nonlinear and parametric oscillators. We show that, contrary to what is usually assumed in the literature, the high harmonic content observed in nonsinusoidal oscillatory dynamics is neither a sufficient nor necessary condition to interpret the associated CFC patterns as epiphenomenal. There is mounting evidence suggesting that the combination of multimodal recordings, specialized signal processing techniques, and theoretical modeling is becoming a required step to completely understand CFC patterns observed in oscillatory rich dynamics of physical and biophysical systems.

23 More

Received 2 August 2020
Revised 29 October 2020
Accepted 3 November 2020

DOI:https://doi.org/10.1103/PhysRevE.102.062401

Physics Subject Headings (PhySH)

Biological neural networks Coupled oscillators Dynamics of networks Neural encoding Neural information Pattern formation Population dynamics Synchronization Synchronization transition

Nonlinear DynamicsPhysics of Living Systems

Authors & Affiliations

Damián Dellavale ^*,†, Osvaldo Matías Velarde^*, Germán Mato , and Eugenio Urdapilleta^‡

Centro Atómico Bariloche and Instituto Balseiro, Comisión Nacional de Energía Atómica (CNEA), Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), Universidad Nacional de Cuyo (UNCUYO), Av. E. Bustillo 9500, R8402AGP San Carlos de Bariloche, Río Negro, Argentina

^*These authors contributed equally to this work.
^†Corresponding author: dellavale@cab.cnea.gov.ar
^‡Corresponding author: urdapile@gmail.com

Article Text (Subscription Required)

Click to Expand

References (Subscription Required)

Click to Expand

Issue

Vol. 102, Iss. 6 — December 2020

Reuse & Permissions

Access Options

Author publication services for translation and copyediting assistance advertisement

Images

Figure 1
Biologically plausible network for unidirectional PAC with external drive representing a slow sensory input entraining fast oscillations underpinning local neural processing in a cortical oscillator (sensory entrainment).
Reuse & Permissions
Figure 2
Synthetic amplitude-modulated signals and derived time series involved in the algorithms for quantification of PAC and harmonicity. Amplitude-modulated signals were computed as described in Appendix pp1-s1 using a sinusoidal modulating at $f_{LF} = 9$ Hz and modulated oscillations at $f_{HF} = 7 \times f_{LF} = 63$ Hz and $f_{HF} = 7.1 \times f_{LF} = 63.9$ Hz for the harmonic PAC and nonharmonic PAC, respectively. The LF signals $[x_{LF} (t)]$ were obtained by filtering the raw signal $x (t)$ using a band-pass filter centered at 9 Hz and a null-to-null bandwidth of 9 Hz. The HF signals $[x_{HF} (t)]$ were obtained by filtering the raw signal $x (t)$ using a band-pass filter centered at 63 Hz and a null-to-null bandwidth of 81 Hz (see Appendix pp1-s5). (a1), (a2) Synthetic amplitude-modulated signals. (b1), (b2) Time series used to compute the TLI metric to quantify spectral harmonicity. Note that the synthetic harmonic and nonharmonic PAC patterns are characterized by $TLI \approx 1$ and $TLI \approx 0$ values, respectively. (c1), (c2) Time series used to compute the PLV and KLMI metrics to quantify PAC. The histograms show the MI from which the KLMI can be computed. The histograms show the distribution of amplitude of HF as a function of the phase of LF.
Reuse & Permissions
Figure 3
Coupling between low and high frequency signals. (a), (c) Upper panels correspond to phase-amplitude coupling via harmonic content: a periodical nonsinusoidal signal, composed by a fundamental sinusoid at 8 Hz plus the first three harmonics with decreasing power, is stereotypically repeated over time [black line in (a)]. By chopping the signal in consecutive segments (red boxes, one prototypical in dark red, others in light red), whose length is the period of the fundamental rhythm, the very same signal is obtained [see (c)]. In (c), the true high frequency signal, i.e., the deterministic (noiseless) signal minus the fundamental oscillatory component (thus avoiding filtering artifacts), is shown as well as the arrow corresponding to each maxima in the chopped high frequency signals (red arrow). These maxima always lie in the very same position compared to the low frequency maxima (black arrow). (d) For comparison purposes, the deterministic fundamental low frequency signal can be observed in (d). (b), (e) Lower panels correspond to a phase-amplitude modulated signal: the phase of a fundamental sinusoid at 8 Hz modulates the intensity of a high frequency signal at 65 Hz. Here, high frequency signals corresponding to chopped segments (blue boxes) does not result in a single trace [see (e)]. Since the modulation is developed only through amplitude, low and high frequency signals are not tightly coupled regarding phase relationships, and each maximum of the high frequency chopped signal (blue arrows) has a distribution (over phases, or relative time) with respect to the low frequency maxima (black arrow). A small level of additive white Gaussian noise [see noisy signals in gray in panels (a) and (b)] does not change conclusions and a concomitant dispersion in the location of high frequency maxima may be observed.
Reuse & Permissions
Figure 4
Hilbert-Filter method for instantaneous frequency estimation. (a) Dynamics of the parametric oscillator (solid black line) generated by simultaneously applying an off-resonance external driving $F_{e}$ and a parametric driving $W_{p}$ tuned at the same frequency $f_{e} = f_{p} = f_{0} / 12 \approx 8.33$ Hz and $θ_{e} = 0$ [see Eqs. (A16) and (A17) in Appendix pp1-s3]. The LF (solid green line) and HF (solid red line) signals where obtained band-pass filtering the raw signal (solid black line) using the BPF whose power responses (i.e., square magnitude) are shown in graph (d) as dotted green and red lines, respectively. The configuration for the parametric oscillator used in these plots is identical to that used in graphs (d) and (e) of Fig. 17. (b) Instantaneous frequency time series computed for the raw signal (solid black line) shown in graph (a). Solid and dotted black lines correspond to the instantaneous frequency computed using the Hilbert-Filter method [Eq. (19)] and the numerical derivative [Eq. (11)], respectively. (c) Instantaneous phase of the frequency time series shown in graph (b), computed via Hilbert transformation. (d) Power spectrum (solid blue line) of the dynamics of the parametric oscillator [solid black line in graph (a)]. The power responses (i.e., square magnitude) of the BPF used to compute the LF and HF signals are shown as dotted green and red lines, respectively. (e) The solid black line represents the power response (i.e., square magnitude) of the band-pass filter $h_{2} (t)$ used to compute the instantaneous frequency time series by means of the Hilbert-Filter method [Eq. (19)]. The band-pass filter $h_{2} (t)$ was implemented as described in Appendix pp1-s5 using a Tukey window in the frequency domain. The dotted black line represents the power response of the low-pass filter $h_{LPF} (t)$ used to compute the instantaneous frequency time series by means of the numerical derivative [Eq. (11)]. (f) Magnitude responses of the band-pass filter $h_{2} (t)$ and the ideal derivator $|ω| f_{s}$ . Regarding the oversampling condition discussed in Appendix pp1-s6, in this case the oversampling ratio is $OSR = f_{s} / f_{LF} = 2000 / 8.33 \approx 240$ .
Reuse & Permissions
Figure 5
Performance of the TLI and ${PLV}_{PPC}$ in quantifying the harmonicity of a synthetic dynamics constituted by the linear superposition of two sinusoidal oscillations at $f_{0} = f_{LF} = 9$ Hz and $f_{HF} = 7 \times f_{LF} = 63$ Hz. In all the cases shown in this figure, we used a sampling rate of $f_{s} = 2000$ Hz and the frequency and amplitude of the LF and HF oscillations were kept unchanged. To obtain all the band-pass filtered signals shown in this figure, we use the BPF as described in Appendix pp1-s5. The bandwidth of the BPF for the LF component (LF BPF) was kept fixed at $B w_{LF} = 9$ Hz. The ${PLV}_{PPC}$ was computed using Eq. (4) with the configuration given by Eq. (A19) and $M = 1$ , $N = 7$ . (a), (d) TLI and ${PLV}_{PPC}$ metrics as a function of the epoch length and taking the level of additive white Gaussian noise (AWGN) as a parameter. The noise level is expressed as the percent of the amplitude of the LF component at $f_{LF} = 9$ Hz scaling the standard deviation $σ$ of the additive white Gaussian noise $N (0, σ)$ . To compute graphs (a) and (d), the bandwidth of the HF BPF was kept unchanged in $B w_{HF} = 99$ Hz. Our implementation of the TLI algorithm (Sec. 2d) requires at least three cycles of the low frequency oscillation ( $f_{LF} = 9$ Hz), which determines the minimum epoch length shown in graphs (a) and (d) ( $3 / f_{LF} \approx 0.3$ s). The maximum epoch length used to compute graphs (a) and (d) was $100 / f_{LF} \approx 11.1$ s. (b), (e) TLI and ${PLV}_{PPC}$ metrics as a function of the HF bandwidth ( $B w_{HF}$ ) corresponding to the BPF used to obtain the HF signal $[x_{HF} (t)]$ , and taking the level AWGN as a parameter. The minimum and maximum $B w_{HF}$ values used to compute the graphs (b) and (e) were 9 and 99 Hz, respectively. To compute the graphs (b) and (e), the epoch length was kept unchanged in $45 / f_{LF} \approx 5$ s. In (a), (b), (d), and (e), the solid lines represent the mean values and the shaded error bars correspond to the standard deviation of 100 realizations at each point. (c) Synthetic dynamics (solid black line) together with the HF and LF signals shown as solid red and green lines, respectively. The synthetic dynamics includes additive white Gaussian noise $N (0, σ)$ with the standard deviation $σ$ corresponding to the $40 %$ of the amplitude of the LF component at $f_{LF} = 9$ Hz. The LF and HF signals where obtained by filtering the raw signal with the band-pass filters whose power responses are shown as dotted green ( $B w_{LF} = 9$ Hz) and red ( $B w_{HF} = 99$ Hz) lines in graph (f), respectively. (f) Power spectrum (solid blue line) of the synthetic dynamics [solid black line in graph (c)] computed using an epoch length of $100 / f_{LF} \approx 11.1$ s. The power responses (i.e., square magnitude) of the BPF used to compute the LF and HF signals are shown as dotted green and red lines, respectively.
Reuse & Permissions
Figure 6
Performance of the ${PLV}_{PPC}$ in quantifying the harmonicity of a synthetic dynamics constituted by the linear superposition of two sinusoidal oscillations. In all the cases shown in this figure, we used a sampling rate of $f_{s} = 2000$ Hz and the amplitude of the LF and HF oscillations were kept unchanged. To obtain all the band-pass filtered signals shown in this figure we use the BPF as described in Appendix pp1-s5. The bandwidth of the BPF for the LF (LF BPF) and HF (HF BPF) components were kept fixed at $B w_{LF} = B w_{HF} = 3$ Hz. The ${PLV}_{PPC}$ was computed using Eq. (4) with the configuration given by Eq. (A19). (A) ${PLV}_{PPC}$ intensity as a function of the frequency ratio $f_{HF} / f_{LF}$ and taking the level of additive white Gaussian noise (AWGN) as a parameter. The LF component was kept fixed at $f_{LF} = 3$ Hz and the frequency of the HF oscillation was varied in the range $2 \times f_{LF} \leq f_{HF} \geq 180 \times f_{LF}$ . The noise level is expressed as the percent of the amplitude of the LF component at $f_{LF} = 3$ Hz scaling the standard deviation $σ$ of the additive white Gaussian noise $N (0, σ)$ . The ${PLV}_{PPC}$ was computed using an epoch length of $45 / f_{LF} \approx 5$ s. (b), (e), (h) Evolution of the ${PLV}_{PPC}$ intensity in-between the harmonic frequency ratios $f_{HF} / f_{LF}$ for three AWGN levels. In (a), (b), (e), and (h), the solid lines represent the mean values and the shaded error bars correspond to the standard deviation of 100 realizations at each point. (c), (f), (i) Power spectrum (solid blue line) of the synthetic dynamics [solid black line in graphs (d), (g) and (j)] corresponding to three frequency ratio values ( $f_{HF} / f_{LF} = 3, 89, 179$ ). The power spectra were computed using an epoch length of $45 / f_{LF} \approx 5$ s. The power responses (i.e., square magnitude) of the BPF used to compute the LF and HF signals are shown as dotted green and red lines, respectively. (d), (g), (j) Synthetic dynamics (solid black line) together with the HF and LF signals shown as solid red and green lines, respectively, corresponding to three frequency ratio values ( $f_{HF} / f_{LF} = 3, 89, 179$ ). The synthetic dynamics includes additive white Gaussian noise $N (0, σ)$ with the standard deviation $σ$ corresponding to the $25 %$ of the amplitude of the LF component at $f_{LF} = 9$ Hz. The LF and HF signals where obtained by filtering the raw signal with the band-pass filters whose power responses are shown as dotted green ( $B w_{LF} = 3$ Hz) and red ( $B w_{HF} = 3$ Hz) lines in (c), (f) and (i), respectively.
Reuse & Permissions
Figure 7
Performance of the TLI metric in quantifying the harmonicity of a synthetic dynamics constituted by the linear superposition of two sinusoidal oscillations. For the sake of comparison purposes, this figure was computed using the same hyperparameters as those used to compute and process the synthetic dynamics shown in Fig. 6.
Reuse & Permissions
Figure 8
Performance of harmonicity (TLI, ${PLV}_{PPC}$ ) and PAC ( ${KLMI}_{PAC}$ ) metrics to characterize an amplitude-modulated synthetic signal. In all the cases shown in this figure, we used a sampling rate of $f_{s} = 2000$ Hz and the frequency and amplitude of the modulating (LF) and the amplitude-modulated (HF) oscillations were kept unchanged. To obtain all the band-pass filtered signals shown in this figure we use the BPF as described in Appendix pp1-s5. The bandwidth of the BPF for the LF component (LF BPF) was kept fixed at $B w_{LF} = 9$ Hz. The ${PLV}_{PPC}$ was computed using Eq. (4) with the configuration given by Eq. (A19) and $M = 1$ , $N = 7$ . The amplitude-modulated signal was synthesized as described in Appendix pp1-s1 using the following hyperpameter values: $c = 1$ (i.e., DSB-C), maximum modulation depth $m = 0$ , $η_{m} = 0$ , we used a sinusoidal modulating $a (t)$ at $f_{0} = f_{LF} = 9$ Hz as given by Eq. (A6), $A_{m} = 1$ , $ϕ_{m} = 0$ , $z_{DBS}$ was set with $f_{HF} = 7 \times f_{LF} = 63$ Hz, $ϕ_{c} = 0$ , $z_{HF} = 0$ , for $z_{h}$ we use $A_{1} = 4$ , $A_{k} = 0 \forall k > 1$ and $ϕ_{k} = 0 \forall k$ . In Eq. (A1), we configured a constant amplitude envelope $E (t) = 1$ . The extrinsic noise level shown in the graphs corresponds to $η (t)$ in Eq. (A1), and is expressed as the percent of the modulating signal $a (t)$ maximum amplitude ( $A_{m}$ ) scaling the standard deviation $σ$ of the additive white Gaussian noise (AWGN) $η \approx N (0, σ)$ . (a), (c), (e) Harmonicity (TLI, ${PLV}_{PPC}$ ) and PAC ( ${KLMI}_{PAC}$ ) metrics as a function of the epoch length and taking the level of AWGN as a parameter. To compute graphs (a), (c), and (e), the bandwidth of the HF BPF was kept unchanged in $B w_{HF} = 99$ Hz. Our implementation of the TLI algorithm (Sec. 2d) requires at least three cycles of the low frequency oscillation ( $f_{LF} = 9$ Hz), which determines the minimum epoch length shown in graphs (a) and (d) ( $3 / f_{LF} \approx 0.3$ s). The maximum epoch length used to compute graphs (a) and (d) was $100 / f_{LF} \approx 11.1$ s. (b), (d), (f) Harmonicity (TLI, ${PLV}_{PPC}$ ) and PAC ( ${KLMI}_{PAC}$ ) metrics as a function of the HF bandwidth ( $B w_{HF}$ ) corresponding to the BPF used to obtain the HF signal $[x_{HF} (t)]$ , and taking the level AWGN as a parameter. The minimum and maximum $B w_{HF}$ values used to compute the graphs (b) and (e) were 18 and 99 Hz, respectively. To compute the graphs (b), (d), and (f), the epoch length was kept unchanged in $45 / f_{LF} \approx 5$ s. In (a)–(f), the solid lines represent the mean values and the shaded error bars correspond to the standard deviation of 100 realizations at each point. (g) Synthetic dynamics (solid black line) together with the HF and LF signals shown as solid red and green lines, respectively. The synthetic dynamics includes additive white Gaussian noise $N (0, σ)$ with the standard deviation $σ$ corresponding to the $40 %$ of the modulating signal $a (t)$ maximum amplitude ( $A_{m}$ ). The LF and HF signals where obtained by filtering the raw signal with the band-pass filters whose power responses are shown as dotted green ( $B w_{LF} = 9$ Hz) and red ( $B w_{HF} = 99$ Hz) lines in graph (h), respectively. (h) Power spectrum (solid blue line) of the synthetic dynamics [solid black line in graph (g)] computed using an epoch length of $100 / f_{LF} \approx 11.1$ s. The power responses (i.e., square magnitude) of the BPF used to compute the LF and HF signals are shown as dotted green and red lines, respectively.
Reuse & Permissions
Figure 9
Temporal evolution of the PAC ( ${PLV}_{PAC}$ ), harmonicity (TLI), and phase clustering ( ${PC}_{LF}$ ) metrics during a synthetic dynamics presenting a transient harmonic PAC pattern. To obtain all the band-pass filtered signals shown in this figure we use the BPF as described in Appendix pp1-s5. (a) Synthetic dynamics (solid black line) together with the HF and LF signals shown as solid red and green lines, respectively. The dynamics (solid black line) was synthesized using Eqs. (A1) and (A4) with the following hyperpameter values: sampling rate $f_{s} = 2000$ Hz, $c = 1$ (i.e., DSB-C), half modulation depth $m = 0.5$ , $η_{m} = 0$ , we used a Gaussian modulating $a (t)$ with the fundamental frequency at $f_{0} = f_{LF} = 3$ Hz as given by Eqs. (A7) and (A8) with $σ \approx 55$ and $A_{m} = 1$ , $z_{DBS}$ was set with $f_{HF} = 89 \times f_{LF} = 267$ Hz, $ϕ_{c} = 0$ , $z_{HF} = 0$ , for $z_{h}$ we use $A_{1} = 4$ , $A_{k} = 1 \forall 2 \leq k \leq 4$ , $A_{k} = 0 \forall k \geq 5$ , and $ϕ_{k} = 0 \forall k$ . The transient harmonic PAC pattern was implemented through the time series envelope $E (t)$ as defined in Eqs. (A2) and (A3), with $α = 0.5$ and $β$ equal to one third of the time series length. Extrinsic noise $η (t)$ was added as shown in Eq. (A1). In this case the noise level corresponds to the $10 %$ of the maximum amplitude of the deterministic part of signal $x (t)$ [i.e., first term of the right-hand side of the Eq. (A1)], scaling the standard deviation $σ$ of the additive white Gaussian noise (AWGN) $η \approx N (0, σ)$ . The LF (solid green line) and HF (solid red line) signals were obtained by filtering the raw signal (solid black line) with the band-pass filters whose power responses are shown as dotted green ( $B w_{LF} = 1$ Hz) and red ( $B w_{HF} = 30$ Hz) lines in graph (e), respectively. (b) Time series showing the temporal evolution of the ${PLV}_{PAC}$ , TLI, and ${PC}_{LF}$ metrics. These time series were computed as described in Sec. 2h using the Algorithm 2 summarized in Table 2, with a sliding window of 20 s in length, i.e., 60 cycles of the slowest oscillatory component at $f_{0} = f_{LF} = 3$ Hz. (c) TLI harmonicity map computed as described in Sec. 2g using a 20-s epoch extracted from the center ( $time \approx 100$ s) of the synthetic dynamics shown in (a). In computing the map, all the TLI values below the significance threshold were set to zero (see Sec. 2g). The pseudocolor scale represents the TLI values ranging from 0 (blue) to 1 (red). (d) Zoom showing two cycles of the synthetic dynamics (solid black line) together with the HF and LF signals shown as solid red and green lines, respectively. The two cycle epoch corresponds to the center ( $time \approx 100$ s) of the synthetic dynamics shown in panel (a). (e) Power spectrum (solid blue line) computed from the synthetic dynamics [solid black line in graph (a)]. The power responses (i.e., square magnitude) of the BPF used to compute the LF and HF signals are shown as dotted green and red lines, respectively. (f) Comodulogram computed as described in Sec. 2g computed from the same epoch used to obtain the harmonicity map (c). In computing the comodulogram, all the $| {PLV}_{PAC} |$ values below the significance threshold were set to zero (see Sec. 2g). The pseudocolor scale represents the $| {PLV}_{PAC} |$ values ranging from 0 (blue) to 1 (red). The harmonicity map (c) and comodulogram (f) were computed using the same BPF (see Appendix pp1-s5).
Reuse & Permissions
Figure 10
Temporal evolution of the PAC ( ${PLV}_{PAC}$ ), harmonicity (TLI), and phase clustering ( ${PC}_{LF}$ ) metrics during a synthetic dynamics presenting a transient harmonic PAC pattern. In this plot we use the same synthetic dynamics and the same set of hyperparameter values to compute the metrics than those described in the caption of Fig. 9, except for the bandwidth of the BPF used to compute the LF component ( $B w_{LF}$ ). In this case, the ${PLV}_{PAC}$ , TLI, and ${PC}_{LF}$ metrics were computed using $B w_{LF} = 13.5$ Hz centered around 7.5 Hz [see the dotted green line in (e)]. This wide BPF produces a nonsinusoidal LF component [see solid green line in (d)], characterized by a nonuniform distribution of phase values producing the increase of the phase clustering ( ${PC}_{LF}$ ) during the dynamics [see solid red line in (b)]. Note the bias in the PAC ( ${PLV}_{PAC}$ ) and harmonicity (TLI) metrics due to the presence of phase clustering ( ${PC}_{LF}$ ). The description of the panels is the same than that given in Fig. 9.
Reuse & Permissions
Figure 11
Comparison of the Algorithms 1 and 2 is summarized in Table 2 aimed to compute time series of CFC metrics. The time series showing the temporal evolution of the ${PLV}_{PAC}$ , TLI and ${PC}_{LF}$ metrics were computed as described in Sec. 2h using a sliding window of 3.33 s in length, i.e., 10 cycles of the slowest oscillatory component at $f_{0} = f_{LF} = 3$ Hz. (a), (c) The time series shown in (a) and (c) were computed from a synthetic dynamics presenting a transient harmonic PAC pattern, as described in Sec. 2h using Algorithms 1 and 2, respectively. The synthetic dynamics used in (a) and (c) was computed using the same set of hyperparameter values as those described in the caption of Fig. 9, with the exception of the extrinsic noise $η (t)$ which in this case was set to $30 %$ of the maximum amplitude of the deterministic part of signal $x (t)$ [i.e., first term of the right-hand side of Eq. (A1)]. (b), (d) The time series shown in (b) and (d) were computed from a synthetic dynamics without PAC, as described in Sec. 2h using Algorithms 1 and 2, respectively. The synthetic dynamics used in (b) and (d) was computed using the same set of hyperparameter values as those described in the caption of Fig. 20, with the exception of the extrinsic noise $η (t)$ which in this case was set to $30 %$ of the maximum amplitude of the deterministic part of signal $x (t)$ [i.e., first term of the right-hand side of Eq. (A1)]. In (a)–(d), the solid lines represent the mean values and the shaded error bars correspond to the standard deviation of 100 realizations at each point.
Reuse & Permissions
Figure 12
Harmonicity-PAC plot computed for a variety of synthetic oscillatory dynamics and taking the amplitude modulation depth as a parameter. The harmonicity-PAC plot was computed as it is described in Sec. 2f. The pseudocolor scale represents the $1 - m$ values ranging from the minimum 0 (blue) to the maximum 1 (red) modulation depth, with $m$ defining the amplitude modulation depth as stated in Eq. (A4). The synthetic dynamics was computed as it is described in Appendix pp1-s1 using sinusoidal modulating signals [Eq. (A6)].
Reuse & Permissions
Figure 13
The TLI is more robust than the ${PLV}_{PPC}$ against changes in the bandwidth of the BPF used to compute the high frequency component. The Van der Pol oscillatory dynamics (solid black line) shown in (a)–(d) were simulated as described in Appendix pp1-s2 using the following hyperpameter values: sampling rate $f_{s} = 2000$ Hz, $μ = 300$ , $ω_{0} = 2 π f_{0}$ , $f_{0} = 10$ Hz, $W_{p} = 0$ , $F_{e} = 0$ , initial conditions $x (0) = 2$ , $\dot{x} (0) = 1$ . With this configuration, we obtain a nonsinusoidal oscillatory dynamics constituted by a fundamental component at $f_{d} = 5.56$ Hz and odd harmonic components at $N \times f_{d}$ with $N = 3, 5, 7, 9, 11, 13, \dots$ . The dynamics was computed without intrinsic noise [ $g_{1} = g_{2} = 0$ in Eq. (A14)] in order to obtain a nonsinusoidal oscillatory dynamics with a constant fundamental period. Extrinsic noise $η (t)$ was added as shown in Eq. (A14). In this case the noise level corresponds to the $20 %$ of the maximum amplitude of the deterministic part of signal [i.e., $x_{1} (t)$ in Eq. (A14)], scaling the standard deviation $σ$ of the additive white Gaussian noise (AWGN) $η \approx N (0, σ)$ . We computed the Van der Pol dynamics for 6-s time interval and then the first 0.1 s (200 samples) of the time series were discarded to remove the transient period of the numerical simulation. For this set of hyperparameter values, (a)–(d) show different realizations of the Van der Pol dynamics. To obtain all the band-pass filtered signals shown in this figure we use the BPF as described in Appendix pp1-s5. The bandwidth of the BPF for the LF component (LF BPF) was kept fixed at $B w_{LF} = f_{d} = 5.56$ Hz [see the dotted green lines superimposed to the power spectra shown in (a)–(d)]. The band-pass filters used to compute the harmonicity metrics were centered at $f_{LF} = 1 \times f_{d}$ Hz (LF BPF) and $f_{HF} = 9 \times f_{d}$ Hz (HF BPF) and, consequently, the ${PLV}_{PPC}$ metric was computed using Eq. (4) with $M = 1$ , $N = 9$ . In (a)–(d) we changed the bandwidth of the BPF for the HF component (HF BPF) whose power response (i.e., square magnitude) is shown as dotted red line superimposed to the power spectra. The resulting band-pass filtered HF signals are shown as solid red lines in the upper graph of (a)–(d). The power spectra and the harmonicity metrics (TLI, ${PLV}_{PPC}$ ) were computed using an epoch length of $\approx 33 / f_{d} \approx 6$ s. (e) Shows the magnitude of the harmonicity (TLI, ${PLV}_{PPC}$ ) and phase clustering metrics ( ${PC}_{LF}$ , ${PC}_{HF}$ ) as a function of the bandwidth of the BPF for the HF component (HF BPF). In (e), the solid lines represent the mean values and the shaded error bars correspond to the standard deviation of 100 realizations at each HF bandwidth value.
Reuse & Permissions
Figure 14
The phase clustering associated to the low frequency component ( ${PC}_{LF}$ ) produces a bias in both TLI and ${PLV}_{PPC}$ metrics. In this figure we use the same synthetic dynamics and the same set of hyperparameter values to compute the metrics than those described in the caption of Fig. 13, except for the configuration of the BPFs used to compute the LF and HF components. In this case, the bandwidth of the BPF for the HF component (HF BPF) was kept fixed at $B w_{HF} = f_{d} = 5.56$ Hz [see the dotted red lines superimposed to the power spectra shown in (a)–(c)]. Besides, the band-pass filters used to compute the harmonicity metrics were centered at $f_{LF} = 5 \times f_{d}$ Hz (LF BPF) and $f_{HF} = 15 \times f_{d}$ Hz (HF BPF) and, consequently, the ${PLV}_{PPC}$ metric was computed using Eq. (4) with $M = 5$ , $N = 15$ . In (a)–(c) we changed the bandwidth of the BPF for the LF component (LF BPF) whose power response (i.e., square magnitude) is shown as dotted green line superimposed to the power spectra. The resulting band-pass filtered LF signals are shown as solid green lines in the upper graph of (a)–(c). The power spectra and the harmonicity metrics (TLI, ${PLV}_{PPC}$ ) were computed using an epoch length of $\approx 33 / f_{d} \approx 6$ s. (d) Shows the magnitude of the harmonicity (TLI, ${PLV}_{PPC}$ ) and phase clustering metrics ( ${PC}_{LF}$ , ${PC}_{HF}$ ) as a function of the bandwidth of the BPF for the LF component (LF BPF). In (d), the solid lines represent the mean values and the shaded error bars correspond to the standard deviation of 100 realizations at each LF bandwidth value.
Reuse & Permissions
Figure 15
Harmonicity-PAC plot computed for the simulated dynamics of the Van der Pol oscillator with intrinsic noise of type additive white Gaussian noise (AWGN). Note that a single nonsinusoidal oscillatory dynamics can produce both harmonic [(d) and (e), (i) and (j)] and nonharmonic [(k) and (l)] PAC patterns. The Van der Pol oscillatory dynamics (solid black line) shown in (a), (d), (f), (i), and (k) were simulated as described in Appendix pp1-s2 using the following hyperpameter values: sampling rate $f_{s} = 2000$ Hz, $ω_{0} = 2 π f_{0}$ , $f_{0} = 10$ Hz, $W_{p} = 0$ , $F_{e} = 0$ , initial conditions $x (0) = 2$ , $\dot{x} (0) = 1$ . To compute the harmonicity-PAC plots shown in (c) and (h), the parameter $μ$ controlling the oscillator nonlinearity was increased from the sinusoidal oscillatory regime [ $μ / ω_{0} \approx 0$ , see (a), (b), (f), and (g)] up to a high nonsinusoidal regime [ $μ / ω_{0} \approx 4.77$ , see (d), (e), (i), (j), (k), (l)]. In (c) and (h), the pseudocolor scale represents the $μ / ω_{0}$ values ranging from $\approx 0$ (blue) to $\approx 4.8$ (red). In (a)–(e), the dynamics of the Van der Pol oscillator was simulated using intrinsic noise of type AWGN applied only on the equation of ${\dot{x}}_{2}$ [ $g_{1} = 0$ and $g_{2} = 0.5$ in Eq. (A14)]. In (f)–(l), the dynamics was simulated by applying the intrinsic noise of type AWGN on the equations of both ${\dot{x}}_{1}$ and ${\dot{x}}_{2}$ [i.e., $g_{1} = g_{2} = 0.5$ in Eq. (A14)]. Therefore, in this case the intrinsic noise components (AWGN) in Eq. (A14) result $η_{1} \approx N (0, 0.5)$ and $η_{2} \approx N (0, 0.5)$ . Extrinsic noise $η (t)$ was added as shown in Eq. (A14). In this case, the noise level corresponds to the $10 %$ of the maximum amplitude of the dynamics $x_{1}$ in Eq. (A14), scaling the standard deviation $σ$ of the additive white Gaussian noise $η \approx N (0, σ)$ . We computed the Van der Pol dynamics for 20-s time interval and then the first 10 s of the time series were discarded to remove the transient period of the numerical simulation. The power spectra [solid blue line in (b), (e), (g), (j) and (l)] were computed using a 10-s epoch from the synthetic dynamics shown in the corresponding graphs [solid black line in (a), (d), (f), (i) and (k)]. In graphs (b), (e), (g), (j), and (l), the power responses (i.e., square magnitude) of the BPF used to compute the LF and HF signals are shown as dotted green and red lines, respectively. To obtain all the band-pass filtered signals shown in this figure we use the BPF as described in Appendix pp1-s5. In all the cases shown in this figure, the bandwidth of the BPF for the LF (LF BPF) and HF (HF BPF) components were set at $B w_{LF} = 12$ Hz centered at 7 Hz and $B w_{HF} = 82.65$ Hz centered at 56.5 Hz, respectively. In (a), (d), (f), (i), and (k), the resulting band-pass filtered LF and HF signals are shown as solid green and solid red lines, respectively. The harmonicity metric (TLI) was computed as it was described in Sec. 2d. The PAC metric ( ${PLV}_{PAC}$ ) was computed using Eq. (4) with the configuration given by Eq. (A20) and $M = N = 1$ .
Reuse & Permissions
Figure 16
Harmonicity-AAC and harmonicity-FFC plots computed for the simulated dynamics of the Van der Pol oscillator with intrinsic noise of type additive white Gaussian noise (AWGN). (a)–(e) To obtain the AAC pattern, the Van der Pol oscillatory dynamics (solid black line) shown in (a) and (d) were simulated as described in Appendix pp1-s2 using Eqs. (A12, A13, A14) with the following hyperpameter values: sampling rate $f_{s} = 2000$ Hz, $ω_{0} = 2 π f_{0}$ , $f_{0} = 10$ Hz, no parametric driving $W_{p} = 0$ , $f_{e} = 5.4$ Hz and $f_{m} = 1.33$ Hz, $A_{m} = 1$ , $c = 1$ (i.e., DSB-C), maximum modulation depth $m = 0$ , initial conditions $x (0) = 2$ , $\dot{x} (0) = 1$ . To compute the harmonicity-AAC plot shown in (c), the external driving amplitude [ $A_{e}$ in Eq. (A13)] was increased from $A_{e} = 0$ (no external driving) up to $A_{e} = 5 \times 10^{4}$ . In (c), the pseudocolor scale represents the $A_{e} / (5 \times 10^{4})$ values ranging from 0 (blue) to 1 (red). (f)–(j) To obtain the FFC pattern, the Van der Pol oscillatory dynamics (solid black line) shown in (f) and (i) were simulated as described in Appendix pp1-s2 using Eqs. (A12, A13, A14) with the following hyperpameter values: sampling rate $f_{s} = 2000$ Hz, $ω_{0} = 2 π f_{0}$ , $f_{0} = 10$ Hz, no external driving $F_{e} = 0$ , $f_{p} \approx 1$ Hz. To compute the harmonicity-FFC plot shown in (h), the intensity of the time variant parameter $W_{p}$ [ $A_{p}$ in Eq. (A12)] was increased from $A_{p} = 0$ (no parametric driving) up to $A_{p} \approx 10$ . In (h), the pseudocolor scale represents the $A_{p} / 34$ values ranging from 0 (blue) to 0.3 (red). For both (c) and (h), the dynamics of the Van der Pol oscillator was simulated using intrinsic noise of type AWGN applied only on the equation of ${\dot{x}}_{2}$ [ $g_{1} = 0$ and $g_{2} = 0.5$ in Eq. (A14)]. Therefore, the intrinsic noise components (AWGN) in Eq. (A14) result $η_{1} = 0$ and $η_{2} \approx N (0, 0.5)$ . Extrinsic noise $η (t)$ was added as shown in Eq. (A14). In this case, the noise level corresponds to the $10 %$ of the maximum amplitude of the dynamics $x_{1}$ in Eq. (A14), scaling the standard deviation $σ$ of the additive white Gaussian noise $η \approx N (0, σ)$ . The Van der Pol dynamics was computed for 20-s time interval and then the first 10 s of the time series were discarded to remove the transient period of the numerical simulation. The power spectra [solid blue line in graphs (b), (e), (g), and (j)] were computed unsing a 10-s epoch from the synthetic dynamics shown in the corresponding graphs [solid black line in graphs (a), (d), (f), and (i)]. In (b), (e), (g), and (j), the power responses (i.e., square magnitude) of the BPF used to compute the LF and HF signals are shown as dotted green and red lines, respectively. To obtain all the band-pass filtered signals shown in this figure we use the BPF as described in Appendix pp1-s5. In computing (c) (AAC), the bandwidth of the BPF for the LF (LF BPF) and HF (HF BPF) components were set at $B w_{LF} \approx 10.3$ Hz centered at 5.4 Hz and $B w_{HF} \approx 17.6$ Hz centered at 27 Hz, respectively. In computing (h) (FFC), the bandwidth of the BPF for the LF (LF BPF) and HF (HF BPF) components were set at $B w_{LF} \approx 10.3$ Hz centered at 5.4 Hz and $B w_{HF} \approx 10.8$ Hz centered at 16.2 Hz, respectively. In (a), (d), (f), and (i), the resulting band-pass filtered LF and HF signals are shown as solid green and solid red lines, respectively. The harmonicity metric (TLI) was computed as it was described in Sec. 2d. For (c), the AAC metric ( ${PLV}_{AAC}$ ) was computed using Eq. (4) with the configuration given by Eq. (A21) and $M = N = 1$ . For (h), the FFC metric ( ${PLV}_{FFC}$ ) was computed using Eq. (4) with the configuration given by Eq. (A24) and $M = N = 1$ .
Reuse & Permissions
Figure 17
Harmonicity-PFC plot computed for the simulated dynamics of the second order parametric oscillator with intrinsic noise of type additive white Gaussian noise (AWGN). Note that two oscillatory dynamics with independent frequencies can produce harmonic PFC patterns [(d), (e) and (f), (g)]. The parametric oscillator dynamics (solid black line) shown in (a), (d), (h), (f) and (k) were simulated as described in Appendix pp1-s3 by simultaneously applying an off-resonance external driving $F_{e}$ and a parametric driving $W_{p}$ tuned at the same frequency and using the following hyperpameter values: sampling rate $f_{s} = 2000$ Hz, $μ = 200$ , $ω_{0} = 2 π f_{0}$ , $f_{0} = 100$ Hz, $f_{p} = f_{e} = f_{0} / 12 \approx 8.3$ Hz, $θ_{e} = 0$ , $A_{e} = 1 \times 10^{5}$ . To compute the harmonicity-PFC plots shown in (c) and (j), the parameter $A_{p}$ controlling the parametric driving intensity was increased from the sinusoidal oscillatory regime [ $A_{p} \approx 0$ , see (a), (b), (h), (i)] up to a high nonsinusoidal regime [ $A_{p} \approx 0.9$ , see (d)–(g) ]. In (c) and (j), the pseudocolor scale represents the $A_{p}$ values ranging from 0 (blue) to 0.9 (red). In (a)–(e), the dynamics of the parametric oscillator was simulated using intrinsic noise of type AWGN applied only on the equation of ${\dot{x}}_{2}$ [ $g_{1} = 0$ and $g_{2} = 0.125$ in Eq. (A18)]. In (f)–(l), the dynamics was simulated by applying the intrinsic noise of type AWGN on the equations of both ${\dot{x}}_{1}$ and ${\dot{x}}_{2}$ [i.e., $g_{1} = g_{2} = 0.125$ in Eq. (A18)]. Therefore, in this case the intrinsic noise components (AWGN) in Eq. (A18) result $η_{1} \approx N (0, 0.125)$ and $η_{2} \approx N (0, 0.125)$ . Extrinsic noise $η (t)$ was added as shown in Eq. (A14). In this case, the noise level corresponds to the $5 %$ of the maximum amplitude of the dynamics $x_{1}$ in Eq. (A18), scaling the standard deviation $σ$ of the additive white Gaussian noise $η \approx N (0, σ)$ . We computed the dynamics of the parametric oscillator for 20-s time interval and then the first 10 s of the time series were discarded to remove the transient period of the numerical simulation. The power spectra [solid blue line in (b), (e), (g), (i), and (l)] were computed unsing a 10-s epoch from the synthetic dynamics shown in the corresponding graphs [solid black line in (a), (d), (f), (h), and (k)]. In (b), (e), (g), (i), and (l), the power responses (i.e., square magnitude) of the BPF used to compute the LF and HF signals are shown as dotted green and red lines, respectively. To obtain all the band-pass filtered signals shown in this figure, we use the BPF as described in Appendix pp1-s5. In all the cases shown in this figure, the bandwidth of the BPF for the LF (LF BPF) and HF (HF BPF) components were set at $B w_{LF} \approx 4.2$ Hz centered at $f_{0} / 12 \approx 8.3$ Hz and $B w_{HF} \approx 179$ Hz centered at $f_{0} = 100$ Hz, respectively. In (a), (d), (f), (h), and (k), the resulting band-pass filtered LF and HF signals are shown as solid green and solid red lines, respectively. The harmonicity metric (TLI) was computed as it was described in Sec. 2d. The PAC metric ( ${KLMI}_{PFC}$ ) was computed using Eqs. (6) and (7) with the configuration given by Eq. (A22). Note that the ${KLMI}_{PFC}$ was normalized with its maximum value in each plot.
Reuse & Permissions
Figure 18
Harmonicity-PAC plot computed for the simulated dynamics of the biologically plausible neural network model shown in Fig. 1 using the threshold linear activation function $S (I_{i})$ [Eq. (2)]. Note that two coupled oscillatory dynamics with independent fundamental frequencies can produce true PAC patterns with high harmonic content via rectification mechanisms [(d) and (e)]. The neural network dynamics (solid black line) shown in (a), (d), (f), (i), (k) and (n) were simulated as described in Sec. 2a using the configuration detailed in Table 1, resulting in an oscillatory dynamics in the gamma band (50 Hz). Besides, we use the following set of hyperpameter values: sampling rate $f_{s} = 20$ kHz, $H_{1} = 0$ , $H_{2} = A_{2} cos (2 π f_{2} t)$ with $f_{2} \approx 3.3$ Hz. To compute the harmonicity-PAC plots shown in (c), (h), and (m), the parameter $A_{2}$ controlling the amplitude of the oscillatory input $H_{2}$ was increased from $A_{2} = 0$ [see (a), (b), (f), (g), (k), (l)] up to $A_{2} = 0.3$ [see (d), (e), (i), (j), (n), (o)]. In (c), (h), and (m), the pseudocolor scale represents the $A_{2}$ values ranging from $\approx 0$ (blue) to $\approx 0.3$ (red). The neural network dynamics was simulated using intrinsic noise $η_{i}$ of type AWGN to the node inputs $I_{i}$ (see Sec. 2a). Panels (c), (h), and (m) were computed with a noise level of $5 %$ , $10 %$ , and $20 %$ of the external input maximum amplitude ( $A_{2} = 0.3$ ) scaling the standard deviation $σ_{i}$ of the additive white Gaussian noise $η_{i} \approx N (0, σ_{i})$ , respectively. We computed the neural network dynamics for 10-s time interval and then the first 5 s of the time series were discarded to remove the transient period of the numerical simulation. The power spectra [solid blue line in (b), (e), (g), (j), (l), and (o)] were computed using a 5-s epoch from the synthetic dynamics shown in the corresponding graphs [solid black line in (a), (d), (f), (i), (k), and (n)]. In (b), (e), (g), (j), (l), and (o), the power responses (i.e., square magnitude) of the BPF used to compute the LF and HF signals are shown as dotted green and red lines, respectively. To obtain all the band-pass filtered signals shown in this figure we use the BPF as described in Appendix pp1-s5. In all the cases shown in this figure, the bandwidth of the BPF for the LF (LF BPF) and HF (HF BPF) components were set at $B w_{LF} \approx 3.3$ Hz centered at $\approx 3.3$ Hz and $B w_{HF} \approx 43.2$ Hz centered at 50 Hz, respectively. In (a), (d), (f), (i), (k), and (n), the resulting band-pass filtered LF and HF signals are shown as solid green and solid red lines, respectively. The harmonicity metric (TLI) was computed as it was described in Sec. 2d. The PAC metric ( ${PLV}_{PAC}$ ) was computed using Eq. (4) with the configuration given by Eq. (A20) and $M = N = 1$ .
Reuse & Permissions
Figure 19
The TLI and ${PLV}_{PPC}$ metrics present a comparable bias when computed on nonharmonically related oscillations. In all the cases shown in this figure, we used a sampling rate of $f_{s} = 2000$ Hz and the bandwidth of the BPF for the LF component (LF BPF) was kept fixed at $B w_{LF} = f_{LF}$ . We use the BPF as described in Appendix pp1-s5. Besides, in all the cases shown in this figure the noise level is expressed as the percent of the amplitude of the LF component at $f_{LF}$ Hz scaling the standard deviation $σ$ of the additive white Gaussian noise $N (0, σ)$ . Panels (a), (b), (d), and (e) were computed using a synthetic dynamics similar to that used in Fig. 5, but in this case it is constituted by two nonharmonic oscillations at $f_{LF} = 9$ Hz and $f_{HF} = 7.2 f_{LF} = 64.8$ Hz. For panels (a), (b), (d), and (e), the ${PLV}_{PPC}$ was computed using Eq. (4) with the configuration given by Eq. (A19) and $M = 1$ , $N = 7$ . (a), (d) TLI and ${PLV}_{PPC}$ metrics as a function of the HF bandwidth ( $B w_{HF}$ ) corresponding to the filter HF BPF used to obtain the HF signal $[x_{HF} (t)]$ , and taking the level AWGN as a parameter. The minimum and maximum $B w_{HF}$ values used to compute the graphs (b) and (e) were 9 and 102.6 Hz, respectively. To compute these graphs, the epoch length was kept unchanged in 5 s. (b), (e) TLI and ${PLV}_{PPC}$ metrics as a function of the epoch length and taking the level of additive white Gaussian noise (AWGN) as a parameter. To compute graphs (b) and (e), the bandwidth of the filter HF BPF was kept unchanged at $B w_{HF} = 102.6$ Hz. Our implementation of the TLI algorithm (Sec. 2d) requires at least three cycles of the low frequency oscillation ( $f_{LF} = 9$ Hz), which determines the minimum epoch length shown in (a) and (d) ( $3 / f_{LF} \approx 0.3$ s). The maximum epoch length used to compute (a) and (d) was $100 / f_{LF} \approx 11.1$ s. (c), (f) The TLI and ${PLV}_{PPC}$ metrics as a function of the epoch length and taking the nonharmonic ratio $R = f_{HF} / f_{LF}$ as a parameter. Panels (c) and (f) were computed using a synthetic dynamics similar to that used in Fig. 5, but in this case it is constituted by two nonharmonic oscillations at $f_{LF} = 3$ Hz and $f_{HF} = R \times f_{LF}$ with $R = 3.2, 13.1, 89.1, 180.1$ . To compute graphs (c) and (f), the bandwidth of the filter HF BPF was kept unchanged at $B w_{HF} = B w_{LF} = f_{LF} = 3$ Hz. The noise level was set to $20 %$ of the amplitude of the LF component at $f_{LF} = 3$ Hz. The minimum and maximum epoch length shown in (c) and (f) are $3 / f_{LF} \approx 1$ s and $100 / f_{LF} \approx 33.3$ s, respectively. In all the panels, the solid lines represent the mean values and the shaded error bars correspond to the standard deviation of 100 realizations at each point.
Reuse & Permissions
Figure 20
Temporal evolution of the PAC ( ${PLV}_{PAC}$ ), harmonicity (TLI), and phase clustering ( ${PC}_{LF}$ ) metrics during a synthetic dynamics without PAC. To obtain all the band-pass filtered signals shown in this figure we use the BPF as described in Appendix pp1-s5. (a) Synthetic dynamics (solid black line) together with the HF and LF signals shown as solid red and green lines, respectively. The dynamics (solid black line) was synthesized using Eqs. (A1) and (A4) with the following hyperpameter values: sampling rate $f_{s} = 2000$ Hz, $c = 1$ (i.e., DSB-C), zero modulation depth $m = 0$ , $η_{m} = 0$ , we used a sinusoidal modulating $a (t)$ with the fundamental frequency at $f_{0} = f_{LF} = 3$ Hz as given by Eq. (A6) with $A_{m} = 1$ , $z_{DBS}$ was set with $f_{HF} = 89 \times f_{LF} = 267$ Hz, $ϕ_{c} = 0$ , $z_{HF} = 0$ , for $z_{h}$ we use $A_{1} = 4$ , $A_{k} = 1 \forall 2 \leq k \leq 4$ , $A_{k} = 0 \forall k \geq 5$ , and $ϕ_{k} = 0 \forall k$ . The transient pattern was implemented through the time series envelope $E (t)$ as defined in Eqs. (A2) and (A3), with $α = 0.5$ and $β$ equal to one third of the time series length. Extrinsic noise $η (t)$ was added as shown in Eq. (A1). In this case the noise level corresponds to the $10 %$ of the maximum amplitude of the deterministic part of signal $x (t)$ [i.e., first term of the right-hand side of Eq. (A1)], scaling the standard deviation $σ$ of the additive white Gaussian noise (AWGN) $η \approx N (0, σ)$ . The LF (solid green line) and HF (solid red line) signals where obtained by filtering the raw signal (solid black line) with the band-pass filters whose power responses are shown as dotted green ( $B w_{LF} = 1$ Hz) and red ( $B w_{HF} = 30$ Hz) lines in graph (e), respectively. (b) Time series showing the temporal evolution of the ${PLV}_{PAC}$ , TLI, and ${PC}_{LF}$ metrics. These time series were computed as described in Sec. 2h using the Algorithm 2 summarized in Table 2, with a sliding window of 20 s in length, i.e., 60 cycles of the slowest oscillatory component at $f_{0} = f_{LF} = 3$ Hz. (c) TLI harmonicity map computed as described in Sec. 2g using a 20-s epoch extracted from the center ( $time \approx 100$ s) of the synthetic dynamics shown in (a). In computing the map, all the TLI values below the significance threshold were set to zero (see Sec. 2g). The pseudocolor scale represents the TLI values ranging from 0 (blue) to 1 (red). (d) Zoom showing two cycles of the synthetic dynamics (solid black line) together with the HF and LF signals shown as solid red and green lines, respectively. The two cycle epoch corresponds to the center ( $time \approx 100$ s) of the synthetic dynamics shown in (a). (e) Power spectrum (solid blue line) computed from the synthetic dynamics [solid black line in (a)]. The power responses (i.e., square magnitude) of the BPF used to compute the LF and HF signals are shown as dotted green and red lines, respectively. (f) Comodulogram computed as described in Sec. 2g computed from the same epoch used to obtain the harmonicity map (c). In computing the comodulogram, all the $| {PLV}_{PAC} |$ values below the significance threshold were set to zero (see Sec. 2g). The pseudocolor scale represents the $| {PLV}_{PAC} |$ values ranging from 0 (blue) to 1 (red). The harmonicity map (c) and comodulogram (f) were computed using the same BPF (see Appendix pp1-s5).
Reuse & Permissions
Figure 21
Temporal evolution of the PAC ( ${PLV}_{PAC}$ ), harmonicity (TLI), and phase clustering ( ${PC}_{LF}$ ) metrics during a synthetic dynamics presenting a transient non harmonic PAC pattern. The synthetic dynamics was synthesized using the same parameter values than those used in Fig. 9, except for the frequency of the carrier in the $z_{DBS}$ signal which in this case was set to $f_{HF} = 88.9 \times f_{LF} = 88.9 \times 3 = 266.7$ Hz. The ${PLV}_{PAC}$ , TLI, and ${PC}_{LF}$ metrics were computed using the same set of hyperparameter values than those used in Fig. 9. The description of the panels is the same than that given in Fig. 9.
Reuse & Permissions
Figure 22
Temporal evolution of the PAC ( ${MVL}_{PAC}$ ), harmonicity (TLI), and phase clustering ( ${PC}_{LF}$ ) metrics during a synthetic oscillatory dynamics constituted by harmonically related rhythms with no PAC. The synthetic dynamics was synthesized using the same parameter values than those used in Fig. 20. The ${MVL}_{PAC}$ [see Eq. (5)], TLI, and ${PC}_{LF}$ metrics were computed using the same set of band-pass filters and hyperparameter values than those used in Fig. 20. The description of the panels is the same as that given in Fig. 20.
Reuse & Permissions
Figure 23
Temporal evolution of the PAC ( ${MVL}_{PAC}$ ), harmonicity (TLI), and phase clustering ( ${PC}_{LF}$ ) metrics during a synthetic oscillatory dynamics constituted by harmonically related rhythms with no PAC. In this plot we use the same synthetic dynamics and the same set of hyperparameter values to compute the metrics than those described in the caption of Fig. 20, except for the bandwidth of the BPF used to compute the LF component ( $B w_{LF}$ ). In this case, the ${MVL}_{PAC}$ [see Eq. (5)], TLI, and ${PC}_{LF}$ metrics were computed using $B w_{LF} = 13.5$ Hz centered around 7.5 Hz [see the dotted green line in (e)]. This wide BPF produces a nonsinusoidal LF component [see solid green line in (d)], characterized by a nonuniform distribution of phase values producing the increase of the phase clustering ( ${PC}_{LF}$ ) during the dynamics [see solid red line in (b)]. Note the bias in the PAC ( ${MVL}_{PAC}$ ) and harmonicity (TLI) metrics due to the presence of phase clustering ( ${PC}_{LF}$ ). The description of the panels is the same than that given in Fig. 20.
Reuse & Permissions
Figure 24
Phase portraits for the simulated dynamics of the Van der Pol oscillator. In this figure we use the same synthetic dynamics and the same set of hyperparameter values as those described in the caption of Fig. 15. In particular, the phase portraits were computed using the dynamics $x_{1}$ in Eq. (A14) which only takes into account the effect of the intrinsic noise, that is, without including the extrinsic (i.e., of observation) noise $η$ . (a) Phase portrait corresponding to the dynamics shown in Fig. 15 (no PAC). (b) Phase portrait corresponding to the dynamics shown in Fig. 15 (no PAC). (c) Phase portrait corresponding to the dynamics shown in Fig. 15 (harmonic PAC). (d) Phase portrait corresponding to the dynamics shown in Fig. 15 (nonharmonic PAC).
Reuse & Permissions
Figure 25
Harmonicity-PAC plot computed for the simulated dynamics of the Van der Pol oscillator with intrinsic noise of type nonadditive white Gaussian noise (NAWGN). In this figure we use the same synthetic dynamics and the same set of hyperparameter values to compute the metrics than those described in the caption of Fig. 15, except for the configuration of the intrinsic noise. In this case, we use nonadditive white Gaussian noise. For the numerical integration of the stochastic differential equation (A14) we use an explicit solver based on the Euler-Heun method [42] using the Stratonovich integral formulation. Importantly, we verified that the harmonicity-PAC plots shown in this figure do not change when computed using the Itô integral formulation. For (a) and (b), the dynamics of the Van der Pol oscillator was simulated using intrinsic noise of type NAWGN applied only on the equation of ${\dot{x}}_{2}$ [ $g_{1} = 0$ and $g_{2} = 0.5 x_{2}$ in Eq. (A14)]. For (c) and (d), the dynamics was simulated by applying the intrinsic noise of type nonadditive white Gaussian noise on the equations of both ${\dot{x}}_{1}$ and ${\dot{x}}_{2}$ [i.e., $g_{1} = 0.5 x_{1}$ and $g_{2} = 0.5 x_{2}$ in Eq. (A14)]. Therefore, in this case the intrinsic noise components (NAWGN) in Eq. (A14) result $η_{1} \approx N (0, 0.5 x_{1})$ and $η_{2} \approx N (0, 0.5 x_{2})$ . Extrinsic noise $η (t)$ was added as shown in Eq. (A14). In this case, the noise level corresponds to the $10 %$ of the maximum amplitude of the dynamics $x_{1}$ in Eq. (A14), scaling the standard deviation $σ$ of the additive white Gaussian noise $η \approx N (0, σ)$ . The harmonicity metric (TLI) was computed as it was described in Sec. 2d. For (a) and (c), the PAC metric ( ${PLV}_{PAC}$ ) was computed using Eq. (4) with the configuration given by Eq. (A20) and $M = N = 1$ . For (b) and (d), we compute the ${KLMI}_{PAC}$ using Eqs. (6) and (7) with the configuration given by Eq. (A20). Note that the ${KLMI}_{PAC}$ was normalized with its maximum value in each plot.
Reuse & Permissions
Figure 26
Phase portraits illustrating the simulated dynamics of the Van der Pol oscillator under an amplitude-modulated external driving. In this figure we use the same synthetic dynamics and the same set of hyperparameter values than those used to compute Fig. 16. In particular, the phase portraits were computed using the dynamics $x_{1}$ in Eq. (A14) which only takes into account the effect of the intrinsic noise, that is, without including the extrinsic (i.e., of observation) noise $η$ . (a) Phase portrait corresponding to the dynamics shown in Fig. 16 for $A_{e} / (5 \times 10^{4}) \approx 0.01$ . (b) Phase portrait corresponding to the dynamics shown in Fig. 16 for $A_{e} / (5 \times 10^{4}) \approx 0.1$ . (c) Phase portrait corresponding to the dynamics shown in Fig. 16 for $A_{e} / (5 \times 10^{4}) \approx 1$ .
Reuse & Permissions
Figure 27
Phase portraits for the simulated dynamics of the second order parametric oscillator. In this figure we use the same synthetic dynamics and the same set of hyperparameter values than those described in the caption of Fig. 17. In particular, the phase portraits were computed using the dynamics $x_{1}$ in Eq. (A18) which only takes into account the effect of the intrinsic noise, that is, without including the extrinsic (i.e., of observation) noise $η$ . (a) Phase portrait corresponding to the dynamics shown in Figs. 17 and 17 (no PFC). (b) Phase portrait corresponding to the dynamics shown in Figs. 17 and 17 (no PFC). (c) Phase portrait corresponding to an intermediate dynamics in-between the cases shown in Figs. 17 and 17 and Figs. 17 and 17. (d) Phase portrait corresponding to the dynamics shown in Figs. 17 and 17 (nonharmonic PFC). (e) Phase portrait corresponding to the dynamics shown in Figs. 17 and 17 (harmonic PFC). (f) Phase portrait corresponding to the dynamics shown in Figs. 17 and 17 (harmonic PFC).
Reuse & Permissions
Figure 28
Harmonicity-PFC plot computed for the simulated dynamics of the second order parametric oscillator with intrinsic noise of type additive white Gaussian noise (AWGN). Note that two oscillatory dynamics with independent frequencies can produce harmonic PFC patterns [(d), (e) and (f), (g)]. In this figure we use the same synthetic dynamics and the same set of hyperparameter values to compute the metrics than those described in the caption of Fig. 17, except for the phase of the external driving $θ_{e} = π / 2$ and the frequency of the parametric and external driving, which were configured as $f_{p} = f_{e} = f_{0} / 11.62 \approx 8.3$ Hz (i.e., $f_{p}$ and $f_{e}$ are nonharmonics of $f_{0}$ ). The harmonicity metric (TLI) was computed as it was described in Sec. 2d. The PFC metric ( ${KLMI}_{PFC}$ ) was computed using Eqs. (6) and (7) with the configuration given by Eq. (A22). Note that the ${KLMI}_{PFC}$ was normalized with its maximum value in each plot.
Reuse & Permissions
Figure 29
Harmonicity-PFC plot computed for the simulated dynamics of the second order parametric oscillator with intrinsic noise of type additive white Gaussian noise (AWGN). Note that two oscillatory dynamics with independent frequencies can produce harmonic PFC patterns [(d), (e) and (f), (g)]. In this figure we use the same synthetic dynamics and the same set of hyperparameter values to compute the filtering and harmonicity metric (TLI) than those described in the caption of Fig. 17, except for the phase of the external driving $θ_{e} = π / 2$ and the frequency of the parametric and external driving, which were configured as $f_{p} = f_{e} = f_{0} / 11.62 \approx 8.3$ Hz (i.e., $f_{p}$ and $f_{e}$ are nonharmonics of $f_{0}$ ). The harmonicity metric (TLI) was computed as it was described in Sec. 2d. In this case, the PFC metric ( ${PLV}_{PFC}$ ) was computed using Eq. (4) with the configuration given by Eq. (A22) and $M = 1$ , $N = 1$ .
Reuse & Permissions
Figure 30
Harmonicity-PAC plot computed for the simulated dynamics of the biologically plausible neural network model shown in Fig. 1 using the softplus activation function $S (I_{i})$ [Eq. (3)]. In this figure we use the same synthetic dynamics and the same set of hyperparameter values to compute the metrics than those described in the caption of Fig. 18, except for the activation function $S (I_{i})$ which in this case was computed using the Eq. (3) with $c = 20$ .
Reuse & Permissions

Physical Review E

covering statistical, nonlinear, biological, and soft matter physics