Distinct Structure of Cortical Population Activity on Fast and Infraslow Timescales

Abstract Cortical activity is organized across multiple spatial and temporal scales. Most research on the dynamics of neuronal spiking is concerned with timescales of 1 ms–1 s, and little is known about spiking dynamics on timescales of tens of seconds and minutes. Here, we used frequency domain analyses to study the structure of individual neurons’ spiking activity and its coupling to local population rate and to arousal level across 0.01–100 Hz frequency range. In mouse medial prefrontal cortex, the spiking dynamics of individual neurons could be quantitatively captured by a combination of interspike interval and firing rate power spectrum distributions. The relative strength of coherence with local population often differed across timescales: a neuron strongly coupled to population rate on fast timescales could be weakly coupled on slow timescales, and vice versa. On slow but not fast timescales, a substantial proportion of neurons showed firing anticorrelated with the population. Infraslow firing rate changes were largely determined by arousal rather than by local factors, which could explain the timescale dependence of individual neurons’ population coupling strength. These observations demonstrate how neurons simultaneously partake in fast local dynamics, and slow brain-wide dynamics, extending our understanding of infraslow cortical activity beyond the mesoscale resolution of fMRI.

A single action potential lasts about a millisecond, and a second suffices for a vast range of sensory-motor and cognitive behaviors, such as recognizing pictures and sounds, getting up or sitting down, or recalling a memory. Accordingly, most neurophysiological research has focused on subsecond timescales. However, several neural processes occur over much longer timescales (Huk et al. 2018). Transitions between sleep and wakefulness and between different stages of sleep occur on timescales of minutes and hours (Weber and Dan 2016;Lecci et al. 2017;Meisel et al. 2017). During wakefulness, changes in arousal can span tens of seconds and minutes, yet they affect performance in subsecond behavioral tasks (Harris and Thiele 2011;Palva and Palva 2012;McGinley et al. 2015). Slow timescale dynamics has been revealed by resting-state fMRI (Raichle 2015), which infers neural activity from the (slow) changes in blood supply to different areas of the brain. However, fMRI monitoring of neural activity is limited to the so called infraslow range of 0.01-1 Hz. Furthermore, both fMRI and other approaches to study mesoscale infraslow cortical dynamics-such as electro-and magneto-encephalography (EEG, ECoG, LFP, MEG, e.g., see Popa et al. 2009;Palva and Palva 2012;Mitra et al. 2018), and intrinsic and voltage-sensitive fluorescent-protein or dye imaging in experimental animals (White et al. 2011;Chan et al. 2015;Kraft et al. 2017)-cannot characterise individual neurons' relationships to infraslow activity.
The relationship of individual neurons to infraslow brain dynamics, and the relationship between a neuron's coupling to infraslow and fast mesoscale cortical dynamics, is thus poorly understood. For example, how much can the firing rate of an individual neuron change over tens of seconds and minutes, and how can these slow dynamics be summarized quantitatively? To what extent are slow changes in firing rate correlated among neurons, and what is the structure of these slow correlations? To what extent are neurons' relationships to slow and fast firing rate fluctuations similar, and might they be driven by the same underlying mechanisms?
Here we addressed these questions by analyzing multihour recordings of neuronal populations in mouse medial prefrontal cortex (mPFC), performed using chronically implanted high-density silicon probes. We found that neuronal spike trains have f 1/ power spectral density (PSD), and that PSD in combination with interspike interval (ISI) distribution suffices for an accurate quantitative model of single neuron spiking dynamics on both fast and slow timescales. Coupling between individual neurons and the population rate (defined as the summed rate of all the spikes recorded by the probe) was timescale-dependent, with many neurons strongly coupled to population rate on fast timescales but weakly coupled on slow timescales, or vice versa. Furthermore, on slow but not fast timescales, neurons' phase preference with respect to the population rate was bimodal. Finally, in frequencies ≤0.1 Hz population rate was highly correlated with arousal as reflected by the pupil area. These results suggest that dynamics on fast and infraslow timescales are distinct processes, and likely regulated by distinct mechanisms at the single neuron level.

