Laminar dynamics of high amplitude beta bursts in human motor cortex

Motor cortical activity in the beta frequency range is one of the strongest and most studied movement-related neural signals. At the single trial level, beta band activity is often characterized by transient, high amplitude, bursting events rather than slowly modulating oscillations. The timing of these bursting events is tightly linked to behavior, suggesting a more dynamic functional role for beta activity than previously believed. However, the neural mechanisms underlying beta bursts in sensorimotor circuits are poorly understood. To address this, we here leverage and extend recent developments in high precision MEG for temporally resolved laminar analysis of burst activity, combined with a neocortical circuit model that simulates the biophysical generators of the electrical currents which drive beta bursts. This approach pinpoints the generation of beta bursts in human motor cortex to distinct excitatory synaptic inputs to deep and superficial cortical layers, which drive current flow in opposite directions. These laminar dynamics of beta bursts in motor cortex align with prior invasive animal recordings within the somatosensory cortex, and suggest a conserved mechanism for somatosensory and motor cortical beta bursts. More generally, we demonstrate the ability for uncovering the laminar dynamics of event-related neural signals in human non-invasive recordings. This provides important constraints to theories about the functional role of burst activity for movement control in health and disease, and crucial links between macro-scale phenomena measured in humans and micro-circuit activity recorded from animal models.

Theories of cortical computation have proposed distinct roles for neural activity in various frequency channels and projections originating from deep and superficial cortical layers ( Adams et al., 2013 ;Arnal and Giraud, 2012 ;Bastos et al., 2018Bastos et al., , 2015Bastos et al., , 2012Donner and Siegel, 2011 ;Fries, 2015Fries, , 2005Friston and Kiebel, 2009 ;Jensen et al., 2015 ;Jensen and Mazaheri, 2010 ;Stephan et al., 2017 ;Wang, 2010 ), but relatively few of these theories address the generation of frequencyspecific activity by inter-laminar dynamics within cortical circuits. Recent computational neural models of beta burst generation ( Law et al., 2019 ;Neymotin et al., 2020 ;Sherman et al., 2016 ) suggest that such bursts can be generated by a short, strong excitatory synaptic input to superficial cortical layers, temporally aligned with a broad, weaker input to deep layers lasting. These synchronized inputs cause current to propagate in opposite directions within a cortical column, resulting in a cumulative dipole with the stereotypical wavelet shape in the time domain as measured by local field potentials (LFPs), electroencephalography (EEG), and magnetoencephalography (MEG; Karvat et al., 2020 ;Kosciessa et al., 2020 ;Little et al., 2019 ;Sherman et al., 2016 ). The prominent peak in the beta burst wavelet shape is generated by the strong supragranular synaptic input that drives current towards deep layers of the cortex. The surrounding tails of the beta burst peak emerge from synaptic input to the deep layers that drive current flow toward superficial layers. Laminar current source density recordings from primary somatosensory cortex in rodents and non-human primates provided initial support for these model-derived predictions ( Sherman et al., 2016 ).
Until recently, testing such circuit-level predictions noninvasively in the human brain was constrained by the low temporal resolution of fMRI and the low spatial precision of EEG and MEG. However, recent developments in high precision MEG (hpMEG; Little et al., 2018 ;Meyer et al., 2017 ;Troebinger et al., 2014b ) have demonstrated sensitivity of the recorded signals to the orientation of cortical columns ( Bonaiuto et al., 2020 ), and the ability to test hypotheses concerning the laminar dominance of frequency-specific induced activity ( Bonaiuto et al., 2018a , b ). We here extended these techniques to develop a temporally resolved laminar analysis which yields estimates of the relative strength of superficial and deep cortical layer activity over the time course of preand post-movement beta bursts. This method thus allows for empirical testing of predictions derived from neural circuit data, non-invasively in humans.
To validate the ability to recover a priori known patterns of laminar activity, we first used a detailed biophysical model of beta burst generation to create simulated MEG data using a generative MEG model. We then tested the ability of the new laminar analysis to recover the simulated strong superficial layer activity at the peak of the burst and deep layer activity surrounding the peak. We then generated predictions for sensor level activity and laminar dynamics from a range of alternative synthetic models with varying dipole laminar locations and current flow directions. Finally, we tested the ability of the analysis to correctly infer the source of laminar activity for a range of simulated co-registration er-ror magnitudes and signal-to-noise ratios. These simulations verify that the new laminar analysis works when the ground truth is known, allow alternative models to be ruled out based on comparison with results from human MEG data, and establish that the human MEG data is of high enough quality to run the analysis.
When applying these analyses to human MEG data, we found that the laminar profile of both pre-and post-movement high amplitude beta bursts in motor cortex conformed to the biophysical model's predictions: activity surrounding the burst peak localized to deep cortical layers, whereas the peak corresponded to activity predominantly in superficial layers. When compared against alternate models of burst generation, we found support only for the model in which bursts are generated by distinct synchronous inputs to deep and superficial cortical layers, which in turn drive current flow in opposite directions. These results thus provide novel demonstration of the laminar dynamics of beta burst generation in human motor cortex, and provide unique linkage between macro-scale phenomena of burst activity measured in humans and their micro-circuit mechanism observed in animal models. More generally, we demonstrate the possibility for uncovering the laminar-resolved dynamics of human burst activity in non-invasive MEG recordings.

Human subject data
Data from eight healthy, right-handed, volunteers with normal or corrected-to-normal vision and no history of neurological or psychiatric disorders were used for our analyses (six male, aged 28.5 ± 8.52 years; Bonaiuto et al., 2018a ;Little et al., 2019 ). The study protocol was in accordance with the Declaration of Helsinki, and all participants gave written informed consent which was approved by the UCL Research Ethics Committee (reference number 5833/001).

