Phase-Amplitude Coupling in Neuronal Oscillator Networks

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 cross-frequency 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 low-frequency (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 $\gamma$ rhythms (>30 Hz) between the cortex and striatum are typically established by cross-regional phase synchronization in the $\theta$ band (4-8 Hz) concurrently with $\theta$-$\gamma$ 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 [29, 26], 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+\mathsf{i}y$), with $\mathsf{i}=\sqrt{-1}$, and $\sigma$ and $\omega$ are real. The SL oscillator is also known as a Hopf oscillator, because, for $\sigma<0$, $z=0$ is globally stable, and for $\sigma>0$ the model undergoes a Hopf bifurcation with $z$ converging to a circular trajectory (limit cycle) of radius $\sqrt{\sigma}$ and angular frequency $\omega$. The SL oscillator can equivalently be described in polar coordinates by letting $z=re^{\mathsf{i}\theta}$, where $r$ and $\theta$ 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 [30, 31, 32, 4, 5], 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 $\sigma$ in Eq. (1) from negative to positive. The oscillation frequency mainly depends on the intrinsic decay time of inhibition in the population [37].

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_{\rm I}<\delta_{\rm f}$, $\sigma(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_{\rm I}>\delta_{\rm f}$, $\sigma(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_{\rm f}$ tends towards $0$, and the fast oscillation may disappear. We illustrate intermittent PAC for large $k_{\rm 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 $\delta_{\rm f}$ in Eq. (2). As pointed out in Ref. [39], since $\sigma$ in Eq. (1) determines the convergence rate to the limit cycle, we then reason that the rate of $r_{\rm f}$ tracking $\sqrt{\sigma(t)}$ is also positively correlated with $\delta_{\rm f}$. Therefore, a larger $\delta_{\rm f}$ means a smaller phase lag. FIG. 2(c) demonstrates that a noticeable phase lag appears for a smaller $\delta_{\rm 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_{\rm f}$ tends to track $\sqrt{\sigma(t)}$=$\sqrt{\delta_{\rm f}}\cdot\sqrt{1+{k_{\rm I}}/{\delta_{\rm f}}\sin(\omega_{\rm I}t)}$, it is clear that a larger ratio ${k_{\rm I}}/{\delta_{\rm f}}$ implies larger fluctuations in $r_{\rm f}$ than a weaker one. That is, the PAC intensity increases with ${k_{\rm I}}/{\delta_{\rm f}}$ (see FIG. 2(d)). Additionally, for a fixed a ratio ${k_{\rm I}}/{\delta_{\rm f}}$, variations of $\delta_{\rm f}$ impact the fluctuations of $r_{\rm f}$: a larger $\delta_{\rm f}$ implies more intense PAC (see FIG. 2(d)).

The input frequency $\omega_{\rm I}$ also has a profound impact on the PAC intensity. A smaller $\omega_{\rm I}$ means a more slowly varying $\sigma(t)$, which is easier for the fast amplitude $r_{\rm f}$ to track. FIG. 2(e) shows that the PAC intensity decreases as $\omega_{\rm I}$ increases. By contrast, the fast frequency $\omega_{\rm 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.

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) become

In polar coordinates, $z_{\rm s}=r_{\rm s}e^{\mathsf{i}\theta_{\rm s}}$ and $z_{\rm f}=r_{\rm f}e^{\mathsf{i}\theta_{\rm f}}$. It can be observed that the terms $\sigma$ and $\omega$ (see Eq. (1)) that describe the amplitude and phase in the slow and fast population read as $\sigma_{\rm s}=\delta_{\rm s}+kx_{\rm f}$, $\omega_{\rm s}=v_{\rm s}+ky_{\rm f}$, and $\sigma_{\rm f}=\delta_{\rm f}+kx_{\rm s}$, $\omega_{\rm f}=v_{\rm f}+ky_{\rm 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_{\rm f}>v_{\rm s}$, one can see that $z_{\rm f}$ oscillates faster than the $z_{\rm s}$. Applying averaging techniques developed in Ref. [43] to the dynamics in Eq. (5), we obtain the approximate guiding system $\dot{\hat{z}}_{\rm s}=(\delta_{\rm s}+\mathsf{i}v_{\rm s}-|\hat{z}_{\rm s}|^{2})\hat{z}_{\rm s}$, which corresponds to the dynamics of the decoupled slow population. The solutions to Eq. (5) and the guiding system satisfy $\|z_{\rm s}(t)-\hat{z}_{\rm s}(t)\|\leq c{v_{s}}/{v_{f}}$ for some positive $c$. Loosely speaking, the slow trajectory ${z}_{\rm s}(t)$ fluctuates in the neighborhood of the trajectory of its decoupled counterpart $\hat{z}_{\rm s}(t)$. Such fluctuations decrease as the two frequencies separate (see Fig. 3(a)). In contrast, Eq. (6) can be rewritten as $\dot{z}_{\rm f}=\big{(}(\delta_{\rm f}+kr_{\rm s}\cos\theta_{\rm s})+\mathsf{i}(v_{\rm f}+r_{\rm s}\sin\theta_{\rm s})-|z_{\rm f}|^{2}\big{)}z_{\rm f}$. Similar to our earlier analysis, the fast amplitude $r_{\rm f}$ tends to track $\sqrt{\delta_{\rm f}+kr_{\rm s}\cos\theta_{\rm s}}$, which implies PAC.

Let the fast population in FIG. 1(c) oscillate in the $\gamma$ band, and the slow population oscillate in the $\theta$ 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_{\rm s}$, which is determined by $\delta_{\rm s}$. Note that the connection strength $k$ also affects PAC intensity. Therefore, one can infer that the PAC becomes more intense as $k\delta_{\rm s}/\delta_{\rm 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 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.

Under the structure depicted in FIG. 1 (d), the dynamics of the neural populations become

where the superscript $p\in\{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 cross-region 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. Let $\bar{r}^{1}_{{\rm f}}$ and $\bar{r}^{2}_{{\rm 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 as ${\rm PLV}=\frac{1}{T}\sum_{t=1}^{T}e^{\mathsf{i}(\Phi{\bar{r}^{1}_{{\rm f}}}(t)-\Phi{\bar{r}^{2}_{{\rm f}}}(t))},$ where $\Phi{r}$ is the phase angle of the amplitude $r$. Let $\Delta\omega_{\rm s}=|v^{1}_{\rm s}-v^{2}_{\rm s}|$ be the natural frequency of the two slow populations. The ratio $k_{c}/\Delta\omega_{\rm s}$ determines the synchrony of the slow populations. FIG. 5(a) shows that the PLV increases with the ratio $k_{c}/\Delta\omega_{\rm s}$. Phase synchronization in the low-frequency 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 approximately

where $p=1,2$. In polar coordinates, $z^{p}_{{\rm s}}=r^{p}_{{\rm s}}e^{i\theta^{p}_{{\rm s}}}$, and Eq. (9) can be rewritten as

For a strong connection $k_{c}$, the system converges to a solution with synchronized frequencies and constant but different amplitudes, i.e., $\dot{\theta}^{1}_{{\rm s}}-\dot{\theta}^{2}_{{\rm s}}=0$, $\dot{r}^{1}_{{\rm s}}=0$ and $\dot{r}^{2}_{{\rm s}}=0$. Letting $\theta^{1}_{{\rm s}}-\theta^{2}_{{\rm s}}=c_{12}$ yields

which has a solution $c_{12}=\arcsin(\frac{v^{1}_{{\rm s}}-v^{2}_{{\rm s}}}{k_{c}\lambda})$ with $\lambda=\frac{r^{2}_{{\rm s}}}{r^{1}_{{\rm s}}}+\frac{r^{1}_{{\rm s}}}{r^{2}_{{\rm s}}}>0$ (the other solution $c_{12}=\pi-\arcsin(\frac{v^{1}_{{\rm s}}-v^{2}_{{\rm s}}}{k_{c}\lambda})$ is unstable). The system converges to a solution wherein the phase $\theta^{2}_{{\rm s}}$ lags behind (resp. leads) $\theta^{1}_{{\rm s}}$ if $v^{1}_{{\rm s}}>v^{2}_{{\rm s}}$ (resp. $v^{1}_{{\rm s}}<v^{2}_{{\rm s}}$), which we illustrate in FIG. 5(b)-(c). Suppose $v^{1}_{{\rm s}}>v^{2}_{{\rm s}}$, then we have that $z^{2}_{\rm s}(t)=r^{2}_{\rm s}\cos(\theta^{2}_{\rm s}(t))+\mathsf{i}\sin(\theta^{2}_{\rm s}(t))=\frac{r^{2}_{\rm s}}{r^{1}_{\rm s}}z^{1}_{\rm s}(t-\frac{c_{12}}{\omega})$ with $\omega$ being the synchronized frequency, which means that $z^{2}_{\rm s}$ lags behind $z^{1}_{\rm s}$ by $\tau={c_{12}}/{\omega}$. 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 cross-frequency phase-amplitude coupling (PAC). A cross-region 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.