Electrophysiological Recordings
All experimental procedures were conducted according to the UK Animals (Scientific Procedures) Act 1986 (Amendment Regulations 2012). Experiments were performed at University College London under personal and project licenses released by the Home Office following institutional ethics review. Adult C57BL/6 mice of both sexes were used.
The experimental procedures for chronically implanting Neuronexus and Neuropixels probes were previously described in (Okun et al. 2016;Jun et al. 2017). Briefly, in an initial surgery under isoflurane anesthesia animals were implanted with a custom built head-plate. Following full recovery and acclimatization to head-fixation, probe implantation was performed under isoflurane anesthesia. The probes were implanted through a craniectomy above mPFC (0.5 mm lateral and 1.8 mm anterior to bregma). Neuronexus probes (A2 × 2-tet with CM16LP connector and Buzsaki32 with CM32 connector) were lowered 1.7 mm, placing the recording sites in the prelimbic cortex (PrL). Neuropixels probes were lowered~3.5 mm (so that the most superficial of the 374 recording sites remained outside of the brain, while the deepest sites were~3.5 mm inside the brain; the recording sites were thus placed in the cingulate, prelimbic, and infralimbic cortices). The probes were oriented approximately parallel to the cortical layers,~0.5 mm lateral offset of the insertion point relative to the midline implied that the probes resided in cortical layers 5 and 6 (which was also confirmed histologically).
Recordings were performed over the course of several months following the probe implantation. For head-fixed recording, mice were placed inside a plastic tube where they could comfortably sit or stand. The recordings lasted 1.5-3 h. In animals implanted with a Neuronexus probe, recordings were performed using OpenEphys (www.open-ephys.org) recording system . Mice with a Neuropixels probe were recorded using SpikeGLX system (github.com/billkarsh/SpikeGLX) developed at Janelia Farm. (Some of the mice were trained and recorded in a behavioral task which the animals would perform for water reward (Lak et al. 2018); the data analyzed here is from recordings of ongoing activity in separate sessions without behavior, performed on days when the animals were not water deprived).
Recordings in freely behaving animals implanted with a Neuronexus probe lasted 4-8 h. Mice were briefly head-fixed to allow attaching the amplifier head-stage to the probe and then released into their home cage, where they were free to engage in any activity of their choice, while being monitored to make sure that the thin cable leading from the amplifier to the OpenEphys box was not entangled.
We have evaluated the quality of spike sorted units using isolation distance metric (Schmitzer-Torbert et al. 2005) and by quantifying the contamination of the refractory periods of the spike autocorrelograms, which was expressed as proportion of the number of spikes in the first 2 ms of the autocorrelogram relative to the autocorrelogram asymptote (Harris et al. 2000). We have limited the analysis to units with isolation distance > 20 and refractory period contamination < 0.2. Our analyses yielded quantitatively similar results when more (and less) stringent criteria were applied.
A possible concern is that our results, instead of reflecting the properties of actual infraslow fluctuations in the firing rates of cortical neuronal populations, are dominated by contamination introduced by unstable recordings. Such concern is not unique to the present work, and was raised in the past regarding estimation of pairwise correlations (Ecker et al. 2010). Here, Neuropixels recordings provided an unprecedented opportunity to detect and monitor drifts, as the recording sites span a contiguous stretch of >3 mm. For spikes detected simultaneously on several contacts, we have computed the vertical location of the "center of mass" of the spike, according to the relative amplitude of the spike waveform on each contact. Changes in these locations over time, particularly for high-amplitude spikes, reveal potential drifts. In the example shown in Fig. S1, multiple drift events are visually apparent. In each event, the vertical locations of high-amplitude spikes at a particular neighborhood of the probe drift~10 μm upwards over the course of~5 s, and over the next~20-40 s return to their original locations. Similar drift pattern occurs~200 μm further down the shank (Fig. S1b,c), which is a strong indication that these drifts are produced by a vertical movement of the probe with respect to cortical tissue, rather than any other cause. In fact, drifts were simultaneously observed at >10 locations across the top 800 μm of the probe. Drifts were not observed when vertical location of low-amplitude spikes was considered. This is consistent with the idea that a vertical movement has a much larger impact on the waveforms of high-amplitude spikes originating from neurons abutting the probe, compared with the low-amplitude spike waveforms of neurons that are more (horizontally) distant.
In the example recording, drifts of high-amplitude spikes were not observed on contacts deeper than~800 μm (Fig. S1d). In recordings from this mouse, all data originating from the top 1 mm was excluded from the analyses. No similar drifts were observed in Neuropixels recordings from the second animal. In recordings performed with Neuronexus probes, the recording sites were located only at the bottom 200 μm of probes which were lowered 1.7 mm into the brain, thus to the extent that drifts of Neuropixels and Neuronexus probes are similar, we do not expect to find vertical drifts in these recordings (since the recording sites were not covering a contiguous interval, the above drift analysis cannot be repeated for Neuronexus recordings).
An additional observation suggesting that our results are not driven by drifts concerns the relationship between amplitude of the different units and their phase preference. If drifts introduce a strong bias into our estimation of phase with respect to population rate (see below), then there might exist some consistent relationship between the phase and spike waveform amplitude of the different units, because drift bias is expected to be stronger for units close to the probe and having high-amplitude spike waveforms. However, no significant correlation between phase and spike amplitude was found in our data. We conclude that drift is an important caveat that has the potential to bias measurements of spiking activity obtained with extracellular probes, however in view of the control analyses explained above, we believe that the phenomena described here are not due to such drifts.

Pupil Tracking
Pupil area was tracked as previously described in Burgess et al. (2017). Briefly, a camera (DMK 21BU04.H or DMK 23U618, The Imaging Source) with a zoom lens (ThorLabs MVL7000) was focused on one of the eyes of the animal. The eye was illuminated by an infrared LED (SLS-0208A, Mightex). Videos of the eye were acquired at ≥30 Hz. In each video frame, excluding frames with blinks, an ellipse was fit to the pupil image, and pupil area was estimated based on this fit.

Single Spike Train Analysis
PSD of individual spike trains (Fig. 2a) was estimated using mtspectrumsegpb function of Chronux toolbox (chronux.org). Note that this function is specifically intended for point processes such as spike trains (Mitra and Bokil 2007). PSD in different frequencies was estimated by breaking the entire recording into segments of appropriate length, and averaging across them. Specifically, segments of length l were used to estimate PSD in l 7/ to l 10/ frequencies (for example, a 2 h recording is broken into 72 000 segments, 0.1 s each, and these segments are used to estimate PSD in 75-100 Hz frequencies. Then the recording is broken into 54 014 segments, 0.133 s each, to estimate PSD in 56-75 Hz frequencies, and so on. Eventually 12 segments of 556 s are used to estimate PSD in 0.013-0.018 Hz frequencies). For frequencies <0.01 Hz the spectrum of the entire recording was computed using mtspectrumpb function without breaking it into segments. For presentation purposes only it was further smoothed using Matlab's smooth function.
PSD of a spike train is closely related to, but distinct from, the PSD of the underlying (continuous) firing rate intensity. If the firing rate λ ( ) t is itself a random process with power spec- , then the power spectrum of the spike train where μ λ is mean firing rate (Lepage et al. 2011). To intuitively see why the μ λ term in the right hand side of the equation is required, consider the simplest case where λ ( )= t const, that is, the case of a constant intensity (homogeneous) Poisson spike train stochastic process. In this case (ω) = λλ S 0 , and (ω) S nn = μ λ . A homogeneous Poisson spike train has constant power in all frequencies for the same reason that white noise has power in all frequencies. The power spectrum is the Fourier transform of the autocorrelation (a result known as the Wiener-Khinchin theorem). The autocorrelation of a Poisson spike train is a delta function at time 0, because of the correlation of each spike with itself, and is zero at all other times since no other spikes are correlated. The Fourier transform of a delta function is a constant function, explaining why the power spectrum of a Poisson spike train is flat.

