A neural mass model of cross frequency coupling

Electrophysiological signals of cortical activity show a range of possible frequency and amplitude modulations, both within and across regions, collectively known as cross-frequency coupling. To investigate whether these modulations could be considered as manifestations of the same underlying mechanism, we developed a neural mass model. The model provides five out of the theoretically proposed six different coupling types. Within model components, slow and fast activity engage in phase-frequency coupling in conditions of low ambient noise level and with high noise level engage in phase-amplitude coupling. Between model components, these couplings can be coordinated via slow activity, giving rise to more complex modulations. The model, thus, provides a coherent account of cross-frequency coupling, both within and between components, with which regional and cross-regional frequency and amplitude modulations could be addressed.


Introduction
Electrophysiological activity of interconnected neurons encompasses multiple oscillatory components [1,2]; these are subject to modulation of frequency [3][4][5] and amplitude [6,7]. Modulation may have various utilities [8,9], such as sequence encoding [10,11], rectification of local neural activity [12,13] and long-rage information transfer [14,15]. The mechanisms behind these effects, however, are still not well-understood. In theory, they could all be manifestations of a single underlying principle. To contemplate this possibility, we propose a model that describes frequency and amplitude modulations as systematic relationships between slow and fast oscillatory components of neural population activity.
Systematic relationships between oscillatory components are collectively known as crossfrequency coupling (CFC). Jensen & Colgin [2] listed six types of CFC of interest to electrophysiology: phase-phase, phase-frequency, phase-amplitude, frequency-frequency, amplitudeamplitude and amplitude-frequency couplings (PPC, PFC, PAC, FFC, AAC, and AFC, respectively). Empirical studies have mostly observed PPC [16][17][18][19][20], PAC [21][22][23] and occasionally AAC [24][25][26]. To our knowledge, PFC, AFC and FFC have not been empirically observed so far. Shortcomings of current CFC detection methods may, in part, be responsible for this. In simulation. In both simulations, the second component has zero mean noise level, thereby preventing any fast oscillations to emerge within the second component. Thus, cross regional PAC and PFC arise via synchronization with the slow oscillations that connect the two components. In the third, fourth, and fifth simulations, FFC, AAC and AFC were generated, respectively, each as a function of the levels of nonzero mean noise input to both components. Based on the results, we propose the present model as an inclusive theoretical framework for the different CFC phenomena.