MRI acquisition
Two MRI scans were acquired prior to MEG scanning with a 3T whole body MR system (Magnetom TIM Trio, Siemens Healthcare, Erlangen, Germany) using the body coil for radio-frequency (RF) transmission and a standard 32-channel RF head coil for reception. The first was a standard T1-weighted scan for creation of each subject's individual head-cast ( Meyer et al., 2017 ), and the second was a high resolution, quantitative multiple parameter map (MPM; Weiskopf et al., 2013 ) for MEG source location.
The T1-weighted protocol used a 3D spoiled fast low angle shot (FLASH) sequence with 1 mm isotropic image resolution, field-of view of 256, 256, and 192 mm along the phase (anterior-posterior, A-P), read (head-foot, H-F), and partition (right-left, R-L) directions, respectively, repetition time of 7.96 ms, and excitation flip angle set to 12°. A single echo was acquired after each excitation to yield a single anatomical image. A high readout bandwidth (425 Hz/pixel) was used to preserve brain morphology, with no significant geometric distortions observed in the images. The acquisition time for this sequence was 3 min 42 s. A 12channel head coil was used for signal reception without either padding or headphones in order to avoid introduction of scalp distortions.
The MPM protocol was comprised of three differentially-weighted, RF and gradient spoiled, multi-echo 3D fast low angle shot (FLASH) sequences and two additional calibration sequences to correct for inhomogeneities in the RF transmit field ( Callaghan et al., 2015 ;Lutti et al., 2012Lutti et al., , 2010, with whole-brain coverage at 800 m isotropic resolution. The FLASH sequences had predominantly proton density (PD), T1, or magnetization transfer saturation (MT) weighting. The PD-and MTweighted volumes used a flip angle of 6°, and the T1 weighted acquisition used a flip angle of 21°. MT-weighting was achieved through the application of a Gaussian RF pulse 2 kHz off resonance with 4 ms duration and a nominal flip angle of 220°prior to each excitation. The field of view was 256 mm head-foot, 224 mm anterior-posterior (AP), and 179 mm right-left (RL). Gradient echoes were acquired with alternating readout gradient polarity at eight equidistant echo times ranging from 2.34 to 18.44 ms in steps of 2.30 ms using a readout bandwidth of 488 Hz/pixel, but only six echoes were acquired for the MT-weighted volume in order to maintain a repetition time (TR) of 25 ms for all FLASH volumes. Partially parallel imaging using the GRAPPA algorithm was employed to accelerate data acquisition, with a speed-up factor of 2 in each phase-encoded direction (AP and RL) with forty integrated reference lines.
To maximize the accuracy of the measurements, inhomogeneity in the transmit field was mapped by acquiring spin echoes and stimulated echoes across a range of nominal flip angles following the approach described in Lutti et al. (2010) , including correcting for geometric distortions of the EPI data due to B0 field inhomogeneity. The total acquisition time for all MRI scans was less than 30 min.
Quantitative maps of proton density (PD), longitudinal relaxation rate (R1 = 1/T1), MT and effective transverse relaxation rate (R2 * = 1/T2 * ) were subsequently calculated ( Weiskopf et al., 2013 ). FreeSurfer (v5.3.0;Fischl, 2012 ) was used to reconstruct pial and white matter surfaces from the MPM volumes for source localization of MEG sensor data. We used a custom FreeSurfer procedure to process MPM volumes, using the PD and T1 volumes as inputs ( Carey et al., 2017 ), resulting in surface meshes representing the pial surface (adjacent to the cerebro-spinal fluid, CSF), and the white/gray matter boundary. We used a custom routine to downsample each surface by a factor of 10 while maintaining the correspondence between surface vertices ( Bonaiuto et al., 2020 ). This yielded two meshes of the same size (same number of vertices and edges), containing about 30,000 vertices each ( M = 30,094.8, SD = 2665.5 over participants).

Head -cast construction
From an MRI-extracted 3D scalp model, we constructed a head-cast that fit between the participant's scalp and the MEG dewar ( Meyer et al., 2017 ;Troebinger et al., 2014b ). Scalp surfaces were first extracted from the T1-weighted volumes acquired in the first MRI protocol using SPM12 ( http://www.fil.ion.ucl.ac.uk/spm/ ). This tessellated surface, along with 3D models of fiducial coils placed on the nasion and the left and right pre-auricular points, was then placed inside a virtual version of the scanner dewar in order to minimize the distance to the sensors while ensuring that the participant's vision was not obstructed. The resulting model (including spacing elements and ficudial coil protrusions) was printed using a Zcorp 3D printer (Zprinter 510). The 3D printed model was then placed inside a replica of the MEG dewar and polyurethane foam was poured in between the model and dewar replica to create the participant-specific head-cast. The fiducial coil protrusions in the 3D model therefore became indentations in the foam head-cast, into which the fiducial coils were placed during scanning. The locations of anatomical landmarks used for co-registration are thus consistent over repeated scans, allowing us to merge data from multiple sessions ( Meyer et al., 2017 ;Troebinger et al., 2014b ).

Behavioral task
Participants completed a visually cued, action-based decision making task in which they responded to an instruction cue projected on a screen by pressing one of two buttons using the index and middle finger of their right hand ( Bonaiuto et al., 2018a ). After a baseline fixation period, a random dot kinematogram (RDK) was displayed for 2 s with coherent motion to the left or to the right. Following a 500 ms delay period, an instruction cue (an arrow pointing to the left or the right), prompted participants to respond by pressing either the left or right button. The level of RDK motion coherence and the congruence between the RDK motion direction and instruction cue varied from trial to trial, but for the purposes of the present study, we collapsed across conditions and analyzed beta activity before and after all button press responses. See Bonaiuto et al. (2018a) for a complete description of the task paradigm and structure.

MEG acquisition and preprocessing
MEG data were acquired using a 275-channel Canadian Thin Films (CTF) MEG system with superconducting quantum interference device (SQUID)-based axial gradiometers (VSM MedTech, Vancouver, Canada) in a magnetically shielded room. A projector was used to display visual stimuli on a screen ( ∼50 cm from the participant), and a button box was used for participant responses. The data collected were digitized continuously at a sampling rate of 1200 Hz. MEG data preprocessing and analyses were performed using SPM12 ( http://www.fil.ion.ucl.ac.uk/spm/ ) using MATLAB R2014a. The data were filtered (5th order Butterworth bandpass filter: 2-100 Hz, Notch filter: 50 Hz) and downsampled to 250 Hz. Eye blink artifacts were removed using multiple source eye correction ( Berg and Scherg, 1994 ). Trials were then epoched from 2 s before the participant's response to 2 s after. Blocks within each session were merged, and trials whose variance exceeded 2.5 standard deviations from the mean were excluded from analysis. Preprocessing code is available at http://github.com/jbonaiuto/meg-laminar .

Burst definition and analysis
Burst were defined using a two-stage process. Firstly, the sensor data time series was inverted onto the subject specific pial surface mesh (note the explicit bias here towards the pial surface is simply to obtain a robust time series estimate so that bursts can be identified) using an Empirical Bayesian beamformer algorithm (EBB; Belardinelli et al., 2012 ;López et al., 2014 ) as implemented in SPM12. Source inversion was completed for a 4 second window of data epoched to the button press. These data were projected into 274 orthogonal spatial (lead field) modes and 4 temporal modes and current density time series were created by multiplying the sensor data by the weighting matrix (M) between sensors and source from the inversion, and the data reduction matrix (U) that specifies the significance of the data modes that map to the cortex. The time series from the vertex closest to the center of the hand knob in the primary motor cortex was then taken forward for burst analysis. The selected current density time series was filtered using a 4th order Butterworth filter in the beta band (13 -30 Hz) and the amplitude determined using the Hilbert function. In order to isolate high amplitude bursts with a similar waveform shape, bursts were defined using an empirically defined threshold ( Little et al., 2019 ) of 1.75 standard deviations above the median beta amplitude. Raw (aside from the broad bandpass and notch filters applied in preprocessing) sensor data were then re-epoched around the peak of the beta burst amplitude ( ± 100 ms). Additionally, a more accurate epoching of bursts was achieved through using the raw time series data and a Woody filter ( Woody, 1967 ) to align each burst epoch to the average burst template (equivalent to cross correlating individual bursts to the average burst and shifting the epoch window according to the lag that maximizes the cross-correlation). We used an error tolerance of 0.1, 100 maximum iterations, and unbiased cross-correlation. This new burst aligned sensor dataset was then used for the time-resolved laminar analyses.
Burst amplitude was compared between pre-and post-movement epochs using a linear mixed model in R (v3.6.3; R Core Team, 2020 ) with the lme4 package (v1.1-21; Bates et al., 2014 ). Minimum beta amplitude in the 50 ms surrounding the burst peak was treated as the de-pendent measure, with epoch as a fixed effect (pre-or post-movement), and subject-specific intercepts as a random effect. Fixed effect significance was tested using a Type II Wald chi-square test.

Localizer source reconstruction
Source inversion for beta burst localization was performed using the empirical Bayesian beamformer algorithm (EBB) as implemented in SPM12. The source inversion was applied to a 200 ms time window, centered on the peak of the average beta burst time course, without Hann windowing. These data were projected into 274 orthogonal spatial (lead field) modes and 4 temporal modes. These inversions used a spatial coherence prior ( Friston et al., 2008 ) with a FWHM of 5 mm. We used the Nolte single shell head model ( Nolte, 2003 ) with the source locations constrained by the vertices of the downsampled cortical surface. We used Bayesian model evidence to compare five methods of estimating cortical columns for constraining source orientations: downsampled surface normal vectors ( Dale and Sereno, 1993 ;Fuchs et al., 1994 ;Hämäläinen and Hari, 2002 ;Hillebrand and Barnes, 2003 ;Lin et al., 2006 ), original surface normal vectors ( Bonaiuto et al., 2020 ), cortical patch statistics ( Lin et al., 2006 ), link vectors ( Dale et al., 1999 ), and variational vector fields ( Fischl and Sereno, 2018 ). The link vectors method yielded the best model fit (pre-movement bursts: 6/8 subjects; post-movement bursts: 7/8 subjects; Fig. S1), and therefore source orientations were constrained according to the link vectors between the pial and white matter surfaces to approximate the orientation of cortical columns ( Bonaiuto et al., 2020 ).
Clusters of vertices with activity above a threshold of 80% of the maximum activity, within a mask of 50 mm centered on the hand knob of the left precentral gyrus were carried forward to the sliding time window source reconstruction.
Induced pre-and post-movement beta activity localizes primarily to deep cortical layers ( Bonaiuto et al., 2018a ). This is also predicted by the biophysical model, which generates beta bursts consisting of two tails reflecting activity in deep layers and one brief peak of superficial layer activity, resulting in a net bias towards deep layers. We confirmed this by performing a Bayes factor comparison between EBB source inversions using pial versus white matter generative models over the entire burst time course, which yielded greater model evidence for the white matter generative model. To account for this predicted net laminar bias and increase sensitivity to potential superficial layer activity, we used the pial surface to localize beta bursts for the subsequent sliding time window source reconstruction.

Sliding time window source reconstruction
We used an adaptive Woody filter to align the beta burst source time series across subjects ( Woody, 1967 ). All subsequent analysis was based on these temporally aligned datasets. The sliding time window source reconstruction was performed using the Multiple Sparse Priors (MSP; Friston et al., 2008 ) algorithm as implemented in SPM12. This was applied to a 40 ms time window of the aligned average beta burst time course, with Hann windowing, with a frequency of interest of 1-256 Hz. As with the localizer source reconstruction, we used a spatial coherence prior with a FWHM of 5 mm, and the Nolte single shell model ( Nolte, 2003 ) with link vector orientations ( Bonaiuto et al., 2020 ). The time window was advanced along the time course of the average beta burst in increments of one time step (4 ms), and the Bayesian model evidence (approximated by free energy) for the generative model was computed within each window. Within each cluster of vertices identified by the localizer source reconstruction, the sliding time window source reconstruction was conducted using all vertices with a link vector angle within 0.1 radians of that of the cluster peak vertex ( Bonaiuto et al., 2020 ). For each selected vertex, the sliding time window MSP inversion was conducted using the pial surface to constrain source locations with the prior set as the cluster vertex, and again using the white matter surface with the prior set as the corresponding vertex on the white matter surface. The difference in the free energy time series between these inversions ( ΔF = F pial -F white matter ) was then averaged over vertices within each cluster, and then across clusters. These free energy difference time series of the aligned burst data were then summed over subjects to yield a fixed effects estimate of laminar dominance dynamics.

Biophysical model
We used the open-source Human Neocortical Neurosolver (HNN) software to simulate our biophysical model of a local neocortical microcircuit under exogenous layer specific synaptic drive ( https://hnn.brown.edu ; Neymotin et al., 2020 ). HNN's model, and the beta burst mechanism, was fully described and validated in prior publications ( Law et al., 2019 ;Sherman et al., 2016 ;Shin et al., 2017 ). HNN's underlying neural model simulates the primary electrical currents in the neocortex that create EEG/MEG signals. The model simulations are based on the biophysical origins of the primary electrical currents (i.e., current dipoles), assumed to be generated by the postsynaptic, intracellular longitudinal current flow in the long and spatially aligned dendrites of a large population of synchronously activated neocortical pyramidal neurons. HNN simulates the primary currents from a canonical model of a layered neocortical column via the net intracellular electrical current flow in the pyramidal neuron dendrites, in the direction aligned with the apical dendrites, multiplied by their length (units nano-Ampere-meters). A scaling factor is applied to the net current dipole output to match the amplitude of recorded data and provides an estimate of the number of neurons contributing to the recorded signal. Simulated multicompartment pyramidal (PN) and single compartment interneurons (IN) are arranged in supra-and infra-granular layers. Neurons receive excitatory synaptic input from simulated trains of action potentials in predefined temporal profiles (see schematic Gaussians in Fig. 3 A) that activate excitatory synapses at the location of the proximal apical/basal and distal apical dendrites of the pyramidal neurons. These two input pathways represent the structure of feedforward (proximal) and feedback (distal) inputs to a laminar cortical circuit. The model is agnostic as to the brain areas providing this input, which could be either thalamic or cortical. In prior studies, we have presumed based on the literature that the thalamus is the source of the proximal and distal beta generating drives in primary cortical areas ( Law et al., 2021 ;Sherman et al., 2016 ). Of note, in Sherman et al. (2016) we examined several alternative beta generating mechanisms and found that those applied here ( Fig. 3 ) provided the closest fit to the human data, with confirmation from invasive laminar electrophysiology in mice and monkeys. The results presented here establish that the source and directionality of the intra-laminar current flow during beta bursts, as predicted in Sherman et al. (2016) and shown in Fig. 3 B, are consistent with laminar resolution source estimation in humans. These results do preclude alternative mechanisms of generation of this pattern of current flow.
In the simulations described in this paper, we used a modified version of the default parameter set distributed with HNN that simulates beta bursts, namely AlphaAndBeta.param file. In brief, this simulation contained 100PN and 35IN per layer and used a stochastic sequence of 10 Hz proximal excitatory synaptic drive, simultaneous with distal 10 Hz excitatory synaptic drive. This pattern of drive was motivated from the fact that 10 Hz alpha rhythmicity is one of the dominant operational models of the thalamus during spontaneous states, and accounted for the alpha/beta SI mu-complex in our earlier work ( Jones et al., 2009 ). The histogram of driving spikes on each cycle of the rhythmic inputs that generated this drive had a Gaussian profile. We have shown that this pattern of input can generate bursts of activity with the characteristic beta event waveform shape that emerges as part of the more continuous somatosensory mu-rhythm ( Jones et al., 2009 ;Sherman et al., 2016 ), but that 10 Hz rhythmicity is not necessary for generation of a single beta burst and the pattern shown in Fig. 3 A is sufficient ( Sherman et al., 2016 ). Here, we examine the model output corresponding to individual bursts only. Individual beta bursts or "events " occur on cycles of the stochastic 10 Hz drive when broad ( ∼100 ms) upward current flow from proximal inputs is synchronously disrupted by downward stronger and faster ( ∼50 ms) current flow from distal inputs, where the ∼50 ms duration of the distal drive creates the dominant beta peak and sets the frequency of the oscillation ( Sherman et al., 2016 ). We used three sets of simulations in order to ascertain the contribution of each cortical layer to the beta waveform shape: 1) proximal drive only, 2) distal drive only; and 3) combined proximal and distal inputs. Each simulation type was run across 50 trials, where a trial (170 ms) refers to a single execution of the model with a defined set of simulation parameters. Results varied across trials with identical parameters due to the stochastic nature of the exogenous proximal and distal drives: on each cycle of the 10 Hz drive, the timing of the synaptic drive is chosen from a Gaussian distribution with mean inter-cycle input time of 100 ms ± 20 ms standard deviation for proximal drive and 100 ms ± 7 ms for distal drive. The mean delay between the proximal and distal drive is 0 ms. Note that the distal drive synaptic weights provided to supra-and infra-granular pyramidal neurons were increased from HNN's default value to 6e-5 S. The distal and proximal current dipole moments were normalized by the maximum of the averaged amplitudes of each, then scaled to match the amplitudes of recorded data. All source code is provided online at https://hnn.brown.edu and https://github.com/jonescompneurolab/hnn , and a parameter file specific to the simulations here will be distributed upon publication.