Single Spike Train Modeling
The spike train model which captures both the fast and slow timescale dynamics of cortical spiking relies on both ISI distribution and PSD of spike trains (Fig. 2). The goal of the model is to generate synthetic spike trains satisfying constraints on both distributions simultaneously.
Let I denote the observed ISI distribution of a spike train. For modeling, I was represented by a histogram with logarithmically spaced bins of all the observed ISIs (32 bins were used to describe ISIs, from 1 ms up to 200 s). Instead of using the PSD of the spike train itself, the model uses the PSD of the underlying continuous firing rate intensity, which we denote by ℘ (the two are closely related but distinct, as described above). For modeling, ℘ was represented by the power of a continuous signal obtained by convolving the observed spike train with a 50 ms FWHM Gaussian (50 parameters were used to represent ℘). With I and ℘ as its inputs (82 parameters in total), the goal of the model is to generate a synthetic spike train whose ISI distribution and PSD are as close as possible to the original spike train. An intermediate step towards this final goal is constructing a continuous firing rate intensity signal ( ) r t for the synthetic spike train. However, we start by constructing a different firing rate intensity signal, ( ) r t 1 , by sampling ISIs from I and convolving the resulting spike train with a 50 ms FWHM Gaussian. Typically ( ) r t 1 will have much less power in the infraslow frequencies than what is required. On the other hand, the straightforward way to generate a signal with power ℘ (i.e., to use inverse Fourier transform) produces a signal whose values are normally distributed with 0 mean, which is inappropriate for a firing rate intensity function. Therefore, we used an iterative algorithm of (Schreiber and Schmitz 1996) to generate a signal ( ) r t with power ℘ and distribution of values of ( ) r t 1 . Once ( ) r t was generated, we sample a spike train n 1 using ( ) r t as a time dependent firing rate signal. In the final step the ISIs of the spike train are adjusted to have the desired distribution I. Specifically, we convert the sequence of ISIs in n 1 into a sequence of ISI ranks, by replacing each ISI with its rank among all the ISIs of n 1 . We build the final output spike train n by sampling from I the same number of ISIs found in n 1 and rearranging them according to the sequence of ranks from n 1 , that is, the ISI rank sequences of n and of n 1 match.
When the model is used to generate an output without an explicit constraint on I (i.e., only ℘ input is provided), it implicitly assumes I has an exponential distribution, with an additional constraint forbidding ISIs <2 ms (representing a hard refractory period).
Although it is possible in principle to consider a spike train model which relies solely on the autocorrelation structure (equivalently the PSD), these models are complex and not necessary for our purposes. The reason we used ISIs to model fast timescale dynamics and PSD for slow timescale dynamics is that ISI distribution provides a simple generative model for fast timescale dynamics of spike trains, whereas methods for constructing spike trains with a given (fast timescale) autocorrelation are substantially more complex (Krumin and Shoham 2009;Macke et al. 2009).

