## Abstract

Brain responses vary considerably from moment to moment, even to identical sensory stimuli. This has been attributed to changes in instantaneous neuronal states determining the system's excitability. Yet the spatiotemporal organization of these dynamics remains poorly understood. Here we test whether variability in stimulus-evoked activity can be interpreted within the framework of criticality, which postulates dynamics of neural systems to be tuned toward the phase transition between stability and instability as is reflected in scale-free fluctuations in spontaneous neural activity. Using a novel noninvasive approach in 33 male human participants, we tracked instantaneous cortical excitability by inferring the magnitude of excitatory postsynaptic currents from the N20 component of the somatosensory evoked potential. Fluctuations of cortical excitability demonstrated long-range temporal dependencies decaying according to a power law across trials, a hallmark of systems at critical states. As these dynamics covaried with changes in prestimulus oscillatory activity in the alpha band (8-13 Hz), we establish a mechanistic link between ongoing and evoked activity through cortical excitability and argue that the co-emergence of common temporal power laws may indeed originate from neural networks poised close to a critical state. In contrast, no signatures of criticality were found in subcortical or peripheral nerve activity. Thus, criticality may represent a parsimonious organizing principle of variability in stimulus-related brain processes on a cortical level, possibly reflecting a delicate equilibrium between robustness and flexibility of neural responses to external stimuli.

**SIGNIFICANCE STATEMENT** Variability of neural responses in primary sensory areas is puzzling, as it is detrimental to the exact mapping between stimulus features and neural activity. However, such variability can be beneficial for information processing in neural networks if it is of a specific nature, namely, if dynamics are poised at a so-called critical state characterized by a scale-free spatiotemporal structure. Here, we demonstrate the existence of a link between signatures of criticality in ongoing and evoked activity through cortical excitability, which fills the long-standing gap between two major directions of research on neural variability: the impact of instantaneous brain states on stimulus processing on the one hand and the scale-free organization of spatiotemporal network dynamics of spontaneous activity on the other.

## Introduction

Neural responses are characterized by remarkable variability, even to identical physical stimuli. This well-known phenomenon has been attributed to fluctuations of the neuronal network's state (Arieli et al., 1996; Sadaghiani et al., 2010), observable via diverse neuronal measures, such as EEG (Rahn and Başar, 1993; Romei et al., 2008; Vanrullen et al., 2011; Forschack et al., 2017; Iemi et al., 2019), BOLD signal (Fox and Raichle, 2007; Becker et al., 2011), local field potentials (Arieli et al., 1996), and single-cell recordings (Azouz and Gray, 1999; Churchland et al., 2010). So far, studies on neuronal variability have mainly focused on the strength of variability (Churchland et al., 2010; Garrett et al., 2013; Dinstein et al., 2015). However, the dynamics of network states over time may provide even further insights into the underlying spatiotemporal organization principles.

In this context, a certain type of fluctuation pattern known as power-law dynamics, which indicates that a signal possesses scale-free properties, is of particular interest. Such power-law relationships represent a hallmark of the (self-)organization of complex systems at a critical state (Sethna et al., 2001; Muñoz, 2018), the point of a phase transition between two distinct system regimens, such as order and disorder (Bak et al., 1987, 1988; Beggs and Plenz, 2003), at which the dynamic range, information processing, and capacity of a system are maximized (Kinouchi and Copelli, 2006; Shew and Plenz, 2013). Figure 1*A* visualizes these system configurations using the Ising model of ferromagnetism (Ising, 1925).

Empirically, power-law dynamics have been found in the size and duration of neuronal avalanches of various species, such as rats (Beggs and Plenz, 2003; Friedman et al., 2012), monkeys (Petermann et al., 2009; Yu et al., 2017), zebrafish larvae (Ponce-Alvarez et al., 2018), and humans (Priesemann et al., 2013; Shriki et al., 2013; Arviv et al., 2015). In the temporal domain, power-law dynamics can be observedin human resting state fMRI networks (Tagliazucchi et al., 2013), as well as in amplitude fluctuations of alpha and beta band activity by way of MEG/EEG (Linkenkaer-Hansen et al., 2001; Palva et al., 2013),associated with long-range temporal dependencies. Thus, criticality may represent a general and appealingly parsimoniousexplanation of neuronal variability inthe brain. However, the functional link between critical dynamics, fluctuations of ongoing neural activity, and cortical excitability remains elusive, particularly for the human brain.

To establish this link, we developed an approach to probe instantaneous cortical excitability on a single-trial level using somatosensory evoked potentials (SEPs) in EEGs on humans in response to electrical median nerve stimuli (Fig. 1*B*). Specifically, the N20 component of the SEP is thought to solely reflect excitatory post-synaptic potentials (EPSPs) of the first thalamocortical volley (Peterson et al., 1995; Wikström et al., 1996; Bruyns-Haylett et al., 2017), generated in the anterior wall of the postcentral gyrus, Brodmann area 3b (Allison et al., 1991). Therefore, the amplitude of this early part of the SEP represents a direct measure of the instantaneous excitability of a well-defined neuronal population in the primary somatosensory cortex. Additionally, to bridge the gap between evoked and ongoing neuronal activity, we evaluated prestimulus oscillations in the alpha band(8-13 Hz) of the same neuronal sources, a classical index of cortical excitability in ongoing neural activity (Klimesch et al., 2007; Romei et al., 2008).

Here we show that fluctuations of cortical excitability are likely to be governed by the same near-critical network dynamics both in ongoing and evoked neural activity.

## Materials and Methods

##### Participants

