Broadband spectral responses in visual cortex revealed by a new MEG denoising algorithm

Currently, non-invasive methods for studying the human brain do not reliably measure spike-rate-dependent signals, independent of other responses such as hemodynamic coupling (fMRI) and subthreshold neuronal synchrony (oscillations and event-related potentials). In contrast, invasive methods – animal microelectrode recordings and human electrocorticography (ECoG) – have recently measured broadband power elevation in field potentials (~50-200Hz) as a proxy for the locally averaged spike rates. Here, we sought to detect and quantify stimulus-related broadband responses using magnetoencephalography (MEG) in individual subjects. Because extracranial measurements like MEG have multiple global noise sources and a relatively low signal-to-noise ratio, we developed an automated denoising technique, adapted from (Kay et al., 2013), that helps reveal the broadband signal of interest. Subjects viewed 12-Hz contrast-reversing patterns in the left, right, or bilateral visual field. Sensor time series were separated into an evoked component (12-Hz amplitude) and a broadband component (60–150Hz, excluding stimulus harmonics). In all subjects, denoised broadband responses were reliably measured in sensors over occipital cortex. The spatial pattern of the broadband measure depended on the stimulus, with greater broadband power in sensors contralateral to the stimulus. Because we obtain reliable broadband estimates with relatively short experiments (~20 minutes), with a sufficient signal-to-noise-ratio to distinguish responses to different stimuli, we conclude that MEG broadband signals, denoised with our method, offer a practical, non-invasive means for characterizing spike-rate-dependent neural activity for a wide range of scientific questions about human brain function. Significance Statement Neuronal activity causes perturbations in nearby electrical fields. These perturbations can be measured non-invasively in the living human brain using EEG and MEG. These techniques have emphasized two kinds of measurements: oscillations and event-related responses. A third type of signal, a stimulus-related increase in power spanning a wide range of frequencies (‘broadband’), is routinely measured in invasive recordings, but not with MEG and EEG. This broadband response is of great interest because unlike oscillations and event-related responses, it is correlated with neuronal spike rates. Here we report quantitative, spatially specific measurements of broadband fields in individual human subjects using MEG. These results demonstrate that a spike-rate-dependent measure of brain activity can be obtained non-invasively from the living human brain.


