Scale-Free Neural and Physiological Dynamics in Naturalistic Stimuli Processing

Abstract Neural activity recorded at multiple spatiotemporal scales is dominated by arrhythmic fluctuations without a characteristic temporal periodicity. Such activity often exhibits a 1/f-type power spectrum, in which power falls off with increasing frequency following a power-law function: P(f)∝1/fβ, which is indicative of scale-free dynamics. Two extensively studied forms of scale-free neural dynamics in the human brain are slow cortical potentials (SCPs)—the low-frequency (<5 Hz) component of brain field potentials—and the amplitude fluctuations of α oscillations, both of which have been shown to carry important functional roles. In addition, scale-free dynamics characterize normal human physiology such as heartbeat dynamics. However, the exact relationships among these scale-free neural and physiological dynamics remain unclear. We recorded simultaneous magnetoencephalography and electrocardiography in healthy subjects in the resting state and while performing a discrimination task on scale-free dynamical auditory stimuli that followed different scale-free statistics. We observed that long-range temporal correlation (captured by the power-law exponent β) in SCPs positively correlated with that of heartbeat dynamics across time within an individual and negatively correlated with that of α-amplitude fluctuations across individuals. In addition, across individuals, long-range temporal correlation of both SCP and α-oscillation amplitude predicted subjects’ discrimination performance in the auditory task, albeit through antagonistic relationships. These findings reveal interrelations among different scale-free neural and physiological dynamics and initial evidence for the involvement of scale-free neural dynamics in the processing of natural stimuli, which often exhibit scale-free dynamics.


Introduction
Many natural stimuli exhibit scale-free temporal or spatial patterns, such that no particular temporal or spatial periodicity predominates (Mandelbrot, 1999). In the spatial domain, it is well documented that natural images follow a P͑f͒ ϰ 1/f ␤ spatial power spectrum, where f is the spatial frequency (Field, 1987). In the temporal domain, scale-free dynamics are characterized by a P͑f͒ ϰ 1/f ␤ temporal power spectrum, where f is the temporal frequency and the power-law exponent ␤ captures the strength of autocorrelation in the signal over time. In dynamics with a larger ␤, trends tend to persist over longer periods of time ( Fig. 1A-C). Time-varying natural images (i.e., natural movies) typically follow a P͑f͒ ϰ 1/f ␤ -type temporal power spectrum (Dong and Atick, 1995). Loudness and pitch fluctuations of natural soundscapes, such as urban and rural environmental noise (De Coensel et al., 2003), speech, and music (Voss and Clarke, 1975), also exhibit 1/f-type temporal power spectra.
In the brain, scale-free dynamics are prominent across multiple observational levels (He, 2014) and manifest in human behavioral output (Gilden, 2001). Two well-studied forms of scale-free neural dynamics are the slow cortical potentials (SCPs; He et al., 2010) and amplitude fluctuations of brain oscillations . The SCPs are the low-frequency (Ͻ5 Hz) component of broadband field potentials that exhibit a 1/f-type power spectrum. Changes in arousal state and task performance alter the power-law exponent in the SCP range recorded by electrocorticography in humans (He et al., 2010). In addition, the SCP correlates with functional MRI (fMRI) signals in both spontaneous fluctuations and stimulus-driven responses (He et al., 2008;He and Raichle, 2009;Kahn et al., 2013;Pan et al., 2013), and the spontaneous fMRI signals exhibit prominent scale-free dynamics (He, 2011).
A parallel line of research has established that amplitude fluctuations of brain ␣ oscillations also follow scalefree dynamics . Its power-law exponent has been shown to vary with task performance (Linkenkaer-Hansen et al., 2004), have a genetic contribution , increase during development (Smit et al., 2011), and differ between patients with Alzheimer's disease and controls (Montez et al., 2009). Moreover, the power-law exponent of ␣-oscillation amplitude predicts that of behavioral fluctuations across individuals (Palva et al., 2013;Smit et al., 2013), indicating a link between neural and behavioral long-range temporal correlation.
Although existing data suggest that SCP phase modulates ␣-oscillation amplitude (He, 2014), whether scalefree dynamics in the SCP and ␣-amplitude fluctuations are related remains unknown. In addition, a rich literature establishes that the healthy human heart is characterized by scale-free dynamics in its interbeat interval, and failure thereof accompanies heart disease or aging (Goldberger et al., 2002). Currently, the literature remains mixed about whether and how scale-free dynamics in the SCP and ␣ amplitude interact with heartbeat dynamics (Palva et al., 2013;Zhigalov et al., 2015). This question is especially intriguing in light of recent data showing that afferent signals from the heart interact with not only the limbic system but also perceptual and cognitive systems in the brain (Park et al., 2014).
We recorded simultaneous magnetoencephalography (MEG) and electrocardiography (ECG) in human subjects in the resting state and while listening to auditory tone sequences whose pitch fluctuations constituted scalefree dynamics with varying degrees of autocorrelation. These auditory stimuli captured second-order statistics in natural stimuli, since their temporal power spectra followed a power-law distribution. Subjects discriminated tone sequences with different degrees of autocorrelation, as captured by the power-law exponent ␤. We investigated the relationships among scale-free neural and heartbeat dynamics across sensors and individuals and over time within an individual. We further examined whether the degree of autocorrelation in scale-free auditory stimuli modulated scale-free neural and heartbeat dynamics, and whether scale-free neural dynamics predicted an individual's ability to discriminate stimuli with different degrees of autocorrelation. We hypothesized, first, that scale-free dynamics in the SCP, ␣-amplitude fluctuations, and heartbeat dynamics are interrelated, and second, that scale-free neural dynamics are involved in the processing of dynamic, scale-free stimuli.