Time Domain Population Coupling on Fast and Slow Timescales
Time domain correlation between spike trains of single neurons and the population rate ( Fig. 3d-f) was computed as previously described in (Okun et al. 2015). Specifically, we computed the inner product between the vectors representing the population rate and single unit spike train at different lags (using Matlab's xcorr), and normalized it by the number of spikes of the single neuron. For fast timescale correlation, the vectors were at 1 ms resolution, and the single neuron spike train was smoothed with Gaussian of 12 ms halfwidth. For slow timescale correlation, the vectors were at 1 s resolution. In both cases the baseline (average values 800-1000 ms away from zero lag for fast timescale correlation, and average values 12-20 s away from zero lag for slow timescale correlation) was subtracted.

Coherence Analysis
For analyzing the relationship between spike trains of individual units and population rate, the latter was obtained by summing all the spikes detected on all the shanks/tetrodes barring the one on which the single unit was recorded. For Neuropixels recordings, where the entire probe consists of one shank (with 374 recording sites over~3.5 mm) this approach was not applicable. Instead, for each unit we have performed our analyses with population rate based on all spikes on the probe (except for those of the unit itself) and with population rate based only on spikes from recording sites >60 μm away from the location of the single unit. All results were almost identical for both conditions. The population rate typically was >100 spikes/s. Coherence between population rate or pupil area and individual units was estimated, together with its confidence interval, using coherencysegpb function of the Chronux toolbox (estimating coherency using coherencysegcpb, where population rate was considered a continuous signal rather than spike count, produced identical results). We used theoretical, asymptotic confidence intervals as computed by Chronux (we have also found that it provides more stringent, i.e., wider, intervals than jackknife). As in the case of PSD estimation, coherence in different frequencies was estimated by breaking the entire recording into segments of appropriate length for each frequency.
There is a non-intuitive correspondence between the coherence of some continuous signal with a spike train and its coherence with the continuous spike intensity signal underlying the spike train (e.g., the neuron's membrane potential). Unlike the more familiar case of coherence between a pair of continuous processes, coherence between a continuous process and a point process (such as a spike train of a neuron) depends on the PSD and the rate of the latter. This (mathematically unavoidable) fact has two implications. First, the f 1/ profile of firing rate PSD implies that coherence of the spike train with population rate falls with frequency, even when coherence between the underlying firing rate intensity and the population rate does not (Fig. S2c,d). Second, because coherence depends on the rate of the spike train, 2 neurons whose firing rate intensities are exactly proportional but unequal do not have the same coherence with population rate. To account for this second issue of rate dependence, we use a correction factor to produce a coherence which reflects a firing rate of 1 spike/s, rather than the actual firing rate of the neuron.
More formally, let ( ) y t be a continuous process, ( ) n t a spike train of a single neuron and λ ( ) t the underlying firing rate, that is, ( ) n t is a doubly stochastic Poisson process with intensity λ ( ) t . It holds that where, C yn and λ C y denote the coherence between ( ) y t and ( ) n t or λ ( ) t , μ is the mean rate of ( ) n t , and λλ S is the power spectrum of λ ( ) t (Aoi et al. 2015). From the above equation, it is clear that 2 spike trains with proportional but unequal rate intensities where α ≠ 1) will have different values for coherence. This is not desirable, therefore instead of reporting the coherence between a spike train and the population rate, we report "rate adjusted coherence" which reflects the coherence that would have been measured if the neuron had a firing rate of 1 spike/s, that is, if its firing intensity was λ( ) μ t / instead of λ( ) t , see Fig. S2a,b for an example. We use a correction factor of ( + (μ − )μ (ω)) − 1 1 / S nn 1/2 , as derived in Aoi et al. (2015), to obtain the rate adjusted coherence. The PSD of the spike train, used for the correction was estimated as described above. The rate adjusted coherence still depends on the PSD of the spike train of the single unit. For instance, it is possible to have 2 intensity functions λ ( ) t 1 and λ ( ) t 2 with equal means of 1 spike/ s, and equal coherence with ( ) y t , but with different power spectra. In this case, Equation (1) implies that if n 1 and n 2 are spike trains with intensities λ 1 and λ 2 , then . Here, we did not attempt to remove this dependence, which would have required an accurate estimate of (ω) λλ S . In practice, λλ S cannot be directly inferred from S nn because the assumption that ( ) n t is a doubly stochastic Poisson process with intensity λ( ) t does not hold. For example, the existence of refractory period in ( ) n t reduces the power in all low frequencies in (ω) S nn (Bair et al. 1994;Rivlin-Etzion et al. 2006). Furthermore, the spiking of actual neurons is driven by changes in the subthreshold membrane potential ( ) V t m which is rather distinct from λ( ) t , as exemplified in Fig. S2. Of note, this discussion primarily applies to high frequencies, whereas in low frequencies μ is significantly lower than (ω) λλ S (or (ω) S nn ) and thus the discrepancy between (ω) C yn and (ω) λ C y is minor (see Equation 1).
To compare the strength of population coupling of 2 simultaneously recorded neurons across all timescales, we have compared their rate adjusted coherences in the following specific frequencies: 0.01, 0.03, 0.1, 0.32, 1, 3.2, 10, 32, and 100 Hz. If the null hypothesis that first neuron has higher rate adjusted coherence in these frequencies could be rejected at P ≤ 0.001 (after using Bonferroni correction for performing 9 comparisons), and the reverse null hypothesis could also be rejected with P ≤ 0.001, the 2 were considered as (a positive) example of a simultaneously recorded pair of neurons where neither neuron dominated the other across all frequencies.

Phase Analysis
Phase of spiking of single units with respect to population rate or pupil area was estimated using the same Chronux toolbox functions used to estimate the coherence (see above). As in the case of PSD estimation, phase in different frequencies was estimated by breaking the entire recording into segments of appropriate length for each frequency. After the phase in each segment was estimated, circular mean and standard deviation were computed. If the distribution of phases (across the segments) had no statistically significant (at P ≤ 0.05) mean, no phase was assigned (e.g., the nonsignificant neurons in the histogram in Fig. 4c,d).
All phases are specified with respect to the population rate, for example, a phase of π − /4 means that the single unit lags behind the population rate, whereas a phase of π /6 means that the unit leads it.
Our MATLAB code for the estimation of PSD, coherence and phase, following the above described procedures, is publicly available on github.com/m-okun/FrequencyDomainPopulationAnalysis.

Linearity and Logarithmicity Indices
), consider the following expression: We define the linearity index of ( ) f x as the value of this expression for = ( + ) k a b /2. Similarly, the logarithmicity index is the expression's value for = k ab . The rationale for these definitions is that for a function changing on a linear scale, total variation in the first and second halves of [ ] a b , is expected to be of comparable magnitude. Thus, linearity index is close to 0 for functions changing on linear scale (the function itself does not have to be linear, e.g., ( ) x sin on any sufficiently long interval), positive for supralinear functions, and negative for sublinear functions. On the other hand, for a function changing on a logarithmic scale, the total variation in [ ] a ab , and [ ] ab b , intervals is expected to be of comparable magnitude, thus its logarithmicity index would be close to 0 (while its linearity index would be negative).
For empirically measured ( ) f x , total variation is contaminated by measurement noise. To avoid this problem, and relying on the fact that phase functions were either monotonic or had just a few extremum points (typical examples shown in Fig. 5), we used the following expression instead of the one given above: ). For each neuron, equation (2) was evaluated using the longest continuous interval of frequencies over which phase was well-defined (i.e., it had no frequencies in which coherency with population rate was not statistically significant). Neurons for which such interval spanned less than an order of magnitude were excluded from the analysis.

