QUANTIFICATION AND STATISTICAL ANALYSIS

Information about the relevant statistical test can be found in the corresponding results section with additional information concerning data preprocessing and selection in the following Data analysis section. Information about sample variables and size is indicated in the results section, and information pertaining to figure panels can be found in the corresponding figure legend. In general, we applied non-parametric statistical tests, thereby avoiding assumptions about the distributions of our empirical data.

All data analysis was performed using custom code and the Fieldtrip toolbox (Oostenveld et al., 2011), both written in MATLAB (MathWorks).

MUA rate was smoothed with a Gaussian (for constant light stimulation: SD = 12.5 ms; for stimulation with pulse trains and sinusoids: SD = 1.25 ms; in each case truncated at ± 2 SD) to obtain the spike density.

To quantify the locking of neuronal responses to optogenetic stimulation, we calculated the Pearson correlation coefficient between MUA spike density and laser intensity as a function of time shift between them.

LFP power spectra were calculated for data epochs that were adjusted for each frequency to have a length of 4 cycles and moved over the data in a sliding-window fashion in 1 ms steps. Each epoch was multiplied with a Hann taper, Fourier transformed, squared and divided by the window length to obtain power density per frequency. For the different stimulation frequencies f, LFP power is shown as ratio of power during stimulation versus pre-stimulation baseline (−0.5 s to −0.2 s relative to stimulation onset).

MUA-LFP locking was quantified by calculating the MUA-LFP PPC (pairwise phase consistency), a metric that is not biased by trial number, spike count or spike rate(Vinck et al., 2010). Spike and LFP recordings were always taken from different electrodes. For each spike, the surrounding LFP was Hann tapered and Fourier transformed. Per spike and frequency, this gave the MUA-LFP phase, which should be similar across spikes, if they are locked to the LFP. This phase similarity is quantified by the PPC as the average phase difference across all possible pairs of spikes. For a given MUA channel, MUA-LFP PPC was calculated relative to all LFPs from different electrodes and then averaged.

Sinusoidal stimulation of different frequencies enabled estimation of neuronal response latencies. This is highly relevant when optogenetic stimulation is used to produce temporal activation patterns at high frequencies. In addition, it validates that the responses we observe are a result of optogenetic stimulation: Neuronal response latencies to optogenetic stimulation are typically on the order of 3–8 ms; By contrast, shorter latency responses are likely to reflect photo-electric artifacts (Cardin et al., 2010). To investigate response latencies, we averaged MUA responses aligned to the peaks of the sinusoids (Figures S5B–S5F). During sinusoidal stimulation, the light was modulated between the respective maximal intensity and nearly zero intensity. Thus, the light crossed the threshold for effective neuronal stimulation at an unknown intensity, and it is not possible to calculate response latencies in the same way as has been done for pulse trains. Therefore, we used a technique of latency estimation that has been developed in the study of synchronized oscillations, and that is based on the slope of the spectrum of the relative phase between two signals (Schoffelen et al., 2005), in our case the light intensity and the MUA. Figure 3C shows this relative-phase spectrum and reveals a strictly linear relationship between relative phase and frequency. A linear frequency-phase relation is a signature of a fixed time lag, because a given time lag translates into increasing phase lags for increasing frequencies (Schoffelen et al., 2005). The slope of this linear relationship allowed us to infer a latency of 5.5 ms, in good agreement with previous reports of neuronal latencies.

The GC spectrum was first estimated separately for each recording site and subsequently averaged over sites. For each trial, we estimated the Fourier transforms of the input (laser) and the output (MUA). Specifically, each trial was segmented into non-overlapping epochs of 500 ms length. Per epoch, the time series of the input and the output were multiplied with a Hann taper, they were zero-padded to a length of 1000 ms, and their Fourier transforms (FTs) were obtained. The FTs were used to calculate the power-spectral densities (PSDs) of the input and of the output, and the cross-spectral density (CSD) between input and output. CSDs and PSDs were averaged over trials and used for the estimation of GC by means of non-parametric spectral matrix factorization (Dhamala et al., 2008). For the example GC spectrum (Figure 4C), the error region was determined by a bootstrap procedure, with 100 iterations, each time randomly choosing 30% of the trials. The shown error boundary is the region containing 95% of the bootstrapped estimates. For the average GC spectrum (Figure 4E), the error region indicates the standard error of the mean across the recording sites.

All inferences were based on the combined data of all animals, for which a given experiment was performed. The resulting inferences are limited to the studied sample of animals, as in most neurophysiological in-vivo studies.

High-resolution spectra of LFP power changes and MUA-LFP PPC were compared between stimulation with blue light and control stimulation with yellow light (Figures 1I and ​and1J).1J). We calculated paired t tests between spectra obtained with blue and yellow light, across recording sites. Statistical inference was not based directly on the t tests (and therefore corresponding assumptions will not limit our inference), but the resulting t-values were merely used as a well-normalized difference metric for the subsequent cluster-based non-parametric permutation test. For each of 10,000 permutations, we did the following: 1) We made a random decision per recording site to either exchange the spectrum obtained with blue light and the spectrum obtained with yellow light or not; 2) We performed the t test; 3) Clusters of adjacent frequencies with significant t-values (p < 0.05) were detected, and t-values were summed over all frequencies in the cluster to form the cluster-level test statistic. 4) The maximum and the minimum cluster-level statistic were placed into maximum and minimum randomization distributions, respectively. For the observed data, clusters were derived as for the randomized data. Observed clusters were considered significant if they fell below the 2.5^th percentile of the minimum randomization distribution or above the 97.5^th percentile of the maximum randomization distribution (Maris and Oostenveld, 2007). This corresponds to a two-sided test with correction for the multiple comparisons performed across frequencies (Nichols and Holmes, 2002).

The neurons in the PING model are Hodgkin-Huxley-like point neurons. The excitatory population consists of a simplified version of model pyramidal neurons introduced by (Traub et al., 1991), the reduced Traub-Miles (RTM). The inhibitory population consists of model basket cells introduced by (Wang and Buzsáki, 1996). The parameters for the model are presented in the tables below, and we refer to the original publication of the model for more details (Börgers, 2017).

PING Neuron parameters:

PING Network parameters:

In order to reproduce the experimentally observed hysteresis effects, we implemented spike frequency adaption in the model pyramidal neurons. The PING+M model is taken from the Adaptation-based, Deterministic Weak PING model from Börgers (Chapter 32 of (Börgers, 2017)). In this model, the previous PING model is modified by the addition of a model M-Current to the pyramidal neurons. Otherwise, the network is identical to the PING model described above.

PING+M Neuron parameters:

In order to investigate the generality of the model results, we next implemented a simple network of leaky-integrate-and-fire neurons. This network was composed of 80% excitatory neurons and 20% inhibitory neurons, coupled via instantaneous synapses. Excitatory neurons were not mutually connected, while the remaining connectivity was all-to-all, with synapse magnitude randomly distributed uniformly between 0 and the respective post-synaptic-potential (PSP) value. Each neuron accumulates postsynaptic potentials until the threshold for spiking is reached. Upon spiking, each neuron transmits to its synaptic partners a post synaptic event and its potential is reset. The membrane voltage of the model LIF neurons is given by: dV/dt=−V/C+I/C, with the membrane timescale tau = R*C, where R is the input resistance of the neuron, C is the membrane capacitance, and I includes both basal and synaptic currents. We drove the network with symmetric single slow sine waves or with white noise. The dynamics of the network were evaluated numerically at a resolution of tau using the Euler method.

LIF Network parameters: