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
Decision letter

Charles E SchroederReviewing Editor; Columbia University College of Physicians and Surgeons, United States
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "A quantitative theory of γ synchronization in macaque V1" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom, Charles E Schroeder (Reviewer #1), is a member of our Board of Reviewing Editors, and the evaluation has been overseen by David Van Essen as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Bard Ermentrout (Reviewer #2).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission. We hope you will be able to submit the revised version within two months.
Summary:
This paper addressed the widely held view that, a stable, shared frequency over time is considered a condition for functional neural synchronization in the gamma range. The findings show that actually, the opposite is true. Instantaneous frequency modulations are critical for adjusting phase relations between different neural ensembles, and thus for synchronization of their oscillatory cycles; this is a key requirement for communication through oscillatory phase synchrony or coherence. The authors show that for local neuron populations oscillating at different gamma frequencies, if similar enough, these frequencies continually attracted and repulsed each other, which enabled preferred phase relations to be maintained in periods of minimized frequency difference. This dynamical interaction is predicted by the theory of weakly coupled oscillators, a fundamental mathematical principle that likely applies widely across brain regions and oscillation frequencies. The paper has a great many strengths, particularly, the outstanding effort at making this topic accessible to the general systems neuroscience audience. However, it could use improvement in several areas.
Essential revisions:
Along with the specific issues noted by the reviewers, there are a few broad concerns that the authors should address: Under reviewer 1 – 1) The issue of a negative control; 2) Negative detuning; and 3) Volume conduction. Reviewer 2's broad concerns were: 1) The authors seem to confuse the PRC and the coupling function throughout the paper; 2) Since they measure Phi, they can now compute the actual interaction function at both sites and this is even more informative than is G; 3) the way that G(φ) is plotted looks very sinusoidal possibly a consequence of the way the phase is extracted; and 4) Since the frequency varies rapidly between the two sites, it may covary between the two populations. Reviewer 3 also wanted to see a control analysis where the full analysis pipeline is applied to nonsimultaneously recorded (and so unsynchronised) trials (either by shuffling trials or combining nonsynchronously recorded probes), to ensure that the SSD step is not picking up some common effects (e.g. stimulus onset) that are inducing the deltaphasedeltaif relationship. Reviewer 3 also asked that the authors give some heuristic guidance for experimental designs to which these techniques can be applied. Presumably they need to systematically modulate both connection strength and frequency difference, as done here? Finally, reviewer 3 requested that the authors deposit the data in a suitable repository such as Dryad to reduce their future administrative burden and to ensure longterm stable access. https://opennessinitiative.org/
Reviewer #1:
This paper addressed the widely held view that, a stable, shared frequency over time is considered a condition for functional neural synchronization in the gamma range. The findings show that actually, the opposite is true. Instantaneous frequency modulations are critical for adjusting phase relations between different neural ensembles, and thus for synchronization of their oscillatory cycles; this is a key requirement for communication through oscillatory phase synchrony or coherence, and idea championed by Varela and von Stein and others, and more recently by Fries. The authors address a specific paradox: cell ensembles oscillating at different gamma range frequencies can somehow communicate. Recording with laminar multielectrodes placed in areas with different parafoveal receptive fields in monkey visual area V1, they show that for local populations running at different gamma frequencies, if similar enough, these frequencies continually attracted and repulsed each other, which enabled preferred phase relations to be maintained in periods of minimized frequency difference. This dynamical interaction is predicted by a physics theory – that of weakly coupled oscillators. This fundamental mathematical principle of synchronization through instantaneous frequency modulations likely applied widely across brain regions and oscillation frequencies. This paper has a great many strengths, particularly, the outstanding effort at making this topic accessible to the general systems neuroscience audience. We also really like the bicyclist analogy, which is also very Dutch. We do not see any serious flaws, however, the paper could use improvement in several areas. Comments on these are below.
First, a few broad issues.
1) Negative control. It would be nice if authors can show a case that can be considered reasonably as a negative control. The reason is that all figures show the cases of positive effects consistent with the theory. Evidence would be stronger if there are no effects like flat G(θ) when there is no interaction. That would also help in arguing the functional role of gamma synchronization later in Discussion. Figure 3—figure supplement 1 shows a negative dependence of interaction strength on cortical distance, which may suggest what a negative control case could be. However, the lowest values are still above 1 Hz. A more desirable control would be to use elevated gamma activity recorded simultaneously in 2 unrelated cortices. One control possibility might be recordings separated by much greater cortical distances. Assuming the authors do not have such cases, recordings during no stimulation or when animals are asleep may work. Even when control conditions are met, interaction strength still might not be zero. It is possible that there may always be some residual amount, as some interactions could remain even for null conditions. However, one might hope to find that interaction strength under control conditions is significantly smaller than values during test conditions. Also, even with a residual interaction strength, preferred phase could deviate from those during test conditions.
2) Negative detuning. Readers might surmise that for a give pair of sites, a pair of δ IFs, one positive and one negative values of same magnitudes, can be obtained, like IF1IF2 or IF2IF1, so that there could be negative detuning as shown in Figure 5. Also, 2 functions of G (θ) can be obtained, and those can be mirror images of each other across the line ΔIF = 0. Is that correct? It may help the authors' case if an example of negative IF is shown in Supplementary Information. If the answer is yes, then I could understand the leftright symmetry of Arnold tongue (as in Figure 5B). However, then why are the cases in Figure 7 not leftright symmetric for physiological data?
3) Results section: "strongly reduce the influence of volume conduction by calculating currentsource density (CSD)" Though the authors' have innovated extensively, the basic laminar electrode approach to data collection in awake monkeys is not new. There are quite a number of studies by a number of groups that used this method to analyze the areal and cellular generators of specific VEP, AEP and EEG components in monkeys between ~1990 and the present. It would be helpful to the broad readership of eLife if the present study were put in the context of that literature. Additionally, effects of volume conduction have been subject of debate over the last several years, and it would be helpful to readers to point to that. This all highlights the particular strength of the CSD approach relative to more typical LFP recordings. Note also, that there is a microscale phenomenon in volume conduction that might be relevant to this paper (Kajikawa J Neurophysiol. 2015).
Now, some more specific concerns
Introduction section: It would be helpful if the authors would be more descriptive introducing concepts of detuning and local stimulus contrast.
Results, first paragraph: Gamma power in the deepest layer is remarkable, but is analyzed rather later as an afterthought. This might be brought forward in the analysis and done in parallel with the other layer compartments.
Figure 1D Supposedly the shading represents SEM? Please state.
Results, second paragraph: It would help to state in the manuscript how the specific components were chosen. According to the Supplementary material, they chose the components with the largest fraction of spectral power in gamma frequencies (2560Hz). They also state that in most cases there was one clear component representing gammaband fluctuations. What was the percentage of cases where there was no clear component? What happens in cases where there was no one clear component?
Figure 2E: Shouldn't mean IF modulation be called mean ΔIF modulation? Do the number of dots in E reflect time samples across 33 trials presented in D? Number of observations used in this analysis should be clearly stated. Also, the relation should be quantified with statistics. Why is this based on so few trials? I understand that monkeys would usually go through several iterations of an experiment. Does this relation hold through different runs?
Results section: "The key to understanding how this dynamic relationship leads to synchronization is that phase relationships associated with lower frequency differences are maintained than phase relationships associated with higher frequency differences" This is unclear to me. Yet, understanding this link is critical. What is meant by "maintained" here?
Results, paragraph two, penultimate sentence: typo – "are (better) maintained"
Results, third paragraph: – Re Figure 3B/Figure 3—figure supplement 1: Tell us the IF for each of the V1 sites before giving us the ΔIF. Also are effects of ΔIF linear x the band?
In the same paragraph: How was the 1Hz modulation amplitude derived?
Figure 3E: this description is at odds with the Figure It looks more like 2x that (~3.5 Hz); it only becomes obvious later that the description refers to 1/2 amplitude modulation. The IF modulation in H looks larger than in E (around 4 vs. 3 respectively). This should be clarified. Also, it would be good to quantify if there is any difference in IF modulation between these two cases.
Results: "The chosen examples were representative for the 805 recorded acrossprobe contact pairs in monkey M1 and 882 pairs in monkey M2." How was this quantified?
Figure 3B, E, H readability of figure would improve if frequency difference (or max and min values used to estimate it) were directly marked on the axes. It would also help to clearly indicate N which was used to obtain these results (same for C, F, I). These results are shown qualitatively. Maybe some quantification using permutations will be beneficial.
Results, fourth paragraph: typo "hence speed is analoguesous to"
The cyclist analogy is outstanding (and very Dutch)!
Subsection “The theory of weakly coupled oscillators (TWCO): A framework for cortical gamma synchronization”: This is a point that occurred at several points throughout the paper. The authors mention horizontal connectivity, but in the communication through synchrony idea (aka CTC of Fries), hierarchical connectivity (e.g., between given RF representations in V1 and V2) is really more relevant lateral connectivity between different eccentricities in V1. Yjis point should be addressed in the Introduction and, while it is mentioned, it should be clarified in the Discussion. Through the authors downplay the notion, the present findings do seem to fit with Prof Singer's Binding Hypothesis.
In the same section: The formula includes only the detuning (Δω) and not the frequency itself. I understand that the authors would like to concentrate on gamma (and they provide clear motivation for that) but it will be good to discuss how this result might be extended to different frequency bands.
Figure 4EG: Aren't the functions supposed to have the shape of a sine wave? Particularly, Figure 4G does not look like that: downward peak look like they have almost twice the width of upward peaks.
Figure 5 legend: typo "we to we"
Subsection “TWCO predicts synchronization properties of V1 cortical gamma rhythms”: "PLV variations over single contact pairs were substantially 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 6 C,D”? It would be helpful to explain why R^{2} 0f 0.18 and 0.32 values are considered substantial.
In the sixth paragraph: Sig s6 is measured data not model result – Should be S5
In the seventh paragraph: Is there any reason for not including noise in the model?
In the same paragraph: "The results show that gamma rhythms with a higher frequency in a pair had the leading phase." This should be unpacked a little. Please state how leading and lagging phases are assigned.
In the tenth paragraph of the same section: "in the (CSDCSD) gamma band (feedforwardfeedback) as a function of detuning and interaction" This should be unpacked a lot. How are feedforward and feedback determined here?
Figure 7 AC: This Figure is not well motivated, it’s not clear why ff vs sf for some and not others. Why are B&C monkeys combined and not in A&B? Does the causality give the same directionality between IF1IF2 and IF2?
Figure 7: Analyses of MUA for Figure 7 needs more explanation. Not clear how gamma phase was defined for spike timing. Did the authors derive/filter spike density functions? And with what kernel/binsize? Also, it is better to shown raw and processed MUA, as shown for CSD in Figure 2A.
Figure 7—figure supplement 4: Does "gamma amplitude" mean the absolute amplitude of gamma band after Hilbert, or the amplitude of real signal? Also, are those amplitudes reflecting the sums of gamma in 2 sites? Please clarify. I may be wrong, but every panel in the figure shows full cycle of amplitude with a peak near zero phase difference, that makes me think it is real amplitude. Is it worth plotting gamma amplitude against interaction strength, ε?
In those plots, peak position is slightly rightward. Does it relate to phase leading/lagging? Same question again as above: If IF1IF2 and IF2IF1 are counted equally, then should the curves be symmetric around Δphase = 0?
Supplementary Information – Correction for CSDinduced phase shifts: "we normalized the phasedifferences for each given contact pair to the condition having the smallest frequency difference," by doing what? Is this parallel translations of points along axes?
In the final paragraph of the subsection: "in the case of mutually anatomically coupled cortical locations, detuning strongly influences the main direction of information flow" As presented, this is questionable.
While grand G(θ) was used to predict PLV, preferred phase difference could differ among pairs according to Figure 6. Were G(θ) values shifted along the axis of phase difference to align the downward peak of the curve to the preferred phase difference for each pair?
Subsection “The Arnold tongue and the regulative parameters of γ synchronization”. Can authors show the range of detuning? Were there cases that can show large deviations in detuning that break down interactions?
In the same subsection: while larger detuning values would decrease opportunities for synchronization across distance in V1, when they are tied strongly to stimulus information in a cell's RF, e.g., the Gabor stimuli used in Ray and Maunsell, there are still opportunities for the neurons in V1 (seeing low contrast) to synchronize with ones in V2 with overlapping RF (also seeing low contrast) and the same for V1/V2 interactions across the different contrast regions of the Gabor. This means basically frequency segregation of channels seeing different "pixels."
Supplementary information: Equation 18 in Supplementary Information duplicates Equation 17. I guess the label "18" was meant for the line: "V(θ)"?
Summary: Likely many of the concerns are simply a result of our misunderstandings. We really enjoyed reading the paper and must say, it really stimulated our thinking. It should be similarly provocative for many others in the field and will likely prove over time to be a fundamental contribution.
Reviewer #2:
I found this to be a very interesting paper and I was quite excited to see how well the theoretical predictions matched the results of the Arnold tongues. I was surprised that the interaction function was so close to a sine wave (a point that I will address shortly below) I think that the paper needs some revision and it appears that the authors are somewhat confused by some of the mathematical terminology in TWCO and I will address that below as well.
Here are my main comments
1) The authors seem to confuse the PRC and the coupling function throughout the paper. They are not the same thing. In the theory of weak coupling, there are three things:
1) The PRC or adjoint function, Z(t) which is the infinitesimal PRC and is hard to compute with real neurons, but east to compute with models (see for example the tutorial: http://www.math.pitt.edu/~bard/bardware/meth3/node15.html#SECTION00052000000000000000on weak coupling of a simple voltage gated model;
2) The interaction function, H(φ) which is defined as:
H(θ_{2}θ_{1}) = (1/T) ∫0TZ(s+θ1). C_{12}(u(s+θ_{1}),u(s+θ_{2})) ds
That is, it is the convolution of the PRC (Z) with the coupling between the oscillators, here denoted by u(t) with 1 = postsynaptic and 2 = presynaptic
3) The odd part of the interaction function which in this paper is G(θ) arises as follows
θ_{1}' = ω_{1} + ε_{12} H(θ_{2} – θ_{1}) Equation 1
θ_{2}' = ω_{2} + ε_{21} H(θ_{1} – θ_{2})
If you let φ = θ_{2} – θ_{1}, the phase difference and if ε_{12}=ε_{21}, then
φ' = Δω + εG(φ) Equation 2
where G(φ) = H(φ) – H(φ) is proportional to the odd part of H(φ)
Equation 2 is what the authors measure; it is not the PRC. I just want to make that clear.
2) Since they measure φ, they can now compute the actual interaction function at both sites and this is even more informative than is G. By measuring θ_{1}' and θ_{2}', and plotting this against φ, equation 1 shows that
θ_{1}' = ω_{1} + ε_{12} H(φ)
θ_{2}' = ω_{2} + ε_{21} H(φ)
thus, in addition to G, they will be able to measure ω_{1} + ε_{12} H(φ) and ε_{21} H(φ). This will give them more information about the direction of coupling and also will give them the full value of H. Indeed, there is nothing here that requires the two H's to even be the same. Thus, I would suggest plotting θ_{1}', θ_{2}' vs φ as well as φ' vs φ.
3) G(φ) is plotted looks very sinusoidal. Is this a consequence of the way the phase is extracted? I suggest taking a model pair of oscillators, such as Hodgkin huxley and coupling them with synapse along with noise and heterogeneity in coupling and in frequency. Then use the method that is used for the data to extract G(φ) and H(φ) and then compare this to the actual determistic values found by setting the noise to zero and computing H and G as in the above. This would assure me that the method of using Hilbert transforms isn’t somehow removing all the Fourier modes beyond the lowest.
4) Since the frequency varies rapidly between the two sites, does this variability covary between the two populations? I ask this because one can take two uncoupled oscillators and apply a broadband correlated signal to them and they will produce a very nice peak in their phasedifference histogram. For example, see Nakao, Arai, Kuramoto PRE 2005, or other papers by Arai & Nakao, or for a neuroscience version, Zhou, Burton, et al. 2013 Frontiers in Comp Neuro.
Reviewer #3:
This study applies the wellknown theoretical model of weakly coupled oscillators to experimental data from macaque V1. This is achieved through an elegant experimental design, which obtains simultaneous paired recordings over a large range of coupling strengths and frequency differences, to allow fitting of the oscillator model. While the idea is simple, the work is technically impressive and has broad important applications for interpretation and analysis of phase relations between brain areas. The writing is extremely clear; the authors have done an excellent job of conveying a technically complex piece of work in a way that should be accessible to a broad audience.
The results are impressive, but I would like to see a control analysis where the full analysis pipeline is applied to nonsimultaneously recorded (and so unsynchronised) trials (either by shuffling trials or combining nonsynchronously recorded probes). This would ensure that the SSD step is not picking up some common effects (e.g. stimulus onset) that is inducing the deltaphasedeltaif relationship e.g. in Figure 2E. (I don't think full permutation inference is necessary for all results, but just some confirmation that indeed the actual experimental data pipeline does produce flat deltaphase / deltaif curves)
The authors propose that the method of fitting TWCO model to data and associated analyses such as Arnold tongue plots might be applicable to other brain regions or frequency band relationships. It would be nice if they could give some heuristic guidance for experimental designs to which these techniques can be applied. Presumably they need to systematically modulate both connection strength and frequency difference, as done here?
The statistical analysis and Figure 6CD show the overall average over all electrode pairs and stimulus conditions. Presumably the experimental contrast already induces a range of PLVs for each specific recording pair, so perhaps they could be evaluated statistically as well (e.g. scatter with pair PLV R2 as a function of distance between the pairs, or something along similar lines).
The authors indicate that data and code are available on request. I would strongly urge them then to deposit the data in a suitable repository such as Dryad to reduce their future administrative burden and to ensure long term stable access. https://opennessinitiative.org/
https://doi.org/10.7554/eLife.26642.025Author response
Reviewer #1:
[…] 1) Negative control. It would be nice if authors can show a case that can be considered reasonably as a negative control. The reason is that all figures show the cases of positive effects consistent with the theory. Evidence would be stronger if there are no effects like flat G(θ) when there is no interaction. That would also help in arguing the functional role of γ synchronization later in Discussion. Figure 3—figure supplement 1 shows a negative dependence of interaction strength on cortical distance, which may suggest what a negative control case could be. However, the lowest values are still above 1 Hz. A more desirable control would be to use elevated γ activity recorded simultaneously in 2 unrelated cortices. One control possibility might be recordings separated by much greater cortical distances. Assuming the authors do not have such cases, recordings during no stimulation or when animals are asleep may work. Even when control conditions are met, interaction strength still might not be zero. It is possible that there may always be some residual amount, as some interactions could remain even for null conditions. However, one might hope to find that interaction strength under control conditions is significantly smaller than values during test conditions. Also, even with a residual interaction strength, preferred phase could deviate from those during test conditions.
The interaction strength (phasedependent IF modulations) reflects the strength of synchronization (even if rhythms are not phaselocked due to high detuning), a measure predicted within the framework of weakly coupled oscillators. If the measure does truly relate to interaction strength, it should depend on anatomical connectivity strength. We did not directly measure anatomical connectivity, but we used the fact that horizontal connectivity decreases as a function of cortical distance in V1. The observation that measured interaction strength decreased with cortical distance was for us a strong indication that the measure has biological relevance. We point this out more clearly now in the text (as suggested by you and reviewer 3) and also include the relevant plot in the main manuscript (Figure 6C):
“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).”
This biological relevance is strengthened further by two analyses suggested by the review (see below).
We agree with the reviewer that a case of no (mutual) connectivity would be beneficial to check whether some interaction strength is remaining. Unfortunately, we did not have such data. However, we have data that get close to this case (Figure 7—figure supplement 3). The laminar probes not only covered the full range of V1 cortical depth, but also a significant part of V2, situated just beneath V1. The receptive fields (RF) of the ‘deep’ V2 were relatively far away from the V1 RF (~3.5°). We therefore analyzed V1V2 pairs with separate RF positions, to measure the resultant IF modulations amplitude. We found that V1V2 pairs had similar IF modulations as observed for V1V1 pairs. The analysis revealed that locations with a large distance (in terms of RF distance) showed much less phasedependent interactions getting close to zero (interaction strength ε= 0.3Hz). Whether some modulations remained for locations without any connectivity, can however not be tested with this data.
As an additional control, we did the same analysis with shuffled trials to test whether stimuluslocked dynamics had an influence on the IF modulations results. We found that shuffling trials did strongly reduce systematic IF modulations. We have included a report of these two additional analyses in the main text:
“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. […] Altogether, these analyses support the conclusion that ε is a biologically meaningful measure of neural interaction.”
2) Negative detuning. Readers might surmise that for a give pair of sites, a pair of δ IFs, one positive and one negative values of same magnitudes, can be obtained, like IF1IF2 or IF2IF1, so that there could be negative detuning as shown in Figure 5. Also, 2 functions of G (θ) can be obtained, and those can be mirror images of each other across the line ΔIF = 0. Is that correct? It may help the authors' case if an example of negative IF is shown in Supplementary Information. If the answer is yes, then I could understand the leftright symmetry of Arnold tongue (as in Figure 5B). However, then why are the cases in Figure 7 not leftright symmetric for physiological data?
For each contact pair we arbitrarily labeled the contact points as IF1 and IF2. However, the assignment for a given pair remained the same across all stimulus contrast conditions, thus inducing a range of frequency differences. A single contact pair therefore had negative as well as positive detuning values. For illustration purposes we showed only positive detuning. We added now a clarifying sentence:
“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 1E).”
We also added further clarifications in the figure legend of Figure 2:
“Note further that ΔIF variations are the result of IF variations occurring simultaneously at the two contact points that together constitute a contact pair.”
Later on in the same figure legend we added that in panel E:
“the ΔIF values tend to be positive, which is related to the sign of the contrast difference and resulting detuning. If for the same pair the contrast difference had been reversed, ΔIF values would have tended to be negative.”
Panel F was newly added in Figure 3—figure supplement 1, and shows ΔIF modulation curves for different detuning values ranging from negative to positive for both monkeys.
Regarding the symmetry of the Arnold tongue: The Arnold tongue of the phaseoscillator model is symmetric in Figure 5, because the underlying interaction function is symmetric (sinusoid with equal negative and positive parts) and because the values were analytically exactly derived, so there were no numerical estimation errors. In Figure 7 we see slight deviations from symmetry because of noise and limited amount of experimental data. Symmetry fits with the experimental observation that interaction function had a positive and negative component and each member of a pair had IF modulations indicating interactions are mutual.
3) Results section: "strongly reduce the influence of volume conduction by calculating currentsource density (CSD)" Though the authors' have innovated extensively, the basic laminar electrode approach to data collection in awake monkeys is not new. There are quite a number of studies by a number of groups that used this method to analyze the areal and cellular generators of specific VEP, AEP and EEG components in monkeys between ~1990 and the present. It would be helpful to the broad readership of eLife if the present study were put in the context of that literature. Additionally, effects of volume conduction have been subject of debate over the last several years, and it would be helpful to readers to point to that. This all highlights the particular strength of the CSD approach relative to more typical LFP recordings. Note also, that there is a microscale phenomenon in volume conduction that might be relevant to this paper (Kajikawa J Neurophysiol. 2015).
We thank the reviewer for this remark and agree fully. We edited the text to more clearly indicate the value of using CSD in the context of existing literature:
“The close positions of recording sites may have led to a contribution of volume condition to synchronization measures. […]The success in reducing volume conduction using CSD favors its use over LFP for spectral analysis at high spatial resolution.”
We did not know the article from Kajikawa and Schroeder and we have read it with a lot of interest. We agree that a deeper understanding of LFP and CSD across layers is an important topic. In our comparison between deep and middlesuperficial layers we have used CSD, which should have limited farreaching crosslaminar contributions present for LFP. We cited the article in this context.
Now, some more specific concerns
Introduction section: It would be helpful if the authors would be more descriptive introducing concepts of detuning and local stimulus contrast.
We improved the text that introduces detuning and local stimulus contrast:
“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 1) was manipulated by varying the sign and amplitude of the spatial variation in contrast.”
And:
“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.”
Results, first paragraph: gamma power in the deepest layer is remarkable, but is analyzed rather later as an afterthought. This might be brought forward in the analysis and done in parallel with the other layer compartments.
We agree that it is very interesting. It does fit observations from (Xing et al., 2012) and (van Kerkoerle et al., 2014) who found also gamma in the deepest layer (we cited them now in the text). We found however that a detailed laminar analysis would go beyond the scope of this study. The main analysis was focused on L24 contacts that represented the largest gamma source (‘compartment’, (Maier et al., 2010)). In Figure 7C, we extended our analysis to deep layer gamma (putatively Layer 6) to show that it did not differ qualitatively from middlesuperficial layers. We believe adding information about the deep layers in Figure 6, and give more detailed description of analysis in text would overload the paper and not add much conceptually. In addition, a more detailed and careful analysis is needed for understanding possible small differences among layers as well as the gammaband interactions occurring within a cortical column. We feel this goes beyond the scope of this paper.
Figure 1D Supposedly the shading represents SEM? Please state.
Thanks. It is corrected in legend of Figure 1.
Results, second paragraph: It would help to state in the manuscript how the specific components were chosen. According to the Supplementary material, they chose the components with the largest fraction of spectral power in gamma frequencies (2560Hz). They also state that in most cases there was one clear component representing gammaband fluctuations. What was the percentage of cases where there was no clear component? What happens in cases where there was no one clear component?
For all trials we chose the SSD component with strongest power in the gamma range (2560HZ). It is true that trials with low gamma power (or low SNR) can lead to unclear SSD components that cover the gamma range. We estimated the percentage of trials where nonoptimal decomposition occurred by evaluating the number of trials with untypical high frequency variation. We estimated the percentage of outlier SSD trials by counting trials with frequency variation exceeding 1.5 times the interquartile range from either the 25^{th} or 75^{th} percentile of the distribution. The resulting estimates of outliers were for M1: 1.46% and for M2: 1.25% of the trials. We describe this now in subsection “Estimation of instantaneous gamma phase, frequency and amplitude” of the Appendix. We did not eliminate these outliers.
Figure 2E: Shouldn't mean IF modulation be called mean ∆IF modulation? Do the number of dots in E reflect time samples across 33 trials presented in D? Number of observations used in this analysis should be clearly stated. Also, the relation should be quantified with statistics. Why is this based on so few trials? I understand that monkeys would usually go through several iterations of an experiment. Does this relation hold through different runs?
Thanks. We corrected ‘mean IF modulation’ to mean ∆IF modulation in Figure 2E
The number of dots in Figure 2E reflects indeed the number of time samples for each of the 33 trials included (1800 x 33 = 59400 points). This now clarified in the legend of Figure 2. We believe we give sufficient statistics in the remainder of the manuscript.
These are trials from one stimulus condition and from one session. We had in total 22 stimulus conditions from which 9 conditions were of relevance for this manuscript. The monkey worked in general for a few hours on this task.
Results section: "The key to understanding how this dynamic relationship leads to synchronization is that phase relationships associated with lower frequency differences are maintained than phase relationships associated with higher frequency differences" This is unclear to me. Yet, understanding this link is critical. What is meant by "maintained" here?
Results, paragraph two, penultimate sentence: typo – "are (better) maintained"
We agree that this sentence was not clear. We changed it to: “are maintained longer over time (slower precession) than”
Results, third paragraph: Re Figure 3B/Figure 3—figure supplement 1: Tell us the IF for each of the V1 sites before giving us the ∆IF. Also are effects of ∆IF linear x the band?
We believe that Figure 2 clearly illustrates the process of how in a single trial we first derive gamma band signals for the different V1 sites (e.g. (see Figure 2A). We believe it is intuitive that from these two signal an instantaneous phase difference can be computed that varies over time, that in turn can lead to a specific instantaneous frequency difference, that also varies over time (Figure 2B,C). The mean ∆IF modulation is critical for making the predictions that are then compared to the theoretical framework, and hence the writing, figures and analysis are focused on that. It is willknown that there are IF variations at single recording sites (e.g. Roberts et al., 2013), but these variations are not relevant for our analysis. Nevertheless, examples of single site V1 IF modulations can be found in the Appendix (Figure 3—figure supplement 1).
In the same paragraph: How was the 1Hz modulation amplitude derived?
It is approximately (maxmin)/2. We derived it using a Fourier approach (using the second and third Fourier component). We added this information in paragraph three of the Results section, and referred for further explanations to the Appendix (section ‘Estimation of interaction strength ε’).
Re Figure 3E: this description is at odds with the Figure It looks more like 2x that (~3.5 Hz); it only becomes obvious later that the description refers to 1/2 amplitude modulation.
We corrected this in the text.
The IF modulation in H looks larger than in E (around 4 vs. 3, respectively). This should be clarified. Also, it would be good to quantify if there is any difference in IF modulation between these two cases.
Notice that the SE shaded area might mislead in coming to a visual judgment. There is actually no difference in IF modulation. In the main text we state for Figure 3H that the amplitude modulation was 1.8Hz, which is 3.6Hz from peak to through. Later on we stated that for Figure 3E:
“the magnitude of the instantaneous frequency modulation did not change”. We now have added between brackets explicitly the value of 1.8Hz in the text.
Results: "The chosen examples were representative for the 805 recorded acrossprobe contact pairs in monkey M1 and 882 pairs in monkey M2." How was this quantified?
We changed the sentence to indicate that the observations, as illustrated in the examples, were characteristic for the whole dataset as shown in the main statistical analysis.
Figure 3B, E, H readability of figure would improve if frequency difference (or max and min values used to estimate it) were directly marked on the axes. It would also help to clearly indicate N which was used to obtain these results (same for C, F, I). These results are shown qualitatively. Maybe some quantification using permutations will be beneficial.
Thanks for the suggestion. We tried to directly mark the frequency difference on the axis, but we found this figure then looked too dense. We would therefore rather like to keep it as is.
This figure was devised as an introductory/illustrative figure. Detailed quantifications are presented in Figure 6 and corresponding text.
Results, fourth paragraph: typo "hence speed is analoguesous to"
Thanks for pointing this out. Corrected.
BTW the cyclist analogy is outstanding (and very Dutch)!
Thanks. There might have been indeed an influence by Dutch culture
Subsection “The theory of weakly coupled oscillators (TWCO): A framework for cortical gamma synchronization”: This is a point that occurred at several points throughout the paper. The authors mention horizontal connectivity, but in the communication through synchrony idea (aka CTC of Fries), hierarchical connectivity (e.g., between given RF representations in V1 and V2) is really more relevant lateral connectivity between different eccentricities in V1. Yjis point should be addressed in the Introduction and, while it is mentioned, it should be clarified in the Discussion. Through the authors downplay the notion, the present findings do seem to fit with Prof Singer's Binding Hypothesis.
We agree with the reviewer that it is an important issue that needs more discussion. We are thankful to the reviewer to promote this discussion.
We question whether there is a clear and founded distinction between ‘binding theories’ and CTC, particularly at the level of synchronization principles involved. Even though anatomical connectivity between cortical areas is different from within a cortical area, we think that coupled oscillator concepts (detuning, IF modulations) likely apply also for betweenarea synchrony. Our ongoing analysis of V1V2 interactions (based on data in Roberts et al., 2013) supports this view as well as published reports that synchrony between cortical areas is observed despite frequency differences (Bosman et al., 2012).
It was not our intention to downplay the Binding by synchrony (BBS) hypothesis as proposed by Prof. Singer. Our data does support a local type of BBS (however not a global type BBS), where feature integration can be supported by gamma synchrony at the level of surround receptive fields. We added a Line to express that we think that gamma do support feature binding at a local spatial level. It reads now:
“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. […]It is therefore not likely that gamma within V1 ‘binds’ whole perceptual objects, but rather binds features locally at the level of surround receptive fields.”
We now explicitly address the extension of interareal interactions in the text:
“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.”
In the same section: The formula includes only the detuning (Δω) and not the frequency itself. I understand that the authors would like to concentrate on gamma (and they provide clear motivation for that) but it will be good to discuss how this result might be extended to different frequency bands.
We completely agree. We added a paragraph stating that the analysis approach can be extended to other cortical/subcortical areas and to different frequencybands. This paragraph reads as follows:
“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 α rhythms (Lakatos et al., 2012; Samaha and Postle, 2015; Schwabedal et al., 2016).”
Figure 4EG: Aren't the functions supposed to have the shape of a sine wave? Particularly, Figure 4G does not look like that: downward peak look like they have almost twice the width of upward peaks.
The underlying interaction function was found to be approximatively a sinusoid. In our original Discussion, we used terms such as “sinusoidallike”. We mention now more explicitly that our analysis might have led to a more smoothed picture of the interaction function than it really is on single trials (due to noise, averaging, CSDrelated phase shifts).
At low detuning, when the IF modulations gets close to zero, the IF modulations get deformed and more asymmetric. We observed this in simulated phaseoscillator model (where the PRC is exactly known), in PING network simulations and in V1 data. At higher detuning the IF modulations converge more and more to the true interaction function shape. That is the reason we estimated the interaction function from IF modulations with sufficient large detuning (>4Hz). It was written in the Supplementary Materials, but we clarified it now also in the main text with the following text:
“Note that the function G(θ) was estimated from data with absolute detuning of more than 4Hz. This was done based on the observation that interaction functions became deformed when detuning was close to (see supplementary materials for further details). Further it avoided smearing due to phase shifts occurring mainly within ± 4Hz.”
Figure 5 legend: typo "we to we"
Thanks. Corrected
Subsection “TWCO predicts synchronization properties of V1 cortical gamma rhythms”: PLV variations over single contact pairs were substantially 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 6 C,D”? It would be helpful to explain why R^{2} 0f 0.18 and 0.32 values are considered substantial.
We agree with the reviewer that the term is vague. However, the values need to be considered in the context that the model predictions are derived from firstprinciple (without direct fitting) as well as that the single contact PLV estimates are noisy. The much higher explained variance values for population means support the view that the model predictions capture an important part of the underlying biological variance.
To address this, we removed the term and replaced it with a sentence to give context to the statistical values. It reads:
‘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 6 C, D), which is considerable taking into account that the model predictions were derived out of first principles and that the single contact data used were noisy.’
In the sixth paragraph: Sig s6 is measured data not model result – should be S5
Thanks. Corrected
In the seventh paragraph: Is there any reason for not including noise in the model?
As shown in Figure 5A, we did include noise in the model. We think the reviewer is referring to the borders (black lines) in Figure 5B, indicating the expected shape of the Arnold tongue (Δω=ε). In the noiseless case, beyond this point oscillators are not completely synchronized anymore and show phase precession. In the presence of noise as is the case in Figure 5B, this sharp transition from complete synchronization to incomplete synchronization does not exist per se. However, this point (or border) is still useful to understand where the synchronization region lies and where PLV values sharply decreases:
“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).”
In the same paragraph: "The results show that gamma rhythms with a higher frequency in a pair had the leading phase.” This should be unpacked a little. Please state how leading and lagging phases are assigned.
A leading phase was defined as the rhythm A having a positive phase difference with rhythm B (phase A –phase B). It is equivalent of rhythm A reaching the oscillation maxima or minima slightly early in time.
For PING networks it meant that neurons of rhythm A spiked earlier than rhythm B (for reviewer’s interest, but not mentioned in paper).
The temporal differences are smaller than the time scale of a full cycle justifying the temporal ordering based on phase differences.
This interpretation was further supported by granger analysis (Figure 8).
Phaseslopeindex (PSI) analysis (not included in the manuscript) did additionally support our interpretation as did crosscorrelation analysis.
We added:
“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.”
In the tenth paragraph of the same section: "in the (CSDCSD) gamma band (feedforwardfeedback) as a function of detuning and interaction” This should be unpacked a lot. How are feedforward and feedback determined here?
The term ‘feedback’ and ‘feedforward’ are indeed difficult terms in this context. We changed them with more descriptive labels (feedforward: X→Y, feedback: X ←Y).
We clarify this in the text and in the main Figure 8.
Figure 7 AC: This Figure is not well motivated, it’s not clear why ff vs sf for some and not others. Why are B&C monkeys combined and not in A&B? Does the causality give the same directionality between IF1IF2 and IF2?
The main finding of the figure is the L24 Fieldfield Arnold tongue mapping, shown in Figure 7A, which is a result that combines the essence of the previous data figures. Our data analysis was focused on L24 contacts, where the main gammaband source is localized. As in previous figures we show the data for the monkeys separately.
B is an extension showing that the same results can be found using spikefield locking. Our motivation was to show that the results can be reproduced by using spiking activity. To limit the number of subplots and in light of all the previous evidence in the paper that the monkeys show comparable results, we opted to show plots combined for both monkeys. The same applies for C.
To address the reviewer’s question on granger causality, we took panel Figure 7D, and separated it into a figure for the two monkeys. The changes in phase difference with detuning (Figure 7, bottom row) were entirely in line with the changes in directed influence computed with the Granger Causality method (Figure 8).
Figure 7: Analyses of MUA for Figure 7 needs more explanation. Not clear how gamma phase was defined for spike timing. Did the authors derive/filter spike density functions? And with what kernel/binsize? Also, it is better to shown raw and processed MUA, as shown for CSD in Figure 2A.
We did our analysis on aggregated smoothed (Gaussian kernel) MUA signals and treated them as continuous signals. That had the advantage to do the exact same analysis pipeline as used for CSD signals. The smoothing kernel is now explicitly mentioned in the main text:
“MUA activity was smoothed with a Gaussian kernel (σ=4ms) and demeaned to obtain a continuous spike density signal that was then analyzed to the same as CSD signals.”
The gammacomponent of the derived smoothed MUA signals was weak compared to CSD signals, which also resulted in much lower PLV estimates. We feel that it is not informative to see raw MUA signals, because it is difficult to visually see the gamma cycles among the noise at single trial level.
Figure 7—figure supplement 4: Does "gamma amplitude" mean the absolute amplitude of gamma band after Hilbert, or the amplitude of real signal? Also, are those amplitudes reflecting the sums of gamma in 2 sites? Please clarify.
The gamma amplitude is the absolute value of the Hilbert transform of the two gamma signals in a pair. In the figure we plotted the mean gamma amplitude as a function of phase difference (populationaveraged). We clarified this in the figure legend of Figure 7—figure supplement 4 as follows: “CSD gamma amplitude modulation from monkey M1 averaged over the two contacts for a given pair with modulation strength expressed in percent change from the mean (yaxis) and as a function of phase difference (xaxis).”
I may be wrong, but every panel in the figure shows full cycle of amplitude with a peak near zero phase difference, that makes me think it is real amplitude.
On the Yaxes of all the panels, we show gammaγ amplitude modulation (% of the mean) as a function of delta phase. This is shown for M1, M2 and he PING model. The different subplots represent different detuning and interaction strength levels.
Is it worth plotting gamma amplitude against interaction strength, ε?
We believe the relationship between gamma amplitude modulation and interaction strength is well illustrated in Figure 7—figure supplement 3.
In those plots, peak position is slightly rightward. Does it relate to phase leading/lagging?
The amplitude modulations were indeed not centered at zero phasedifference. It is likely to relate to detuning and the established preferred phaserelation.
Similar properties were observed for PING network simulations, which convinced us that these properties probably tell something about how networks interact. However, we believe that this should be investigated by its own systematically in a future study.
Same question again as above: If IF1IF2 and IF2IF1 are counted equally, then should the curves be symmetric around Δphase = 0?
The detuning values were not absolute (not rectified), so we selected values from one side of the Arnold tongue. Therefore, e.g. for averaging of pairs with detuning of 4Hz, we averaged over conditions where for IF1IF2 or IF2IF1 the second contact had frequency 4Hz higher than the first one. That is the reason the curves do not need to be symmetric
Supplementary Information – Correction for CSDinduced phase shifts: "we normalized the phasedifferences for each given contact pair to the condition having the smallest frequency difference," by doing what? Is this parallel translations of points along axes?
Yes, it is a parallel translation. We subtract the phase difference from the condition of smallest frequency difference with the phase differences from all other conditions. It shifts all points along the axis and does not change any other distribution properties. This was clarified in the Appendix, as follows:
“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).”
In the final paragraph of the subsection: "in the case of mutually anatomically coupled cortical locations, detuning strongly influences the main direction of information flow" As presented, this is questionable.
We agree with the reviewer that the sentence should be written more carefully. We have done so and now provide also a bit more context, which now read as follows:
“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”.
We argue that for mutually coupled cortical location, like between V1 locations, the phaselag is dependent on detuning. Further, as shown for the granger analysis, but also for network simulations published by us and others (Besserve et al., 2015; Buehlmann and Deco, 2010; Cannon et al., 2014; Lowet et al., 2016a), it will shape the information flow between the networks as measured by information theoretical tools (e.g. transfer entropy). But we agree that much more work is needed and these results are only indicative.
While grand G(θ) was used to predict PLV, preferred phase difference could differ among pairs according to Figure 6. Were G(θ) values shifted along the axis of phase difference to align the downward peak of the curve to the preferred phase difference for each pair?
To average the grand G we selected conditions with sufficient detuning such that they had similar phase differences (was written in the Supplementary Materials) to avoid smearing. It is now mentioned in the text:
“Note that the function G(θ) was estimated from data with absolute detuning of more than 4Hz. This was done based on the observation that interaction functions became deformed when detuning was close to (see supplementary materials for further details). Further it avoided smearing due to phase shifts occurring mainly within ± 4Hz.”
Subsection “The Arnold tongue and the regulative parameters of gamma synchronization”. Can authors show the range of detuning? Were there cases that can show large deviations in detuning that break down interactions?
The range of detuning values can be seen from Figure 6 and Figure 7.
In Figure 7 conditions of high detuning (Δω=+6Hz) and low interaction strength (ε=1Hz) led to PLV values close to zero (strong desynchronization).
We added two descriptive sentences:
“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.”
In the same subsection: while larger detuning values would decrease opportunities for synchronization across distance in V1, when they are tied strongly to stimulus information in a cell's RF, e.g., the Gabor stimuli used in Ray and Maunsell, there are still opportunities for the neurons in V1 (seeing low contrast) to synchronize with ones in V2 with overlapping RF (also seeing low contrast) and the same for V1/V2 interactions across the different contrast regions of the Gabor. This means basically frequency segregation of channels seeing different "pixels."
This is a statement we agree on. In fact, we are working on a modeling paper that deals with implications for V1 and V2 interactions and the shape of V2 RF.
In our Plos comp bio paper (Lowet et al., 2015), we argue that detuning can help in local feature integration/ figureground segregation. Image areas of low heterogeneity ‘bind’ because of low detuning and image regions of high heterogeneity decouple because of high detuning (e.g. around object borders). This we describe in, as quoted here:
“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.”
We believe that the selforganization of gamma rhythmic activity will have implication for how V2 neurons integrate V1 information and also the relation between V2 neural populations: “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.”
Future studies will need to investigate better how V1V2 gammaband interactions operate (in both directions) and the implications for visual computation.
Supplementary information: Equation 18 in Supplementary Information duplicates Equation 17. I guess the label "18" was meant for the line: "V(θ)"?
Thanks for pointing it out.
Reviewer #2:
[…] 1) The authors seem to confuse the PRC and the coupling function throughout the paper. They are not the same thing. In the theory of weak coupling, there are three things:
1) the PRC or adjoint function, Z(t) which is the infinitesimal PRC and is hard to compute with real neurons, but east to compute with models (see for example the tutorial: http://www.math.pitt.edu/~bard/bardware/meth3/node15.html#SECTION00052000000000000000on weak coupling of a simple voltage gated model;
We thank the reviewer for providing the link.
2) The interaction function, H(φ) which is defined as:
H(θ_{2}θ_{1}) = (1/T) ∫0TZ(s+θ1). C_{12}(u(s+θ_{1}),u(s+θ_{2})) ds
That is, it is the convolution of the PRC (Z) with the coupling between the oscillators, here denoted by u(t) with 1 = postsynaptic and 2 = presynaptic
3) The odd part of the interaction function which in this paper is G(θ) arises as follows
θ_{1}' = ω_{1} + ε_{12} H(θ_{2} – θ_{1}) Equation 1
θ_{2}' = ω_{2} + ε_{21} H(θ_{1} – θ_{2})
If you let φ = θ_{2} – θ_{1}, the phase difference and if ε_{12}=ε_{21}, then
φ' = Δω + εG(φ) Equation 2
where G(φ) = H(φ) – H(φ) is proportional to the odd part of H(φ)
Equation 2 is what the authors measure; it is not the PRC. I just want to make that clear.
We agree and thank the reviewer for his comments. We indeed did not take into account a coupling function in the equations that capture synaptic dynamics. This implies that we cannot be sure to have measured the underlying PRC.
We corrected this and clarified the relationship between the measured interaction function G to the PRC and the coupling function. Given that the synaptic dynamics that underlie gamma rhythms are relatively fast, we think that the (oddpart) PRC will be likely similar to the measured interaction function G. We clarified this in our text:
“It is worth noting that the interaction function G(θ) is not identical with 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.”
Since they measure φ, they can now compute the actual interaction function at both sites and this is even more informative than is G. By measuring θ_{1}' and θ_{2}', and plotting this against φ, equation 1 shows that
θ_{1}' = ω_{1} + ε_{12} H(φ)
θ_{2}' = ω_{2} + ε_{21} H(φ)
thus, in addition to G, they will be able to measure ω_{1} + ε_{12} H(φ) and ε_{21} H(φ). This will give them more information about the direction of coupling and also will give them the full value of H. Indeed, there is nothing here that requires the two H's to even be the same. Thus, I would suggest plotting θ_{1}', θ_{2}' vs. φ as well as φ' vs. φ.
We thank the author for pointing to this important issue. In this manuscript we have focused on the bidirectional interaction function G, being the central function of the equation used, for predicting experimental synchronization data.
Prior to submission, we have tried to obtain (trustable) unidirectional interaction functions H for experimental data. However, the difficulty is to estimate the intrinsic frequency ω_{1} or ω_{2} which are hidden variables. Without a good estimate of the intrinsic frequencies, it cannot be determined which part of the single H is negative or positive. Future work is required to resolve these issues.
As we show in supplementary materials we can however obtain an estimate of the general shape of single H and their amplitudes for V1 data (see examples in Figure 3—figure supplement 1)
We estimated the single H of PING networks, which was highly asymmetric (mainly positive). This was possible because for PING networks the intrinsic frequencies ω_{1} and ω_{2} can be obtained by simply decoupling the networks.
The detuning variable is also a hidden variable. Taking the assumption that both rhythms have similar interaction strengths, we found that the detuning is the mean of the measured interaction function. This is also true even if the underlying H is asymmetric. We tested the validity of this approach using phaseoscillator simulations as well as PING simulations where the true detuning value was known.
3) G(φ) is plotted looks very sinusoidal. Is this a consequence of the way the phase is extracted? I suggest taking a model pair of oscillators, such as Hodgkin huxley and coupling them with synapse along with noise and heterogeneity in coupling and in frequency. Then use the method that is used for the data to extract G(φ) and H(φ) and then compare this to the actual determistic values found by setting the noise to zero and computing H, G as in the above. This would assure me that the method of using Hilbert transforms isn’t somehow removing all the Fourier modes beyond the lowest.
To show that the obtained G is not trivial, we simulated phase oscillators with different types of true G. As we show in Figure 5—figure supplement 1 the resultant G reflected well the true G for different SNR.
We need however to state that the noise in our measurements might have smoothed our interaction function estimates. Also that we applied CSD might have led to further variability and averaging out of more transient components. We state this problem more explicitly now:
“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) γband interaction functions.”
4). Since the frequency varies rapidly between the two sites, does this variability covary between the two populations? I ask this because one can take two uncoupled oscillators and apply a broadband correlated signal to them and they will produce a very nice peak in their phasedifference histogram. For example, see Nakao, Arai, Kuramoto PRE 2005, or other papers by Arai & Nakao, or for a neuroscience version, Zhou, Burton, et al. 2013 Frontiers in Comp Neuro.
It concerns the question whether synchronization is due to direct mutual connections between local cortical locations or due to coordination by common input modulation.
Testing the idea directly would require measuring the feedforward synaptic inputs originating from the retinageniculate pathway and see whether input correlation strength among V1 locations does influence gammaband synchronization. This unfortunately goes beyond the scope of this study.
However, we think that shared broadband input as explanation for synchronization patterns found here is unlikely for the following reasons:
Empirically (Bosman et al., 2009; Brunet et al., 2013; Lowet et al., 2016b; Ray and Maunsell, 2010; Roberts et al., 2013) gamma frequency fluctuations in V1 and between V1 and V2 have been shown to have positive correlation, however with a strong dominance in the lower frequency bands (deltatheta 15Hz band). Low frequency modulations appear unlikely to be good candidate for phase synchronization of gamma rhythms.
A plausible scenario is that two cortical locations share directly some common input. This would imply that neurons from the locations have partially overlapping receptive fields (RFs). However, cortical V1 locations with clearly distinct RF positions still showed substantial phase coordination.
Previous studies have suggested that gammaband synchronization in V1 has a spatial scale that fits well the known spatial spread of horizontal connectivity (Eckhorn et al., 2001; Palanca and DeAngelis, 2005). This was confirmed by our own data. This fits also well with observations that gammaband properties depends on surround receptive field activation (Gieselmann and Thiele, 2008). Further, recent work in mice have shown that V1 gammaband synchronization depends on particular interneurons (SOM+) that are important for horizontal cortical interactions (Veit et al., 2017).
Furthermore, frequency correlation per se can be cause as well as consequence of synchronization. For example, PING networks with mutual interactions with unshared input variability will exhibit increased frequency correlation with stronger connectivity. That is the reason is very important to measure the correlation at the input level to properly test the idea.
We wrote on this important topic in the Discussion:
“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). 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).”
Reviewer #3:
[…] The results are impressive, but I would like to see a control analysis where the full analysis pipeline is applied to nonsimultaneously recorded (and so unsynchronised) trials (either by shuffling trials or combining nonsynchronously recorded probes). This would ensure that the SSD step is not picking up some common effects (e.g. stimulus onset) that is inducing the deltaphasedeltaif relationship e.g. in Figure 2E. (I don't think full permutation inference is necessary for all results, but just some confirmation that indeed the actual experimental data pipeline does produce flat deltaphase / deltaif curves)
First, we apologize that we forgot to mention the exact time window used for our analysis. Our trials consisted of 1sec baseline fixation and 2sec visual stimulation. For our gammaband analysis we focused on the stimulation period 0.2sec2sec. This was fixed in the text, which read as follows:
“For analysis the stimulation period (0.22sec) was used, not including the first 200ms to avoid stimulusevoke transients.”
The authors propose that the method of fitting TWCO model to data and associated analyses such as Arnold tongue plots might be applicable to other brain regions or frequency band relationships. It would be nice if they could give some heuristic guidance for experimental designs to which these techniques can be applied. Presumably they need to systematically modulate both connection strength and frequency difference, as done here?
Thanks for this suggestion. We wrote a paragraph in the Discussion section about how detuning and interaction strength could be manipulated in future studies:
“Nevertheless, future studies are required to test to what extent weakly coupled oscillator principles apply to different frequency bands across brain regions. […] 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.”
The statistical analysis and Figure 6CD show the overall average over all electrode pairs and stimulus conditions. Presumably the experimental contrast already induces a range of PLVs for each specific recording pair, so perhaps they could be evaluated statistically as well (e.g. scatter with pair PLV R2 as a function of distance between the pairs, or something along similar lines).
We computed for each single contact pair separately how much variability over the 9 contrast conditions could be explained by the model predictions. The variation for each single contact pair represent variation by contrast conditions /detuning. The results are PLV R^{2} of M1: 0.27 ± 0.0002, n=802 and M2: 0.1 ± 0.0001, n=882. For phase differences we found that R^{2} of M1: 0.52 ± 0.0001, n=802 and M2: 0.44 ± 0.0004, n=882. The R2 values can differ from the main analysis, because the variability over single contacts (representing differences in cortical distance/interaction strength) is missing. Notice also that the values need to be interpreted in the context of noisy single contact pair estimation and that the model was not directly fitted to the data. We included the analysis in the text for PLV and for mean phase difference
The authors indicate that data and code are available on request. I would strongly urge them then to deposit the data in a suitable repository such as Dryad to reduce their future administrative burden and to ensure long term stable access.https://opennessinitiative.org/
We agree with the reviewers that data sharing is an important issue in neuroscience and beyond. We will upload (after the acceptance of the manuscript, on dryad as suggested) the singletrial CSD data of each monkey allowing the reader to assess the reported observations and to test new ones. This way the final version of the manuscript will contain a link to the project files. Full access to the (raw) data is given upon reasonable request, given the size and complexity of the data. To enhance other researchers’ understanding of the observed phenomena and TWCO, we also added Matlab codes for phaseoscillator simulations and PING network simulations. Further, we will add codes for Singular spectrum decomposition (SSD).
https://doi.org/10.7554/eLife.26642.026Article 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).
Reviewing Editor
 Charles E Schroeder, Columbia University College of Physicians and Surgeons, United States
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.
Metrics

 2,583
 Page views

 463
 Downloads

 42
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Neuroscience
Deciphering patterns of connectivity between neurons in the brain is a critical step toward understanding brain function. Imagingbased neuroanatomical tracing identifies areatoarea or sparse neurontoneuron connectivity patterns, but with limited throughput. Barcodebased connectomics maps large numbers of singleneuron projections, but remains a challenge for jointly analyzing singlecell transcriptomics. Here, we established a rAAV2retro barcodebased multiplexed tracing method that simultaneously characterizes the projectome and transcriptome at the single neuron level. We uncovered dedicated and collateral projection patterns of ventromedial prefrontal cortex (vmPFC) neurons to five downstream targets and found that projectiondefined vmPFC neurons are molecularly heterogeneous. We identified transcriptional signatures of projectionspecific vmPFC neurons, and verified Pou3f1 as a marker gene enriched in neurons projecting to the lateral hypothalamus, denoting a distinct subset with collateral projections to both dorsomedial striatum and lateral hypothalamus. In summary, we have developed a new multiplexed technique whose paired connectome and gene expression data can help reveal organizational principles that form neural circuits and process information.

 Neuroscience
Blindness affects millions of people around the world. A promising solution to restoring a form of vision for some individuals are cortical visual prostheses, which bypass part of the impaired visual pathway by converting camera input to electrical stimulation of the visual system. The artificially induced visual percept (a pattern of localized light flashes, or ‘phosphenes’) has limited resolution, and a great portion of the field’s research is devoted to optimizing the efficacy, efficiency, and practical usefulness of the encoding of visual information. A commonly exploited method is noninvasive functional evaluation in sighted subjects or with computational models by using simulated prosthetic vision (SPV) pipelines. An important challenge in this approach is to balance enhanced perceptual realism, biologically plausibility, and realtime performance in the simulation of cortical prosthetic vision. We present a biologically plausible, PyTorchbased phosphene simulator that can run in realtime and uses differentiable operations to allow for gradientbased computational optimization of phosphene encoding models. The simulator integrates a wide range of clinical results with neurophysiological evidence in humans and nonhuman primates. The pipeline includes a model of the retinotopic organization and cortical magnification of the visual cortex. Moreover, the quantitative effects of stimulation parameters and temporal dynamics on phosphene characteristics are incorporated. Our results demonstrate the simulator’s suitability for both computational applications such as endtoend deep learningbased prosthetic vision optimization as well as behavioral experiments. The modular and opensource software provides a flexible simulation framework for computational, clinical, and behavioral neuroscientists working on visual neuroprosthetics.