Generation of simulated MEG datasets
All simulations were based on a dataset acquired from one human participant, which was used to determine the sensor layout, sampling rate (1200 Hz, downsampled to 250 Hz), number of samples (51), and simulated dipole locations and orientations for the simulations. In each simulation, we specified spatially distributed source activity centered at a single vertex on the pial surface and the corresponding vertex on the white matter surface, approximating the deep and superficial ends of a cortical column ( Bonaiuto et al., 2020 ). The orientations of the simulated dipoles were defined by the vector connecting the two vertices (link vector) with opposite polarities: the pial dipole pointed toward the white matter dipole and vice versa. The time course of simulated activity at the pial vertex was given by the cumulative dipole moment from 50 beta bursts generated by the biophysical model run with only the distal drive applied, and that of the white matter vertex by the cumulative dipole moment from 50 bursts generated by the model run with only the proximal drive applied. The amplitudes of these time courses were scaled to yield beta bursts matching those seen in the human data, with the pial and white matter vertex dipoles having mean magnitudes of 8 and 6 nAm, respectively. We then spatially smoothed these simulated dipole time courses with a Gaussian kernel (FWHM = 5 mm), to obtain two patches of spatially distributed activity. We then used a single shell forward model ( Nolte, 2003 ) to generate a synthetic dataset from the simulated source activity. Unless otherwise specified, Gaussian white noise was added to the simulated data and scaled in order to yield a per-burst amplitude SNR level (averaged over all sensors) of − 20 dB.
Simulated data from the simplified synthetic model (Fig. S2) were generated in the same way as with the biophysical model, but the pial and white matter vertex dipoles had Gaussian activity time courses, centered within the epoch. The pial vertex Gaussian had a width of 10 ms and a magnitude of 6 nA m, and the white matter vertex Gaussian had a width of 25 ms and a magnitude of 4.5 nAm.

