Clinical Validation of the Champagne Algorithm for Evoked Response Source Localization in Magnetoencephalography

Magnetoencephalography (MEG) is a robust method for non-invasive functional brain mapping of sensory cortices due to its exceptional spatial and temporal resolution. The clinical standard for MEG source localization of functional landmarks from sensory evoked responses is the equivalent current dipole (ECD) localization algorithm, known to be sensitive to initialization, noise, and manual choice of the number of dipoles. Recently many automated and robust algorithms have been developed, including the Champagne algorithm, an empirical Bayesian algorithm, with powerful abilities for MEG source reconstruction and time course estimation (Wipf et al. 2010; Owen et al. 2012). Here, we evaluate automated Champagne performance in a clinical population of tumor patients where there was minimal failure in localizing sensory evoked responses using the clinical standard, ECD localization algorithm. MEG data of auditory evoked potentials and somatosensory evoked potentials from 21 brain tumor patients were analyzed using Champagne, and these results were compared with equivalent current dipole (ECD) fit. Across both somatosensory and auditory evoked field localization, we found there was a strong agreement between Champagne and ECD localizations in all cases. Given resolution of 8mm voxel size, peak source localizations from Champagne were below 10mm of ECD peak source localization. The Champagne algorithm provides a robust and automated alternative to manual ECD fits for clinical localization of sensory evoked potentials and can contribute to improved clinical MEG data processing workflows.


Introduction
Planning resective brain surgery, whether for the removal of structural lesion or of seizure onset zone, requires the estimation of the location of functional regions of cortex in order to plan a strategy that maximizes benefit while minimizing risk of postsurgical functional deficit. Magnetoencephalography (MEG) is particularly well-suited for functional preoperative brain mapping prior to surgery of tumor or vascular lesion because its results do not depend primarily on blood flow that can be altered by tumor growth and does not suffer from susceptibility artifacts and vascular confounds seen, e.g., with fMRI (Kreidenhuber et al. 2019). Like other neurophysiologic methods, it has exquisite temporal resolution, and with the choice and application of appropriate source imaging techniques, can provide results with a high degree of spatial resolution. This is one of several papers published together in Brain Topography on the "Special Issue: Computational Modeling and M/EEG. The clinical standard for MEG source localization of interictal epileptiform discharges (IEDs, or "spikes") as well as for functional landmarks from sensory evoked responses is the equivalent current dipole (ECD) localization algorithm (Burgess et al. 2011). Despite its widespread use, the technique of manual ECD fitting is subjective, labor-intensive, and sensitive to noise. Many semi-and fully automated algorithms have been elaborated to counter these difficulties. These have included, for example, the application of the coherent Maximum Entropy of the Mean (cMEM) approach as described in Chowdhury et al. (2013) to the localization of IEDs, a method termed distributed magnetic source imaging or dMSI, and found to be successful in comparison to ECD fitting by Pelligrino et al. (2018). Recently, the Bayesian multidipole iterative Monte Carlo approach "SESAME" described by Sommariva and Sorrentino (2014) was evaluated in comparison to standard ECD methods by Luria et al. (2020) for the localization of IEDs, likewise showing excellent performance. These methods were also noted to be advantageous because results were more operator-independent.
Over the last decade we have developed and gained experience with an empirical Bayesian algorithm with powerful abilities for MEG source reconstruction and time course estimation (Wipf et al. 2010;Owen et al. 2012); in its current form, it is called Champagne. Here, we evaluate the performance of the Champagne algorithm in preoperative brain tumor patients for localization of functional cortices (sensory and auditory) where there was minimal failure in localizing sensory evoked responses using standard ECD methods. Our aim is to establish equivalence between the clinical standard of manual ECD fitting and automated Champagne in this scenario, and also to illustrate several cases where Champagne may provide specific benefit in the case of ECD failure.

