Mapping language from MEG beta power modulations during auditory and visual naming

Evaluation of language dominance is an essential step prior to epilepsy surgery. There is no consensus on an optimal methodology for determining language dominance using magnetoencephalography (MEG). Oscillatory dynamics are increasingly recognized as being of fundamental importance for brain function and dysfunction. Using task-related beta power modulations in MEG, we developed an analysis framework for localizing and lateralizing areas relevant to language processing in patients with focal epilepsy. We examined MEG responses from 29 patients (age 42±13 years, 15M/14F) during auditory description naming (ADN) and visual picture naming (PN). MEG data were preprocessed using a combination of spatiotemporal filtering, signal thresholding, and ICA decomposition. Beta-band 17-25Hz power decrements were examined at both sensor and source levels. Volumetric grids of anatomical source space were constructed in MNI space at 8mm isotropic resolution, and beta-band power changes were estimated using the dynamic imaging of coherent sources beamformer technique. A 600ms temporal-window that ends 100ms before speech onset was selected for analysis, to focus on later stages of word production such as phonologic selection and motor speech preparation. Cluster-based permutation testing was employed for patient- and group-level statistical inferences. Automated anatomic labeling atlas-driven laterality indices (LIs) were computed for 13 left and right language- and motor speech-related cortical regions. Group localization of ADN and PN consistently revealed significant task-related decrements of beta-power within language-related areas in the frontal, temporal and parietal lobes as well as motor-related regions of precentral/premotor and postcentral/somatomotor gyri. A region-of-interest analysis of ADN and PN suggested a strong correlation of r=0.74 (p<0.05, FDR corrected) between the two tasks within the language-related brain regions, with the highest spatial overlap in the prefrontal areas. Laterality indices (LIs) consistently showed left dominance (LI > 0.1) for most individuals (93% and 82% during ADN and PN, respectively), with average LIs of 0.40±0.25 and 0.34±0.20 for ADN and PN, respectively. Source analysis of task-related beta power decrements appears to be a reliable method for lateralizing and localizing brain activations associated with language processing in patients with epilepsy.