Subjects
The experiment was approved by the Institutional Review Board of the National Institute of Neurological Disorders and Stroke. All subjects were right-handed and neurologically healthy with normal hearing. Nineteen subjects between 19 and 30 years old (mean age 24.7; 12 females) participated in a ϳ3-h long MEG session with simultaneous ECG recording. Two subjects did not have ECG data because one subject was tested before the ECG recording was implemented, and the other subject's ECG data was too noisy to reliably extract R peaks. All subjects provided written informed consent.

Stimulus creation
We created auditory tone sequences whose pitch fluctuations had five levels of autocorrelation strength, spanning from fractional Gaussian noise (fGn) to fractional Brownian motion (fBm). We used a circulant embedding algorithm (Helgason et al., 2011) to create fGn time series with Hurst exponents of 0.5, 0.6, and 0.75 (corresponding to power-law exponent ␤ ϭ 0, 0.2, and 0.5, respectively, where ␤ ϭ 2H -1), as well as fBm time series with Hurst exponents of 0.25 and 0.5 (corresponding to power-law exponent ␤ ϭ 1.5 and 2, where ␤ ϭ 2H ϩ 1). We created six unique 600-element long series for each level of ␤, where each element x j of x ␤,i is taken to represent the pitch of the jth tone in the sequence. We verified that each synthesized x ␤,i indeed had the desired ␤ by computing the autocorrelation function and performing power spectral analysis for each sequence ( Fig. 1A and 1B, respectively). Each auditory sequence i was unique to ensure that subjects would respond to statistical properties of the sequence rather than memorizing particular features.
After verifying the autocorrelation properties of the sequences, we translated and scaled each x ␤,i so that its elements ranged from log(220) to log(880) and discretized the series such that each element took on one of 25 values evenly spaced on the log scale. Let us refer to the scaled, translated, and discretized series as p ␤,i . Each p ␤,i thus represents a time series of tone pitches, where pitch varies in semitone steps between 220 and 880 Hz and exhibits autocorrelation prescribed by ␤. This range of pitch values was chosen to span the isoloudness region of human hearing (i.e., with identical amplitude, subjective loudness varies minimally with changing pitch in this range).
To produce an auditory stimulus for each tone sequence, we first computed the time series of tone frequencies as f _␤, i ϭ exp͑p _␤, i ͒, in Hz. For each f j ʦf _␤, i , we constructed a sinusoidal sound wave of duration 300 ms at a 44,100-Hz sampling frequency (thus yielding 13,230 samples for each tone) according to y _j,s ϭ Acos͓2f j ͑s /SR͒ ϩ j ͔, 1 Յ s Յ 13, 230 , where s denotes sample number, f j denotes tone frequency, SR denotes the sampling frequency of 44,100 Hz, and A ϭ 1. The amplitude of the tones was kept identical throughout the sequence, i.e., the tones were not amplitude-modulated. The 300-ms duration was chosen for ease of listening (Patel and Balaban, 2000). Because Sequences of two different overall range were presented (two left columns, large fluctuation range; two right columns, small range), to demonstrate that trend strength (i.e., autocorrelation) is independent from overall range. D, Group-average (n ϭ 19) conditional probability of behavioral response given stimulus ␤, color-coded by the proportion of response ␤ at each stimulus-␤ level. Each row sums to 1. each tone was 300 ms and each sequence contained 600 tones, each sequence had a total duration of 180 s.
Cosine waves y _j, s for each tone j were concatenated, such that there was no silence period between consecutive tones. For the first tone j ϭ 1, the phase was set to 0. For subsequent tones j Ͼ 1, the absolute value of was set to the arccosine of the final sample of tone j -1. The sign of for tone j was set such that the first-order time derivative was continuous across the transition from the end of tone j -1 to the start of tone j. This ensured that there was a smooth and continuous transition between cosine waves in the junctions where tone frequency changed.
Auditory sequences were presented using the Psych-PortAudio function of the Psychophysics Toolbox (Brainard, 1997) in Matlab (Mathworks, Natick, MA). The audio was delivered through specialized ear tubes that were custom fitted to work within the MEG scanner. We used Etymotic ER-3 Insert Headphones, in which the frequency response is flat to 5 kHz. The plastic tubing from the transducer to the earpiece had a speed-of-sound delay of around 10 ms, which was corrected in MEG data analyses.

Experimental design
After presentation of each auditory tone sequence (3 min long), subjects were asked to judge its "trend strength" on a scale of 0 to 4. Higher "trend strength" was explained to the subjects as the tendency for a trend to persist over longer period of time, which captured the strength of autocorrelation in the time series. Ascending entries on the "trend strength" scale corresponded to ascending levels of stimulus ␤. Visual performance feedback was presented after every trial to assist subjects in learning how to use the scale to accurately characterize stimulus ␤. The feedback indicated what trend strength rating had been entered by the subject, what the true trend strength of the sequence was, and whether the subject's trend strength rating was correct, close to correct (off by one level of trend strength), or incorrect (off by two or more levels of trend strength). To ensure adequate task performance, all subjects were trained during an initial behavioral session that took place at least a few days previously in which shortened versions of stimuli sequences were used. Subjects who performed adequately (behavioral Ͼ 0.2, see below) were invited back for the MEG/ECG experiment (82.4% of all subjects tested). At the start of the initial behavioral session as well as the main MEG/ECG experiment, subjects were visually presented with shortened examples of each level of trend strength, as exemplified in Fig. 1C, which shows that trend strength was independent of overall range. Subjects also completed one practice block at the start of the main experiment.
Stimulus sequences were initially created with the intention of having ␤ values of 0, 0.5, 1.01, 1.5, and 2, and feedback presented to the subject was consistent with this scheme. However, because of an aliasing artifact, stimuli with an intended ␤ ϭ 1.01 had an empirical ␤ closer to 0.2 (Fig. 1A,B), and so we treat these stimuli as having ␤ ϭ 0.2 in all analyses. Because of this complication, the implicit mapping between trend strength rating and stimulus ␤ communicated to subjects via visual feedback was erroneous for the ␤ ϭ 0.2 case. Nonetheless, subjects demonstrated an ability to accurately detect that pitch autocorrelation was weaker for stimulus ␤ ϭ 0.2 than for ␤ ϭ 0.5 (Fig. 1D). This suggests that subjects were able to accurately classify stimulus ␤ in spite of the feedback error and further justifies our treatment of these stimuli as ␤ ϭ 0.2 rather than ␤ ϭ 1.01. Additionally, we verified that no analysis presented herein yielded different statistical inference if using the originally intended value of stimulus ␤ ϭ 1.01 rather than the empirically derived value of ␤ ϭ 0.2.
We also included a rest condition (3-min long trials) in which no auditory stimulus was presented. Throughout the experiment-during both auditory task and resting state-subjects were asked to keep their gaze fixed on a cross presented at the center of the screen to minimize eye movement. Each subject completed 36 trials in total, which included six trials per condition (five stimulus ␤ levels plus rest condition). The stimulus sequences differed across trials within the same ␤ level. The trials were grouped in blocks of three, resulting in 12 blocks in total. Stimulus ␤ was randomized across blocks, whereas rest trials were evenly dispersed throughout the experiment (always presented as the second trial of even numbered blocks). The head position of the subject was measured with respect to the MEG sensor array using coils placed on the left and right preauricular points and the nasion. Before each block subsequent to the first block, the subject self-corrected the head position to the same position recorded at the start of the first block using a custom visual-feedback program to minimize head displacement across the experiment. Video monitoring of the subject during the experiment ensured that subjects stayed alert and did not close their eyes for extended periods of time. After excluding trials that had failure of behavioral response or drowsiness as shown by eye closure, 14 subjects had all 36 trials, two subjects had 35 trials, two subjects had 33 trials, and one subject had 30 trials.