SNR and co -registration error simulations
Beta bursts were generated using the simplified model and Gaussian white noise was added to the simulated data and scaled in order to yield per-trial amplitude SNR levels (averaged over all sensors) between − 50 and − 20 dB to generate synthetic datasets across a range of realistic SNRs (typically ranging from − 40 to − 20 dB; Goldenholz et al., 2009 ). This was repeated 50 times per SNR level and the minimum free energy difference during the tails of the burst ( − 40 -− 10 ms and 10 -40 ms) and maximum free energy difference during the burst peak ( − 10 -10 ms) was computed.
Similarly, we simulated between-session co-registration error by introducing a linear transformation of the fiducial coil locations in random directions (0 mm translation and 0°rotation up to 2 mm translation and 2°rotation) to generate a realistic range of uncertainty concerning the location of the brain relative to the MEG sensors ( Adjamian et al., 2004 ;Gross et al., 2013 ;Ross et al., 2011 ;Singh et al., 1997 ;Stolk et al., 2013 ;Whalen et al., 2008 ). In these simulations, the per-trial amplitude SNR was set to − 20 dB. This was repeated 50 times per co-registration error level, and the minimum free energy difference during the tails of the burst and maximum free energy difference during the burst peak was computed.

Consistent waveform shapes of beta bursts in human motor cortex
We first identified bursts of precentral beta activity in MEG data from human participants performing a visually cued action decisionmaking task involving button press responses made using the right hand ( Bonaiuto et al., 2018a ). Bursts were identified at the source level as periods when beta amplitude in the hand area of left motor cortex exceeded an empirically defined threshold (burst duration: M = 94.08, SD = 29.19 ms; Little et al., 2019 ). We segmented each participant's preprocessed (not beta bandpass filtered) dataset into 200 ms epochs, centered on the peak of each burst, and split each dataset into bursts that occurred prior to ( M = 326.8 SE = 133.3 bursts per subject; 2614 total bursts), or after ( M = 722.5, SE = 155.0 bursts per subject; 5780 total bursts) the button press response in each trial. After alignment, the average time series of pre-movement beta bursts matched the stereotyped wavelet-like shape previously described in somatosensory cortex ( Sherman et al., 2016 ), and appeared as a single dipolar field pattern centered over left sensorimotor cortex with a transient reversal in direction around the burst peak ( Fig. 1 A). This mean pre-movement burst waveform shape was consistent across all participants ( Fig. 1 B). Beta bursts that occurred post-movement had the same spatial and temporal features as pre-movement bursts ( Fig. 1 C), but had a greater peak amplitude ( X 2 (1) = 95.77, p < 0.001; Fig. 1 D). This suggests that a common mechanism may underlie the generation of pre-and post-movement motor beta bursts.