Patient Characteristics
A retrospective analysis of our clinical MEG database identified 21 patients who underwent preoperative mapping of auditory and somatosensory cortices with magnetic source imaging between March 2018 and May 2018 before tumor resection at the University of California, San Francisco. Please refer to Table 1 for the demographic, clinical, and pathological characteristics of each patient. The final cohort Tasks During the MEG procedure, each individual lay supine on a comfortable bed. For the auditory task, a 1 kHz tone was presented binaurally for 400-msec duration with an interstimulus interval of 2000-msec and a jitter between 0 to 100-msec at random for a total of 120 trials. Data was collected in epochs of 600-msec with a 100-msec pre-stimulus interval for each trial. For the somatosensory task, lip and index fingers were painlessly stimulated using pneumatically driven pulses (20-25 psi, 40 msec duration) from balloon diaphragms. Epochs of 400-msec duration with a 100-msec pre-stimulus interval were collected for 512 trials, comprising 256 trials at each stimulus site: right and left lip sites (RLip, LLip) and right and left second digit (RD2, LD2). 512 trials of stimulus with an inter-stimulus interval of 500msec and a jitter of 50-msec were randomly presented in one block. Two somatosensory blocks were performed in order to test a total of four somatosensory sites (RD2, RLip, LD2, and LLip). The location of the patient's head position relative to the MEG sensors was determined at the beginning and ending of the collection by means of three small fiducial marks placed at nasion, left preauricular, and right preauricular points. The pre-and post-trial locations were used to measure head movement during each of the tasks (Kirsch et al. 2007;Nagarajan et al. 2008;Traut et al. 2019).

MEG Preprocessing
MEG sensor data for each task were averaged across trials after removing any trials with eye blink, muscle activity, or other obvious artifact. As a result, each subject had one averaged auditory evoked field dataset and four averaged somatosensory evoked field datasets (RD2, RLip, LD2, LLip).
All MEG datasets were 3rd order gradient denoised and detrended. Next, the clinical MR images from a 3-tesla unit (GE Medical Systems) were transferred to MEG workstations and used in conjunction with MEG localization information for coregistration. Coregistration was done using CTF software (CTF Omega 2000, CTF MEG, Coquitlam, BC, Canada) and coregistration errors passed a maximum threshold of 0.5 cm.

Equivalent Current Dipole
The equivalent current dipole localization algorithm utilizes an iterative least squares minimization to compute the strength and location of dipoles in a single spherical volume of uniform conductivity that could account for the sensory data. The auditory evoked field data were filtered with a high pass between 1 and 4 Hz and for most cases filtered with a low pass of 40 Hz. The somatosensory evoked field data were filtered with a high pass between 1 and 15 Hz and for most cases filtered with a low pass of 40 Hz. The primary auditory cortex in each hemisphere was localized based on the M100 auditory response peak using ECD fit (Pross et al. 2015;Chang et al. 2016). The primary somatosensory cortex corresponding to each of the tactile stimulation sites was localized based on the earliest response peak ~ 45 msec poststimulus onset seen in the contralateral hemisphere. Dipole fits were accepted based on a goodness-of-fit values > 0.95 and having 95 % confidence volume of reconstructions < 0.1 cm 3 .

