Generation of Phase-amplitude Coupling of Neurophysiological Signals in a Neural Mass Model of a Cortical Column

Phase-amplitude coupling (PAC), the phenomenon where the phase of a low-frequency rhythm modulates the amplitude of a higher frequency, is becoming an important neurophysiological indicator of short- and long-range information transmission in the brain. Although recent evidence suggests that PAC might play a functional role during sensorimotor, and cognitive events, the neurobiological mechanisms underlying its generation remain imprecise. Thus, a realistic but simple enough computational model of the phenomenon is needed. Here we propose a neural mass model of a cortical column, comprising fourteen neuronal populations distributed across four layers (L2/3, L4, L5 and L6). While experimental studies often focus in only one or two PAC combinations (e.g., theta-gamma or alpha-gamma) our simulations show that the cortical column can generate almost all possible couplings of phases and amplitudes, which are influenced by connectivity parameters, time constants, and external inputs. Furthermore, our simulations suggest that the effective connectivity between neuronal populations can result in the emergence of PAC combinations with frequencies different from the natural frequencies of the oscillators involved. For instance, simulations of oscillators with natural frequencies in the theta, alpha and gamma bands, were able to produce significant PAC combinations involving delta and beta bands.


Introduction
It has been hypothesized that phase-amplitude coupling (PAC) of neurophysiological signals plays a role in the shaping of local neuronal oscillations and in the communication between cortical areas (Canolty and Knight, 2010).PAC occurs when the phase of a low frequency oscillation modulates the amplitude of a higher frequency oscillation.A classic example of such phenomenon was demonstrated in the CA1 region of the hippocampus (Bragin et al., 1995), where the phase of the theta band modulates the power of the gamma-band.Computational models of the theta-gamma PAC generation in the hippocampus have been proposed (Kopell et al., 2010) and are based on two main types of models.The first model consists in a network of inhibitory neurons (I-I model) (White et al., 2000) whereas the second model is based on the reciprocal connections between networks of excitatory pyramidal and inhibitory neurons (E-I model) (Kopell et al., 2010;Tort et al., 2007).In such models fast excitation and the delayed feedback inhibition alternate, and with appropriate strength of excitation and inhibition, oscillatory behavior may continue for a while.When the gamma activity produced by the E-I or I-I models is periodically modulated by a theta rhythm imposed by either an external source or theta resonant cells within the network (White et al., 2000), a thetagamma PAC is produced.Recently, the generation of theta-gamma PAC was studied (Onslow et al., 2014) using a neural mass model (NMM) proposed by Wilson and Cowan (Wilson and Cowan, 1972).In NMMs spatially averaged magnitudes are assumed to characterize the collective behavior of populations of neurons of a given type instead of modeling single cells and their interactions in a realistic network (Jansen and Rit, 1995;Wilson and Cowan, 1972).Specifically, the Wilson and Cowan model consists of an excitatory and inhibitory populations mutually connected.While the models mentioned above have improved our understanding of the physiological mechanism that give rise to theta-gamma PAC, we lack modeling insights into the generation of PAC involving other frequency pairs.This is critical since experimental studies have shown that the PAC phenomenon is restricted neither to the hippocampus nor to theta-gamma interactions.In fact, PAC has been detected in pairs involving all possible combinations of low and high frequencies: deltatheta (Lakatos et al., 2005), delta-alpha (Cohen et al., 2009;Ito et al., 2013), delta-beta (Cohen et al., 2009;Nakatani et al., 2014), delta-gamma (Florin and Baillet, 2015;Gross et al., 2013;Lee and Jeong, 2013;Nakatani et al., 2014;Szczepanski et al., 2014), theta-alpha (Cohen et al., 2009), thetabeta (Cohen et al., 2009;Nakatani et al., 2014), theta-gamma (Demiralp et al., 2007;Durschmid et al., 2013;Florin and Baillet, 2015;Lakatos et al., 2005;Lee and Jeong, 2013;McGinn and Valiante, 2014;Wang et al., 2011), alpha-beta (Sotero et al., 2013), alpha-gamma (Osipova et al., 2008;Spaak et al., 2012;Voytek et al., 2010;Wang et al., 2012), andbeta-gamma (de Hemptinne et al., 2013;Wang et al., 2012).Furthermore, although experimental studies usually focus in one or two PAC combinations, most of the combinations mentioned above can be detected in a single experiment (Sotero et al., 2013).This suggest a diversity and complexity of the PAC phenomenon that haven't been grasped by current computational models.
In this work we propose a neural mass model of a cortical column that comprises 4 cortical layer and 14 neuronal populations, and study the simultaneous generation of all PAC combinations mentioned above.The neuronal populations modeled have natural frequencies in the theta, alpha and gamma bands.However, due to the effective connectivity between them, oscillations at the delta and beta bands appear and result in PAC involving these frequencies.We then focus on five combinations: delta-gamma, theta-gamma, alpha-gamma, delta-beta, and beta-gamma, and explore how changes in model parameters such as strength of the connections, time constants and external inputs, strengthen or weaken the PAC phenomenon.