Biophysical computational modeling predicts a specific laminar mechanism for beta burst generation
We next sought to determine the putative layer-resolved cortical mechanisms behind the generation of beta bursts in motor cortex. For this purpose, we extended previously developed hpMEG techniques ( Bonaiuto et al., 2018a( Bonaiuto et al., , 2018b to yield a temporally resolved estimate of laminar activity. This new analysis involved construction of generative MEG models of oriented current patches from white matter and pial surfaces extracted from each subject's MRI volume ( Fig. 2 A), representing deep and superficial cortical layers. Each subject's burst-aligned data were then averaged, and source inversion was conducted on the average beta burst time series to localize the bursts in source space ( Fig. 2 B). After determining the source space locations of precentral beta bursts, we then proceeded to use a sliding time window analysis to compare the relative strength of current flow in deep and superficial layers over the time course of the bursts. At each pial surface vertex selected by the localizer inversion, we identified the corresponding white matter surface vertex in the direction of the estimated orientation of the cortical column at that location using the link vectors method ( Bonaiuto Fig. 1. Sensor-level pre-and post-movement beta bursts. A) Single subject sensor-level broadband data aligned to the peak of premovement motor beta bursts. The mean spatial topography (top) indicates a field reversal over left sensorimotor cortex centered on the burst peak. Individual pre-movement burst waveforms (middle), taken from the preprocessed (not beta bandpass filtered) sensor data indicated by a black circle in the spatial topography (MLP34), exhibit a stereotyped, wavelet-like shape that becomes more apparent when averaged over all pre-movement bursts (shaded area indicates the standard error of the burst waveform across all bursts). B) Mean pre-movement burst waveforms from all subjects ( n = 8; 2614 total bursts) after alignment with a Woody filter (top) and averaged over subjects (bottom). The shaded area indicates the standard error of the mean burst waveform across subjects. Post-movement motor beta bursts have the same spatial and temporal pattern as pre-movement bursts at the single subject (C) and group level (D; 5780 total bursts). et al., 2020 ), as this method resulted in the best model fit compared to other methods tested (Fig. S1). These vertex pairs were then used as Bayesian priors in a sliding time window source inversion ( Fig. 2 C). This involved determining the more likely generative (pial and white matter) model, given the burst activity, within a small time window of the average burst waveform, using Bayesian model evidence (approximated by free energy) comparison. Following each time window comparison, we advanced the time window to obtain a time series of the free energy difference between the models (approximating the Bayes factor; Fig. 2 C). The resulting time series provides an estimate of the likelihood that the electrical currents generating the field signal at each time point were stronger either in deep or superficial layers (as approximated by our pial and white matter surfaces).

Laminar -resolved activity can be identified from simulated MEG data
Before applying this analysis to human MEG data, we first verified that this approach could reliably identify the time course of laminar activity in simulated data for which the ground truth was known. We therefore simulated a set of beta bursts using the biophysical model, as in our previous work ( Fig. 3 A, 19). In this model, the synaptic drive is generated by trains of action potentials in predefined temporal patterns, whose Gaussian shaped histograms are schematically depicted in Fig. 3 A, that activate fast excitatory synapses at layer specific locations as shown. This drive generates intracellular current flow up or down within the cortical pyramidal neuron dendrites, and the cumulative current dipole moment is estimated from the net longitudinal intracellular currents across the pyramidal neurons, multiplied by their length (see SI Appendix, Biophysical model for further details; Neymotin et al., 2020 ;Sherman et al., 2016 ). When only the superficial distal drive is present, the cumulative current dipole exhibits a downward deflecting current into the cortex ( Fig. 3 B, top), and when only the deep proximal drive is present, the cumulative current dipole exhibits an upward deflecting current out of the cortex ( Fig. 3 B, bottom). The cumulative dipole moment from the combination of these inputs results in a dipole moment that resembles the recorded beta burst waveform ( Fig. 3 B, middle).
We then created a generative MEG model based on the anatomy of a human subject and simulated a deep and a superficial source at the location identified by the localizer source inversion for that subject. Both sources were oriented according to the estimated cortical column orien-

Fig. 2.
Sliding time window laminar source inversion. A) Pial and white matter surfaces are extracted from quantitative maps of proton density and T1 times from a multi-parameter mapping MRI protocol. B) Source inversion over the entire burst time course was used to localize the average beta burst sensor signal (inlay). The peak pial surface vertex and the corresponding vertex on the white matter surface were used as priors in the following sliding window inversion. C) Sliding time window source inversion was performed using a 40 ms wide window. For each iteration, source inversion was run using a pial generative model with the pial vertex from the localization inversion as a prior, and using a white matter generative model with the corresponding white matter vertex as a prior. The difference in free energy between the two models (F pial -F white matter ) was used to determine the laminar locus of dominant activity as the window advanced along the average time series of the beta burst (example data shown from a single human subject). tation at that location ( Bonaiuto et al., 2020 ). Following the biophysical model prediction, the superficial current pointed in the direction of the white matter surface (i.e. into the cortex), whereas the deep current flow pointed in exactly the opposite direction (out of the cortex). The time course of each current in the generative MEG model was determined by the cumulative dipole moment of the biophysical model with only the distal (superficial dipole) or proximal (deep layer dipole) synaptic drive applied. This generative model was then used to create simulated MEG sensor data which had the same spatial and temporal features observed in the human MEG data, providing support for the model predictions ( Fig. 3 D). This process was repeated 10 times to yield a dataset of virtual subjects. Source inversion on these simulated datasets revealed the same wavelet-like burst waveform shape in source space current density ( Fig. 3 E). As expected, the sliding time window source inversion identified the bilaminar pattern of simulated activity used to generate the data ( Fig. 3 F): with significantly more evidence for the white matter (deep) model at the beginning and end of the burst ( − 50 -− 18 ms and 6 -50 ms) and for the pial (superficial) model around the peak of the burst ( − 14 -2 ms). These results generalized beyond the particular time series generated by the biophysical model, as we obtained the same results using a highly simplified model with deep and superficial currents given by scaled Gaussian time series (Fig. S2), as well as using the biophysical model with pink rather than white noise (Fig. S3).
Having validated that the new temporally resolved laminar analysis could correctly reconstruct the time course of the relative strength of deep and superficial layer activity in simulated data, we next wanted to determine if the analysis could rule out a range of alternative synthetic generative models which differed in terms of their deep and superficial dipole moment waveforms and polarities. We first tested a synthetic model in which the distal drive consisted of a temporally broad and weak signal, directed into the cortex, whilst the proximal drive consisted of a stronger and briefer signal, directed out of the cortex ( Fig. 4 A). This resulted in a simulated sensor dataset with an oppositely oriented field and reversed polarity waveform to that observed in the human subject MEG data ( Fig. 4 B). Importantly, the temporally resolved laminar analysis was able to correctly identify that this dataset was generated by a pattern of activity that was stronger in superficial layers at the beginning and end of the burst, and deep layers around the peak of the burst ( Fig. 4 C). We next tested a model in which the deep and superficial layer dipoles pointed away rather than towards each other ( Fig. 4 D). The simulated sensor data from this model also had oppositely oriented spatial and temporal features compared to those observed in human subject MEG data ( Fig. 4 E), and the temporally resolved laminar analysis again correctly identified the time course of laminar activity ( Fig. 4 F). We tested a range of other synthetic models including two models containing single deep or superficial dipoles in isolation, and models where dipole currents were oriented in the same direction. In each of these alternative models, the analysis correctly identified the time course of dominant laminar activity (Fig. S4). The only model whose simulated dataset matched that of the human subject data in terms of field direction and waveform shape was the original biophysical model in which strong superficial layer activity was temporally aligned with broad, weaker activity in deep layers. The superficial layer synaptic activity drove net current flow toward the deep layers and vice versa for the deep layer activity. It should be noted that these alternate synthetic models are not intended to realistically explain the generation of beta bursts in sensorimotor cortex, and the examples in Fig. 4 were generated without using the biophysical model. For example, it is likely biophysically impossible for deep layer synaptic inputs to drive current flow away from the superficial layers. However, these simulations demonstrate that the temporally resolved laminar analysis can, in principle, distinguish between different patterns of laminar activity in human MEG data.

