Abstract
Detection of statistical irregularities, measured as a prediction error response, is fundamental to the perceptual monitoring of the environment. We studied whether prediction error response is associated with neural oscillations or asynchronous broadband activity. Electrocorticography was conducted in three male monkeys, who passively listened to the auditory roving oddball stimuli. Local field potentials (LFPs) recorded over the auditory cortex underwent spectral principal component analysis, which decoupled broadband and rhythmic components of the LFP signal. We found that the broadband component captured the prediction error response, whereas none of the rhythmic components were associated with statistical irregularities of sounds. The broadband component displayed more stochastic, asymmetrical multifractal properties than the rhythmic components, which revealed more self-similar dynamics. We thus conclude that the prediction error response is captured by neuronal populations generating asynchronous broadband activity, defined by irregular dynamic states, which, unlike oscillatory rhythms, appear to enable the neural representation of auditory prediction error response.
SIGNIFICANCE STATEMENT This study aimed to examine the contribution of oscillatory and asynchronous components of auditory local field potentials in the generation of prediction error responses to sensory irregularities, as this has not been directly addressed in the previous studies. Here, we show that mismatch negativity—an auditory prediction error response—is driven by the asynchronous broadband component of potentials recorded in the auditory cortex. This finding highlights the importance of nonoscillatory neural processes in the predictive monitoring of the environment. At a more general level, the study demonstrates that stochastic neural processes, which are often disregarded as neural noise, do have a functional role in the processing of sensory information.
- auditory cortex
- broadband response
- mismatch negativity
- prediction error
- rhythmic components
- multiscale multifractal analysis
Introduction
Detection of novel sensory information enables adaptive interaction with the surrounding environment (Clark, 2013a; Whitmire and Stanley, 2016). In the predictive coding framework of brain functioning, this interaction is characterized by a reciprocal loop between sensory predictions and prediction error signals (Friston and Kiebel, 2009; Bastos et al., 2012). Neural mechanisms of prediction error are typically studied by presenting a series of “standard” stimuli with intermittently occurring deviant stimuli, also called “oddballs,” and by contrasting brain responses between these stimulus categories (Chennu et al., 2013; Parras et al., 2017; Lumaca et al., 2019). This way, event-related potentials (ERPs) and a range of neural oscillations have been identified as neural markers of prediction error. The most widely studied deviance in ERPs is the auditory mismatch negativity (MMN), a negative deflection of electrical event-related potential recorded on the scalp or using intracranial electrodes (Näätänen et al., 1978, 2007; Halgren et al., 1995; Rosburg et al., 2005). MMN originates from the primary auditory cortex (Alho, 1995; Alain et al., 1998; Edwards et al., 2005), and it peaks at ∼150–200 ms in humans, while the peak latencies <100 ms are typically reported in monkeys (Javitt et al., 1992; Komatsu et al., 2015; Camalier et al., 2019). In addition to MMN, prediction error responses are observed in a variety of frequency ranges including theta (3–8 Hz; Fuentemilla et al., 2008; Hsiao et al., 2009; Ko et al., 2012; Choi et al., 2013; MacLean and Ward, 2014), alpha (8–12 Hz; Ko et al., 2012; MacLean and Ward, 2014), beta (14–30 Hz; Haenschel et al., 2000; MacLean and Ward, 2014), and gamma (>30 Hz; Marshall et al., 1996; Haenschel et al., 2000; Edwards et al., 2005; Eliades et al., 2014; MacLean and Ward, 2014; Dürschmid et al., 2016) ranges.
Several interpretations could be formulated aiming to explain the abundance of prediction error responses in the frequency dimension. First of all, there could be multiple independent neural mechanisms sensitive to stimulus deviance. This suggestion, however, does not explain why there would be so many distinct mechanisms with an identical functional role. Alternatively, frequency-specific detectors of prediction error might be only partially independent, forming hierarchical cross-frequency interactions. For instance, rhythms of different frequency bands could drive each other (e.g., delta phase could modulate theta amplitude and theta phase could modulate gamma amplitude in the auditory cortex; Lakatos et al., 2005). Yet another possibility, which we pursue in the present study, is that a broad frequency range of deviance responses, including theta, alpha, beta, and gamma bands, points to a broadband prediction error response, which is not restricted to a particular frequency band, but instead is driven by an arrhythmic or asynchronous neural signal across a wide frequency range. In fact, a large number of studies reported deviance effects to run across several frequency bands (Haenschel et al., 2000; Hsiao et al., 2009; Ko et al., 2012; MacLean and Ward, 2014; Chao et al., 2018), arguably alluding to arrhythmic processing of unexpected stimuli.
The electrophysiological signal recorded by scalp EEG or local field potentials (LFPs) is a summed activity of both postsynaptic and action potentials. Postsynaptic potentials contribute to the rhythmic oscillations of different frequency bands (Buzsáki et al., 2012), reflecting neural synchrony at specific timescales. Contrary to this, empirical data analysis and modeling suggest that the average neuronal firing rate produces asynchronous, broadband changes across a wide frequency range (Medvedev and Kanwal, 2004; Henrie and Shapley, 2005; Mukamel et al., 2005; Hwang and Andersen, 2011; Buzsáki et al., 2012; Li et al., 2019; Guyon et al., 2021). Such rhythmic and broadband components of the LFP signal can be decomposed using spectral principal component analysis (PCA; Miller et al., 2009a,b, 2017; Miller, 2019), in this way separating synchronous and asynchronous neural activity. The broadband component of the LFP power spectrum is commonly characterized by a power law function (Freeman and Zhai, 2009; He, 2014; Hermes et al., 2019), which reflects the lack of any specific temporal beat (e.g., 10 Hz) in the signal. Contrary to this, rhythmic components produce frequency-specific spectral peaks that deviate from the power law. In fact, the electrocorticography (ECoG) power is characterized by at least three different power law regions of which the transitions vary across individuals and recordings in humans (He et al., 2010; Chaudhuri et al., 2018) and nonhuman primates. The functional relevance of this heterogeneous scaling is discernible as, for instance, levels of arousal across a gradual progression from awake to anesthesia (Gifani et al., 2007) or deep sleep (Ma et al., 2006; Weiss et al., 2009) can manifest selectively within power law changes at different timescales. Such complex dynamics across different LFP timescales can be characterized by multiscale multifractal analysis (MMA; Gierałtowski et al., 2012), which was developed to analyze signal fluctuations on a wide range of timescales like those observed in LFP signals.
In the present study, we aimed to assess whether such broadband neural dynamics rather than frequency-specific rhythms underlie prediction error in the auditory cortex in the primate brain. Furthermore, we sought to contrast the multiscale dynamics of the broadband LFP component to that of the rhythmic components. Importantly, while several previous LFP and EEG studies have linked MMN to gamma-range activity (Marshall et al., 1996; Haenschel et al., 2000; Edwards et al., 2005; Eliades et al., 2014; MacLean and Ward, 2014; Dürschmid et al., 2016), and often referred to it as a “broadband” neural signal, the spectral power of the raw signal conflates genuinely asynchronous broadband activity with gamma oscillatory activity (for a difference in visual processing, see Hermes et al., 2019). Thus, it remains uncertain whether the reported gamma band correlates of MMN are driven by gamma oscillations or nonoscillatory broadband activity. Addressing this issue, we used spectral PCA analysis, which decoupled broadband activity from the oscillatory components of the neural signal and enabled us to assess their distinct functional associations with MMN.
Our key hypotheses were based on the predictive coding explanation for the mismatch negativity, or indeed any violation or omission response. In brief, in predictive coding formulations there are two neurocomputational mechanisms in play. First, mismatches between ascending sensory afferents and descending predictions are registered by prediction error units (usually considered to be superficial pyramidal cells; Bastos et al., 2012; Shipp, 2016). These prediction error responses are then passed to higher hierarchical levels of processing to drive Bayesian belief updating, optimize descending predictions, and resolve prediction errors in a continuous process of recurrent hierarchical message passing. According to this framework, a change in the predicted frequency would be reported by prediction error units tuned to both the new and old (i.e., oddball and standard) frequencies. Physiologically, these responses are mediated by changes in neuronal firing rates; either over broad frequency bands (e.g., as measured by multiunit activity) or, on some accounts, high gamma frequencies (Bastos et al., 2015a,b).
The second neurocomputational mechanism of predictive coding corresponds to an optimization of the gain or excitability of prediction error populations, immediately following an unpredicted sensory input. Computationally, this corresponds to increasing the precision of prediction errors, so that they have a greater influence on belief updating (Feldman and Friston, 2010; Clark, 2013b). In engineering, this corresponds to a change in the Kalman gain (Rao and Ballard, 1999). Physiologically, this is manifested in terms of an increased postsynaptic gain or loss of self-inhibition. Electrophysiologically, this change in inhibition–excitation balance would be expressed in terms of a (deviant-induced) dynamic instability, of the sort characterized by scale invariance. Please see the study by Friston et al. (2012) for a mathematical analysis in terms of Lyapunov and scaling exponents (see their Fig. 7 for a detailed simulation of induced auditory responses).
On the predictive coding account, the first mechanism rests on top-down predictions based on a generative model of the auditory stream. This provides a formal description of the model adjustment hypothesis (Garrido et al., 2009). The second mechanism depends on changing postsynaptic sensitivity, which offers a formal explanation for sensory-specific adaptation and related phenomena (e.g., spike rate adaptation). From our perspective, both mechanisms make a clear prediction: responses to deviant stimuli will be expressed in broadband (nonoscillatory) response, reflecting the activity of prediction error populations. Furthermore, such broadband signals should show evidence of increasing dynamic instability compared with the oscillatory responses as measured in terms of scaling exponents. Based on the current knowledge reviewed above about the two mechanisms of predictive coding, we make the following predictions: (1) auditory oddball responses to deviant stimuli will be expressed in broadband (nonoscillatory) component of an electrophysiological/ECoG signal, reflecting the activity of prediction error neural populations; and (2) these broadband responses will show increased dynamic instability relative to that of frequency-specific rhythms, as measured in terms of scaling exponents at several timescales. We tested these two hypotheses using PCA-based spectral decoupling PCA and multiscale multifractal analysis. Regarding broadband MMN response, while predictive processing entails both local predictions across adjacent stimuli as well as global predictions across different sequences of stimuli (Chennu et al., 2013, 2016), we focused on the relatively fast local prediction errors in the present study.
Materials and Methods
Subjects
We tested three adult male common marmosets (Callithrix jacchus) that weighed 320–380 g. Monkeys were implanted with an ECoG electrode array under general anesthesia, and all efforts were made to minimize suffering. All surgical and experimental procedures were performed under the National Institutes of Health Guidelines for the Care and Use of Laboratory Animals and were approved by the RIKEN Ethical Committee (No. H26-2–202). ERP data of one monkey (Fr) was reported previously (Komatsu et al., 2015), whereas datasets of monkeys Go and Kr are new.
Implantation of ECoG arrays
Chronically implanted, customized multichannel ECoG electrode arrays (Fig. 1B–E; Cirtech) were used for neural recordings (Komatsu et al., 2015, 2017). We implanted 32 (the left hemisphere of monkey Fr), 64 (the right hemisphere of monkey Go), and 64 (the right hemisphere of monkey Kr) electrodes in the epidural space. For the 32-electrode array, each electrode contact was 1 mm in diameter and had an interelectrode distance of 2.5–5.0 mm (Komatsu et al., 2015). For the 64-electrode array, each electrode contact was 0.6 mm in diameter and had an interelectrode distance of 1.4 mm in a bipolar pair (Komatsu et al., 2017). The electrode array covered the frontal, parietal, temporal, and occipital lobes. The additional four electrodes of monkey Fr covered part of the right frontal lobe. The animals were initially sedated with butorphanol (0.2 mg/kg, i.m.), and surgical anesthesia was achieved with ketamine (30 mg/kg, i.m.) and medetomidine (350 μg/kg, i.m.). The animals were then positioned in a stereotaxic frame (Narishige) and placed on a heating pad during surgery. Vital signs were monitored throughout surgery. Implantation of the electrode arrays involved the removal of a bone flap (∼2 cm along the anterior–posterior axis, ∼1 cm along the mediolateral axis) over the parietal cortex. The array was advanced into the epidural space. After positioning the electrode array, connectors were attached to the bone using dental acrylic and titanium (size, 1.0 × 0.1 mm) or PEEK (size, 1.4 × 2.5 mm) screws. The reference electrodes were placed in the epidural space, and the ground electrodes in the episkull space. The anti-inflammatory corticosteroid dexamethasone (1.25 mg/kg, i.m.) was administered after surgery to prevent brain swelling. The animals were given antibiotics and analgesics daily for 5 d after surgery. Following the recovery of the animals, the position of each electrode in the arrays was identified based on computer tomography, and then coregistered to a template T1-weighted anatomic magnetic resonance image (MRI; http://brainatlas.brain.riken.jp/marmoset/; Hikishima et al., 2011; monkey Fr) or preacquired MRI (monkeys Go and Kr) using MRIcron software (http://www.mricro.com; Rorden et al., 2007). In all monkeys, the electrode array covered the frontal, parietal, occipital, and temporal cortices, including the primary auditory area (Fig. 1F–I).
Stimuli and task
We adopted a roving oddball paradigm (Cowan et al., 1993; Haenschel et al., 2005; Garrido et al., 2008). The trains of 3, 5, or 11 repetitive single tones of 20 different frequencies (250–6727 Hz with intervals of one-quarter octave) were pseudorandomly presented. Tones were identical within each tone train but differed between tone trains (Fig. 1A). Because tone trains followed on from one another continuously, the first tone of a train was considered to be an unexpected deviant tone, because it was of a different frequency than that of the preceding train. The final tone was considered to be an expected standard tone because it was preceded by several repetitions of this same tone. To avoid analytical artifacts stemming from differences in the number of standard and deviant stimuli, we considered only the last tone of a train as standard. There were 240 changes from standard to deviant tones in a single recording session. Pure sinusoidal tones lasted 64 ms (7 ms rise/fall), and stimulus onset asynchrony was 503 ms. Stimulus presentation was controlled by MATLAB (MathWorks) using the Psychophysics Toolbox extensions (Brainard, 1997; Pelli, 1997; Kleiner et al., 2007). Tones were presented through two audio speakers (Fostex) with an average intensity of 60 dB SPL around the ear of the animal.
ECoG recording and preprocessing
ECoG recordings were taken in the passive listening condition while monkeys were awake. In each recording session, the monkey Fr was held in a drawstring pouch, which was stabilized in a dark room, and the monkeys Go and Kr sat on a primate chair in a dimly lit room. The length of a single session was ∼15 min: the first 3 min of data were used for another auditory paradigm consisting of standard tones (data not shown), and the remaining 12 min of data were used for the roving oddball sequences. For monkey Fr, data from three sessions were used for analysis, which resulted in 720 (240 × 3) standard and deviant presentations. For monkeys Go and Kr, data from six sessions were used for analysis, which resulted in 1440 (240 × 6) standard and deviant presentations.
ECoG signals were recorded using a multichannel data acquisition system (Cerebus, Blackrock Microsystems) with a bandpass filter of 0.3–500 Hz and then digitized at 1 kHz. In the signal preprocessing, those signals were rereferenced using an average reference montage and high-pass filtered >0.5 Hz using a sixth-order Butterworth filter. We segmented datasets from −903 to 400 ms relative to the onset of the unexpected tone, so that each segment would include a pair consisting of a deviant and a standard immediately preceding the deviant, as well as a baseline of 400 ms preceding the standard tone (Fig. 1A). The segments were then further divided into standard epochs and deviant epochs (−400 to 400 ms, with 0 ms indicating the onset of tone). While we intended to carry out the main analyses on shorter epochs in the present study (−100 to 350 ms; see below), we initially created wider epochs that were used for different analyses not reported here. For the spectral decoupling, ECoG segments of the longest 11-tone trains were created [−200 to 5500 ms, with 0 ms indicating the onset of the first tone of the sequence (deviant tone), see below]. Parts of the dataset are shared in the public server Neurotycho.org (http://neurotycho.org/; Nagasaka et al., 2011).
Functional localization of electrodes-of-interest: event-related potentials and time–frequency analysis
ECoG electrode-of-interest was identified functionally by contrasting LFP between standard and deviant stimuli (0–350 ms), separately for each electrode (Fig. 2A). For functional localization of MMN waveforms, and only for it, a low-pass filter of 40 Hz was applied, using a sixth-order Butterworth filter. ECoG recordings were rereferenced with respect to the common average reference across all electrodes. Data were then epoched around the onset of tones (−100 to 350 ms), and baseline correction was applied by subtracting the mean of the 100 ms period before the stimulus onset. MMN was assessed by comparing the standard ERP and deviant ERP, and one electrode with the largest MMN amplitude was selected for each monkey for further analyses. In all three monkeys, the identified electrode-of-interest was located in the auditory cortex (Fig. 1G–I). In addition, given that MMN is associated with high-gamma response (Marshall et al., 1996; Haenschel et al., 2000; Edwards et al., 2005; Eliades et al., 2014; MacLean and Ward, 2014; Dürschmid et al., 2016), we inspected whether the electrodes with peak MMN amplitude would also show the highest gamma power difference between standard and deviant tones. We plotted time–frequency charts between standard and deviant stimuli (0-350 ms), separately for each electrode. Epochs were first bandpass filtered in multiple successive 10-Hz-wide frequency bins (from 1 to 250 Hz) using a zero phase shift noncausal finite impulse filter with 0.5 Hz roll-off. Next, for each bandpass-filtered signal, we computed the power using a standard Hilbert transform. Finally, the resulting time–frequency charts were z score transformed from −100 to 0 ms. As expected, the largest gamma effect was indeed observed in the electrodes with the largest MMN amplitude (Fig. 2A).
Decoupling the cortical spectrum to isolate Broadband and Rhythmic spectral components
To extract the course of broadband spectral activity, we conducted the spectral decoupling of raw LFP signal (Miller et al., 2009a,b, 2017; Miller, 2019). As for the ERP analysis, ECoG potentials were rereferenced with respect to the common average reference across all electrodes. For the selected electrodes-of-interest (see above), discrete epochs of power spectral density (PSD) were calculated using the long epochs corresponding to the 11-tone sequences [−200 to 5500 ms around the first tone (deviant tone) of each sequence]. For each trial (i.e., each 11-tone sequence), individual PSDs were normalized with an element-wise division by the average power at each frequency, and the obtained values were log transformed. The eigenvectors [principal spectral components (PSCs)] from this decomposition revealed distinct components of cortical processing, with the first three PSCs explaining most variance: the broadband PSC (i.e., no predominant peak present in the PSD); the first rhythmic PSC (alpha rhythm, ∼ 10 Hz); and the second rhythmic PSC (delta rhythm, ∼2 Hz; Figs. 2C, 3A,H,O).
To identify components of stimulus-related changes in the PSD (Fig. 3), an inner product matrix of the normalized PSDs previously performed on the 11-tone sequences was diagonalized with a singular value decomposition; and was then applied to individual epochs of −100 to 350 ms in length, with 0 ms indicating the stimuli onset for both standard and deviant tones (Fig. 3B,G,I,N,P,U), and to the 11-tone sequences themselves (Fig. 3C–E,J–L,Q,S). Following the study by Miller et al. (2009b), the entire time series was z scored per trial (reported in intuitive units, because this measure is approximately normally distributed) and exponentiated, and then a value of 1 was subtracted (setting the mean closer to 0). Then, to make both conditions comparable by setting the baseline period to 0, we further performed a baseline correction on the prestimulus period by subtracting the mean value per trial between −100 and 0 ms. The first PSC allowed us to obtain the “broadband time course,” which has been shown to reflect a power law in the cortical PSD (Miller et al., 2009a), and the second and third PSCs uncovered the “rhythmic time courses” with distinct frequency peaks.
Cross-individual decoding
To assess the generalizability of our findings in the temporal domain, we conducted between-subjects cross-decoding, known to be a more conservative method of validation compared with within-subjects decoding (Gevins et al., 1998; Visell and Shao, 2016). A univariate temporal decoding model was applied on each PSC time course on the selected auditory cortex electrodes, aiming to decode the stimuli expectancy categories (i.e., standards vs deviants; Fig. 4). The ADAM-toolbox was used on the broadband and rhythmic PSC time courses with epochs from −100 to 350 ms (Fahrenfort et al., 2018). Crucially, and for each individual neural component, we trained a linear discriminant analysis (LDA) classifier in one monkey and tested it in a separate monkey for obtaining cross-individual decodability of stimuli expectancy category (i.e., standard vs deviant trials). Next, a backward decoding algorithm, using either stimulus expectancy category was applied, training the classifier on one monkey and testing on another monkey. An LDA was used to discriminate between stimulus classes (e.g., deviant vs standard trials) after which classification accuracy was computed as the area under the curve (AUC), a measure derived from signal detection theory. This allowed us to determine whether standard could be discriminated from deviant better than chance. Instead of relying on a hypothetical chance-level AUC, we nonparametrically computed the distribution of actual chance-level accuracies for every time point by running 500 iterations, randomly permuting the class labels on every iteration. This in turn allowed us to obtain a p value for every time sample. Each p value was computed by dividing the number of times that the AUC value under random permutation exceeded the actually observed AUC value and dividing this value by the total number of random permutations. These nonparametrically obtained p values were corrected for multiple comparisons over time, using a false discovery rate (FDR) correction (p < 0.05, q = 0.05; Benjamini and Hochberg, 1995; Benjamini and Yekutieli, 2001), thus yielding multiple comparison-corrected time windows of above-chance accuracy. This procedure allowed us to assess whether the decodability of one brain with a classifier trained on another would succeed at a meaningful time window of interest (i.e., ∼50–100 ms when MMN was observed; Fig. 2A). Notably, given that decoding was carried only on one electrode that showed the strongest effect, spatial heterogeneity between monkeys was removed entirely, which effectively left only one dimension (i.e., time) in which animals could display differences.
Dynamic characterization of the scaling behavior
To characterize the scaling properties of the neural activities of all PSCs of all monkeys, we first quantified the relationship between log(power) and log(frequency) of all individual trials (trains of 11 tones of standards and deviants combined).
Continuous power spectral densities.
The power spectral density (band: 1–100 Hz) of each trial (baseline removed) for the PSCs studied was computed by applying the modified Welch periodogram method as implemented in the MATLAB pwelch() function. We used 50% overlapping Hann windows of 1.024 s, with 250 ms removed from the beginning/end of each trial to avoid edge artifacts. Subsequently, we combined all individual trials for each PSC after removing their 200 ms baseline; the resulting series had a length of 1,327,920 (Monkey Fr) and 2,655,840 samples (Monkeys Kr and Go).
The functional repertoire of neural networks hinges on the multiscale temporal organization of their coordinated activity. We aimed at obtaining a “summary” metric of the macroscopic fluctuations of neural activity that could tease apart the dynamic behavior of the extracted ECoG components. The slope of the double logarithmic relationship between the frequency (f) and the power spectrum [S(f); Fig. 5A]—the spectral exponent (β)—can be conceived as such a measure. However, although, for example, a linear composition of stationary oscillatory signals the power spectrum accurately represents the signal, for macroscopic neural activity, which is inherently nonstationary and nonlinear (Paluš, 1996), spectral analysis can be ambiguous and may yield identical results to dynamically distinct processes (Stanley et al., 1999). A scaling analysis allows quantifying regularities in data with a wide-band spectrum—without being tied to the rigidity of an oscillatory/“clock-like” code—by using the concept of self-affinity (loosely, self-similarity; Fig. 2D). Self-affine signals are characterized by having a statistically similar structure at different levels of magnification (i.e., a level of persistence or “memory”). Persistence manifests in a signal as “random noise” superposed on a background with many cycles. These cycles are, however, not periodic (i.e., they cannot be extrapolated with time; Mandelbrot, 1983). The degree of persistence is quantified by the scaling exponent, higher persistence signifies that the signal is less irregular at all scales than the ordinary Brownian motion (uncorrelated noise; Mandelbrot, 1983). Thus, lower values of the exponent are synonymous with higher irregularity and unpredictability. Physiologic signals often display different degrees of persistence depending on the range of timescales considered. In the PSC data, this was apparent from the spectrum analysis, which revealed changing points (crossovers) in the decay of power with frequency (i.e., the decay was only piecewise linear), suggesting that the spectral densities of the PSC were not purely self-affine. Furthermore, neural signals can be more heterogeneous. The degree of persistence may differ with the moment of the time series analyzed (Kantelhardt, 2011; i.e., while traditional scaling analysis quantifies the scaling relationship in second-order statistics of a time series [(mono)fractality], a generalized analysis looks for scaling relationships in qth-order statistics, where q can take any arbitrary value. Quantifying the persistence of multiple interwoven fractal subsets (scalings) within the signal is commonly referred to as multifractal analysis. A signal is multifractal if it requires different scaling exponents, dependent on the order of statistics studied, for its description; conversely, the signal is monofractal if its scaling behavior is equal for all the moments analyzed (Fig. 2E, bottom).
Multiscale multifractal analysis.
Since our goal was to study and contrast the multifractal properties of the different ECoG components and our preliminary spectral analysis revealed the location of crossovers (changes of scaling) varied across individuals and PSCs, a predefinition of the timescales of interest was impossible. Thus, we used a method called MMA (Fig. 2E; Gierałtowski et al., 2012), a data-driven scaling analysis robust to nonstationarity. MMA is an extension of the detrended fluctuation analysis (DFA; Peng et al., 1995), an established method to quantify the monofractal scaling behavior of nonstationary signals, robust to some extrinsic trends (Hu et al., 2001). DFA is essentially a modified root mean square analysis of a random walk (Peng et al., 1995). Briefly, for a given time series,
The multifractal DFA extension (Kantelhardt et al., 2002) permits assessment of multifractality in the signals through generalizing the average fluctuation function F(s) in Equation 1, in particular, the variance
Equation 1 corresponds to the case q = 2 and assesses how the SD changes with scale. This situation corresponds to the standard DFA, which is equivalent to the Hurst (1951) exponent for stationary random processes. By varying q, we can study shorter and longer fluctuations i.e., how the qth root of the qth-order variance changes with scale (s). Typically,
We applied MMA to the PSCs of the three marmosets within a range of
Multifractal and singularity spectrum.
After obtaining h(q,s), we also computed the corresponding singularity spectrum. This representation yields a quantification of the different scaling relationships present (quantified by h above) paired with their respective proportion (density or dimension) in the signal. The scaling relationships are designated singularity strengths (α, the Lipschitz–Hölder exponent) and their densities,
One can connect the exponents obtained with MMA with the mass exponent
The exponent
There are known caveats in the computation of
We aimed at obtaining several indices (Fig. 6A, graphic overview). It is noteworthy that the singularity spectrum typically has a convex upward shape and its left branch (right branch) corresponds to the
A narrow width implies the signal is monofractal (self-affine), the length of the width gives an indication of the deviation from monofractality (i.e., of the degree of multifractality; Ihlen, 2012). Remarkably,
All parameters obtained from the singularity spectrum were computed for all marmosets and different timescales for each of the neural components, and are presented individually or are averaged across scales.
Surrogate data.
We created 50 shuffled surrogates by randomly permuting in temporal order the samples of the original time series of the PSCs of each marmoset. If the shuffling procedure yields time series exhibiting simple random behavior (
Statistics
For ERP MMN (Fig. 2A), pairwise comparisons were used by comparing the mean amplitude of pairs of adjacent standard (i.e., the last tone of the N train) and deviant (i.e., the first tone of the N + 1 train) stimuli. Similarly, for the spectrally decoupled time series (Fig. 3F,M,T), we performed separate repeated-measures ANOVA (RANOVA) of mean amplitude for each monkey between PSC (Broadband, Rhythmic 1, Rhythmic 2) and stimulus (standards, deviants), using Bonferroni correction for post hoc comparisons. Mean amplitude was always calculated over ±40 ms interval around the maximum peak detected after averaging all stimulus/component categories. In the case of the MMA (Fig. 5E), Hurst exponents were compared using RANOVA for each monkey among PSCs (Broadband, Rhythmic 1, and Rhythmic 2), and post hoc comparisons were Bonferroni corrected. Statistical analyses were performed using open-source statistical software jamovi (version 0.9; https://www.jamovi.org).
Results
Using epidurally implanted electrodes, we recorded ECoGs from three awake common marmosets, who passively listened to the stream of varying tones (Fig. 1A). By contrasting neural responses to standard and deviant tones, which were physically matched across expectancy conditions, we first identified an MMN deviance response from the auditory cortex (Fig. 2A). Afterward, we decomposed the raw LFP signal into broadband and rhythmic spectral components (Miller et al., 2009a,b, 2017; Miller, 2019; Fig. 2C). Spectral decomposition allowed us to assess whether MMN is driven by the broadband rather than oscillatory components of the LFP signal. In the following, we report a single-trial analysis that was conducted separately for each monkey (referred to as Fr, Go, and Kr), using electrodes located in the auditory cortex (Fig. 1G–I).
Experimental design and ECoG electrode arrays. A, Using a roving oddball paradigm, 20 different single tones were presented in the trains of 3, 5, or 11 identical stimuli. Any two subsequent trains consisted of different tones. This way, while the adjacent standard (depicted in black) and deviant (depicted in green) tones deviated in the frequency because of the transition between the trains, the two expectancy conditions were physically matched, as the first and the last tones of the same train were treated as deviant and standard tones in the analysis of the adjacent stimuli pairs. B, The complete 32-electrode array and connector viewed from the front. G, Ground electrode; R, reference electrode. C, A fitting example of the 32-electrode array on a model brain. D, The complete 64-electrode array and connector viewed from the front. E, A computed tomography image of the implanted 62-electrode array registered to the MRI of monkey Kr (two electrodes were cut during implantation). F, An enlarged view of the temporal area of the marmoset, showing the core (A1, R, and RT) and belt areas. The borders between each auditory area are estimated by overlaying the Common Marmoset Brain Atlas (http://brainatlas.brain.riken.jp/marmoset/modules/xoonips/listitem.php?index_id=66) on the standard brain. A1, Primary auditory cortex; R, area R (rostral auditory cortex); RT, area RT (rostrotemporal auditory cortex). G, Locations of the 32 electrodes in the monkey Fr. The red circle indicates the electrode used for the spectral decomposition and MMA analyses. H, Locations of the 64 electrodes in the monkey Go. The red circle indicates the electrode used for the spectral decomposition and MMA analyses. I, Locations of the 62 electrodes in the monkey Kr. The red circle indicates the electrode used for the spectral decomposition and MMA analyses.
ECoG data analyses. A, Time courses of ERP waveforms of the standard (black) and deviant (green) stimuli conditions. The 0 ms time point indicates the onset of a given tone. Error shades represent the SEM, calculated across all trials at each time point. Each subplot represents a different monkey. B, Time–frequency charts showing a spectral power response to the auditory stimuli, expressed as the difference between the standard and the deviant tones. The 0 ms time point indicates the onset of tones. C, Spectral decoupling. Temporally adjacent raw LFP segments of the standard tone (i.e., the last stimulus of the previous train) and the deviant tone (i.e., the first stimulus of the subsequent train) were extracted for the spectral PCA. First, the fast Fourier transform was used to calculate log power (1–250 Hz) of the raw LFP signal, which was afterward normalized across all trials within a given expectancy condition. Normalized spectral snapshots were input into spectral PCA, which separated broadband and rhythmic components. Principal spectral components were reconstructed back to the time series for the subsequent contrast between the expected and unexpected stimuli conditions. D, Two views of regularity in temporal variability. Some simple machines display mechanistic behavior that can be described by decomposing their output variables into multiple frequency-dependent oscillators or combinations of these (spectral analysis). By contrast, brains display complex behavior with output variables that appear erratic, intermittent, and nonstationary, and are the result of an inextricable interdependence of processes at many temporal scales. In this case, a measure that is scale-free and a “summary” of the whole activity is more adequate to characterize the regularities of the signal. This can be thought of as obtaining the spectrum of the different fractal/scaling exponents/singularities hidden in the signal (multifractal analysis). E, Overview of MMA (see Results, Materials and Methods).
Auditory ERP in the raw LFP signal
First, we confirmed that perturbation of auditory cortex with a deviant tone compared with a preceding standard tone increased the mean amplitude of auditory evoked potentials in the MMN time window (Fr: t(1,719) = −7.37, p < 0.001, Cohen's d = 0.275; Go: t(1,1439) = −4.60, p < 0.001, Cohen's d = 0.121; Kr: t(1,1439) = −9.27, p < 0.001, Cohen's d = 0.244; Fig. 2A). The latency of ERP peaks (58–66 ms) was consistent previous MMN studies of nonhuman primates (Javitt et al., 1992; Komatsu et al., 2015).
Auditory evoked responses reconstructed with broadband and rhythmic components
Aiming to differentiate the broadband component of LFP signal from rhythmic sources, we conducted spectral principal component analysis that decouples the PSD into components reflecting different underlying neural dynamics (Fig. 2C; see Materials and Methods). Using this technique, a broadband component can be identified by a uniform power increase (i.e., a component without clear peaks in the PSD) across a wide range of frequencies (Fig. 3A,H,O, red lines). In addition to broadband spectral component, the technique also reveals a diverse set of narrow-band oscillatory components, revealed by peaks in the PSD (Fig. 3A,H,O, blue and black lines).
Event-related broadband response of stimulus expectancy. A, H, O, Element magnitude of the major PSCs in the frequency domain (1–240 Hz). In this and other subplots, the Broadband PSC is depicted in red, the Rhythmic 1 PSC (alpha) in blue, and the Rhythmic 2 PSC (delta) in black. B, I, P, A narrow window of back-reconstructed time series of the broadband and rhythmic PSCs, locked to the onset of tones (0 ms). Standard and deviant stimuli are averaged together. C–E, J–L, Q–S, Back-reconstructed time series of the Broadband and Rhythmic PSCs along a sequence of 11 identical tones. The 0 ms value indicates the onset of the deviant tone. F, M, T, ANOVA results of the stimulus expectancy (standard, deviant) and the spectral component (Broadband, Rhythmic 1, Rhythmic 2) contrast. Significant main effects were observed for the PSC (Fr: F(2,1438) = 341.70, p < 0.001, η2 = 0.322; Kr: F(2,2878) = 113.00, p < 0.001, η2 = 0.073; Go: F(2,2878) = 78.60, p < 0.001, η2 = 0.052) and the stimulus expectancy (Fr: F(1,719) = 14.1, p < 0.001, η2 = 0.01; Kr: F(1,1439) = 23.60, p < 0.001, η2 = 0.016; Go: F(1,1439) = 4.81, p < 0.029, η2 = 0.003) factors, and the interaction between the PSC and the stimulus expectancy (Fr: F(2,1438) = 17.20, p < 0.001, η2 = 0.02; Kr: F(2,2878) = 20.80, p < 0.001, η2 = 0.014; Go: F(2,2878) = 15.49, p < 0.001, η2 = 0.011). Error bars indicate the standard error of the mean (SEM). =, The main effects; x, the interaction. G, N, U, Stimuli-locked waveforms show post hoc comparisons between the standard and deviant stimuli in the broadband and rhythmic PSCs, which revealed larger amplitude for the deviant stimuli in the Broadband PSC contrast (Fr: t = 6.96, pBc < 0.001; Kr: t = 7.84, pBc < 0.001; Go: t = 5.48, pBc < 0.001), but not in the Rhythmic 1 (Fr: t = 0.378, pBc = 1.00; Kr: t = 0.612, pBc = 1.00; Go: t = 0.397, pBc = 0.99) or the Rhythmic 2 (Fr: t = 0.812, pBc = 1.00; Kr: t = −0.033, pBc = 1.00; Go: t = −1.567, pBc = 1.00) PSC contrasts. pBc = p values after Bonferroni correction. Note: For visualization purposes, 11-tone time series were smoothed with an 80 ms Gaussian envelope (SD = 80 ms) and single-tone time series were smoothed with a 2 ms Gaussian envelope (SD = 20 ms). However, all statistical analyses were conducted on nonsmoothed data.
This way, three major PSCs, one representing a Broadband component and two representing Rhythmic components, were identified from the auditory LFP signal. PSCs were highly consistent across three monkeys (Fig. 3A,H,O), matching tightly with the original depiction of spectral PCA (Miller et al., 2009b, their Fig. 1A). To assess which of these three major PSCs encode auditory prediction error response, components were back-projected to the time dimension.
We found that the Broadband PSC carried a characteristic auditory event-related broadband (ERBB) response, reminiscent of auditory ERP, compared with the relatively low amplitude of responses derived from the rhythmic PSCs with alpha (Rhythmic 1, ∼9–11 Hz) and delta (Rhythmic 2, ∼1–3 Hz) peaks. Arguably, the Rhythmic 1 component represents endogenous alpha oscillations generated in the auditory cortex (Haegens et al., 2015; Keitel and Gross, 2016; Billig et al., 2019), whereas the Rhythmic 2 component reflects temporal structure of exogenous stimulation (i.e., the stimulus-onset asynchrony of 503 ms). The ERBB response was evident in the average of individual—standard and deviant—responses (Fig. 3B,I,P) as well as along the whole sequence of 11 identical tones, and a response was diminished in the tone sequences reconstructed from the Rhythmic components (Fig. 3C–E,J–L,Q,S). Repeated-measures ANOVA of mean amplitude between the PSC (Broadband, Rhythmic 1, Rhythmic 2) and the stimulus expectancy (standard, deviant) factors revealed the main effects for the PSC and the stimulus expectancy, and the interaction between the PSC and the stimulus expectancy factors (Fig. 3F,M,T). Post hoc comparisons showed that the ERBB response locked to the deviant tones had a larger amplitude compared with the ERBB response locked to the standard tones in the Broadband PSC contrast, but not in the Rhythmic 1 PSC or the Rhythmic 2 PSC contrasts (Fig. 3G,N,U). We thus conclude that the MMN response recorded by the ECoG of the auditory cortex is driven by broadband rather than rhythmic components of the LFP signal.
Cross-individual decoding of stimulus expectancy with broadband and rhythmic components
While the single-subject results of ERBB response were highly consistent across all three monkeys (Fig. 3), we wanted to establish whether the broadband prediction error response of an individual monkey can be extrapolated to other individuals of the same species. This would indicate that the prediction error information generated in the auditory cortex is implemented similarly across monkeys. We thus assessed the cross-individual generalizability of the ERBB response by decoding the stimuli expectancy using the Broadband and Rhythmic PSCs. Using all trials of a respective PSC of one monkey, we trained an LDA classifier to learn stimulus categories (standard vs deviant) in the auditory cortex electrode (Fig. 1G–I). Afterward, we decoded stimulus categories using the same PSC in a different monkey. Using Broadband PCS, we obtained significant decodability in all six pairs of comparisons (i.e., cross-individual decoding among three monkeys; Fig. 4A). The time windows of significant decoding above chance level (50% AUC) were consistent with MMN and BRBB responses (Figs. 2A, 3G,N,U). Contrary to this, no significant cross-individual decodability was observed using Rhythmic 1 and 2 PCSs (Fig. 4B). These findings confirm the cross-individual generalizability of the broadband PSC encoding of stimulus expectancy.
Decoding of stimulus expectancy across monkeys with Broadband and Rhythmic components. Classification of stimulus expectancy conditions (standard, deviant) was conducted in one of the monkeys (plotted in green). Afterward, the LDA classifier was tested on the other two monkeys (plotted in purple). Time points of FDR-corrected significant decoding (AUC) of stimulus categories above are depicted in gray (for details, see Material and Methods). A, Decoding was successful in all six pairs using Broadband PSC. B, Decoding did not exceed the chance level using the Rhythmic 1 and Rhythmic 2 PSCs.
Multiscale multifractal analysis of broadband and rhythmic neural components
We hypothesized that during the evoked response, the multiscale temporal organization of the broadband differs from that of the rhythmic components. In particular, we sought to characterize and compare the scale-free temporal properties of the segregated neural components. These properties, both in terms of uniform scaling and intermittency (scale-dependent changes in scaling), are associated with perceptual and cognitive processing (Werner, 2010; He, 2014; Papo, 2014). We further hypothesized that the broadband component—the neural signal subserving oddball detection—has a more stochastic multiscale temporal organization that allows greater dynamic flexibility. The scale-free nature of the neuronal population firing rate, manifested in the broadband PSC (Manning et al., 2009; Miller, 2010), is usually estimated by determining the slope of the log-log function of PSD (power vs frequency), also referred to as 1/f (fractal) scaling. However, often the PSD is not characterized by a single exponent and may show a scale dependence (Miller, 2010; Chaudhuri et al., 2018) and/or different scaling, depending on the statistical moment, and hence exhibit multifractality (Nagy et al., 2017). Indeed, the single-trial auditory responses (standards and deviants) revealed a piecewise linear decay of power with frequency in each marmoset (Fig. 5A), suggesting that the dynamics of the underlying processes may have scale-free properties but also a heterogeneous scaling dependent on frequency (timescale). This is noticeable by the different slopes that characterize the 1/f-like PSD depending on the frequency range (Fig. 5A), precluding the fitting of a unique line to estimate the slope across the whole spectrum. Thus, to fully characterize the scale-free properties of the three components, we sought to test for the presence of scale-dependent multifractality in the series of increments of neural activity in the marmoset auditory cortex.
Multiscale multifractal analysis of the broadband and rhythmic dynamics. A, Double-logarithmic plots of the power spectral densities of the Broadband (blue), Rhythmic 1 (pink), and Rhythmic 2 (green) components during all trials (trains of standard and deviant tones) reveal a piecewise, approximately linear decay of power with frequency. The average scaling (fractal) properties of the power spectral densities (last column and black lines in other columns) are distinct across frequencies, spectral components, and marmosets. Notably, while spectral power densities are plotted in the current figure, weights of spectral PCA across frequencies are plotted in Figure 3A,H,O. B, MMA method. Top, Log-log plots of the fluctuation functions
Multifractality requires the presence of different scaling exponents (h) of different moments of the fluctuations (q) over a wide range of timescales (s; Kantelhardt et al., 2002). To quantify the multifractal behavior, we applied MMA (Gierałtowski et al., 2012; see Materials and Methods; Figs. 2E, 5B). In a nutshell, the raw fluctuations of the spectral components (Fig. 2C,E) were integrated (Fig. 2E, profile) and split into windows of 10–600 ms (Fig. 1A, ∼ISI) from which the quadratic trend (Fig. 2E, colored lines, profile) was removed. The scaling exponent [Hq(s)] characterizes the relationship between the average fluctuation of the integrated and detrended signal associated with a given moment (q) and the window of observation [scale (s); Figs. 2E, 5B]. For q = 2, SD changes with scale are quantified; for q < 0, short-term variability are quantified; and for q > 0, longer fluctuations are quantified (see Eq. 2; see also Materials and Methods). Hq(s) for a range of q and s values is represented by a surface; if the surface is approximately flat, the signal is monofractal; and if Hq(s) varies substantially with q, the signal is multifractal. Using MMA, we found that all three PSCs show considerable variability in their h(q,s) values (Hurst surfaces; Fig. 5C). The Rhythmic components displayed similar surfaces, distinct from the nonlinear profile of h across q
To determine whether the multifractality, depicted in the Hurst surfaces (Fig. 5C), is caused by the temporal correlations of the signal distribution, we created a distribution of shuffled surrogates (i.e., copies of the original data with identical mean, variance, and histogram distribution but no temporal structure). While the mean Hurst surfaces of the surrogate distribution showed for all monkeys a decrease in multifractality (p < 0.001; Fig. 5C), the averaged Hurst exponent values indicated that the neural dynamics approached randomness (h = 0.5) for all monkeys (Fig. 5E). Therefore, the multifractality is caused mostly by the temporal correlations but also by a fat-tailed probability distribution. We subsequently computed the multifractal spectrum f(α). Analogous to a Fourier analysis (i.e., the decomposition of a signal into a sum of components with fixed frequencies), f(α) can be understood as the decomposition of a signal into a set of exponents, α (Mandelbrot, 2003; Figs. 2D, 6A). Their relative presence in the signal is weighted by the f(α) function. The Broadband activity interweaved more densely sets of singularities that are less self-similar than those of the Rhythmic components and displayed a lower degree of multifractality and a more asymmetrical f(α) value (Fig. 6C,D), suggesting that its dynamics differ from simple multiplicative cascades. The shape of the multifractal spectra for the Broadband activity also displayed a right truncation (Fig. 6C), which is expected because of the leveling of h(q,s) for q < 0 (Ihlen, 2012).
Contrast of the singularity spectra-based parameters of the Broadband and Rhythmic neural activity. A, Illustration of the singularity spectrum [
To sum up, MMA revealed that the generalized scale-dependent Hurst exponent h(q,s) and the derived f(α) curves of the dynamics of Broadband and Rhythmic components show multifractality as well as marked differences of this property. Importantly, the Broadband component more closely approached a stochastic asymmetrical multifractal distribution.
Discussion
In the present study, we compared two alternative views of prediction error processing, namely whether LFP oscillatory or broadband components of neural activity encode deviant sensory stimuli. We first replicated previous research by showing that auditory MMN peaks in the auditory cortex. Afterward, we separated the LFP signal into Broadband and Rhythmic components and repeated MMN analyses separately for each component. While the main two Rhythmic components present in the data were not able to distinguish between the standard and deviant tones, the Broadband component indexed the stimuli difference in the auditory cortex, confirming our first hypothesis. The findings were highly consistent across all three marmosets, and the cross-individual decoding in the temporal domain successfully classified stimuli category (standard or deviant) when data were trained on one monkey and tested on a different one. While the decoding accuracy was modest (i.e., AUC range, 0.52–0.60), cross-individual decoding is more conservative because of individual differences such as the curvature of cortical gyri. Importantly, significant decoding was observed only with the Broadband PSC, as the decoding was unsuccessful with the Rhythmic PSCs. Our findings suggest that auditory MMN response, a classical marker of prediction error, is primarily driven by the broadband component of LFP signal.
Contrary to the oscillatory narrow-band signals, such as alpha rhythm, broadband activity recorded with ECoG primarily reflects asynchronous neuronal firing (Henrie and Shapley, 2005; Mukamel et al., 2005; Hwang and Andersen, 2011; Buzsáki et al., 2012; Li et al., 2019; Guyon et al., 2021). Likewise, our earlier modeling of power-law changes in Broadband PCS suggests that they reflect a change in mean population-averaged firing rate (Miller et al., 2009a; Miller, 2010). The association between broadband signals and neuronal spikes has been demonstrated both during passive tasks (Manning et al., 2009; Guyon et al., 2021) as well as in response to sensory stimulation (Medvedev and Kanwal, 2004; Henrie and Shapley, 2005; Hwang and Andersen, 2011; Li et al., 2019). While most of the evidence comes from visual neuroscience (Henrie and Shapley, 2005; Hwang and Andersen, 2011), a direct link between broadband activity and neuronal firing also has been demonstrated in auditory cortices (Medvedev and Kanwal, 2004; Mukamel et al., 2005), although late broadband high-frequency response in A1 may also reflect dendritic processes separable from action potentials (Leszczyński et al., 2020). Overall, given that broadband PSC primarily reflects the mean firing rate of neuronal populations, and that neuronal spiking correlates with the high-frequency LFP in the auditory cortex, including poststimulus responses, our findings suggest that prediction error responses are encoded through more flexible, dynamic, unstable broadband codes than relatively more stable oscillatory codes in the poststimulus time window.
While it has been argued that a stimulus-driven phase reset of slow frequencies in the range of delta and theta oscillations may underlie prediction error response (Fuentemilla et al., 2008; Ko et al., 2012; Arnal et al., 2015), we found that the Rhythmic 2 component with a distinctive delta peak and a considerable contribution from theta range activity (Fig. 3A,H,O) did not discriminate between standard and deviant tones. Likewise, the Rhythmic 1 component representing alpha range activity did not encode prediction error response, although several previous studies linked MMN to alpha-band power (Ko et al., 2012; MacLean and Ward, 2014). Yet, despite these discrepancies, our findings do not imply that low-frequency counterparts of neural activity do not contribute to predictive coding: long-term dependencies are relevant in sensory prediction in the auditory cortex (Rubin et al., 2016) and low-frequency neural oscillations are instrumental for enhancing (Schroeder and Lakatos, 2009) and gating information in the auditory cortex (Lakatos et al., 2013). In particular, given that we did not separate evoked and induced neural activity, it is possible that low-frequency ongoing oscillations could constrain predictive error responses during the prestimulus period. As our analyses focused solely on the poststimulus responses, follow-up studies should investigate the role of Broadband and Rhythmic PSCs in the evoked versus induced components of neural signals at both prestimulus and poststimulus time windows.
Confirming the second hypothesis, our results demonstrate that prediction error processing is subsumed by an asynchronous broadband activity with dynamic properties that are very distinct from those of the rhythmic components. Importantly, this difference is unveiled when a multiscale approach is used to characterize fluctuations with several degrees of resolution (multiple fractal hierarchies), and it is patent in the surfaces and multifractal spectrum; the difference is shown to be equivocal by simply observing the power spectral densities or doing a classic Hurst analysis. The broadband component is distinctive from the other components by its lower level of self-similarity and multifractality and also by its asymmetric multifractal spectrum. Importantly, all the broadband and rhythmic components displayed multifractal fluctuations. We interpret this finding as multifractality reflecting a generic feature of neuronal networks, with cognition operating concurrently with modulations of this property (Papo, 2014). Arguably, spike trains represent information with a multifractal temporal coding (Fetterhoff et al., 2015), and the integrated multifractal spectrum permits inference of the tuning curve of spiking activity in primates (Fayyaz et al., 2019). This could be a more effective dynamic fingerprint of how information is encoded in neuronal assemblies than the one provided by oscillatory rhythms. This hypothesis is bolstered by ideas that synchronization per se only arises in collective states where no new information can be created. In contrast, adaptive behavior emerges from more subtle forms of coordination (e.g., through the metastability or asynchronous coupling of spatiotemporal patterns of neural activity; Friston, 2000; Tognoli and Kelso, 2014). The multifractality present in the recordings reveals how the macroscopic neural dynamics are intermittent and its spectral density changes with time, which has been hypothesized to be a facet of temporal metastability (Friston, 1997; Tognoli and Kelso, 2014); at the core of metastability is the broken symmetry of spatiotemporal patterns (Kelso, 1995), which was present only in the broadband activity. Thus, it is perhaps this facet of multifractality that is relevant. In fact, the more asymmetrical multifractal spectra of the broadband activity suggest that this feature may be a proxy of a dynamic regime that allows the breakdown of symmetry, which is characteristic of systems that can perceptibly or meaningfully react to afferent inputs (Freeman and Vitiello, 2006).
Furthermore, the prediction error processing by neural assemblies in the auditory cortex is sustained by an irregular broadband component with certain dynamic behavior. Specifically, tone and aftertone processing are characterized by small fluctuations lying in a tight range of the nonergodic dynamic regime (h just above 1), which has been proposed as an explanation for the 1/f noise of cognitive processes (Grigolini et al., 2009), and large fluctuations with apparent weak correlations. The latter finding aligns with sensory discrimination of visual stimuli by monkeys being more accurate if network activity in the interstimulus period is more desynchronized (Beaman et al., 2017). Moreover, it has been substantiated by theoretical accounts that used simulated recurrent neural networks to show that correlated presynaptic input and weakly correlated network spiking activity can coexist (Renart et al., 2010) and that asynchronous irregular states are responsive to afferent stimuli while synchronized oscillatory states are not (Zerlaut and Destexhe, 2017). Arguably, the arrhythmic components of the signal, while often discounted as “noise,” enable neural dynamics to be ready to jump from one state to another quickly. In such a case, the background neural noise can continuously reorganize itself by staying in an active predictive state that is constantly ready to change, depending on the incoming stimuli and internal prediction error responses.
Our findings were enabled by a novel approach to quantify these complex dynamics of neural systems, the brain's so-called “stochastic chaos” (Freeman et al., 2001). Future studies are anticipated to extend MMA to wider frequency ranges (>100 Hz), with a fine-grained resolution to arguably uncover the spike tuning underlying sensory-state discrimination (Fayyaz et al., 2019). The broadband prediction error response should be studied further using hierarchical auditory prediction paradigms that can discriminate sensory and top-down prediction error responses (Bekinschtein et al., 2009; Chennu et al., 2013, 2016). Developed in human studies, such paradigms have recently been successfully applied in the common marmosets (Chao et al., 2018). Furthermore, while the marmoset model of MMN is deemed successful and very stable, as indicated by cross-individual decoding, the current study should be replicated using LFP recordings in humans.
Importantly, our findings reveal a meso-scale neural dynamics of MMN and contribute to a unifying framework for the micro-level to macro-level neural mechanisms of the prediction error response. While most of the auditory MMN studies are conducted at the macro level using scalp EEG recordings or at the meso-level LFP, auditory prediction error responses also have been identified using single-neuron recordings (Ulanovsky et al., 2003, 2004; Pérez-González et al., 2005; Solomon and Kohn, 2014; Nieto-Diego and Malmierca, 2016; Parras et al., 2017). In particular, individual neurons located in the primary auditory cortex increase spiking rate following the presentation of oddball stimuli, which has been observed in different mammal species, including cat (Ulanovsky et al., 2003, 2004), rat, and mouse (Nieto-Diego and Malmierca, 2016; Parras et al., 2017). Similar responses also have been identified in subcortical neurons (Pérez-González et al., 2005; Parras et al., 2017). In particular, a subclass of neurons located in the dorsal and external cortices of the inferior colliculus of the rat respond selectively to novel auditory stimuli, while muting their response to repetitive stimuli (Pérez-González et al., 2005). A recent study of single-neuron activity recorded from different auditory centers in rats and mice suggests that prediction error response is organized hierarchically along the nonlemniscal auditory pathway, which is composed of inferior colliculi, medial geniculate bodies, and the primary auditory cortex with sensitivity to the deviant tones increasing along the pathway (Parras et al., 2017). MMN-like deviance sensitivity of firing rate increases further in the nonprimary regions of the auditory cortex (Nieto-Diego and Malmierca, 2016). How do such micro-level single-neuron responses relate to the MMN potentials recorded with ECoG and/or EEG? Are there different neuronal mechanisms at different levels of measurement, such as single-neuron spiking rate versus neuronal oscillations recorded using ECoG/EEG?
Our study indicates that increased neuronal firing rate may underlie prediction error responses not only at the micro level of single-neuron recordings, but also at the higher meso-level LFP measurements. In particular, we show that MMN prediction error response is driven by the Broadband component of the meso-level LFP signal. Given that the Broadband PSC largely reflects the stochastic neuronal firing rate, as suggested by previous modeling studies (Miller et al., 2009a; Miller, 2010), our findings are consistent with the account that auditory prediction error response is indeed encoded at a single action potential level within neuronal populations, which generate broadband signal at meso-level electrophysiology and most likely at macro-level electrophysiology. Broadband LFP activity provides indirect access to the total spiking output of neurons, as shown by a growing number of experiments and simulations (Freeman, 2004; Rasch et al., 2008; Crone et al., 2011). Thus, the reported Broadband activity in this study provides a “proxy” for investigating the neuronal mechanisms underlying auditory prediction error. As such, the mesoscopic information of the Broadband LFP component represents a crucial link between macroscopic-level EEG and the microscopic-level spiking activity of neural populations (Buzsáki et al., 2012).
How could our LFP-based broadband results be reconciled with the abundant literature on frequency-specific MMN results, mostly derived from EEG experiments that do not find broadband MMN response across all frequencies? First, the low-frequency range of broadband effects can be obscured by coincident changes in specific rhythmic phenomena (Miller, 2010). Second, too often classical frequency bands are loosely equated to specific rhythms (Lopes da Silva, 2013), and the views of collective neural network activity as oscillations lend too much emphasis to “rhythmicity” (Cole and Voytek, 2017) when, in reality, in those narrow-band analyses perhaps no characteristic frequency oscillation was present and/or may even be spurious and caused by filtering (de Cheveigné and Nelken, 2019). Third, EEG artifacts may decrease signal-to-noise ratio in certain segments of the broadband signal, in which case only relatively clean segments would survive as significant detectors of prediction error response, especially when the effect sizes are small. For instance, blink artifacts may distort neural signals in the delta and theta frequency ranges (Gasser et al., 1992), whereas muscular artifacts are likely to interfere with activity in the beta and gamma ranges (van de Velde et al., 1998). Given that unexpected stimuli modulate blinking rate (Bonneh et al., 2015) and facial expression (Reisenzein and Studtmann, 2007), an extreme case of which is a startle reflex (Brown et al., 1991), different levels of ocular and muscular EEG artifacts would be expected in response to standard and deviant tones. Furthermore, MMN amplitude and latency are modulated by drowsiness (Winter et al., 1995; Sallinen and Lyytinen, 1997; Nashida et al., 2000). Given that spontaneous fluctuations of alertness level are likely to occur during passive oddball paradigms, episodes of increased drowsiness could interfere with neural processing in the theta and alpha frequency ranges (Noreika et al., 2020a,b). Thus, depending on the experimental demands, the selection and training of participants, and the data preprocessing steps, certain segments of the broadband signal may be occluded by artifactual or irrelevant signals when contrasting standard and novel stimuli, yielding narrow-band deviance responses that in fact originate from a scale-free broadband component of the neuronal signal. The speculative role of EEG artifacts in the preclusion of broadband response should be tested in the future using simultaneous EEG and LFP recordings.
To conclude, we show that in a well studied paradigm of auditory oddballs, oscillations do not constitute a means to temporally constrain information processing of auditory prediction error. They are perhaps the tips of the iceberg, the latter being an arrhythmic broadband component with asymmetric multifractal stochastic properties at several timescales. Our article establishes the relevance of the broadband activity to encode relatively low-level auditory patterns and provides a theoretical background and empirical tools to probe which predictive values lie under the “noisy” surface in other paradigms and sensory modalities. While the present study focused exclusively on the local prediction error response in the auditory cortex, follow-up studies should investigate whether broadband or oscillatory signals underlie the top-down frontotemporal predictions of incoming stimuli (Garrido et al., 2008; Wacongne et al., 2011; Chennu et al., 2016).
Footnotes
The authors declare no competing financial interests.
This article is dedicated to the memory of Prof. Walter J. Freeman (1927–2016) whose pioneering work on neurodynamics has inspired countless meaningful insights during the execution of this project. We thank Dr. Daniel Bor, Dr. Laura Imperatori, Dr. Robin Ince, and Viktor Sadílek for contributing to valuable discussion, and Prof. Karl J. Friston for helping to sharpen theoretical aspects of the study.
- Correspondence should be addressed to Andrés Canales-Johnson at afc37{at}cam.ac.uk or Valdas Noreika atv.noreika{at}qmul.ac.uk