A neural mass model of a cortical column
Figure 1 shows the proposed model obtained by distributing four cell classes in four cortical layers (L2/3, L4, L5, and L6).This produced 14 different neuronal populations, since not all cell classes are present in every layer (Neymotin et al., 2011).Excitatory neurons were either regular spiking (RS) or intrinsically bursting (IB), and inhibitory neurons were either fast-spiking (FS), or low-threshold spiking (LTS) neurons.Each population performs two operations.Post-synaptic potentials (PSP) are converted into an average density of action potentials (AP) using the sigmoid function: where the variable x represents the PSP, and parameters e 0 , v 0 and r, stand for maximal firing rate, the PSP corresponding to the maximal firing rate e 0 , and the steepness of the sigmoid function, respectively.The second operation is the conversion of AP into PSP, which is done by means of a linear convolution with an impulse response ݃ሺ‫ݐ‬ሻ given by: ݃ሺ‫ݐ‬ሻ ൌ ‫݁ݐ݇ܩ‬ ି௧  (2) where G controls the maximum amplitude of PSP and k is the sum of the reciprocal of membrane average time constant and dendritic tree average time delay (Jansen and Rit, 1995).The convolution model with impulse response (2) can be transformed into a second order differential equation (Jansen and Rit, 1995;Sotero et al., 2007).Then, the temporal dynamics of the average PSP in each neuronal population ‫ݔ‬ can be obtained by solving a system of 14 second order differential equations: where n = 1,…,14 and m = 1,…,14.The populations are numbered from 1 to 14 following the order: [L2RS, L2IB, L2LTS, L2FS, L4RS, L4LTS, L4FS, L5RS, L5IB, L5LTS, L5FS, L6RS, L6LTS, L6FS].Notice we labeled layer 2/3 simply as L2.As can be seen in (3), neuronal populations interact via the connectivity matrix Γ .Inputs from neighboring columns are accounted via ‫‬ which can be any arbitrary function including white noise (Jansen and Rit, 1995).Thus, (3) represents a system of 14 stochastic differential equations.The 'damping' parameter ܾ critically determines the behavior of the system.If we set to zero the connections between the populations ሺΓ ൌ 0, ݊ ് ݉ሻ then for ܾ 1 (overdamped oscillator) and ܾ ൌ 1 (critically damped oscillator) each neuronal population returns to steady-state without oscillating.If ܾ ൏ 1 (underdamped oscillator) each population is capable of producing oscillations even if the inter-population coupling is set to zero.
Note that the case ܾ ൌ 1 corresponds to the Jansen and Rit model (Jansen and Rit, 1995) which has been extensively used in the literature (David and Friston, 2003;Grimbert and Faugeras, 2006;Sotero and Trujillo-Barreto, 2008;Sotero et al., 2007;Ursino et al., 2010;Valdes-Sosa et al., 2009;Zavaglia et al., 2006).Thus, in those models an individual population is not capable of oscillating, and is the presence of inter-population connections (nonzero Γ , ݊ ് ݉) that produces oscillatory behavior that mimics observed Electroencephalography (EEG) signals.However, realistic models of only one population are able to produce oscillations (Wang and Buzsaki, 1996).To account for this possibility we introduced the parameter ܾ with values ܾ ൏ 1. Tables 1 and 2 presents the parameters of the model and their interpretation. .

Computation of phase-amplitude coupling
The analytic representation yሺ‫ݐ‬ሻ of a filtered signal xሺ‫ݐ‬ሻ can be obtained by means of the Hilbert transform ‫ܪ‬൫xሺ‫ݐ‬ሻ൯: where aሺ‫ݐ‬ሻ and φሺ‫ݐ‬ሻ are the instantaneous amplitude and phase, respectively.As a measure of PAC we will use the flow of information from the phase to the amplitude.Recently it was shown (Liang, 2014) that the information flow from a time series to another one (in our case from φሺ‫ݐ‬ሻ to aሺ‫ݐ‬ሻ) can be calculated as: . This is a provisional file, not the final typeset article where ‫ܥ‬ ఝ is the sample covariance between aሺ‫ݐ‬ሻ and φሺ‫ݐ‬ሻ, ‫ܥ‬ ,ௗ is the covariance between aሺ‫ݐ‬ሻ and ݀aሺ‫ݐ‬ሻ/݀‫,ݐ‬ ‫ܥ‬ ఝ,ௗ is the covariance between φሺ‫ݐ‬ሻ and ݀aሺ‫ݐ‬ሻ/݀‫,ݐ‬ and ‫ܥ‬ and ‫ܥ‬ ఝఝ are the variances of the instantaneous amplitudes and phases, respectively.The units of ܶ ఝ՜ given by equation ( 5) are in nats/s.A nonzero value of the information flow ܶ ఝ՜ indicates PAC.A positive ܶ ఝ՜ means that φሺ‫ݐ‬ሻ functions to make aሺ‫ݐ‬ሻ more uncertain, whereas a negative ܶ ఝ՜ indicates that φሺ‫ݐ‬ሻ tends to stabilize aሺ‫ݐ‬ሻ (Liang, 2014).
A significance value can be attached to ܶ ఝ՜ using a surrogate data approach (Penny et al., 2009), where we shuffle the amplitude time series aሺ‫ݐ‬ሻ and calculate (using ( 5)), 1000 surrogate values.From this surrogate dataset we first compute the mean μ and standard deviation ߪ, and then compute a zscore as: Values satisfying |ܼ| 1.96, are significant with ߙ ൌ 0.05.Significant Z values are set to 1 and nonsignificant ones to zero.Then, Z values are multiplied to ܶ ఝ՜ .

Results
Simulated data were generated by the numerical integration of the system (3).For this, the local linearization method (Jimenez and Biscay, 2002;Jimenez et al., 1999) was used with an integration step of 10 -4 s.The values of the parameters are shown in tables 1 and 2. Five seconds of data were simulated and the first two seconds were discarded to avoid transient behavior.Thus, subsequent steps were carried out with the remaining three seconds.Figure 3 shows the spectrum of the signals presented in Figure 2. The six excitatory populations have their main spectrum peak in the alpha band, but they also present energy in the delta and theta band.
Slow inhibitory populations have the highest peak in the theta band, but also have energy in delta, alpha and beta bands.Fast inhibitory populations were set to present a peak in the gamma band but due to the interaction with other populations, they present significant peaks in other frequencies as well, especially in the theta and alpha bands.This is evident from Figure 4  To test the existence of PAC, we filtered each PSP in Figure 2 into six frequency bands: delta (0.1-4 Hz), theta (4-8 Hz), alpha (8-12 Hz), beta (12-30 Hz), and gamma (30-80 Hz).To this end, we designed finite impulse response (FIR) filters using Matlab's signal processing toolbox function firls.m.To remove any phase distortion, the filters were applied to the original time series in the forward and then the reverse direction using Matlab's function filtfilt.m.The Hilbert transform was then applied and instantaneous phases and amplitudes for each frequency band and each neuronal population were obtained.Ten different PAC combinations between a low-frequency phase and a higher-frequency amplitude were computed: delta-theta, delta-alpha, delta-beta, delta-gamma, thetaalpha, theta-beta, theta-gamma, alpha-beta, alpha-gamma, and beta-gamma.Each PAC combination consisted of a matrix of 14x14 PAC values representing all possible interactions between the 14 neuronal populations.
Results shown in Figure 5 include nine out of the ten PAC combinations.The theta-alpha PAC combination was not included since no significant values were obtained for the set of parameters used.The strongest PAC values were found for the delta-theta, alpha-beta, and delta-beta combinations.Specifically, negative information flow from the phase of alpha in L5RS, L5LTS and L5FS to the amplitude of beta in L4FS.This means that the phase of alpha in those three populations stabilizes the amplitude of beta in L4FS.On the other hand, the strong positive information flow from the phase of delta in L4FS to the amplitude of theta in L5RS and L6RS means that L4FS tends to make the activity in L5RS and L6RS more uncertain.As seen in the figure, theta-gamma, alphabeta, and alpha-gamma PAC, presented only negative information flow from phases to amplitudes.
Delta-gamma presented only positive flow and it was from L6RS to L2FS.Delta-theta, delta-alpha, delta-beta, and theta-beta, presented both positive and negative information flow.