Temporally resolved estimates of lamina -specific currents require high precision MEG data
The spatial precision of MEG data is classically limited by its signalto-noise ratio (SNR) and error in localization of the brain relative to the Fig. 3. A biophysical model of beta burst generation predicts a bilaminar time course of beta burst dynamics. A) Beta bursts were generated by a model with a broad proximal excitatory synaptic drive temporally aligned with a strong distal synaptic drive. B) The model was run with just the distal drive (top) and just the proximal drive (bottom), and the resulting cumulative dipole moments were used as source signals to generate the simulated MEG sensor data. Each line shows the time series from a single burst simulation ( n = 50) and the average is shown as a black line. The middle panel shows the cumulative dipole moment generated from the model with both drives, exhibiting the same waveform and time-frequency (inlay) features observed in the human MEG data. C) The generative model included two oppositely oriented currents positioned at corresponding locations on the pial and white matter surfaces in motor cortex, which represent the superficial and deep cortical layers, with source activity given by the model run with the distal and proximal synaptic drives, respectively. D) Simulated sensor data generated by the model has the same spatial and temporal features as beta bursts observed in the human subject MEG data used to determine the simulated current locations and orientations (inlay). E) Time course of source current density resulting from the localizer source inversion on the simulated sensor datasets. Each simulation ( n = 10) is shown as a colored line and the black line corresponds to the average over all simulations. F) The sliding window source inversion correctly identifies that the simulated bursts were generated by activity predominately in deep layers at the beginning and end of the burst, with stronger superficial layer activity at the peak of the burst. The dashed lines (at ΔF = ± 3) show the point at which one laminar model is 20 times more likely than the other model. MEG sensors (co-registration error; Hillebrand and Barnes, 2002;Troebinger et al., 2014b ). We used subject-specific 3D printed head-casts to reduce within-session movement and reduce between-session head position variability, allowing us to reduce co-registration error and increase SNR by averaging across sessions ( Meyer et al., 2017 ;Troebinger et al., 2014b ). Before applying the new analysis to human MEG data, we first wanted to determine if the SNR and co-registration accuracy achieved by this approach were within the range required for accurate estimates of lamina-specific currents. We therefore used the simplified beta burst model (Fig. S2) to generate simulated MEG datasets with a range of SNR and co-registration error levels, and examined the resulting laminar biases during the tails and peak of the simulated bursts (see Methods). At low SNR levels, no significant laminar biases could be detected, but as SNR increased above − 30 dB, the deep bias at the tail ends of the burst and superficial bias at the peak of the burst could be resolved ( Fig. 5 A). Moreover, we saw that the ability to resolve this laminar bias slowly decreased with increasing co-registration error, and deep and superficial biases could no longer be reliably detected when this error increased beyond 2 mm and 2° ( Fig. 5 B). However, the SNR and co-registration error of our human subject data were well within this bounds (premovement burst SNR: M = 299.72, SE = 2.23 dB; post-movement burst SNR: M = 295.15, SE = 1.82 dB; between-session movement variability < 0.6 mm; within-session movement maximum = 0.22 mm; Meyer et al., 2017 ).

The laminar activity pattern of motor beta bursts in human MEG data conforms to predictions from biophysical modeling
After validating the temporally resolved laminar analysis on simulated data, and estimating the required SNR and co-registration accuracy, we then proceeded to apply this analysis to pre-and postmovement beta bursts from human MEG data. High amplitude premovement beta bursts from motor cortex at the source level had the same wavelet-like current density waveform observed in the sensor data ( Fig. 6 A). As predicted by the biophysical model, activity at the beginning and end of pre-movement bursts predominately occurred in deep cortical layers, whereas activity at the peak of the bursts was biased superficially ( Fig. 6 B). Prior to the burst peak, the free energy difference reached a minimum of − 34.73, indicating that the deep layer model was more strongly supported by the MEG data than the superficial layer model (a free energy difference of ± 3 indicates that one model is approximately 20 times more likely than the other). After the burst peak, the Fig. 4. Alternate synthetic models with reversed proximal and distal dipole moment waveforms or polarities generate different predictions. A) Synthetic model with reversed distal (top) and proximal (bottom) dipole moment waveforms from which simulated MEG data were generated. B) The synthetic model with reversed dipole moment waveforms generated simulated sensor data with oppositely oriented spatial and temporal features, compared to beta bursts observed in the human MEG data used to determine the simulated current locations and orientations (inlay). C) The sliding window source inversion correctly identifies that the bursts generated by the reversed dipole moment waveforms model were generated by activity predominately in superficial layers at the beginning and end of the burst, with stronger deep layer activity at the peak of the burst. D) Alternative synthetic model with reversed dipole polarities. E) Similar to the reversed dipole moment waveforms model, the reversed dipole polarity model generated simulated sensor data with oppositely oriented spatial and temporal features, compared to bursts in human MEG data (inlay). F) The sliding time window source inversion can determine the laminar source of dominant activity in the reversed dipole polarity model.