Results
To examine the intrinsic spiking dynamics of single cortical neurons on timescales extending to tens of seconds and minutes we used multisite silicon probes to record the activity of neuronal populations in the frontal cortex. All the recordings were performed using chronically implanted probes (Okun et al. 2016;Jun et al. 2017). We used 16-and 32-channel Neuronexus probes in 6 animals and 374-channel Neuropixels probes in 2 additional animals. The recordings lasted 1.5-3 h in head-fixed mice, standing or sitting in a plastic tube (22/26 recordings, 730/775 neurons), and 4-8 h in freely behaving mice residing in their home cage (4/26 recordings, 45/775 neurons).

Single Neurons Show Dynamics at Multiple Timescales
On fast timescales the spiking dynamics of cortical neurons can be summarized by the ISI distribution. The characteristic irregular firing of cortical neurons results in ISIs varying by several orders of magnitude (Softky and Koch 1993), with some neurons also exhibiting ISI histogram peaks indicating rhythmicity at particular frequencies (Fig. 1a). A neuron's ISI distribution was not on its own sufficient to account for the structure of its spike train at long timescales. Indeed, visual inspection shows that cortical firing rates typically fluctuate on timescales of minutes or more (Fig. 2b), longer than almost all ISIs of neurons with firing rates ≥1 spike/ s. Synthetic spike trains created by randomly reshuffling the original ISIs did not have this slow timescale dynamics (Fig. 2c). The discrepancy between actual and ISI-shuffled data could be summarized by the spike count Fano factor: the variance divided by the mean of spike counts in bins of prescribed temporal duration (Fig. 2d). For bins of short duration (1-100 ms), Fano factors were close to 1, and the Fano factors of the original and shuffled data were similar. However, for bins of 1-100 s, the Fano factors of actual data were several-fold higher. Across all analyzed neurons (n = 775), the Fano factor of the number of spikes in 1024 ms bins was 1.6 times higher in the actual data compared with ISI-shuffled trains, and with 16 384 ms bins it was 4.8 times higher (for a summary across all bin sizes see the error of ISI model in Fig. 1e, below).
Although the ISI distribution could not alone capture the infraslow porion of a cell's spiking dynamics, the combination of ISI distribution and spike train PSD provided a good approximation. Because our recordings lasted several hours, we were able to compute power spectra down to very low frequencies, where the PSD values were much higher than for fast frequencies, indicative of infraslow dynamics (Fig. 1a). We devised an algorithm that generates synthetic spike trains with prespecified PSD and ISI distributions (see Materials and Methods). The slow-timescale firing dynamics of these synthetic spike trains was visually similar to the original data (compare Fig. 1b with Fig. 2b) and closely matched the observed Fano factors over multiple timescales (Fig. 1c), as expected from the analytical relationship between Fano factor and autocorrelation of a stationary spike train ). The full model was better than models that used either PSD or ISI distribution independently (Fig. 1d). At slow timescales (e.g., 1 min; Fig. 1e), Fano factors predicted from PSD alone are significantly closer to the actual values than ISI-based predictions, but still not as accurate as the full model. For fast timescales (e.g., 100 ms), ISIs predict spike count accurately, but the PSD alone is insufficient (Fig. 1e). The full model respects both constraints, and as a result provides predictions that are significantly better than either ISI distribution or PSD alone (e.g., for 1-min bins its average error is 47%, compared with 1072% and 305% for ISI only and PSD only models, Fig. 1d,e). These results also demonstrate the major contribution of slow dynamics to spiking variability in the cortex.
Cortical neurons are diverse in their intrinsic dynamics. This diversity is well characterized at short timescales by differences in spike regularity (Maimon and Assad 2009) and by the differing propensity of neurons to emit complex-spike bursts (McCormick et al. 1985;de Kock and Sakmann 2008), but dynamical diversity at slow timescales is largely unexplored. To address this question, we observed that the PSD of most neurons had power-law profile over the 0.01-1 Hz range, with the exponent significantly different between neurons (Fig. S3a). Fitting spike train power with a β const/f function over 0.01-1 Hz revealed that the power-law exponent β covered a range of 0.39 ± 0.19 (mean and standard deviation for n = 775 neurons), and was conserved when fit separately for the first and second halves of each recording ( Fig. 1f; R 2 = 0.69 overall; median R 2 of individual recordings = 0.63; P <0.05 in each recording). The power law exponent was unrelated to mean firing rate and weakly related to bursting (Fig. S3b-d). We conclude that cortical neurons are diverse in the strength of their infraslow firing rate fluctuations, and that the structure of these fluctuations can be summarized, to first approximation, by the PSD slope β.

