A quantitative theory of gamma synchronization in macaque V1
Abstract
Gammaband synchronization coordinates brief periods of excitability in oscillating neuronal populations to optimize information transmission during sensation and cognition. Commonly, a stable, shared frequency over time is considered a condition for functional neural synchronization. Here, we demonstrate the opposite: instantaneous frequency modulations are critical to regulate phase relations and synchronization. In monkey visual area V1, nearby local populations driven by different visual stimulation showed different gamma frequencies. When similar enough, these frequencies continually attracted and repulsed each other, which enabled preferred phase relations to be maintained in periods of minimized frequency difference. Crucially, the precise dynamics of frequencies and phases across a wide range of stimulus conditions was predicted from a physics theory that describes how weakly coupled oscillators influence each other’s phase relations. Hence, the fundamental mathematical principle of synchronization through instantaneous frequency modulations applies to gamma in V1 and is likely generalizable to other brain regions and rhythms.
https://doi.org/10.7554/eLife.26642.001Introduction
Synchronization, the ability of oscillators to mutually adapt their rhythms (Pikovsky et al., 2002; Winfree, 1967), is a ubiquitous natural phenomenon. Neural synchronization in the gammarange has been reported both in subcortical structures (Akam et al., 2012; Steriade et al., 1993; Zhou et al., 2016) and in cortical areas (Fries, 2015; Gray and Singer, 1989; Gregoriou et al., 2009). Gamma rhythms emerge in activated neural circuits in which fastspiking inhibitory neurons play a central role (Cardin et al., 2009; Tiesinga and Sejnowski, 2009; Traub et al., 1996). A prime example is the emergence of gamma rhythms in the early visual cortex during visual stimulus processing (e.g. Brunet et al., 2015; Gail et al., 2000; Gray and Singer, 1989; Hermes et al., 2015; Ray and Maunsell, 2010; Roberts et al., 2013). Gamma synchronization has been related to the formation of neural assemblies within (Gail et al., 2000; Gray and Singer, 1989; Havenith et al., 2011; Vinck et al., 2010) and across brain areas (Bosman et al., 2012; Gregoriou et al., 2009; Grothe et al., 2012; Jia et al., 2013a; Roberts et al., 2013; Sirota et al., 2008; Zhou et al., 2016).
The precise temporal coordination of presynaptic spikes increases their effectiveness on postsynaptic targets (Fries et al., 2001; Tiesinga et al., 2004) and can thereby modulate the effectiveness of neural communication (Börgers et al., 2005; Cannon et al., 2014; Womelsdorf et al., 2007), as shown between V1 and V4 during visual attention (Bosman et al., 2012; Grothe et al., 2012). Temporal coordination in terms of spike timing (phase code) might be an efficient and robust mechanism for information coding (Havenith et al., 2011; Jensen et al., 2014; Maris et al., 2016; Tiesinga and Sejnowski, 2009; Vinck et al., 2010). Further, gamma rhythmic inhibition might increase coding efficiency through sparsening (Chalk et al., 2015; Jadi and Sejnowski, 2014; Vinck and Bosman, 2016) and normalization (Gieselmann and Thiele, 2008; Ray et al., 2013) of neural activity. These network consequences of gamma have led to influential hypotheses about the function of gamma in sensation and cognition (Buehlmann and Deco, 2010; Buzsáki and Wang, 2012; Eckhorn et al., 2001; Fries, 2015; Gray and Singer, 1989; Maris et al., 2016; Miller and Buschman, 2013), including a role in perceptual grouping (Eckhorn et al., 2001; Engel et al., 1999; Gray and Singer, 1989) and in visual attention (Bosman et al., 2012; Fries, 2015; Gregoriou et al., 2009; Miller and Buschman, 2013).
Surprisingly, in spite of important scientific advances, it is not well understood how gamma rhythms synchronize and what the underlying principles of synchronization are. For example, recent experimental observations of large variability in gamma oscillation frequency have raised doubts about the robustness and functionality of gamma synchronization in the brain. It has been observed that frequency fluctuates strongly over time (Atallah and Scanziani, 2009; Burns et al., 2010; 2011) and that different cortical locations can express different preferred frequencies at a single moment in time (Bosman et al., 2012; Ray and Maunsell, 2010). That these observations have led to doubts on the functionality of gamma synchronization reveals a stationary view of synchronization, which assumes that the underlying oscillatory dynamics are stable at a fixed phaserelation and shared frequency. This is also reflected in the widespread use of stationary methods to assess gamma synchronization, of which spectral coherence is a prime example (Carter et al., 1973). From a dynamic systems perspective, however, synchronization is primarily a nonstationary process (Izhikevich, 2007; Izhikevich and Kuramoto, 2006; Kopell and Ermentrout, 2002; Pikovsky et al., 2002; Winfree, 1967), because oscillators mutually adjust their rhythms through phase shifts (i.e. through changes in the instantaneous frequency).
Here, by using a combination of theoretical and experimental techniques, we studied the dynamical principles of gamma synchronization in monkey visual area V1. We simultaneously recorded gammarhythmic neural activity at different V1 cortical locations and studied their synchronization properties while using local stimulus contrast (Ray and Maunsell, 2010; Roberts et al., 2013) to modulate the frequency difference averaged over time (detuning). Strikingly, even when the mean frequencies did not match (detuning > 0), we often observed that gamma rhythms synchronized. This was achieved by continuously varying their instantaneous frequency difference, which permitted the temporary maintenance of a preferred phase relationship during reoccurring periods of minimized instantaneous frequency difference. The interplay between the detuning and the amount of instantaneous frequency modulations regulated the phaselocking strength and the preferred phaserelation between V1 locations. Furthermore, to achieve a principled understanding of our observations, we applied the theoretical framework of weakly coupled oscillators to our data (Ermentrout and Kleinfeld, 2001; Hoppensteadt and Izhikevich, 1998; Kopell and Ermentrout, 2002; Kuramoto, 1991; Pikovsky et al., 2002). We found that a single differential equation accounted well for the nonstationary frequency modulations and further allowed for precise predictions of how the phaselocking and the phaserelation between gamma rhythms changed across conditions.
Results
Frequency differences regulate the phase synchronization process between local monkey V1 gamma rhythms
We first asked how synchronization within V1 was influenced by mean frequency differences, and by the distance between recording sites. To this aim, we recorded simultaneously from two to three laminar probes (each with 16 recording contacts spaced along the recording shaft, see Supplementary Materials for alignment procedure) in cortical area V1 of two macaques (M1 and M2) (Figure 1A). We used distances between probes of 1–6 mm, matching approximately the extent of V1 horizontal connectivity (Stettler et al., 2002). Notably, horizontal connectivity strength declines strongly with distance between cortical locations (Stettler et al., 2002), so that increasing interprobe distance indexes decreasing horizontal connectivity strength. The monkeys fixated centrally while a fullscreen static squarewave grating with spatially varying contrast was shown (Figure 1B). The local contrast varied periodically over visual space such that different contrasts were presented to different cortical locations. The magnitude of contrast difference (ranging from 0% to ~43%, see Table S1) was manipulated by varying the sign and amplitude of the spatial variation in contrast. The stimulus gratings induced gamma power in layers 2–4 and in the deepest layer (Figure 1C, Figure 1—figure supplement 1; van Kerkoerle et al., 2014; Xing et al., 2012). The gamma frequency increased systematically with higher local contrast (linear regression, single contact level, M1: R^{2} = 0.38, M2: R^{2} = 0.27, both p<10^{−10}, Figure 1D, Figure 1—figure supplement 2). The range of the frequency shift in our data (~5 Hz) was smaller than in Roberts et al., 2013, reflecting a narrower contrast range used here; from ~20% to ~60% (Table 1). The tight relationship between contrast and gamma frequency allowed us to induce different mean frequencies in nearby cortical locations separated by as little as 1–6 mm (e.g., Figure 2A).
The close positions of recording sites may have led to a contribution of volume condition to synchronization measures. The LFP, despite being local in comparison to extracranial electrical field measure like EEG, still might integrate signals over a scale of up to 1 cm horizontally (Kajikawa and Schroeder, 2011; Lindén et al., 2011; Xing et al., 2009), which may affect the interpretability of layerdependent analysis (Kajikawa and Schroeder, 2015). Using laminar probes enabled us to reduce the influence of volume conduction by calculating currentsource density (CSD), as the second spatial derivative of LFP signals measured along each probe (Mitzdorf and Singer, 1977; Schroeder et al., 1991a; Vaknin et al., 1988). The success in reducing volume conduction using CSD favors its use over LFP for spectral analysis at high spatial resolution. Next, we used a singular spectrum decomposition technique (SSD, (Bonizzi et al., 2014)) to extract gamma components from the CSD. From these singletrial gamma signals, we estimated the instantaneous frequency and phase at individual recording sites and the instantaneous phase difference between sites. In the example shown of a single pair of recording sites (Figure 2A), the stimulus induced a gamma frequency of 36 Hz at one probe and 32 Hz at the other, because different contrasts appeared in the respective receptive fields. As shown in the raw trace of the instantaneous phase differences in Figure 2B, the gamma phase difference was not constant over time, but continuously exhibited modulations and shifts. Sometimes, the phase difference changed slowly and at other times, it changed faster. The change of phase difference over time is called ‘phase precession’ (Pikovsky et al., 2002). Note that this should be distinguished from a phenomenon of the same name: the precession of preferred spiking phase in the theta cycle observed in rats moving through hippocampal place fields (Skaggs et al., 1996). In the present study, the ‘rate’ or ‘speed’ of precession is expressed as the instantaneous frequency difference in Hz (Figure 2C–D). We found that the observed modulations in phase difference were not random as would be expected if different frequencies precluded synchronization. Instead, the instantaneous phase difference was related to the instantaneous frequency difference. In Figure 2E, we plotted the instantaneous frequency difference as a function of the instantaneous phase difference. Each black point represents one momentary observation and the blue line the average using binning of 0.25 rad width. The plot shows that the instantaneous frequency difference (ΔIF) tends to be lower at certain phase differences than at others. We observed that the average frequency difference was close to 0 for phases differences between 0 and 2 radians, but was much higher at other phase relationships. The key to understanding how this dynamic relationship leads to synchronization is that phase relationships associated with lower frequency differences are maintained longer over time (slower precession) than phase relationships associated with higher frequency differences. This can be readily appreciated by the higher density of dots in Figure 2E between 0 and 2 radians.
Three key examples from our results nicely illustrate the dynamics of the relationship between instantaneous frequency difference and phase difference (Figure 3). These examples were derived from our experimental design, in which we varied the cortical distance between probes in a pair (varying horizontal connectivity strength), and in which we systematically varied for each pair the contrast difference (9 levels), and hence the mean gamma frequency difference. We show in these examples positive frequency differences for illustration, but negative differences were also present for single contact pairs in our data, depending on the sign of the contrast difference (see Figure 3—figure supplement 1). In the first example (Figure 3 column 1), we show two cortical locations separated by a relatively large distance of ~5 mm, presented with a visual contrast difference of 17% (Figure 3A). This yielded an overall mean frequency difference of 5 Hz (Figure 3B). If this frequency difference were constant, the phase difference would advance at a phase precession rate of 2π every 200 ms, which would preclude synchronization. However, the frequency difference was not constant. Instead, the instantaneous frequency difference changed as a function of phase difference (Figure 3B, Figure 3—figure supplement 1) with a modulation amplitude of ~1 Hz (approximately (maxmin)/2; see Appendix). At the smallest frequency difference (4 Hz, yellow point), the phase precession was slowest (2π every 250 ms). As a result, the probability distribution of phase differences over time (Figure 3C) was nonuniform giving a phaselocking value (Lachaux et al., 1999) (PLV) of 0.11. The peak of the distribution, the ‘preferred phase’, was at 1.3 rad, in line with the minimum of the instantaneous frequency modulation shape. In the second example, we chose a pair with a similar frequency difference of 4.8 Hz but a reduced distance (~2.5 mm, Figure 3D). The instantaneous frequency modulation was larger with a modulation amplitude of 1.8 Hz (Figure 3E) and a minimum around 3 Hz at the preferred phase. Because a lower minimum frequency difference corresponds to slower phase precession at the preferred phase than in the previous example, the preferred phase was maintained for longer. This resulted in a narrower phase difference distribution, indicating higher synchrony (PLV = 0.32, Figure 3F). The peak of the distribution was centered at a smaller phase difference (0.78 rad). In the third example, the cortical distance remained the same as in Figure 3D but the frequency difference was reduced (2.8 Hz) by eliminating the contrast difference (Figure 3G and Figure 1—figure supplement 2). Compared to example 2, the magnitude of the instantaneous frequency modulation did not change (modulation amplitude 1.8 Hz, Figure 3H), but showed a lower mean difference and a minimum close to zero (1 Hz, Figure 3H). Thus, the associated phase difference (0.48 rad) could be maintained for even longer periods and the phase difference probability distribution became even more pronounced and narrower (PLV = 0.51, Figure 3I). The three examples illustrate how the mean frequency difference and cortical distance (a proxy of the strength of horizontal interactions) determine the dynamic relationship between the instantaneous frequency difference and the phase difference during synchronization. In the following sections, we will show how these observations were characteristic of the whole dataset comprising 805 recorded acrossprobe contact pairs in monkey M1 and 882 pairs in monkey M2.
Before framing the relationship between instantaneous frequency modulations and the phase difference distribution in a mathematical manner (next section), we illustrate that relationship by an analogy with two cyclists on a circular circuit. Their speed is calculated as the number of circuits they complete in a given time, hence speed is analogous to frequency. Phase is equivalent to position on the track and phase difference is equivalent to the distance between the cyclists. Phase locking therefore is analogous to the amount of time they spend at a consistent distance from each other. The phase precession rate is analogous to the speed with which the distance between the cyclists changes. If the cyclists maintain constant speeds as they go around the track, the distance between them will vary at a constant rate, and they will only maintain a consistent distance (phase difference) if they both cycle at the same speed. Hence with stable instantaneous speed (i.e. frequency), phase locking is either absent or complete. If the cyclists do vary their speed, more complex patterns become possible. If they vary their speeds independently, the phasedifference distribution will be flat. However, suppose the slower cyclist can travel faster in the slipstream of the faster cyclist, then the amount of time the cyclists travel close to each other will be greater than the time they spend far away. They might cycle around the whole circuit this way (complete phase locking), but more likely the faster cyclist will get away and the cyclists will travel at their natural speed until they come together again allowing the slower cyclist to speed up. Hence, all phase differences are represented, but some are overrepresented – specifically phases where the faster cyclist is just in the lead. In our experiment, we measured the instantaneous gamma frequency, similar to looking on the speedometer of each bike, and the instantaneous phase difference, corresponding to the distance between the cyclists. This allowed us to understand the resultant probability distribution of phase differences (yielding the phase locking value and the average phase difference). We found that the average speed difference and the speed modulation strength defined the probability distribution.
The theory of weakly coupled oscillators (TWCO): A framework for cortical gamma synchronization
We now show how the observed synchronization behavior can be accounted for within the mathematical framework of the theory of weakly coupled oscillators (Ermentrout and Kleinfeld, 2001; Hoppensteadt and Izhikevich, 1998; Kopell and Ermentrout, 2002; Kuramoto, 1991; Pikovsky et al., 2002; Winfree, 1967). Many oscillatory phenomena in the natural world represent dynamic systems with a limitcycle attractor (Winfree, 2001). Although the underlying system might be complex (e.g. a neuron or neural population), the dynamics of the system can be reduced to a phasevariable if the interaction among oscillators is weak. If interaction strength is weak, amplitude changes are relatively small and play a minor role in the oscillatory dynamics. In this way, V1 neural populations can be approximated as oscillators, ‘weakly coupled’ by horizontal connections (Figure 4A). The manner in which mutually coupled oscillators adjust their phases, by phasedelay and phaseadvancement, is described by the phase response curve, the PRC (Brown et al., 2004; Canavier, 2015; Izhikevich, 2007; Kopell and Ermentrout, 2002; Schwemmer and Lewis, 2012). The PRC is important, because if the PRC of a system can be described, the synchronization behavior can be understood at a more general level and hence predicted across various conditions.
According to the theory, the synchronization of two coupled oscillators can be predicted from the forces they exert on each other as a function of their instantaneous phase difference. The amount of force is here defined as interaction strength, which is modulated as a function of phase difference by an interaction function that is closely related to the PRC (for a detailed discussion of the relationship between the two functions, please see TWCO predicts synchronization properties of V1 cortical gamma rhythms). In addition, each oscillator has an intrinsic (natural) frequency and its own source of phase noise, making the oscillators stochastic. Hence, the phase precession of two oscillators is given by:
where $\stackrel{.}{\theta}$ is the time derivative of the phase difference θ (the rate of phase precession), ∆ω the detuning (the intrinsic frequency difference), ε the interaction strength (scalar function), G(θ) the interaction function (mutual PRC), and η the combined phase noise, where $\eta \sim N\left(0,\sqrt{2}{\sigma}^{2}\right)$, see (Figure 4B). Phase noise is defined here as variation that is unrelated to interaction, which occurs for neural oscillators due to inherent instabilities of the generative mechanism (Atallah and Scanziani, 2009; Burns et al., 2010). This type of variation is distinct from measurement noise, which is unrelated to the dynamics of the system. We express ω, ε and η in units of Hz (1Hz = 2π*rad/s). The time derivative $\stackrel{.}{\theta}$ is also expressed in Hz (instantaneous frequency, IF).
Note that here, detuning ∆ω is the intrinsic or natural frequency difference between two oscillators, which is the frequency difference oscillators would have without any interaction. The measured detuning can differ from the intrinsic detuning ∆ω if the oscillators exhibit synchronization. In model simulations or while solving analytical equations, intrinsic frequencies and frequency differences are known, whereas in empirical data the intrinsic detuning ∆ω needs to be estimated from the measured detuning. Likewise, whereas ε and η are variables that can be set in analytical equations or simulations, they are not directly given in empirical data and need to be estimated. The issue of estimation is treated in the next section. Note that throughout the text, the symbols ω, ε and η are used to refer to known variables in analytical or modeling contexts, and to estimates of those variables in the description of our empirical data.
Below, we discuss the results of solving Equation 1 analytically (see Appendix for more information), which allowed us to study changes in the phasedifference probability distribution as a function of detuning ∆ω and interaction strength ε. The phasedifference probability distribution was characterized by the PLV and the mean (preferred) phase difference. The analytical solutions as a function of detuning ∆ω and interaction strength ε can be understood more easily by first considering the noisefree case. In the noisefree case (σ = 0), one can solve the equation for zeropoints (equilibrium points), meaning that the phase precession is zero ($\stackrel{.}{\theta}=0$, i.e. zero frequency difference). To reach equilibrium, the detuning ∆ω and the interaction term εG(θ) need to be counterbalanced, and three cases can be considered. First, when detuning is smaller than the interaction strength (∆ω<=ε), there is a particular phase difference at which an equilibrium can be reached. At equilibrium, there is no phase precession and thus PLV equals 1 (full synchronization). Second, when interaction strength is zero (ε = 0), the asynchronous oscillators display continuous linear phase precession and have zero PLV, with the exception of zero detuning. Third, when detuning is larger than a nonzero interaction strength (∆ω>ε, ε >0), oscillators exhibit nonlinear phase precession over time, characteristic for the intermittent synchronization regime (Ermentrout and Rinzel, 1984b; Izhikevich, 2007; Pikovsky et al., 2002, Figure 4C). The phase precession rate (instantaneous frequency difference) is determined by the detuning ∆ω, the modulation shape G(θ), and the modulation amplitude ε. Around the preferred phaserelation, the instantaneous frequency difference is reduced (‘slow’ precession in Figure 4C), whereas away from the preferred phaserelation, the instantaneous frequency is larger (‘fast’ precession in Figure 4C). For a given ∆ω and ε, a characteristic relationship can be predicted between ΔIF and Δphase (Figure 4D), indicative of the interaction function G(θ). Note that in the noiseless regime, a PLV between 0 and 1 can be obtained, varying between intermittent and full synchronization. However, including phase noise (σ > 0) has important effects on the synchronization behavior (Izhikevich, 2007; Pikovsky et al., 2002). The noise flattens the phaserelation distribution and can induce full cycles of phase precession (phase slips) that also lead to instantaneous frequency modulations. For noisy oscillators, the intermittent synchronization regime is the default regime for a large parameter range.
To show the applicability of the theory, we first reproduced the three empirical examples shown in Figure 3 by numerical simulations of Equation 1 and by varying detuning ∆ω and interaction strength ε. We assumed a sinusoidal G(θ) (see Kuramoto model, Breakspear et al., 2010; Kuramoto, 1991) and a phase variability of SD = 18 Hz (similar to our experimental data). As we did also with empirical data (see Figures 5 and 6), detuning was estimated here from the mean frequency difference at which the instantaneous frequency difference (∆IF) modulations were centered, whereas the interaction strength was estimated from the amplitude of the modulations (Figure 4D). As shown in Figure 4E–J, our simulations showed the same relation between the instantaneous frequency difference modulations and the properties of the phase difference probability distribution as observed for V1 gamma (Figure 3B–I).
To test whether the same synchronization properties could be reproduced by simulation data from a more biologically plausible model, we constructed a model consisting of two mutually coupled pyramidalinterneuron gamma network (PING) networks (Figure 3—figure supplement 2). The PING network captures essential biophysical properties of cortical gamma rhythmicity (Börgers et al., 2005; Fries, 2015; Tiesinga and Sejnowski, 2009) and can be considered a biologically plausible instantiation of an oscillator in V1. As excitatory input drives gamma frequency (Buia and Tiesinga, 2006; Jia et al., 2013b; Llinás et al., 1991; Ray and Maunsell, 2010; Roberts et al., 2013), detuning was manipulated by independently varying the excitatory synaptic input strength to the two networks. The interaction strength was manipulated by changing the crossnetwork excitatory synaptic strengths. Using this more detailed model, we fully reproduced the synchronization properties obtained with the Kuramoto oscillator model (Figure 3—figure supplement 2). This shows that the latter model, despite its simplicity, captures essential aspects of neural synchronization.
Estimating the underlying parameters and function of TWCO in observed data
To demonstrate the value of TWCO for understanding V1 gamma synchronization, we first assessed the ability of the theory to accurately predict monkey V1 recording data quantitatively (Figure 5A). Second, we tested whether we could reconstruct the Arnold tongue, which is a central prediction of the theory. The Arnold tongue describes the synchronization region in the parameter space of detuning and interaction strength (Figure 5B) and provides a general intuitive description of the gamma synchronization behavior.
To achieve the first goal, the theoretical parameters of Equation (1) need to be estimated. This equation can then be solved to predict the expected phasedifference probability distribution. Here, we were interested in two key properties of the distribution, the phaselocking value (PLV) and the mean phase difference. The theory predicts that the phasedifferencedependent modulation of the instantaneous frequency difference (∆IF(θ)) is determined by the detuning ∆ω and the interaction term εG(θ). As shown above, we consistently observed modulations in ∆IF(θ) in our experimental datasets (Figure 3). Importantly, the timeaveraged modulation of the instantaneous frequency $\mathrm{\Delta}\overline{\mathrm{I}\mathrm{F}}(\theta )$ directly relates to the deterministic term ∆ω+εG(θ), as noise is averaged out (see more in the Appendix). Based on this relation, the two parameters (∆ω and ε) as well as the shape of function G(θ) were estimated from the experimentally observed modulation of ∆IF(θ) (Figure 5A, Figure 5—figure supplement 1). Equation (1) contains a white noise process η determined by variance σ^{2} (mean = 0). The variance was determined by estimating the overall observed frequency variability in our gammaband signal (taking SNR into account, see Appendix).
Based on these theoretical considerations, we estimated ∆ω and ε separately for each contact pair between probes in each experimental condition. The interaction strength ε was estimated by the modulation amplitude of the averaged modulation in the intrinsic frequency difference $\mathrm{\Delta}\overline{\mathrm{I}\mathrm{F}}(\theta )$. The detuning ∆ω was estimated by the average of the intrinsic frequency difference $\mathrm{\Delta}\overline{\mathrm{I}\mathrm{F}}(\theta )$ computed over the full range of instantaneous phase differences [π π]. By contrast, we estimated a single G(θ) function and σ value from each monkey separately, therefore assuming stability of underlying PRCs and of the noise sources. The function G(θ) was estimated by the normalized $\mathrm{\Delta}\overline{\mathrm{I}\mathrm{F}}(\theta )$ modulation shapes. We validated the approach using phaseoscillator simulations (Figure 5—figure supplement 1). Note that the function G(θ) was estimated from data with absolute detuning of more than 4 Hz. This was done based on the observation that interaction functions became deformed when detuning was close to (see for more in Appendix). Further, it avoided smearing due to phase shifts occurring mainly within ±4 Hz. Given G(θ) and the value σ, the equation could be mathematically (analytically) solved for any values of detuning ∆ω and interaction strength ε. This means that for each contact pair and condition, we could derive precise predictions of differences in instantaneous frequency, phase relation, and phase locking (PLV) for comparison with the observed data.
A potential problem is that SNR influences both the PLV and the interaction strength estimate (problem of circularity). Further, the variables detuning and interaction might not be completely independent, due to factors like SNR. We therefore did not directly use the individual interaction strength values for comparison, but first binned contact pairs according to cortical distance (±0.25 mm). For each cortical distance, we then computed the averaged interaction strength. All contact pairs within a cortical distance bin were then assigned the same interaction strength. This step circumvented the problem of circularity and dependence of variables, but it also limited the maximum prediction accuracy that could be achieved.
To achieve the second goal, reconstructing the Arnold tongue, we mapped the observed PLV and mean phase differences as a function of detuning and interaction strength (using cortical distance binned as above) to obtain the Arnold tongue (Pikovsky et al., 2002). To demonstrate the expected shape of the synchronization region (Figure 5B), we mapped the analytically derived PLV and mean phase difference from TWCO Equation (1) in the ∆ωε parameter space. We observed a triangular synchronization region (Figure 5B) described as the Arnold tongue. This reflects the fact that stronger interaction strengths ‘tolerate’ larger detuning (∆ω<=ε). Further, a clear phase gradient along the detuning dimension can be observed. The oscillator with a higher frequency led the oscillator with a lower frequency in terms of their phases.
TWCO predicts synchronization properties of V1 cortical gamma rhythms
We then assessed whether the theory predicted the experimental gammaband PLV values recorded from V1 (focusing on layers 2–4) using the estimation procedure as described above. We estimated for each contact pair and stimulus condition their detuning values (ranging from about −6 Hz to 6 Hz) as well as their interaction strength. The phase noise parameter and the interaction function G(θ) were estimated for the two monkeys separately.
The interaction function G(θ) was estimated as being approximately a sinusoidal function (Figure 6B) with relatively symmetric negative and positive components (Akam et al., 2012). This means that phase precession was accelerated (increase in frequency) or reduced (decrease in frequency) depending on the precise phasedifference. This type of interaction function allows for robust phaselocking for negative as well as positive detuning values (see symmetric Arnold tongue below). This is because negative detuning can be counterbalanced by the positive component of G(θ) and the positive detuning by the negative component of G(θ). It is worth noting that the interaction function G(θ) is not identical with the PRC. This is because the interaction function G(θ) is the convolution of the PRC with the coupling function (Ermentrout, 1996). In the present data, the exact form of the underlying (mainly synaptic) coupling function underlying V1 gamma synchronization was out of reach, and we only estimated here its overall strength ε. This contrasts with modeling data where synaptic coupling strengths are known and the coupling function can be computed. Nevertheless, the synaptic/electrical dynamics that underlie gamma rhythms are relatively fast, and we therefore expect that the interaction function G(θ) is closely related to the PRC. Hence, whenever we use the terms PRC and G(θ) in the context of our empirical data analysis, we keep their conceptual distinction in mind while considering them similar for practical purposes.
The phase noise parameter σ was found to be relatively large (M1: σ = 19 Hz, M2: σ = 20 Hz) indicating substantial frequency variability not explained by the interaction function (likely due to inherent noise and interactions with other cortical locations). The detuning ∆ω was positively correlated with the local contrast difference (linear regression, M1: R^{2} = 0.28, M2: R^{2} = 0.25, both p<10^{−10}) and with MUA rate difference between probes (linear regression, M1: R^{2} = 0.53, M2: R^{2} = 0.36, both p<10^{−10}) in line with Ray and Maunsell, 2010. The interaction strength ε was found to be inversely correlated with the cortical distance between probes (linear regression, M1: R^{2} = 0.41, M2: R^{2} = 0.29, both p<10^{−10}, Figure 6C), in line with the known decrease of V1 horizontal connectivity with distance (Stettler et al., 2002).
To test further the idea that the interaction strength ε is a biologically meaningful measure of neural interaction more thoroughly, we repeated the analysis of interaction strength ε over cortical distance between probes with trialshuffled data. A large interaction strength ε surviving the shuffling may reveal an influence of a stimuluslocked component on ε. This permutation analysis led to populationaveraged IF modulation curves that were nearly flat, with values on average of ε = 0.31 Hz ± 0.002 in M1 and ε = 0.28 Hz ± 0.006 in M2. This is much lower than the ε values of 1–2 Hz observed without shuffling (Figure 6C). This may have been due to the fact we had only ~30 trials to shuffle per condition. This likely was not enough to obtain optimal randomization. Indeed, applying the same procedure to phaseoscillator simulations with 30 simulation trials also led to a remaining value of ε = 0.2 Hz ± 0.009. Furthermore, the higher the trial number, the closer the value got to zero (100 trials = 0.1 Hz ± 0.004, 500 trials = 0.05 Hz ± 0.002, 1000 trials = 0.03 Hz ± 0.004). Nevertheless, we cannot exclude that the small remaining nonzero value of ε after shuffling to some extent reflected a minor contribution of stimulusdependent dynamics in our data. In an attempt to empirically test interaction strength in a case where no or weak anatomical connectivity is expected, we analyzed additional V1V2 pair recordings with farremoved RFs in monkey M1 (Fig S9, GI). The interaction strength we observed (ε = 0.3 Hz) was very small, not different from shuffled trials, in line with the expected weak connectivity between involved recordings sites. Altogether, these analyses support the conclusion that ε is a biologically meaningful measure of neural interaction. Having estimated detuning ∆ω, interaction strength ε, the interaction function G(θ), and the phase noise η, we were in a position to predict the properties of synchronization for each contact pair by solving the Kuramoto equation (Figure 5A).
We found that the gamma PLV variations over single contact pairs were significantly captured by the analytical predictions as a function of ∆ω and ε (model accuracy: M1: R^{2} = 0.18, n = 7245, M2: R^{2} = 0.32, n = 7938, Figure 6D,E). This is particularly striking, given that the model predictions were derived out of first principles and single contact data were noisy. We also tested whether the model predicted variation of PLV evaluated for each single contact pair separately, where variation is induced mainly by detuning (model accuracy: M1: R^{2} = 0.27 ± 0.0002, n = 802, M2: R^{2} = 0.1 ± 0.0001, n = 882). The population means, defined as the averaged PLV values of contacts pairs with a similar detuning and cortical distance (bin size: ±0.35 Hz, ±0.3 mm), were very well predicted (model accuracy: M1: R^{2} = 0.83, M2: R^{2} = 0.86, both n = 638). To illustrate this, we plotted in Figure 6F,G the population means and the predictions for different detuning values for a single, medium interaction strength bin (M1: ε = 1.7, M2: ε = 1.6). The observed PLVs (dots) corresponded very well to the predictions (red line).
We also analyzed the mean phase difference (preferred phaserelation). A positive phase difference (phase X – phase Y) means that contact X leads (precedes in time) contact Y in terms of the phase of its oscillatory activity. Note that the temporal differences were smaller than the time scale of a full cycle, justifying the use of phase differences to indicate temporal ordering. The phase difference ranged nearly between –pi/2 to pi/2 in both M1 and M2. Again, single contact pair data was substantially captured by the analytical predictions as a function of ∆ω and ε (model accuracy: M1: R^{2} = 0.56, n = 7245, M2: R^{2} = 0.3, n=7938 Figure 6H,I). Furthermore, we tested whether the model predicted variability of phase difference evaluated for each single contact pair separately. This variability mainly represents variability induced specifically by detuning (model accuracy: M1: R^{2} = 0.52 ± 0.0002, n = 802, M2: R^{2} = 0.44 ± 0.0004, n = 882). The observed population means for different ∆ω and ε values followed the analytical predictions precisely (model accuracy: M1: R^{2} = 0.92, M2: R^{2} = 0.88, both n = 638). In Figure 6J, K, we plotted the population means and the predictions, but this time as a function of a range of detuning values for a medium interaction strength, further illustrating a good correspondence. The gamma rhythm with the higher frequency in a pair had the leading phase and the mean phase difference increased with increased detuning. To our knowledge, this is the first demonstration that phase locking values and preferred phase differences in primate cortex can be quantitatively predicted based on theoretical principles and limited knowledge of the system.
To further test the ability of TWCO to predict observed neural synchronization behavior, we plotted the observed CSDCSD gamma PLVs in V1 as a function of ∆ω and ε for both M1 and M2. In this manner, we tested whether we would observe an Arnold tongue in the V1 data, which is a synchronization region with the shape of an inverted triangle defined by its regulative parameters ∆ω and ε and a core prediction of TWCO (See Figure 5B). Figure 7A shows the observed PLV (colorcoded) plotted as a function of ∆ω and ε, revealing a structure that fitted the predicted Arnold tongue in both monkeys. As predicted, conditions of high interaction strength and low detuning showed strong gamma synchronization, whereas conditions of low interaction strength and high detuning yielded weak gamma synchronization. Notably, a model consisting of two coupled PING networks, in which interaction strength was manipulated by changing synaptic connectivity between the networks, and detuning by imposing differential excitatory drive, also yielded the Arnold tongue (Figure 7—figure supplement 1).
Using the estimated parameters, we also predicted the borders of the Arnold tongue analytically (black lines), which captured the outline of the observed Arnold tongue well. Due to intrinsic frequency variability (phase noise), the PLV values were not expected to decrease as sharply as expected from noiseless coupled oscillators (see Figure 5B). Further, in both monkeys (Figure 7A bottom), the map of mean phase difference showed a clear phase gradient across the detuning dimension as expected from the TWCO (Figure 5B). The results show that gamma rhythms with a higher frequency in a pair had the leading phase. Furthermore, for a given detuning, stronger interaction strength led to a reduction of the phase difference (see also yellow dots in Figure 3B,E,H and Figure 4E,F,G).
As an additional test of the robustness of our findings and their applicability to neural spiking data, we replicated our analysis in spikeCSD coupling measurements (see for more in the Appendix). We computed the PLV and mean phase difference between multiunit activity (MUA) recorded from a contact of one probe and the CSD recorded from a contact of another probe. MUA activity was smoothed with a Gaussian kernel (σ = 4 ms) and demeaned to obtain a continuous spike density signal that was then analyzed similar to CSD signals. As shown in Figure 7B, we observed a similar Arnold tongue structure for spikeCSD measurements. The same analysis using SpikeSpike measurements also resulted in a similar Arnold tongue structure (Figure 7—figure supplement 2).
We have thus far confined analysis to pairs in middle and superficial layers. We therefore further separately investigated interactions between deep layers (Figure 7C). CSDCSD analysis between deep contacts (L56) confirmed a similar Arnold tongue structure showing that the Arnold tongue properties do apply across the cortical layers. Our laminar probes reached also cortical area V2 lying beneath of V1 (Figure 1—figure supplement 1). We tested for V1V2 pairs whether they exhibited similar phasedependent instantaneous frequency modulations. We found that this was indeed the case (Figure 7—figure supplement 3).
The systematic variation of the phase difference between contact pairs by detuning indicates that detuning can affect the information flow between gamma rhythms (Besserve et al., 2015; Buehlmann and Deco, 2010; Cannon et al., 2014; Lowet et al., 2016). This is because spikes from a neural rhythm that leads another neural rhythm in time are more effective (Buehlmann and Deco, 2010; Cannon et al., 2014; Fries, 2015). To test this further, we mapped the main direction of Granger causal influence (see more in the Appendix) in the (CSDCSD) gamma band (X→Y vs X←Y) as a function of detuning and interaction strength. We observed that a change in the sign of detuning and phase difference was linked to a change in the direction of strongest granger causality (Figure 8).
We observed one property of synchronization that was not accounted by the model equations. We found that gamma (instantaneous) amplitude (the absolute of analytical signal) varied weakly or moderately as a function of phase difference (Figure 7—figure supplement 4) in our experimental V1 data These amplitude variations were replicated also in simulation data of two mutually coupled PING spiking networks. The gamma amplitude variation became stronger with interaction strength. It has been shown before that increased mutual entrainment of synchronizing local gamma rhythms can enhance their amplitudes (Womelsdorf et al., 2007). However, TWCO remains highly predictive even in conditions of weaktomoderate amplitude variations as long these variations do not strongly change the phase trajectory (Izhikevich, 2007; Kopell and Ermentrout, 2002; Pikovsky et al., 2002).
In Figure 9, we summarize schematically our main findings of how gamma synchronization between cortical locations is determined by their interaction strength and detuning and how it relates to the theory of weakly coupled oscillators, exemplified by the Arnold tongue. We propose that anatomical coupling is an important factor defining the interaction strength, however by itself is not sufficient to fully predict the amount of functional gammaband interactions. Critical in addition is the amount of detuning that can functionally couple or decouple anatomically connected cortical locations. The crucial combined contribution of detuning Δω and anatomical connectivity (related to ε) to synchronization is illustrated in three specific cases (Figure 9A–C), two of which yielding very low synchronization (with Δω, ε coordinates falling just outside the Arnold tongue), and one of which yielding strong synchronization (with a Δω, ε coordinate falling inside the Arnold tongue) (Figure 9D–F). Furthermore, Figure 9A–C (see arrows) illustrate that in the case of mutually anatomically coupled cortical locations, detuning influences the temporal relationship and possibly the direction of information flow between synchronized gamma rhythmic neural assemblies.
Discussion
The present study shows that gamma synchronization in awake monkey V1 adheres to theoretical principles of weakly coupled oscillators (Ermentrout and Kleinfeld, 2001; Hoppensteadt and Izhikevich, 1998; Kopell and Ermentrout, 2002; Kuramoto, 1991; Pikovsky et al., 2002; Winfree, 1967), thereby providing insight into the synchronization regime of gamma rhythms and its principles. Given the generality of the synchronization principles, they are likely to also apply to other brain regions and frequency bands.
Intermittent synchronization: the role of nonstationary frequency modulations
Our findings reveal the importance of phasedependent frequency modulations for synchronizing V1 gamma rhythms. These modulations show that a fixed and common frequency is not required for phase coordination. To the contrary, stronger nonstationary frequency modulations led to stronger synchronization, and thus to more reliable phase coordination. Frequency modulations arise naturally in the intermittent synchronization regime (Ermentrout and Rinzel, 1984b; Izhikevich, 2007; Pikovsky et al., 2002), when oscillators cannot remain in a stable equilibrium due to detuning and noise. Given the variable nature of gamma rhythms in vivo (Atallah and Scanziani, 2009; Burns et al., 2010; Ray and Maunsell, 2010; Roberts et al., 2013), intermittent synchronization is the most likely regime for their phase coordination. Although complete synchronization is not achieved in this regime, phase coordination remains sufficiently robust to influence the strength and directionality of information flow, by rendering particular phaserelations more likely than others (Battaglia et al., 2012; Buehlmann and Deco, 2010; Fries, 2015; Maris et al., 2016). The observation of nonstationary frequency modulations also has methodological implications. Gamma rhythms are often studied with stationary methods, for example spectral coherence or stationary Granger measures, yet our findings are not in line with the (weaksense) stationarity assumption (Lachaux et al., 1999; Lowet et al., 2016). Timeresolved nonstationary methods are therefore more appropriate to study the dynamics underling gamma synchronization (Bonizzi et al., 2014; Huang, 2005; Lachaux et al., 1999).
The Arnold tongue and the regulative parameters of gamma synchronization
Previous studies have established diversity in the phaselocking (Eckhorn et al., 2001; Gray and Singer, 1989; Ray and Maunsell, 2010) and in the phaserelations (Maris et al., 2016; Vinck et al., 2010) of gamma rhythms in the primate visual cortex. However, how this observed diversity in phaserelation and phaselocking is regulated was not well established. Here, we show that mainly two parameters determined gamma synchronization: the detuning Δω and the interaction strength ε. This was highlighted in the mapping of the Arnold tongue, offering a graphical understanding of how these parameters shape gammaband synchronization. Detuning represents a desynchronization force, whereas the interaction strength represents a synchronization force. In our experiment, the former was modulated by input drive differences associated with different local contrasts, and the latter by changes in connectivity strength associated with horizontal cortical distances between electrodes. Their interplay defined the resultant phaselocking strength and the preferred phaserelation between gamma rhythms. The observed role of detuning is in agreement with a previous study in the rat hippocampus (Akam et al., 2012), in which optogenetic entrainment strength and phase of gamma rhythms were dependent on the frequencydetuning. The results also agree with theoretical concepts of oscillatory interactions (Ermentrout and Kopell, 1984a; Hoppensteadt and Izhikevich, 1998; Sancristóbal et al., 2014; Tiesinga and Sejnowski, 2010). We suggest that small detuning values (mainly <∆10 Hz) reported in the present study and much larger shifts in the gamma frequencyrange (25–50 Hz to 65–120 Hz) as reported in the rat hippocampus (Colgin et al., 2009) represent different but complementary mechanisms for controlling gamma synchronization. On the one hand, only a small difference in gamma frequency will leave a possibility for synchronization while a large difference will preclude synchronization. So, large shifts in detuning open or close opportunities for synchronization. On the other hand, at small levels of detuning that offer opportunities for synchronization, small changes in instantaneous frequency will modulate the exact strength and direction of the gammamediated information flow. Hence, instantaneous frequency modulations, which define the interaction strength, reflect the overall ability of two cortical locations to engage in gammaband synchronization. These modulations are mediated by anatomical connectivity and further modified by oscillation amplitude. Hence, an important source of instantaneous V1 gamma frequency modulations is the underlying network (intermittent) synchronization process, which means that variations in gamma frequencies do not argue against a functional role of gamma synchronization (see Bosman et al., 2009; Burns et al., 2011; 2010; Roberts et al., 2013). Furthermore, we show that the shape of the instantaneous frequency modulations reflects the underlying interaction function G(θ), which in our recording data likely is closely related to the PRC (Hoppensteadt and Izhikevich, 1998; Kopell and Ermentrout, 2002; Kuramoto, 1991; Pikovsky et al., 2002; Winfree, 1967). The interaction function describes how the oscillators advance or delay each other’s phase development to coordinate their phaserelation. We observed approximately symmetric sinusoidallike functions in V1 gamma that resemble the basic function of the widelyused Kuramotomodel (Breakspear et al., 2010). This is in agreement with the biphasic PRC of gamma rhythms observed in the rat hippocampus (Akam et al., 2012) and fits with our observed symmetric Arnold tongues (Izhikevich, 2007; Kopell and Ermentrout, 2002; Pikovsky et al., 2002). Importantly, here we estimated the bidirectional interaction function G(θ). This function can be symmetric despite the presence of asymmetric individual (unidirectional) PRCs (Cannon and Kopell, 2015; Wang et al., 2013), as long as the rhythms interact approximately equally strongly, which is a plausible assumption between V1 locations.
The interaction functions we estimated here might be smoother than they really are due to limitations of our analysis arising from noise, averaging, and steps taken to reduce volume conduction. Future studies are required to characterize in more detail the (unidirectional/bidirectional) gammaband interaction functions. Unidirectionally connected neural groups, for example between certain cortical areas, might have asymmetric interaction functions and an asymmetric Arnold tongue. In this situation, a frequency difference between cortical areas (Bosman et al., 2012; Cannon et al., 2014) might be favorable for optimal information transmission.
We found smalltomoderate phasedependent variations of oscillation amplitude, which were not accounted for by the model equations. They were observed both in V1 data and PING simulations, indicating they are of biological origin. Future work is necessary to better understand their relevance. In addition, we assumed that synchronization between V1 locations emerged due to mutual horizontal interactions, yet common input fluctuations might further shape V1 gamma synchronization (Wang et al., 2000; Wiesenfeld and Moss, 1995; Zhou et al., 2013), especially for neurons with similar receptive fields. Although we did not investigate the possible effects of common input, the observation that gamma synchronization occurred between V1 locations with distinct receptive fields and with a dependence on cortical distance as expected from anatomical connectivity (Gail et al., 2000; Gieselmann and Thiele, 2008; Palanca and DeAngelis, 2005; Ray and Maunsell, 2010; Stettler et al., 2002) indicates that crosscolumnar gammaband synchronization depends strongly on direct mutual horizontal interactions (Veit et al., 2017).
Role of V1 gamma synchronization for visual processing and broader relevance
In our experiment, detuning was dependent on the local contrast difference (Ray and Maunsell, 2010; Roberts et al., 2013), known to change neural excitation in V1 (Sclar et al., 1990), while the interaction strength was dependent on the underlying horizontal connectivity strength, here varied by cortical distance (Stettler et al., 2002). Gamma synchronization is therefore informative about the sensory input (Besserve et al., 2015) and about the underlying structure of connectivity. Indeed, the frequency of gamma rhythms is modulated by various sensory stimuli (Fries, 2015) and by cognitive manipulations (Bosman et al., 2012; Buzsáki and Wang, 2012; Fries, 2015) suggesting that frequency control is a potential avenue for modulating functional gammaband coordination and information transfer (Besserve et al., 2015; Buehlmann and Deco, 2010; Lowet et al., 2016). Further, as phase lag is dependent on detuning, detuning may influence the direction of information flow among mutually coupled oscillators. This is in line with granger causality analysis in our paper (Figure 8), but also with network simulations published by us and others (Besserve et al., 2015; Buehlmann and Deco, 2010; Cannon et al., 2014; Lowet et al., 2016) showing that detuning will shape the information flow between model networks as measured by information theoretical tools (e.g. transfer entropy). Nevertheless, much more work is needed to explore the influence of detuning on directionality of information flow, and the results in the present paper are only suggestive.
The effect of detuning on synchronization was strongly modulated by interaction strength, which we demonstrated to relate strongly to the strength of horizontal connectivity. Horizontal connectivity in V1 is not only local, but also exhibits remarkable tuning to visual features, orientation being a prime example (Stettler et al., 2002). Hence, innate and learned connectivity patterns likely affect the interaction strength and hence the synchronization patterns of gamma rhythms within V1. These properties suggest V1 gamma as a functional mechanism for early vision (Eckhorn et al., 2001; Gray and Singer, 1989) by temporally coordinating local neural activity as a function of sensory input and connectivity. In agreement with previous studies (Eckhorn et al., 2001; Palanca and DeAngelis, 2005; Ray and Maunsell, 2010), V1 gamma synchronization was found to be mainly local and limited to a narrow range of frequency differences. It is therefore not likely that gamma within V1 ‘binds’ whole perceptual objects. Instead, it is more likely to bind features locally at the level of surround receptive fields. Furthermore, recent studies on the gammaband response during natural viewing (Brunet et al., 2015; Hermes et al., 2015) have found variable levels of synchronization power for different natural images. In accordance with these observations, the revealed Arnold tongue of V1 gamma implies that natural image parts with high input/detuning variability (heterogeneity) will induce no or weak synchronization, whereas parts with low input/detuning variability (homogeneity) will induce stronger synchronization. This is also in line with proposals linking gamma synchronization with surround suppression/normalization (Gieselmann and Thiele, 2008; Ray et al., 2013) and predictive coding (Vinck and Bosman, 2016). Our findings and theoretical interpretation shed new light onto the operation of gamma synchronization in the brain and will permit new and more detailed descriptions of the mechanisms by which synchronization is regulated by cognitive and sensory inputs.
Finally, we propose that the mechanism we have described for gamma synchronization in V1 also holds outside the visual cortex. Gamma synchronization across cortical areas have been observed in spite of frequency differences (Bosman et al., 2012; Gregoriou et al., 2009), which is further supported by our additional analysis of V1V2 interactions. Together, this suggests that similar principles likely operate for gammaband interareal interactions. Further, the instantaneous gamma frequency fluctuations that we have shown to be instrumental in regulating synchronization, have also been observed in the rat hippocampus by Atallah and Scanziani (2009). Their analysis suggested that these fluctuations, which reflected rapid phase shifts due to changes in excitationinhibition balance, might be critical for gammamediated information flow. Likewise, Nguyen et al. (2009) observed instantaneous frequency modulations during ripples in rodent hippocampus, revealing dynamics that may be indicative of processes related to learning and memory. These findings support our proposal that cyclebycycle modulations in frequency that regulate gamma synchronization also happen in other frequency bands and in other brain regions or structures. Nevertheless, future studies are required to test to what extent weakly coupled oscillator principles apply to different frequency bands across brain regions. Importantly, as long as the instantaneous phase of a neural rhythm can be determined, the methods used in this study can be applied. Instantaneous phase extraction has been for example applied to theta rhythms (Belluscio et al., 2012; Buzsáki, 2002) or alpha rhythms (Lakatos et al., 2005; Samaha and Postle, 2015; Schwabedal et al., 2016). In future studies, optogenetic tools (Boyden et al., 2005; Fenno et al., 2011; Zhang et al., 2007) will be highly useful to modify oscillation properties like detuning in a precise manner. Variation of sensory or cognitive variables can also be a powerful and natural way of modulating network states if enough is known about the system (e.g., Bosman et al., 2012). Interaction strength could be estimated from anatomical knowledge or manipulated by optogenetics (e.g. by targeting celltypes involved in a specific type of anatomical connectivity and varying oscillation amplitude). Aside of emerging new technological possibilities for network state modulation, a tight combination of experimental and dynamic systems theory will be critical for fruitful analysis and interpretation of neural oscillatory data.
In summary, the present paper offers the first predictive theory of synchronization, which we suggest can be used to assess the mechanisms of synchronization in various frequency bands, and to assess their contribution to diverse forms of cognition.
Materials and methods
Species used and surgical procedures
Request a detailed protocolTwo adult male rhesus monkeys were used in this study. A chamber was implanted above early visual cortex, positioned over V1/V2. A head post was implanted to headfix the monkeys during the experiment. All the procedures were in accordance with the European council directive 2010/63/EU, the Dutch ‘experiments on animal acts’ (1997) and approved by the Radboud University ethical committee on experiments with animals (Dier Experimenten Commissie, DEC).
Recording methods
Request a detailed protocolV1 recordings were made with 2 or 3 Plexon Uprobes (Plexon Inc.) consisting of 16 contacts (150 µm intercontact spacing). We recorded the local field potential (LFP) and multiunit spiking activity (MUA). For the main analysis, we used the currentsource density (CSD, (Vaknin et al., 1988)) to reduce volume conduction. We aligned the neural data from the different laminar probes according to their cortical depth and excluded contacts coming from deep V2. Layer assignment was based on the stimulusonset CSD profile (Schroeder et al., 1991a) and the interlaminar coherence pattern (Maier, 2010b). Receptive field (RF) mapping was achieved by presenting at fast rate highcontrast black and white squares pseudorandomly on a 10 × 10 grid (Roberts et al., 2013). For RF mapping we used CSD signals and spikes.
Task and visual stimuli
Request a detailed protocolThe monkeys were trained for headfixation and were placed in a Faradayisolated darkened booth at a distance of 57 cm from a computer screen. Stimuli were presented on a Samsung TFT screen (SyncMaster 940bf, 38°x30° 60 Hz). During stimulation (2 s) and prestimulus time (1 s) the monkey maintained a central eye position (measured by infrared camera, Arrington, 60 Hz sampling rate). The monkey's task was to passively gaze on a fixation point while a stimulus was shown. The monkey was rewarded for correct trials. The local stimulus contrast was manipulated in a full screen static squarewave grating (2 cycles/degree, presented at two opposite phases randomly interleaved). Contrast was varied smoothly over space such that different RFs had different contrast values. The direction of the contrast difference was parallel to the arrangement of RFs and orthogonal to the orientation of the grating. The stimulus was isoluminant with the prestimulus grey screen. We presented 9 different contrast modulation conditions (Table.S1). Cortex software (http://dally.nimh.nih.gov/index.html) was used for visual stimulation and behavioral control.
Data analysis
Request a detailed protocolWe analyzed gamma rhythms in the visual stimulation period (0.2 s  2 s). We discarded the first 200 ms to avoid stimulusonset transients. To investigate dynamical changes in the gamma phase and frequency over time, we estimated the instantaneous gamma phase and frequency using the singular spectrum decomposition of the signal (SSD [Bonizzi et al., 2014]) combined with HilbertTransform or waveletdecomposition. The phaselocking value (PLV) was estimated as the mean resultant vector length (Lachaux et al., 1999) and the preferred phaserelation as the mean resultant vector angle. For experimental data, we estimated the signaltonoise ratio (SNR) to reduce the influence of measurement noise on estimates. Phase flipping due to CSD computation was corrected.
Theoretical and computational modeling
Request a detailed protocolUsing the theory of weakly coupled oscillators, we investigated the phaselocking as well as the mean phase difference of two mutually coupled noisy phaseoscillators with variable frequency difference (detuning) and interaction strength. The stochastic differential equation was solved analytically (Pikovsky et al., 2002). The analytical results correctly predicted the numerical simulations.
Statistics
The accuracy of the theoretical predictions for the experimental data was quantified as the explained variance R^{2}.
Data availability
Request a detailed protocolExperimental data sets, modeling and analysis tools are available to all interested researchers upon request from the corresponding author. For singular spectrum decomposition visit https://project.dke.maastrichtuniversity.nl/ssd/.
Appendix
Surgical procedures
Two adult male rhesus monkeys (Macaca mulatta) were used in this experiment. Two chambers were implanted above early visual cortex, one positioned over V1/V2 and the second over V4. For the experiment reported here we used data from the V1/V2 chamber only. A head post was implanted to headfix the monkey during the experiment. All the procedures were in accordance with the European council directive 2010/63/EU, the Dutch ‘experiments on animal acts’ (1997) and approved by the Radboud University ethical committee on experiments with animals (Dier‐Experimenten‐Commissie, DEC).
Recording techniques
V1 recordings were made with Plexon Uprobes (Plexon Inc.) consisting of 16 contacts (10 µm diameter, 0.5–1 mΩ impedance, and 150 µm intercontact spacing). Three probes were inserted through a sharp guide tube, which was lowered through granulation tissue to just above the level of the dura surface. The probes were arranged in a linear manner separated from each other by ~2–3 mm. The probes were then advanced by separate microdrives (Nan Instruments LTD.). The probes were connected to headstages of high input impedance, and data were acquired via the Plexon ‘Multichannel Acquisition system’ (MAP, Plexon Inc.). The measured extracellular signal was filtered online between 150 Hz and 8 kHz to extract spiking activity and filtered between 0.7 Hz and 300 Hz to obtain the’ local field potential’ (LFP). The signal was amplified and digitized with 1 kHz for the LFP and 40 kHz for the spike signal. The data was converted from Plexon to Matlab file format and cut into trials from fixation onset to stimulus offset using the fieldtrip toolbox (Oostenveld et al., 2011). For the LFP data, the line noise was removed using the fieldtrip toolbox dft filter, which fits a sine and a cosine at 50, 100 and 150 Hertz and subtracts these components from the data. We collected 7 recording sessions in monkey M1 and 6 sessions in M2. Each recording session had on average ~590 trials in M1 and ~718 trials in M2.
Current source density (CSD)
First for extrapolating the CSD to the outermost contacts of our probes, at the top and bottom of the probe, a replica of the LFP of respectively the first and last contact was appended (Vaknin et al., 1988). The LFP was then smoothed with a Gaussian (zerophase) filter of a SD of 1.2 and range of 5 (effectively weighting signals around the centre electrode by 24% in the centre, 20% immediate neighbours, 12% 2 contacts away, 5% 3 contacts away). Then the standard CSD algorithm was applied for each contact position x, our intercontact spacing h of 150 µm and a conductivity C of 0.3 S/m:
We used CSD signals for the main analysis to reduce effects of volume conduction (see also section MUACSD and MUAMUA analysis).
Receptive field mapping
Receptive fields (RFs) were mapped using both spiking and LFP information as described in (Roberts et al., 2013). Briefly, monkeys fixated centrally while highcontrast black and white squares of sizes 0.11degree were presented pseudorandomly on a 10 × 10 grid. The locations where the spiking or the LFP response exceeded the 75th percentile of the response distribution were defined as the RF. Other than in Roberts et al. (2013), the LFP response was also used based on the envelope of the broadband gamma power (30–150 Hz) in the CSD, which we found to produce a localized result in line with spiking RFs. CSDs were computed as described above, but with a smaller Gaussian filtering of SD 0.6 and a filter range of 2, meaning that only the two neighbouring electrodes of the centre electrode had some iAcademic Pressnfluence on the RF estimate of a given contact. This was done to avoid mislocalization of RF shifts in size or position that are indicative of a shift to a different column or to V2 (Gattass et al., 1981) (see Figure 1—figure supplement 1B, rightmost plots for an example of CSD and spiking RFs with such a shift). To obtain estimates of cortical distance (in mm) between the probes we took advantage of the wellknown retinotopy of V1. We measured the distance between RF centres and calculated the cortical distance by converting differences in visual degrees using a cortical magnification factor (CMF, [Schwartz, 1980; Sereno et al., 1995]). The CMF was estimated individually for each monkey where we used the measured the physical distance between the laminar probes (fixed to the holder) before insertion into cortex (M1:~2.7 mm/deg, M2:~2.5 mm/deg).
Laminar alignment
We inserted the laminar probes on each recording day. The exact laminar positions of the probes (Figure 1—figure supplement 1) differed within and between sessions and hence we depthaligned the probes based on their stimulusevoked response and interlaminar coherence characteristics (Maier, 2010b). For depthalignment (to assign each contact a particular cortical depth value) we used the following procedure:
We computed the CSDVEP response. The different sinksource profiles were aligned using a paralleltempering technique (Frenkel and Smit, 2001). This is an iterative procedure that minimizes the squared error between all probes, shifting the position of one probe by one position on each iteration. Central to the parallel tempering algorithm is the parallel start of the procedure at multiple ‘temperatures’, each of which in our case starts with a different initial, random offset in the probes. Higher temperatures accept higher increases in error with a shift in the position of a probe. If a procedure running at a high temperature achieves a lower error than another temperature (overcoming a local minimum), it swaps the achieved shift vector with a lower temperature to find the new minimum around it. Similar to (Godlove et al., 2014) (using a genetic algorithm), we implemented a lenient maximum shift constraint between electrodes (allowing by shifts of 4 channels upwards and downwards, which for any two probes enforces a minimal overlap of 50%) to prevent trivial solutions. For our data, we used 3000 iterations at 4 different temperatures and different error tolerances per temperature (log spaced between zero and 1). The procedure showed asymptotic behaviour (no further decrease in error) at <= 1500 iterations. Note that the optimal number of iterations required for this algorithm will depend on the number of probes/sessions entered.
We then computed the within laminar probe LFP coherence matrix (Carter et al., 1973). It has been shown that there is sharp decrease in coherence around the L4/L5 border (Maier, 2010b). We chose this to refine the depth alignments of step 1 using the coherence matrix and again parallel tempering with the initial values defined by the output of step 1. An advantage of the coherence matrix is that it is a robust feature and insensitive to possible gain differences among contacts.
We manually checked for outliers of which none were found in this dataset.
Layer assignment
Channels were labelled as supragranular, granular and infragranular based on the location of the initial sinksource reversal (as established by the position of the reversal in the aligned grand average) in relation with known anatomy. We consider the position of the sinksource reversal to correspond to the edge of layers 4 and 5 (Mitzdorf and Singer, 1977; Schroeder et al., 1991a). Specifically, given our intercontact spacing of 150 micrometres and about 500 micrometres width generally used per layer (van Kerkoerle et al., 2014; Maier, 2010b), channels from this border to 450 micrometre below it were labelled infragranular, channels up to 450 micrometre above as granular, and channels 600 above it and higher as supragranular. Data were averaged within supra and granular layers or infragranular layers in agreement with the two separable sites of gammapower synchronization as indicated in the text.
Definition of the V1White MatterV2 borders
The depth probes often collected signals beyond the lower V1 layer 6 border and often reached the deep V2 infragranular layers. When the probes reached deep V2 the RFs shifted abruptly several degrees as expected form V1V2 retinotopy (Figure 1—figure supplement 1B, rightmost plots) (Gattass et al., 1981). The white matter situated between the two areas appeared relatively thin, often comprising 1–2 contacts (150300 microns).
To estimate the lower V1 Layer 6 boundary, we first used spiking RFs to determine the transition. We computed a RF centre distance measure, referenced to L4L5 border, to determine at which contact the transition to deep V2 occurred. Before the transition, often 1 or 2 contacts did not show spike RFs at all and were thus likely to represent white matter. V1 Layer 6 border was then defined as the contact with the last low RF centre distance (threshold <0.5 deg). In probes with low spiking quality; we used CSD signals (filtered in the gamma range (30–150 Hz) for determining the V1 L6 border.
Singlesession RF and CSD evaluation
For each session and probe, the CSD from fullscreen checkerboard flashes (37), the task and RF data were plotted sidebyside. CSDs from flashes and the grating onset were very similar in the initial response (data not shown). The taskdata from a single, highcontrast condition was split in an early and a later half to detect any changes in depth over the session and also compared with flash CSDs before and after the task (where available). Recordings were stable in depth according to this measure. The RF mapping was used to detect changes in the size or location of RFs over depth and to ascertain that there were no gradual drifts in RF location, indicative of a probe not inserted fully orthogonal to cortex. In cases were noticeable shifts were observed, the affected deeper channels were removed from the analysis. The final cutoff between deep V1 and white matter/V2 was determined based on the distance from the layer 4/5 reversal (see Layer assignment). This border, 450 microns below the 4/5 reversal, was typically above the level where RF shifts were observed, leading to removal of further deep channels from the analysis.
Visual stimulation paradigm
The monkeys were trained to accept headfixation and were placed in a Faradayisolated darkened booth at a distance of 57 cm from a computer screen. Stimuli were presented on a Samsung TFT screen (SyncMaster 940bf, 38° x 30° 60 Hz). The screen was calibrated to linearize luminance as function of RGB values. During stimulation and prestimulus time the monkey maintained eye position (measured by infrared camera, Arrington 60 Hz sampling rate) within a square window of 2 × 2°. This window was relatively large to allow for noise associated with the camera, recording with a second highspeed highresolution camera showed that eye position was generally held more stable than the window required. The monkey was rewarded if for keeping gaze within the eye window during the whole trial.
We aimed to manipulate gamma frequency differences between three recorded locations in V1 each separated by ~2–3 mm, corresponding to receptive fields (RF) separated by ~1 degree in visual space. The probes were arranged linearly either perpendicularly or parallel to the lunate sulcus, thus receptive fields were arranged respectively either horizontally or vertically. To manipulate gamma frequency differences, we manipulated local stimulus contrast differences in a large squarewave grating (2 cycles/degree, presented at two opposite phases randomly interleaved). Contrast was varied smoothly between the three locations. The direction of the contrast difference was parallel to the arrangement of RFs and orthogonal to the orientation of the grating. To avoid that the contrast manipulation would attract exogenous and endogenous attention (possibly appearing as an object or object boundary), we manipulated contrast differences in a repeating symmetric pattern over the entire screen. Additionally, the stimulus was isoluminant at all points and was isoluminant with the prestimulus grey screen. The contrast at the location of the centre RF, was constant over all conditions. We presented 8 levels of contrast difference and one stimulus where contrast was the same at all points. The exact contrasts differed slightly between the two monkeys since we used different screens (of the same type) which had somewhat different luminance levels. Contrast levels are given in Table 1.
We aimed to align the stimulus so that receptive fields at the three cortical locations would align with the highest, lowest and midpoint of one cycle of the contrast variation. However, RFs did not always fall exactly as we wished and there was often some variability in RFs within each probe. To get the best alignment that we could on a given session, we placed the stimulus such that receptive fields from the upper portion of the central probe fell on the midpoint between the peak and trough of the contrast variation. We then selected a stimulus where the distance between the peak and trough best matched the distance between RFs from the flanking probes. In most cases this lead to a peaktotrough distance of 2 degrees. In some cases, we used a distance of 1 or of 3 degrees. In some sessions we recorded with only two probes in V1. In those cases, the stimulus was aligned so that the midpoint was midway between the RFs of the two probes.
Most analysis was based on the measured gamma frequency rather than the stimulus contrast and so any mismatch between the stimulus contrast a particular RF received and the contrast we planned to present did not affect our conclusions. Where statistical analysis (see sections below ‘Effects of visual contrast and eccentricity on gamma frequency’) was based on stimulus contrast we took the stimulus contrast which was present at the centre of the measured RF of each single electrode contact. For Figure 1 and Figure 1—figure supplement 2 the data is shown binned by stimulus contrast values for illustration.
Effects of visual contrast and eccentricity on gamma frequency
Local stimulus contrast had a significant effect on the V1 gamma frequency (linear regression, M1: R^{2} = 0.38, n = 1179 M2: R^{2} = 0.27, n = 1134, both p<10^{−10}, slope: ~0.15 Hz/contrast, see Figure 1—figure supplement 2) in both monkey M1and M2 confirming previous studies of monkey and human visual cortex (Hadjipapas et al., 2015; Hall et al., 2005; Jia et al., 2013a; Ray and Maunsell, 2010; Roberts et al., 2013; Self et al., 2016). Stimulus contrast lead to a monotonic increase of the frequency, here measured as the mean of the instantaneous gamma frequency (similar results were obtained using the conventional frequency of the power spectral peak). Both LFP and CSD gamma gave the same result. The frequency increase was approximately linear in the range tested, however it might deviate from linearity if the whole contrast range is considered. Further, in comparison to prior studies (Hadjipapas et al., 2015; Jia et al., 2013a; Ray and Maunsell, 2010; Roberts et al., 2013), we used here wholefield gratings with local spatially varying contrast. The MUA spike rate also significantly increased with stimulus contrast (linear regression, M1: R^{2} = 0.14, n = 1179, M2: R^{2} = 0.12, n = 1134, both p<10^{−10}), which has been well established by previous work (Contreras and Palmer, 2003; Sclar et al., 1990). It suggests that the frequency change is due to a change of network excitation (Tiesinga and Sejnowski, 2009; Traub et al., 1996). We inserted laminar probes acutely into the visual cortex and the different probes, depending on their arrangement, recorded from cortical location coding for different visual eccentricities. There was also variation across sessions. It has been shown in previous work that the V1 gamma frequency is modulated by eccentricity (Lima et al., 2010; van Pelt and Fries, 2013). We confirmed these observations. The gamma frequency significantly decreased with visual eccentricity (linear regression, M1: R^{2} = 0.12, n = 1179, M2: R^{2} = 0.15, n = 1134, both p<10^{−10}). We also observed that the MUA spike rate decreased with visual eccentricity (linear regression, M1:R^{2} = 0.04, n = 1179, M2: R^{2} = 0.08, n = 1134, both p<10^{−10}) similarly to gamma frequency. Frequency differences (detuning) between all V1 pairs were here a function of both stimulus contrast, being the strongest factor, and visual eccentricity (multiple linear regression, M1: ∆contrast, R^{2} = 0.28, ∆eccentricity, R^{2} = 0.09, n = 7245; M2: ∆contrast, R^{2} = 0.25, ∆eccentricity, R^{2} = 0.11, n = 7938, all p<10^{−10}). We observed that the frequency difference was closely related to MUA spike rate difference among probes (linear regression, M1: R^{2} = 0.53, n = 7245, M2: R^{2} = 0.36, n = 7938, both p<10^{−10}) indicating that gamma frequency differences (and hence detuning) between locations are related to excitability differences. The lower excitability in more eccentric locations could reflect network differences or that stimulus, with a spatial frequency of 2 cycles/degree, was better suited to more foveal sites.
Estimation of instantaneous gamma phase, frequency and amplitude
For quantifying the phaselocking value and the preferred phase difference we relied on the reconstruction of the instantaneous phase (Picinbono, 1997). Methods based on the instantaneous phase deal better with nonstationary dynamics (than e.g. spectral coherence), which were present in the gammaband signals investigated here. The main challenge is to decompose the often complex, multicomponent measured LFP/CSD signal, into a welldefined gamma oscillatory component from which the instantaneous phase can be extracted (i.e., after a HilbertTransform or directly from a timefrequency representation (TFR),[Le Van Quyen et al., 2001]). We used a method based on the singular spectrum decomposition of the signal (SSD, see https:// project.dke.maastrichtuniversity.nl/ssd/) (Bonizzi et al., 2014). SSD is a recently proposed method for the decomposition of nonlinear and nonstationary time series (Bonizzi et al., 2012; Bonizzi et al., 2014) in a completely datadriven manner. The method originates from singular spectrum analysis (SSA), which is a nonparametric spectral estimation method used for analysis and prediction of time series. For a given signal x(t) we applied SSD for each trial separately to extract the gamma oscillatory components (SSD_{γ}). Here a short overview is presented. For more information see (Bonizzi et al., 2014). The following steps were implemented to retrieve the gamma oscillatory component SSD_{γ} (Bonizzi et al., 2012), where each iteration reproduces one component. The iteration stopped when 10 components were extracted or only 1% residual variance remained.
1. The signal x(t) is embedded giving a trajectory matrix X:
Particular to the SSD approach, the embedding dimension M is automatically estimated in a completely datadriven manner as 1.2*Fs/fmax, with fmax being the dominant frequency in the power spectral density (PSD) of x(n), and Fs the sampling frequency. The factor 1.2 allows M to cover a time span 20% larger than the average period of the wanted component (to account for a variable period).
2. The singular value composition (SVD) of the trajectory matrix X is then computed:
3. Out of the M principal components of X, an approximated version of X is obtained by selecting those principal components with a dominant frequency in the range [fmax  δf; fmax + δf], where the width of the dominant peak δf is estimated by means of a Gaussian interpolation of the power spectral density of the time series x(t). Then signal is then reconstructed by diagonal averaging. The reconstructed component signal is subtracted from the original signal and a new iteration of steps is started.
The SSD procedure results in a set of components representing rhythmic variation of the signal with different dominant frequencies. We were interested in the component which represented the gammaband. We therefore selected the component which had the largest fraction of spectral power in gamma frequencies [25 Hz60Hz] for each single trial. In the large majority of cases, there was a single dominant component representing gammaband fluctuations in the LFP/CSD signals. To get an estimate of the percentage of outlier SSD trials, we counted trials with an instantaneous frequency variation exceeding 1.5 times the interquartile range from either the 25th or 75 percentile of the distribution. We found that according to this criterion, 1.46% in M1% and 1.25% in M2 of the SSD decomposed trials could be considered as outliers, with high frequency variations indicating that for those trials the SSD decomposition was likely not optimal. As the percentage of outliers was small, we did not remove them from analysis.
For deriving the instantaneous phase of a SSD component, the Hilbert transform (HT) was applied using the Matlab implementation.
where HT(SSD_{γ}) is the HilbertTransform of the selected SSD gamma component. The HT of a realvalued signal is added as imaginary component to the realvalued signal itself to obtain the analytic signal. SSDα_{γ} is the analytical signal of the SSD_{γ}. The instantaneous phase φ can then easily be derived from the analytic signal:
Arg is the argument of the complex value SSDa_{γ}. The instantaneous frequency (IF) can be determined as the derivative of the instantaneous phase. The phases need to be unwrapped before applying the derivative. However, the IF might exhibit strong outliers if the signal is noisy. We used therefore a Savitzkygolay filter (Schafer, 2011) to smooth the phase trajectory (and hence the IF) using a polynomial fitting approach (kernel = 31 ms).
The HT is a standard approach for reconstruction of the instantaneous phase, however a problem of HT is its sensitivity to low SNR. We therefore used another approach for estimating instantaneous phase that is more robust against noise, but remains valid (Lowet et al., 2016). We approximated the instantaneous phase by using the timefrequency representation (TFR) of the signals using Morlet wavelets (Le Van Quyen et al., 2001). This approach was used mainly for estimating phaselocking strength (PLV). Morlet wavelet approach was defined as follows:
where ${W}_{SS{D}_{\gamma}}\left(t,\omega \right)$ is the wavelet coefficient of the gamma SSD component and Ψ^{*}_{t,ω} is the complex conjugate of the Morlet wavelets, both as a function of time t and frequency ω. Morlet wavelets were defined as:
Where σ defines the width of the wavelet which also defines the number of cycles (nc = 6fσ). Here we used 6 cycles. The argument of the complex wavelet coefficients gives the instantaneous phase for each frequencytime point:
Estimation of phaselocking strength and mean phase difference
The mean phase difference was defined as the mean circular phase difference between two oscillations (averaged in the complex domain), where θ = ϕ_{1} ϕ_{2}:
with a range of [π, π]. Arg is the argument function and θ is the instantaneous phase difference derived from the Hilbert transform. For estimating phaselocking we computed the phaselocking value (PLV, [Lachaux et al., 1999]) based on the instantaneous phase derived from the wavelet TFR. The PLV was computed by averaging the complex values with unit amplitude:
where T is total number of time points (trials were concatenated). The frequencies ω_{1} and ω_{2} were chosen based on the frequency of the gamma spectral power peaks of the respective contacts. The PLV ranges from 1, corresponding to full phase consistency, to 0, corresponding to fully random. Importantly, the PLV measure allows that oscillations have different frequencies (a form of crossfrequency coupling measure, [Lowet et al., 2016]). Both, HTPLV or wavelet TFRPLV gave similar results. However, the wavelet TFRPLV is more robust for SNR changes over different probes or sessions and we chose this as our preferred method for the main analysis. The main results were also not dependent on applying SSD and similar results could be obtained by combining filtering and HT or wavelet TFR on raw signals. For the MUA signals analysed (see below) we used wavelet TFRPLV on the raw MUA signals.
Correction for CSDinduced phase shifts
When applying CSD on a laminar probe the resultant signal from a given contact will likely show a constant (artificial) phase shift relative to the phase of the original LFP. This is because the CSD computes the difference among nearby LFP contacts which can change the polarity. For statistical analysis on single contact level these shifts are not problematic (as they are constant for a given contact pair) nor for the directionality measures, but it would give a scrambled picture for the Arnold tongue mapping, where all contact pairs are needed for analysis. To reduce the effect of the phase shifts, we normalized the phasedifferences for each given contact pair to the condition having the smallest frequency difference (this corresponds to a parallel translation). Hence, for CSD the phasedifference is by definition 0 at frequency difference (detuning) zero. This was done because gamma oscillations had zero phase difference at zero frequency difference shown by 1) LFPLFP analysis 2) confirmed by MUAMUA analysis. An alternative correction of the CSD phase difference using the estimated timelags from the PSI gave similar results.
SNR estimation and SNRcorrection
For experimental data it is important to consider (external) measurement noise. Measurement noise is noise that adds to the biological signal and is completely unrelated to the underlying dynamics. The amount of measurement noise is often expressed as the signaltonoise ratio (SNR). Despite the fact that the SNR from invasive LFP or MUA measurements is higher than noninvasive EEG/MEG measures, the SNR is still a limiting factor and needs to be considered for a better interpretation of the data. At low SNR, the PLV is largely underestimated. For example, a SNR of 3 can reduce the PLV more than half. Further, it also important for separating effects of true gamma amplitude from effects by SNR. A further important motivation for considering SNR correction was to be able to compare experimental PLV to the analytical predictions from the coupled oscillator equations which are SNR free.
In the data the exact amount of biological signal and external noise is unclear and needs to be approximated. We approximated the gammaband SNR by using the fact that most of the gamma power is induced by stimulation. We therefore compared gamma power during stimulation to gamma power during baseline period. The power spectra in the baseline period looked similar to 1/f indicating that the approximation is plausible. The gamma SNR was defined as follows:
To obtain PLV values that are relatively SNR independent, we simulated artificial oscillatory synchronization data using phaseoscillator equations (Lowet et al., 2016) for different SNR levels. We applied the exact same PLV estimation procedure as used for experimental data and quantified how SNR level does change the PLV estimate. The PLV estimates were compared to analytical derived expected PLV by solving the phaseoscillator equations. From these analyses we derived a SNR inverse function which gives a correction factor for the PLV measured at a particular SNR.
In addition, we performed the same procedure for the estimation of the interaction strength ε in experimental data which is also sensitive to SNR. At low SNR, the interaction strength ε is underestimated. Also here we computed a correction factor based on simulated data with different level of SNR.
Theory of weakly coupled oscillators
Detailed reviews and mathematical descriptions of the theory, also its extensions and limitations, can be found in a number of publications (Ermentrout and Kleinfeld, 2001; Hoppensteadt and Izhikevich, 1998; Kopell and Ermentrout, 2002; Kuramoto, 1991; Pikovsky et al., 2002; Winfree, 1967). According to the theory of weakly coupled oscillators, the phase evolution of two given cortical V1 locations is reduced to:
where φ_{1,2} is the phase, $\stackrel{.}{\phi}}_{1,2$ its temporal derivative, ω_{1,2} is the preferred frequency, ε_{12} and ε_{21} are the interaction strengths, H_{12} and H_{21} are the single PRCs and η_{1,2} is a phasenoise term with η_{1,2} ~ N(0, σ^{2}) N being the normal distribution. The two equations, as given in the main text Equation 1, can be further simplified to:
where θ = φ_{1} –φ_{2} is the phase difference, ∆ω = ω_{1} – ω_{2} the detuning, εG(θ)= ε_{21}H_{21}(θ)ε_{12}H_{12}(θ) the combined interaction term with ε being the interaction strength and G(θ) the mutual PRC (oddparts) and η= η_{1} –η_{2} the phase noise with η ~N(0, $\sqrt{{2\sigma}^{2}}$).
Analytical derivation of phaselocking and mean phase difference
Equation 15 is a stochastic differential equation (Langevin equation) and was solved as described in (Pikovsky et al., 2002). Equation 15 can be rewritten in the form of a FokkerPlanck equation that has been developed to give an analytical solution for the evolution of a probability distribution P of a particle influenced by a drag force (first term on the right side of the equation) and a random Gaussian noise process (second term). The drag force is here the combined systematic force of detuning ∆ω and the interaction function εG(θ):
The stationary (timeindependent) solution $\overline{P}\left(\theta \right)$ of the FokkerPlanck equation which is:
where C is a normalization constant defined by$\underset{0}{\overset{2\pi}{\int}}\overline{P}\left(\theta \right)d\theta =1$. V(θ) represents the influence of systematic force as a function of phase difference. $\overline{P}\left(\theta \right)$ is the phase difference probability distribution and describes how likely a particular phase difference is to occur. A uniform distribution means that every phase difference is equally likely and the oscillator are hence asynchronous. If the distribution approximates a delta distribution (meaning only one phase difference has nonzero probability), then the oscillators are completely synchronized. All other distributions in between signify intermittent (partial) synchronization (also called cycle slipping or phase walkthrough, [Izhikevich, 2007; Pikovsky et al., 2002]). To quantifying the narrowness of the distribution, we use the phaselocking value (the mean resultant vector length, [Lachaux et al., 1999]) defined here as:
Further, we were also interested in the mean phase difference, also described as the preferred phase difference, defined here as:
A phase difference between oscillators in neural networks implies spike timing differences. It has been shown that spiketiming is an important characteristic in addition to spike synchrony (Dan and Poo, 2004; Heitmann et al., 2013; London and Häusser, 2005; Markram et al., 2012; Masquelier et al., 2009; Tiesinga et al., 2008).
Biophysical modelling of coupled gammagenerating neural networks
To demonstrate that the results from the phaseoscillator equations are generalizable to more biophysically realistic neuronal network oscillations (see Figure 3—figure supplement 1), we simulated two coupled excitatoryinhibitory spiking neural networks generating pyramidalinterneuron gamma (PING, [Tiesinga and Sejnowski, 2009]) oscillations.
The neural voltage dynamics v were of the Izhikevichtype (Izhikevich, 2003) and defined as follows:
The coupled differential equations were numerically solved using the Euler method (1 ms step size). The networks were both composed of two types of neurons: 200 regular spiking neurons RS (a = 0.02, b = 0.2, c=65mV, d = 8) and 50 fastspiking interneuron FS (a = 0.1, b = 0.2, c=65mV, d = 2). RS were excitatory neurons and FS inhibitory neurons (ratio 4:1). The neural networks were alltoall synaptically connected. Synapses were modelled as exponential decaying functions, reset to 1 after the presynaptic neurons fired. Synaptic connection values had a maximum synaptic connection strength (max syn). The synaptic strengths were chosen from a random uniform distribution defined between the 0 and the maximal connection strength.
Within a network, RS neurons projected excitatory synaptic AMPA (decay constant = 2 ms) connections onto FS neuron (max syn = 0.45) and among themselves (max syn = 0.05). FS neurons projected synaptic GABAA (decay constant = 8 ms) connections onto RS neurons (max syn = −0.35) and among themselves (max syn = −0.2). For crossconnections between the networks, we included RS→FS connections (E→I, max syn(default)=0.015) and RS→RS connections (E→E, max syn(default)=0.007)1450011. We did not include internetwork FS→FS or FS→ RS connections to reflect that V1 horizontal connectivity is dominated by excitatory connections originating from pyramidal cells(Angelucci and Bullier, 2003; Angelucci et al., 2002; Bosking et al., 1997; Boucsein et al., 2011; Stettler et al., 2002).
The input drive to RS neurons was composed of a fixed input current to each neuron (=10), unique Gaussian input noise for a given neuron (SD ±3) and Gaussian input noise shared among neurons (SD ±1) of the same network. Thus each network received Gaussian input noise to RS neurons with the effect of inducing instantaneous frequency variation in the network over time (similar to intrinsic phase noise in the phaseoscillator model). For FS neurons, each received a fixed input current (=4) and Gaussian input noise (SD ±3). FS neurons received further excitatory drive from the RS neurons. For estimating the instantaneous phase, phase difference and frequency of the network oscillation we used a population signal defined as mean membrane voltage of all RS neurons of a given network. We simulated in total n = 697 conditions (17 coupling and 41 detuning conditions) to compare it to analytical predictions.
Instantaneous frequency modulations by phasedifference
Synchronization counteracts the phase precession by either accelerating or decelerating the precession depending on the form of the phaseresponse curve (PRC). Hence, phase difference dependent frequency modulations are expected from synchronization theory. To quantify the phase difference dependent frequency modulation in simulation/experimental data, we first computed for each pair of oscillations their instantaneous phases and their derivative (instantaneous frequency). To estimate the modulation, we computed the mean instantaneous frequency for a given instantaneous phase difference. For this, we binned the instantaneous phase difference data into equal bin sizes (bin size = 0.1 rad), and for each bin we estimated the mean instantaneous frequency, here for contact 1.
where IF the instantaneous frequency, T_{θ} is the maximal number of time points having phase difference θ, t_{θ} are individual time points with phase difference θ.
Estimation of detuning value ∆ω
The intrinsic frequency, the frequency an oscillator would have without interactions with other oscillators, could not be directly measured experimentally. The simple mean (emergent) frequency difference between oscillations will change as a function of synchronization. The stronger the synchronization, the closer the (emergent) frequency difference will be become; up to the point they are complete synchronized (common frequency).
Yet, the intrinsic frequency can be approximated from the phase difference dependent instantaneous frequency fluctuations. If there are no interactions among oscillators, the measured frequency is equal to the intrinsic frequency. However, if the oscillators synchronize, the instantaneous frequency will fluctuate as a function of the phase difference. At the preferred phase difference, the IF difference between oscillators is minimal, whereas at the antipreferred phase it is maximal. Importantly, if both the interaction strength and the PRC are similar for both oscillators, then the mean of $\mathrm{\Delta}\overline{IF}\left(\theta \right)$ will be equal to the detuning. Hence to derive the detuning, we first assumed that the interaction functions between oscillations were symmetric, which seems plausible, considering the isotropic horizontal connectivity properties in V1 (Stettler et al., 2002). The detuning value was then defined as follows:
assuming that (1) ε_{12} ≈ ε_{21} and (2) H_{12}≈ H_{21.} The validity of the approach was tested using phaseoscillator model as well as the coupled PING network model. In the former, the true detuning was a given parameter and in the latter the detuning could be measured by decoupling the two PING networks. Both modelling types showed that the detuning could be robustly retrieved, if interaction strengths were approximately symmetric. The Equation 26 can be adapted to deal with cases of interaction asymmetry (ε_{12} ≠ ε_{21}) using the individual interaction strengths ε_{12} and ε_{21} for a weighted averaging. In our experimental V1 data we observed covariation and similarity of individual interaction strengths ε_{12} and ε_{21} that is in line with the isotropic connectivity structure of V1. We observed systematic deviations from symmetry in cases when the PING networks or V1 contacts had large amplitude differences. Oscillation amplitude can influence the interaction strength ε, because a network that can send a larger amount and more synchronized spikes to another network will have a stronger influence(Fries, 2015; Tiesinga et al., 2004; Womelsdorf et al., 2007).
Estimation of interaction strength ε
A straightforward method is for each contact pair to estimate the modulation amplitude ε as the (minmax)/2 of the modulation. Even though the method works in many cases, especially for the PING simulation data, it is not very robust against SNR and has a tendency of overestimating the interaction strength as tested with phaseoscillator simulations where the true interaction is known. We therefore used another approach (used in the main analysis) based on the Fourier transform (FFT) of the modulation function:
where ω is frequency. The first Fourier coefficient is the mean offset of the modulation. Since we observed an approximately sinusoidal shape of the frequency modulation that was periodic over a phase differences of 2π, the amplitude of the modulation is captured in the second Fourier coefficient. We also included the third Fourier coefficient to capture to some extent the asymmetries observed in the modulation shape. The higher Fourier coefficients should mainly represent noise. For estimating modulation strength ε we summed the second and third Fourier coefficient and subtracted the estimated noise. This noise was assumed to be uniform across all Fourier components. It was therefore estimated as the mean amplitude of the second quadrant of N Fourier coefficients (N defined by number of phase bins of the modulation function).
This gave much more robust estimates for lower SNR data and reduced the tendency of overestimation for lower interaction strengths. Rather it had a weak tendency of underestimation (especially if the modulation shape is more asymmetric). However, for both methods the estimation of interaction strength ε systematically decreased with lower SNR.
For the experimental data, we scaled the ε values for the analytical predictions to account for the SNR in macaque V1 data. For this, as described above, we estimated the known interaction strengths of the phaseoscillator data and added different level of external noise to mimic SNR seen in monkey data. This yielded a curve giving the accuracy of the interaction strength estimates as a function of SNR. The inverse of the curve gave us the rescaling values to compensate for SNR. Based on the estimated SNR of monkey M1 and M2 data, we rescaled the estimate of ε to make it approximately SNR independent. For the main analysis conditions, with larger detuning were chosen to estimate the modulation function of instantaneous frequency by phase difference because data points were more equally distributed over different phases for these conditions. For each contact pair we used all conditions with a detuning value larger than 4 Hz and took the mean of those estimated ε values.
In our experimental data we always had cases of large detuning, owing to our experimental manipulation of contrast difference. In cases where only small detuning values are recorded, we suggest that estimations of ε can be based on the absolute instantaneous frequency differences. This will avoid cancellation of fluctuations around zero which would give severe underestimation of the true underlying interaction strength. The price will be a tendency of overestimation for very low interaction strengths.
Estimation of the G(θ)
Given that the detuning value (∆ω) and the interaction strength (ε) were estimated from the mean instantaneous frequency modulation by phasedifference ($\mathrm{\Delta}\overline{IF}\left(\theta \right)$), then G(θ) can be simply estimated by following equation:
This approach was tested using the phaseoscillator model (Figure 5—figure supplement 1), where G(θ) was a known function. Using the described approach any shape of G(θ) could be estimated from the simulation data assuming the oscillators were mutually connected with approximately similar interaction strengths (hence being symmetric).
The G(θ) describes how the rate of phase precession (equivalent to instantaneous frequency difference) between oscillator is altered as a function of phasedifference, whereas the single PRC H(θ) describes how the phase evolution of a single oscillator (=instantaneous frequency) is altered as a function of phasedifference. The H(θ) in the PING networks were asymmetric with a stronger positive (advancing) component and a weak negative (delaying) component, hence being more of the socalled PRC Type 1 (Cannon and Kopell, 2015; Schwemmer and Lewis, 2012). When PING networks were unidirectional coupled, they exhibited asymmetric Arnold tongues due to the asymmetric H(θ). Whether these applies to unidirectional coupled brain areas needs to be tested. For the main analysis, the PING networks were mutually connected with the same strength. The resultant G(θ) had therefore symmetric negative and positive components.
For making predictions of the PLV or mean phasedifference we assumed that the underlying interaction properties (the shape of the PRCs) did not change over the different conditions. We therefore used one estimation of G(θ) function for the whole dataset for each monkey 14500111450011or for all PING simulations. We assumed that the underlying interaction properties (the shape of the PRCs) did not change over the different conditions.
To obtain a G(θ) population estimation for a whole dataset we averaged single absolute G(θ) from all contact pairs that had a sufficient level of detuning (∆ω>4 Hz). We took the absolute to make G(θ) independent of sign and so avoid cancelling each other out during averaging. Taking only conditions where ∆ω>4 Hz was necessary first to assure low synchronization and therefor to have a more uniform phasedifference probability distribution to ease the estimation of the instantaneous frequency difference for all phasedifferences. Second, the minima of the G(θ) shifted in its mean (preferred) phase difference mainly within the range of −4 Hz to 4 Hz. This would lead to smearing of the obtained population estimation. Restricting to conditions of ∆ω>4 Hz ensured that the individual minima approximately overlapped.
While we used the population G(θ) for the main analysis, using a G(θ) from a single contact pair led to good prediction values of the whole dataset in many cases. However, there were also individual cases that deviated from the norm. Generally, estimation of G(θ) in contact pairs with low detuning that show very high level of synchronization is difficult as the oscillators remain constantly around their preferred phaserelation. This was the case especially in PING simulation with low phase noise levels. In this case, perturbation techniques might be more appropriate. Further, in cases of strong amplitude differences between contacts, we observed asymmetries at the level of single PRC leading to different G(θ) properties. Further, we observed small to moderate amplitude modulations as a function of phase difference (Figure 7—figure supplement 4) that might have affected the shape of G(θ) of contact pairs with strong interaction values. The use of CSD for reducing volume conduction led to additional noise due to artificial phase shifts. Given all these considerations, single G(θ) could be noisy for experimental single contact pairs. The population average G(θ) was a better representation of the interaction properties of V1 horizontal connections.
The use of IF(θ) for estimating single PRC or G(θ) needs careful considerations. However, our work shows how much important information these modulations can contain about the underlying synchronization process. Our approach can be applied to other brain regions and frequencybands to improve understanding of the underlying synchronization properties.
Estimation of noise variance σ^{2}
The estimation of the phase noise variance σ of the noise process η(t) from data is not trivial due to measurement noise and general complexity of LFP/CSD signals. The phase noise is intrinsic to the oscillatory process (dynamic noise) and relevant for the understanding of the dynamics and therefore distinct to external measurement noise. Phase noise implies variability in the instantaneous frequency of oscillators (see Equation 1315) and the overall variability of instantaneous frequency should scale with the noise variance σ^{2}. We approximated the noise variance σ^{2} by determining in the phaseoscillator model what σ^{2} value would produce the same observed instantaneous frequency difference distribution as observed in the PING or experimental V1 data. This approach was more robust than estimating the phase noise variance around the interaction function. It is important to note that the observed frequency variance is not the same as the (intrinsic) variance going into Equation 1315. This is because synchronization also counteracts the intrinsic variability. The procedure involved two main steps:
Estimate the (population average) standard deviation of the observed instantaneous frequency difference distribution of SSD gamma.
Using phaseoscillator equations find the value for σ that can reproduce the observed standard deviation of the observed instantaneous frequency difference distribution giving the observed signalnoiselevel.
Analytical predictions for PING and experimental V1 gamma
Using the Equations 17 and19 and the estimated G(θ) and noise variance σ from the data we could make predictions for any value of detuning Δω and interaction strength ε. The phase difference probability distribution $\overline{P}\left(\theta \right)$ was analytically predicted, from which we quantified the phaselocking strength (PLV, see Equation 20) and the mean phase difference (Equation 21). The predicted PLV and mean phase difference was compared to the observed PLV and mean phase difference from the data with the same detuning and interaction strength.
For PING networks we modulated the interaction strength ε by changing the internetwork connectivity strength and the detuning Δω by given different excitatory input drive to the two networks. For each simulation we had an estimate of interaction strength ε and detuning Δω.
For experimental V1 gamma, interaction strength ε was modulated by cortical distance and detuning Δω by local contrast and to a weaker extent eccentricity (see above). For each contact pair, we had their interaction strength ε and their detuning Δω.
Mapping of the Arnold tongue
For PING data, we mapped the data corresponding to their detuning and internetwork synaptic connectivity strength. For data corresponding to particular connectivity strength, we estimated the interaction strength ε and used these values for the rescaling of the yaxis, because the interaction strength ε was the parameter we wanted to compare to the theoretical model.
For experimental V1, we binned the contact pairs according to their detuning (±bin size = 0.35 Hz, bin steps = 0.2 Hz) and cortical distance (±bin size = 0.6 mm, bin steps = 0.3 mm). This average data binned according to detuning and cortical distance is what is termed population data. For data corresponding to a particular bin, we estimated the interaction strength ε and use these values for the rescaling of the yaxis. This was done to make sure that the interaction strength ε dimension was independent of the detuning ∆ω dimension, because binning directly using interaction strength ε and detuning had a potential risk of inducing dependencies between dimensions (e.g., due to SNR fluctuation) as both were based on estimation of ∆IF(θ).
Evaluation of prediction accuracy of analytical model
We estimated the accuracy of the model predictions using the coefficient of determination, R^{2}, for phaselocking strength (PLV) and the mean phase difference. Notice that here we evaluate the model accuracy without optimizing the parameters to enhance fitting.
SS_{res} is sum of square of the prediction error, the residuals of the difference between observed data and the predicted data, and the SS_{Tot} is sum of square of the demeaned observed data.
For the PING networks we observed that for both PLV (R^{2} = 0.93, n = 697) and mean phase difference (R^{2} = 0.94, n = 697) the model predictions explained a large significant part of the variance.
For experimental V1 gamma data we observed that the model predictions captured also a large significant part of the PLV (population level: M1:R^{2} = 0.88, n = 638, M2: R^{2} = 0.9, n = 638; level of single contact pairs: M1: R^{2} = 0.18, n = 7245, M2: R^{2} = 0.32, n = 7938) and mean phase difference variance in both monkeys (population level: M1: R^{2} = 0.94, n = 638, M2: R^{2} = 0.88, n = 638; level of single contact pairs: M1: R^{2} = 0.56, n = 7245, M2: R^{2} = 0.27, n = 7938). The populationlevel data represent the binned and averaged singlecontact pairs and conditions according to detuning and cortical distance (see section Mapping of the Arnold tongue).
Analysis of L2L4 and L5L6 gammaband synchronization
For the main analysis of synchronization, we limited the analysis to data recorded from L2L4 representing most the gamma power in V1 (Buffalo et al., 2011; van Kerkoerle et al., 2014; Maier, 2010b; Roberts et al., 2013; Xing et al., 2009). The lowest gamma power was observed around the L4L5 border. We observed a second gamma peak around L5L6 (van Kerkoerle et al., 2014; Xing et al., 2012) and gamma power going into deep V2. To distinguish L6 from deep V2 we used marked receptive fields shifts (as described above) as indicator for the transition from V1 to V2.
We did the exact same analysis for quantifying synchronization between pairs of L56 gamma as used for L24 gamma (Figure 7C). We could confirm the observation of an Arnold tongue in terms of PLV and mean phase difference also for the deep gamma showing that the observed synchronization properties can be generalized over different laminar compartments. We propose that calculating the PRC and Arnold tongue between various cortical locations would be a fruitful way to understand the connectivity between brain networks.
MUACSD and MUAMUA analysis
We also computed gamma PLV and the mean phase difference using multiunit activity (MUA) by computing both CSDMUA locking and MUAMUA locking (Figure 7B and Figure 7—figure supplement 2). The MUA represent a local population spike rate signal and it is thought to reflect more the ‘output’ of the network, whereas LFP/CSD represent the synaptic input of the network (Buzsáki and Schomburg, 2015; Buzsáki et al., 2012). Further, in the main analysis we estimated the synchronization gamma behaviour by using currentsource density (CSD) signals derived from our V1 16contact laminar probes. The important advantage of CSD compared to the local field potential (LFP) is the strong reduction of volume conduction which would substantially bias the PLV as well as the mean phase difference the closer the laminar probes are. The local second spatial derivation of nearby contact on the laminar probes for deriving CSD (Einevoll et al., 2013; Vaknin et al., 1988) reduced the effect of far electrical fields. However, application of CSD can likely not completely eliminate the influence of volume conduction of very near probes. Therefore, we used the more local MUA signal to test whether we can confirm a similar gamma synchronization behaviour as observed with CSD. A disadvantage of MUA signal in our recording data was its much lower SNR than LFP/CSD signals. We analysed the aggregate MUA signals of all L2L4 contacts of a single laminar probes, converted the spikes into spike densities smoothed using a Gaussian filter (σ = 4 ms). In Figure 7—figure supplement 2 the results of the MUACSD and MUAMUA analysis are illustrated showing that the Arnold tongue in terms of PLV and mean phasedifference could be observed in MUACSD as well MUAMUA signals. The results show that similar gamma synchronization behaviour in V1 can be observed at the level of CSD, representing mainly synaptic inputs, and spiking data, representing neural output. It also shows that volume conduction, already minimized for CSD signals, cannot be an influential determinant.
Instantaneous amplitude modulation by phase difference
We also investigated whether instantaneous gamma amplitude changed as a function of phase difference (Figure 7—figure supplement 4). In both, the PING networks as well as V1 gamma oscillations, we observed small to moderate amplitude modulation (up to ~15% modulation from mean amplitude). The modulations observed in the PING model looked strikingly similar to the V1 gamma amplitude modulations (compare Figure 7—figure supplement 4A/B with Figure 7—figure supplement 4C). These modulations are not expected from the weakly coupled oscillator theory, but as mentioned above (section Estimation of detuning value ∆ω), oscillation amplitude can influence the interaction strength ε as more synchronized spikes are more effective in influencing receiving neurons. This might affect synchronization behaviour in terms of phaselocking and mean phase difference, especially if amplitude modulations become more substantial. Future work should explicitly address the effect of amplitude on synchronization (Aronson et al., 1990).1450011
Granger causal analysis
We investigated Granger causal interactions (Granger, 1969) by fitting a full multivariate autoregressive model.
Where a the value of an each discrete time series $X}_{i$ at time bin t is predicted based on a linear combination of the past K time intervals (i.e. the maximum lag of the model) of itself ($n=i$) and the past of the other simultaneously recorded signals ($n\ne i$) as well as a constant offset term ($\alpha}_{i$). Finally the residuals are captured in a noise term (${\u03f5}_{i1}$). The coefficients (${\alpha}_{i},{\beta}_{ink}$) were fitted to the data by way of least squares regression.
Besides the full model in Equation 31, a second (restricted) model was fitted, where the past of signal j was left out.
Signal $X}_{j$ is said to have a Granger causal interaction with $X}_{i$ if the full model fits statistically better than the restricted model, since this suggests that the past of signal $X}_{j$ affects the present of signal $X}_{i$. The two fitting operations gave two residual noise signals, ${\eta}_{i1}\left(t\right)$ and ${\eta}_{i2}\left(t\right)$. If either of these models fitted the data significantly better, than the variance of its respective residual time series should be significantly lower than the other.
For all our analyses of Granger causal interactions of the SSD gamma components we used a maximum lag (K) of 10 time bins, that is, 10 ms. Figure 7D show the mean differences between the Fvalues calculated in a feedforward or feedback direction (i.e. ${F}_{j\to i}{F}_{i\to j}$) for signals i and j as a function of detuning and coupling strength (i.e. cortical distance).
Data availability

Data from: A quantitative theory of gamma synchronization in macaque V1Available at Dryad Digital Repository under a CC0 Public Domain Dedication.
References

Reaching beyond the classical receptive field of V1 neurons: horizontal or feedback axons?Journal of PhysiologyParis 97:141–154.https://doi.org/10.1016/j.jphysparis.2003.09.001

Circuits for local and global signal integration in primary visual cortexJournal of Neuroscience 22:8633–8646.

Amplitude response of coupled oscillatorsPhysica D: Nonlinear Phenomena 41:403–449.https://doi.org/10.1016/01672789(90)90007C

Dynamic effective connectivity of interareal brain circuitsPLoS Computational Biology 8:e1002438.https://doi.org/10.1371/journal.pcbi.1002438

Crossfrequency phasephase coupling between θ and γ oscillations in the hippocampusJournal of Neuroscience 32:423–435.https://doi.org/10.1523/JNEUROSCI.412211.2012

Singular spectrum analysis improves analysis of local field potentials from macaque V1 in active fixation taskConference proceedings : ... Annual International Conference of the IEEE Engineering in Medicine and Biology Society. IEEE Engineering in Medicine and Biology Society. Annual Conference 2012:2945–2948.https://doi.org/10.1109/EMBC.2012.6346581

Singular spectrum decomposition: A new method for time series decompositionAdvances in Adaptive Data Analysis 06:1450011.https://doi.org/10.1142/S1793536914500113

Orientation selectivity and the arrangement of horizontal connections in tree shrew striate cortexJournal of Neuroscience 17:2112–2127.

A microsaccadic rhythm modulates gammaband synchronization and behaviorJournal of Neuroscience 29:9471–9480.https://doi.org/10.1523/JNEUROSCI.119309.2009

Millisecondtimescale, genetically targeted optical control of neural activityNature Neuroscience 8:1263–1268.https://doi.org/10.1038/nn1525

Generative models of cortical oscillations: neurobiological implications of the kuramoto modelFrontiers in Human Neuroscience 4:190.https://doi.org/10.3389/fnhum.2010.00190

Optimal information transfer in the cortex through synchronizationPLoS Computational Biology 6:e1000934.https://doi.org/10.1371/journal.pcbi.1000934

Attentional modulation of firing rate and synchrony in a model cortical networkJournal of Computational Neuroscience 20:247–264.https://doi.org/10.1007/s1082700663580

Is gammaband activity in the local field potential of V1 cortex a "clock" or filtered noise?Journal of Neuroscience 31:9658–9664.https://doi.org/10.1523/JNEUROSCI.066011.2011

The origin of extracellular fields and currentsEEG, ECoG, LFP and spikesNature Reviews Neuroscience 13:407–420.https://doi.org/10.1038/nrn3241

What does gamma coherence tell us about interregional neural communication?Nature Neuroscience 18:484–489.https://doi.org/10.1038/nn.3952

Mechanisms of gamma oscillationsAnnual Review of Neuroscience 35:203–225.https://doi.org/10.1146/annurevneuro062111150444

Phaseresetting as a tool of information transmissionCurrent Opinion in Neurobiology 31:206–213.https://doi.org/10.1016/j.conb.2014.12.003

The leaky oscillator: Properties of inhibitionbased rhythms revealed through the singular phase response curveSIAM Journal on Applied Dynamical Systems 14:1930–1977.https://doi.org/10.1137/140977151

Neurosystems: brain rhythms and cognitive processingEuropean Journal of Neuroscience 39:705–719.https://doi.org/10.1111/ejn.12453

Estimation of the magnitudesquared coherence function via overlapped fast Fourier transform processingIEEE Transactions on Audio and Electroacoustics 21:337–344.https://doi.org/10.1109/TAU.1973.1162496

BookNeural oscillations as a signature of efficient coding in the presence of synaptic delaysCold Spring Harbor Labs Journals.

Response to contrast of electrophysiologically defined cell classes in primary visual cortexJournal of Neuroscience 23:6936–6945.

Modelling and analysis of local field potentials for studying the function of cortical circuitsNature Reviews Neuroscience 14:770–785.https://doi.org/10.1038/nrn3599

Temporal binding, binocular rivalry, and consciousnessConsciousness and Cognition 8:128–151.https://doi.org/10.1006/ccog.1999.0389

Type I membranes, phase resetting curves, and synchronyNeural Computation 8:979–1001.https://doi.org/10.1162/neco.1996.8.5.979

Frequency Plateaus in a Chain of Weakly Coupled Oscillators, ISIAM Journal on Mathematical Analysis 15:215–237.https://doi.org/10.1137/0515019

Beyond a pacemaker’s entrainment limit: phase walkthroughThe American Journal of Physiology 246:R102–106.

The development and application of optogeneticsAnnual Review of Neuroscience 34:389–412.https://doi.org/10.1146/annurevneuro061010113817

BookUnderstanding molecular simulation: from algorithms to applicationsAcademic Press.

Rapid feature selective neuronal synchronization through correlated latency shiftingNature Neuroscience 4:194–200.https://doi.org/10.1038/84032

Visual topography of V2 in the macaqueThe Journal of Comparative Neurology 201:519–539.https://doi.org/10.1002/cne.902010405

Switching neuronal inputs by differential modulations of gammaband phasecoherenceJournal of Neuroscience 32:16172–16180.https://doi.org/10.1523/JNEUROSCI.089012.2012

Synchrony makes neurons fire in sequence, and stimulus properties determine who is aheadJournal of Neuroscience 31:8570–8584.https://doi.org/10.1523/JNEUROSCI.281710.2011

A dendritic mechanism for decoding traveling waves: principles and applications to motor cortexPLoS Computational Biology 9:e1003260.https://doi.org/10.1371/journal.pcbi.1003260

Stimulus Dependence of Gamma Oscillations in Human Visual CortexCerebral Cortex 25:2951–2959.https://doi.org/10.1093/cercor/bhu091

HilbertHuang Transform Its Appl1–26, Introduction to the HilbertHuang Transform and its related mathematical problems, HilbertHuang Transform Its Appl.

Encyclopedia of Mathematical Physics48–53, Weakly coupled oscillators, Encyclopedia of Mathematical Physics.

Simple model of spiking neuronsIEEE Transactions on Neural Networks 14:1569–1572.https://doi.org/10.1109/TNN.2003.820440

Dynamical Systems in Neuroscience: The Geometry of Excitability and BurstingDynamical Systems in Neuroscience: The Geometry of Excitability and Bursting.

Temporal coding organized by coupled alpha and gamma oscillations prioritize visual processingTrends in Neurosciences 37:357–369.https://doi.org/10.1016/j.tins.2014.04.001

Generation of field potentials and modulation of their dynamics through volume integration of cortical activityJournal of Neurophysiology 113:339–351.https://doi.org/10.1152/jn.00914.2013

Chapter 1 Mechanisms of phaselocking and frequency control in pairs of coupled neural oscillatorsHandbook of Dynamical Systems 2:3–54.https://doi.org/10.1016/S1874575X(02)800224

Collective synchronization of pulsecoupled oscillators and excitable unitsPhysica D: Nonlinear Phenomena 50:15–30.https://doi.org/10.1016/01672789(91)90075K

Measuring phase synchrony in brain signalsHuman Brain Mapping 8:194–208.https://doi.org/10.1002/(SICI)10970193(1999)8:4<194::AIDHBM4>3.0.CO;2C

An oscillatory hierarchy controlling neuronal excitability and stimulus processing in the auditory cortexJournal of Neurophysiology 94:1904–1911.https://doi.org/10.1152/jn.00263.2005

Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchronyJournal of Neuroscience Methods 111:83–98.https://doi.org/10.1016/S01650270(01)003727

Synchronization dynamics in response to plaid stimuli in monkey V1Cerebral Cortex 20:1556–1573.https://doi.org/10.1093/cercor/bhp218

Dendritic computationAnnual Review of Neuroscience 28:503–532.https://doi.org/10.1146/annurev.neuro.28.061604.135703

Distinct superficial and deep laminar domains of activity in the visual cortex during rest and stimulationFrontiers in System Neuroscience 4:31.https://doi.org/10.3389/fnsys.2010.00031

Diverse phase relations among neuronal rhythms and their potential functionTrends in Neurosciences 39:86–99.https://doi.org/10.1016/j.tins.2015.12.004

Spiketimingdependent plasticity: a comprehensive overviewFrontiers in Synaptic Neuroscience 4:2.https://doi.org/10.3389/fnsyn.2012.00002

Oscillations, phaseoffiring coding, and spike timingdependent plasticity: an efficient learning schemeJournal of Neuroscience 29:13484–13493.https://doi.org/10.1523/JNEUROSCI.220709.2009

Cortical circuits for the control of attentionCurrent Opinion in Neurobiology 23:216–222.https://doi.org/10.1016/j.conb.2012.11.011

Laminar segregation of afferents to lateral geniculate nucleus of the cat: an analysis of current source densityJournal of Neurophysiology 40:1227–1244.

Measuring instantaneous frequency of local field potential oscillations using the Kalman smootherJournal of Neuroscience Methods 184:365–374.https://doi.org/10.1016/j.jneumeth.2009.08.012

FieldTrip: Open source software for advanced analysis of MEG, EEG, and invasive electrophysiological dataComputational Intelligence and Neuroscience 2011:1–9.https://doi.org/10.1155/2011/156869

On instantaneous amplitude and phase of signalsIEEE Transactions on Signal Processing 45:552–560.https://doi.org/10.1109/78.558469

Synchronization: A Universal Concept in Nonlinear ScienceAmerican Journal of Physics 70:655.https://doi.org/10.1119/1.1475332

Strength of gamma rhythm depends on normalizationPLoS Biology 11:e1001477.https://doi.org/10.1371/journal.pbio.1001477

The Speed of AlphaBand Oscillations Predicts the Temporal Resolution of Visual PerceptionThe Speed of AlphaBand Oscillations Predicts the Temporal Resolution of Visual Perception.

Role of frequency mismatch in neuronal communication through coherenceJournal of Computational Neuroscience 37:193–208.https://doi.org/10.1007/s1082701404957

What is a savitzkygolay filter?IEEE Signal Processing Magazine 28:111–117.https://doi.org/10.1109/MSP.2011.941097

Alphawave frequency characteristics in health and insomnia during sleepJournal of Sleep Research 25:278–286.https://doi.org/10.1111/jsr.12372

Regulation of spike timing in visual cortical circuitsNature Reviews Neuroscience 9:97–107.https://doi.org/10.1038/nrn2315

Inhibitory synchrony as a mechanism for attentional gain modulationJournal of PhysiologyParis 98:296–314.https://doi.org/10.1016/j.jphysparis.2005.09.002

Mechanisms for phase shifting in cortical networks and their role in communication through coherenceFrontiers in Human Neuroscience 4:196.https://doi.org/10.3389/fnhum.2010.00196

Analysis of gamma rhythms in the rat hippocampus in vitro and in vivoThe Journal of Physiology 493 (Pt 2):471–484.https://doi.org/10.1113/jphysiol.1996.sp021397

Cortical gamma band synchronization through somatostatin interneuronsNature Neuroscience 20:951–959.https://doi.org/10.1038/nn.4562

Gammaphase shifting in awake monkey visual cortexJournal of Neuroscience 30:1250–1257.https://doi.org/10.1523/JNEUROSCI.162309.2010

Hippocampal CA1 pyramidal neurons exhibit type 1 phaseresponse curves and type 1 excitabilityJournal of Neurophysiology 109:2757–2766.https://doi.org/10.1152/jn.00721.2012

Biological rhythms and the behavior of populations of coupled oscillatorsJournal of Theoretical Biology 16:15–42.https://doi.org/10.1016/00225193(67)900513

Spatial spread of the local field potential and its laminar variation in visual cortexJournal of Neuroscience 29:11540–11549.https://doi.org/10.1523/JNEUROSCI.257309.2009

Impact of neuronal heterogeneity on correlated colored noiseinduced synchronizationFrontiers in Computational Neuroscience 7:113.https://doi.org/10.3389/fncom.2013.00113
Article and author information
Author details
Funding
Nederlandse Organisatie voor Wetenschappelijk Onderzoek (45109025)
 Mark Jonathan Roberts
Nederlandse Organisatie voor Wetenschappelijk Onderzoek (45304002)
 Peter De Weerd
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank N Kopell, W Singer, A Bastos, P Fries, C Micheli, F Smulders, J.v.d. Eerden, J Karel, P Bonizzi, A Hadjipapas, AA.v.d. Berg for discussion. Supported by NWO VICI grant 45304002 to PDW and NWO VENI grant 45109025 to MJR. All data are stored at the Department of Psychology and Neuroscience, Maastricht University, The Netherlands. We thank the Radboud University Nijmegen for hosting our experiments, and staff of the Central Animal Facility (CDL) for expert assistance.
Ethics
Animal experimentation: All the procedures were in accordance with the European council directive 2010/63/EU, the Dutch 'experiments on animal acts' (1997) and approved by the Radboud University ethical committee on experiments with animals (Dier‐Experimenten‐Commissie, DEC).
Version history
 Received: March 8, 2017
 Accepted: August 21, 2017
 Accepted Manuscript published: August 31, 2017 (version 1)
 Version of Record published: September 25, 2017 (version 2)
 Version of Record updated: November 8, 2017 (version 3)
Copyright
© 2017, Lowet et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.