Fig. 5.
High precision MEG data enables localization of laminar dynamics. A) At SNR levels below − 30 dB (dashed vertical line), it is impossible to detect laminar biases during the tail ends and peak of simulated bursts (shaded area shows standard error). Above − 30 dB the full time course of laminar activity can be resolved. B) When co-registration error is less than 2 mm and 2°(dashed vertical line), it is possible to resolve the deep and superficial biases during the tails and peak of simulated beta bursts (shaded area shows standard error). Above 2 mm and 2°, laminar bias detection is unreliable. minimum free energy difference was − 31.55. By contrast, the free energy difference reached a maximum of 49.15 around the peak of burst, suggesting that the superficial layer model was more strongly supported by the data than the deep layer model. High amplitude post-movement beta bursts had the same wavelet-like current density waveform shape at the source level ( Fig. 6 C), and had the same bilaminar pattern activity as pre-movement bursts and as predicted by the biophysical model (preburst minimum ΔF = − 71.32; post-burst minimum ΔF = − 48.02; burst peak ΔF = 51.47; Fig. 6 D). These results were robust to choice of anal-ysis parameter values (Figs. S5-6), individual subject differences (Fig.  S7), and cortical column orientation (Fig. S8).
We ran a control analysis to ensure that these results were spatially specific, using the same burst-averaged data as the main analysis, but with random vertices on the pial surface (and corresponding vertices on the white matter surface) as Bayesian priors in the sliding time window source inversion. This control analysis was designed to ensure that the bilaminar pattern of activity characterizing precentral beta bursts was specific to the cortical location in which these bursts were occurring. Fig. 6. Pre-and post-movement motor beta bursts exhibit the predicted bilaminar dynamics. A) The aligned pre-movement beta burst source level current density time courses averaged across subjects (shaded area shows standard error). B) Activity at the beginning and end of pre-movement bursts localized to deep cortical layers, whereas activity at the peak of the bursts localized superficially, in agreement with the results from the biophysical model (inset, reproduced from Fig. 3 F). This was not true for spatially randomized (blue) surrogate data which yielded a flat Bayes factor time course that was not biased to either surface (shared area shows standard error). Post-movement beta bursts had similar source level temporal (C) and laminar (D) characteristics as pre-movement bursts (inset demonstrates how ΔF is computed).
The control analysis yielded flat estimates of laminar activity, without any significant bias toward the deep or superficial layers, for both preand post-movement bursts ( Fig. 6 C,D). This was also true for a temporal control analysis, which involved shuffling the bursts in the time domain before averaging them, prior to running the localizer and sliding time window source inversions (Fig. S9).