Population Coupling Strength is Unrelated at Fast and Slow Timescales
To understand how the slow dynamics of individual neurons was related to that of the entire population, we started by considering how individual neurons are related to the population rate-the summed rate of all spikes detected on the probe. At short timescales, neurons vary continuously in the strength of their coupling with population rate (Okun et al. 2015). To characterise the relation of neurons to the population across multiple timescales, we extended this analysis into the frequency domain. Analysis in frequency domain does not suffer from the inherent ambiguity of time domain analysis, where correlations computed using a time bin of a specific duration reflect co-modulation not just on the scale of the bin, but also on all slower timescales (Brody 1999).
The PSD of population rate had 1/f profile, similar to the profile of PSDs of single neuron spike trains. In high frequencies the value of the PSD of any spike train is dominated by the firing rate term (see the section on single spike train analysis in Materials and Methods), therefore in high frequencies PSD of population rate was close to the sum of PSDs of all the individual spike trains that together comprise the population rate. However, in all frequencies <1 Hz, the population rate PSD was several fold higher than the sum of PSDs of the individual spike trains (Fig. 3a, Fig. S4). As the PSDs of independent neurons would add linearly, this result indicates that infraslow fluctuations in firing rates of neurons were correlated in all these frequencies. To estimate the coherence between population rate and spike trains of specific neurons, we considered the former as a continuous function of time, and computed its "rate adjusted coherence" (Aoi et al. 2015) with the spike train of each neuron, which accounts for differences in mean firing rate between neurons (Fig. S2a,b; see Materials and Methods). To verify that this method provided a reliable measure, we estimated coherence separately in two halves of single recordings, obtaining similar estimates for coherence on both slow and fast timescales for most neurons (for 0.1 and 10 Hz, correspondingly: overall R 2 = 0.63 and 0.75, median R 2 of individual recordings = 0.59 and 0.61, P <0.05 in 19/26 and 22/26 recordings; Fig.  3b,c, see also Fig. S7).
Coherence analysis revealed widely diverse relationships to population rate, both between neurons, and between timescales within individual neurons. In all cases, coherence decayed to zero with increasing frequency. This behavior is a by-product of the point processes nature of spike trains and need not signify a decay in coherence between population rate and the membrane potential of individual neurons (Fig. S2c,d; see the section on coherence analysis in Materials and Methods). Importantly, the manner of this decay varied greatly between neurons (Fig. 3d-f). In some cases the relative strength of different neurons' population coupling was conserved across frequencies (e.g., the red neuron in Fig. 3d has consistently larger coherence than the green or blue neurons). However, simultaneously recorded neurons often showed different rates of coherence decay: in 25% of simultaneously recorded pairs each neuron had a significantly stronger coherence in a subset of frequencies (Fig. 3e,f). Furthermore, some neurons' coherence with population rate was nonmonotonic (15% of cells, e.g., green neuron in Fig. 3f). On average across neurons, rate adjusted coherence with population rate remained <0.5 in all our recordings, even in frequencies as low as 0.1-0.01 Hz (Fig. 3b,  Fig. S5).

Phase of Population Coupling Differs Across Timescales
Coherence is an indication of a constant phase relationship between 2 processes. Thus, if a neuron has high coherence with population rate at a given frequency, this means it fires at a reliable phase with respect to the population-but does not imply that this phase is zero. Phase analysis showed that most neurons had a stable phase preference with respect to population rate across halves of the recording on both slow and fast timescales (Fig. 4a,b, see also Fig. S7). It also revealed a major difference between phases of population coupling on slow and fast timescales, with out of phase activity several-fold more likely in the infraslow range (Fig. 4a-e). Specifically, at ≥10 Hz just 5% of cells had phase closer to π than to 0, whereas at ≤ 0.3 Hz this was the case for ≥ 28% of the cells (P < 10 −28 , Z-test for equality of 2 proportions). This however did not completely summarize a neuron's phase preference: even within a mode, there remained significant correlation in a neuron's precise phase from one half of the recording to the other (Fig. 4a), and at high frequencies phases also had reliable nonzero values across the 2 halves of the dataset (Fig. 4b,d). The preferred phase distribution in the infraslow frequencies was not symmetric: more neurons had phases between π /4 and π 3 /4 (i.e., leading the population rate) than between π −3 /4 and π − /4 (lagging the population rate, e.g., 16% vs. 10% at 0.1 Hz, P <0.001). The fact that neurons show an asymmetric phase distribution relative to their summed activity might seem contradictory, but was possible because neurons which lagged the population had higher firing rates compared with neurons which led it (4.5 ± 4.8 vs. 3.1 ± 4.8 spikes/s, for π −3 /4 to π − /4 vs. π /4 to π 3 /4, P <0.0001). In contrast to the phase preference of individual cells, the relative phase of population rate on different shanks or tetrodes (or on different segments of the Neuropixels probe) was close to 0 in all frequencies (Fig. S6), suggesting that variations in population coupling phase mainly differ within, rather than between local populations.
The phase spectra of individual neurons were diverse, and could show a complex dependence on frequency. According to the phase distribution histogram (Fig. 4e) one expects to find neurons' phase preference to be close to either 0 or π at infraslow frequencies, and close to 0 at high frequencies. While many neurons conformed to this pattern (Fig. 5a, Fig. S7a), neurons that were anticorrelated with infraslow population rate differed in the dependence of phase preference on frequency: it was discontinuous, with clear subdomains and a drop in coherence in some neurons, but gradual in others (Fig. 5b,c, Fig. S7b,  c). We also observed neurons whose phase preference did not fit the overall pattern, for example, having phase preference of π/2 in infraslow frequencies (Fig. 5d, Fig. S7d) or exhibiting altogether different behaviors (Fig. 5e, Fig. S7e).
Phase modulation seemed to occur on logarithmic rather than linear scale (Fig. 5c-e). To assess the rate of phase changes, we devised indices which quantify the linearity and logarithmicity of the phases (see Materials and Methods). The linearity index is 0 for phase changing on linear scale, and it is positive (negative) for phase changing on supra-(sub-) linear scale. Similarly, the logarithmicity index is 0 for phase changing on logarithmic scale, and it is positive (negative) for phase changing on supra-(sub-) logarithmic scale. We found that phase spectra were overwhelmingly changing sublinearly (just 4% had positive linearity index), whereas the logarithmicity index values were about equally distributed around 0 (logarithmicity index of 59% of the neurons was positive, Fig. 5f). We conclude that the phase between single neuron and population rates predominantly changes on a logarithmic scale with frequency.
The logarithmic rate of phase preference change implies that phase in nearby frequencies is similar. In other words, when only a small range of frequencies is considered (e.g., on linear scale), the phase is approximately constant, and thus to a first approximation single cell neuronal dynamics with respect to population rate is independent of timescale. To test this prediction we compared how well constant phase and linear phase models fit the phase preference in nearby frequencies (0.1 Hz vs. 0.03 Hz or 0.32 Hz). The former model corresponds to timescale-independent dynamics, the latter to a lead or lag by a fixed time interval between an individual neuron and the population rate. As predicted by the logarithmic rate of phase change, the constant phase model fit the data substantially better than the linear model (23% explained variance vs. no explained variance for 0.32 Hz, and 28% vs. 17% explained variance for 0.03 Hz, Fig. S8).