Generation of theta-gamma and alpha-gamma PAC
Some of the parameters presented in Table 1 were taken from the neural mass literature (Jansen and Rit, 1995;Wendling et al., 2000), and the ones with no equivalent in the literature were assigned physiologically reasonable values.Thus, it is necessary to explore how the change in the parameters affect the PAC values.In this section, for the sake of simplicity, we focus on two PAC combinations which involve the gamma rhythm and have been of great interest in the literature: theta-gamma, and alpha gamma.
Figure 6 shows results for the theta-gamma PAC combination when varying external input and time constant parameters values.Each panel present a different simulation, in which we varied one of the parameters while maintaining the others with the same values as in Table 1.Panels from A) to D) present PAC between the 14 populations for different values of the external input ‫‬ ହ to the L4RS population: 0, 150, 300, and 1000, respectively.Panels from E) to H) present PAC values when the external input ‫‬ to the L4FS population was changed to: 0, 300, 500, and 800, respectively.Finally, panels from I to L present PAC values when time constants for the L2/3RS and L4RS, ݇ ଵ , ݇ ହ , populations where changed to: ݇ ଵ ൌ ݇ ହ ൌ 40, 100, 200, 300 ‫ݏ‬ ିଵ , respectively.
Figure 7 shows results for the theta-gamma PAC combination when connectivity parameters are varied.Panels from A) to C) show results when multiplying the entire connectivity matrix C with a factor of 1.5, 0.5 and 10, respectively.Panels from D) to I) present results when changing only one connection in the connectivity matrix: D) ‫ܥ‬ ଵଵ ൌ 40, E) ‫ܥ‬ ଷଷ ൌ െ40 , F) ‫ܥ‬ ଵ ൌ 0, G) ‫ܥ‬ ଵସଵସ ൌ െ45 , H) ‫ܥ‬ ଵଶହ ൌ 0, I) ‫ܥ‬ ହ ൌ 0. Figures 8 and 9 present the same simulations as in Figures 6 and 7, but for the alpha-gamma PAC.
The strongest positive PAC of all the simulations explored was found in the theta-gamma combinations for: ݇ ଵ ൌ ݇ ହ ൌ 200 ‫ݏ‬ ିଵ (Figure 6K) and it corresponded to the connection from L2FS to L2LTS.The strongest negative PAC was also found in Figure 6K and corresponded to the connection from L4RS to L4LTS.Thus, layers 2/3 and 4 received the strongest information flow from the phases of all populations in the cortical column.Overall, changes in external input and time constants produced theta-gamma PAC higher than alpha-gamma PAC, whereas changes in connectivity produced higher alpha-gamma PAC values.