Methods Overview
As illustrated in Fig 1, the present model consists of two identical components labelled as Node 1 and Node 2. Each component is a four-populations NMM node which simulates electrophysiological activity [36]. The four populations represent, respectively, pyramidal neuron, excitatory interneuron, slow inhibitory interneuron, and fast inhibitory interneuron of the cortex. In each node, interaction between the pyramidal neurons and slow inhibitory interneurons generates a steady-state activity, which is considered as its main (intrinsic) oscillation (e.g. alpha band activity in visual cortex in [33]. Unique to our four-populations model, the fast-inhibitory interneuron population has dynamic self-feedback. The feedback enables the fast inhibitory interneurons to generate oscillations in the gamma range. This fast activity can Each node is identical to the NMM in [36]. The nodes are mutually connected by the dynamics of the pyramidal neurons. Two independent Gaussian noise sources drive the NMM nodes. https://doi.org/10.1371/journal.pone.0173776.g001 occur simultaneously with the slow activity, which is a prerequisite for CFC in the model. The two four-populations nodes are mutually connected via two between-node pyramidal neuron populations. The emergence and stability of fast activity in each node depends on two parameters, P 1 and P 2 , which represent the mean levels of noise input to the nodes [36]. Small values deactivate the fast oscillations (damping). Values above the damping regime determine the type of modulation between the slow and fast oscillation. For parameter values below a critical threshold, resonance occurs, allowing frequency modulation, namely, PFC, AFC, and FFC; above the threshold the fast oscillations reach a steady-state, in which the frequency is invariant, and amplitude modulation results, allowing for PAC and AAC. For an analytical account of this mechanism, see [44]). We systematically vary the values of P 1 and P 2 in order to demonstrate that PFC, PAC, FFC, AAC and AFC appear as a function of the parameters P 1 and P 2 , without changing model structure or connectivity.

Model
Each of the four neural population units of a node, pyramidal neurons, excitatory interneurons, slow inhibitory interneurons, and fast inhibitory interneurons, has two variables: mean firing rate fr(t) which is a function of mean membrane potential v m (t), and mean post synaptic potential v(t). The mean membrane potential is obtained from the weighted summation of inputs. The mean membrane potential is converted to the mean firing rate via a sigmoid function: Here v max is the maximum firing rate of the population of neurons, v θ is the value of the potential for which a 50% mean firing rate is achieved, and r is the slope of the sigmoid at v θ ; The parameter v θ will serve as mean firing threshold. The sigmoid function is parametrized equally for all populations in the model. A second-order differential equation called the dendritic transfer function relates the mean firing rate, fr(t), of a population to the un-weighted mean post synaptic potential, v(t), of that population via the dendritic transfer function: where G and ω characterize the mean strength and speed, respectively. Different population units are distinguished by their time constant T ¼ 1 o ½sec and the low frequency gain of the dendritic transfer function G o ½mV:sec. The dendritic transfer function of each population unit is denoted as H lx (S), in which l 2 {p: pyramidal neuron population, q: excitatory interneuron population, f: fast inhibitory interneuron population, s: slow inhibitory interneuron popula-tion} and x 2 {1,2} is the node number, e.g. H s2 (S) denotes the dendritic transfer function of the slow inhibitory interneuron population in the second node.
Four interconnected population units constitute a single node. Inputs to each population unit are the weighted outputs from the other units, reflecting the mean excitatory or inhibitory post synaptic potentials of an excitatory or inhibitory population (EPSP or IPSP). The input synaptic potentials are weighted by constant connectivity strengths which are called synaptic gains. A synaptic gain between populations is represented by C αβ where α, β 2 {p: pyramidal neuron population, q: excitatory interneuron population, f: fast inhibitory interneuron population, s: slow inhibitory interneuron population}. Here α and β refer to target and source population, respectively; e.g. the constant C ps is the gain of the synaptic connection from the slow inhibitory population to the pyramidal neuron population.
The pyramidal neuron population unit excites all three interneuron population units within the node and in turn receives input from them. The slow inhibitory interneuron population unit inhibits the fast inhibitory interneuron population unit. The fast inhibitory interneuron population has an extra state variable, v ff (t), for the dynamic self-inhibition. The dynamics of the self-feedback in the fast inhibitory population is of first order and serves as a low pass filter with cut-off frequency f c ¼ 1 ð2pt f Þ . As shown in Fig 1, Nodes 1 and 2 have identical structure. The two nodes are mutually connected via the pyramidal populations. The synaptic gain of Node 1 to Node 2 is represented by K 21 , while that in the inverse direction is given as K 12 . The dendritic transfer functions of pyramidal neurons between the nodes are labeled as H p21 (S) and H p12 (S). For each node, independent Gaussian noise excites both the fast inhibitory interneuron and pyramidal neuron units. The noise represents the mean firing rate of unmodeled external neural populations. The noise is passed through the dendritic transfer function of pyramidal neurons between the nodes, i.e., H p21 (S) and H p12 (S), in order to convert it to mean post synaptic potentials. The potentials are weighted by factors K p and K f in order to excite the pyramidal neurons and fast inhibitory interneurons, respectively. For each node, the mean membrane potential of the pyramidal population is the main output. The pyramidal population unit has been the main output (and input) unit in previous cortical neural mass models [29,35]. The mean membrane potential of the neurons represents population activity of the node; thus, this output is comparable with macroscopic electrophysiological brain signals, such as, LFP, ECoG, and EEG. The dynamics of the full model is described as follows.
Equations for Node 1. Pyramidal neuron population: Excitatory interneuron population: Slow inhibitory interneuron population: Fast inhibitory interneuron population: External inputs: Equations for Node 2. Pyramidal neuron population: Excitatory interneuron population: Slow inhibitory interneuron population: Fast inhibitory interneuron population: External inputs: Parameters. Table 1 shows the parameter values that are common to all current CFC simulations. The two nodes are identically parametrized, except for the time constant of the dynamic self-feedbacks in the fast inhibitory interneurons, τ f1 and τ f2 . As we will explain further in the subsection Choice of Slow and Fast Oscillation Frequencies, these parameters were fixed at different values for the sake of obtaining FFC, AAC and AFC. The remaining parameter values are those suggested by [33] and [30], same as in our previous study [36]. Jansen and Rit [33] assumed that synaptic connections are active in equal proportions between interacting units, and proposed to the fix ratios of synaptic gains. The set of parameters was modified by Wendling et al., [30] to model fast inhibitory interneurons. The modified set of parameters has been the basis set for fast oscillations in many neural mass modeling studies, including the current one.
Noise sources and fast oscillations. Table 2 shows how the values of the main parameters P 1 and P 2 were systematically varied in order to enable different types of CFCs. Two independent noise sources, NðP 1 ; s 2 1 Þ and NðP 2 ; s 2 2 Þ, excite both the pyramidal neurons and the fast

Parameter Interpretation Notation Value
Synaptic gain from excitatory interneurons to pyramidal neurons C qp1 , C qp2 135 Synaptic gain from pyramidal neurons to excitatory interneurons C pq1 , C pq2 108 Synaptic gain from slow inhibitory interneurons to pyramidal neurons C sp1 , C sp2 33.75 Synaptic gain from pyramidal neurons to slow inhibitory interneurons C ps1 , C ps2 33.75 Synaptic gain from fast inhibitory interneurons to pyramidal neurons C fp1 , C fp2 40.5 Synaptic gain from pyramidal neurons to fast inhibitory interneurons C pf1 , C pf2 27 Synaptic gain from fast inhibitory interneurons to slow inhibitory interneurons C fs1 , C fs2 10.8 Synaptic gain of fast inhibitory interneurons self-feedback C ff1 , C ff2 135 Noise excitation weight for pyramidal neurons K p1 , K p2 40 Noise excitation weight for fast inhibitory interneurons K f1 , K f2 108 Average time constant of pyramidal neurons membrane potential, the inverse divided by 2π is equivalent to low cut-off frequency [sec] Average time constant of between node pyramidal neurons membrane potential, the inverse divided by 2π is equivalent to low cut-off frequency [sec] inhibitory interneurons of Nodes 1 and 2, respectively. The noise sources represent the unmodeled mean firing rates, which originate from various cortical and subcortical regions. As long as the mean noise level P is too small, fast oscillations fail to arise due to damping. Higher levels of noise trigger fast oscillations, for which P serves as a bifurcation parameter. We will refer to the bifurcation value as limit-cycle threshold. When a P parameter value is below the limitcycle threshold, the fast inhibitory population with feedback serves as a resonance filter. Such a filter is selective to a specific frequency in the input and the output contains the resonance frequency component. In this regime, the frequency of the fast oscillation is sensitive to the slow oscillations in the modulating input. Hence, the resonance is suitable as a frequency modulation mechanism. For PFC and FFC simulations, the parameter was set in the resonance range. When P value is above the threshold, a limit cycle emerges corresponding to the fast oscillation. In this type of oscillation, the amplitude of the fast oscillation is sensitive to modulatory input. Therefore, the limit cycle oscillation is a suitable mechanism for amplitude modulation. For PAC and AAC simulations, thus, the P parameters are set to obtain a limit-cycle osillation. For AFC, P values were set to obtain a limit cycle in one node and resonance oscillations in the other. Slow oscillations, generation and synchronization. As in previous NMM studies [29,45], the model generates slow oscillations through the interactions between the pyramidal and slow inhibitory interneuron populations within each node. When nodes are connected, the slow oscillations synchronize, desynchronize, or show intermittent synchronization behavior, depending on the coupling strengths K 12 and K 21 between them [46]. Synchronization will not arise for small values of K 12 and K 21 . We used K 12 = K 21 = 40 throughout this study, which assures that the slow oscillations are synchronized.
Choice of slow and fast oscillation frequencies. For frequency modulation, intuitively, the condition f fast >> f slow should be met where f fast is the frequency of the fast oscillation and f slow is that of the slow one. We set the slow oscillation in the delta band and the fast oscillation in the gamma band, the slowest and the fastest commonly observed in EEG. To generate slow oscillation in the delta band (~3Hz), the time constants of the pyramidal neurons in both nodes were set to o p1 ¼ o p2 ¼ 10 rad s . This value differs by a factor 10 from that used in previous studies to obtain alpha oscillations [29,30,36]. In order to keep the stationary conditions of the model unaffected, the low frequency gain ( G o ) was kept in the same proportion, meaning G was changed from 3.2 (for alpha) to 0.32 (for delta). The frequency of the fast oscillations is tuned through the time constant of the dynamic feedback, τ f1 and τ f2 . These were set to 0.5 ms and 1ms, respectively, resulting in 51Hz and 42Hz gamma activity in the first and the second node, respectively. These frequencies were chosen to be different in order to show that for FFC, AAC, and AFC in general the fast oscillations do not need to have the same frequencies.
For amplitude modulation, the condition should be met that f fast f slow > 2 [27]. The fast and slow frequencies chosen above meet these conditions.

Evaluation of Cross-Frequency Coupling (CFC)
To assess the CFC generated by our model, the output of each node was decomposed into fast-oscillation (FO) and slow-oscillation (SO) signals, respectively, by high-pass or low-pass filtering at a 15Hz cut-off frequency (f c ). For filtering we used the eeglab13.4.4b Matlab toolbox [47]. FO and SO signals were always normalized to zero mean and unit variance. To evaluate frequency modulation, the average number of zero crossings of the FO signal was calculated for each positive and negative phase of the SO signal in the same node. The average number of zero crossings was obtained by calculating the zero-crossing rate according to the following equation: in which T is the time interval of SO in negative (from sine(-π) to sine(0)) or positive (from sine(0) to sine(π)) phases and f s is the sampling frequency. The indicator function I{A} is 1 if its argument is true and 0 otherwise. The average over T does not reflect small fluctuations, which do not cross zero. As a result, the number zcr tends to be lower than the actual frequency, i.e., the average numer of zero crossing is not identical to the average frequency, but may serve as an index is of it. The power spectral density (PSD) of the output signal was calculated using the [48] method. A Hanning time window of length equal to 16348 data point was moved stepwise with 25% overlap. Amplitude modulation was evaluated by visual inspection, comparing the envelope of the FO signals with the SO signal.

Simulation specifications
The model was implemented and ran in Matlab/Simulink 8.0 (R2012b), using a 4 th order Runge-Kutta method with fixed time step of 0.0001 seconds. The number of iterations was chosen sufficiently large to assure stationarity of the dynamics. Transient dynamics were omitted from data analyses and graphical presentations of the results. The model is available at Zenodo website [49].

Phase-frequency coupling
For P 1 a value of 4.5, lower than the limit cycle threshold, was selected (Table 2), which assures that the node operates as a resonance filter circuit in the gamma band. As shown in Fig 2, for every positive phase of the slow oscillation in Node 1, the average zero crossings of fast oscillation in Node 1 was increased and for every negative phase it was decreased. Amplitude of the fast oscillation was slightly modulated as well, but far less pronounced than frequency. This illustrates that in the resonance regime, the fast oscillation is more sensitive to frequency than to amplitude modulation. The value of P 2 was set to 0 in order to deactivate the fast oscillation in the second node. The PFC arises within Node 1. The slow oscillation in the first node was synchronized with that in the second node (Fig 2). As a consequence, the frequency of the fast oscillation in Node 1 also correlates with the phase of the slow oscillations in Node 2, resulting in cross-node PFC. Note that cross-node PFC was generated only through the settings of mean noise input levels, P 1 and P 2 , while connection strengths remained fixed. The cross-node PFC arises through cross-node phase-phase coupling of the slow oscillations, combined with local phase-frequency coupling. Although the fast oscillation in Node 2 appears to be oscillating at around 18Hz, this is due to the cut-off frequency of the filters. As we will report in the Power

Phase-amplitude coupling
The mechanism for implementing cross-node PAC is similar to above, except that in the first node, the value of parameter P 1 was set to 7, which is higher than the limit cycle threshold which triggers a steady-state gamma band oscillation. P 2 again was set to 0, in order to deactivate the fast oscillation in the second node ( Table 2). The amplitude of the fast oscillation in Node 1 was modulated by the slow activity via the coupling between the slow and the fast inhibtory interneurons. Through the synchronization of the slow oscillations across the nodes, the slow oscillation in Node 2 also correlates with the amplitude changes of the fast oscillation in Node 1, resulting in cross-node PAC (Fig 3).

Frequency-frequency coupling
With FFC, two distinct fast oscillations change their frequencies over time. Thus, two frequency-modulated signals are required. To enable this, P 1 and P 2 were both set to 4.5, i.e. in the resonance regime. The resulting FFC is shown in Fig 4. Fast oscillations in both nodes show increase and decrease in average number of zero crossings within, respectively, the negative and positive phases of the slow oscillations. The synchronized slow oscillations modulate Amplitude-amplitude coupling AAC refers to the correlation between two amplitude signals from two distinct fast oscillations. Therefore, AAC requires a pair of amplitude-modulated signals. To achieve this, P 1 and P 2 were both set to 7, well above the limit cycle threshold. Time constants for the fast interneuron populations self-feedback were set to obtain a fast oscillation of about 51Hz in Node 1, and of about 42Hz in Node 2. As shown in Fig 5, the amplitudes of the fast oscillations were modulated by synchronization between the slow oscillations. As a result, the amplitude of the fast oscillation in Node 1 correlated with that of the fast oscillation in Node 2. AAC is thus the product of two local PACs of which the slow signals are synchronized between nodes.

Amplitude-frequency coupling
AFC is the most complex type of CFC; amplitude of one fast oscillation correlates with the fluctuations in the frequency of the other. To obtain AFC, PAC is generated in Node 1, i.e., the P 1 is set to 7, well above the limit-cycle threshold value, and P 2 is set to 4.5, in the resonance regime. PAC within Node 1 and PFC within Node 2 were combined via the cross-node phase synchronization. The synchronized slow oscillations between nodes simultaneously modulate the amplitude of the fast oscillation in Node 1 and the frequency of the fast oscillation in Node 2, resulting in AFC between the fast oscillations of Nodes 1 and 2 (Fig 6). Power Spectrum Densities (PSDs) of model output Fig 7 shows PSDs of outputs from Nodes 1 and 2 in each CFC simulation. Peaks appear in the delta (< 4Hz) and gamma (>30Hz) bands, which indicates stationary oscillatory activity. In case of frequency modulation (Panels a and c), the gamma activity shows wider peaks than in case of amplitude modulation (Panels b and d). In Panel e, where both modulations occur, a narrow (Node 1) and a wide (Node 2) bands can be observed in the gamma activity. In the PFC and PAC results (Figs 2 and 3), the fast oscillation in Node 2 appeared to have oscillatory activity around 18Hz. However, no peak at 18Hz exists in Panels a and b. Thus, the apparent '18Hz oscillation' is shown to be an artifact of filtering. (Note that numbers of zero crossings of fast oscillations in Figs 2-6 are slightly lower than that of their actual frequency due to the estimation method of the zero crossing)

Discussion
As suggested by Jensen and Colgin [2], coupling between slow and fast neural activity could take one of six forms: PFC, PAC, FFC, AAC, AFC or PPC. The first five of these were obtained in a model consisting of two mutually connected four-populations NMM nodes. For PAC or PPC, models have previously been developed with greater neurological detail than the present one [50][51][52]. Instead, we aimed to understand the systematic relationship among the cross-frequency couplings (CFCs). Neural mass models afford this based on their relative simplicity while preserving sufficient neural plausibility for simulating population level electrophysiological signals [30,33,45]. The neural mass modeling approach has been used for both quantitative estimation of the system parameters in cortical and subcortical regions [43,[53][54][55] and qualitative simulation for neural systems behavior [31,[56][57][58]. Here, all results are qualitative, as they serve the purpose of introducing a plausible common mechanism for the emergence of different types of CFC. Each node of the model represents a patch of cortical tissue, anywhere between a few millimeters and a centimeter in diameter, in which oscillations in various frequencies can be observed [41]. Whereas nodes were connected symmetrically to represent interconnected brain electrical activity sources within a cortical region, asymmetric inter-node connectivity is more likely to represent hierarchical structures, such as top-down/bottom-up connectivity in visual cortex [59]. Addressing this issue will require further elaboration of the inter-node connectivity of the model.
To generate different types of CFC, the model structure and connectivity were not changed, rather, the level of input noise of the nodes (parameters P 1 and P 2 ) was varied. The noise level determines the stability of the fast activity: when P is in the resonance regime, local PFC occurs, while the parameter is set above the limit-cycle threshold, local PAC emerges. The results showed that PFC and PAC are two phases of a CFC which are determined by the stability of the fast inhibitory activity. The modeled fast activity represents interneuron gamma (ING). Thus, the simulation results predict that, as the mean noise input level to a cortical region changes, stability of ING will change, and CFC will change in appearance between frequency and amplitude modulation. The slow activity of the nodes was synchronized via the inter-node connection. The local fast activities are modulated via the synchronized slow activity. Our results show how local ING can be modulated across brain regions via a common intrinsic slow activity, e.g., alpha band activity in visual cortices. Cross-regional PAC was reported even between bilateral motor cortices in monkeys [60]. For the inter-hemispheric PAC, other regions, such as, the subthalamic nucleus, could also contribute as a pacemaker. A pacemaker region could replace or complement the between-node pyramidal population [43,61]. This would not affect to the framework of cross-regional PAC/PFC generation, as long as the slow activities of the individual nodes remain synchronized in the same frequency band. Across brain regions of which the main frequencies differ, e.g., alpha rhythms in visual cortices and theta in hippocampus, m:n synchrony (where m6 ¼n) is needed for modulation of local fast activities. Neural mass modeling studies have been addressing the issue of cross-frequency phase-phase coupling (PPC), in which coupling strength between the nodes is the major parameter of interest [46]. The current study, however, does not address the issue since our main parameters were the mean noise input levels.
The model also successfully implemented AAC, FFC and AFC. In all simulations, amplitude and/or frequency of fast activity was modulated by synchronized slow activity. The assumption of a common modulator is in accordance with empirical studies. For example, AAC measurement provided evidence that two fast activities were modulated by a common slow activity in LFP signals [24][25][26]. To our knowledge, PFC, FFC and AFC have not yet been reported in brain signals. The model predicts FFC if slow activity in two regions is phase synchronized while stability of ING in both is intermediate, i.e., in the resonance regime. Likewise, AFC is expected when the slow activity is synchronized while ING is stable in one region but intermediately stable in the other.
Our PFC, FFC and AFC results provide a possible explanation why frequency modulations are scarcely reported. Frequency modulation occurred when the model was in the resonance regime, where amplitude modulation also occurrs to some extent. This suggests that in the brain signal, PFC and PAC may co-exist. Since PFC measures are still under-developed compared to PAC measures [6,[62][63][64][65], PFC may fail to be observed in a noisy signal such as EEG, while the weaker PAC is detected.
The frequency modulation could be more than a mere transient state to amplitude modulation. Empirical studies have shown that fast inhibitory interneurons fire at a specific phase of the gamma band oscillation [37,66]. The gamma frequency, in this context, is synonymous to the average number of spiking. As the gamma frequency increases/decreases according to the slow phases, the inhibitory population changes in excitability. In the current model, the range between the maximum and minimum frequencies is larger in the frequency modulation than in the amplitude modulation. Such a wide range of excitability is beneficial, for example, for implementing a gating function for information transfer.
In our model, frequency was determined by the fixed time constants of the neural populations (see also [67]). Frequency modulation was realized by through the effect of mean input noise level on the gain of the sigmoid block. The product of the gain and the time constant determines the resonance frequency. In principle, the time constants could also be considered as function of noise level [40]. This, however, would not change the principle behind the current approach. We conclude that two fundamental forms of CFCs in cortex are PFC and PAC. They emerge as modulations of ING. The couplings can either be kept local or spread across regions via PPC, depending on the synchronization of slow intrinsic oscillation between these regions.