Infraslow Dynamics Correlates With Pupil Diameter
Head-fixed mice, such as those we recorded here, show fluctuations in alertness levels over time. To address the degree to which the infraslow dynamics we observed could relate to alertness, we monitored the animals' pupil area in a subset of head-fixed recordings (Fig. 6a). At 0.03 Hz, 65% (350/541) of the neurons were significantly coherent with the pupil area signal, and the magnitude of this coherence was consistent when estimated from separate halves of the recording ( Fig. 6b; P <0.01 in 10/13 recordings, the median percentage of variance in one half of the data explained by the other half across recordings: 41%). Phase preferences were similarly stable, and the phase distribution had 2 clear peaks~π rad apart (Fig. 6c), consistent with the existence of 2 populations positively and negatively coupled to arousal (Stringer et al. 2018).
Next, we considered how individual neurons' coupling to the pupil and to the local population rate are related. Visual inspection of population rate and the pupil area signals suggested the two are similar in the infraslow range (Fig. 6a), which was confirmed by coherency analysis showing that the two signals were highly coherent in frequencies ≤0.1 Hz (Fig.  6d). Correspondingly, individual neurons' coherence with pupil area closely matched their coherence with population rate (Fig.  6e). Phases with respect to the population rate and the pupil area were also closely matched; a consistent gap between the two (which at 0.03 Hz constituted 0.78 ± 0.51 rad) indicates that neuronal spiking preceded the pupil signal (Fig. 6f).
Coupling to pupil fluctuations was related to the infraslow dynamics of firing of individual neurons. We observed that the power-law exponent β obtained from fitting the spike train power with a β const/f function (Fig. 1f) correlated with pupil coherence across the recorded neurons (r = 0.43, P <10 −9 ). Interestingly, this correlation was highly significant only for neurons firing in phase with the pupil  with the rate adjusted coherence with pupil (at 0.03 Hz), for neurons whose phase with respect to pupil fluctuation was close to 0 (the left peak in the bimodal histogram in f, shown in orange; r = 0.64, P < 10 −9 ). (h) No significant correlation was found between the PSD power-law slope and pupil coherence for neurons whose phase with respect to the pupil was close to π (the right peak in the bimodal histogram in f, shown in purple; P = 0.66). (i, j) Coherence (phase) of individual neurons with pupil signal on slow timescale (0.03 Hz) and their coherence (phase) with population rate on fast timescale (10 Hz) are uncorrelated (P = 0.85, Spearman correlation in i, P = 0.31, circular correlation in j). (k) Fast timescale phase was uncorrelated with coherence with the pupil (P = 0.77). (l) A significant correlation between phase with pupil at 0.03 Hz and coherence with population rate at 10 Hz (P < 0.001) was observed, where neurons anticorrelated with the pupil were more coherent with population rate on fast timescales. fluctuations, and insignificant for the antiphase neurons (Fig. 6g,h).
Population coupling on fast and slow timescales were largely unrelated: we observed no relationship between coherence with population rate on fast timescales and coherence with the pupil signal (P = 0.85, Fig. 6i; P >0.2 in each individual recording, Spearman correlation) and no relationship between the phases (P = 0.31, Fig. 6j, P >0.17 in each individual recording, circular correlation). Similarly, no significant relationship between fast timescale phase and coherence with the pupil was observed (P = 0.77, Fig. 6k), although phase with pupil did show a statistically significant relation to fast timescale coherence with the population (Fig. 6l, this might be due to some subclasses of cortical neurons differing in their fast and slow timescale population coupling properties). Considered together, these observations are consistent with the idea that population coupling on slow timescales is controlled by separate mechanisms from the local synaptic inputs driving fast timescale population coupling.