Discussion
We combined temporally resolved laminar inverse analysis and biophysical modeling to show that movement-related high amplitude beta bursts in human motor cortex arise from distinct and temporally aligned deep and superficial layer excitatory synaptic drives that generate alternating dominant dipole currents in deep and superficial layers. The deep layer excitatory synaptic input initially drives current toward the superficial layers (i.e., out of the cortex), and this is transiently interrupted by the superficial layer input which drives current in the opposite direction to generate the dominant peak in the beta burst waveform. This pattern of activity is in accordance with predictions from a biophysical model of somatosensory beta burst generation ( Jones et al., 2009 ;Law et al., 2019 ;Neymotin et al., 2020 ;Sherman et al., 2016 ), suggesting that motor and somatosensory beta bursts share a common mechanism.
In addition to identifying the mechanism of high amplitude beta burst generation in human motor cortex, we validate the use of high precision MEG for temporally resolved estimates of the relative strength of deep and superficial cortical activity. Previously developed techniques were only able to localize the laminar source of induced and temporally averaged activity over a relatively wide time window ( Bonaiuto et al., 2018a , b ;Troebinger et al., 2014a ), and were only able to determine relative laminar biases due to differences in SNR in the beta and gamma band signals analyzed ( Bonaiuto et al., 2018a ). Here we have extended these methods to resolve the time course of laminar neural dynamics, analyzed putative deep and superficial layer signals with comparable SNR, and taken advantage of the ability of hpMEG to infer the orientation of cortical columns ( Bonaiuto et al., 2020 ) in order to anatomically constrain current orientations and polarity in laminar MEG generative models. The resulting method yields estimates of the time series of laminar dominance reflective of synaptically induced currents that are in-line with depth electrode recordings from rodents and nonhuman primates ( Sherman et al., 2016 ), providing further support for this family of laminar hpMEG analyses.
We focused this analysis on high amplitude, short duration, beta bursts, which have a stereotyped waveform shape ( Little et al., 2019 ;Sherman et al., 2016 ), in order to maximize SNR for our sliding time window analysis of the average beta burst time course. However, previous studies have demonstrated bursts of beta activity can also have longer duration and possibly lower amplitude ( Feingold et al., 2015 ;Jasper and Penfield, 1949 ;Kilavik et al., 2012 ;Murthy and Fetz, 1992 ;Saleh et al., 2010 ;Salenius et al., 1997 ;Sanes and Donoghue, 1993 ;Seedat et al., 2020 ). The amplitude threshold used here was chosen to maximize the correlation between the number of bursts detected and the time-averaged beta amplitude ( Little et al., 2019 ). Moreover, lagged coherence, a measure of signal rhythmicity ( Fransen et al., 2015 ), demonstrates that the phase of beta band activity is unpredictable over three cycles ( Little et al., 2019 ). These high amplitude bursts therefore account for the majority of beta power in our data.
One interesting possibility is that lower amplitude and longer busts of beta activity have a distinct function and generation mechanism from short, high amplitude bursts. In support of this idea, traveling waves of beta events in macaque motor cortex occur in distinct spatiotemporal patterns depending on their amplitude ( Rubino et al., 2006 ), with longer burst events sometimes exhibiting multiple patterns in which velocity is tightly linked to amplitude ( Denker et al., 2018 ). Using a hidden Markov model rather than an amplitude threshold for burst detection, the highest amplitude and shortest duration (250-300 ms) beta bursts are found in sensorimotor cortex ( Seedat et al., 2020 ), however when the number of states in the model was increased, the duration of bursts detected was reduced. This may be because the extra states allowed the model to account for more diversity in detected beta bursts. We used a threshold on the amplitude envelope of beta bandpass filtered data to detect bursts, which was therefore not sensitive to differences in burst peak frequency and therefore may have averaged over laminar biases at the beginning and end of the bursts. However more refined methods could be used to detect and cluster bursts varying in peak frequency, amplitude, and duration ( Echeverria-Altuna et al., 2021 ;Zich et al., 2020 ). One such promising technique is empirical mode decomposition ( Huang et al., 1998 ), which could also be used to extract and align burst waveform dynamics using phase-aligned instantaneous frequency ( Fabus et al., 2021 ;Quinn et al., 2021 ), rather than the Woody filter used here.
Motor and sensory cortex differ in terms of their afferent projections from other regions ( Arikuni et al., 1988 ;Felleman and Van Essen, 1991 ;Friedman et al., 1986 ;Jones et al., 1978 ;Tokuno and Tanji, 1993 ), but have the same layer-specific input patterns ( Kuramoto et al., 2009 ). In line with this, our results suggest a common mechanism for the generation of beta bursts across sensorimotor cortex (though this does not imply that all cortico-basal beta bursts are cortically generated). However, beta bursts in these regions have distinct relationships with behavior ( Little et al., 2019 ;Shin et al., 2017 ), for example, pre-and post-movement bursts in motor cortex are related to different aspects of motor responses ( Little et al., 2019 ), and may therefore have distinct functions. Identifying the source of the deep and superficial layer projections to motor cortex that drive pre-and post-movement beta bursts will help to shed light on their functional role in motor behavior. Future extensions to the technique developed here could identify these sources by incorporating measures of time-varying functional connectivity ( Astolfi et al., 2008 ;Hesse et al., 2003 ;Youssofzadeh et al., 2016 ), allowing for inference of lamina-specific connectivity.
We restricted our analyses to a predefined ROI within contralateral motor cortex. We found support for a model of beta burst generation initially developed to explain activity in somatosensory cortex, but transient bursts of beta activity have also been found in prefrontal regions were they have been implicated in working memory and actionstopping Jana et al., 2020 ;Lundqvist et al., 2016 ;Wessel, 2020 ). Our results suggest that the beta frequency peak in time frequency decomposition of sensorimotor activity is largely a consequence of the waveform shape of high amplitude "beta " bursts rather than the amplitude of oscillatory activity. Prefrontal beta bursts might therefore be characterized by a different waveform shape that could be generated by distinct mechanisms. A comparison of burst generation mechanisms across cortical regions exhibiting bursts of beta activity is an interesting avenue for future research.
There are some limitations to the temporally resolved laminar analysis presented here. Like previously developed laminar hpMEG techniques ( Bonaiuto et al., 2018a , b ), the new analysis can only determine the relative bias of deep versus superficial layer activity; it cannot recover the time course of activity in each layer. However, coupling the analysis with a biophysical model and testing a range of alternative synthetic models allows inference of lamina-specific activity time courses by determining the directionality of deep and superficial currents. Future extensions involving generative models with surfaces representing intermediate cortical layers, and higher SNR data afforded by cryogen-free MEG sensors (optically-pumped magnetometers; OPMs; Boto et al., 2016Boto et al., , 2018Holmes et al., 2018;Iivanainen et al., 2017Iivanainen et al., , 2019Knappe et al., 2014 ), may be able to use an approach similar to current source density methods ( Pettersen et al., 2008 ;Rappelsberger et al., 1981 ;Schroeder et al., 1998Schroeder et al., , 1991Schroeder et al., , 1990 to increase the spatial specificity of laminar inference. An interesting question is whether the same study could have been conducted without head-casts but with appropriate head-tracking. We do not know the answer, but our simulations showing co-registration tolerances of 2 mm and 2° ( Fig. 5 B) suggest that this would have been challenging.
Although we found deep and superficial biases at different periods over the time course of beta bursts, these results are in fact consistent with the finding that temporally averaged pre-and post-movement induced beta activity predominates in deep cortical layers ( Bonaiuto et al., 2018a ). The biophysical model adopted here generates beta bursts consisting of two tails of deep layer activity and a transient peak of superficial layer activity, resulting in an overall bias towards deep layers. We compensated for this predicted net deep layer bias by using a pial surface localizer, but future work could overcome potential depth biases by adopting an approach using localization with an intermediate layer surface or temporal basis functions generated by biophysical model predictions. More generally, these results question the notion that activity in various frequency channels is categorically segregated into different cortical layers. Instead, they suggest a view in which the generation of frequency-specific activity has a temporal and laminar profile, and only appears to be restricted to deep or superficial layers when measuring the net bias through the lens of temporally averaging induced activity. More generally, our work bridges between macro-scale phenomena measured in humans and micro-circuit activity recorded from animal models, and therefore facilitates probing the mechanisms behind beta bursts. The combination of biophysically principled neural circuit modeling together with precise characterization of cortical activity in humans provides novel ways to connect data from human and animal models.
In earlier work ( Troebinger et al., 2014a ), we showed that underestimation of the true patch extent will tend to bias estimates towards deeper layers, whereas an overestimation of patch extent will tend to bias layer estimates more superficially. In our parameter sensitivity analyses (Fig. S6), however, we found laminar inference to be robust over a range of realistic patch sizes. Here we had a specific assumption that activity is driven from a combination of deep and superficial sources with equal patch extents. However, we should mention that an alternative hypothesis, which we did not explore, is that the bursts arise from a single layer within an expanding and contracting local coherence in current flow.
Our analysis aligned MEG sensor data to beta amplitude peaks and treated these signals as event-related fields (ERFs). Sensory and motor ERFs (and event-related potentials; ERPs) have waveforms which are composed of a series of dynamic components generated by temporal patterns of activity in both deep and superficial cortical layers ( Mehta et al., 2000a( Mehta et al., , 2000bOlson, 2001 ;Schall et al., 2020 ;Schroeder et al., 1992Schroeder et al., , 1991Szymanski et al., 2011 ), but studies combining laminar electrode and ERP recordings from nonhuman primates are rare ( Woodman, 2012 ). The sliding time window inversion technique presented here can also be used to non-invasively test hypotheses, potentially constrained by biophysical modeling ( Jones et al., 2009( Jones et al., , 2007Neymotin et al., 2020 ), concerning the generation of event-related neural activity by interlaminar dynamics.
The finding that temporally transient deep and superficial layer excitatory synaptic drives generate movement-related beta bursts in human motor cortex provides relevant mechanistic constraints to theories of the functional relevance of pre-and post-movement beta activity ( Alegre et al., 2008( Alegre et al., , 2002Cao and Hu, 2016 ;Cassim et al., 2001 ;Classen et al., 1998 ;Engel and Fries, 2010 ;Fetz, 2013 ;Jasper and Penfield, 1949 ;Murthy and Fetz, 1992 ;Pfurtscheller et al., 1996 ;Reyns et al., 2008 ;Roelfsema et al., 1997 ;Saleh et al., 2010 ;Spitzer and Haegens, 2017 ;Tan et al., 2016Tan et al., , 2014. Such theories are often framed in terms of slowly modulated beta oscillations, and generally predict sustained beta activity arising from either sustained inputs or recurrently maintained activity. On the contrary, the transitory confluence of deep and superficial layer synaptic drives in beta burst generation suggests a dynamic functional role for motor beta bursts ( Little et al., 2019 ), such as the integration of different sensorimotor information signals during pre-and post-movement periods.  J Bonaiuto). SL was funded by a clinical research training grant from the Wellcome Trust ( 105804/Z/14/Z ). SAN was funded by an NIDCD (R01DC012947-06A1) and an Army Research Office Grant (W911NF-19-1-0402). SRJ and SAN were supported by NIH R01EB022889 and NIH R01MH106174. The Wellcome centre for Human Neuroimaging is supported by a strategic award from Wellcome (091593/Z/10/Z). The funders had no role in the preparation of the manuscript.

Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.neuroimage.2021.118479 .