67
The time-varying electric and magnetic fields near neural tissue provide an indirect but rich source 68 of information about the activity of neural populations (reviewed by Buzsaki et al., 2012). These 69 signals include rapid, 'evoked' responses that are time-locked to stimulus events (Norcia et al.,70 2015), oscillatory responses (Berger, 1929), and non-oscillatory, broadband signals (Miller et al.,71 2007; Miller et al., 2009c). Broadband signals associated with sensory or motor tasks have been 72 widely observed in human electrocorticography, or 'ECoG', (Miller et al., 2014) and animal 73 microelectrode recordings (Henrie and Shapley, 2005). The broadband signal is an elevation in 74 spectral power, typically spanning 50 to >200Hz (Miller et al., 2009b), and has attracted a great 75 deal of attention for several reasons. 76 First, the broadband signal is correlated with the level of neural activity (multi-unit spiking), and 77 hence provides a way to study population-level spiking activity in a cortical region (Liu and 78 Newsome A second challenge in measuring broadband extracranially is that the response is most evident in 94 high frequencies (> 60Hz), and the signal amplitude at these frequencies is low. While intracranial 95 recordings have relatively high signal-to-noise ratios (SNR) even at these higher frequencies (Miller 96 et al., 2014), EEG and MEG do not (Hämäläinen et al., 1993). Broadband signals can extend to lower 97 frequencies (Harvey et al., 2013;Winawer et al., 2013), but oscillatory processes in lower frequency 98 bands often mask broadband measures in this range (Miller et al., 2009c). Here, we sought to measure broadband signals quantitatively in the human brain using a non-107 invasive method (MEG). In order for this important, spike-dependent signal to be useful, it is 108 necessary to measure it reliably in individual subjects, with a high SNR. A high SNR is essential if 109 this signal will be widely used to study differences across stimuli, tasks, or groups. We developed a 110 novel, automated MEG denoising algorithm adapted from prior fMRI work . Our 111 experiments were designed to elicit spatially localized neural responses in visual cortex, and eye 112 movements were measured in a subset of subjects to test for possible confounds from non-neural 113 sources. 114 Subjects provided written informed consent. The experimental protocol was in compliance with the 122 safety guidelines for MEG research and was approved by the University Committee on Activities 123

Methods
involving Human Subjects at New York University and by the ethics committee of the National 124 Institute of Information and Communications Technology (NICT). 125

Display 126
Stimuli were generated using MATLAB (MathWorks, MA) and PsychToolbox (Brainard, 1997 The stimuli were contrast-reversing dartboard patterns (12 square wave contrast reversals per  137 second), windowed within either a half circle (left or right visual field) or full circle (bilateral visual 138 field) aperture, with a diameter of 22 degrees at NYU (26 degrees at CiNet). Mean luminance gray 139 (206 cd/m² (NYU), 83 cd/m² (CiNet)) was used as background color for the dartboards and was 140 shown in the full field during blank trials between stimulus periods ( Figure 1). 141

Experimental design 142
One run consisted of six seconds flickering 'on' periods, alternated with six seconds 'off' mean 143 luminance periods, repeated 6 times (72 seconds). The order of the left-, right-or both-visual field 144 apertures was random. There was a fixation dot in the middle of the screen throughout the run, 145 switching between red and green at random intervals (averaging 3 seconds). The subjects were 146 instructed to maintain fixation throughout the run and press a button every time the fixation dot 147 changed color. The subjects were asked to minimize their blinking and head movements. After 148 every 72-second run, there was a short break (typically 30-s to 1 minute). Each subject participated 149 in 15 runs.

156
Within a run, the order of both-, left-, and right-field flickering periods was randomized. Fifteen runs were obtained per 157 subject, so that there were 30 repetitions of each stimulus type across the 15 runs. The fixation dot is increased in size for 158 visibility. Actual fixation dot was 0.17 degrees in radius (6 pixels).

MEG signal acquisition 160
Data for the main experiment were acquired continuously with a whole head Yokogawa MEG 161 system (Kanazawa Institute of Technology, Japan) containing 157 axial gradiometer sensors to 162 measure brain activity and 3 orthogonally-oriented reference magnetometers located in the dewar 163 but away from the brain area, used to measure environmental noise. The magnetic fields were 164 sampled at 1000Hz and were filtered during acquisition between 1Hz (high pass) and 200Hz (low 165 pass). 166 In a subset of subjects (S6-S8), eye movements were recorded by an EyeLink 1000 (SR Research  167 Ltd., Osgoode, ON, Canada). Right eye position data were continuously recorded at a rate of 1000Hz. 168 Calibration and validation of the eye position was conducted by having the subject saccade to 169 locations on a 5-point grid. Triggers sent from the presentation computer were recorded by the 170 EyeLink acquisition computer. The same triggers were recorded simultaneously by the MEG data 171 acquisition computer, allowing for synchronization between the eye-tracking recording and MEG 172 recording. 173 The 4 data sets acquired with an Elekta Neuromag at CiNet and were pre-processed in MATLAB 174 (MathWorks, MA, USA) using the identical code and procedure. The CiNet data were acquired as 175 102 pairs of planar gradiometer signals (204 sensors). Data were analyzed from each of the 204 176 gradiometers separately and paired into 102 locations for mesh visualization (e.g., the broadband 177 signal-to-noise-ratio for sensor 121 and 122 out of 204 would be averaged to show one signal-to-178 noise-ratio in the position of sensor 61 out of 102 MATLAB. Using either the environmentally-denoised data or raw data, the signals were divided into 200 short epochs. Each stimulus type (left-, right-, or both-hemifield, or blank) was presented in 6-s 201 blocks, and these blocks were divided into 6 non-overlapping 1-s epochs. We discarded the first 202 epoch of each 6-s block to avoid the transient response associated with the change in stimulus. 203 After epoching the data, we used a simple algorithm to detect outliers. We first defined a 'data 204 block' as the 1-s time series from one epoch for one sensor. So a typical experiment consisted of 205 ~170,000 data blocks (157 sensors x 1080 1-s epochs). We computed the standard deviation of the 206 time series within each data block, and labeled a block as 'bad' if its standard deviation was more 207 than 20 times smaller or 20 times larger than the median standard deviation across all data blocks. 208 The time series for bad data blocks were replaced by the time series spatially interpolated across 209 nearby sensor (weighting sensors inversely with the distance). Further, if more than 20% of data 210 blocks were labeled bad for any sensor, then we removed the entire sensor from analysis, and if 211 more than 20% of data blocks were bad for any epoch, then we removed the entire epoch from 212 analysis. Typically, two to seven sensors and 2%-4% of the epochs were removed per session for 213 the NYU data. For the CiNet datasets, almost no sensors or epochs were removed (one sensor and 214 one epoch across all data sets). These preprocessing steps were implemented with 215 dfdPreprocessData.m. 216

Computation of stimulus-locked and broadband responses 217
Data were summarized as two values per sensor and per epoch: a stimulus-locked and a broadband 218 power value. These calculations were done by first computing the Fourier transform of the time 219 series within each epoch (Figure 2A,B). 220 The stimulus-locked signal was then defined as the amplitude at the stimulus-locked frequency 221 (12Hz). The broadband response was computed as the geometric mean of the power across 222 frequencies within the range of 60-150Hz, excluding multiples of the stimulus-locked frequency 223 (see also Figure 2 AB). The geometric mean is the exponential of the average of the log of the signal. 224 We averaged in the log domain because log power is better approximated by a normal distribution 225 than is power, which is highly skewed. These two calculations converted the MEG measurements 226 into a broadband and a stimulus-locked summary metric, each sampled once per second ( Figure  227 2C). The two summary metrics were computed by the functions getstimlocked.m and 228 getbroadband.m. 229 We then bootstrapped across epochs to compute confidence intervals on the signal estimates (per 230 sensor and per condition). For each of 1000 bootstraps, we sampled n epochs with replacement, 231 where n is the total number of epochs in the experiment. We then computed the average response 232 across epochs for each stimulus condition, minus the average across blank epochs. This provided 233 one summary measure for each of the three stimulus conditions and each of the two dependent 234 measures (broadband and stimulus-locked) for each of the 1000 bootstraps. Finally, we took the 235 median across bootstraps as the estimate of signal and half of the 68% confidence interval across 236 bootstraps as the estimate of the noise ( Figure 2D,E). For some analyses, the ratio of these values 237 was defined as the signal-to-noise ratio (SNR). 238  We used the stimulus-locked signal to identify the noise pool because this signal had a very high 289 SNR, and could easily by measured prior to running our denoise algorithm, and because we 290 assumed (and confirmed by inspection) that sensors with broadband responses also had stimulus-291 locked responses. 292 For most subjects, most of the sensors in the noise pool were located over the front of the head (see 293 for example Figure 3A). 294

Filtering of time series 295
As described above, the broadband summary metric was derived from power at a limited range of 296 temporal frequencies (60-150Hz, excluding multiples of the stimulus frequency). After defining the 297 noise pool, the time series of all sensors in all epochs were filtered to remove signal at all 298 frequencies not used to compute the broadband signal. Hence the remaining time series contained 299 power only at frequencies defining the signal of interest. This step was important because the noise 300 pool, though selected for a low stimulus-locked SNR, could nonetheless have contained a small, 301 residual stimulus-locked signal. This residual signal would have been correlated with the 302 experimental design (larger when stimuli were present than absent) and hence projecting it out of 303 the data could have caused a systematic bias (see the script denoisingProjectingInVariance.m). 304

PCA 305
Following filtering, the next step in the algorithm was principal component analysis (PCA). This 306 identified the common components of the time series across the sensors in the noise pool. PCA was 307 computed separately for each 1-s epoch ( Figure 3C). This means that denoising occurred at the 308 same temporal scale (1 second) as the computation of the summary metrics. This differs from some 309 denoising algorithms, in which noise regressors are identified over a much longer time period, e.g., 310 several minutes (Vigario, 1997). Denoising at a short-time scale can be advantageous if the spatial 311 pattern of the noise responses is not consistent across the entire experiment. As a control 312 comparison, we also ran our algorithm by identifying PC time series on the entire duration of the 313 experiment (~20 minutes) rather than epoch by epoch. (See Results, 'Control analyses for MEG 314 Denoise algorithm'.) 315 Projecting out PCA components 316 The first one to ten principal components (PCs) in each epoch were projected out of the time series 317 for all sensors, using linear regression. This resulted in ten new data sets: One with PC 1 projected 318 out, one with PC 1 and 2 projected out, etc. up to 10 PCs projected out ( Figure 3D). After projecting 319 out the noise components, we summarized the data into a stimulus-locked and broadband 320 component as described in Figure 2.

Statistical comparisons 329
To assess the effect of the MEG Denoise algorithm on the broadband SNR, we compared the 330 broadband SNR after applying MEG Denoise to the broadband SNR either without denoising or after 331 applying other denoising algorithms. To make these comparisons, we first identified 10 sensors of 332 interest from each subject. These sensors of interest were the 10 with the highest SNR in any of the 333 three stimulus conditions, either before or after denoising, excluding sensors from the noise pool. 334 For each of the three stimulus conditions, we then took the average SNR from these 10 sensors 335 without denoising or after applying MEG Denoise or another denoising algorithm. Finally, we 336 conducted two-tailed t-tests, paired by subject (n=8), between the broadband SNR after MEG 337 Denoise to the broadband SNR without denoising (or with another algorithm). The t-tests were 338 conducted separately for each of the three stimulus conditions (both-hemifield, left-hemifield, and 339 right-hemifield). 340

Control analyses 341
To investigate the validity of our algorithm, we ran multiple control analyses. In particular, it is 342 important to rule out the possibility that the denoising algorithm produces significant results even 343 when the data contains no sensible signal. To test this, we compared the difference in SNR of 344 denoised data with the following controls: (1) phase-scrambling the PC time series, and (2) using all 345 sensors to define the noise with PCA rather than only a subset of sensors that have little to no 346 stimulus-locked signal. We also assessed the effect of identifying and projecting out PC time series 347 equal in length to the entire experiment (~20 minutes), rather than PC time series matched in 348 length to our analysis epochs (1-s). This comparison tested the assumption that denoising in 349 shorter epochs was advantageous, possibly due to the pattern of noise sources differing over the 350 course of the experiment. three NYU subjects (S6-S8) whether there was a difference in rate between the 'off' baseline 355 periods and 'on' stimulus periods, and within the three stimulus (both-, right-, left-hemifield) 356 conditions. Microsaccades were identified as changes in position with above a relative velocity 357 threshold (6º/s) and a minimum duration of 6 ms, as reported in Engbert & Mergenthaler (2006) to 358 analyze rate and direction of microsaccades as well as separating MEG data into epochs that did and 359 did not contain microsaccades. 360

361
A large field 'on-off' stimulation experiment was used to characterize the stimulus-locked (steady 362 state evoked field, 'SSVEF') and broadband responses in visual cortex measured with MEG. The two 363 measures are reported below, both prior to and after applying our new denoising algorithm. 364

365
In each stimulus condition (left-, right-, and both-hemifield), the stimulus contrast reversed 12 366 times per second, so the stimulus-locked signal was measured at 12Hz and harmonics. Because the 367 largest component was at 12Hz, we defined the stimulus-locked signal for a particular stimulus 368 condition as the amplitude at 12Hz, averaged over all 1-second epochs with that stimulus (typically

384
Both the stimulus-locked and broadband signals were largest in medial, posterior sensors, as 385 expected from activations in visual cortex (Seki et al., 1996). For the stimulus-locked signal, the 386 both-hemifield condition tended to produce broadband signals in bilateral posterior sensors, 387 whereas the single-hemifield conditions produced responses that were lateralized, with higher SNR 388 contralateral to the stimulus. This pattern could be seen in an example subject and in the average 389 across subjects ( Figure 5). The lateralization of the stimulus-locked signal was less clear in the 390 average across subjects due to imperfect alignment of the sensors showing the largest differential 391 response to the left-and right-hemifield stimuli. In each of the 8 individual subjects and in each of 392 the 3 conditions, the stimulus-locked response was evident, with the signal at least 10x above the 393 noise (data not shown). 394 395 Figure 5. Topographic map of stimulus-locked and broadband responses. Data from subject S1 (left) and averaged 396 across subjects S1-S8 by sensor (right). The top 3 rows show data from the 3 stimulus conditions (both-, left-, and right-397 hemifield) compared to blank, and the lower row shows data as the left-only minus right-only conditions. The dependent 398 variable plotted for the single subject data is the signal-to-noise ratio at each sensor, computed as the mean of the 399 contrast (stimulus minus blank) across bootstraps divided by the standard deviation across bootstraps (bootstrapped 400 over epochs). For the group data, the signal-to-noise ratio is the mean of the subject-specific SNRs at each given sensor.

401
The same scale bar is used for all stimulus-locked plots. For the broadband plots, one scale bar is used for the first three 402 rows, and a different scale bar with a smaller range is used for the fourth row. Made with dfdMakeFigure5.m.

403
The spatial pattern of broadband signals was qualitatively similar to the spatial pattern of the 404 stimulus-locked signal, with bilateral posterior responses in the both-hemifield condition, and 405 lateralized responses in the single-hemifield conditions ( Figure 5, individual example and group-406 averaged data). However, the broadband responses had much lower signal-to-noise than the 407 stimulus-locked responses, and in many of the individual subjects, broadband was not evident in 408 one or more conditions (data not shown). The broadband responses were less reliable for the left-409 and right-hemifield conditions than for the both-hemifield conditions. 410 The fact that broadband responses were evident in a few subjects in some conditions indicates that 411 it is possible to measure broadband fields with MEG. However, if this signal cannot be measured 412 reliably in many subjects and many conditions, then the practical value of measuring broadband 413 with MEG is limited. This motivated us to ask whether denoising the MEG data could unmask 414 broadband signals, making it more reliable across subjects and stimulus conditions. 415

417
The MEG data were denoised using a new algorithm as described in detail in the Methods section. 418 In brief, for each subject a subset of sensors that contained little to no stimulus-locked responses 419 were defined as the noise pool. Once the noise-pool was defined, the time series in each sensor and 420 in each epoch was filtered to remove all signals not contributing to the broadband measurement.

435
We first illustrate the effect of denoising with an example from a single sensor in one subject 436 ( Figure 6). This sensor showed a broadband response both prior to, and after, denoising. The 437 benefit of denoising was not evident when comparing the mean power spectra before and after 438 denoising ( Figure 6A). Denoising did not reduce the variability in power across frequencies, nor did 439 it increase the separation in the spectra for the contrast stimulus and the blank. Instead, the effects 440 of denoising are better appreciated by examining the variability across epochs rather than across 441 frequencies ( Figure 6B). The biggest effect is that the broadband power estimates became less 442 variable across epochs, both for the blank condition and the stimulus condition. This is indicated by 443 the narrower distributions in the response amplitudes for the two conditions ( Figure 6B, main 444 panels) and for the difference between conditions ( Figure 6B, insets). The standard deviation of the 445 difference distributions decreased more than two-fold (from 0.79 to 0.35) as a result of denoising. 446 There are two other secondary patterns evident in these distributions. First, the mean broadband 447 power of both the blank and stimulus condition decreased as a result of denoising (for the both-448 hemifield condition, 35.8 versus 26.1, prior to versus after denoising; for the blank, 28.7 versus 449 21.4). This was expected because projecting out signal reduces power. Second, the contrast 450 between the two conditions (difference between the means) reduced: 7.0 prior to denoising versus 451 4.8 after denoising. The combination of these two effects was that the percent difference was little 452 changed, with broadband power from the contrast-stimulus about 25% more than for the blank 453 before and after denoising. Hence denoising did not increase the estimate of the percent signal 454 change. 455 It is important to consider how these effects interact. Because the reduction in variability across 456 epochs was the biggest effect of denoising (more than 2-fold), there was more than a doubling of 457 SNR, computed as the median divided by the variability of the difference distribution. In sum, the 458 spectral plots show that the variability in power across frequencies was little affected by denoising 459 (Figure 6A), whereas the distribution plots show that the variability in total broadband power 460 across epochs was reduced considerably ( Figure 6B). 461 We now consider the effect of denoising across sensors, subjects, and stimulus conditions. 462 Projecting out noise PCs substantially increased the signal-to-noise ratio of the broadband 463 measurement in visually responsive sensors. For example, in the both-hemifield condition for 464 subject S1, the median SNR of the 10 most visually responsive sensors increased from 5 to 10 after 465 denoising ( Figure 7A, blue solid line), similar to the example sensor shown earlier ( Figure 6B). In 466 contrast, the SNR of the 75 sensors in the noise pool was relatively unaffected by denoising ( Figure  467 7A, blue dashed line). This was expected because sensors in the noise pool were unlikely to 468 distinguish stimulus from blank. Across the 8 subjects in the both-hemifield condition, taking the 469 mean of the 10 most visually responsive sensors for each subject, the SNR increased about 3-fold 470 (from 1.6 to 5.0), with a numerical increase in every subject ( Figure 7B). Because the SNR stabilized 471 in all subjects with 10 or fewer PCs projected out, in subsequent analyses, for simplicity we report 472 the effects of denoising with exactly 10 PCs. A comparison of the SNR before denoising (0 PCs 473 projected out) and after (10 PCs projected out) summarized across all subjects and the three 474 stimulus conditions shows increases in SNR for every subject in all conditions ( Figure 7C) 475 (p=0.0001, p=0.0007, p=0.0022 for two-tailed t-tests, 0 v 10 PCs, for both-, left-, and right-hemifield 476 conditions, respectively). 477

486
In principle, the SNR increases could have arisen from increased signal, decreased noise, or both. To 487 distinguish among these possibilities, we compared the signal level alone and the noise level alone 488 before and after denoising. As in prior results, the signal was defined as the difference in broadband 489 power between the contrast pattern and the blank (median across bootstraps), and the noise was 490 defined as the variability of this difference metric (half of the 68% confidence interval across 491 bootstraps). For all three stimulus conditions in most subjects, the signal was largely unaffected by 492 denoising, staying at a similar level or decreasing slightly, while the noise level went down 493 substantially (Figure 8). These analyses indicate that the increase in SNR from denoising ( Figure 7 Expressed as a percentage increase over baseline, the broadband response to the both-hemifield 497 stimulus after denoising was ~10.9±1.7% averaged across the top 10 sensors in each subject (mean 498 ± sem across subject), and 12.6%±1.6% for the top 5 sensors. This contrasts with the much larger 499 stimulus-locked response, which was a nearly 8-fold increase over baseline even prior to denoising 500 (678%±226% increase over baseline for the top 5 sensors).

504
The effect of denoising the broadband signal was not uniform across the sensor array. In general, 505 sensors where we expected visual activity (over the posterior, central part of the head) showed 506 increased SNR following denoising. In the example subject S1 as well as the average across subjects, 507 the denoised broadband response was observed in bilateral sensors for the both-hemifield 508 condition, and with a contralateral bias (relative to the midline) in the two lateralized conditions 509 (Figure 9). For the both-hemifield stimulus, broadband responses were evident in sensors over the 510 posterior, middle of the head in most individual subjects ( Figure 10).

521
To validate the assumptions in our denoising algorithm, we ran three control analyses. In one 522 control analysis, we concatenated all epochs to derive noise regressors from the whole 523 experimental time series (Figure 11, 2 nd bar, compared to using the default of 1-s epochs to derive 524 noise regressors -1 st bar). The elevation in broadband SNR was significantly less when we 525 concatenated all epochs (p = 0.0016, p = 0.0023 and p = 0.0447, for the three stimulus conditions 526 respectively). In the second control analysis, the noise pool included all sensors rather than only 527 those sensors that were not visually responsive. Here, the noise regressors included some signal as 528 well as noise, and hence should be of less benefit. This expectation was confirmed, in that there was 529 no increase in SNR when the algorithm was run with the omission of the noise-pool-selection step 530 (Figure 11, 3 rd bar, p = 0.0014, p = 0.0015 and p = 0.0020 for the three stimulus conditions 531 respectively). In a 3 rd control analysis, we phase-scrambled each of the epoch-by-epoch noise time 532 series. The phase-scrambled regressors were temporally uncorrelated with the actual time series in 533 the noise. As a result, we found no change in SNR levels ( Figure 11,  MEG sensors which face away from the head and measure environmental rather than physiological 550 fields. By design, these algorithms project out time series from the subspace spanned by the 551 reference sensors, thereby reducing environmental noise, but not physiological noise. Applying 552 either one of these two algorithms alone to the 8 data sets reported above increased the broadband 553 signal-to-noise ratio, evident in the group-averaged sensor plots ( Figure 12A, columns 3-4 versus 554 column 2), and the increased SNR in the 10 most responsive sensors ( Figure 12B, 2 nd and 3 rd bar 555 versus 1 st bar in each plot). 556 In planned comparisons, we evaluated the SNR increase of each algorithm or combination of 557 algorithms to the increase from MEG Denoise alone. The increase from each of the two 558 environmental algorithms alone was significantly less than that from our new MEG Denoise 559 algorithm ( Figure 12A, column 5 versus columns 3-4; Figure 12B, 4 th bar versus 2 nd and 3 rd ). 560 Applying two algorithms in sequence, first either CALM or TSPCA, followed by MEG Denoise, also 561 resulted in a large increase in broadband SNR ( Figure 12A, columns 6 and 7). For all three stimulus 562 conditions, the combination of MEG Denoise and CALM resulted in the largest gain in SNR, 563 significantly larger than MEG Denoise alone for two out of the three conditions ( Figure 12B, 5 th 564 versus 4 th bars). This indicates that the MEG Denoise algorithm and an environmental algorithm 565 captured some independent noise. 566

577
In a separate analysis, we ran the MEG Denoise algorithm to evaluate its effect on the stimulus-578 locked signal. The methods were identical to those used to denoise the broadband signal except for 579 the omission of one step, the step in which we filtered the time series to remove temporal 580 components that do not contribute to the broadband signal. Denoising modestly increased the 581 stimulus-locked SNR for all stimulus conditions for most subjects (Figure 13, top). The SNR 582 increased numerically in all subjects (n=8) and in all stimulus conditions, although the percentage 583 increases were smaller than those for denoising the broadband signal, ~20% increase compared to 584 two-fold. As in the case of denoising the broadband signals, the main contribution to the increase in 585 SNR for the stimulus-locked signal was a decrease in variability across epochs (Figure 13, bottom), 586 rather than an increase in the signal level (Figure 13, middle). 587

593
To test whether the findings reported above generalize to other instruments and experimental 594 environments, we conducted the same experiment using a different type of MEG system, an Elekta 595 360 Neuromag at CiNet. The CiNet system contains paired planar gradiometers, in contrast to the 596 axial gradiometers used in the Yokogawa MEG at NYU, and the scanner is situated in a different 597 physical environment, with potentially very different sources of environmental noise. The pre-598 processing pipeline at this imaging center often includes a denoising step based on temporally 599 extended signal source separation (tSSS) (Taulu and Simola, 2006;Taulu and Hari, 2009 The identical experiments were conducted with 4 new subjects. As expected, all three stimulus 605 types led to a large stimulus-locked response in the posterior sensors, with a peak SNR of more 606 than 10 in the group averaged data ( Figure 14A, column 1). A modest, spatially specific broadband 607 signal was measured from the undenoised data for each stimulus type ( Figure 14A column 2), with 608 a peak SNR of 1-2 in the group-average data for all three conditions. Unlike the NYU data, in the 609 CiNet data the MEG Denoise algorithm on the raw data did not generally result in an increase in the 610 broadband SNR (group data, Figure 14A, columns 2 and 3; individual subjects, Figure 14B, left side 611 of each subplot). However, when the raw data were pre-processed with the tSSS algorithm ( Figure  612 14A, column 4), application of MEG Denoise increased the SNR in all 3 stimulus conditions for 3 out 613 of 4 subjects, and in 2 out of 3 stimulus conditions for the 4th subject. Together, the MEG Denoise 614 algorithm increased the SNR by 2-3 fold, similar to the NYU data (both-hemifield: 2.8 to 5.6; left-615 hemifield: 0.8 to 2.4; right-hemifield: 2.01 to 4.4; means across subjects 1-4, top 10 sensors each, for 616 the tSSS data and the MEG Denoised tSSS data). Just as with the NYU MEG data set, the combination 617 of an algorithm tailored to find environmental noise (tSSS) and our algorithm produced the most 618 robust results, indicating that MEG Denoise and the environmental denoising algorithm removed at 619 least some independent sources of noise. 620

662
We separated the MEG signal into two components, one time-locked and one asynchronous with 663 the stimulus. The stimulus-locked component was clearly visible in all subjects with minimal 664 preprocessing. The asynchronous signal, spanning 60-150Hz, was visible with little preprocessing 665 in some subjects and the mean across subjects. However, the SNR was low compared to the 666 stimulus-locked component. With our denoising algorithm, SNR more than doubled, resulting in 667 reliable, spatially specific broadband signals in all individuals. We showed in a subset of subjects 668 that the broadband signals could not be explained by systematic biases in the pattern of fixational 669 eye movements, supporting the interpretation that the broadband fields arise from neural activity 670 rather than artifacts associated with eye movements. 671 These results are qualitatively consistent with intracranial measurements . 672 However, it has proven difficult to measure extracranial broadband signals arising from neural 673 activity. Below, we discuss the significance of broadband responses, challenges in measuring them 674 extracranially, and the generalizability of our denoising algorithm. 675

676
A century ago Berger and others described oscillations in surface EEG between 10 and 25Hz 677 (Berger, 1929;Adrian and Matthews, 1934). More recently, using intracranial recordings, Crone 678 and colleagues (1998) reported power increases in higher frequencies (75-100Hz) associated with 679 motor movements. Subsequently, this high frequency response was interpreted as a broadband 680 (not oscillatory) signal, thought to reflect increased neuronal activity, rather than increased 681 synchrony (Miller et al., 2007;Miller et al., 2009c). In support of this interpretation, studies showed 682 that broadband signals correlate with single (Manning et al., 2009) and multiunit spike rates (Ray 683 and Maunsell, 2011), and are observed throughout cortex ( Figure 16). Under some conditions, 684 broadband is also correlated with the BOLD signal (  The broadband response occurs in many brain areas for many types of stimuli, spanning at least 50-704 150Hz, and can extend to lower and higher frequencies (Miller et  First, MEG sensors pool over a large area, so baseline power reflects activity from a large fraction of 734 the brain, whereas visually driven broadband responses likely come from confined regions 735 (Krusienski et al., 2011). In contrast, both baseline and visual responses in ECoG electrodes arise 736 from the same cortical patch. Second, the amplitude depends not only on pooling area, but also 737 phase coherence. If broadband signals arise from incoherent neural activity, and stimulus-locked 738 signals from coherent (synchronous) activity, then the former will grow with the square-root of the 739 number of sources, and the latter with the number of sources. Since MEG pools over much larger 740 populations than ECoG, the ratio of incoherent signal strength (broadband) to coherent (stimulus-741 locked) will be much lower. This logic is supported by modeling (Linden et al., 2011) and empirical 742 studies with intra-and extracranial measures, which found that the most coherenet intracranial 743 signals were best transmitted outside the head (Pfurtscheller and Cooper, 1975 and EEG sensors. Many of these noise sources are spectrally broad and hence particularly 749 problematic when investigating neural broadband signals. 750 Although spike fields generated from eye movements can be mistaken for broadband neural 751 activity (Yuval-Greenberg and Deouell, 2009), it is unlikely that our spatially-specific broadband 752 measures were substantially contaminated by eye movement artifacts. This was confirmed by 753 analyses of eye movement data, and the fact that middle posterior sensors where we observed 754 broadband are not usually associated with MEG spike field artifacts . A second eye 755 movement confound, the electromagnetic fields arising from movement of the retina-to-cornea 756 dipole, causes low frequency artifacts (4-20Hz; (Keren et al., 2010)) and therefore is unlikely to 757 have affected our broadband measures (60-150Hz). 758 Head muscles can also cause spectrally broadband contaminants (Muthukumaraswamy, 2013), as 759 can external noise sources, e.g., nearby electrical equipment. However, these noise sources are 760 unlikely to be confined to occipital sensors and to co-vary with stimulus condition, and hence do 761 not explain our broadband observations. Moreover, it is likely that these noise sources, if present, 762 were included in our noise pool, and hence MEG Denoise would have reduced their effects. 763 MEG Denoise and other denoising algorithms 764 MEG Denoise uses PCA on a subset of sensors to remove noise. In principle, it can capture any noise 765 source contributing to the noise pool, including environmental, oculomotor, muscular, and neural. 766 This differs from algorithms designed to remove environmental noise. Hence MEG Denoise is 767 complementary to these methods. We found that the most effective analysis was either MEG 768 Denoise alone, or MEG Denoise following an environmental denoising algorithm. Our algorithm has 769 much in common with ICA denoising (Vigario, 1997), with some important differences. First, PCA, 770 unlike ICA, ranks components by variance explained. Second, MEG Denoise explicitly identifies 771 noise sensors. These features enable the algorithm to be fully automated, making it easy to denoise 772 at the time scale of individual events (e.g., >1,000 one-second epochs). If the spatial pattern of the 773 PCs vary over time, then deriving the components independently within short epochs is more 774 effective ( Figure 11). 775 To use MEG Denoise for other experimental designs, analyses, or scanners, one would need to 776 change some input parameters. In addition to the experimental design matrix and data, required 777 inputs include the experiment-specific functions to summarize the MEG responses. In our 778 experiments, one function computed the stimulus-locked signal and was used to define the noise 779 pool. For most of our analyses, a second function computed the broadband power to evaluate the 780 results. In principle, one could use a single function to define the noise pool and evaluate the data 781 (as we did for denoising the stimulus-locked signal). For other experiments, one might use a 782 function that computes the amplitude or latency of an evoked response, or the power in a limited 783 temporal frequency band, or any measure relevant to the experiment. Alternatively, one could run a 784 separate localizer experiment to identify a pool of potential sensors of interest and a pool of noise 785 sensors, and then manually enter the list of noise sensors to denoise the main experiment. There 786 are several other optional inputs, such as the method to identify the noise pool, the accuracy metric 787 (SNR/R 2 ). Here, we used the defaults for all optional inputs. 788

789
Stimulus-driven broadband brain responses can be quantitatively characterized in individual 790 subjects using a non-invasive method. Because we obtain high SNR measures from short 791 experiments, the broadband signal can be used to address a wide range of scientific questions. 792 Access to this signal opens a window for neuroscientists to study signals associated with neuronal 793 spike rates non-invasively at a high temporal resolution in the living human brain. 794