Champagne Source Reconstruction
The Champagne algorithm is a robust empirical Bayesian method to estimate brain source localization and the time course of multiple neural sources of MEG data (Wipf et al. 2010;Owen et al. 2012). To apply Champagne, first the clinical MR images were spatially normalized to the Montreal Neurological Institute (MNI) atlas template using SPM (http:// www. fil. ion. ucl. ac. uk/ spm). Next, we calculated a three-component lead field for each voxel (e.g. see Huang et al. 1999) using a forward model of sensor activity at a spatial resolution of 8mm spanning the entire brain by using the NUTMEG open-source analysis toolbox in MATLAB (Dalal et al. 2008;Dalal et al. 2011;Owen et al. 2012). We used the NUTMEG software for visualization of source time-course, overlaid on coregistered MRIs. The time-series and leadfield calculations were estimated for voxels in the individual MRI space (co-registered to the MEG data) and then transformed to MNI coordinates. Spatiotemporal neural activity from both hemispheres following stimulus presentation was reconstructed. The activation patterns were then displayed on the MNI brain template for each subject.
In order to determine activation patterns, Champagne was deployed on every averaged dataset across each subject. The input to the automated Champagne algorithm were: (1) the MEG data, selected for the active time-window of interest; (2) the lead field matrix calculated from the forward model; and (3) baseline MEG data (pre-stimulus) used to calculate the model noise covariance. The only free parameter was the number of iterations used for the algorithm convergence, typically set to 100, based on experience with both simulations and with many real datasets. The algorithm output is the time course of source activity for each orientation, and source variance (power) for each voxel.
For comparison with ECD, we recorded the maximum of the estimated source activity within the active time window of interest. For the auditory evoked fields, we used a baseline time-window of − 100 to −5 ms pre-stimulus onset and an active time-window of 25-250 ms post-stimulus onset.For the somatosensory evoked fields, we used a baseline timewindow of − 100 to − 5 ms pre-stimulus onset and an active time-window of 5-105 ms post-stimulus onset.

Comparison of Results Across Algorithms
Both the Champagne peaks and ECD localizations were chosen based on waveforms at consistent latencies yielding successful source localization results in the perirolandic region of the brain for the somatosensory evoked responses and in the auditory cortex for the auditory evoked responses. If a peak or dipole could not be found in these regions it was considered a failed localization. We calculated the success rate of each pipeline for each task across the cohort. Here, the success rate is defined as the number of successfully localized sources as described above out of the number of potential localizable sources for the auditory and somatosensory tasks across the cohort using each analysis.
For each stimulus, tri-planar MNI coordinates of the corresponding dipole fit and of the corresponding Champagne peak were recorded. As expected, the auditory task resulted in two sources from each hemisphere, while the somatosensory task resulted in a single source in the hemisphere contralateral to the stimulus presentation. Ideally, datasets for each subject thus included 12 total MRI coordinates, one from ECD fit and one from Champagne, for each stimulus: right AEF, left AEF, RD2, LD2, RLip, and LLip. Next, the Euclidean distances were calculated between ECD-Champagne coordinate pairs for each stimulus across subjects. The effectiveness of the Champagne algorithm in localizing source activity in this clinical population was determined based on how close the Champagne peaks were to the corresponding dipole fits for each stimulus. The cases in which a dipole location or Champagne peak could not be found were excluded from this analysis.

Results
The overall success rate for source localization was greater for Champagne (96.8 %) than for ECD analysis (93.7 %). Table 2 provides a breakdown of the success rate of source localization for each stimulus type and analysis method across the cohort.
Across stimulus types and modalities, the average distance between Champagne peak and ECD location was 9.3 ± 4.8 mm. Figure 1 a violin plot showing the distances between Champagne peaks and ECD locations across each stimulus, provides a more detailed breakdown of the performance of Champagne analysis compared to ECD analysis for each category. Across the entire cohort, the average distance between Champagne peak and ECD location was 9.1 ± 5.7 mm for the right auditory evoked field and 10.1 ± 5.3 mm for the left auditory evoked field. The average distance between Champagne peak and ECD location was 11.9 ± 8.4 mm for the right index (RD2) somatosensory evoked field and 8.5 ± 4.1 mm for the left index (LD2) somatosensory evoked field. The average distance between Champagne peak and ECD location was 10.7 ± 5.7 mm for the right lip (RLip) somatosensory evoked field and 7.1 ± 3.3 mm for the left lip (LLip) somatosensory evoked field task. Overall, across these various stimulus types and modalities, the average distance between Champagne peak and ECD location was 9.3 ± 4.8 mm.
We present cases that had high concordance between locations of Champagne peaks and respective ECD locations for each category of stimulus. These cases are shown in Figs. 2 and 3, and 4. Figure 2 is for Patient 14 who had concordant Champagne and ECD localization of auditory evoked fields (AEFs). Figure 3, for Patient 13, shows index finger somatosensory evoked field (SEFs) (top row  Figure 4, for Patient 2, shows concordance between Champagne and ECD fit for lip SEFs (top row right lip stimulus and bottom row left lip stimulus); note that Champagne was successful in spite of a relatively noisy background in baseline signal.

Discussion
We compared the performance of Champagne algorithm to the clinical standard for MEG source localization of sensory evoked fields, equivalent current dipole (ECD) fit. The Champagne peaks were on average < 10 mm away from corresponding ECD localization. We found strong concordance between Champagne and ECD for both auditory and somatosensory tasks throughout our clinical cohort. Despite the success of Champagne algorithm, there are certain drawbacks that we predict might result in algorithm failure. Specifically, it is possible that the noise is not adequately modelled in the baseline due to artifactual data collection. In addition, there may be large co-registration errors that may result in faulty localization. We have circumvented such problems through standard troubleshooting steps that mitigate the effects of such issues.
In a few of the cases, Champagne was successful in localizing a source despite ECD failing to localize a source: this was true for three left and one right auditory evoked responses, and two left lip somatosensory responses. This provided insight into Champagne's ability to provide localization even when ECD was unable to perform and provided some preliminary indication that Champagne algorithm was more robust than ECD fit in some cases. Further testing is still required to determine the sensitivity and specificity of Champagne against direct electrocortical stimulation mapping, a validated functional localization method that is considered to be clinical gold standard.
A tool like Champagne could also be especially useful in improving efficiency in MEG clinical data processing workflows. Current data processing workflows, especially with presurgical mapping data from brain tumor patients, can be taxing, labor intensive, and sensitive to error. Due to its automated procedures, we hope that the Champagne algorithm can be a useful and effective resource for source for a wide range of users. For example, since Champagne is an automated tool, it can be used as a first-pass source localization of sensory evoked fields. Champagne algorithm can also provide localization where ECD fitting is difficult or where it fails due to complexity of source activity, noisy data, or error of model. As noted above, as other groups have found for IEDs, algorithms that are robust in the face of technical limitations and that are usable for a wider range of users can provide practical solutions to clinical MEG source localization. Fig. 1 Distance between Champagne and ECD localization. Violin plot showing the difference in millimeters (y-axis) between Champagne peak activity and ECD locations for each stimulus type (x-axis). R_AEF right auditory evoked field; L_AEF left auditory evoked field; RD2 right second digit somatosensory evoked field; LD2 left second digit somatosensory evoked field; RLip right lip somatosensory evoked field; LLip left lip somatosensory evoked field and right (bottom) hemispheres. The peak activity is at roughly 100 msec, consistent with the M100 auditory cortical response. Column C shows snapshots of brain MR slices in coronal and axial planes with the Champagne peak activity shown in red, corresponding to the peaks marked with red vertical cursors in column B. Column D overlays both Champagne peak (red circles) and ECD fit (green squares) on brain MR slices in coronal and axial plane . Column A shows the time series of each sensor's activity across time (ms). The top graph in column A overlays all left hemisphere sensors and the bottom graph in column A overlays all right hemisphere sensors. Column B displays the time course of cortical activity at the specific voxel corresponding to peak auditory cortical signal in the left (top) and right (bottom) hemispheres. The peak activity is between 53 and 57 ms, consistent with an index finger somatosensory cortical response. Column C shows snapshots of brain MR slices in the axial plane with the Champagne peak activity shown in red, corresponding to the peaks marked with red vertical cursors in column B. Column D overlays both Champagne peak (red circles) and ECD fit (green squares) on brain MR slices in coronal and axial plane Column B displays the time course of cortical activity at the specific voxel corresponding to peak somatosensory cortical signal in the left (top) and right (bottom) hemispheres. The peak activity is between 21 and 23 msec, consistent with a lip somatosensory cortical response. Column C shows snapshots of brain MR slices in the axial plane with the Champagne peak activity shown in red, corresponding to the peaks marked with red vertical cursors in column B. Column D overlays both Champagne peak (red circles) and ECD fit (green squares) on brain MR slices in coronal and axial plane 1 3

Conclusion
In this study, we examine Champagne's performance for source localization of sensory evoked fields in a clinical cohort. We provide strong evidence that the Champagne algorithm is equivalent to ECD fit. We also highlight Champagne as a viable alternative to the commonly clinically used ECD algorithm.
Funding This work was funded in part by the following Grants: UCOP-MRP-17-454755, NIH grants R01EB022717, R01NS100440, R01DC176960, R01DC017091, R01AG062196, and a research contract from Ricoh MEG Inc.
Availability of data and material deidentified neurophysiologic data available upon request.
Code availability: the software described is freely available as part of the open-source Neurodynamic Utility Toolbox for Magnetoencephaloand Electroencephalography (NUTMEG).

Conflict of interest None.
Ethics approval all procedures were approved by the UCSF Institutional Review Board.

Consent to participate all subjects provided consent for participation in research.
Consent for publication all subjects provided consent for publication of deidentified information.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.