Introduction
Language mapping is a critical step before epilepsy and brain tumor surgery. Careful delineation of cortical regions associated with language processing can minimize the risk of developing post-surgical functional deficits in this patient group. While language task-related magnetoencephalographic (MEG) responses have been studied by many groups in patients undergoing presurgical evaluations for epilepsy or brain tumors (Doss et al., 2009;Hirata et al., 2010;Papanicolaou et al., 2004;Raghavan et al., 2017), there is no consensus on an optimal methodology for delineating language areas within the cortex or determining hemispheric language dominance. This may in part be due to the widely distributed and multi-faceted nature of cortical networks associated with language function (Binder et al., 1997;Fedorenko and Thompson-Schill, 2014;Friederici, 2011;Price, 2000), limitations of MEG responses or source modeling methods in localizing all of the regions that are engaged by a task, and the neurophysiological effects of epilepsy on the language network (Raghavan et al., 2017).
MEG is a measure of neural activity with excellent temporal resolution (<1 ms), moderate spatial precision (mm), and the ability to measure cortical rhythms across a wide range of frequencies (0.1 to >100 Hz). The neural processing of words is composed of various subprocesses including acoustic, phonetic, phonologic, semantic, and syntactic analysis, and the temporal properties of these processes can potentially be resolved by time-sensitive methods like MEG (Helenius et al., 2002;Indefrey and Levelt, 2004). Functional magnetic resonance imaging (fMRI) has been widely used to localize language areas in the cortex, however, this technology lacks the temporal resolution to elucidate the transient dynamics of language processes. Although MEG has been used to assess hemispheric language dominance as an alternative to the invasive Wada test (Wada and Rasmussen, 1960), precise localization and lateralization of language areas in individual patients remain a challenge using current MEG methodologies (Kuchcinski et al., 2015;Papanicolaou et al., 2004;Rodin et al., 2013).
Classic neurophysiological studies have used time-locked electromagnetic changes at stimulus onset, collectively referred to as evoked responses to study language comprehension processes (Kutas and Van Petten, 1988). In particular, the N400 component of the evoked response has been extensively studied as a measure of brain activity related to semantic processing (Baggio and Hagoort, 2011;Brown and Hagoort, 1993;Koelsch et al., 2004;Lau et al., 2008). However, attempts to localize later components (e.g., P600 generators) have failed to provide reliable source models (Friederici et al., 2003;Friederici, 2011). While evoked responses mainly represent neural activity that is time-and phase-locked to the stimulus, responses that have significant jitter in their trial-to-trial latency are better detected by averaging power changes in different frequency bands across trials (Pfurtscheller and Lopes da Silva, 1999;Tallon-Baudry and Bertrand, 1999). There is now a compelling body of empirical evidence from electrocorticography (ECoG), local field-potential (LFP) recordings, MEG, and EEG that local high-frequency activity in the gamma band (>60 Hz) accompanied by a broader field of power-decreases in the alpha (8-12 Hz) and beta (13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30) bands are a signature of cortical engagement during tasks (Crone et al., 2006(Crone et al., , 1998Eulitz et al., 1996;Miller et al., 2007;Pfurtscheller, 1991;Pfurtscheller and Andrew, 1999;Singh et al., 2002). Broad-band high gamma activity is related to the firing rates of local neuronal populations (Manning et al., 2009;Nir et al., 2007;Ray and Maunsell, 2011), is more focally expressed than concomitant decreases in alpha and beta power (Crone et al., 1998;Edwards et al., 2005;Miller et al., 2007), and also appears to represent the EEG/MEG correlate of the fMRI BOLD signal (Hall et al., 2014;Logothetis et al., 2001;Mukamel et al., 2005;Singh, 2012;Singh et al., 2002). Although imaging of high-gamma activity would be ideal to localize cortical areas engaged by a language task, detectability of high-gamma activity outside primary sensory cortical regions in EEG and MEG remains challenging, and only a few recent studies have reported success (Dansereau et al., 2014;Hashimoto et al., 2017;Lee et al., 2010). This likely derives from factors such as signal loss due to the mixing of out-of-phase rhythms at the scalp-electrode/MEG-sensor (Pfurtscheller and Cooper, 1975) and the presence of a substantial myogenic noise floor at these frequencies (Jerbi et al., 2009;Whitham et al., 2007). Consequently, many MEG studies have instead exploited power decrements at lower frequencies, especially in the beta band, to localize cortical engagement during language tasks (Eulitz et al., 1996;Fisher et al., 2008;Grabner et al., 2007;Kadis et al., 2008;Liljestr€ om et al., 2005;Pang et al., 2011;Passaro et al., 2011).
A few MEG studies have used beta power decrements for lateralizing and localizing brain areas relevant to language production and/or comprehension in patients being evaluated for neurosurgical procedures (Findlay et al., 2012;Fisher et al., 2008;Hirata et al., 2010Hirata et al., , 2004. These findings support the value of beta power decrements as a marker of cortical engagement during language processing. Despite the strong association of beta-band power modulation with task-related activity, further evidence is needed to support its use as a clinical tool for language mapping. In the present work, we developed a beta power decrement mapping framework to provide localization and lateralization of language areas relevant to two different naming tasks in patients with focal epilepsy. We designed the mapping procedure to have power at both the single-patient level, necessary for clinical utility, and the group level, for validity and generalizability. For source modeling, we used the Dynamic Imaging of Coherent Sources beamformer (DICS-BF) technique, which enables the study of cortical sources of oscillatory activation in the frequency domain (Gross et al., 2001). DICS-BF has been successfully used for language lateralization in a small sample of healthy adults during a continuous recognition memory task (Passaro et al., 2011). Using DICS-BF, we examined beta-modulation effects during auditory description naming and visual picture naming in 29 patients undergoing investigations for epilepsy surgery. These are cognitively complex tasks that require The auditory description naming task (ADN) included 110 trials of naming in response to a brief spoken description. Stimulus presentation was paced by the vocal responses: after the patient responded to the current stimulus, a keyboard key controlled by the MEG engineer was used to present the next stimulus. (B) The picture naming (PN) task included 100 trials of naming in response to line drawings of everyday concrete objects. Stimuli were visually presented with a random ISI of 2-3s. stimulus recognition, concept selection, phonologic word-form selection, and motor speech preparation before speech initiation (Indefrey and Levelt, 2004). We elected to focus here on later stages of these processing cascades, such as word-form selection and motor speech preparation, which are hypothesized to be similar across the two tasks. This was accomplished by analyzing data within a 600-ms time window preceding and time-locked to the vocal response on each trial. This method has the added benefit of minimizing inter-trial variability, which should enhance the overall accuracy of localization and lateralization analyses.

Patients
MEG data were analyzed from 29 patients diagnosed with focal epilepsy. Their mean age was 42 years (standard deviation AE 13, range 21-68). There were 15 men and 14 women; all patients except one were right-handed. Prior video-EEG showed right-hemispheric seizure onset in 15 cases (52%), left-hemispheric onset in 13 (45%), and bilateral onset in 1 case (3%). Patients were enrolled prospectively, and written informed consent was obtained from all participants. Study procedures were approved by the Institutional Review Board (IRB) of the Medical College of Wisconsin.

Language tasks
During the data acquisition, participants completed two overt response language experiments: auditory description naming (ADN) and visual picture naming (PN).
The ADN task consisted of 110 trials in which a brief spoken description of an object was presented (e.g., "Stick used to hit a ball" or "Insect that makes honey") via MEG-compatible headphones (TIP-300, Nicolet Biomedical, Madison, WI). Patients were asked to respond verbally with the name of the described object, e.g., "bat" or "bee". A total of 86 different stimuli were used. To obtain an adequate number of trials, 24 of these stimuli were presented twice, for a total of 110 trials; see Supp. Table 1 Fig. 2. MEG language mapping framework. (A) MEG preprocessing. MEG data were filtered using the tSSS algorithm, bandpass filtered in a frequency range of 0.1-40Hz, and trials containing artifacts (SQUID jumps, eye-blinks, head movement, muscle) were rejected. Trials or sensors with a variance that exceeded 3 Â 10 À24 T, kurtosis larger than 15, or z-score larger than 4 were removed. Cardiac artifacts were removed from preprocessed signals using independent component analysis. (B) Average time-frequency response (TFR) analysis. Average of all sensor-level event-related power changes relative to baseline were utilized to support the selection of frequency of interests, findings from an individual patient. (C) Response-locked analysis. Left panel, trials were redefined based on vocal onsets. A 600 ms windows of post-stimulus responses, conditioned to be 100 ms before vocal onset, were selected to contrast against pre-stimulus responses. a fixed time interval of 300 ms. The right panel shows the average TFR of response-locked activations of an individual patient.
for a complete list. Stimuli were selected from a set of previously normed description naming cues (Hammeke et al., 2005). For the current study, the stimuli were audio recordings produced by a native English speaker and presented at normal listening levels (60 dB above normal hearing levels). Stimulus duration ranged from 1.3 to 2.6 s (mean AE SD ¼ 1.8 AE 0.24 s). Similar ADN tasks have been used in PET and fMRI studies (Balsamo et al., 2002;Bookheimer et al., 1998), in cortical stimulation mapping (Hamberger et al., 2003), and lesion localization studies (Pillay et al., 2017).
Stimuli were manually advanced by the clinical MEG engineer after hearing the vocal responses from the participant (Fig. 1A). After each speech response, the experimenter hit a key, followed by a programmed delay of 1200 ms (1000 ms interval of a small fixation cross and a 200 ms interval without the fixation mark) before presentation of the next stimulus. This assured that there would be no temporal overlap of cortical activation from the participant's vocalization with the following trial. A sound guide was used to direct participant vocal responses outside of the magnetically shielded room to a microphone for recording with the acquisition. To minimize artifact from eye movements, patients were asked to focus on a marked fixation point while responding to the auditory cues. Stimuli were presented using Psychtoolbox-3 (Brainard, 1977) in MATLAB 2017b (The Mathworks, Inc.) on a PC (Windows XP).
The PN task consisted of 100 images of everyday concrete objects selected from Snodgrass and Vanderwart (1980), presented with a randomly jittered interstimulus interval between 3 and 5 s. Similar investigations have been conducted using MEG (Bowyer et al., 2004). Pictures were randomly presented for 1sec followed by a blank screen. PsychoPy v3.0 (Peirce, 2007) was used to design the experiment and display pictures on a screen via a Panasonic DLP projector (PT-D7700U-K). Patients were instructed to name the presented pictures as soon as an object was displayed on the screen (Fig. 1B). Spoken responses during both PN and ADN tasks were recorded using a microphone connected to a transducer placed outside the magnetically shielded room.
The language tasks were performed as part of a typical epilepsy surgery evaluation, which included a supine resting scan to localize epileptiform activity. Patients were requested to be sleep deprived the night before and were encouraged to sleep during the initial portion of the MEG recording, a 40-60-min continuous acquisition. A small break was taken while the MEG was converted from supine to vertical positioning for language testing. Task performance was monitored by the MEG technologist to make sure the patient remained alert and responded to the target trials. Response accuracy was not measured during the two experiments, but reaction time (vocal onset relative to stimulus onset) was used to partially evaluate the overall performance of patients.

Data acquisition
Patients underwent EEG, MEG, and MRI in this study. Only MEG and MRI data were analyzed. MEG data were acquired using a 306-channel (204 planar gradiometers and 102 magnetometers) whole-head biomagnetometer system (Vectorview™, Elekta-Neuromag Ltd., Helsinki, Finland) in a magnetically shielded room (ETS-Lindgren, Eura, Finland) located at Froedtert Hospital, Milwaukee, WI, USA. The raw data were acquired with a sampling rate of 1 kHz and a high-pass filter with a 0.03Hz cutoff frequency. The position of the participant's head relative to the sensors was determined using four head-position indicators coils attached to the scalp surface. Three anatomical landmarks (nasion and pre-auricular points) and the digitized head shape were digitized using a Polhemus Fastrak system (Polhemus; Colchester, VT) for alignment with the anatomical MRI.
The anatomical MRI of each patient was acquired with a GE Healthcare Discovery MR750 3-T MR system. The high-resolution T1 image was acquired with a matrix size of 512 Â 512 Â 192 and a spatial resolution of 0.5 Â 0.5 Â 1 mm.

MEG data analysis
Response-locked beta-band power changes from 110 trials during ADN and 100 trials during PN were analyzed. Prestimulus periods were used for baseline correction.

Preprocessing and time-frequency response analysis
A temporal variant of signal space separation (tSSS) using MaxFilter software v2.2 (Elekta-Neuromag, Helsinki, Finland) was applied to remove external magnetic interferences and discard noisy sensors (Taulu and Simola, 2006). Data were epoched from À1 to 3s relative to the onset of stimuli, and bandpass filtered (Butterworth with an order of 4) from 1 to 40Hz. Trials containing artifacts (SQUID jumps, eye-blinks, head movement, or muscle) were removed by thresholding based on a variance exceeded 3 Â 10 À24 T, kurtosis larger than 15, and z-score larger than 4. Cardiac artifacts were removed using independent component analysis (ICA) based on the infomax algorithm (Bell and Sejnowski, 1995). See Fig. 2A for a summary. To improve computational efficiency the algorithm was limited to 20 components. Across all patient data, the cleaning algorithm resulted in an average of five noisy channels, ten bad trials, and two artifactual ICA components per dataset.
Voice signals were preprocessed using a 4th order 70-Hz Butterworth high-pass filter. Vocal onsets were marked using a Hilbert transform (to obtain audio envelope fluctuations) followed by an abrupt change detector (Killick et al., 2012). Marked onsets were visually inspected for false positives. Trials corresponding to voice responses faster than 400 ms (empirically determined) after stimulus onset were discarded.
Sensor-level time-frequency response (TFR) analysis was conducted to confirm the presence of beta-band modulations and verify consistency with previous literature Laaksonen et al., 2012;Pang et al., 2011;Piai et al., 2012;Youssofzadeh and Babajani-feremi, 2019). Multitaper TFR analysis was performed on all sensor locations with the following specifications: a time range of À0.5 to 2.5s and frequency range of 1-40Hz; a frequency-dependent Hanning taper with a length of three cycles (ΔT ¼ 3/f, f is the frequency of interest) and a spectral smoothing of ΔF ¼ 0.8 (f was applied to estimate Fourier representations of the trial signals). The TFRs were baseline-corrected using 0.3s pre-stimulus onset data. Before source modeling, the sensor level TFRs were used to select the frequency of interest (FOI) specific to the response data, by inspecting peak power decrease effects in a frequency range of 17-25Hz. The average of all sensor-level event-related power changes relative to baseline was used to support the selection of frequencies of interest. TFR analysis of ADN and PN task responses form a representative patient is shown in Fig. 2B.
To capture the dynamics of the language production network, a response-locked analysis of MEG trial activation was carried out. As shown in Fig. 2C (left), responses were identified from the microphone signal, and trials were re-defined using a 600 ms window extending from 700 ms before voice onset to 100 ms. Beta-power during this response period was contrasted against a 300 ms window before stimulus presentation. A TFR analysis conducted for ADN and PN task responses from a representative patient is shown in Fig. 2C (right).

MEG beamforming source analysis
MEG data were coregistered to the T1-weighted anatomical MRI using common fiducial markers and head shape digitization. T1 images were segmented using the SPM8 toolbox to extract brain tissues. A singleshell semi-realistic head model constructed from the patient's segmented MR scans was used to estimate the forward model and leadfield matrix (Nolte, 2003). To account for inter-individual variability in brain structure, atlas-based Montreal Neurological Institute (MNI) -aligned grids in individual headspace were created. The source model was defined on a regular 3D grid in MNI space; individual grids were volumetrically warped to a template grid (MNI 152 with 8-mm resolution) using a nonlinear transformation. A grid with 4464 fixed grid points inside the brain was obtained, per patient.
Beta-band source power was estimated using a frequency-domain spatial filtering technique called Dynamic Imaging of Coherent Sources (DICS) (Gross et al., 2001). DICS is a linearly constrained minimum variance beamformer in the frequency domain. It estimates the covariance matrix to calculate the spatial filter using the sensor-level cross-spectral density (CSD) matrix and applies the filter to the sensor-level CSD to reconstruct the source-level CSDs of pairwise voxel activations, providing coherence measures between the source pairs, enabling both functional connectivity and power source mapping analyses (van Vliet et al., 2018).
In this work, DICS was used for power mapping alone. Sensor-level single-trial CSDs were estimated in a beta frequency range of 17-25Hz using a fast Fourier transform (FFT) and multitapering with multiple tapers from discrete prolate spheroidal sequences (DPSS) sequence (Slepian and Pollak, 1961). Center frequency was selected at 21Hz with a spectral smoothing window of 4Hz (to have sufficiently smooth spectra) and a zero-padding of 4 s (to interpolate baseline and active windows to the same/high spectral resolution). We believe, the multitapering, padding, and smoothing effectively improve spectral power-estimatesand minimize any difference in the quality of PSDs that may arise from unequal lengths of the baseline and test data-segments, thereby reducing any related bias to beamforming source modeling. The center frequency was adjusted slightly based on the desynchrony effects observed in TFR responses of individual patients (see, 2.4.1).
The patient-specific leadfields and CSDs were used to estimate the inverse spatial filters (beamformers) for all the grid points from the combined "post" and "pre" stimulus data conditions, the so-called "common filter". The common spatial filter effectively minimizes the possibility of condition differences that may lead to biased filter estimates (Popov et al., 2018). A lambda regularization parameter for the common filter estimation was specified as 10%, to reduce sensitivity to noise and increase consistency of the spatial maps across individuals. This step resulted in a spatial filter for each volumetric source location. To estimate the source power, a fixed dipole orientation was used. The largest of the three dipole directions per spatial filter was obtained using singular value decomposition (SVD), and power corresponding to the orientation that maximizes the power was reported. Signals from both gradiometers and magnetometers were used for the source modeling. However, the different order of magnitude of the signals affects the covariance matrix used for the beamformer spatial filter computation. To account for rank deficiency of the covariance matrix due to MaxFilter, a truncation parameter, Kappa, was set to the SVD edge discontinuity of combined pre-and post-stimulus responses (range 60-70). Using the already computed spatial common filter and regularization and truncate parameter settings, "pre-" and "post-stimulus" source trials were estimated.
Statistical significance at the level of individual and group was determined using cluster-based non-parametric permutation testing (Maris et al., 2007). A one-tailed (only power decrement effect) dependent-sample t-test was applied to statistically quantify the differences in mean neural DICS beta source power between "post" and "pre" interval responses. Monte Carlo permutation testing was performed with a cluster alpha level of 0.05, 5000 randomizations, and alpha cluster level of 0.05 and extreme statistics to avoid multiple comparisons. This step resulted in t-maps for every patient. The cluster-based correction aims to improves the sensitivity of statistics. A group of neighboring sources that show similar behavior are grouped and the randomization distribution of the clusters is tested to reject the null hypothesis. Clusters in FieldTrip toolbox are identified based on spatial contiguity (inspired by the Matlab image toolbox function bwlabeln) of a 3-dimensional source volume that is binarized by the cluster-level threshold. The maximum statistics (i.e., the maximum of summed t-values across random partitions, 'maxsum') are used for controls for the expected proportion of false positives, i.e., to avoid the issue of multiple comparisons.
For a group-level analysis, patient-specific t-statistics were subjected to an independent samples t-test against a null hypothesis. Here, the goal was to determine whether the global null hypothesis can be rejected, i.e. that there is at least one cluster across all source regions which makes data non-exchangeable across patients. Similar to the patient-level statistics, a Monte Carlo permutation testing, an alpha cluster of 0.05, 5000 randomizations, and a cluster level of 0.05 with extreme statistics to avoid multiple comparisons was used for group-level statistics.

Laterality and ROI analyses
A volumetric automated anatomical labeling (AAL) atlas (Tzourio-Mazoyer et al., 2002) consisting of 116 subdivisions of cortical, subcortical, and cerebellar regions (parcels) was employed to combine and summarize power source measures (t-values) across regions. Parcellation was conducted at patient and group levels, to report laterality indices (LIs) and the most strongly activated regions (ROIs), respectively.
A laterality or asymmetry index, ranging from À1 to 1 was calculated using the formula, LI ¼ (L À R)/(L þ R), where L and R are t-statistic source power values obtained for 13 left and right hemisphere regions, respectively. Source power values for a given ROI were obtained by averaging the t-values of grid points within the ROI. Following the approach suggested by previous MEG studies (Findlay et al., 2012;Papanicolaou et al., 2004;Raghavan et al., 2017;Tanaka et al., 2013), LI indices greater than 0.10 were considered left hemispheric dominance, less than À0.10 right hemispheric dominance, and intermediate values as symmetric activation.
A Pearson linear correlation analysis was conducted to assess the consistency of localized sources in the 26 (13 Â 2) selected parcels between the ADN and PN tasks. Pearson correlation was examined with an FDR threshold value of q ¼ 0.05. Laterality indices of PN and ADN were also examined using the Pearson correlation, and ROIs with highest laterally scores were reported. Fig. 3. Language-specific parcels. 13 frontotemporal, parietal and sensorymotor cortical regions were selected for language localization and laterality analysis. Random colored areas represent parcellated brain regions according to the AAL anatomical labeling scheme (Tzourio-Mazoyer et al., 2002  Localized sources corresponding to ADN and PN tasks identified by the DICS beamformer source analysis (post-vs. pre-stimulus intervals). t-values were derived from a between-subject independent sample t-test against a null hypothesis. Suprathresholded t-values between 4 and 6 are shown. (B) AAL group source parcellation. For comparison purposes, source activations are displayed on an inflated template (MNI-152) with dark representing sulci and gray representing gyri. Suprathresholded values with t-values greater than 5 are shown.

Results
Patients tolerated the MEG language testing well, and no adverse events were reported. Patients' response times to the ADN and PN tasks were 2230 AE 556 and 897 AE 648 msec (mean AE SD), respectively. Time interval from vocal responses to presentation of the next stimulus for all trials during the ADN task was, 2147 AE 383 msec (mean AE SD) ranging 1351-2998 msec, and during the PN task was, 3074 AE 397 msec (mean AE SD) ranging 1944-3825 msec. These support no overlap of neural activity from the prior response into the following trial. Fig. 4 shows group localized sources for ADN and PN tasks based on the beta (17-25Hz) power changes, as identified by the response-locked DICS beamformer analysis and permutation statistics. Findings revealed decreases in beta power predominantly in canonical left frontal language, left superior temporal, and left sensorimotor (pre-and postcentral gyri) regions during both tasks, and in right sensorimotor regions during PN.
A whole-brain AAL parcellation analysis revealed strong modulation of beta in the left prefrontal cortex areas (IFG pars opercularis t ¼ 5.8 and IFG pars triangularis, t ¼ 5) during ADN, whereas PN showed widespread power decrement in left temporal, parietal and frontal cortices (supramarginal gyrus, t ¼ 5.1, superior temporal gyrus t ¼ 4.7, postcentral t ¼ 5.2, and IFG pars triangularis, t ¼ 4.9) as well as right precentral gyrus, t ¼ 5.1. Regions with t-values > 4 are reported in Table 1.
A Pearson correlation analysis was conducted to assess the consistency of beta power decrement within the 26 ROI parcels across the two tasks (see section Laterality and ROI analyses). The analysis showed a strong correlation of r ¼ 0.72 (p < 0.05, FDR corrected) between ADN and PN tasks (Fig. 5A). Activation peaks (highest source values across 26 parcels for each patient's data) projected on a cortical surface template supported an overall consistency in localization between ADN and PN task responses, as shown in Fig. 5B.

Discussion
We developed an analysis pipeline that characterizes the changes in beta-band oscillatory activity within language-related brain circuits. Using a frequency-resolved beamforming source analysis in MEG, our pipeline showed prominent beta power decrements in the left hemisphere during auditory and visual naming tasks in patients undergoing evaluations for epilepsy surgery. The beta power decrements localized to canonical language and motor speech production areas at the individual patient level and showed a consistent left-hemispheric dominance at the group-level. Source-modeling of task-related beta power decrements may, therefore, be suitable for presurgical evaluations of hemispheric language dominance.
Beta-band power decrement has been reported to be predominantly observed in the left inferior frontal and middle frontal cortices during overt and covert language tasks, e.g., auditory word/verb generation, spoken sentence comprehension, and visual picture/object naming, in both healthy and clinical populations (Armeni et al., 2019;Findlay et al., 2012;Fisher et al., 2008;Hirata et al., 2004;Lewis and Bastiaansen, 2015;Liljestr€ om et al., 2015;Piai et al., 2012;Wang et al., 2012;Weiss and Mueller, 2012). Consistent with previous findings, our analysis of ADN and PN task responses showed decrements of beta frequency power within canonical language-related areas of the left frontal lobe and also within premotor, somatosensory, and auditory regions linked to motor speech production (Tourville et al., 2019), as in Fig. 4 and Table 1. The sources of the beta-band power decrease estimated in these frontal and temporal regions have specifically been linked to word production operations (Baggio and Hagoort, 2011;Bastiaansen and Hagoort, 2015;Lewis and Bastiaansen, 2015;Wang et al., 2012). Our time-frequency analyses at the individual patient level (e.g., Fig. 2C) also support the presence of robust power decrements in the range of beta-band power changes. As for the contributions of extracanonical areas, decreased power of the motor-cortical beta rhythm has been observed during the activation of the motor system (which is also necessary for vocalization during our naming tasks), and during its secondary involvement by language stimuli such as action verbs (Saarinen et al., 2006;van Elk et al., 2010;Weiss and Mueller, 2012).
The two experiments used in our study were not identical, thus differences in the patterns of activations were expected. The ADN task uses auditory descriptions of objects as stimuli, whereas the PN task uses pictures of everyday objects, resulting in a substantial difference in the timing of presentation of stimuli across the tasks (see section Language tasks). This difference was reflected in the time-frequency responses (e.g. Fig. 2C) as well as in the level of overlap between language and motor speech relevant areas (Fig. 4). Both tasks showed a significant decrease in beta power in the left IFG, although these effects were somewhat stronger for the ADN task. These findings agree with several M/EEG and fMRI studies supported strong contributions by the IFG (BA44/45, pars opercularis/pars triangularis) during semantic or orthographic verbal fluency or word generation tasks (de Guibert et al., 2010;Friederici, 2011;Le Bihan et al., 2012;Leh ericy et al., 2000;Youssofzadeh et al., 2018), and with numerous fMRI studies showing the involvement of the left inferior frontal cortex in word and concept selection (Badre et al., 2005;Novick et al., 2005;Thompson-Schill et al., 1997;Wagner et al., 2001). The somewhat stronger prefrontal activation for ADN compared to PN may reflect the greater working memory and other executive demands involved in processing the ADN stimuli, which consisted of strings of words that must be individually recognized and then combined semantically and syntactically. In contrast, the PN responses showed more widespread beta modulations in bilateral sensorimotor regions associated with motor speech control. Although the reason for this difference is unclear, we speculate that in the case of picture naming, motor speech processes may be more tightly time-locked to speech output than in the ADN task. That is, the much slower stimulus presentation in the ADN task might allow for a slower and less uniform time course of activation of motor control processes across trials, resulting in a weaker average signal in these regions. Despite these differences, the group activation patterns for the two tasks were highly correlated (r ¼ 0.74), Fig. 5A. This similarity was likely enhanced by our use of a response-locked time window for analysis, which emphasizes later stages of naming and minimizes sensitivity to earlier activity in sensory processing systems. Both tasks ultimately engage modality-independent language regions for semantic processing, word retrieval, and articulatory planning related to the naming task (Price, 2012).
Our analysis pipeline has some key advantages over conventional language mapping techniques. First, as mentioned above, the responselocked analysis captures the dynamics of language production-related brain activities and minimizes contributions from prelinguistic sensory systems. The suitability of response-locked analysis has been demonstrated for MEG auditory verb-generation and EEG picture-naming tasks (Findlay et al., 2012;Piai et al., 2014;Ri es et al., 2013). Second, we used a beamforming frequency-based source modeling method, DICS-BF, to obtain source estimates of oscillatory changes related to language processing. DICS-BF is especially suitable for analyzing induced activity and oscillatory dynamics (Liljestr€ om et al., 2005). Third, almost all MEG source estimates are vulnerable to the unavoidable phenomena of volume conduction (forward model) and instantaneous field spread (inverse model) (Drakesmith et al., 2013;Palva et al., 2018;Wang et al., 2018). To reduce the effect of volume conduction, we employed a semi-realistic shaped volume conductor model (Nolte, 2003) and computed volume conduction models based on patients' anatomical information to improve the reliability of source estimates. We also analyzed the data at the level of the underlying sources by computing the beamformer filter. Source-space beamformer effectively mitigates the interpretational difficulties introduced by electromagnetic field spread (Gross et al., 2013;Schoffelen and Gross, 2009). Additionally, our mapping approach relied on changes in oscillatory strength between pre-and (response-locked) post-stimulus data conditions, rather than estimating the changes in one condition alone. This differential analysis is less sensitive to background noise and does not underestimate deep sources (Quraan et al., 2011). Fourth, most previous MEG studies using beamformer methods have reported group activation results based on simple averages without any attempt to form a statistical analysis of voxel and cluster significance. However, statistical methods are essential if meaningful comparisons are to be made between the results for different participant cohorts or modalities. In our analysis pipeline, we employed both individual-level and group-level statistics to minimize the confounding effects that could arise from power differences between data conditions and individual variability. Overall, our analysis pipeline is data-driven, optimal time and frequency parameters are selected based on the participant's time-frequency responses, and "power maps" are derived from oscillatory activations supported by statistical measures, thus potentially providing better reliability and validity for language mapping than conventional methods.
There is no consensus on an optimal MEG methodology for assessing language dominance during presurgical evaluations for epilepsy or tumor surgery. The conventional dipolar source-modeling method still in common use is suboptimal for modeling cortical activity from spatially extensive language networks (Huang et al., 2016;Raghavan et al., 2017). Source mapping of evoked activations may not always lead to a full delineation of the language network (Doesburg et al., 2015;Hirata et al., 2004;Laaksonen et al., 2012;Piai et al., 2014), perhaps because the oscillatory changes are not taken into account, and the averaging of evoked responses is ineffective if responses are not phase-locked across trials (Piai et al., 2014). Our group laterality analysis yielded a consistent left-hemispheric dominance for the ADN and PN task responses (Fig. 5), in agreement with previous MEG studies that found beta power decrements to be a robust marker for language lateralization (Findlay et al., 2012;Fisher et al., 2008;Hirata et al., 2010), lending further support to the beta modulation approach for determining language dominance in MEG clinical applications.
Although our framework has been designed for clinical investigations, it can be easily generalized to characterize language network dynamics in healthy controls. Specifically, the approach can be used to interrogate any arbitrarily defined brain network for patterns of connectivity (information flux within a theoretically derived network defined by anatomical boundaries), to gain a better understanding of network function. Some recent MEG studies using functional connectivity measures have emphasized the importance of dynamic interactions between different regions of the brain during language comprehension and production to explain the function of these language networks (Schoffelen et al., 2017;Youssofzadeh et al., 2017). For example, oscillatory activity underlying language comprehension has been suggested to involve reciprocal interactions originating from frontal and parietal regions, with interactions peaking at the beta-band frequency (Schoffelen et al., 2017). MEG patterns from a network connectivity analysis revealed regions involving in a covert verb generation task in children, together with increased and left-lateralized connectivity with age (Youssofzadeh et al., 2017).
It is important to note that the DICS-BF approach requires the specification of both time and frequency ranges for computation. In our study, we epoched our data to a time window thought to be relevant to word production during naming and restricted analyses to the beta frequency band. The choice of temporal and spectral window to assess other aspects of language function is not trivial, especially for sentence-level language comprehension. With additional experience, we anticipate proposing guidelines so researchers can make informed decisions about choosing parameters for optimal language mapping. One concern with our methodology is the choice of unequal window lengths, i.e., comparison of a longer 600 ms active time window against a shorter 300 ms baseline period. While multitapering, padding, and smoothing techniques effectively improve spectral factorization (smooth PSD estimates with same/ high spectral resolution) of the data pieces, there remains the possibility that the PSD variance-estimates for two windows are unequal. Opting for identical window lengths is, therefore, recommended for similar investigations in the future. As another limitation of the current study, we did not validate findings against other modalities. Validation against fMRI language mapping should be attempted in the future. By measuring the level of concordance with fMRI and possibly comparing against Wada testing, we hope to provide improved sensitivity and specificity of presurgical mapping of language functions in patients with epilepsy. Finally, our analysis framework was designed for volumetric source analysis, but it can be extended to surface-based source analysis. As a preliminary investigation, a surface-based DICS-BF (implemented in Brainstorm toolbox, see section Code and data availability) with sources that are constrained to the individual's cortical surface showed a high consistency of spatial localization with the volumetric sources (implemented using Fieldtrip toolbox). A sample result from a patient completing the ADN task is shown in Fig. 6, supporting the generalizability of the method.
To summarize, source analysis of MEG beta power changes is a powerful method for lateralizing and localizing task-related cortical engagement. Findings based on a DICS-BF source analysis revealed regions involved in language and motor speech processing and showed consistent left-hemispheric dominance of these networks in patients with focal epilepsy during two naming tasks. The proposed pipeline may be suitable for mapping the language cortex in presurgical evaluations.

Ethics statement
The experimental protocol was approved by the Institutional Review Board (IRB) of the Medical College of Wisconsin (MCW). All participants gave written informed consent, and the approval process of the IRB complies with the declaration of Helsinki.
Code and data availability MEG and MRI data were analyzed in FieldTrip toolbox v20190419 (fieldtriptoolbox.org), activation patterns were visualized using CONN toolbox ver.17c (web.conn-toolbox.org), All analyses were completed in MATLAB2018a (The Mathworks, Inc.). A surface-based version of the pipeline was implemented in Brainstorm toolbox (http://neuroimage. usc.edu/brainstorm/) and is publicly available at, https://github.co m/vyoussofzadeh/DICS-beamformer-for-Brainstorm.Anonymized MEG and structural MRI data, as well as custom code, are available on demand by emailing the corresponding author.