Spatial and temporal correlations in human cortex are inherently linked and predicted by functional hierarchy, vigilance state as well as antiepileptic drug load

The ability of neural circuits to integrate information over time and across different cortical areas is believed an essential ingredient for information processing in the brain. Temporal and spatial correlations in cortex dynamics have independently been shown to capture these integration properties in task-dependent ways. A fundamental question remains if temporal and spatial integration properties are linked and what internal and external factors shape these correlations. Previous research on spatio-temporal correlations has been limited in duration and coverage, thus providing only an incomplete picture of their interdependence and variability. Here, we use long-term invasive EEG data to comprehensively map temporal and spatial correlations according to cortical topography, vigilance state and drug dependence over extended periods of time. We show that temporal and spatial correlations in cortical networks are intimately linked, decline under antiepileptic drug action, and break down during slow-wave sleep. Further, we report temporal correlations in human electrophysiology signals to increase with the functional hierarchy in cortex. Systematic investigation of a neural network model suggests that these dynamical features may arise when dynamics are poised near a critical point. Our results provide mechanistic and functional links between specific measurable changes in the network dynamics relevant for characterizing the brain’s changing information processing capabilities.

It sounds like you have just randomly permuted the time series data. This will not maintain the power distribution as you claim. Reordering the time series points will change the power distribution, which you can easily verify by computing power spectrums before and after the random permutation. Also consider two extreme examples: ordering the points in descending order versus ordering the points based on every other point being descending from the max value or ascending from the min value. One will include only low frequencies, and one will include high frequencies. It is unclear whether your method is valid but the explanation is lacking or whether your chosen method has fundamental flaws.
Standard methods exist for creating surrogate data with a power spectrum that matches given data. On option is IAAFT, Schreiber & Schmitz (1996)  We thank the reviewer for their comment. We agree that if one shuffles the raw time series before applying the power spectral analysis the power distribution is not preserved. However, within the context of our manuscript, we want to kindly reiterate that this is not how the surrogate data were defined. Specifically, in our manuscript we first perform the power spectral analysis and extract the high gamma power (one value every 0.125 seconds), thus creating a new power time series of high gamma power. Only thereafter we shuffle this power time series of high gamma power in each 120 second segment. Thus, the distribution of high gamma power does not change as all values for high gamma power are preserved in this time series, but the timely order of these high gamma power values is lost and thus the temporal correlations are eliminated. This approach does not need any further assumptions or coefficient estimation and hence is very robust and computationally efficient. We apologize for this misunderstanding and have expanded on the explanation of how the surrogate data analysis performed to clarify our approach: "Specifically, after extracting the broadband γ-power time series, this power time series was then randomly shuffled within each 2-min segment. As the time series shuffling only happens after the power spectral density is computed the power distribution in the data is preserved but the temporal correlations in the data are eliminated." Thank you for helping us to clarify this more.

46
An essential ingredient for information processing is thought to be the ability of neural 47 circuits to integrate information over appropriate periods of time and across different 48 cortical areas. Characterization of neural systems according to their temporal and spatial 49 correlations (TC and SC, respectively) has consequently provided important insight into their 50 information processing capabilities. For example, in decision-making and working memory 51 tasks 1-4 , the ability to integrate over extended periods of time may increase the signal-to-52 noise ratio and afford to maintain some memory of past activity. In non-human primates, 53 temporal correlations have been found to increase along the functional hierarchy providing 54 a unifying principle for information integration across different timescales 1,5,6 . Similarly, the 55 ability to integrate information in space, across functionally specialized regions, is 56 considered essential for normal brain functioning 7-9 . Consequently, effective propagation 57 and integration of information in space has been shown to depend on vigilance state where 58 it is maximized during wake and breaks down during slow-wave sleep 10 . 59 60 While temporal and spatial correlations are thus promising in providing a principled 61 approach to characterizing information integration, essential questions related to them 62 remain only incompletely understood. First, temporal and spatial correlations have mostly 63 been studied independently. An understanding if and how temporal and spatial correlations 64 are related is still missing. Second, previous studies have largely been limited in duration, 65 focusing on certain tasks or windows of interest, thus providing only an incomplete picture 66 of the variability and fluctuations of these indices, their potential dependence on vigilance 67 state and interventions, such as medication treatment. Third, it is still unclear whether 68 concepts like the hierarchical ordering of temporal scales also apply to human brains. 69 Although there is strong evidence for a hierarchical ordering of temporal correlations in non-70 human primates 6 , evidence in humans has been limited to comparisons between core and 71 periphery networks in MEG 11 or fMRI 12 . Thus, whether the concept of hierarchical ordering 72 of temporal correlations, as measurable in iEEG, along processing pathways similarly applies 73 to human brains remains currently unresolved 13 . 74 75 Here we use long-term invasive EEG data to comprehensively map temporal and spatial 76 correlations according to cortical topography, vigilance state and drug dependence over 77 extended periods of time. We show that temporal and spatial correlations in cortical 78 networks are intimately linked, break down during slow-wave sleep (SWS) and decline under 79 antiepileptic drug (AED) application. Further, we report that temporal correlations in human 80 electrophysical signals follow a gradient which aligns with the functional hierarchy. Finally, 81 we study a neuron network model, that reproduces the tight interconnection between 82 spatial to temporal correlations, their hierarchical ordering and decline under AED action. 83 Our results provide novel mechanistic and functional links between specific measurable 84 changes in the network dynamics relevant for characterizing the brain's changing 85 information processing capabilities.