EEG data were recorded from 33 male human subjects. Two subjects were excluded because no clear SEPs were visible in the single-trial analysis, probably because of suboptimal placement of the stimulation electrodes and a low signal-to-noise ratio (SNR) of the EEG. The remaining sample of 31 subjects had a mean ± SD age of 26.9 ± 5.0 years. All participants were right-handed (lateralization score, mean ± SD, 92.9 ± 11.7), as assessed using the Edinburgh Handedness Inventory (Oldfield, 1971), and did not report any neurologic or psychiatric disease. All participants gave informed consent and were reimbursed monetarily. The study was approved by the local ethics committee.

##### Stimuli

Somatosensory stimuli were applied using electrical stimulation of the median nerve. A noninvasive bipolar stimulation electrode was positioned on the left wrist (cathode proximal). The electrical stimuli were designed as squared pulses of a 20 µs duration. The stimulus intensity was set to 1.2 × motor threshold, leading to clearly visible thumb twitches for every stimulus, as individually determined by a staircase procedure before the experiment. Stimuli were applied using a DS-7 constant-current stimulator (Digitimer).

##### Procedure

During the experiment, participants were seated comfortably in a chair their hands extended in front of them in the supinate position on a pillow. Electrical stimuli were presented in a continuous sequence with interstimulus intervals ranging from 663 to 763 ms (randomly drawn from a uniform distribution; ISI_{average} = 713 ms). In total, 1000 stimuli were applied, divided into two blocks of 500 stimuli each with a short break in between. Participants were instructed to relax and fixate their gaze on a cross on a computer screen in front of them while receiving the stimuli.

##### Data acquisition

EEG data were recorded from 60 Ag/AgCl electrodes at a sampling rate of 5000 Hz using an 80-channel EEG system (NeurOne, Bittium) with a bandwidth of 0.16-1250 Hz. Electrodes were mounted in an elastic cap (EasyCap) at the international 10-10 system positions FP1, FPz, FP2, AF7, AF3, AFz, AF4, AF8, F7, F5, F3, F1, Fz, F2, F4, F6, F8, FT9, FT7, FT8, FT10, FC5, FC3, FC1, FC2, FC4, FC6, C5, C3, C1, Cz, C2, C4, C6, CP5, CP3, CP1, CPz, CP2, CP4, CP6, T7, T8, TP7, TP8, P7, P5, P3, P1, Pz, P2, P4, P6, P8, PO7, PO3, PO4, PO8, O1, and O2. Four additional electrodes were placed at the outer canthus and the infraorbital ridge of each eye to record the electro-oculogram. During recording, the EEG signal was referenced to FCz, and POz served as ground. All impedances were kept to <10 kΩ. For source reconstruction, EEG electrode positions were measured in 3D space individually for each subject using the Patriot motion tracker (Polhemus). Additionally, the compound nerve action potential (CNAP) of the median nerve was measured using two bipolar electrodes, positioned at the inner side of the left upper arm.

Structural T1-weighted MRI scans (MPRAGE) of all participants but two were obtained during a different testing date (within the same year of the experiment or up to 3 years earlier), on a 3T Verio, Skyra, or Prisma scanner (Siemens).

##### EEG preprocessing

Stimulation artifacts were cut out and interpolated between −2 and 4 ms relative to stimulus onset using Piecewise Cubic Hermite Interpolating Polynomials. The EEG data were bandpass filtered between 30 and 200 Hz, sliding a fourth-order Butterworth filter forward and backward over the data to prevent phase shift. With this filter, we specifically focused on the N20-P35 complex of the SEP, which emerged from frequencies of >35 Hz as indicated by time-frequency analyses on pilot data, and omitted contributions of later (slower) SEP potentials of no interest. Furthermore, this filter effectively served as baseline correction of the SEP since it removed slow trends in the data, reaching an attenuation of 30 dB at 14 Hz, thus ensuring that fluctuations in the SEP did not arise from fluctuations with slower frequencies (e.g., alpha band activity). (The relationship between decibels [dB] and magnitude is defined as ). Bad segments of the data were removed by visual inspection, resulting in 989 trials on average per participant. The data were then rereferenced to an average reference. Eye artifacts were removed using Independent Component Analysis. For analysis of SEPs, the data were segmented into epochs from −100 to 600 ms relative to stimulus onset. EEG preprocessing was performed using EEGLAB (Delorme and Makeig, 2004), and custom-written scripts in MATLAB (The MathWorks).

##### Single-trial extraction using canonical correlation analysis (CCA)

Single-trial SEPs were extracted by applying a variant of CCA, as previously proposed by Waterstraat et al. (2015). CCA is used for finding weights *w _{x}* and

*w*that mutually maximize the correlation between two signals

_{y}*X*and

*Y*, so that:

For extracting single-trial SEPs, we constructed *X* as a two-dimensional matrix (time × channel) containing all single-trial epochs (concatenated in the time domain), whereas *Y* contained the average SEP, concatenated as often as there were epochs (also concatenated in the time domain). The resulting weight matrix *w _{x}* represents spatial filters that, in combination with

*w*, maximize the correlation between single-trial activity (

_{y}*X*) and the average SEP (

*Y*). By considering all concatenated single trials at once, latent spatial components are identified which are shared across single trials and best approximate the spatiotemporal profile of the overall SEP. This approach, however, does not require a match in the amplitude between averaged and single-trial SEPs, but only the temporal profile of the spatially projected single trials being similar (i.e., maximally correlated) to the average SEP. To particularly focus on the early portion of the SEP, the spatial filters

*w*were trained using shorter segments from 5 to 80 ms after stimulus but applied to the entire epochs from −100 to 600 ms. We derived a number of spatially distinct components by applying the spatial filters to the single-trial matrix, here denoted as CCA components as follows:

_{x}To characterize the CCA components in more detail, their spatial patterns were computed as follows:
and components were visually identified that showed a tangential spatial pattern over the central sulcus as is typical for the N20-P35 complex (referred to as tangential CCA components). Furthermore, components were identified that showed a peak in the activity time course at 15 ms (referred to as thalamic CCA components; only in a subset of the sample). This procedure was performed individually for every subject for the first four CCA components, as sorted by their canonical correlation coefficients. Since CCA is insensitive to the polarity of the signal, the resulting tangential CCA components were standardized so that the N20 always appeared as a negative peak in the SEP (i.e., by inverting their spatial filters *w _{x}*, if necessary). Furthermore, CCA is insensitive to the order of trials. Thus, the same spatial filters

*w*are obtained when permuting the order of single-trial SEPs, and it is therefore not possible that CCA influences the temporal structure of SEP amplitudes across trials.

_{x}##### SEP peak amplitudes and prestimulus oscillatory activity

N20 peak amplitudes were defined as the minimum value in single-trial SEPs of the tangential CCA components ±2ms around the latency of the N20 in the within-subject average SEP. P35 peak amplitudes were defined accordingly as the maximum around the latency of the P35 in the within-subject average SEP.

To extract prestimulus alpha band activity, the data were first cut into segments from −500 to −5 ms relative to stimulus onset. Next, the data segments were mirrored to both sides (symmetric padding) and bandpass filtered between 8 and 13 Hz (fourth-order Butterworth filter forward and backward applied). By segmenting the data before filtering, we ensured that there was no contamination by poststimulus activity. Mirroring the data segments served to minimize filter-related artifacts at the edges of the prestimulus segments. To make a direct comparison with the SEP possible, we applied the spatial filter corresponding to the tangential CCA component to the prestimulus data. Subsequently, the prestimulus alpha envelope was measured by taking the absolute values of the signals processed with the Hilbert transform. To derive one prestimulus alpha metric for every trial, amplitudes of the alpha envelope were averaged in the prestimulus time window of interest between −200 and −10 ms.

##### EEG source reconstruction

