Phase-Amplitude Coupling in Neuronal Oscillator Networks

Phase-amplitude coupling (PAC) describes the phenomenon where the power of a high-frequency oscillation evolves with the phase of a low-frequency one. We propose a model that explains the emergence of PAC in two commonly-accepted architectures in the brain, namely, a high-frequency neural oscillation driven by an external low-frequency input and two interacting local oscillations with distinct, locally-generated frequencies. We further propose an interconnection structure for brain regions and demonstrate that low-frequency phase synchrony can integrate high-frequency activities regulated by local PAC and control the direction of information flow across distant regions.

Oscillatory activity at multiple frequency bands is widely observed in the brain, among other natural, biological and technological systems. It is thought that synchronization of brain oscillations within the same frequency band supports effective neural communication [1], and extensive studies have investigated this phenomenon [2][3][4][5]. Yet, there is an emerging body of evidence revealing that brain oscillations in distinct frequency bands are not independent, and they are instead found to interact through cross-frequency coupling [6,7].
As one of the most prevalent representations of crossfrequency coupling [8], phase-amplitude coupling (PAC) occurs when the power (or amplitude) of a high-frequency rhythm is locked to (and often modulated by) the phase of a low-frequency rhythm (see FIG. 1(a)). Brain signals recorded by various techniques, e.g., local field potential (LFP), electroencephalograph (EEG), and magnetoencephalograph (MEG), have revealed that PAC emerges in numerous brain regions, including auditory and prefrontal cortices [9,10], nucleus [11], and hippocampus [12], and plays a crucial role in motor functions [13] and cognitive processes such as working memory [14], attention [15], and learning [16].
Despite extensive empirical work on PAC, its underlying mechanisms remain largely unknown. Modeling PAC can help us uncover the conditions under which this phenomenon emerges, shedding light on the working principles of numerous cognitive functions. Most of the existing models focus on the microscopic scale [17][18][19][20], but they fail to capture the behavior of larger neural populations. Some attempts have been made to describe PAC at the population level with the aid of neural mass models [21][22][23][24]. Two architectures have been accepted to generate PAC [25]: 1) a neural population oscillating at high-frequency (fast) is modulated by an external lowfrequency (slow) input (see FIG. 1(b)); 2) two populations that respectively generate fast and slow oscillations interact with each other (see FIG. 1(c)). Yet, a model capable of assimilating the aforementioned instances with large-scale synchronization of multiple brain regions is still critically missing.
PAC usually does not emerge in the brain as an isolated phenomenon. In fact, regions exhibiting local PAC may also be phase-coupled in low frequency bands. The combination of cross-region phase synchronization and local cross-frequency PAC can thus effectively integrate local information in distant brain regions, which is processed and regulated by high-frequency oscillations [8]. For instance, correlated γ rhythms (>30 Hz) between the cortex and striatum are typically established by crossregional phase synchronization in the θ band (4-8 Hz) concurrently with θ-γ PAC [26]. Finally, while synchronization phenomena enable effective communication, the question of which regions send and which ones receive such information often remains unanswered. Determining the directionality of information flow is of paramount importance for revealing the spatio-temporal hierarchy of brain processes. For instance, it has been shown that during working memory maintenance the prefrontal cortex receives information from the sensory cortices [27] and sends information to the inferior temporal cortex [28].
In this Letter, we propose a neurobiologically plausible model based on the Stuart-Landau (SL) equation that allows for the emergence of cross-frequency coupling in both architectures depicted in FIG. 1(b)-(c). Further, inspired by previous studies showing that low-frequency oscillations are better suited for establishing long-distance interactions than high-frequency ones [26,29], we provide an interconnection structure (see FIG. 1(d)) in which cross-region interconnections exist only between slow populations. This structure creates a framework that integrates cross-frequency coupling with same-frequency phase synchronization, and constitutes a building block that can be easily scaled to form large networks. Under this structure, we demonstrate the important role that long-distance same-frequency phase synchrony, together with regional PAC, plays in the coordination of high-frequency local activity and in information routing. We start by introducing the SL equation: where z is a complex state (z = x + iy), with i = √ −1, and σ and ω are real. The SL oscillator is also known as a Hopf oscillator, because, for σ < 0, z = 0 is globally stable, and for σ > 0 the model undergoes a Hopf bifurcation with z converging to a circular trajectory (limit cycle) of radius √ σ and angular frequency ω. The SL oscillator can equivalently be described in polar coordinates by letting z = re iθ , where r and θ are referred to as the amplitude and phase angle of the oscillator, respectively. While phase-only models such as the Kuramoto model have been widely used to study synchrony in brain oscillations [4,5,[30][31][32], the SL model, like many other phase-amplitude models [33,34], can capture richer behaviors than phase-only ones. We use the SL equation to describe the collective dynamics of a single neural population. The SL equation can be derived from the well-established Wilson-Cowan model of interaction between excitatory and inhibitory neural populations [35], and from the reduction of a network of inhibitory integrate-and-fire neurons [36]. Despite its apparent simplicity, the SL equation is a very powerful approximation of neural dynamics around the Hopf bifurcation. In fact, local synaptic coupling can synchronize the firing activity of neurons [25], and then bring the population from a stationary to an oscillatory regime by shifting the parameter σ in Eq. (1) from negative to positive. The oscillation frequency mainly depends on the intrinsic decay time of inhibition in the population [37]. Emergence of PAC by exogenous inputs -We first show how PAC can arise in a neural population oscillating at a high frequency subject to a slow exogenous input (i.e.,

Input
r Oscillatory activity the architecture depicted in FIG. 1(b)). Consider the dynamics of the neural population with an input u where The input u can represent various sensory signals from sensory cortices [25,38]. We let u be a sinusoidal signal of strength k I , i.e., u = k I sin(ω I t). Consequently, Then, σ in Eq.
(1) becomes time-dependent, i.e., σ(t) = δ f + k I sin(ω I t), and the fast amplitude r f tends to converge to σ(t), yielding PAC. Notice that, depending on how the fast oscillation behaves throughout the slow cycle, PAC can be categorized into continuous and intermittent. In the former, the fast oscillation is constantly active, while in the latter, the fast oscillation only appears in a certain phase interval of the slow cycle [25]. The input magnitude determines which type of PAC will occur. For a weak input, where k I < δ f , σ(t) is always positive and the fast population remains in the oscillatory regime, which implies that continuous PAC arises (see FIG. 2(a)). Conversely, for a strong input, where k I > δ f , σ(t) can be negative in some interval of the slow cycle. In these intervals, the fast population is driven out of the oscillatory regime. Thus, r f tends towards 0, and the fast oscillation may disappear. We illustrate intermittent PAC for large k I in FIG. 2(b).
One can observe from FIG. 2(a)-(b) that the peaks and valleys of the fast amplitude oscillation appear almost simultaneously as those of the slow input oscillation, respectively. However, this situation does not always happen, and it highly depends on δ f in Eq. (2). As pointed out in Ref. [39], since σ in Eq. (1) determines the convergence rate to the limit cycle, we then reason that the rate of r f tracking σ(t) is also positively correlated with δ f . Therefore, a larger δ f means a smaller phase lag. FIG. 2(c) demonstrates that a noticeable phase lag appears for a smaller δ f in comparison with FIG. 2(a).
After revealing the emergence of PAC, we investigate how this phenomenon depends on the model parameters. More precisely, we utilize the Modulation Index (MI) [40] to measure PAC intensity (see definition in the Supplementary Material). As we have reasoned that the fast amplitude r f tends to track σ(t)= √ δ f · 1 + k I /δ f sin(ω I t), it is clear that a larger ratio k I /δ f implies larger fluctuations in r f than a weaker one. That is, the PAC intensity increases with k I /δ f (see FIG. 2(d)). Additionally, for a fixed a ratio k I /δ f , variations of δ f impact the fluctuations of r f : a larger δ f implies more intense PAC (see FIG. 2(d)).
The input frequency ω I also has a profound impact on the PAC intensity. A smaller ω I means a more slowly varying σ(t), which is easier for the fast amplitude r f to track. FIG. 2(e) shows that the PAC intensity decreases as ω I increases. By contrast, the fast frequency ω f has little influence on the PAC intensity (see FIG. 2(f)).
So far we have focused on the emergence of PAC in the architecture shown in FIG. 1(b). That is, a slowly oscillating input to the fast population modulates the fast amplitude. However, some in vitro studies have shown that besides high-frequency rhythms, low-frequency ones can also be generated locally in a single brain region [17,41]. These heterogeneous rhythms then interact as in the architecture depicted in FIG. 1(c). We now turn to the study of this architecture. PAC by local interactions: shaping of a brain region -As depicted in FIG. 1(c), a fast population may locally interact with a slow one. The two population dynamics arė where z s = x s + iy s , z f = x f + iy f , and δ s , δ f > 0. The oscillators' natural frequencies are determined by v s and v f , with v f > v s . The complex function f (·) describes how the two populations are interconnected. If the connection is unidirectional from the slow population to the fast one (a case studied in Ref. [42]), i.e., f (z f ) = 0, then Eq. (4) reduces to Eq. (2) with input u = f (z s ). In this case, our analysis above concludes that PAC emerges if f (z s ) is a function of the slow phase. We next turn our attention to the more interesting case wherein the interaction between the two populations is bidirectional. In this Letter we simply consider f (z) = kz. Then, Eq. (3) and Eq. (4) becomė  . 1 (c), where the slow frequency is 6.5 Hz (θ band), and the fast one is 30Hz (γ band). (c) PAC intensity (measured by the modulation index (MI)) increases with kδs/δ f . Given kδs/δ f , PAC becomes more intense for larger δ f .
In polar coordinates, z s = r s e iθs and z f = r f e iθ f . It can be observed that the terms σ and ω (see Eq. (1)) that describe the amplitude and phase in the slow and fast population read as σ s = δ s + kx f , ω s = v s + ky f , and σ f = δ f + kx s , ω f = v f + ky s , respectively. Interestingly, the amplitudes and phases of the two populations become dependent upon one another due to the interconnection. Notice that such dependence is asymmetric due to the frequency difference. We find that the slow oscillation remains relatively independent from the fast one. Since v f > v s , one can see that z f oscillates faster than the z s . Applying averaging techniques developed in Ref. [43] to the dynamics in Eq. (5), we obtain the approximate guiding systeṁ z s = (δ s +iv s −|ẑ s | 2 )ẑ s , which corresponds to the dynamics of the decoupled slow population. The solutions to Eq. (5) and the guiding system satisfy z s (t)−ẑ s (t) ≤ cv s /v f for some positive c. Loosely speaking, the slow trajectory z s (t) fluctuates in the neighborhood of the trajectory of its decoupled counterpartẑ s (t). Such fluctuations decrease as the two frequencies separate (see Fig. 3(a)). In contrast, Eq. (6) can be rewritten aṡ z f = (δ f + kr s cos θ s ) + i(v f + r s sin θ s ) − |z f | 2 z f . Similar to our earlier analysis, the fast amplitude r f tends to track √ δ f + kr s cos θ s , which implies PAC. Let the fast population in FIG. 1(c) oscillate in the γ band, and the slow population oscillate in the θ band (a choice compatible with empirical observations). FIG. 3(b) shows that PAC emerges from the interaction of these two populations. The PAC intensity here also depends on the slow amplitude r s , which is determined by δ s . Note that the connection strength k also affects PAC intensity. Therefore, one can infer that the PAC becomes more intense as kδ s /δ f increases (see FIG. 3(c)).
We point out that the architecture depicted in FIG. 1(c) constitutes the basic organization of a brain region containing heterogeneous populations. This architecture can be easily extended to account for situations in which more than two frequencies coexist in a single Following previous empirical findings [44], in this simulation we have fixed the slow and fast frequencies to 10 and 80 Hz, respectively. We have also set δs = 0.217 and δ f = 0.038 in the SL model (in blue), while the coupling k is kept unitary. The integration of Eq. (5)-(6) was performed on Matlab with standard variable-step ODE45 solver.
brain region.
To corroborate the claim that our proposed architecture can reproduce PAC that resembles empirical signals in individual brain regions, we resort to in silico experiments. As a common recording technique, local field potentials (LFP) are typically recorded by an electrode placed in the interested brain region to measure the local voltage variation that results from charge separation in the extracellular space. LFP signal analyses revealed that PAC takes place both in the striatum and the hippocampus [44]. Following the method in Ref. [40] (which we elaborate in the Supplementary Material), we generate synthetic LFP signals with diverse PAC intensity. We found that the model proposed in Eq. (3)-(4) in the architecture depicted in FIG. 1(c) can be tuned to be almost indistinguishable, after an initial transient, from synthetic LFP signals. We present an example in FIG. 4.
Because numerous brain functions require synchronous activation of different brain regions, we now extend the single-region architecture discussed above (FIG. 1(c)) to multiple coupled regions. Phase synchrony governs PAC across distant regions -As a building block for more complex network structures, we consider a two-region clique as depicted in FIG. 1(d), wherein interaction across brain regions is established by oscillations in the low-frequency range. Then, we show how cross-region phase synchrony contributes to integrate local high-frequency activities across long distances and to control the direction of information flow between regions.
Under the structure depicted in FIG. 1 (d), the dynamics of the neural populations becomė where the superscript p ∈ {1, 2} is the region index (−p = 2 if p = 1 and −p = 1 otherwise). Notice that, differently from the within-region connections across frequencies, the cross-region connection within the same frequency band is diffusive, which permits synchrony of slow oscillations. Here, k p and k c are the within-and crossregion connection strength, respectively.
We find that the cross-region connection strength k c determines the phase synchrony between the two slow populations. One can reason that if the slow phases are synchronized, the fast amplitudes also behave coherently, provided there is local PAC in each region. Letr 1 f andr 2 f be the average amplitudes in the two regions. To measure the correlation of these two amplitudes, we calculate the phase locking value [40] defined where Φr is the phase angle of the amplitude r. Let ∆ω s = |v 1 s − v 2 s | be the natural frequency of the two slow populations. The ratio k c /∆ω s determines the synchrony of the slow populations. FIG. 5(a) shows that the PLV increases with the ratio k c /∆ω s . Phase synchronization in the lowfrequency range coordinates local high-frequency activities regulated by PAC, which is consistent with empirical studies such as Ref. [26].
Finally, we put forth that our model may shed new light on the directionality of information flow across brain regions. Specifically, we show how natural frequency differences affect the lead-lag relationship of oscillatory rhythms between brain regions. It suffices to investigate the lead-lag direction in the low-frequency range, as we already know that slow oscillations govern communications across regions. Moreover, within each region the slow population remains relatively independent from the fast one. Thus it is sufficient to study the two diffusively coupled slow populations in FIG. 1(d), whose dynamics are approximatelẏ where p = 1, 2. In polar coordinates, z p s = r p s e iθ p s , and Eq. (9) can be rewritten aṡ For a strong connection k c , the system converges to a solution with synchronized frequencies and constant but different amplitudes, i.e.,θ 1 s −θ 2 s = 0,ṙ 1 s = 0 andṙ 2 s = 0. Suppose v 1 s > v 2 s , then we have that z 2 s (t) = r 2 s cos(θ 2 s (t)) + i sin(θ 2 s (t)) = r 2 s r 1 s z 1 s (t − c12 ω ) with ω being the synchronized frequency, which means that z 2 s lags behind z 1 s by τ = c 12 /ω. This behavior may suggest that the former acts as the information sender, and the latter acts as the receiver, since the lead-lag relationship typically encodes the directionality of the information flow [45,46].
To summarize, we postulate that synchronization supports effective communication, and that the natural frequency differences determine the information flow directionality. We conjecture that the brain is capable of controlling this directionality effectively by manipulating the natural frequencies of neural populations. Further experimental studies are needed to investigate these conjectures. We remark that our findings provide some useful insights regarding the inference of directionality in information flow. As low-frequency oscillations are much more engaged in cross-region communication than the high-frequency counterparts, the latter may be filtered out from recorded signals when analyzing information flow between regions.
In this Letter, we have focused on modeling crossfrequency phase-amplitude coupling (PAC). A crossregion structure for interaction of multiple brain regions has also been proposed. The combination of this interconnection structure and our model may pave a way to study other cross-frequency coupling apart from PAC, same-frequency phase synchrony, and their interplay.