94
We analysed invasive electroencephalography (iEEG) recordings from 23 patients to 95 characterize spatial and temporal correlations (STC) as functions of vigilance state, 96 antiepileptic drug (AED) action, and of functional cortical hierarchy. Broadband -power (56-97 96 Hz) was used as an index of population firing rate near an electrode 14-18 to derive spatial 98 and temporal correlations and based on prior comparative work across broad frequency 99 ranges 62 . In line with these prior studies, we observed this high-frequency domain activity to 100 be best suited to resolve the correlations within cortical dynamics compared to other, lower 101 frequency bands (see Supplementary Figures S2 and S3 for TC and SC respectively).  102 Specifically, we measured the decay speed of the autocorrelation function of this signal as 103 an index for temporal correlations (TC; 5,19 ) and the decay of the pairwise cross-correlation 104 function over distance as an index of spatial correlations (SC; 20-22 ; Figure 1). 105 106 Spatial and temporal correlations are tightly interconnected 107 Previous research on STC as a metric for information integration in cortex has primarily 108 focused on task-related activity or investigation of relatively short periods, thus providing 109 only limited insight into fluctuations of STC over time, their co-variation, and dependence on 110 external factors, such as medications. We observed TC and SC to fluctuate significantly over 111 the course of a full day (Figure 2 A, B). Importantly, SC and TC were highly correlated, i.e., 112 when TC increased, SC tended to increase as well and vice versa. dosage normalized by the defined daily dosage and did not differentiate between different 146 AED due the high variability of AED given across patients (see Supplementary Table S2). Even 147 though different AED have different mechanisms of action, including ion-channel blockers, 148 e.g., Phenytoin, GABAergic drugs, e.g., Clobazam, which increase inhibition, previous work 149 has suggested that their influence on STC may be similar 19 . Comparison between full days 150 with highest and lowest AED load in each patient indicated both TC and SC to be higher 151 during the low AED days ( Temporal correlations increase with cortical functional hierarchy 166 Previous work in non-human primates has indicated that cortical areas exhibit a hierarchical 167 ordering in their timescales of their temporal correlations 6 . The sparse spatial sampling with 168 only a few electrodes per regions only allowed the evaluation of TC but not SC as a function 169 of functional hierarchy in our data. To determine the existence of a hierarchical ordering of 170 TC also in human cortex, we identified the five cortical regions along the visual pathway 171 which are similar in function to the regions used for the analysis in non-human primates in 6 . 172 The areas are colour-coded in Figure  days and nonSWS. We focused this analysis primarily on low AED load and nonSWS episodes 179 in an attempt to more closely resemble unperturbed, non-sleep related brain activity. TC 180 tended to increase from regions with low hierarchy to regions higher in hierarchy. The 181 average slope of increase was 0.4 ± 0.2 sec./area ( <0.05; Wilcoxon signed-rank test). No 182 increase with functional hierarchy was observed in the time-shuffled surrogate data ( Figure  183 3 D, grey bars). inhibitory synaptic strength by inh , for example via GABAergic drugs. As the results for both 213 mechanisms are qualitatively similar in the model, we only report them for the former (for a 214 direct comparison see also 19 ). 215 216 In the absence of AED action, collective dynamics exhibited a phase transition from low 217 activity to a high-activity phase when connection strength was increased. STC increased 218 when λ was increased to the critical value of λ=1 (  for the covariation of TC and SC could be given by slow changes of an external drive. 292 However, this might seem unlikely based on investigations over 48 hours and the systematic 293 changes observed under AED load or with cortical hierarchy. In particular the ability to tune 294 a system in a way that leads to a collapse of signatures of criticality has been argued to be a 295 strong indicator for critical dynamics emerging from network interactions in the unperturbed 296 state 45,46 . In our analysis, the increased AED loads can be regarded as tuning mechanisms of 297 the system, i.e., decreasing the network connectivity, which lead to decreased STC as 298 measures of the distance to criticality. 299 We here report to our knowledge for the first time in humans that temporal 300 correlations within highly resolved electrophysical signals increase along the functional 301 hierarchy. This is in agreement with observations in non-human primates using single-unit 302 activity 6 and results from a fMRI study in humans 12 . It is also in alignment with other recent 303 studies showing that more complex tasks, often associated with regions higher in the 304 functional hierarchy, might profit from longer temporal correlations 47,48 . Together, these 305 findings may thus indicate that long temporal correlations could be beneficial for 306 information integration especially in regions with higher functional hierarchy. The neural 307 network model studied here may exhibit such a gradient of temporal correlations by tuning 308 only one parameter in contrast to previous models explaining such gradients 1 . Specifically, 309 when tuning the model closer to criticality, dynamics exhibited longer spatio-temporal 310 correlations and thus longer memory within the signal 38. Extended memory (or also self-311 affinity) in a signal is also often quantified by the Hurst exponent which is 0.5 for 312 uncorrelated signals and grows up to 1 for correlated, self-affine signals. Our model closely 313 matched this expected behavior as it showed Hurst exponents close to 0.5 for simulations 314 far from criticality (λ≠1) and Hurst exponents around 1 at the critical value of the absolute 315 value of the largest eigenvalue (λ=1). Our model, however, consists only of one network, and 316 may thus only mirror one cortical region at a time. Therefore, it does not account for region-317 to-region interactions which, besides other effects like neural transmission speeds, might 318 play an important role leading to the correlation structure observed in the data. Hence, 319 further research is warranted to validate the mechanism leading to the spatial and temporal 320 correlation gradients we observed. 321 Long spatio-temporal correlations along with other capabilities 21,32,36,49-52 may 322 provide some benefits for information processing as information can be stored for longer 323 periods and inputs from more sites can be integrated. Conversely, when spatio-temporal 324 correlations decline, as under AED or during SWS, this may potentially impair information 325 processing. Impairments of cognition and perception are widely observed side effects under 326 AED treatment 53-57 and can be observed for neurological disorders as well 58,59 . Spatio-327 temporal correlations may therefore potentially serve as biomarkers for these deficits in 328 future research. 329 Even though we accounted for many confounding factors, including AED load, SWS, 330 seizures and subclinical events, removal of all channels with interictal epileptiform 331 discharges, it needs to be noted that the data still stem from patients with epilepsy which 332 may limit transfer of these concepts to healthy subjects. For instance, cortical dynamics may 333 potentially take multiple days to relax back to "normal" after surgery for electrode 334 implantation. As the AED load often correlates with the time after surgery this could 335 potentially be another, at least partial contributor to the observed changes in STC under 336 changing levels of AED. In our data, 17 out of 23 patients had their higher medication day 337 closer to the surgery date than the lower medication day. A sub-analysis of the six patients 338 with the low AED day closer to the surgery, however, did not show a significant opposite 339 effect in the other direction, which might have been an indication of the surgery effect. 340 Thus, while the current data does not suggest time to surgery to be a main driving factor, 341 further research is needed to better delineate the effects of AED load and relaxation of 342 cortical dynamics after surgery. Furthermore, these effects should also not affect changes of 343 STC during SWS and along the cortical hierarchy, as these were compared within patients 344 and did not have a timely delay to surgery. 345 Furthermore, the sleep staging used in this manuscript is only able to detect SWS. 346 Even though more elaborate algorithms which are also able to detect other sleep stages 347 from iEEG are in principle available they are not applicable to our data set due to the 348 sampling frequency 60 or data set duration 61 . While our research highlights the importance 349 to study temporal and spatial correlations over extended periods of time to capture their 350 variability, the sampling across cortical areas is naturally limited by electrode positioning 351 determined by the clinical need. This leads to a relative under-sampling for example of the 352 ACC in comparison to other regions. This sparse spatial sampling also limited us in calculating 353 spatial correlations for individual regions along the hierarchy. Empirical evidence indicative 354 of a strong co-variation between spatial and temporal correlations along with conceptual 355 insights from our model, however, strongly suggest that also spatial correlations should 356 follow a gradient along the functional hierarchy, similar to temporal correlations. The authors have no conflicts of interest to declare. 366 367 368

369
Preprocessing of invasive electroencephalography data 370 Multi-day invasive electroencephalographic (iEEG) recordings of 23 patients with epilepsy 371 undergoing presurgical evaluation at the Epilepsy Center of the University Hospital of 372 Freiburg, Germany, were analysed (12 female, age 28±13 years, mean ± standard deviation). 373 The data set was made available in the Epilepsiae database, their use for research was 374 approved by the ethics committee of the University of Freiburg and writteninformed 375 consent that the clinical data might be used and published for research purposes was given 376 by all patients 62 . Detailed data on the included patients can be found in the Supplementary 377 Table S2. 378 379 Data from two full days (midnight to midnight each), the day with highest and lowest anti-380 epileptic drug (AED) loads, were analysed in each patient. The day with highest and lowest 381 AED loads were determined by the summed prescribed drug dosages normalized by their 382 individual defined daily dosages across drugs. If more than one day had the same AED load, 383 highest and lowest AED days were chosen so that the time between days was maximized. 384 Data from each day was processed in consecutive 1 h segments. First, power line noise and 385 its first harmonic were removed using a notch filter at 50 and 100 Hz. Second, iEEG data, 386 which was sampled at either 256, 512 or 1024 Hz, was down-sampled to the common 387 frequency 256 Hz by applying an anti-aliasing filter and then decimating the signal. After 388 visual examination occasional 1 h data segments with additional frequency-restricted 389 artefacts were removed from further analysis. All epileptic seizures, including a 10-min 390 preictal and 10-min postictal period, were excluded for the purpose of this study. 391 Furthermore, all subclinical seizures, including 2 min before and after, as well as all channels 392 containing interictal discharges as labelled in the Epilepsiae database, were excluded from 393 the analysis 62 . 394 395 Assignment of electrodes to brain regions 396 The number of electrodes varied between patients and included both surface and depth 397 electrodes. Electrode placement was solely determined by clinical considerations. Individual 398 brain surfaces were constructed from MRI images and warped onto a common standard 399 MNI-152 template using the software package AFNI 63 . As a result, electrode positions were 400 assigned to a certain region if they fell within 9.5 mm of it (Euclidean distance), which is a 401 trade-off between positional accuracy and maximising the number of electrodes for each 402 region. Generally, the regions are far apart allowing for such a coarse assignment and an 403 observation of posterior to frontal changes of our measures. Electrodes were excluded when 404 their positions were not differentiable. Five human brain areas along the visual pathway 405 were chosen, functionally corresponding to the five regions investigated in non-human 406 primates in 6 . Hence, we adopted the names medial temporal (MT) and lateral intraparietal 407 (LIP) area in visual cortex, lateral prefrontal (LPFC), orbitofrontal (OFC) and anterior cingulate 408 cortex (ACC). The function associated with LIP in non-human primates can be more likely 409 associated with the medial part of the intra-parietal sulcus in humans 64 . Therefore, we 410 investigated this brain area but kept the name LIP for comparability to 6 . All regions were 411 drawn as set of regions constructed in 65 . In Figure 1 all regions except ACC can be seen. The 412 full list of the sets for each region can be found in the Supplementary Table S1. 413 414

Quantification of temporal correlations 415
Power fluctuations in the high -band have been shown to provide a local, spatiotemporal 416 estimate of population spike rate variations near an electrode 14-18 . We therefore evaluated 417 the fluctuations in the median power between 56-96 Hz and, for simplicity, refer to this 418 frequency band as broadband -power. Following previous work, for each iEEG channel, the 419 time series of broadband -power fluctuations were obtained by calculating the power every 420 125 ms (Welch's method, Hanning window; 5,19 . We observed the median broadband -421 power to be approximately log-normally distributed (see Supplementary Figure S1) and 422 therefore applied the logarithm with base 10 to the time series (see middle panel in Figure 1  423 for a typical time series). 424 425 Next, autocorrelation functions were calculated from consecutive, 75% overlapping 2-min 426 data segments of these individual-channel power time series (see top right panel in Figure  427 1). Following previous work, temporal correlations (TC) were defined as the time lag when 428 the autocorrelation function dropped for the first time below half of the difference between 429 the value at the first lag and the baseline (red dashed line in Figure 1) Additionally, note that the minimal value of the TC measure is one time lag, i.e., 0.125 435 seconds, which leads to a non-zero TC measure even in surrogate data (defined below). 436 437 Quantification of spatial correlations 438 Spatial correlations (SC) were calculated on the same broadband -power time series as 439 temporal correlations. Specifically, for all possible pairs of channels, the Pearson cross-440 correlation was calculated for each consecutive, non-overlapping 2-min segments. 441 Correlation values from channel pairs were averaged according to their Euclidean distance 442 with a bin size of 1 mm to get a distance dependent cross-correlation function. SC were then 443 defined as the average cross-correlation within the 7 to 79 mm interval (grey shaded area 444 within the cross-correlation panels in Figure 1). The interval was chosen as most of the 445 patients had electrodes within these distances (Figure 2 F). However, results did not change 446 when other distances were used, see Supplementary Figure 5. As this is a direct measure of 447 the cross-correlation function it is expected to be zero in surrogate (time-shuffled) data, as 448 can be seen in Figure 2F. 449  450  451  452  453  454  455  456  457  458  459  460  461  462  Surrogate Data  463 As additional controls, SC and TC were also calculated in the same way from randomly time-464 shuffled power time series (see lower panels in Figure 1). Specifically, after extracting the 465 broadband γ-power time series, this power time series was then randomly shuffled within 466 each 2-min segment. As the time series shuffling only happens after the power spectral 467 density is computed the power distribution in the data is preserved but the temporal 468 correlations in the data are eliminated. This was done for each electrode separately and, by 469 consequence, eliminated spatial correlations as well. Using these time-shuffled surrogate 470 data, TC and SC were calculated. Comparison of our data results with these surrogate data 471 allowed to rule out that our findings simply stem from fluctuations in the power but were 472 indeed caused by changes in the correlations in space and time. 473 474 475 Staging of vigilance states 476 To evaluate spatial and temporal correlations (STC) as a function of vigilance state, 477 specifically during and outside of slow-wave sleep (SWS), we calculated a vigilance index 478 following previous work 23 which is defined on 30-sec windows as the band power ratio 479 Θ + + high + spindle .
(1) 480 481 Segments were classified as SWS whenever the respective vigilance index was one standard 482 deviation above the mean vigilance index for the given full day 23 . Full 24-hour days were 483 scored individually for each patient to decrease the effects of possible multi-day rhythms 484 and potential impact under changing AED levels. A 2-min segment was finally classified as 485 SWS only if all 30-sec segments in it were classified as SWS and as nonSWS otherwise. 486 487 Neural network model 488 We studied STC as a function of network state in a parsimonious neural network model. In 489 the model = 1600 neurons were distributed on an equidistant 2-dimensional grid with 490 periodic boundary conditions. The weights of an all-to-all connected adjacency matrix were 491 first assigned random strengths from a uniform distribution which then were adjusted by a 492 Gaussian profile 493 − 2 2 2 , (2) 494 495 with being the Euclidian distance between neuron and neuron and scaling the width 496 of the profile similar to 66,67 . Self-connections were omitted and for all calculations shown 497 was set to = 4 as a trade-off between a fully connected and a not connected network as in 498 Ref. 67 . A total of = 20% of the neurons was randomly set to be inhibitory, i.e., their 499 weights are multiplied by −1 and the remaining 80% of neurons were excitatory, i.e., they 500 had positive outgoing connections. To model the influence of AED the weights of inhibitory 501 and excitatory neurons in the weight matrix can be separately scaled by a factor inh and 502 exc , respectively. The unperturbed network is given for inh = exc = 1 and increasing only 503 inh leads to stronger inhibition whereas decreasing exc leads to decreased excitable output 504 of the excitatory neurons. 505 Each neuron can be in either of two states ( ) ∈ {0,1} at any timepoint . Here, ( ) = 0 506 corresponds to the neuron not firing and ( ) = 1 to the neuron firing. The firing 507 probability in the next time step ( + 1) is determined by the sum over all inputs and 508 takes the form 509

719
Table S2 Meta data of the patients from the EPILEPSIAE database 3 . Drug loads are given in fractions of the 720 defined daily dosage, see Table S3 for full names of the drugs. Total number of electrodes are given with the 721 number of electrodes showing interictal epileptiform discharges (IED) in parenthesis.

722
Table S3 Names of drugs, their abbreviations as used in Table S2 and their defined daily dose.