To reconstruct the sources of the EEG signal, we estimated lead field matrices based on individual brain anatomies and individually measured electrode positions. The structural T1-weighted MRI images (MPRAGEs) were segmented using the Freesurfer software (http://surfer.nmr.mgh.harvard.edu/), and a 3-shell boundary element model was constructed which was used to compute the lead field matrix with OpenMEEG (Kybic et al., 2005; Gramfort et al., 2010). For 2 subjects, a template brain anatomy (ICBM152; Fonov et al., 2009) was used as no individual MRI scans were available. For 1 subject, standard electrode positions were used instead of individually measured positions. The lead field matrices were inverted using eLORETA (Pascual-Marqui, 2007), and sources were reconstructed for the spatial patterns of the tangential CCA component of every subject. Next, individual source spaces were transformed into a common source space based on the ICBM152 template using the spherical coregistration with the FSAverage atlas (Fischl et al., 1999) derived from Freesurfer, to average the obtained sources of the CCA components across subjects. The calculation of the individual head models and visualization of the sources were performed using Brainstorm (version 3.4; Tadel et al., 2011). The MATLAB implementation of the eLORETA algorithm was derived from the MEG/EEG Toolbox of Hamburg (METH).

##### Processing of peripheral electrophysiological data (median nerve CNAP)

Analogously to the EEG data, stimulation artifacts were cut out and interpolated between −2 and 4 ms relative to stimulus-onset using Piecewise Cubic Hermite Interpolating Polynomials. Next, the data were high-pass filtered at 70 Hz to extract the short-latency CNAP peak of only a few milliseconds duration, applying a fourth-order Butterworth filter forward and backward to prevent phase shift. Additionally, notch filters (fourth-order Butterworth) were applied from 48 to 52 Hz and 148 to 152 Hz, respectively. Epochs were extracted from −100 to 600 ms relative to stimulus onset.

##### Detrended fluctuation analysis (DFA)

Power-law dynamics in the fluctuations of early SEPs as well as of prestimulus alpha band activity were quantified using DFA (Peng et al., 1994; Kantelhardt et al., 2001; Hardstone et al., 2012). DFA calculates the fluctuation (i.e., SD) of a cumulative signal on different time scales and tests whether its distribution follows a power law: , where F denotes the fluctuation function, τ the signal length (or window size), and α the power-law exponent. The DFA exponent α quantifies the extent of power-law dynamics of a signal, with α > 0.5 indicating persistent auto-correlations; whereas α = 0.5 is expected for a signal without a correlated temporal structure (i.e., white noise).

DFA was chosen since it requires less strict assumptions about the stationarity of the signal than the classical auto-correlation function or power spectral density (Linkenkaer-Hansen et al., 2001; Hardstone et al., 2012; Sangiuliano Intra et al., 2018), and its power-law exponent estimates are robust even when removing large portions of the data or when analyzing data that have been stitched together (Chen et al., 2002). Furthermore, the choice of DFA allowed direct comparisons with the existing literature on scale-free behavior in neuronal data.

We analyzed power-law dynamics in the fluctuation of SEP and prestimulus alpha amplitudes across trials with window sizes ranging from 7 to 70 trials, which correspond to time windows of ∼5 to 50 s. The same temporal window sizes were selected for the DFA of continuous alpha band activity. To examine whether a power-law distribution was an appropriate model for the observed dynamics, we performed model comparisons as described by Ton and Daffertshofer (2016), similarly to suggestions by Botcharova et al. (2014) and Clauset et al. (2009). In case of the existence of a power law, the relation between the log-transformed amplitude fluctuation and the log-transformed window size should be linear. This linear log-log model was compared with a number of alternative models, comprising higher-order polynomial, exponential, and logarithmic functions, as well as splines with two or three linear sections (for more details see Ton and Daffertshofer, 2016). These comparisons were performed at every latency of the early SEP, fitting the models by maximum likelihood estimation and evaluating the model fit by the Bayesian Information Criterion.

##### Evaluation of SNR

The SNR of the single-trial SEP, as measured by the tangential CCA component, was quantified as the quotient of the root-mean-square signal in the time range of the SEP (10-50 ms) and a prestimulus baseline (−50 to −10ms), so that .

The same procedure was applied to estimate the SNR of the CNAP and of the thalamic CCA component. For the CNAP, we chose time windows from 5 to 8 ms and −8 to −5ms, and for thalamic activity 12 to 18 ms and −18 to −12ms, to estimate signal and noise, respectively.

##### Simulation of the relationship between SNR and DFA exponent

Signals with DFA exponents systematically varying in the range from α = 0.5 to α = 0.8 were generated by filtering white noise with IIR filters whose coefficients depended on the desired DFA exponents as described by Schaworonkow et al. (2015) according to the algorithm of Kasdin (1995). The length of these time series was set to 1000 data points corresponding to our empirical data from the SEP fluctuation across trials. These time series were mixed with white noise, that is, stochastically independent time series with DFA exponents of α = 0.5. The time series with varying DFA exponents were mixed with the noise at varying SNRs ranging from 0.001 to 6, defined as . This procedure was repeated 100 times to account for the variance in the generation of random time series. Subsequently, DFA exponents of the mixed time series were measured, and the average DFA exponent of the simulated signal was identified for which the SNR and DFA exponent of the mixed time series corresponded to our empirical analysis of SEP fluctuations.

##### Simulation of the influence of temporal filtering on DFA exponents

To confirm that our temporal filtering did not cause the DFA exponent increases in the early SEP, we applied the same filtering to surrogate data with stochastically independent SEP fluctuations. SEP fluctuations across trials were simulated by decreasing or increasing an average SEP time course by a randomly generated factor for every trial. These signals were superimposed on continuous pink noise, which was bandpass filtered between 30 and 200 Hz (fourth-order Butterworth filter applied forward and backward), using a SNR of 2, a typical value for empirical data. Subsequently, DFA was applied across trials for every time sample of the simulated SEP, corresponding to above-described DFA analyses of the empirical SEPs.

##### Statistical analyses

We compared the empirical DFA exponent time courses with surrogate data and applied cluster-based permutation tests to assess whether, and at which latencies, DFA exponents were significantly higher than it would be expected for stochastically independent fluctuation (i.e., white noise). First, we determined the expected DFA exponents for stochastically independent fluctuation by shuffling the trial order of our data and applying DFA to it. To account for variability because of random shuffling, this step was repeated 1000 times, and DFA exponents of these iterations were averaged, yielding an average surrogate DFA exponent time course for every subject. (Averaged across all samples and subjects, the mean DFA exponent was α = 0.512, thus slightly increased compared with the theoretical DFA exponent of white noise of α = 0.5. This small empirical deviation may have been caused by the asymptotic behavior of DFA for small window sizes.) Next, the DFA exponents of the data with intact trial order were compared with the average DFA exponents of the surrogate data, using a two-sample *t* test, resulting in a *t* value for every comparison over the time course of the SEP. To obtain clusters of increased DFA exponents, *t* values were thresholded at *p _{pre}* = 0.001. Within clusters,

*t*values were summed up to cluster

*t*values

*t*. The same procedure was repeated 1000 times for the surrogate data, always comparing one surrogate dataset with the average surrogate data, which provided us with the distribution of cluster

_{cluster,empirical}*t*values under the null hypothesis. Next, a cutoff value

*t*was defined at the 99.9th percentile corresponding to a cluster threshold of

_{cluster,crit}*p*= 0.001. Finally,

_{cluster}*t*of all clusters in the empirical DFA exponent time course was compared with

_{cluster,empirical}*t*to identify clusters of significantly increased DFA exponents. With this procedure, we controlled for the number of samples over the SEP time course, intersubject variability, and the distribution of amplitude values of the SEP from which DFA exponents were derived.

_{cluster,crit}Analogously, DFA exponents of prestimulus alpha band activity were statistically tested using a *t* test on group level, comparing them with the average DFA exponents of the subject-specific null distributions, which were calculated from 1000 surrogate datasets with shuffled trial order each. The statistical significance of the DFA exponents of continuous alpha band activity was compared with a *t* test to the average DFA exponent of 1000 surrogate datasets that were derived from generating white-noise signals of comparable length as the EEG recordings (3,600,000 data points), filtering them in the alpha band range and extracting alpha amplitudes in the same way as for the empirical data (corresponding to the procedure suggested by Nikulin and Brismar, 2004).

To test the relationship of SNR and DFA exponents, we correlated the average SNR of single-trial SEPs with the average DFA exponents between 10 and 50 ms after stimulus, across participants using Spearman correlation.

Furthermore, we assessed the relationship between single-trial N20 peak amplitudes and prestimulus alpha amplitudes using a linear-mixed-effects model of the following form:
with subject as random factor, estimating the fixed effect as well as the random slope of the predictor prestimulus alpha amplitude with the dependent variable N20 peak amplitude (intercepts were included both for the fixed and random effects). Additionally, the relationship between N20 peak amplitudes and prestimulus alpha amplitudes was assessed with a permutation-based approach in which we compared their Spearman correlation coefficients with those from surrogate prestimulus alpha amplitudes with the same auto-correlated structure but shuffled phases generated by Adjusted Amplitude Fourier Transform (Theiler et al., 1992), as suggested by Schaworonkow et al. (2015). Empirical correlation coefficients were averaged across subjects after Fisher's *z* transformation and compared with the null distribution of 10,000 averaged correlation coefficients from the surrogate analyses to obtain the corresponding *p* value.

The covariation of N20 and P35 peak amplitudes was tested using a random-slope linear-mixed-effects model with P35 peak amplitude as dependent variable, N20 peak amplitude as independent variable, and subject as random factor.

Finally, we correlated DFA exponents of alpha activity (both prestimulus and continuous) and DFA exponents of the SEP. To account for the variability of the DFA exponent time course in the early SEP, we calculated the arithmetic mean of DFA exponents in four subsequent time windows, 20–25, 25–30, 30–35, and 35-40 ms, and computed their Spearman correlation coefficients with the DFA exponents of alpha activity (prestimulus and continuous), respectively. To control for the multiple testing across subsequent time windows, we applied Bonferroni correction.

For all statistical analyses, the significance level was set to *p* = 0.05. Correlation analyses as well as permutation-based statistics were performed in MATLAB (version 2017b, The MathWorks). For all correlational relationships between different parameters, we used Spearman's correlation coefficient as a measure robust to deviations of the data from the normal distribution. The DFA model comparisons were performed using publicly available code provided by Ton and Daffertshofer (2016). The linear-mixed-effects models were calculated in R (version 3.5.1, R Core Team, 2018) using the *lmer* function of the *lme4* package (version 1.1-21; Bates et al., 2015) with the default options (including Restricted Maximum Likelihood estimation of unstandardized fixed-effect coefficients). To derive a *p* value for the fixed-effect coefficients, the denominator degrees of freedom were adjusted using Satterthwaite's method (Satterthwaite, 1946) as implemented in the R package *lmerTest* (version 3.1-1; Kuznetsova et al., 2017).

##### Data and code availability

The data that support the findings of this study are available on request from the corresponding author (T.S.). The data cannot be made publicly available because of the privacy policies for human biometric data according to the European General Data Protection Regulation.

The custom-written code that was used for data processing and statistical analyses is publicly available at https://osf.io/jzqdt/?view_only=dcb94617841445859adc5496d33c3cee.

## Results

### SEPs and neuronal generators

The SEPs, averaged across all participants and trials, are shown in Figure 2*A*. The N20 component is visible as a negative peak at ∼20 ms at electrodes contralateral to the stimulation site and posterior to the central sulcus. Furthermore, the scalp topography at 20 ms shows a tangential dipole centered over the central sulcus (Fig. 2*B*), consistent with the assumption of neuronal generators located in the anterior wall of the postcentral gyrus.

To extract single-trial SEPs, we used a variant of CCA in which spatial filters were trained based on a pattern matching between average SEP and single trials (Fedele et al., 2013; Waterstraat et al., 2015). With this method, we obtained a set of spatially distinct CCA components for every individual subject. In all subjects, a prominent CCA component was identified that displayed the pattern of the typical N20 tangential dipole (Fig. 2*C*) and showed a clear peak at ∼20 ms after stimulus (Fig. 2*A*). Furthermore, subsequent source reconstruction of the spatial pattern revealed that the strongest generators were located in the anterior wall of the postcentral gyrus (Fig. 2*D*,*E*). We focus on this CCA component in the following analyses and refer to it as the tangential CCA component.

Single-trial SEPs retrieved from the tangential CCA component are displayed in Figure 3*A* for an exemplary subject. It is apparent that the amplitude of the early SEP fluctuates over trials, however without a clear trend (Fig. 3*B*).

### Temporal dynamics in single-trial SEP amplitude fluctuations

To evaluate the characteristics of SEP fluctuations across trials, we applied DFA (Kantelhardt et al., 2001; Hardstone et al., 2012). The DFA exponent α quantifies the extent of power-law dynamics of a signal, with α > 0.5 indicating persistent auto-correlations, whereas α = 0.5 would suggest a signal without a temporal structure (i.e., white noise). DFA was performed on the amplitudes at every latency relative to the stimulus onset, across time windows of 7–70 trials (i.e., equivalent to ∼5–50 s; exemplarily illustrated in Fig. 3*B*,*C*). Applying DFA at all latencies provided a DFA exponent time course for every subject (Fig. 3*D*) that indicates which portions of the early SEP show power-law dynamics (i.e., DFA exponents > 0.5). Subsequently, DFA exponent time courses were averaged across participants (Fig. 3*E*). The average explained variance of the power-law relationships at all latencies of the time range between 10 and 50 ms was *R*^{2} > 0.99, indicating a near-perfect fit of the DFA method for these data.

Increased DFA exponents were observed particularly in the early part of the SEP with an onset around the latency of the N20 component, whereas surrogate data generated by shuffling the trial order yielded DFA exponents close to α = 0.5. Two prominent peaks in the DFA exponent time course are visible from Figure 3*E*, with DFA exponents of α = 0.575 and α = 0.577, at latencies of 25 and 33 ms after stimulus, respectively. However, the absolute value of DFA exponents highly depends on the SNR of the signal, as is further examined in simulations below, which suggest DFA exponents of at least α = 0.63 when the SNR bias is taken into account. The observation of two prominent DFA exponent peaks was statistically confirmed as two main clusters were found around these two peaks by cluster-based permutation tests (*p*s_{cluster} < 0.001). The DFA exponents were characterized by a similar, yet not identical, time course compared with the magnitude of the SEP (Fig. 3*E*). Although the first significant DFA exponent cluster emerged together with the peak of the N20 component, the first DFA exponent peak appeared slightly later. This suggests that long-range temporal dependencies were not most pronounced at the N20 peak but rather while the potential returned back to baseline. Yet, this is not contradicting the notion of power-law dynamics in fluctuations of cortical excitability since recent evidence from pharmacological studies suggests that excitatory processes of the N20 component dominate even until the rising flank of the P35, the next component after the N20 in the SEP (Bruyns-Haylett et al., 2017).

The second DFA time course peak co-occurred with the second prominent peak of the SEP, the P35 component. This second DFA peak most likely reflects activity propagated from the N20 component to the P35 as these two components moderately covaried in our data, β_{fixed} = –0.378, *t*_{(29.800)} = −9.342, *p* < 0.001, as tested by a random-slope linear-mixed-effects model with P35 peak amplitude as dependent variable, N20 peak amplitude as independent variable, and subject as random factor. Similarly, the two smaller clusters of increased DFA exponents at ∼44 and 68 ms (Fig. 3*E*) may reflect the propagation of dynamics in earlier SEP components to later processing stages.

Temporal filtering in the preprocessing of the data cannot have caused these long-range temporal dependencies since (1) no DFA exponent increases were observed during the prestimulus baseline of the SEP and (2) additional control analyses did not show increased DFA exponents when applying the same preprocessing to stochastically independent SEP fluctuations (as tested with simulated data).

Furthermore, model comparisons based on maximum likelihood estimation (Ton and Daffertshofer, 2016) indicated that a linear model best fitted the log-log relationship between amplitude fluctuation and window size at 81.59% of the latencies in the time range from 10 to 50 ms after stimulus, as indicated by the median Bayesian Information Criterion on group level. Taking into account that the finite data length and measurement noise in single-trial amplitudes have presumably degraded the fit of the linear log-log model for the underlying neural signals (Botcharova et al., 2014), we therefore conclude that our data are largely consistent with the hypothesis that amplitude fluctuations across time are distributed according to a power law.

### Do power-law dynamics originate from the neuronal fluctuations in the periphery or at the thalamic level?

To investigate whether the observed temporal dynamics in cortical SEPs may arise from fluctuations in peripheral nerve excitability, we applied the same procedure as described above to the CNAP of the median nerve measured at the inner side of the upper arm. As expected, the nerve potential peaked at ∼6 ms after stimulus and fluctuated over trials (Fig. 4*A*). However, no increased DFA exponents were observed (Fig. 4*B*).

In addition, a CCA component was identified in 13 of the 31 participants that contained SEP activity already at 15 ms (Fig. 4*C*,*D*), most likely reflecting the P15 component of the SEP, which is thought to represent thalamus-related activity (Albe-Fessard et al., 1986). Also, the spatial pattern of this CCA component suggested a deeper and more medial source than the tangential CCA component (Fig. 4*E*). Importantly, the DFA exponents of this subcortical activity did not show any increase in the range of the P15 component (Fig. 4*D*), thus being in contrast with the DFA exponent increase for early cortical potentials.

### DFA exponents and SNR

Since it is known from previous studies that the SNR highly affects the measurement of power-law dynamics (Blythe et al., 2014), we investigated the relationship between DFA exponents and SNR in single-trial SEPs. On average, across all participants, the SNR of the tangential CCA component was , and showed a positive rank correlation with average DFA exponents in the time range from 10 to 50 ms after stimulus (*r* = 0.358, *p* = 0.048).

Additionally, we further clarified this relationship with simulations: We mixed signals expressing different DFA exponents with white noise (DFA exponent α = 0.5), for a range of SNRs, and measured the DFA exponent of these mixed signals. As is visible from Figure 4*F*, DFA exponents of the mixed signals are attenuated toward α = 0.5 when lowering the SNR. Given an SNR of 1.68 and an empirical DFA exponent of α = 0.575, as was the case for the tangential CCA component in the present study, our simulations suggest an underlying source with a DFA exponent of α ≈ 0.63. Yet this value most likely still underestimates the “true” power-law dynamics of the system as the signal term contained in the empirical estimate of the SNR is not noise-free but a mixture of both signal and noise. This leads to an overestimation of the SNR and in turn to an underestimation of the degrading impact of noise on the scaling exponent.

To relate this simulation also to the other measures for which we calculated DFA exponents, we calculated the SNR of the CNAP at the upper arm and thalamic CCA components. Here, we found SNRs of 2.20 (SD = 0.85) and 1.33 (SD = 0.11), respectively, suggesting that our signal quality was sufficient to detect DFA exponent increases if they had been there since the SNR of the CNAP was even higher than that of the SEP and the SNR of the thalamic CCA component was just slightly lower.

### Power-law dynamics in alpha band activity and its relation to the SEP

Since previous studies on cortical excitability in MEG/EEG focused on oscillatory activity in the alpha band, we investigated both its correspondence to the early part of the SEP as well as its DFA exponents.

To test the relationship between alpha oscillatory activity and SEP amplitude, we performed a regression analysis between the mean alpha amplitude in a prestimulus window from −200 to −10ms and the peak amplitude of the N20 component. Alpha activity was extracted from the same neuronal sources as the SEP by applying the spatial filter of the tangential CCA component. A significant negative relationship was found on group level using a random-slope linear-mixed-effects model (β_{fixed} = −0.034, *t*_{(25.095)} = −4.895, *p* < 0.001). Thus, higher prestimulus alpha activity was associated with more negative N20 peak amplitudes (Fig. 5*A*).

To control for spurious covariation caused by the auto-correlated structure of both signals, we additionally ran permutation tests using surrogate data with comparable temporal structure as suggested by Schaworonkow et al. (2015). Aggregated on group level, these tests confirmed the negative relationship between prestimulus alpha activity and N20 peak amplitude (*r*_{group-level} = −0.035, *p* < 0.001).

Next, we investigated the DFA exponents of mean prestimulus alpha amplitude across trials. Averaged across subjects, we observed a mean DFA exponent of α = 0.60, which significantly differed from DFA exponents for shuffled trial order (*t*_{(30)} = 6.627, *p* < 0.001). Also, DFA exponents in continuous, ongoing alpha activity were significantly higher than chance level across subjects (α = 0.66, *t*_{(30)} = 9.5327, *p* < 0.001).

To further test the relationship between the fluctuation patterns of prestimulus activity and of the SEP, we correlated DFA exponents of prestimulus alpha amplitude and DFA exponents of the SEP across participants (using Spearman correlation). DFA exponents of the SEP time series were averaged across data samples in four consecutive time windows of 5ms each, between 20 and 40 ms after stimulus. DFA exponents of prestimulus alpha activity were correlated with the DFA exponents of the first time window from 20 to 25 ms (*r* = 0.486, *p* = 0.025, Bonferroni-corrected; Fig. 5*B*). However, this relationship did not emerge for any other time window between 25 and 40 ms (all *p* > 0.3). Similarly, DFA exponents of continuous alpha activity were correlated with DFA exponents of the SEP in the first time window from 20 to 25 ms after stimulus (*r* = 0.500, *p* = 0.019, Bonferroni-corrected), but not with DFA exponents in the subsequent 5 ms time windows between 25 and 40 ms (all *p* > 0.3).

Notably, the SNR of the SEP cannot explain the relation between DFA exponents of alpha activity and DFA exponents of the SEP, as no relationship was found between SNR of the SEP and DFA exponents of prestimulus alpha activity (*r* = 0.208, *p* = 0.261).

Together, both amplitude and temporal structure of oscillatory activity in the alpha band thus relate to the corresponding characteristics of the early SEP responses, establishing alink between these two measures of instantaneous cortical excitability.

## Discussion

In the present study, we investigated the temporal dynamics of neuronal excitability in the human primary somatosensory cortex by short-latency SEPs. Fluctuations of excitability demonstrated power-law dynamics across trials, extending the previous notion that neuronal systems operate close to a critical state (Linkenkaer-Hansen et al., 2001; Beggs and Plenz, 2003; Poil et al., 2012; Palva et al., 2013; Priesemann et al., 2013) to variability in stimulus-evoked responses in the human sensory system. In addition, fluctuations in prestimulus alpha band activity and initial cortical excitation were related through their amplitudes as well as their temporal structure. For the first time, these findings thus link critical dynamics in ongoing and evoked activity as measured noninvasively in the human EEG, and directly associate the observed power-law dynamics with variability in cortical excitability.

### What do temporal dynamics in SEPs tell about the functioning of the neural system?

Given that the SEP in the time range of the N20 component reflects initial excitatory cortical processes (Peterson et al., 1995; Wikström et al., 1996; Bruyns-Haylett et al., 2017), the observed power-law dynamics indicate that instantaneous excitability does not vary stochastically independently over time (i.e., like white noise) but is characterized by long-range temporal dependencies. This means that cortical excitability at a given moment is related to its fluctuation history and contains information about subsequent dynamics up to 50 s later.

Such power-law dynamics have often been interpreted within the hypothesis that the underlying system is poised near a critical state (Bak et al., 1987, 1988; Sethna et al., 2001; Beggs and Plenz, 2003; Kitzbichler et al., 2009), at which the dynamic range, information processing, and memory capacity of a system are maximized (Kinouchi and Copelli, 2006; Shew and Plenz, 2013). Hence, it may be beneficial for neural systems to be tuned to this close-to-critical state, although the variability inherent to it may impair an exact mapping of stimulus features and neuronal activity in sensory processing.

Although well corresponding to previous reports of power-law dynamics in oscillatory activity (Linkenkaer-Hansen et al., 2001; Palva et al., 2013), our SNR-adjusted power-law exponents of α ≈ 0.63 deviated from α = 1, the expected DFA exponent for the critical state (Poil et al., 2012). This may be explained by neuronal dynamics not always being poised exactly at the critical state but rather residing in a near-critical regime (Priesemann et al., 2013), resulting from transient rebalancing of order and disorder in neural dynamics to meet environment-imposed demands (Ponce-Alvarez et al., 2018).

Typically, the power-law relationships in models of complex systems as well as in empirical neuronal avalanche recordings have been measured in the spatial domain, such as the distribution of size and duration of neuronal avalanches. Yet, critical systems should also express power-law dynamics in the temporal domain as shown for the Ising model (Zhao et al., 2017). Thus, the observed long-range temporal dependencies in cortical excitability in our data might indeed correspond to near-critical dynamics in the spatial domain as observed by Beggs and Plenz (2003) in their seminal study on the scale-free behavior of spontaneous neuronal avalanches in slices of the rat somatosensory cortex.

Specifically, the temporal power-law dynamics in our data could reflect that neuronal excitability spatially differs across the network, in agreement with the notion of scale-free neuronal avalanches, which in turn leads to the observed amplitude fluctuations in the early SEP over time. Whereas a recent approach to measure neuronal avalanches from thresholded broadband EEG data supported the notion of critical dynamics also after stimulus presentation (Arviv et al., 2015), our findings of power-law dynamics in early SEPs demonstrate near-critical dynamics, for the first time, in primarily stimulus-related processes in the human brain and suggest fluctuations of cortical excitability to be the driving underlying mechanism.

### Dissociation of temporal dynamics in the cortex from peripheral and subcortical variability

To corroborate the notion of the observed power-law dynamics being a cortical phenomenon, we examined their origin in more detail.

First, we measured power-law dynamics at the onset of the N20 component (not earlier), giving no cause to assume generators of power-law dynamics in stimulus processing upstream of the primary somatosensory cortex. Second, source reconstruction confirmed that the strongest generators of the N20 lay in the anterior wall of the postcentral gyrus (Fig. 2*E*) as is expected from the literature for the N20 component (Allison et al., 1991). Third, no power-law dynamics were present in peripheral variability as measured from the CNAP of the median nerve. Fourth, thalamus-related activity reflected in the P15 component of the SEP in a subsample of 13 subjects did not show any power-law dynamics either, suggesting that neuronal variability is stochastically independent, even at final subcortical processing stages.

Hence, we conclude that the observed power-law dynamics most likely are of cortical origin.

### Relationship between prestimulus alpha activity and initial cortex excitation

Following the idea that oscillatory activity in the alpha band reflects cortical excitability (Klimesch et al., 2007; Romei et al., 2008; Sauseng et al., 2009; Zrenner et al., 2018), we tested whether this measure in a prestimulus window was related to the initial cortex excitation as assessed by the N20 amplitude. In our data, lower amplitudes of alpha band activity, hypothesized to reflect a state of increased excitability, were associated with smaller (less negative) N20 amplitudes. Although at a first glance this finding seems to contradict the hypothesis of low alpha activity being associated with higher excitability (Klimesch et al., 2007; Jensen and Mazaheri, 2010), it may be explained in a straightforward manner by the neurophysiological basis of EEG generation.

The scalp EEG reflects relative changes in collective charge distributions resulting from neuronal activation manifested in primary post-synaptic currents (PSCs; Kandel et al., 2000; Lopes da Silva, 2004; Ilmoniemi and Sarvas, 2019). The magnitude of an EEG potential *U* emerging through synchronous activity of a well-specified neuron population, as is assumed for the N20 component of the SEP, should adhere to the following general relationship:
where I denotes the sum of local primary PSCs, the number of involved neurons, and the lead field coefficient projecting source activity to the electrodes on the scalp. Since and *LF* can be assumed to be constant for the N20 component in a sequence of unchanging median nerve stimuli, we believe that primarily I, reflecting excitatory PSCs, contributed to the amplitude variability of the early part of the SEP.

Now, assuming that states of higher neuronal excitability are associated with membrane depolarization on a cellular level, the electrical driving force for further depolarizing inward trans-membrane currents is decreased and less current is needed to reach the threshold potential for excitatory responses (Castro-Alamancos, 2009). This leads to decreased PSCs at high excitability states (Deisz et al., 1991) and would result in lower amplitudes in the EEG. Hence, one should rather expect decreased N20 components following low prestimulus alpha activity, as was the case in our data.

The relationship between prestimulus alpha activity and early SEP amplitudes is further corroborated by their corresponding degree of power-law dynamics in the time range of the SEP from 20 to 25 ms after stimulus (Fig. 5*B*). Intriguingly, the idea that this link is established via cortical excitability manifested in prestimulus membrane potentials is consistent with a recent study, which indeed demonstrated near-critical dynamics in membrane potential fluctuations in the turtle visual cortex (Johnson et al., 2019).

Although this notion should be treated with some caution as it is based on the strong assumption that mechanisms observed in single-cell recordings (in animals) can be generalized to cell populations in the human cortex, the present findings could represent the missing link between power-law dynamics on micro (single-cell) and macro (cell population) scale and would relate findings of criticality in neuronal avalanches to noninvasively measured EEG potentials in humans.

### Implications for the perspective on neural variability

Why neuronal systems express large variability, particularly in perceptual processes, has been an enduring question for many years. In this study, we show that the temporal structure of such variability provides further insights into organizing principles of large-scale neuronal dynamics in the somatosensory system. Not only is stimulus-evoked activity dependent on prestimulus neuronal states, but even the instantaneous network states themselves seem to fluctuate in a systematic manner, in our data manifested in near-critical dynamics of instantaneous cortical excitability.

This fluctuation pattern of cortical excitability corroborates the notion of a delicate equilibrium between robustness and flexibility to external stimuli (Muñoz, 2018), which is maintained through balancing excitatory and inhibitory driving forces (Beggs and Plenz, 2003; Poil et al., 2012). In line with this, neurotransmitters affecting the excitation/inhibition ratio of cortical neurons may tune these near-critical dynamics, as has been shown for alpha band activity (Pfeffer et al., 2018). This may in turn enable neural networks to adaptively adjust their excitability during stimulus processing by shifting their dynamics closer to or further away from the critical state.

In conclusion, proximity to criticality may represent a compellingly parsimonious explanation of moment-to-moment fluctuations in neural responses with potential benefits for information processing in the brain.

## Footnotes

The authors declare no competing financial interests.

V.V.N. was supported in part by the HSE Basic Research Program and the Russian Academic Excellence Project 5-100. We thank Sylvia Stasch for participant recruitment and help with data collection; and Alice Hodapp for supporting the source reconstruction of the EEG.

- Correspondence should be addressed to Tilman Stephani at stephani{at}cbs.mpg.de or Vadim V. Nikulin at nikulin{at}cbs.mpg.de