Discussion
We used frequency domain analysis to examine the activity of neuronal populations in mPFC across frequencies spanning 4 orders of magnitude (0.01-100 Hz). Our findings point to a fundamental difference between fast and infraslow timescale cortical dynamics. The strengths of a neuron's population coupling at fast and slow timescales were unrelated; furthermore, at fast timescales nearly all neurons fired at preferred phases close to 0 relative to population rate, whereas at slow timescales the phase distribution was bimodal, with preferred phases of~30% of the neurons closer to π . Population coupling in infraslow, but not fast frequencies reflected coupling to brain-wide arousal signal (pupil area). While these general rules held for most neurons, a great diversity of fine-detailed behaviors was seen within local populations, for example regarding the slow-timescale dynamics as captured by a neuron's power spectrum, and the way its coherency with population rate depended on frequency.
The difference between population coupling at fast and slow timescales likely indicates different mechanisms driving these types of coupling. Fast timescale dynamics reflects local synaptic activity (Haider and McCormick 2009), and the fasttimescale population coupling of individual neurons is correlated with the number of the local synaptic connections they receive (Okun et al. 2015). In contrast, infraslow dynamics correlates with global, brain-wide phenomena related to arousal, which are controlled at least in part by neuromodulatory inputs (McGinley et al. 2015;Reimer et al. 2016); a similar mechanism has been suggested for the global component of resting state fMRI measurements (Scholvinck et al. 2010;Wong et al. 2013;Turchi et al. 2018). The fact that a neuron's population coupling on fast and slow timescales were uncorrelated therefore suggests that the degree to which a neuron's firing is controlled by global brain states is unrelated to its local connectivity; for instance, a neuron weakly affected by neuromodulatory tone could still be strongly wired into the local network. The hypothesis that fast and slow population coupling arise through different mechanisms is supported by observations of neurons whose phase with population rate had discontinuous subdomains in high and low frequencies, as one would expect to see if the slow but not fast fluctuations are produced by a mechanism that suppresses the firing of these neurons while elevating the population rate (Fig. S9). This hypothesis is also supported by the fact that only slow-timescale population coupling phases were bimodal, an observation consistent with previous work in mouse primary visual cortex (V1). Specifically, in a previous study in V1 we observed neurons with weak fasttimescale population coupling but only very few with negative fast-timescale coupling (Okun et al. 2015). Furthermore, finding two populations of neurons in V1 that couple negatively as well as positively to arousal (on slow-timescales) was reported by Vinck et al. (2015) and Stringer et al. (2018). Whereas the recording conditions in these works were not identical to ours (mice were recorded on a wheel, and in  arousal was increased by a sensory stimulus), the combined evidence from these studies seems to suggest that the dynamical properties we have described are not unique to the deep layers of mPFC. On the other hand, based on analysis of Utah array recordings in primates, it was suggested that the structure of neuronal activity in cortical networks is preserved across timescales (Kiani et al. 2015); it is unclear if this discrepancy is due to difference in species, cortical areas, or other factors.
Most of our present day knowledge on infraslow cortical dynamics comes from fMRI studies of resting state activity (Buckner et al. 2013;Raichle 2015;Foster et al. 2016). Because fMRI provides a blood-oxygen-level dependent (BOLD) signal rather than a direct measure of neural activity, it is limited to measurements on infraslow timescales. Multiple studies have shown that the BOLD signal correlates with population rate (Logothetis et al. 2001;Ma et al. 2016;Mateo et al. 2017), although disagreements on the BOLD signal's interpretation remain (e.g., see Winder et al., 2017). Our study provides an account of how individual neurons' activities combine to produce infraslow fluctuations in population rate, and hence in BOLD (to the extent the two are correlated). The low-frequency power of the population rate was 2-5 times larger than it would be if cells were independent of each other (Fig. 3a, Fig. S4). Because the recorded populations were spread over hundreds of micrometers, this increase would likely have been even higher if the recorded populations were concentrated in a smaller volume. We observed that the contribution of single neuron activity to the mesoscale signal is limited in two ways. First, for majority of neurons their coherence with population rate remained relatively low (typically between 0.2 and 0.4) even in the 0.01-0.1 Hz range of frequencies (Fig. S5), and for some neurons coherence in this range was found to be even lower than for higher frequencies (Fig. 3f). Second, the infraslow fluctuations in firing rate of many neurons were partially or completely out of phase with the population (Fig. 4). The present study is limited to measurements of spiking activity in one cortical area (including in cases where several population rate signals from different shanks or parts of the silicon probe were compared), the coherency relationship between neurons and population rates of distal cortical areas remains to be elucidated, for example, see Mitra et al. (2018). Another potential caveat with linking the present work to resting state fMRI studies is the degree to which the activity we observed is in the pure resting state regime (Logothetis et al. 2009;Winder et al. 2017). While it is possible to find short intervals during which a mouse does not move, this is not the case for intervals longer than a few seconds (e.g., typically mice move their eyes every 5-10 s). Thus our results should not be viewed as describing pure spontaneous activity (an ideal which is impossible to achieve in practice for infraslow timescales in awake subjects), but pertain to the actual cortical dynamics which is partially driven by intrinsic behaviors.
At the single neuron level, power spectral analysis was consistent with scale-free dynamics in the infraslow frequency range (Fig. 1a, Fig. S3). Such dynamics is typical of neuronal activity on various spatial scales, from fMRI measurements to ion channels, e.g., (Lowen et al. , 1999He 2011), and elsewhere in biology, for example, organization of heart beats (Bassingthwaighte et al. 1994). Similar spectra have been reported for retinal and thalamic cells recorded in anesthetized cats Lowen et al. 2001), and more recently in resting humans (Nir et al. 2008), with a mean power-law exponent of 0.45, close to the 0.39 value observed here (Fig. 1f,  Fig. S3d). For an intuitive interpretation of this value consider that for a signal with power spectrum proportional to 1/f 0.4 , 66% of the total power slower than any chosen frequency ω is concentrated in frequencies ≤ω/2. As a result of these slow changes in firing rate of individual neurons, their spike count variance in minute bins was on average~10-fold higher than what fast-timescale spiking dynamics alone (i.e., the ISI model) would predict (Fig. 1e). The interpretation and causes of such scale-free dynamics are controversial. Here we find that in neurons positively (but not negatively) correlated with arousal the power-law exponent is highly correlated with the strength of pupil coupling (Fig. 6g,h), implicating brain-wide neuromodulatory mechanisms. Another suggested possibility is that on slow timescales it is inherent to recurrently connected networks (Chaudhuri et al. 2018). Yet another possibility is that scale-free behavior could be caused by single-cell intrinsic mechanisms such as firing rate adaptation (Marom 2010;Xu and Barak 2017), which can result in a frequency-independent lead of~0.2 rad of the output spiking over sinusoidal input currents with periods <1 Hz (Lundstrom et al. 2008;Pozzorini et al. 2013). These effects can build up across more complex networks: for example, when rat whiskers were stimulated by white noise on top of which sinusoidal modulation with 0.3-0.03 Hz frequency was added, barrel cortex neurons preceded the sinusoidal stimulus envelope by~0.8 rad on average, while thalamic neurons were leading by less than half as much (Lundstrom et al. 2010). An alternative, recently proposed possibility is that infraslow firing rate fluctuations are driven by slow changes in ion concentrations (Krishnan et al. 2018). While the contribution of each of these mechanisms remains to be elucidated, it is likely that their effect on cortical dynamics is particularly complex on intermediate timescales (~1 Hz) where they interact with fast-timescale local synaptic activity.