# Phase-Amplitude Coupling in Neuronal Oscillator Networks

Yuzhen Qin<EMAIL_ADDRESS>Department of Mechanical Engineering, University
of California, Riverside  Tommaso Menara  Department of Mechanical
Engineering, University of California, Riverside  Danielle S. Bassett
Department of Physics and Astronomy, Department of Bioengineering, Department
of Electrical & Systems Engineering, Department of Neurology, Department of
Psychiatry, University of Pennsylvania & the Santa Fe Institute  Fabio
Pasqualetti<EMAIL_ADDRESS>Department of Mechanical Engineering,
University of California, Riverside

###### Abstract

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 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.

Figure 1: An overview of PAC. (a) Illustration of PAC. (b)-(c) Different
architectures for the emergence of PAC: (b) a fast population with an external
slow input (e.g., sensori stimulus); (c) two interacting populations with
distinct local oscillations (which also constitute a simplified brain region).
(d) Cross-region interconnection structure.

We start by introducing the SL equation:

$\displaystyle\dot{z}=(\sigma+\mathsf{i}\omega-|z|^{2})z,$ (1)

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].

_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., the architecture depicted in FIG. 1(b)). Consider the dynamics of
the neural population with an input $u$

$\displaystyle\dot{z}_{\rm f}=(\delta_{\rm f}+\mathsf{i}\omega_{\rm f}-|z_{\rm
f}|^{2})z_{\rm f}+z_{\rm f}u,$ (2)

where $z_{\rm f}=x_{\rm f}+\mathsf{i}y_{\rm f}$, and $\delta_{\rm f},v_{\rm
f}>0$. The input $u$ can represent various sensory signals from sensory
cortices [38, 25]. We let $u$ be a sinusoidal signal of strength $k_{\rm I}$,
i.e., $u=k_{\rm I}\sin(\omega_{\rm I}t)$. Consequently,

$\displaystyle\dot{z}_{\rm f}=(\delta_{\rm f}+k_{\rm I}\sin(\omega_{\rm
I}t)+\mathsf{i}\omega_{\rm f}-|z_{\rm f}|^{2})z_{\rm f}.$

Then, $\sigma$ in Eq. (1) becomes time-dependent, i.e., $\sigma(t)=\delta_{\rm
f}+k_{\rm I}\sin(\omega_{\rm I}t)$, and the fast amplitude $r_{\rm f}$ tends
to converge to $\sqrt{\sigma(t)}$, yielding PAC.

Figure 2: PAC by exogenous inputs. (a)-(b) A small $k_{\rm I}$ yields
continuous PAC, while a large one yields intermittent PAC ($\delta_{\rm f}=6$,
$k_{\rm I}$ is 3 for panel (a) and 18 for panel (b)). (c) $\delta_{\rm f}$ in
panel (a) is reduced to $0.3$. (d)-(f) The dependence of PAC intensity
(measured by the modulation index (MI)) on the input strength (d), the input
frequency (e), and the fast frequency (f). (If not specified otherwise,
$\delta_{\rm f}=0.3$, $k_{\rm I}=3$ $\omega_{f}=3$ Hz, and $\omega_{\rm
I}=0.5$ Hz.)

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.

_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 are

$\displaystyle\dot{z}_{{\rm s}}=(\delta_{\rm s}+\mathsf{i}v_{\rm s}-|z_{\rm
s}|^{2})z_{\rm s}+z_{s}f(z_{\rm f}),$ (3) $\displaystyle\dot{z}_{{\rm
f}}=(\delta_{\rm f}+\mathsf{i}v_{\rm f}-|z_{\rm f}|^{2})z_{\rm
f}+z_{f}f(z_{\rm s}),$ (4)