Generation of delta-beta, delta-gamma, and beta-gamma PAC
As shown in the figures 3 and 4, although the populations in the cortical columns were modeled as oscillators with natural frequencies in the theta, alpha and gamma band, due to connections between them, activity in the delta and beta band emerged.This produces PAC involving delta phases and beta phases and amplitudes as shown in Figure 5.In this section we focus on three PAC combinations: delta-beta, delta-gamma, and beta-gamma, and explore how their strength change as The strongest PAC value was positive and was found between the phase of delta in L4RS and the amplitude of beta in L2LTS (Figure 11G).This was the strongest PAC across all the 10 PAC combinations computed in this paper for the range of parameters explored.Another interesting results was that for the beta-gamma combination, when varying the connectivity parameters all values were found to be negative (Figure 15), meaning that beta oscillations tended to stabilize the gamma oscillations.

Discussion
We have proposed a neural mass model that captures the phase-amplitude coupling between layers in a cortical column.The model comprises fourteen interconnected neuronal populations distributed across four cortical layers (L2/3, L4, L5 and L6).Excitation and inhibition that impinge on a population are represented by the average PSP that are elicited in dendrites belonging to a population.
Phase-amplitude coupling between the neuronal populations was modeled using a measure of the information flow between two times series, proposed recently by (Liang, 2014).

Simplifications and assumptions
We omitted layer 1, because it does not include somas (Binzegger et al., 2004).Based on experimental reports on the strength of the inputs to each layer (Binzegger et al., 2004;Jellema et al., 2004), we considered external inputs to the RS and FS populations in layer 4, thus neglecting possible external inputs to other layers.
In this work, we used the information flow (or information transfer) to measure the causation between the phase of a low frequency rhythm and the amplitude of a higher frequency rhythm, and used it as an index of PAC.We chose this measure over transfer entropy (Schreiber, 2000) because the latter is difficult to evaluate, requiring long time series (Hlavackova-Schindler et al., 2007).
Moreover, recent studies have shown that it is biased as its values depends on the autodependency coefficient in a dynamical system (Runge et al., 2012).On the other hand, the information flow (Liang, 2014) is calculated using a simple analytical expression, equation ( 5), that depends only on sample covariances, which is obtained under the assumption that the two time series are components of a general two dimensional linear system.Then, our use of equation ( 5) to measure PAC relies on the assumption that the relationship between phases and amplitudes of different frequencies is well described by a linear model.Although this seems restrictive at first glance, we should remember that a linear model (and in general any model) of cross-frequency coupling (CFC) is necessarily the result of a nonlinear model at the level of the neuronal activities (Chen et al., 2008).