Data acquisition and preprocessing
Experiments were conducted in a whole-head 275channel CTF MEG scanner (VSM MedTech, Coquitlam, BC, Canada). MEG data were collected with a sampling rate of 600 Hz and an anti-aliasing filter at Ͻ150 Hz. Analyses were performed on 271 sensors after excluding four malfunctioning sensors. The Fieldtrip package (Oostenveld et al., 2011) implemented in Matlab was used for data preprocessing, and analyses were conducted using Fieldtrip and custom-written code. We used independent component analysis to remove artifacts related to eye blinks, eye movements, heartbeat, breathing, and slow movement drift. Empty-room recording was collected in a prior experiment to verify that instrument noise was orders of magnitude lower than the signal we analyzed.

Measuring scale-free parameters
Various scale-invariance measures are mathematically related and thus can be reasonably compared (Eke et al., 2002). Power-law exponents of SCP and ␣-oscillation amplitude fluctuations were estimated using the common power spectral analysis. To estimate the power-law exponent ␤ of the SCP, a fast Fourier transform was applied to the MEG signal from each sensor in each trial (3 min long) to compute its power spectrum. The power spectrum was plotted in double-logarithmic scale ( Fig. 2A). Because a power spectrum following P͑f͒ ϰ 1/f ␤ can be rewritten as log͓P͑f͔͒ Ϫ ␤log͑f͒, the negative slope in the log-log plot provides a convenient measure of the powerlaw exponent ␤. In line with previous studies (He et al., 2008(He et al., , 2010, the SCP ␤ was estimated in the range of 0.005-5 Hz (power spectra were calculated on each trial lasting 180 s; thus the lowest frequency visible was 0.0056 Hz).
To estimate the power-law exponent ␤ of ␣-oscillation amplitude fluctuations, we first extracted ␣ oscillations from continuous MEG signal in each trial using a thirdorder Butterworth filter between 6.7 and 13.3 Hz and computed its instantaneous amplitude envelope by applying the Hilbert transform. A fast Fourier transform was then applied to the ␣-oscillation amplitude fluctuation to create a power spectrum for each sensor in each trial. The ␣-oscillation amplitude ␤ was estimated from the log-log plot of power spectrum using the 0.1-to 1-Hz range (Fig.  2B).
Following previously established methodology, we applied detrended fluctuation analysis (DFA) on heartbeat dynamics (Goldberger et al., 2002;Hardstone et al., 2012). This computation was applied to the time variation in interbeat interval, measured as the interval in seconds between adjacent R peaks ( Fig. 2C, left). DFA analysis was conducted as follows. First, the interbeat interval time series from each trial was integrated, and the mean was subtracted. We then estimated the local trend in nonoverlapping windows of equal length using a least-squares fit and determined the fluctuation as variance on the local trend for a given window. Five different window lengths were used: 4, 8, 16, 32, and 64. The log-log plot of mean fluctuation (F) against window length (l) constitutes the DFA plot (Fig. 2C, right), and its slope estimates the DFA exponent ␣, following the relation: F͑l͒ ϰ l ␣ , or log͓F͑l͔͒ ϰ ␣log͑l͒. Theoretically, the power-law exponent ␤ is related to the DFA exponent ␣, following ␤ ϭ 2␣ Ϫ 1 (Eke et al., 2002).

Interrelations between scale-free neural and heartbeat dynamics
To qualitatively assess the relationship between SCP ␤ and ␣-oscillation amplitude ␤ across sensors, we first generated grand-average scalp topography for SCP ␤ and ␣-amplitude ␤ across subjects and trials ( Fig. 2A,B, Figure 2. Characterization of neural and physiological dynamics. A, Left, MEG signal power spectra from an example subject in each task condition (averaged across all sensors and trials within a condition). Middle, example MEG signal power spectrum from a single trial in one subject. The red lines indicate the frequency range for extracting SCP power-law exponent ␤: 0.005-5 Hz. Right, grand-average topographic map of SCP ␤ across the scalp. B, As in A, but for ␣-oscillation amplitude. ␣ Oscillation was filtered in the 6.7-to 13.3-Hz range, and its amplitude time series was extracted using the Hilbert transform. ␣-Amplitude ␤ was extracted in the 0.1-to 1-Hz range (red lines in the middle panel). C, For analysis on heartbeat dynamics, the interbeat interval was calculated as the difference in time (s) between adjacent R-peaks in the ECG recording (left panel). The interbeat interval time series was subjected to DFA to extract the DFA exponent ␣, which describes the power-law relationship between fluctuation magnitude and the length of observation in scale-free dynamics (right panel). right). No formal statistical test was performed for the relationship between SCP ␤ and ␣-amplitude ␤ across all sensors due to the difficulty in accurately accounting for the degrees of freedom associated with the sensors. However, a scatterplot of SCP ␤ and ␣-amplitude ␤ across sensors is included for descriptive purposes (Fig.  3A).
To investigate the relationship between SCP ␤ and ␣-amplitude ␤ across subjects, we first defined two clusters of sensors based on the topographic map of SCP ␤ (Fig. 3C) and then averaged SCP ␤ and ␣-amplitude ␤, respectively, across sensors within each cluster. Pearson's correlation was calculated between SCP ␤ and ␣-amplitude ␤ across subjects for each of the two clusters.
To assess the relationship between any two of our three scale-free parameters (SCP ␤, ␣-amplitude ␤, or ECG ␣) over time, we calculated Pearson's correlation between them across all 36 trials within each subject. Pearson's r values were transformed into Fisher's z-values, which were subjected to a one-sample t-test across subjects. Statistical significance was assessed by a cluster-based nonparametric permutation test (Nichols and Holmes, 2002;Maris and Oostenveld, 2007). To this end, we shuffled one variable across trials for 1,000 iterations. For each iteration, Pearson's correlation between the two variables was recomputed, and the r-value was transformed into Fisher's z-value and submitted to a onesample t-test across subjects as with the original data. Clusters were defined for both the original and shuffled data as contiguous groups of sensors with p-values less than 0.05 and t-values of the same sign. Summing the t-values created a summary measure of each cluster. To build the null distribution, the absolute value of the summed t-statistic of the largest magnitude was extracted from each iteration. Finally, the absolute magnitude of the summed t-statistic in each cluster from the original data was compared to the null distribution. A cluster survived cluster-based correction if 2.5% or fewer of observations in the null distribution surpassed the absolute value of the cluster's summed t-statistic (corresponding to p Ͻ 0.05 in a two-tailed test). Correlation between MEG signal power and ␤ was tested similarly. Interrelations between scale-free neural and physiological dynamics. A, Scatterplot of SCP ␤ against ␣-amplitude ␤ across all MEG sensors (averaged across subjects and stimulus conditions). The inset shows the grand-average topographical plots of SCP ␤ and ␣-amplitude ␤, reproduced from Fig. 2. B, Top, nested-frequency analysis between SCP phase and ␣ amplitude. Groupaverage MI z-score topography plot is shown (middle), along with the phase-amplitude histogram for the two dominant clusters (see insets). Bottom, scatterplot between empirically measured SCP ␤ and simulated ␣-oscillation amplitude ␤ across all sensors. C, Left, scatterplot across subjects between SCP ␤ and ␣-amplitude ␤ in a posterior sensor cluster (see inset). Right, scatterplot across subjects between SCP ␤ and ␣-amplitude ␤ in an anterior sensor cluster (see inset). D, Correlation between SCP ␤ and SCP power across trials within an individual. Pearson correlation values were Fisher's z transformed and subjected to a group-level one-sample t-test. t Values are plotted in the topography plot, with a single cluster encompassing all sensors surviving cluster-based correction. The group-level mean and SEM of Fisher z values (averaged across all sensors) are shown in the bar graph to the right (p Ͻ 0.001, one-sample t-test across subjects). E, Similar to D, but for within-subject, across-trial correlation between SCP ␤ and ECG ␣. Two clusters of sensors survived cluster-based correction. The group-level mean and SEM of Fisher z values (averaged across all significant sensors) are shown in the bar graph to the right (p Ͻ 0.001, one-sample t-test across subjects).

SCP: ␣-oscillation nested-frequency analysis and simulation
Given previous EEG and electrocorticography findings showing a nested-frequency relationship between SCP phase and ␣-oscillation amplitude (Vanhatalo et al., 2004;Monto et al., 2008;He et al., 2010;He, 2014), a natural question is whether this nested-frequency relationship produces any correlation between SCP ␤ and ␣-amplitude ␤. To address this question, we performed simulations to reveal what kind of relationship between SCP ␤ and ␣-amplitude ␤ would be expected if it were driven entirely by the nested-frequency relationship between them.
First, we quantified the nested-frequency relationship between SCP phase and ␣-oscillation amplitude using the well-established modulation index (MI;Tort et al., 2008;He et al., 2010). We extracted the SCP phase time series by using a third-order Butterworth filter between 0.005 and 1 Hz (using a 0.005-to 5-Hz filter yielded nearly identical results) and then applying the Hilbert transform. The ␣-oscillation amplitude time series was derived as described above. SCP ␣ nested-frequency plots (for each subject, sensor, and condition) were generated by binning the SCP phase time series into 20 evenly spaced phase bins and averaging the ␣-oscillation amplitude within each phase bin. The MI was computed based on this nestedfrequency plot, which uses an inverted entropy measure to quantify its deviation from a uniform distribution. For statistical testing, the MI value was converted into an MI z-score by comparison with a null distribution generated by shuffling the phase time series using five equal-length segments, following a previously described method (He et al., 2010). A preliminary analysis suggested that stimulus condition did not modulate MI z-score, and thus different conditions were combined in the subsequent simulation.
We next used the empirically measured SCP phase series in conjunction with the nested-frequency relationship between SCP phase and ␣-oscillation amplitude at each sensor to construct simulated ␣-oscillation amplitude time series. To this end, the SCP-␣ nested-frequency plot as described above was averaged across subjects and conditions for each sensor. This distribution was smoothed using a five-bin-wide moving average. For each sample of the SCP phase time series, we used spline interpolation of the nested-frequency distribution at that sensor to determine what ␣-oscillation amplitude would be predicted by the SCP phase. Finally, we computed the ␤ of the empirical SCP time series and the simulated ␣-oscillation amplitude series, averaged these ␤ values over trials and subjects, and examined their Pearson correlation across sensors. This result reveals the contribution of the nested-frequency relationship between SCP and ␣ oscillations to the correlation between SCP ␤ and ␣-oscillation amplitude ␤.

Stimulus modulation of scale-free neural and heartbeat dynamics
We probed whether the strength of autocorrelation in the auditory stimulus (as captured by its power-law ex-ponent ␤) modulated scale-free neural and heartbeat dynamics. For each sensor in each subject, stimulus ␤ was correlated with SCP ␤, ␣-oscillation amplitude ␤, or ECG ␣ across the 30 task trials using Spearman's rank correlation. Spearman's rho values were transformed into Fisher's z-values (Fieller et al., 1957) and submitted to a one-sample t-test across subjects at each sensor. Statistical significance was established using a cluster-based nonparametric permutation test as described earlier.

Behavioral performance assessment and correlation with scale-free neural dynamics
To evaluate subjects' behavioral performance in the stimulus ␤ discrimination task, we first visualized the conditional probability of behavioral response ␤ (i.e. the ␤ corresponding to the subject's trend strength rating) at a given stimulus ␤ (Fig. 1D). The behavioral performance of each subject was assessed by Spearman's rank correlation between stimulus ␤ and response ␤ across all task trials. The Spearman's , or behavioral , captures a subject's behavioral performance for the entire experiment while allowing for some leniency in which subjects could be close to the right answer but not exact. To test whether scale-free neural dynamics predicted behavioral performance on a subject-by-subject basis, we computed Pearson's correlation between behavioral and either SCP ␤ or ␣-oscillation amplitude ␤ at each sensor across subjects.

Behavioral performance
The across-trial Spearman's correlation between stimulus ␤ and response ␤ ("behavioral ") was significant for every subject (n ϭ 19, p ranged from p Ͻ 0.0001 to 0.048 [ Table 1, line a]), indicating that all subjects could perform the stimulus ␤ discrimination task significantly above chance level. Behavioral ranged from 0.36 to 0.88, with an average value of 0.66. Fig. 1D shows the groupaverage conditional probability map of behavioral response given stimulus ␤. The accuracy of subjects' behavioral responses is reflected by the concentration of the probability distribution along the diagonal, where response ␤ is equal to stimulus ␤. Note that incorrect behavioral responses tended to be close to the correct response (i.e. off-diagonal elements that are closer to the diagonal have higher probability than those that are farther away). The Spearman's correlation coefficient thus provides a more informative and natural measure for quantifying overall behavioral performance than proportion correct. In particular, Spearman's correlation, but not the proportion of correct responses, takes into account the magnitude of response error. Interestingly, subjects were better at discriminating between stimulus ␤ in the fBm (␤ ϭ 1.5 or 2) than the fGn range (␤ ϭ 0, 0.2 or 0.5; Fig. 1D). We note that stimuli in the fBm range tended to sound more melodic than those in the fGn range.

Scale-free dynamics in neural activity and heart rate variability
Consistent with earlier reports (Dehghani et al., 2010;He, 2014), the power spectrum of MEG signals followed power-law scaling, with peaks at discrete frequencies corresponding to various brain oscillations (see Fig. 2A, left, for example power spectra from a single subject). We extracted power-law exponent ␤ for SCP from the 0.005to 5-Hz frequency range ( Fig. 2A, middle). The grand average (across stimulus conditions and subjects) topographical distribution of SCP ␤ across MEG sensors is shown in Fig. 2A (right). SCP ␤ ranged from 0.87 to 1.12 across sensors (mean 0.97), and exhibited an anteriorposterior gradient with frontal sensors displaying higher ␤, and hence longer temporal autocorrelation in the SCP. This finding is consistent with previous MEG (Dehghani et al., 2010) and fMRI (He, 2011) observations. We extracted instantaneous amplitude fluctuations of ␣ oscillations (filtered in the 6.7-to 13.3-Hz range), and computed its power spectrum (see Fig. 2B, left, for result from an example subject). In line with previous reports , ␣-oscillation amplitude fluctuations exhibit power-law scaling in the power spectrum. We extracted the power-law exponent ␤ of ␣ amplitude in the range of 0.1-1 Hz (Fig. 2B, middle). ␣-Amplitude ␤ ranged from 0.32 to 0.61 across sensors (mean 0.51) and displayed an anterior-posterior gradient opposite to that of SCP, such that posterior sensors had higher ␤, and accordingly, stronger autocorrelation (Fig.  2B, right).
Finally, following established procedures for analyzing heartbeat dynamics (Goldberger et al., 2002), we defined R-peaks from ECG recordings and constructed an interbeat interval time series, which was subjected to DFA analysis to extract the DFA exponent ␣ (Fig. 2C). Theoretically, the DFA exponent ␣ is directly related to the power-law exponent ␤ (see Materials and Methods), both of which capture the strength of autocorrelation in a time series, one using a time-domain approach (DFA exponent ␣), the other a frequency-domain approach (power-law exponent ␤). Across 17 subjects with simultaneous ECG-MEG recordings, ECG ␣ ranged from 0.57 to 1.2, with a mean of 0.82 across subjects, consistent with previous reports (Goldberger et al., 2002).

Anticorrelation between SCP and ␣-oscillation amplitude power-law exponents
We next explored the relationship between scale-free dynamics in SCP and ␣-oscillation amplitude fluctuations across the scalp and subjects. Because stimulus condition minimally modulated power-law exponent of SCP or ␣-oscillation amplitude (see below), for this analysis, we pooled data across all conditions (including five stimulus ␤ levels and resting condition).
The scalp topography of SCP ␤ and ␣-oscillation amplitude ␤ ( Fig. 2A,B, right) display opposite anterior-posterior gradients, indicating that as the strength of autocorrelation in SCP increases, the strength of ␣-oscillation amplitude fluctuations tends to decrease. To qualitatively assess this relationship, we plotted the two measures, each averaged over 19 subjects, against each other across all sensors (Fig. 3A). This revealed a negative relationship between SCP ␤ and ␣-oscillation amplitude ␤ across MEG sensors. We then assessed whether there might also be an anticorrelation between SCP ␤ and ␣-amplitude ␤ across subjects. To this end, we first defined two clusters of sensors based on the scalp topography of SCP ␤, distributed over frontal and posterior regions, with relatively high and low ␤, respectively. We then assessed across-subject correlation between SCP ␤ and ␣-oscillation amplitude ␤ for each cluster of sensors. In the posterior cluster, we observed a significant negative correlation between SCP ␤ and ␣-amplitude ␤ across subjects ( Fig. 3C; r ϭ - . Could the significant anticorrelation between SCP ␤ and ␣oscillation amplitude ␤ in the posterior cluster across subjects be driven by a relation between their respective power? Two pieces of evidence suggest that this is not the case. First, SCP power and ␣-oscillation power were found to correlate positively across subjects (r ϭ 0.91, p ϭ 8.5 ϫ 10 -8 [Table 1, line d]; but note that this finding in itself could be due to measurement variation across subjects). Second, a partial correlation analysis revealed that after controlling for the effects of SCP and ␣-oscillation power, the anticorrelation between SCP ␤ and ␣-amplitude ␤ across subjects in the posterior sensor cluster was unchanged (r ϭ -0.48, p ϭ 0.05 [ Table 1, line e]).
The above results reveal an intriguing negative relationship between SCP ␤ and ␣-amplitude ␤ across the scalp and subjects, such that stronger autocorrelation in the SCP is accompanied by weaker autocorrelation in the amplitude fluctuations of ␣ oscillations. In light of previous observations of a nested-frequency relationship between SCP phase and ␣-oscillation amplitude (Vanhatalo et al., 2004;He, 2014), these findings raise a natural question: is the anticorrelation between SCP ␤ and ␣-amplitude ␤ driven by the nested-frequency relationship between them? To test this hypothesis, we quantified the nestedfrequency pattern between SCP phase and ␣-oscillation amplitude in each MEG sensor (Fig. 3B, top) and simulated ␣-oscillation amplitude time series for each sensor in each subject, based on the sensor-specific nestedfrequency pattern and the empirically measured SCP phase time series. We then computed the power-law exponent ␤ of the simulated ␣-oscillation amplitude time series and plotted it against the measured SCP ␤ across all sensors (Fig. 3B, bottom). This simulation reveals a positive relationship between SCP ␤ and ␣-amplitude ␤, suggesting that the negative relationship observed in the empirical data cannot be explained by the nestedfrequency relationship between SCP and ␣ oscillations.
Finally, we investigated whether the amount of power in the SCP or ␣ range was related to their respective ␤ across time within an individual (see Methods). We found a robust positive correlation between SCP power and ␤: after cluster-based correction for multiple comparisons, all MEG sensors across the entire scalp demonstrated a significant positive correlation (Fig. 3D). By contrast, ␣-oscillation power had no significant correlation with its amplitude ␤ after cluster-based correction. This result is consistent with previous findings showing that SCP ␤ changes by modulating the power in the lowest frequency ranges, thereby causing a positive correlation between its power and ␤ (He et al., 2010).

Relationship between scale-free neural and physiological dynamics
We further investigated whether the strength of autocorrelation in scale-free neural and heartbeat dynamics comodulated across time within an individual. To this end, we computed correlations between SCP ␤ and the DFA exponent ␣ of heartbeat dynamics measured by ECG. This analysis revealed two significant clusters, one over the left central cortex, and the other over the right central cortex extending into frontal areas (Fig. 3E). No significant correlation was found between ␣-amplitude ␤ and ECG ␣ after correction for multiple comparisons.

Stimulus condition modulates scale-free dynamics in ␣-oscillation amplitude
Does the strength of autocorrelation in the stimulus sequence ("stimulus ␤") modulate the strength of autocorrelation within scale-free neural or physiological dynamics? To answer this question, we computed Spearman's rank correlation between stimulus ␤ and the autocorrelation parameter from neural or heartbeat dynamics (respectively, SCP ␤, ␣-oscillation amplitude ␤, and ECG ␣) across all trials during the auditory task (30 trials in total) for each subject.
We found that as stimulus ␤ increased, ␣-amplitude ␤ progressively decreased in a posterior sensor cluster overlying visual cortex at a corrected p ϭ 0.023 (Table 1, line f; Fig. 4A, top). For the sensors within this cluster, the mean ␣-amplitude ␤ across subjects in each condition is plotted in Fig. 4A (bottom). Interestingly, white noise input (␤ ϭ 0) enhanced ␣-amplitude ␤ compared to the rest (p ϭ 0.0065 [ Table 1, line g], paired t-test across subjects), and there was a trend effect of stimuli with strong autocorrelation (␤ ϭ 2) reducing ␣-amplitude ␤ compared to the rest (p ϭ 0.0774 [ Table 1, line h]). A control analysis indicated that stimulus ␤ did not significantly influence MEG signal power in the ␣ range. No significant correlation between stimulus ␤ and SCP ␤ was found, nor between stimulus ␤ and ECG ␣, suggesting that the strength of autocorrelation within SCP and heartbeat dynamics was robust to the range of scale-free stimuli used in this experiment.

Scale-free neural dynamics predicted behavioral performance
Finally, we investigated whether scale-free neural or physiological dynamics predicted a subject's behavioral performance in the auditory task. Behavioral performance was assessed by Spearman's correlation between stimulus ␤ and response ␤ as described above (behavioral ). We found that behavioral performance correlated negatively with SCP ␤ (averaged across all conditions) in a group of sensors distributed over frontocentral areas ( Fig.  4B; r ϭ -0.62, p Ͻ 0.0044 [ Table 1, line i]), and positively with ␣-amplitude ␤ in a group of sensors over left frontotemporal cortices ( Fig. 4C; r ϭ 0.64, p ϭ 0.0030 [Table 1, line j]). These results suggest that higher autocorrelation within ␣-oscillation amplitude fluctuations, and lower au-tocorrelation within the SCP, predict better discrimination performance on scale-free auditory stimuli. Importantly, neither SCP power nor ␣-oscillation power correlated with behavioral performance, suggesting that signal power was not a mediating factor between scale-free neural dynamics and behavioral performance. A control analysis further suggested that the above effects are regionally specific: ␣-amplitude ␤ in the frontocentral area (Fig. 4B

Discussion
In this study, we investigated the interrelations among scale-free dynamics in the SCP, ␣-oscillation amplitude fluctuations, and heartbeat dynamics across MEG sensors and subjects and over time within an individual. We further explored their modulation by scale-free dynamic stimuli and tested whether an individual's scale-free neural dynamics predicted the ability to tell scale-free stimuli apart based on autocorrelation property. Below, we summarize our findings in turn and discuss their implications.

Interrelations among scale-free neural and physiological dynamics
Across the scalp, a qualitative pattern emerged such that sensors exhibiting stronger autocorrelation (hence, larger power-law exponent ␤) in the SCP tended to have weaker autocorrelation in the ␣-oscillation amplitude fluctuations. In addition, SCP ␤ and ␣-amplitude ␤ were anticorrelated across subjects within a large posterior sensor cluster. This anticorrelation could not be explained by the nested-frequency coupling between SCP and ␣ oscillations, as a control analysis based on simulation suggested that the phase-amplitude coupling between SCP and ␣ oscillations contributes to a positive correlation between their power-law exponents instead. Together, these results suggest that not only do scale-free dynamics exist within both arrhythmic brain activity and amplitude fluctuations of brain oscillations, but these different scale-free neural dynamics are related and follow a systematic antagonistic pattern. Functionally, this anticorrelation may be important for preventing excessive longrange temporal correlation in the brain, such that strong autocorrelation in one type of neural signals impedes the generation of strong autocorrelation in another. Because proper brain functioning requires a balance of sufficient order and flexibility, such anticorrelation may be evidence of a negative feedback mechanism whereby selforganized brain activity is regulated across levels to avoid excessive regularity or overly random fluctuation. The current study does not address the mechanism giving rise to the anticorrelation between SCP ␤ and ␣-amplitude ␤ across subjects. In particular, it remains unknown whether these two measures have a common mechanism or different mechanisms under common influence, or alternatively, whether one measure influences the other directly or indirectly. Nonetheless, developmental and genetic contributions that have been shown to influence ␣amplitude ␤ Smit et al., 2011) indicate possible starting points for future investigations to probe.
We further observed that SCP power positively correlated with SCP ␤ across time within an individual, indicat- . Stimulus modulation of scale-free neural dynamics and prediction of behavioral performance. A, Spearman rank correlation was calculated between stimulus ␤ and ␣-amplitude ␤ across the 30 task trials for each subject, and the group average is plotted for all sensors (top left panel). A posterior sensor cluster survived cluster-based correction at p Ͻ 0.05 (top right panel). For this significant sensor cluster, ␣-amplitude ␤ averaged across sensors was plotted for each stimulus ␤ level and rest condition (bottom panel), which shows the mean and SEM across subjects. B, Left, Pearson's correlation value between behavioral performance (measured as behavioral ) and SCP ␤ (averaged across all conditions) across subjects, thresholded at a p Ͻ 0.05 level. Nonsignificant sensors are shown as a uniform green background. Right, SCP ␤ averaged across significant sensors is plotted against behavioral for all subjects. C, Same as B, but for the correlation between ␣-amplitude ␤ and behavioral performance.
ing that higher ␤ in the SCP is a result of higher power in the lowest frequencies. This resonates with a previous finding on the pattern of power spectral changes in this frequency range during task performance (He et al, 2010).
Our results reveal a novel relationship between scalefree neural and physiological dynamics, with the strengths of autocorrelation in the SCP and heartbeat dynamics (captured by SCP ␤ and ECG ␣, respectively) positively comodulating across trials within an individual. Because neither measure was influenced by stimulus condition, this relationship is due to their intrinsic fluctuations over time. By contrast, we did not observe a significant relationship between ␣-amplitude ␤ and ECG ␣; such a relationship was reported in Palva et al. (2013) but failed to be reproduced in a later study from the same group (Zhigalov et al., 2015), who reported a correlation between the scaling exponents of ␦-oscillation amplitude and heartbeat dynamics; however, they tested many frequency bands and brain systems without correcting for multiple comparisons. In addition, neither of these two previous studies directly recorded ECG, but rather used independent component analysis-extracted component from the MEG recording to substitute for heartbeat signal. Our result with direct ECG recording suggests that fluctuations in slow, arrhythmic neural activity coordinate with heart signals, although the directionality of this influence remains unknown at present. We speculate that a tight correlation between SCP ␤ and ECG ␣ may be because the time scales at which SCP and heartbeat dynamics fluctuate are comparable, both taking place on the order of many seconds (SCP, 0.2-200 s; heartbeat, 4 -64 s; compared with ␣ amplitude, 1-10 s). Together, the anticorrelation between SCP ␤ and ␣-amplitude ␤ and the positive correlation between SCP ␤ and ECG ␣ may suggest the SCP as a central link that connects scale-free neural and physiological dynamics across scales and systems. Nonetheless, neuroanatomical interpretation for the spatial distribution of sensors whose SCP ␤ correlate with ECG ␣ (Fig.  3E) would be better informed by future investigations using invasive recordings and/or source reconstruction.
More broadly, it has been shown that the brain exerts strong autonomic influence on and receives feedback from the heart (Craig, 2002;Gray et al., 2007). Previous studies suggest that scale-free heartbeat dynamics may be adaptive, with its long-range temporal correlation serving as a self-organizing mechanism for highly intricate processes that generate fluctuations across wide timescales (Ivanov et al., 1996). Indeed, highly periodic or rigid behaviors may narrow functional responsiveness, as shown by the observation that the breakdown of scalefree heart dynamics and appearance of excessive regularity often accompany pathologies such as severe congestive heart failure (Goldberger et al., 2002).

Stimulus modulation of scale-free neural dynamics
We observed a systematic modulation of ␣-oscillation amplitude dynamics by scale-free auditory stimuli, such that ␣-amplitude ␤ decreased with increasing stimulus ␤ in a posterior sensor cluster overlying occipital cortex. Our stimuli captured a range of stationary and nonstationary patterns, from fractional Gaussian noise to fractional Brownian motion. This result suggests that listening to stimuli that exhibit strong autocorrelation reduces autocorrelation in ␣ amplitude fluctuations in visual regions. A control analysis further suggested that stimulus ␤ had no effect on MEG signal power in the ␣ range. Why should an auditory task affect scale-free neural dynamics in visual regions? Although the underlying mechanisms of this phenomenon require future investigation, a speculative possibility is that higher stimulus ␤ translates into lower ␣-amplitude ␤ in visual regions owing to cross-modality interaction carried by an inhibitory pathway from auditory cortex to visual cortex (Iurilli et al., 2012).
In contrast, we did not observe a significant correlation between stimulus ␤ and SCP ␤ after cluster-based correction. This negative finding could have several reasons. It is possible that SCP reflects a backbone of brain network structure that remains unperturbed by changes in arousal state (He et al., 2008) or the range of scale-free stimuli used herein. Yet at present, we cannot rule out the possibility that the sample size in the current study was insufficient for detecting an effect in the SCP.

Prediction of behavioral performance
We found intriguing evidence suggesting that baseline characteristics of scale-free dynamics in the SCP and ␣-amplitude fluctuations predicted an individual's performance in discriminating between scale-free auditory stimuli exhibiting different levels of autocorrelation. Better performance correlated with higher ␣-amplitude ␤ and lower SCP ␤. Moreover, neither SCP nor ␣-oscillation power correlated with behavioral performance, suggesting specific behavioral relevance of scale-free parameters. Previous studies have shown that ␣-amplitude ␤ correlates with long-range temporal correlation in behavioral fluctuations across normal subjects (Palva et al., 2013;Smit et al., 2013). Yet, it is unclear whether longer or shorter autocorrelation in behavioral fluctuations is adaptive. On the other hand, discriminating natural stimuli based on their time-aggregate statistics should confer behavioral advantage in an ecologically natural environment. To our knowledge, this is the first study demonstrating that properties of scale-free neural dynamics predict behavioral performance across normal individuals.
Why should lower SCP ␤ and, conversely, higher ␣-amplitude ␤ predict better behavioral performance? One possibility is that SCP fluctuations include frequencies an order of magnitude lower than ␣-amplitude fluctuations (0.005-5 vs. 0.1-1 Hz). Thus, this pattern of result is consistent with the idea that there may be an optimal range of autocorrelation that is most conducive to performing this task: relatively weak autocorrelation in the very long time scales encompassed by the SCP, and relatively high autocorrelation in the comparatively shorter time scales spanned by ␣-amplitude fluctuations. Higher ␣-amplitude ␤ might also suggest a state closer to criticality with higher information-processing capacity (Shew et al., 2011;Poil et al., 2012). Tantalizing clues supporting the existence of an optimal range of scale-free neural dynamics exist from studies of clinical populations. For example, breakdown of long-range temporal correlation inand ␣-oscillation amplitude fluctuations has been observed in depression (Linkenkaer-Hansen et al., 2005), Alzheimer's disease (Montez et al., 2009), andschizophrenia (Nikulin et al., 2012). On the other hand, abnormally high long-range temporal correlation in ␤-band amplitude fluctuations is found in seizure-onset areas (Monto et al., 2007).
Finally, in our experiment, the auditory stimuli were constructed such that the scale-free statistic, embodied in the power-law exponent ␤, is the only difference between categories of auditory tone sequences. All other statistics, including tone duration, pitch range, sequence length, and higher-order statistics (which are random), were identical across stimulus categories (␤ ϭ [0, 0.2, 0.5, 1.5, 2]). Moreover, behavioral discrimination was carried out on power-law exponent ␤, not the specific sequence presented; this was ensured by presenting six unique sequences at each ␤ level and asking subjects to make discrimination about ␤ only. Hence, in our task, subjects' ability to discriminate different auditory tone sequences was specifically related to their ability to process the scale-free statistic of the stimuli, and our findings establish the role of scale-free brain activity in processing scale-free statistics of naturalistic stimuli. On the other hand, our results do not imply that the function of scale-free brain activity is specific to the processing of scale-free or natural stimuli. It is possible that similar correlations may be observed for tasks that do not explicitly require the evaluation of scale-free stimulus statistics. Future studies investigating performance in such tasks will delineate the functional specificity (or generality) of scale-free neural activity.
In summary, we observed novel relationships among scale-free dynamics in distinct components of neural and physiological activity, including the SCP, ␣-oscillation amplitude fluctuations, and heartbeat dynamics. We further demonstrate that scale-free neural dynamics can be systematically perturbed by scalefree dynamical stimuli that capture second-order statistics (i.e. autocorrelation or power spectrum) of the natural environment. Moreover, the baseline characteristics of scale-free neural dynamics in an individual predict their ability to discriminate scale-free dynamical stimuli based on their autocorrelation property. These results shed light on the complex interrelations among scale-free neural and physiological dynamics at different levels and how they may contribute to adaptive behavior in the natural environment.