where $z_{\rm s}=x_{s}+\mathsf{i}y_{\rm s},z_{\rm f}=x_{\rm
f}+\mathsf{i}y_{\rm f}$, and $\delta_{\rm s},\delta_{\rm f}>0$. The
oscillators’ natural frequencies are determined by $v_{\rm s}$ and $v_{\rm
f}$, with $v_{\rm f}>v_{\rm s}$. The complex function $f(\cdot)$ 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_{\rm f})=0$, then Eq. (4) reduces to Eq. (2) with input $u=f(z_{\rm s})$.
In this case, our analysis above concludes that PAC emerges if $f(z_{\rm 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) become

$\displaystyle\dot{z}_{\rm s}=\big{(}(\delta_{\rm s}+kx_{\rm
f})+\mathsf{i}(v_{\rm s}+ky_{\rm f})-|z_{\rm s}|^{2}\big{)}z_{\rm s},$ (5)
$\displaystyle\dot{z}_{\rm f}=\big{(}(\delta_{\rm f}+kx_{\rm
s})+\mathsf{i}(v_{\rm f}+ky_{\rm s})-|z_{\rm f}|^{2}\big{)}z_{\rm f}.$ (6)

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.

Figure 3: PAC by local interactions. (a) $D_{s}/\delta_{\rm s}$ with
$D_{s}=\sup_{t}||z_{\rm s}(t)-\hat{z}_{\rm s}(t)||$ measures the derivation of
$z_{\rm s}(t)$ from the solution of its decoupled counterpart. As $v_{\rm
f}/v_{\rm s}$ increases, the derivation shrinks dramatically ($k=5,\delta_{\rm
f}=0.3$). (b) PAC emerges by local interaction between two populations as in
FIG. 1 (c), where the slow frequency is 6.5 Hz ($\theta$ band), and the fast
one is 30Hz ($\gamma$ band). (c) PAC intensity (measured by the modulation
index (MI)) increases with $k\delta_{\rm s}/\delta_{\rm f}$. Given
$k\delta_{\rm s}/\delta_{\rm f}$, PAC becomes more intense for larger
$\delta_{\rm f}$.

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.

Figure 4: The SL model precisely reproduces synthetic LFP signals. 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
$\delta_{\text{s}}=0.217$ and $\delta_{\text{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.

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 become

$\displaystyle\dot{z}^{p}_{{\rm s}}=(\delta^{p}_{\rm s}+\mathsf{i}v^{p}_{\rm
s}-|z^{p}_{\rm s}|^{2})z_{\rm s}+k^{p}z^{p}_{\rm s}z^{p}_{\rm
f}+k_{c}(z^{-p}_{{\rm s}}-z^{p}_{{\rm s}}),$ (7)
$\displaystyle\dot{z}^{p}_{{\rm f}}=(\delta^{p}_{\rm f}+\mathsf{i}v^{p}_{\rm
f}-|z^{p}_{\rm f}|^{2})z_{\rm f}+k^{p}z^{p}_{\rm f}z^{p}_{\rm s},$ (8)

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].

Figure 5: Roles of Low-frequency phase synchrony. (a) PLV of the fast
amplitudes in the two regions against the ratio $K_{c}/\Delta\omega_{\rm s}$.
The line indicated the mean of 1000 random realizations contained in the
shaded area. It can be inferred that the fast populations becomes more
correlated as the synchrony of the slow oscillations increases. (b)-(c) Lead-
lag relationship between the two signals: $z^{1}_{\rm s}$ precedes $z^{2}_{\rm
s}$ if $v^{1}_{\rm s}>v^{2}_{\rm s}$ and vice versa.

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

$\displaystyle\dot{z}^{p}_{{\rm s}}=(\delta^{p}_{\rm s}+\mathsf{i}v^{p}_{\rm
s}-|z^{p}_{\rm s}|^{2})z_{\rm s}+k_{c}(z^{-p}_{{\rm s}}-z^{p}_{{\rm s}}),$ (9)

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

$\displaystyle\dot{r}^{p}_{{\rm s}}=(\delta^{p}_{\rm s}-({r^{p}_{{\rm
s}}})^{2})r^{p}_{{\rm s}}+k_{c}(r^{-p}_{{\rm s}}\cos(\theta^{-p}_{{\rm
s}}-\theta^{p}_{{\rm s}})-r^{-p}_{{\rm s}}),$
$\displaystyle\dot{\theta}^{p}_{{\rm s}}=v^{p}_{{\rm
s}}+k_{c}\frac{r^{-p}_{{\rm s}}}{r^{p}_{{\rm s}}}\sin(\theta^{-p}_{{\rm
s}}-\theta^{p}_{{\rm s}}).$

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

$\displaystyle v^{1}_{{\rm s}}-v^{2}_{{\rm s}}-k_{c}\left(\frac{r^{2}_{{\rm
s}}}{r^{1}_{{\rm s}}}+\frac{r^{1}_{{\rm s}}}{r^{2}_{{\rm s}}}\right)\sin
c_{12}=0,$

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.

###### Acknowledgements.

This work was supported in part by awards ARO W911NF-18-1-0213, ARO
71603NSYIP, and NSF NCS-FO-1926829.

## References

  * Palva and Palva [2012] S. Palva and J. M. Palva, Trends in Cognitive Sciences 16, 219 (2012).
  * Tass et al. [2003] P. Tass, T. Fieseler, J. Dammers, K. Dolan, P. Morosan, M. Majtanik, F. Boers, A. Muren, K. Zilles, and G. Fink, Physical Review Letters 90, 088101 (2003).
  * Nicosia et al. [2013] V. Nicosia, M. Valencia, M. Chavez, A. Díaz-Guilera, and V. Latora, Physical Review Letters 110, 174102 (2013).
  * Qin et al. [2020] Y. Qin, M. Cao, B. D. Anderson, D. S. Bassett, and F. Pasqualetti, arXiv preprin arXiv:2003.07405 (2020).
  * Menara et al. [2020] T. Menara, G. Baggio, D. S. Bassett, and F. Pasqualetti, IEEE Transactions on Control of Network Systems 7, 302 (2020).
  * Canolty et al. [2006] R. T. Canolty, E. Edwards, S. S. Dalal, M. Soltani, S. S. Nagarajan, H. E. Kirsch, M. S. Berger, N. M. Barbaro, and R. T. Knight, Science 313, 1626 (2006).
  * Jensen and Colgin [2007] O. Jensen and L. L. Colgin, Trends in Cognitive Sciences 11, 267 (2007).
  * Canolty and Knight [2010] R. T. Canolty and R. T. Knight, Trends in Cognitive Sciences 14, 506 (2010).
  * Lakatos et al. [2005] P. Lakatos, A. S. Shah, K. H. Knuth, I. Ulbert, G. Karmos, and C. E. Schroeder, Journal of Neurophysiology 94, 1904 (2005).
  * Voloh et al. [2015] B. Voloh, T. A. Valiante, S. Everling, and T. Womelsdorf, Proceedings of the National Academy of Sciences 112, 8457 (2015).
  * Cohen et al. [2009] M. X. Cohen, N. Axmacher, D. Lenartz, C. E. Elger, V. Sturm, and T. E. Schlaepfer, Journal of Cognitive Neuroscience 21, 875 (2009).
  * Axmacher et al. [2010] N. Axmacher, M. M. Henseler, O. Jensen, I. Weinreich, C. E. Elger, and J. Fell, Proceedings of the National Academy of Sciences 107, 3228 (2010).
  * Yanagisawa et al. [2012] T. Yanagisawa, O. Yamashita, M. Hirata, H. Kishima, Y. Saitoh, T. Goto, T. Yoshimine, and Y. Kamitani, Journal of Neuroscience 32, 15467 (2012).
  * Roux and Uhlhaas [2014] F. Roux and P. J. Uhlhaas, Trends in Cognitive Sciences 18, 16 (2014).
  * Lakatos et al. [2008] P. Lakatos, G. Karmos, A. D. Mehta, I. Ulbert, and C. E. Schroeder, Science 320, 110 (2008).
  * Tort et al. [2009] A. B. Tort, R. W. Komorowski, J. R. Manns, N. J. Kopell, and H. Eichenbaum, Proceedings of the National Academy of Sciences 106, 20942 (2009).
  * White et al. [2000] J. A. White, M. I. Banks, R. A. Pearce, and N. J. Kopell, Proceedings of the National Academy of Sciences 97, 8128 (2000).
  * Tort et al. [2007] A. B. Tort, H. G. Rotstein, T. Dugladze, T. Gloveli, and N. J. Kopell, Proceedings of the National Academy of Sciences 104, 13490 (2007).
  * Wulff et al. [2009] P. Wulff, A. A. Ponomarenko, M. Bartos, T. M. Korotkova, E. C. Fuchs, F. Bähner, M. Both, A. B. Tort, N. J. Kopell, W. Wisden, et al., Proceedings of the National Academy of Sciences 106, 3561 (2009).
  * Fontolan et al. [2013] L. Fontolan, M. Krupa, A. Hyafil, and B. Gutkin, The Journal of Mathematical Neuroscience 3, 16 (2013).
  * Akam et al. [2012] T. Akam, I. Oren, L. Mantoan, E. Ferenczi, and D. M. Kullmann, Nature Neuroscience 15, 763 (2012).
  * Onslow et al. [2014] A. C. Onslow, M. W. Jones, and R. Bogacz, PloS one 9, e102591 (2014).
  * Jedynak et al. [2015] M. Jedynak, A. J. Pons, and J. Garcia-Ojalvo, Frontiers in computational neuroscience 9, 14 (2015).
  * Chehelcheraghi et al. [2017] M. Chehelcheraghi, C. van Leeuwen, E. Steur, and C. Nakatani, PloS one 12, e0173776 (2017).
  * Hyafil et al. [2015a] A. Hyafil, A.-L. Giraud, L. Fontolan, and B. Gutkin, Trends in Neurosciences 38, 725 (2015a).
  * von Nicolai et al. [2014] C. von Nicolai, G. Engler, A. Sharott, A. K. Engel, C. K. Moll, and M. Siegel, Journal of Neuroscience 34, 5938 (2014).
  * Lara and Wallis [2014] A. H. Lara and J. D. Wallis, Nature neuroscience 17, 876 (2014).
  * Daume et al. [2017] J. Daume, T. Gruber, A. K. Engel, and U. Friese, Journal of Neuroscience 37, 313 (2017).
  * Von Stein and Sarnthein [2000] A. Von Stein and J. Sarnthein, International Journal of Psychophysiology 38, 301 (2000).
  * Montbrió and Pazó [2018] E. Montbrió and D. Pazó, Physical Review Letters 120, 244101 (2018).
  * Qin et al. [2019] Y. Qin, Y. Kawano, O. Portoles, and M. Cao, arXiv preprint arXiv:1906.01065 (2019).
  * Menara et al. [2019] T. Menara, G. Baggio, D. S. Bassett, and F. Pasqualetti, in _IEEE Conf. on Decision and Control_ (Nice, France, 2019), pp. 4697–4704.
  * Salova and D’Souza [2020] A. Salova and R. M. D’Souza, arXiv preprint arXiv:2006.06163 (2020).
  * Matheny et al. [2019] M. H. Matheny, J. Emenheiser, W. Fon, A. Chapman, A. Salova, M. Rohden, J. Li, M. H. de Badyn, M. Pósfai, L. Duenas-Osorio, et al., Science 363 (2019).
  * Schuster and Wagner [1990] H. G. Schuster and P. Wagner, Biological Cybernetics 64, 77 (1990).
  * Brunel and Hakim [1999] N. Brunel and V. Hakim, Neural Computation 11, 1621 (1999).
  * Börgers and Kopell [2005] C. Börgers and N. Kopell, Neural Computation 17, 557 (2005).
  * Mazzoni et al. [2008] A. Mazzoni, S. Panzeri, N. K. Logothetis, and N. Brunel, PLoS Comput Biol 4, e1000239 (2008).
  * Provansal et al. [1987] M. Provansal, C. Mathis, and L. Boyer, Journal of Fluid Mechanics 182, 1 (1987).
  * Tort et al. [2010] A. B. L. Tort, R. Komorowski, H. Eichenbaum, and N. Kopell, Journal of Neurophysiology 104, 1195 (2010).
  * Whittington and Traub [2003] M. A. Whittington and R. D. Traub, Trends in Neurosciences 26, 676 (2003).
  * Hyafil et al. [2015b] A. Hyafil, L. Fontolan, C. Kabdebon, B. Gutkin, and A.-L. Giraud, elife 4, e06213 (2015b).
  * Sanders et al. [2007] J. A. Sanders, F. Verhulst, and J. Murdock, _Averaging methods in nonlinear dynamical systems_ , vol. 59 (Springer, 2007).
  * Tort et al. [2008] A. B. Tort, M. A. Kramer, C. Thorn, D. J. Gibson, Y. Kubota, A. M. Graybiel, and N. J. Kopell, Proceedings of the National Academy of Sciences 105, 20517 (2008).
  * Nolte et al. [2008] G. Nolte, A. Ziehe, V. V. Nikulin, A. Schlögl, N. Krämer, T. Brismar, and K.-R. Müller, Physical Review Letters 100, 234101 (2008).
  * Moon et al. [2015] J.-Y. Moon, U. Lee, S. Blain-Moraes, and G. A. Mashour, PLoS Comput Biol 11, e1004225 (2015).

##  I Supplemental Material

Modulation Index (MI). To measure the intensity of phase-amplitude coupling
(PAC), we employ the Modulation Index (MI) [40]. In the following, we
summarize the steps to compute the MI.

We first derive the phase-amplitude distribution $P$ of two complex signals –
a slow one $z_{\rm s}(t)=x_{\rm s}(t)+\mathsf{i}y_{\rm s}(t)$, and fast one
$z_{\rm f}(t)=x_{\rm f}(t)+\mathsf{i}y_{\rm f}(t)$ – as follows:

  1. (1)

We extract the time series of the phases of $z_{\rm s}(t)$. For each time
instant $t$, the phase is $\theta_{\rm s}(t)$ in $z_{\rm s}(t)=r_{\rm
s}(t)e^{\mathsf{i}\theta_{\rm s}(t)}$.

  2. (2)

Extract the time series of the amplitude envelope of $z_{\rm f}(t)$.
Analogously, the amplitude at each time instant $t$ is $r_{\rm f}(t)$ in
$z_{\rm f}(t)=r_{\rm f}(t)e^{\mathsf{i}\theta_{\rm f}(t)}$.

  3. (3)

Divide the phases $\theta_{\rm s}(t)$ into $n$ bins, where $n$ can be any
positive number (we use 18 in the our text). For each bin $i$, calculate the
mean of $r_{\rm f}$, denoted by $\langle r_{\rm f}\rangle_{i}$.

  4. (4)

Normalize the mean amplitude in each bin by

$\displaystyle P(j)=\frac{\langle r_{\rm f}\rangle_{i}}{\sum_{j}^{n}\langle
r_{\rm f}\rangle_{j}}.$

Notice that the phase-amplitude distribution satisfies $P(i)\geq 0$ for any
$i$ and $\sum_{j=1}^{n}P(j)=1$, akin to a probability mass function. The
intensity of PAC can be inferred intuitively by visual inspection of the plot
of $P$ over all the bins. Nevertheless, to provide a quantitative measure of
PAC intensity, the Modulation Index is defined as

$\displaystyle{\rm MI}=\frac{D_{KL}(P,U)}{\log(n)}$

where $D_{KL}(P,U)=\log(n)-H(P)$, $H(P)=-\sum_{j=1}^{n}P(j)\log(P(j))$, and
$U$ represent the uniform distribution. Loosely speaking, the Modulation Index
is a scaled Kullback-Leibler distance of the phase-amplitude distribution from
the uniform distribution, satisfying $0\leq\text{MI}\leq 1$. $\text{MI}=1$
means that $P$ is a Dirac-like distribution, while $\text{MI}=0$ is equivalent
to a uniform distribution. Thus, for any signal, the larger the MI is, the
more intense the PAC.

Generating synthetic local field potential (LFP). Following [40], we generate
synthetic LFP signals as

$\displaystyle x(t)=r_{\rm f}(t)\sin(2\pi\omega_{\rm f}t)+r_{\rm
s}\sin(2\pi\omega_{\rm s}t)+W(t),$

where $r_{\rm f}(t)$ is the amplitude envelope of the fast oscillation,
$r_{\rm s}$ is a constant representing the amplitude of the slow oscillation,
and $W(t)$ is a Gaussian white noise. To ensure that the synthetic signal
contains phase-amplitude coupling, $r_{\rm f}(t)$ is defined as

$\displaystyle r_{\rm f}(t)=\bar{r}_{\rm f}\frac{\mu\sin(2\pi\omega_{\rm
s}t)+2-\mu}{2},$

where $\bar{r}_{\rm f}(t)$ is a constant determining the maximum amplitude of
the fast oscillation, and $\mu$ is the fraction of the amplitude envelope
modulated by the slow oscillation. The PAC intensity is thus controlled by the
parameter $\mu$. LFP signals with diverse PAC intensity can be generated by
manipulating $\mu$.