On the generation of PAC combinations involving delta and beta rhythms
Our simulations suggest that the interconnection of neuronal populations can produce PAC combinations with frequencies different from the natural frequencies of the oscillators involved.For instance, our model of oscillators with natural frequencies in the theta, alpha and gamma bands, was able to produce significant PAC involving delta and beta rhythms: delta-theta, delta-alpha, delta-beta, delta-gamma, theta-beta, and alpha-beta.The strength of the coupling was affected by the strength of the connections between the populations, the inhibitory and excitatory time constants and the strength of the external input to the model.Interestingly, the strongest PAC value found over the range of the parameters explored in this paper, was not in the theta-gamma or alpha-gamma combination, but in .CC-BY 4.0 International license peer-reviewed) is the author/funder.It is made available under a The copyright holder for this preprint (which was not .http://dx.doi.org/10.1101/023291doi: bioRxiv preprint first posted online Jul. 27, 2015; Phase-amplitude coupling in a cortical column model the delta-beta combination when varying the connectivity parameters.This demonstrates that effective connectivity values are critical in the emergence of the PAC phenomenon.

Realistic vs neural mass models of PAC generation
The first computational models of PAC generation were realistic models of the theta-gamma coupling in the hippocampus (Kopell et al., 2010).These models considered networks of hundreds of interconnected neurons which were individually modeled by either a single compartment (White et al., 2000) or realistically represented by multiple compartments for the soma, axon, and dendrites (Tort et al., 2007).A practical disadvantage of this approach is that it needs high computational power.But more important than that, the use of such realistic models produces hundreds or thousands of variables and parameters, making it difficult to determine their influence on the generated average network characteristics.This is especially critical if we are interested in analyzing the link of PAC and other mesoscopic phenomenon like functional magnetic resonance signals (Wang et al., 2012).The analysis of multiple PAC combinations as done in this paper would be even more difficult with realistic metworks.By comparison, our model of one cortical column comprised only 14 second-order (or 28 first-order) differential equations, which can be easily solved in any modern personal computer.

The generated dynamics is not restricted to PAC
The neural mass model presented in this paper is relatively simple but it can generate a rich temporal dynamics.Studies of the dynamics generated by the Jansen and Rit model, which is the basis for our model, can be found elsewhere (Faugeras et al., 2009;Grimbert and Faugeras, 2006;Spiegler et al., 2010).In this paper we focused on PAC, which is only one type of the more general phenomenon of CFC which is the result of the nonlinearities in the brain dynamics.Then, it is not unexpected to find other types of CFC in the signals generated by our model (see for instance the temporal dynamics of L6RS in Figure 2 which corresponds to frequency modulation).Thus, in addition to PAC, other five types of CFC could be explored (Jirsa and Muller, 2013): amplitude-amplitude coupling (AAC), phase-phase coupling (PPC), frequency-frequency coupling (FFC), phase-frequency coupling (PFC), and frequency amplitude coupling (FAC).All these types of couplings can be easily calculated using equation (5) after replacing aሺ‫ݐ‬ሻ and φሺ‫ݐ‬ሻ with the appropriate time series.This allows to compute all types of CFC with the same equation, and easily compare them since all results will be in the same unit of measurement (i.e.nats/s).Zavaglia, M., Astolfi, L., Babiloni, F., and Ursino, M. (2006).A neural mass model for the simulation of cortical activity estimated from high resolution EEG during cognitive or motor tasks.J Neurosci Methods 157, 317-329.

Figure 2
Figure 2 presents the temporal evolution of the average PSP in each neuronal population.Time series colored in red correspond to excitatory populations (L2RS, L2IB, L4RS, L5RS, L5IB, L6RS) whereas inhibitory populations (L2LTS, L4LTS, L5LTS, L6LTS) are presented in green.As seen in the figure, the generated signals present the characteristic 'waxing and waning' (i.e, amplitude modulation) observed in real EEG signals.
where we show the spectrum of the population when all values in the connectivity matrix are set to zero: Γ ൌ 0. Peaks in Figure4correspond to the natural frequency of oscillation of the populations: L2RS (9

.
figures 14 and 15 shows the results for the beta-gamma PAC.

Figure 1 .
Figure 1.Proposed neural mass model of the cortical column.A) Layer distribution of the four neuronal types.The excitatory populations are the intrinsically bursting (IB), and the regulatory spiking (RS).The inhibitory population are the low-threshold spiking (LTS) and fast spiking (FS).B) Connectivity matrix values used for coupling the 14 populations modeled.Negative values correspond to inhibitory connections.

Figure 2 .
Figure 2. Simulated temporal evolution of the postsynaptic potentials of all populations for one realization of the model.Excitatory populations are depicted in red and inhibitory ones in green.

Figure 3 .
Figure 3. Frequency spectrum of the postsynaptic potentials shown in Figure 2. Excitatory populations are depicted in red and inhibitory ones in green.

Figure 4 .
Figure 4. Frequency spectrum of the neuronal populations when the connectivity matrix is set to zero.Excitatory populations are depicted in red and inhibitory ones in green.

Figure 5 .
Figure 5. Phase-amplitude coupling corresponding to the simulation presented in Figure 1.Nosignificant values were set to zero and are depicted in white.A) delta-theta B) delta-alpha C) deltabeta D) delta-gamma E) theta-alpha F) theta-beta G) theta-gamma H) alpha-beta I) alpha-gamma.

Figure 7 .
Figure 7. Theta-gamma phase-amplitude coupling: changes in the connectivity parameters.A)