Some individuals may have both hands amputated, in which case they have no hand with which to grasp a force-sensing device. In such cases, the system can make assumptions about the finger flexion forces at certain time points. For example, when the amputee is asked to concentrate the force on a specific finger, the system can assume an arbitrarily higher flexion force for that finger and a minimal force on the other fingers during that time frame. Amputees also can be prompted to vary the total grip force of a grasp where in some time frames the amputee can be prompted to exert a minimal total grasp force, and in other time frames the amputee can be prompted to exert a high grasping force. For each time frame, the system can make an assumption about the total grasp force and assume how the force would be distributed across the fingers during a grasp such as a regular closed fist depending on previous data from able-bodied individuals. For each case, the opposing thumb force can be assumed to be approximately equal to the summation of the four finger forces. An advantage of this approach is that it requires the same effort across subjects for exerting specific grasping forces because the assumption made of the grasping force for each time frame is independent of the magnitude of EMG signal (because of normalization) and of the overall muscular strength of the subject.

## CONCLUSION

We have demonstrated proof of principle in the use of five EMG electrodes for predicting hand pose and individual finger forces using a one-against-one SVM and RFR, respectively. The present system has potential for myoelectric control of dexterous hand

### REFERENCES

prostheses. Future work should explore additional methods of feature selection, signal filtering, machine learning classification, Kalman filtering, and training. New myoelectric control systems should be adjustable with small amounts of new training data without the need to retrain the entire system so that drifts in the EMG signal over time do not decrease classification accuracy.

### ETHICS STATEMENT

Subjects gave informed written consent to participate in the experiments. The experiments were approved by the Institutional. Review Board at Arizona State University (Protocol: #1201007252) and were in accordance with the Declaration of Helsinki. Prior the experiment, participants were given time to read a consent form that described the experimental procedure. After reading the form, participants were asked if they had any questions. After any questions were answered, participants were asked to sign and date the consent form. Next, the experimenter signed and dated the consent form. No participants from vulnerable populations were used in this study.

### AUTHOR CONTRIBUTIONS

AG designed the prosthetic control system, designed most of the experimental protocol, carried out the experimental protocol, collected data, analyzed data, and made data plots. PA provided expertise in the area of myoelectric control of robotic devices, thereby playing a role in design of the prosthetic control system and experimental protocol design. MS provided expertise in motor control and neuroscience. He contributed to the experimental design and data analysis methods. All authors contributed to manuscript preparation.

### ACKNOWLEDGMENTS

The authors thank Dr. Justin Fine for assistance with statistical analysis.

### FUNDING

Research reported in this publication was supported by the Eunice Kennedy Shriver National Institute of Child Health and Human Development of the National Institutes of Health under Award Number R21HD081938 and the Grainger Foundation.

**Discalimer:** The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer AA and handling editor declared their shared affiliation and the handling editor states that the process nevertheless met the standards of a fair and objective review.

*Copyright © 2017 Gailey, Artemiadis and Santello. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Prediction of STN-DBS Electrode Implantation Track in Parkinson's Disease by Using Local Field Potentials

Ilknur Telkes 1 , Joohi Jimenez-Shahed2 , Ashwin Viswanathan3 , Aviva Abosch4 and Nuri F. Ince1 \*

*1 Clinical Neural Engineering Lab., Biomedical Engineering Department, University of Houston, Houston, TX, USA, 2 Department of Neurology, Baylor College of Medicine, Houston, TX, USA, 3 Department of Neurosurgery, Baylor College of Medicine, Houston, TX, USA, 4 Department of Neurosurgery, University of Colorado, Aurora, CO, USA*

Optimal electrophysiological placement of the DBS electrode may lead to better long term clinical outcomes. Inter-subject anatomical variability and limitations in stereotaxic neuroimaging increase the complexity of physiological mapping performed in the operating room. Microelectrode single unit neuronal recording remains the most common intraoperative mapping technique, but requires significant expertise and is fraught by potential technical difficulties including robust measurement of the signal. In contrast, local field potentials (LFPs), owing to their oscillatory and robust nature and being more correlated with the disease symptoms, can overcome these technical issues. Therefore, we hypothesized that multiple spectral features extracted from microelectrode-recorded LFPs could be used to automate the identification of the optimal track and the STN localization. In this regard, we recorded LFPs from microelectrodes in three tracks from 22 patients during DBS electrode implantation surgery at different depths and aimed to predict the track selected by the neurosurgeon based on the interpretation of single unit recordings. A least mean square (LMS) algorithm was used to de-correlate LFPs in each track, in order to remove common activity between channels and increase their spatial specificity. Subband power in the beta band (11–32 Hz) and high frequency range (200–450 Hz) were extracted from the de-correlated LFP data and used as features. A linear discriminant analysis (LDA) method was applied both for the localization of the dorsal border of STN and the prediction of the optimal track. By fusing the information from these low and high frequency bands, the dorsal border of STN was localized with a root mean square (RMS) error of 1.22 mm. The prediction accuracy for the optimal track was 80%. Individual beta band (11–32 Hz) and the range of high frequency oscillations (200–450 Hz) provided prediction accuracies of 72 and 68% respectively. The best prediction result obtained with monopolar LFP data was 68%. These results establish the initial evidence that LFPs can be strategically fused with computational intelligence in the operating room for STN localization and the selection of the track for chronic DBS electrode implantation.

Keywords: local field potentials, subthalamic nucleus, microelectrode recordings, least mean square algorithm, LDA classification

#### Edited by:

*Paolo Bonifazi, Tel Aviv University, Israel*

#### Reviewed by:

*Alessandro Stefani, University of Rome, Italy J. Luis Lujan, Mayo Clinic, USA*

> \*Correspondence: *Nuri F. Ince<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Neuroprosthetics, a section of the journal Frontiers in Neuroscience*

Received: *29 December 2015* Accepted: *21 April 2016* Published: *09 May 2016*

#### Citation:

*Telkes I, Jimenez-Shahed J, Viswanathan A, Abosch A and Ince NF (2016) Prediction of STN-DBS Electrode Implantation Track in Parkinson's Disease by Using Local Field Potentials. Front. Neurosci. 10:198. doi: 10.3389/fnins.2016.00198*

## INTRODUCTION

Deep brain stimulation (DBS) of the subthalamic nucleus (STN) is an effective therapy for the treatment of the motor symptoms of Parkinson's disease (PD) (Herzog et al., 2004; Hariz, 2012). DBS surgery involves localization of the motor territory of the STN, for permanent implantation of a DBS electrode at this site. Although the exact mechanism of DBS remains to be elucidated, STN stimulation is well-tolerated and improves all of the cardinal symptoms of PD (Levy et al., 2002). However, STN stimulation can result in side effects arising from the spread of stimulation to structures surrounding the STN (Richardson et al., 2009). Moreover, sub-optimal positioning of DBS electrodes accounts for up to 40% of cases of inadequate efficacy of stimulation postoperatively (Okun et al., 2005). Thus, developing quantitative electrophysiological methods to define the optimal site of stimulation may help optimize DBS outcomes.

The task of the neurosurgeon is to place the DBS electrode within the motor territory of the STN, and well within the STN borders such that current does not spread to the surrounding structures, thereby resulting in stimulation-limiting side effects (Richardson et al., 2009). Although the surgical procedure varies somewhat between medical centers, targeting of the STN during DBS surgery generally includes preoperative stereotactic imaging (MRI), used in conjunction with stereotactic atlases. This step is followed by intraoperative electrophysiological techniques consisting of the conversion of neural activity, in the form of single-unit neuronal activity (SUA) recorded at the microelectrode tip, into audio and visual signals. This procedure is experience-based and depends critically on the neurosurgeon's and neurophysiologist's ability to recognize entry into the STN, based on a variety of cues.

In order to obtain a three-dimensional map of the STN and surrounding structures, multiple microelectrode recording tracks (typically up to five) (Benabid et al., 2009) are carried out, either sequentially or simultaneously. Determination of the optimal track for DBS implantation is a key component to successful therapeutic outcome. Optimal track selection is primarily based on microelectrode recording of single unit activity (MER/SUA), which is used to identify cells with firing characteristics consistent with STN neurons and response characteristics confirming the motor sub-territory of the STN (Falkenberg et al., 2006). Despite the common usage of MER/SUA during stereotactic surgery for PD, limitations of this technique include difficulties interpreting complex signal patterns to localize the anatomical borders of the STN, highly overlapping spiking characteristics of single neurons around the target structure, recording SUA from a very small region, sensitivity of SUA to noise, susceptibility of SUA to small amounts of blood or edema within the microelectrode track, and the binary nature of SUA (unlike local field potentials; LFP), all of which may affect the accuracy of STN localization in PD (Chen et al., 2006; Gross et al., 2006; Novak et al., 2011). The caliber of single-unit recordings can be easily diminished due to drift of the recorded unit away from the electrode tip, as a consequence of transmitted pulsations of the brain and other environmental conditions (Sanghera and Grossman, 2004).

Interpretation of SUA recordings with computational intelligence was proposed as a new approach to help clinical decision making in the operating room (Wong et al., 2009). However, such approaches are still susceptible to the challenges of isolating single neurons in the operating room. LFPs represent the aggregate activity of neuronal populations, and are particularly sensitive to synchronous and oscillatory firing patterns (Priori et al., 2004; Gross et al., 2006). Recent studies indicate that LFPs in PD correlate with both motor and nonmotor symptoms of the disease, and their signals are more robust than SUA (Priori et al., 2013; Thompson et al., 2014). Importantly, LFPs are an objective and quantitative metric while MER/SUA is more qualitative and subject to inter-practitioner variability. Although, the functional role of LFPs during DBS surgery is not fully established, we propose that they can be used to contribute to target localization in PD. In the present study, for the purpose of assisting with clinical decision making, we aimed to develop an automated approach by processing LFPs from multiple tracks to localize the dorsal border of STN and predict the macroelectrode implantation track identified by the neurosurgeon based on SUA interpretation. In the next sections first we describe our data collection methods and then detail our signal processing and classification techniques. We study the role of LFP sub-bands in prediction of the location of the dorsal border of STN. Moreover, we explore different decision criteria fused with LFP sub-band features toward prediction of optimal track selected by the neurosurgeons. We show experimental results obtained from 22 patients and discuss our results and demonstrate that LFPs can be used effectively in the operating room for clinical decision making.

## METHODS

### Patients and Surgical Procedure

This is a multicenter study in which patients were recruited in either University of Minnesota Medical Center or Baylor St. Luke's Medical Center. The experimental protocol was approved by the Institutional Review Boards of the University of Minnesota and Baylor College of Medicine. All patients provided written informed consent to participate in the study. Intraoperative LFPs were recorded from 22 patients (14 men, 8 women), who were diagnosed with idiopathic PD, and exhibited typical motor symptoms which were tremor, rigidity, and bradykinesia. Disease duration ranged between 1 and 20 years, with a mean of 10.55 years (standard deviation of 4.7 years) (**Table 1**). All patients discontinued short-acting Parkinson's medications at 12 h prior to surgery, and long-acting medications at least 24 h prior to surgery. As per standard clinical protocol, target coordinates and trajectory to the STN, were identified by preoperative stereotactic MRI, which was fused to a stereotactic computed tomography (CT), on a neuro-navigational platform (StealthStation, Medtronic Corp, MN). Then, again based on standard clinical protocol, three simultaneous tracks were performed in each subject (**Figure 1A**). The superior and inferior borders of STN, along with the optimal depth for positioning the DBS electrode, were determined by the clinical team via electrophysiological mapping using MER/SUA, and the DBS

*† UPDRS* = *Unified Parkinson's Disease Rating Scale.*

*‡ Pre-Operative Medication OFF-to-ON UPDRS Scores: Total Improvement.*

electrode was implanted by the neurosurgeon based on these spatial data, followed by macro stimulation to confirm electrode location based on benefit and side effect profile—i.e., location within motor territory of STN, but not so close to border with adjacent internal capsule or medial lemniscus, that lowthreshold stimulation- induced side effects were detected followed by confirmatory intra-operative imaging modalities. Surgeries were performed in awake patients, under the benefits of local anesthesia. In 3 of 22 patients, microelectrode mapping of right and left STN occurred on different surgical dates, as the surgical procedures were staged for clinical reasons. Therefore, these recordings were counted as separate, enabling 25 individual STN microelectrode recordings for LFP-based optimal track prediction.

#### Intra-Operative Recordings and Track Selection

Following standard stereotactic techniques, and insertion of three brain cannulas and microelectrodes (Abosch et al., 2002), MER/SUA recording was carried out using a Microguide system (AlphaOmega Inc., USA) at 12 kHz. Simultaneous LFPs were recorded with an XLTEK-EMU128FS system (Natus, San Carlos, California) at 2 kHz with 16 bit A/D resolution or gHIAmp (gTec Inc., Graz, Austria), 38.4 kHz with 24 bit A/D resolution. The LFP recordings were obtained from a 1 mm wide stainless steel contact which is 3 mm (NeuroProbe, AlphaOmega Inc., USA) or 1 mm (MicroTargetingTM, FHC Inc., USA) above the SUA recording tip and referenced to the cannula (**Figure 1B**). All microelectrodes were advanced toward the estimated target using a NeuroDrive (AlphaOmega Inc., USA) with micrometer resolution. In order to synchronize the SUA and LFP recordings, the digital depth information of the NeuroDrive is transmitted from the MicroGuide system to LFP recording system using a TCP/IP connection. Initial recordings began 20 mm above the intended final location of the electrode tip ("target") as determined by direct targeting methods and proceeded until the electrode reached 3 mm below the MER-determined target, within the substantia nigra. Electrodes were lowered in 1 mm steps until 10 mm above "target," and then in 0.5 mm steps. Duration of recordings at each depth was 15–30 s.

At each depth, the subjects sequentially rested and after a certain depth (<10 mm) executed limb movements for 10–15 s period. The neurosurgeon used standard clinical techniques for localizing the STN, via real-time auditory and visual analysis of the recorded SUA. The dorsal, ventral, and posterior borders of the STN were identified by noting increased background noise and cell firing rate, and the STN neurons were examined for movement-responsive receptive fields. In particular, the superior border of STN was determined when the background activity increased and border cells were first observed among one of the tracks in MER-SUA. This position was used as the target value in STN border identification. Among three tracks, the track with the longest span of bursting cell firing and movement responsive fields was selected for the chronic DBS electrode implantation. This track was labeled as the "optimal track" and used as the target variable in LFP based track number prediction. The neurosurgeons were blinded to the LFP recordings in the operating room and the identification of STN dorsal border and the selection of optimal track for chronic DBS electrode implantation was not influenced by the LFP recordings. The STN border and the track selected by the neurosurgeon are predicted with the LFP data using the signal processing methods detailed below.

#### Signal Processing

All data were analyzed offline. Before any processing, all recorded signals were visualized with a custom in-house developed software and annotated to distinguish artifact and/or epochs of resting, active, and passive movements. Based on these annotations, resting state data were extracted into MATLAB (Mathworks, Natick, Massachusetts) and the data recorded by XLTEK system and gHIAmp system were downsampled to 1 and 1.2 kHz, respectively, for analysis.

A schematic diagram of our signal processing pipeline is given in **Figure 2**. As an initial step the raw signals were visualized and it was observed that tracks were difficult to distinguish, due to a high amount of common activity masking spatially localized activity and/or artifacts resulting from abrupt movements of the patient and other environmental factors. In order to eliminate the common activity among tracks, but still preserve the trackspecific neural activity, the LFP data on each track were decorrelated using a least mean square (LMS) algorithm with a steepest descent update. The general formula for the decorrelation method is as follows:

Initialization: w(0) = 0.5, n = 0, 1, 2, . . .

$$\mathbf{y}^{\prime}(n) = \mathbf{w}^{T}(n)\mathbf{x}^{\prime}(n) \tag{1}$$

$$e\left(n\right) = d\left(n\right) - \left.\chi\left(n\right)\right|\tag{2}$$

$$\hat{e}\left(\eta\right) = \begin{cases} \text{sign}\left(e\left(\eta\right)\right) \* \text{ 20} & \text{if } |e\left(\eta\right)| > \text{20} \\ e\left(\eta\right) & \text{otherwise} \end{cases} \tag{3}$$

$$
\omega(n+1) = \omega(n) + \mu \hat{e}(n) \ge (n) \tag{4}
$$

where y(n) is the filter output, eˆ(n) is the residual which is the de-correlated signal, d(n) is the desired signal, µ is the step size, and w(n) is time varying filter coefficient (Hayes, 1996). In the current method, each channel, d(n), was predicted by using a

linear weighted combination of other two channels, x(n). LFP activity from three tracks were recorded continuously during the entire surgery while the microelectrodes were traveling to the estimated target. Consequently, the signal characteristic varied over depths. Since in each depth the signal was recorded for 15– 30 s, temporal variability exists in the signal. Therefore, the filter coefficients, wT (n), were updated on a sample by sample basis recursively to make the system to adapt to depth and time varying signal properties. At each iteration, the error, e(n), was calculated and this residual was used as the de-correlated LFP data in future steps for feature extraction and visualization. At 20 mm above the estimated target, all three tracks showed very similar signal characteristics indicating that they were in the white matter. Therefore, the initial filter coefficients were selected as the average of two channels with equal weights with the initialization of the filter coefficients w (n) = 0.5. By using this adaptive approach, we aimed to eliminate the common activity across tracks and suppress localized artifacts caused by patient movements and environmental factors. In order to prevent the system from being affected by high amplitude artifacts and to preserve the robustness, the error was saturated by using a 20µV threshold (Equation 3). This threshold was determined experimentally and we observed that the system recovered from localized artifacts pretty fast even if the artifact amplitude was too large.

Due to differences in spatial correlation of low and high frequency bands, the monopolar data were, first, decomposed into two frequency bands which were 8–200 and 200–450 Hz by using a 2nd order Butterworth IIR filter (**Figure 2**). The LMS algorithm was individually applied to these subbands with step size of µ = 0.0002. Each track was de-correlated by using LFPs on the other two tracks. The algorithm was applied to each depth by transferring filter coefficients to the next depth. In this way, filter coefficients were not required to start from 0.5 at each new depth so that the algorithm would adapt faster and can use both temporal and spatial information of the past. Decomposed and de-correlated data were re-merged and spectral analysis was performed. In this regard, a modified Welch periodogram method with a robust statistics was used (Telkes et al., 2014). Specifically, a fast Fourier transform (FFT) was computed with a 1024 samples long Hanning window and the window was shifted with 50% overlap. Since some artifacts destroy or dominate the power spectrum estimate obtained with mean operator, the median of the spectra of all sliding windows was calculated to eliminate localized artifacts in the spectrum. The method was repeated for each depth and all spectra were combined to visualize depth-varying power spectrum of LFPs on multi tracks. Generated depth-frequency maps (DFMs) were resampled with a 0.25 mm depth resolution and linearly interpolated to obtain equidistance depth values. The maps were smoothed using a Gaussian kernel filter to suppress noise and to reveal beta and high frequency band oscillations (HFOs). Then, DFMs were normalized with the average baseline of three tracks and transformed into log scale using the Equations (5) and (6). The tracks were not normalized by their own baseline but by the mean of all three tracks in order to compare the signal power between them. The baseline used for normalization was selected as the highest depths which assumed to be in the white matter. Therefore, the baseline was determined as top 5 depths (20– 15 mm above the estimated target) in 22 recordings. However, in rest of the three recordings, since the analysis started from lower than 20 mm (such that 18 mm) due to artifacts, the baseline segment was kept shorter and selected as top 3 depths. The purpose of using higher depths was to avoid from including any thalamic activity in normalization segment. The baseline normalization formula is noted below:

$$
\overline{b}\_{\text{avg}} = \frac{\left(\overline{b}\_1 + \overline{b}\_2 + \overline{b}\_3\right)}{3} \tag{5}
$$

$$
\overline{n}\_{\rm dfm} = \ 20 \times \log\_{10} \left( \frac{\overline{r}\_{\rm dfm}}{\overline{b}\_{\rm avg} + \Phi(f)} \right) \tag{6}
$$

where b1, b2, b3 are the baseline spectrum of each 3 track, bavg is the average baseline power, rdfm indicates the DFM, 8(f)

is a small regularization parameter which is applied for each frequency f and ndfm is the normalized DFM.

In order to observe the depth-varying frequency content of LFPs, DFMs of the patients were visualized. We noted that when the electrodes reached the STN border identified by the neurosurgeon, generally there was also an excessive activity in the beta and HFO range. In order to identify the most beneficial track along with the dorsal border of the STN, the sub-band power was extracted from all tracks and normalized by using a subject-specific average baseline. Based on the distribution of neural activity on the tracks, the sub-band frequencies were designated 11–32 Hz for beta band and 200–450 Hz for HFOs. The distribution of power in the STN among all tracks and the distribution of power only on the selected track inside and outside of STN (above the dorsal border of STN) were investigated by box and whisker plots. Student's t-test with twosample was used to check if the distributions were significantly different or not.

#### Classification

After sub-band power features were normalized between zero and one with a Max-Min normalization method for inter-subject comparison, a linear discriminant analysis (LDA) was used for classification. The principle of LDA is to maximize the separation of classes while keeping the within class densities small by using linear combination of features, −→v · −→z (Alpaydin, 2010). The linear discriminant function:

$$\begin{aligned} g\_i \left( z \middle| \nu\_i, \ \nu\_{i0} \right) &= \nu\_i^T z + \nu\_{i0} \\ &= \sum\_{j=1}^d \nu\_{ij} z\_j + \nu\_{i0} \end{aligned} \tag{7}$$

where gi (z) is the discriminant function for the input features zj with sum of the weights vj and threshold values vi0.

#### Localization of the Dorsal Border of STN

In the present study, the dorsal border of STN identified by clinical team is predicted from the depth varying LFP data by using the decision distance of a linear classifier as shown in **Figure 3A**. First an LDA classifier was trained by contrasting the LFP sub-band features coming from inside and outside of STN (above the dorsal border of STN). This classifier was evaluated at each depth and the returned decision distance was used as a measure of confidence. The depth with the highest confidence for IN-STN decisions was marked. Then we traced the decision distances above this depth and found the location where the LDA classifier voted for OUT-STN. This point where the classifier made IN vs. OUT decision transition was finally chosen as the predicted dorsal border of STN. The difference between prediction and the STN border identified by MER-SUA was calculated in each patient and the root mean square (RMS) of the prediction errors was used to quantify the performance of the classifier. Further, statistical analysis by using Student's t-test and F-test was conducted in order to compare the mean and the variance of predictions obtained by different subband features, respectively.

#### Prediction of the Optimal Track

The optimal track selection among three tracks is done by the neurosurgeons through interpreting the excessive single cell firings within the STN. Consequently, for the prediction of the optimal track using LFP data, an LDA classifier was trained by contrasting the LFP subband features of selected track vs. unselected tracks below the dorsal border of STN. This classifier was evaluated at all depths as in the STN border prediction and the returned decision distance was used as a measure of confidence. The distance returned by the linear classifier

optimal track. (A) Decision strategy in the prediction of dorsal border of STN. The classifier trained by contrasting the LFP sub-band features coming from inside and outside of STN gives a decision at each depth from top to bottom. These decision distances were used as a measure of confidence. Depth with the highest confidence for IN-STN decisions was marked. Then decisions were traced back where the LDA classifier voted for OUT-STN. This point where the classifier made IN vs. OUT decision transition was chosen as the predicted dorsal border of STN. (B) Decision strategy in the optimal track prediction. Classifier trained by contrasting the LFP subband features of selected track vs. un-selected tracks below the dorsal border of STN returned decision distances at all depths. The decision distances were used in three different scenarios: (1) Below the SUA-based STN border, (2) below the one standard deviation from the average of STN border, and (3) below the LFP-based STN border obtained from prediction of dorsal border of STN. The decision distance with the longest span was used as a decision criterion to predict the optimal track.

was used in three different scenarios for final decision making (See **Figure 3B**). In the first scenario, the optimal track was predicted below the STN border provided by the neurosurgeon for that specific test subject based on the SUA interpretation. This represents the setup in which we fuse SUA- and LFPbased information. In the second scenario, no SUA information about the STN border of the test subject was used, and the decisions were given below the one standard deviation from the average of STN border estimated from all training subjects. In the third scenario, the optimal track decisions were made below the STN border which was derived solely from the LFP data. Specifically, here we explored whether or not the LFP could predict the optimal track without any SUA-based interpretation. We studied the classification performance below and above the STN border in these three different scenarios. A schematic diagram related to this process is given in **Figure 3B**. Depthvarying LFP subband features of each track were classified using the trained LDA and a label and related decision distance were generated by the classifier for each depth. We classified one of the tracks as the optimal one based on the longest span of decision distances voting for optimal track within the STN. Note that the longest span is a common approach used intraoperatively by neurosurgeons for MER-SUA-based optimal track selection. Also note that, if the track selected by neurosurgeon in the operating room did not match with the decision of the algorithm, the decision was counted as a misclassification.

The prediction of optimal trajectory was investigated using individual sub-band powers, beta and HFO, and their combination. To explore the benefit of the LMS algorithm over monopolar signals (raw signals), the same classification procedure was carried out with the raw data.

Finally, in order to assess the efficiency and reproducibility of the classification, a leave-one-subject-out method was used. In each step, one subject was used for testing, whereas the other subjects were used for training the LDA classifiers for STN dorsal border and optimal track prediction. The procedure was repeated until the whole sample was classified. In addition, this procedure was performed separately for individual beta and HFO subbands of LFP and their combinations to examine their efficacy in classification performance.

In order to explore a relationship between classification results and post-operative simulation parameters used for the initial programming 6 months after the surgery were compared in correctly classified and misclassified groups. The distribution of stimulation amplitude was investigated by box and whisker plots. Student's t-test with two-sample was used to check if the distribution of simulation amplitude, frequency, and pulse width were significantly different or not.

## RESULTS

### De-Correlation of LFP Data from Multiple Tracks

We analyzed LFP data derived from 75 MER tracks in patients with PD who were undergoing STN DBS electrode placement. Typical raw SUA and LFP data coming from various depths were shown in **Figure 4**. The red dashed line indicates the dorsal border of STN. The correlation between the two modalities is clearly seen. When the electrode enters the STN, both single cell firing and oscillatory activity increase.

In **Figures 5A,B**, the mean correlation matrices of raw LFP data filtered in beta band (11–32 Hz) and HFO band (200– 450 Hz) were shown. The correlation matrix in the beta band (**Figure 5A**) explicitly shows that the spatial correlation between tracks is high whereas the correlation between tracks in HFO range is small (**Figure 5B**). The small amount of correlation at HFO band in raw data also shows that oscillations at higher frequencies are more localized than the oscillations at lower frequencies. For these reasons, the LFP data were de-correlated with LMS separately in these frequency bands. It was found that the correlation between tracks is reduced after the LMS-based preprocessing step (**Figure 5C**) which helped more to distinguish the tracks.

**Figure 6** demonstrates the effect of LMS algorithm by comparing it to the raw LFP signals. As it can be seen in **Figures 6A–C**, the raw LFP data have a high amount of common activity across all tracks at various depths which masks the spatially and temporally distinguishing patterns during targeting. In **Figure 6B**, DFMs indicate that the high-energy low band activity among tracks masks other oscillations. The common activity across three tracks and the high energy low band oscillations can be also seen in the power spectrum shown in **Figure 6C** which was generated from the LFP data below the dorsal border of STN. On the other hand, target specific oscillations are clearly seen on de-correlated LFP data (**Figure 6D**). In particular, the energy in the first track is much higher than the other two tracks and it is easier to observe the track differences and the estimated STN border depth for the target localization. The DFMs of these tracks shown in **Figure 6E** demonstrate that the first track contains LFPs with higher energy in low and high frequency bands below the dorsal border of the STN which is marked with a white dashed line. Furthermore, the power spectrum shown in **Figure 6F** demonstrates that not all three tracks show excessive beta activity. There is an increase in the gamma band (35–55 Hz) and great enhancement in HFO range (200–400 Hz) in the first track compared to other tracks. The LMS algorithm not only reveals the pathological beta oscillations but also the HFOs having lower energy.

#### Spatio-Spectral Patterns of Multitrack LFP

In order to provide a sense of the depth-varying frequency content of multitrack LFPs, we demonstrated representative normalized DFMs of de-correlated LFP data of all three tracks from four patients in **Figure 7**. In each map, the dorsal STN border is marked with a white dashed line. The excessive beta oscillations can be clearly seen in the first subject dominantly in the center track and localized to certain depths (**Figure 7A**). The power of beta oscillations in the posterior and medial track is weak, yet it can be still observed. Furthermore, there is a strong and track-specific HFO around 350 Hz which is well aligned with the low band activity. On the other hand, in the subject presented in **Figure 7B**, beta oscillations are observed in all tracks along with HFOs. Although the excessive LFP activity occurs below the STN dorsal border as for the patient presented in the **Figure 7A**, the excessive depth varying spectral patterns are pretty track and region specific. The HFO in the center track sits at 350 Hz while it is located at 250 Hz in the posterior track. The lateral track shows wider but weaker oscillations. DFMs in **Figure 7C** demonstrate similar LFP characteristics to the first subject (**Figure 7A**) with dominant beta oscillations and HFOs in the center track. Similarly, these oscillations are well aligned below the dorsal border and highly stronger than the beta oscillations in the other two tracks. Distinctly, strong oscillations at higher depths are observed above the dorsal border in the posterior track which might be related to thalamic activity. Note that we observed high frequency activity localized at higher depths above the dorsal border of STN in at least one of the un-selected tracks in 56% of recordings. The number of the un-selected tracks with the observed oscillations were: 8 posterior, 5 center, and 1 medial. Similar to the HFO activity seen in **Figure 7C**, these oscillations were noted from 11.5 ±

specific oscillations are visible on de-correlated LFP data. The red dashed line shows the dorsal border of STN determined by MER-SUA. (E) DFM of de-correlated LFPs in each microelectrode track. The high energy oscillations in distinct low and high frequency bands inside the STN (below the white dashed line) can be clearly seen in the first track. (F) Power spectrum of de-correlated LFPs. Tracks do not show excessive and correlated beta activity anymore. There is an increase in the gamma band (35–55 Hz) and great enhancement in HFO range (200–400 Hz) in the first track compared to other tracks.

2.6 to 5.7 ± 2.4 mm (average values) above the estimated final location of the electrode tip. The tracks having higher-depth HFOs do not include strong beta activity. These oscillations have a longer spatial span with lower power. It is likely that these oscillations rise from thalamic structures (Hutchison et al., 1998; Falkenberg et al., 2006), and given their spatial distribution in relation to the beta band activity, they might be used as markers for STN localization. The fourth representative subject shown in **Figure 7D** introduces a different LFP characteristic compared to others. None of the tracks are associated with strong, long span of beta oscillations. Specified border is not aligned with weak beta oscillations but the short lasting excessive one in the lateral track. All tracks demonstrate spatially different weak-to-minor HFOs. The overlap in the LFP activity between tracks, the weak activity across tracks and thalamic oscillations are some of the factors contributing to the challenges to the prediction of dorsal border of STN and selection of optimal track.

#### Intra-Track and Inter-Track Differences of LFP Spectra

For the neurosurgeon selected track, the distribution of beta and HFO subband powers above and below the dorsal border of STN are given in **Figure 8A**. The analysis shows that there is a significant difference between the power inside and outside the STN region (above the STN dorsal border) (t = 44.72, p < 0.001; t = 34.89, p < 0.001) in the selected track. As seen from the box-plot in **Figure 8A**, the sub-band power is much higher

inside the STN. When the subband power was compared between the optimal and other tracks (**Figure 8B**), the distributions were found significant as well (t = 16.47, p < 0.001; t = 15.17, p < 0.001). The significance is consistent at the beta band and HFO band in both distributions. The variance of HFO power in the un-selected tracks is higher than the variance in the selected track. Based on our previously mentioned findings, we postulate that thalamic activity in un-selected tracks might contribute to increased variance of HFO power when the distribution includes entire track.

#### Localization of the Dorsal Border of STN

The progression of prediction of STN dorsal border for representative subjects and the average results estimated from the entire patient population is shown in **Figure 9**. In **Figures 9A,B**, decision distances returned by the classifier voting either for IN-STN or OUT-STN are shown for two representative subjects. The decision distance curves were obtained from the fused beta and HFO features. Note that the predicted STN border is shown with an arrow corresponds to the position where we find the maximum confidence point associated with IN-STN and trace back to the depth crossing zero. The dorsal border of STN provided by the neurosurgeon based on SUA interpretation is shown with a dashed vertical line. **Figure 9A** shows a late prediction of the dorsal border (e = −0.75 mm) while **Figure 9B** indicates an early border prediction (e = +1 mm). **Figure 9C** demonstrates the average border decisions with the standard deviation coming from all test subjects by using individual subband powers, beta and HFO, and fused features. The overall localization error of the dorsal border of STN was quantified by calculating RMS of the error between target values and LFP predictions across all subjects. The red and blue lines show the decisions obtained with beta and HFO band features indicating an RMS error of 1.98 mm and 1.18 mm, respectively. The mean value of prediction error for beta band features was 0.83 ± 1.84 mm while the mean of error for HFO band features was −0.23 ± 1.18 mm. The decisions obtained through the fused beta and HFO band features had an RMS error of 1.22 mm with mean of 0.24 ± 1.22 mm. In **Figure 9D**, the distribution of prediction errors are shown for each studied subband and their fusion. Student's two sample t-test analysis indicated that the difference between mean values of beta-based prediction error and HFO-based prediction error was significantly different (t = 2.22, p = 0.0322) while no statistically significant difference

FIGURE 9 | Localization of dorsal border of STN. Panels (A,B) show the decision distances returned by the classifier voting either for IN-STN or OUT-STN for two representative subjects. Predicted STN border is shown with an arrow corresponds to the position where we find the maximum confidence point voted for IN-STN and trace back to the depth crossing zero. The black dashed lines indicate the dorsal border of STN obtained by MER/SUA. The predicted border is 0.75 mm lower than the SUA-based border in the first subject (A) and 1 mm higher in the second (B). Panel (C) demonstrates the average border decisions with the standard deviation coming from all test subjects by using individual sub-band powers, beta (red) and HFO (blue), and fused features (green). (D) Comparison of prediction errors. The error (RMS) of the predictions obtained with beta and HFO band features are 1.98 and 1.18 mm. The RMS of the error in predictions obtained from fused features is 1.22 mm. The variance of RMS values were compared by *F-test*. The difference between beta-based and HFO-based prediction error is marginally significant (*p* = 0.054) with α = 0.05. +indicates the outlier error value.

was found neither between beta-based prediction error and the error of fused features (t = 1.23, p = 0.2244) nor HFObased prediction error and the error of fused features (t = −1.25, p = 0.2185). When the variances of these distributions were compared by using an F-test, the analysis showed that the difference between beta-based and HFO-based border prediction was only marginally significant [F(1,2) = 2.42, p = 0.054] while there was no statistically significant difference between the variance of individual sub-band powers (beta and HFO) and fused power [F(1,3) = 2.26, p = 0.075; F(2,3) = 0.94, p = 0.88, respectively].

#### Prediction of the Optimal Track

We studied the optimal track classification in three different scenarios using the LFP data: (i) below the SUA-based STN border, (ii) below one standard deviation from the average STN border obtained from training data, and (iii) the LFP-based STN border. We trained the LDA classifier using individual subband powers and their combination. Our results toward the prediction of optimal track from LFP data is given in **Table 2**. We note that the best results were obtained from the combined subband power features and consistently in all these scenarios, the optimal track prediction accuracy was 80% (shown in bold type) indicating that the classifier can predict the track targeted to the STN in 20/25 recordings. These results show that prediction of optimal track can be performed independently from single unit recordings. When the beta and HFO subband features were used individually, the classification accuracy dropped to 72 and 68% respectively. When the procedure was repeated with the raw data, the prediction rate was poor. In particular, the classification accuracy was 64% for beta band power and 68% for the HFO and fused features which supports the observation that HFOs obtained in monopolar configuration are already highly de-correlated among different tracks.

Despite the spatially localized thalamic oscillations, the classification results obtained above the STN border were quite poor. The prediction accuracy was found to be 40% when the classification was computed above SUA-based or LFP-based STN border. Decision accuracy with average STN border was even lower at 36% by using fused sub-band power. The results indicate that the LDA classifier trained with the LFP features above the STN does cannot predict the optimal track with a reasonable accuracy and was close to chance level.

The progression of classification over depths for three representative subjects are shown in **Figures 10A–C**. In each plot, the STN border location provided by the neurosurgeon based on the SUA interpretation is also represented with a vertical dashed line. The decision distances in both selected and un-selected tracks returned by the classifier are close to each other down to the dorsal border of STN. Since the spectral characteristics of LFPs change inside the STN compared to higher depths, we observe a sudden change between the decision distances as well. If only one of the tracks deviates from the others and reaching the highest confidence level, it is easily determined as the optimal track by the classifier. If more than one track are voted for the optimal track (below the zero line in **Figure 10**) with high confidence levels, algorithm gives the optimal track decision by computing the longest span of the selected track votes. The progression of the classification for a misclassified subject is given in **Figure 10C**. The average decisions for the optimal track of all subjects with the associated standard deviation are given in **Figure 10D**. A clear separation is observed in decision distances between selected and un-selected tracks indicating a high percentage of correct prediction among the subjects.

#### Distribution of Selected Tracks

**Table 3** shows the frequency of selected tracks based on MER-SUA interpretation and LFP processing. As per standard clinical protocol, the initial trajectory to target STN is determined by preoperative stereotactic imaging. Three tracks are selected by the neurosurgeon based on the initial planning for microelectrode recordings. The initial expectation is that the center track will hit the STN while other tracks account for possible targeting

#### TABLE 3 | Frequency of track selection in the present study.

#### Comparison of SUA-based and LFP-based track selection

*The selection frequency based on LFPs are shown in bold type.*

error. Then based on the MER-SUA recordings the optimal track is selected among these three trajectories. Although the image based planning aims to hit the STN with the center track, **Table 3** demonstrates that intraoperative MER-SUA-based decisions among 25 recordings is not biased toward the center track. We note that the selection frequency is higher in anterior track based on MER-SUA mapping. In addition, the posterior track is not selected at all. When the selection frequency based on LFPs is studied, it can be seen that both MER-SUA and LFP decisions (shown in bold) match with a high percentage but LFP based prediction was more in favor of the center track. Overall, our results indicate that stereotactic planning does not perfectly correlate with intraoperative electrophysiology based track selection and highlight the variance in track selection.

#### Post-Operative Programming Parameters

We explored whether there exist any difference in programming parameters between the correctly classified and misclassified patients. In particular we investigated the post-operative simulation parameters such as voltage, frequency and pulse width which were selected during the programming 6 months after the surgery. The distribution of selected stimulation voltages are presented in **Figure 11**. We note that the average stimulation voltage used in correctly classified group is 1.72 ± 0.63 V while it is 2.12 ± 0.69 V in misclassified group. Student's two sample t-test analysis indicates that the difference in voltages between two groups is not statistically significant (t = −1.16, p = 0.2595). However, one of the misclassified subjects we observed has high level beta and HFO activity in both selected and un-selected track and this patient is stimulated with 1 V. This indicates that both tracks could be viable. When this subject is removed from the misclassified group, we note that the mean post-operative stimulation voltage increases to 2.4 V ± −0.46 for the misclassified population. The difference in post-operative stimulation voltages between correctly classified and misclassified groups without this outlier becomes marginally significant (t = −1.92 p = 0.0685). No significant difference is found either in the frequency (183.1 ± 5.8 Hz) or in the pulse width (90 ±

the un-selected track-2 is passing zero line and reaching the highest confidence level, it is determined as the optimal track by the algorithm. Panel (D) shows the average optimal track decisions with the standard deviation coming from all test subjects. It shows a clear separation in decision distances between selected (red) and un-selected tracks (blue) indicating a high percentage of correct prediction among the subjects.

26.5µs) between groups (t = −0.74, p = 0.4692 and t = 0.96, p = 0.3477, respectively).

## DISCUSSION

Significant variability exists in the axial and coronal orientation of the STN in humans (Patel et al., 2008), and the motor territory of the STN is small, measuring ∼4–6 mm extent from dorsal to ventral. These factors combined with brain shift between preoperative stereotactic imaging and intraoperative electrode brain penetration can lead to targeting errors in the operating room. Pre- and intra-operative clinical imaging methods alone are suboptimal for accurate placement of a DBS electrode; they are subject to distortion, and visualization of a clear differentiation between the STN and surrounding structures can be difficult. In this regard, our study also indicates a considerable amount of variance in track selection. Under the assumption of hitting the STN through center track by image based planning, track selection was not found to be biased toward the center track in the intraoperative MER/SUAbased decisions indicating that stereotactic planning does not perfectly correlate with intraoperative electrophysiology based track selection.

Accurate localization of STN motor territory through intraoperative electrophysiology is a crucial step for DBS electrode implantation (Zonenshayn et al., 2000; Gross et al., 2006). As recently as 2013, an international survey of highvolume DBS implanting sites revealed that 83% of centers use microelectrode recording indicating that the most commonly used electrophysiological mapping method remains MER-SUA recordings (Abosch et al., 2013). However, the method has several limitations in practice as subjective interpretation of complex signal patterns to localize the anatomical borders of the STN, being less stable and more "susceptible to technical (e.g., impedance) and physiological (e.g., cerebrospinal fluid and blood) fluctuations" (Thompson et al., 2014). As Gross et al. indicated, the number of groups using solely macroelectrode/DBS mapping to target the STN without any microelectrode recording is high (Chen et al., 2006; Gross et al., 2006; Telkes et al., 2014). Although there are advantages of using macroelectrode/DBS technique alone in STN targeting like carrying less amount of risk for intracranial hemorrhage since there is no multiple trajectories and due to the blunt-tip of the macroelectrode (Xiaowu et al., 2010), the drawbacks like microlesion effect which might limit the clinician's ability to test or therapeutic effectiveness in the operating room and poor spatial resolution should be taken into consideration (Rezai et al., 2006; Wang et al., 2014). Although asleep, MRI-based non-MER-guided surgery is gaining sway (Starr et al., 2014), the possibility of clinically testing a DBS electrode prior to permanent

implantation does not exist as yet in the context of such a procedure.

Earlier investigations have documented that excessive beta oscillations in certain basal ganglia structures, especially the STN, represent a pathophysiological feature of PD (Weinberger et al., 2006; Kane et al., 2009; Lopez-Azcarate et al., 2010; Oswal et al., 2013). Excessive beta band (8–30 Hz) activity is detected when the electrodes enter into the STN (Levy et al., 2002; Kühn et al., 2008; Brittain and Brown, 2014). Similarly, excessive oscillations at very high frequency ranging from 200–400 Hz are also observed (Priori et al., 2004; Lopez-Azcarate et al., 2010; Özkurt et al., 2011). Even though these high frequency oscillations are considered to represent a pro-kinetic state, and appear with dopaminergic medication and/or induced movement (Foffani et al., 2003, 2006; Trottenberg et al., 2006), others have demonstrated that HFOs (>200 Hz) can still be observed in the STN during the medication OFF state or at rest (Lopez-Azcarate et al., 2010). In our study, all patients discontinued with their short and long acting medication before the surgery and were in OFF state. As others, we observed HFOs in the resting state and increased band power along with entry into the STN. Existence of excessive beta band and HFOs within the STN in PD can be used in target localization. However, the variability and patientspecific characteristics of spatial distribution of excessive beta and HFOs should be taken into consideration (Chen et al., 2006; Weinberger et al., 2006; Wang et al., 2014).

Despite a few publications using intraoperative microelectrode LFPs for STN localization (Michmizos et al., 2008; Holdefer et al., 2010; Wang et al., 2014), to our knowledge, no studies exist on the functional use of LFPs recorded from multiple microelectrodes for the selection of the optimal trajectory targeting the STN in PD. The present study demonstrates that using spectral features of LFP to identify the optimal track without any decorrelation technique provide suboptimal results due to widely distributed neural signals and/or artifacts masking the spatially and temporally distinguishing patterns during targeting. Therefore, the LMS algorithm is used as an efficient technique to decorrelate the tracks by keeping localized activities in each. The adaptive LMS algorithm is widely used in the biosignal processing field since early 80s for signal enhancement due to its efficiency and low complexity (Widrow et al., 1975; Ferrara and Widrow, 1982; Chen et al., 1990). Since the decorrelation is being done recursively without violating the causality constraint, where each channel is predicted by the current samples of other data channels, the algorithm can be easily executed on standard PC architectures in real-time. Since it is an adaptive technique, the time and depth varying parameters allows tracking time and depth varying LFP activity and does not suffer from the cross talk as much as in the common average based derivation. It should be noted that since it estimates current signal by using a linear combination of other signals, the LMS algorithm cannot fully eliminate the high amplitude artifacts if they are not distributed among the tracks, which constitutes a major drawback of the algorithm. One way of reducing the effect of large artifacts and keep the system stable is to use an error threshold with upper and lower boundaries. Another important factor influencing the benefit of the algorithm is the learning parameter µ. It should be investigated by considering signal properties such that adaptation of the system should be neither very slow nor very fast.

Spectral analysis showed that beta oscillations are getting stronger as the electrodes approach the STN. Not only beta oscillations but also strong HFOs can be observed in the STN area well aligned with beta oscillations. This strong relation is noted in the tracks selected by neurosurgeons in 17 recordings out of 25. In rest of the eight recordings, HFO was either weak or fully absent or they were noted only in one of the un-selected tracks. The energy changes above and in the STN were used to localize the dorsal border. The RMS error of prediction for the dorsal border of STN is obtained from 1.18 to 1.98 mm when the different features are used. The minimum prediction error is found with the power of HFOs (1.18 mm) indicating that, despite the unknown functional role of these high frequency components, they can still be eligible in STN targeting. The LFP is a continuous process and does not suffer from SUA isolation challenge while the target variable is a SUA driven information which is also prone to interpretation error and isolation challenge. Considering the dorsa-ventral size of STN, 1.22 mm prediction error in depth with the fused features may represent 11% in DBS electrode 3387 or 16% in 3389 difference which can be easily compensated with the multiple contacts of the DBS electrode.

The features computed above the STN border provided poor results in prediction of optimal track. We note that the optimal track can be predicted with higher accuracy with the features obtained below the dorsal border of STN. Analysis manifest that 14 recordings out of 25 (56% of entire dataset) indicate spatially distinct HFOs together with beta activity above the dorsal border of STN (see **Figures 7C,D**) in at least one of the un-selected tracks. It can be assumed that these relatively weak oscillations located away from the dorsal border of STN are recorded from thalamic structures. To our knowledge, this considerable amount of thalamic oscillations in PD are not well studied phenomena. These findings presented here might be used as spatial markers in STN localization and might form the basis of further investigation into PD pathophysiology from a spatio-spectral perspective.

It should be noted that our classification technique could not predict the selected track in 20% of the subjects. We did not observe a gender difference between these five misclassified subjects. Specifically, three of them were men and two of them were women. None of the misclassified patients were tremor dominant. Three of these patients were typical PD and the other two were bradykinetic/rigid. The mean age and disease duration were 62.2 ± 13.6 and 12 ± 4.6, respectively and were not significantly different from other correctly classified patients (62.1 ± 8.3 and 10.1 ± 4.6, respectively). The misclassification in the 20% of the patients occurred due to many different factors in LFP signal including weak activity or similar activity between tracks. During these recordings, we did not use any sedation. Therefore, weak activity cannot be related to anesthesia. In one patient with typical PD phenotype, the LFP signal was weak across all tracks. We noted that the beta and HFO activity started to develop toward the bottom border of the planned target deeper than the other patients. We believe that in this particular patient the weak activity across all tracks can be described with the electrode positions. Our observations indicate that the three tracks just started to enter the STN and did not fully went through it. In another misclassified case, the LFP activity was quite strong and similar in two out of three tracks. Therefore, the classifier output was very close for these two tracks. In the other three patients, the HFO activity in the SUA selected track was weak compared to LFP selected track. Studies hypothesize that maximum beta band (13–32 Hz) and gamma band (48–220 Hz) power is highly correlated with stimulation programming parameters in DBS chronic electrode (Ince et al., 2010). When a particular contact pair on the electrode shows strong beta and gamma oscillations, it's assumed that the electrode is closer to the source so that lower stimulation would provide better symptom improvement and less side effects. The present study supports these results. We observed higher stimulation voltages in those patients where the LFPs did not correlate with MER/SUA selected tracks. Despite the lack of statistical significance, the stimulation voltages in the 6-mo-programming of implantable pulse generator (IPG) indicate lower values in the patient group having stronger LFPs in beta and/or HFO bands. A study with larger sample size would be needed to test the validity of this observation.

#### REFERENCES

Abosch, A., Hutchison, W. D., Saint-Cyr, J. A., Dostrovsky, J. O., and Lozano, A. M. (2002). Movement-related neurons of the subthalamic nucleus in patients with Parkinson disease. J. Neurosurg. 97, 1167–1172. doi: 10.3171/jns.2<IP_ADDRESS>7

#### CONCLUSION

The present study describes an automated approach for electrophysiological localization of STN, using microelectroderecorded LFPs acquired during DBS surgery simultaneous to MER. This work is novel, in that it is the first study to combine different sub-band features derived from beta (11– 32 Hz) and HFOs (150–450 Hz) of LFPs in order to (1) estimate the optimal track for DBS implantation, and (2) identify the dorsal STN border, with high accuracy. This work also contributes to knowledge about the neurophysiology of PD by describing the spatial localization of HFOs. Because recording LFPs simultaneous with MER/SUA does not prolong the total duration of surgery, using this technique online in the operating room would increase the chance of optimal placement of the DBS macroelectrode within the motor subterritory of the STN, without an appreciable downside. Fused with existing mapping techniques, automated online LFP analysis may increase the accuracy of the DBS macroelectrode placement. This might contribute to the efficacy of DBS by reducing the stimulation voltage and associated side effects. Since the electrode placement is guided by LFP activity, the current technique could also be useful to monitor the LFP events which are capable to fine tune the future DBS settings in a closed loop paradigm (Ince et al., 2010; Rouse et al., 2011; Priori et al., 2013).

Further prospective investigations regarding the clinical outcomes using this technique of optimal track selection are warranted.

### AUTHOR CONTRIBUTIONS

NI and IT designed the study, collected the data, wrote the manuscript, and conducted the analyses. AA and AV performed the surgeries, contributed to data collection, and interpretation of the results. JJ performed the behavioral tests during surgery, evaluated the condition of subjects, and contributed to interpretation of the results. All authors reviewed the manuscript and approved the final version of the manuscript.

#### ACKNOWLEDGMENTS

We wish to thank all the patients who participated in this study. We are grateful to Akshay Gupte for his help with clinical data collection. We are also thankful to Ibrahim Onaran, Anish Sen, and Jiang Tianxiao for their help with data collection and visualization. This research was supported by National Science Foundation, award CBET-1343548.

Abosch, A., Timmermann, L., Bartley, S., Rietkerk, H. G., Whiting, D., Connolly, P. J., et al. (2013). An international survey of deep brain stimulation procedural steps. Stereotact. Funct. Neurosurg. 91, 1–11. doi: 10.1159/000343207

Alpaydin, E. (2010). Introduction to Machine Learning, 2nd Edn. Cambridge, MA: The MIT Press.

mechanism in Parkinson's disease. J. Neurosci. 30, 6667–6677. doi: 10.1523/JNEUROSCI.5459-09.2010

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Telkes, Jimenez-Shahed, Viswanathan, Abosch and Ince. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# EEG-Based Quantification of Cortical Current Density and Dynamic Causal Connectivity Generalized across Subjects Performing BCI-Monitored Cognitive Tasks

Hristos Courellis 1, 2 \*, Tim Mullen1 , Howard Poizner 3 , Gert Cauwenberghs 2, 3 and John R. Iversen1

1 Swartz Center for Computational Neuroscience, University of California, San Diego, San Diego, CA, United States, 2 Department of Bioengineering, University of California, San Diego, San Diego, CA, United States, 3 Institute for Neural Computation, University of California, San Diego, San Diego, CA, United States

#### Edited by:

Timothée Levi, University of Bordeaux 1, France

#### Reviewed by:

David Thomas Bundy, Kansas University of Medical Center Research Institute, United States Alireza Mousavi, Brunel University London, United Kingdom

> \*Correspondence: Hristos Courellis<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Neuroprosthetics, a section of the journal Frontiers in Neuroscience

Received: 10 January 2017 Accepted: 20 March 2017 Published: 17 May 2017

#### Citation:

Courellis H, Mullen T, Poizner H, Cauwenberghs G and Iversen JR (2017) EEG-Based Quantification of Cortical Current Density and Dynamic Causal Connectivity Generalized across Subjects Performing BCI-Monitored Cognitive Tasks. Front. Neurosci. 11:180. doi: 10.3389/fnins.2017.00180 Quantification of dynamic causal interactions among brain regions constitutes an important component of conducting research and developing applications in experimental and translational neuroscience. Furthermore, cortical networks with dynamic causal connectivity in brain-computer interface (BCI) applications offer a more comprehensive view of brain states implicated in behavior than do individual brain regions. However, models of cortical network dynamics are difficult to generalize across subjects because current electroencephalography (EEG) signal analysis techniques are limited in their ability to reliably localize sources across subjects. We propose an algorithmic and computational framework for identifying cortical networks across subjects in which dynamic causal connectivity is modeled among user-selected cortical regions of interest (ROIs). We demonstrate the strength of the proposed framework using a "reach/saccade to spatial target" cognitive task performed by 10 right-handed individuals. Modeling of causal cortical interactions was accomplished through measurement of cortical activity using (EEG), application of independent component clustering to identify cortical ROIs as network nodes, estimation of cortical current density using cortically constrained low resolution electromagnetic brain tomography (cLORETA), multivariate autoregressive (MVAR) modeling of representative cortical activity signals from each ROI, and quantification of the dynamic causal interaction among the identified ROIs using the Short-time direct Directed Transfer function (SdDTF). The resulting cortical network and the computed causal dynamics among its nodes exhibited physiologically plausible behavior, consistent with past results reported in the literature. This physiological plausibility of the results strengthens the framework's applicability in reliably capturing complex brain functionality, which is required by applications, such as diagnostics and BCI.

Keywords: brain-computer interface (BCI), electroencephalography (EEG), causality analysis, source localization, LORETA, independent component analysis, motor activity, spatial reach and saccade

## INTRODUCTION

Dynamic intra-cortical interactions offer a wealth of information that can be used to develop an understanding of brain functionality and quantify the effects of brain pathologies on brain functionality. Traditional static functional characteristics of individual brain regions, and dynamic interactions within a network of brain regions (Moeller et al., 2015) offer limited insight into brain function as they typically do not examine the causality among these interactions. The advent of Grangercausal analysis provided a framework to quantify the asymmetric causal interactions between regions of the cortex through specific frequency bands, which hold functional significance for those regions (Barnett and Seth, 2014; Seth et al., 2015), and has shed light on a large number of network-related disorders, such as Alzheimer's and Parkinson's Disease (Pievani et al., 2014). However, this connectivity information is not limited in use to diagnostics, as an increasing number of brain-computer interface (BCI) applications are utilizing connectivity features for intent discrimination when conducting BCI-related tasks (Zhang et al., 2012; Kabbara et al., 2016).

While a potential wealth of information is available in intra-cortex interactions, a systematic algorithmic approach that allows for reliable generalization of cortical network dynamics across test subjects in a subject group is not readily available. Many sophisticated tools exist (Linder et al., 2011; Niso et al., 2013; Barnett and Seth, 2014), which allow for both parametric and non-parametric estimation of dynamic causal interactions among neural activity signals within an individual. However, these tools exhibit various limitations, such as strict channel-space applicability or lack of accounting of common sources, thus, giving rise to spurious causal interactions within the chosen cortical network (Trongnetrpunya et al., 2015). Additional limitations exist in the source-localization algorithms applied to reconstruct cortical-source activity from the measured channel space data. Many popular algorithms stemming from Independent Component Analysis (ICA) are not able to achieve reproducible source localization outcomes across test subjects (Bell and Sejnowski, 1995). The inability to control where these algorithms localize the sources for decomposition on separate subjects presents a large problem since the sources represent nodes in a cortical network, and an ICA-based algorithm could localize a different network for every subject in a group, making group-level network comparisons difficult. Independent Component (IC) ordering is a highly non-trivial task (Hyvärinen and Oja, 2000) for which very few automated approaches exist barring manual selection, a prohibitively tedious and time consuming task for increasingly large datasets.

In this paper, we develop an algorithmic and computational framework that addresses the aforementioned issues through a multi-step, open-source pipeline of analytic tools (Delorme et al., 2011; Iversen et al., 2014). The pipeline combines the deconvolution power of ICA; the systematic modeling of the geometry of gray matter, thereby creating a common source space across subjects with the use of individually warped boundary element models (BEMs); the selectivity afforded by Cortical Low Resolution Electromagnetic Tomography (cLORETA)—a source-space based, distributed source localization algorithm; and the robustness of the Shorttime direct Directed Transfer Function (SdDTF)—a causality measure, computed from Multivariate Autoregressive (MVAR) model coefficients, that addresses the spurious connectivity affecting Granger Causal methods (Korzeniewska et al., 2003).

A similar approach to source localization and connectivity analysis was proposed by Mullen et al. (2015) as a means of demonstrating that the connectivity features considered in this framework can be applied with great success in real-time brain state identification. In order to demonstrate the power afforded by this framework in group-level analysis, a BCIrelated task originally introduced by Park et al. (2014) was employed, requiring subjects to reach and saccade to presented visual targets while their cortical activity was recorded using electroencephalography (EEG) as the neuroimaging modality. The approach utilized for BCI-related task is a passive-BCI approach, also described as an implicit human-computer interaction, wherein passive EEG equipment that senses the user is used to evaluate cognitive states of the user (Schmidt, 2000; George and Lecuyer, 2010). While, in this case, the closedloop feedback elements of the task simply involve auditory and visual confirmation of reach/saccade task completion, the applicability of this approach to more active real-time closedloop feedback BCI systems is also discussed. Although the motor activity present in the human cortex during simple planned motor actions has been carefully documented with respect to the activity of individual cortical regions (Johnson et al., 1996; Filimon, 2010; Heider et al., 2010; Glover et al., 2012), the dynamic interactions among these cortical regions during planning and execution of the task have not been documented as extensively. Park et al. were able to differentiate between spatial movement of the arm and eye by employing channel-space empirical mode decomposition, but did not characterize connectivity between cortical regions underlying these processes. The algorithmic approach implemented in this paper was first described by Iversen et al. (2014), where the method was applied to a smaller subject group. In this paper, we introduce modifications to the original pipeline in Iversen et al. (2014) including a new maximum-power heuristic that captured previously undetected connectivity involving the Supplementary Motor Area and Precuneus, regions expected to exhibit connectivity given the motor-related BCI task, and non-parametric statistical thresholding of the computed causal connectivity values. Through application of our new pipeline, described fully in the Methods section that follows, we demonstrate in the Results section the strength and specificity of the inference that can be made regarding the dynamic causal interactions present among regions of the cortex during reach/saccade planning and execution. The Discussion section illustrates the physiological plausibility of the findings reported in the Results section with reference to previous findings reported in the literature. The Conclusion section summarizes the features, the scalability, extensibility, and the applicability of the proposed algorithmic and computational framework.

#### METHODS

#### Overview

Electrophysiological data were collected from human participants conducting a BCI-monitored task, and processed according to the methodology and computational tools outlined in **Figure 1**, using MATLAB as the computational platform. The workflow in **Figure 1** was employed to extract cortical source information from a cortical network defined for all of the participants, and, subsequently, identify statistically significant spectro-temporal interactions among the nodes (sources) of this network.

### Participants and Task

Ten healthy, right-handed individuals (age: 20.8 ± 2.6 years) participated in this study. The study was approved by the

FIGURE 1 | General outline for cross-subject quantification of cortical connectivity. Initial data processing and artifact rejection is conducted in EEGLAB (Delorme and Makeig, 2004) (A,B). ROIs are selected through IC clustering, also conducted in EEGLAB (C). Head model construction and Electric Lead Field calculation are conducted in MoBILAB (Ojeda et al., 2014) (D). Source localization is conducted in BCILAB (Delorme et al., 2011) through an implementation of cLORETA that iteratively estimates regularization hyperparameters (E). Connectivity estimation is conducted using SIFT, an extension of EEGLAB (Delorme et al., 2011) (F,G). All of the above software and their internal dependencies are run in MATLAB.

Human Subjects Institutional Review Board of the University of California, San Diego. Written informed consent was obtained from all subjects prior to the start of the experiment.

The task involved a time-constrained eye movement and reach using the subject's dominant hand (right) from a target in the center of a touch-screen to a laterally-offset target as depicted in **Figure 2** (Park et al., 2014). Each trial was either a control condition (lift), where the participant moved the stylus off of the central fixation dock after a brief target onset period (500–700 ms), or a test condition (reach/saccade), where the participant shifted their visual focus on and reached with the stylus to a new target that appeared in the upper right or lower left of the screen. Trials were defined as a combination of the planning period, 0 to 500 ms after the end of the target onset window, during which the subject planned either the reach/saccade or lift, and subsequent execution from 500 to 1000 ms after the target onset window, during which the subject moved their hand and reached for the new target. Each subject completed 256 trials, 96 of which were control condition trials, with the remaining 192 being test condition trials. The position of the new target randomly varied by a small angular amount (± 5 degrees) about a 30 degree direction centered at the position of the initial target—the center of the screen. Detailed description of the cognitive task performed by the participants can be found in Park et al. (2014).

#### EEG Recording

Scalp electroencephalographic (EEG) activity was recorded with a sampling rate of 512 Hz from 64 electrodes positioned on the scalps of the participants in an extended International 10–20 system configuration. The electrode array was grounded/referenced on a standard Driven Right Leg/Common Mode Sense (DRL/CMS) reference (Biosemi Inc., Amsterdam, The Netherlands), and each electrode was individually impedance-tested to ensure low impedance across the entire electrode array (below 5 k for each electrode). The 3D position of each electrode was digitized in order to construct individualized head models for each participant (see Park et al., 2014 for details on the digitization system).

#### Preprocessing and Artifact Rejection

Scalp electroencephalographic (EEG) channel space data preprocessing was conducted using EEGLAB (Delorme and Makeig, 2004) through a procedure outlined in **Figure 1A**. EEG data from each of the 10 participants were band-pass filtered,

was monitored by 64-channel EEG. Participants held stylus over a white fixation point in the center of the screen (Top Left) and were prompted for each trial to either lift and retouch the stylus (Control Condition) or to reach and saccade (Test Condition) to a randomly generated green target (Bottom Left). Epoch structure (Bottom Right) was identical for control and test conditions, with planning conducted while the participants remained on the fixation point with the green target present, and execution beginning when the fixation point was removed. Analysis window includes both planning and execution periods of reach/saccade or lift. Direction of the reach (upper right vs. lower left) was not discriminated since our focus was on the general causal dynamics involved in reaching rather than the encoding of spatial directionality.

retaining frequencies from 1 to 55 Hz, and re-referenced to the channel average since the data was recorded reference-free. The data were epoched, extracting the time period of interest, which encompassed both the planning phase and the execution phase of every trial. No distinction was made between reach directions since the emphasis for this analysis was placed on the general process of reach/saccade planning and execution rather than directional discrimination and, as such, both upward and downward angled reach trials were pooled and analyzed together. The control trials (Lift) and test trials (Reach/Saccade) were analyzed separately, using a cortical network defined over the same ROIs. Differences between the control and test conditions were examined by considering differences in network dynamics between the control and test conditions during the planning and execution phases. Noise and artifact laden trials and channels were identified and removed in a semi-automated manner. The kurtosis of each channel or trial was used as the identification metric, with channels or trials exhibiting a kurtosis Z-score > +5 or < −5 being flagged and subsequently visually inspected before removal.

The epoched channel space data were decomposed into sets of maximally independent components, where independence was achieved by minimizing mutual information between components using Infomax ICA (Bell and Sejnowski, 1995) (process outlined in **Figure 1B**). These components represented a combination of putative effective cortical sources, muscular artifacts, ocular artifacts, and electrical activity from the heart, among other noise sources. Through a semi-automated procedure, cortical sources were defined by employing independent components that both visually, by identifying power spectra and topographical maps with dipolar cortical source character, and computationally, through dipole-fitted residual variance, corresponded to dipoles (each fitted to a dipole with low residual variance). This procedure is particularly well suited for EEG artifact removal since the electrophysiological activity of interest present in the cortex manifests detectably in the form of an electrical dipole. The linear mixing of the electric fields generated by these dipoles, and corresponding unmixing through ICA allows for very precise retention of cortical data and rejection of non-cortical artifacts. Using the weight matrix computed during ICA, the retained cortical independent components were re-projected to channel space, creating cleaned channel space electrical signals that were further downsampled to 128 Hz for use in subsequent processing.

#### ROI Selection Methods-Clustering

Cortical regions of interest (ROIs) for source localization were determined by spatially and functionally clustering the ICs aggregated across all datasets/participants into a common cortical volume, as outlined in **Figure 1C**. The clustering feature vector for each IC was constructed by combining information from: the topographical scalp map of the IC, the event related potential (ERP) associated with the IC and its power spectrum, the spectral power content of a rolling window Fourier Transform (FT) across the IC time-series data before ERP-type averaging, known as the Event-Related Spectral Perturbation (ERSP) of the IC, the associated phases from this rolling FT window that give the Inter-Trial Coherence (ITC) of the IC, the ERP-image constructed from concatenation of all trials for the epoched event, and the three dimensional location of the equivalent current dipole fitted to the IC (Makeig, 1993). These features were selected to ensure that the IC clusters were both spatially tight and functionally homogeneous, thus promoting the creation of clusters from the same cortical region that demonstrated similar spectro-temporal dynamics during the conducted task. PCA was used to reduce both the dimensionality of each of the above features, as well as the final feature vector for each IC to 10 dimensions corresponding to the 10 highesteigenvalues, and, subsequent, K-means clustering of the ICs was conducted with K = 14 groups, whose clustering outcome yielded more spatially tight and functionally homogeneous clusters than the outcomes for K = 9 and K = 16 groups. The centroid of each cluster was computed and anatomically labeled using the Talairach atlas (Lancaster et al., 2000).

#### Head Model Construction

The (BEMs) that defined the geometry of the cortical source space for current localization were constructed using MoBILAB (Ojeda et al., 2014), as outlined in **Figure 1D**. The BEM approximation for the cortex of each participant included a standardized head model—the MRI MNI Colin27 three-layer BEM comprised of scalp, skull, and cortical layers. The cortical layer, where the source localization n occurs, was discretized into a triangular mesh containing 4825 vertices, which were automatically aggregated into anatomical regions based on the MNI standard gray matter atlas, used for ROI selection. A customized head model was constructed for each of the 10 participants based on the digitized electrode positions on EEG cap conforming to the participant's scalp during the experiment. Using standard electrical conductance values for the three BEM layers (0.33 S/m Scalp, 0.022 S/m Skull, 0.33 S/m Cortex Oostendorp et al., 2000), an electrical lead field matrix was derived, relating the conduction of electric fields generated by current dipoles from the cortical mesh (bottom layer) to the top of the skin (outer layer) where the field potential is recorded by scalp electrodes. This lead field matrix, along with the discrete surface Laplacian of the cortex layer of the BEM, were computed using the OpenMEEG toolbox (Gramfort et al., 2010) which was also used later for computing the current density during source localization.

### Source Localization-Cloreta

Cortically Constrained Low Resolution Electromagnetic Tomography (cLORETA) (Pascual-Marqui et al., 1994) was used for source localization due to the fixed source-space it required for current localization, allowing for cross-participant ROI consistency. The source-space was the cortical BEM constructed in the previous section that accounted for the conduction of electric fields generated in the cortex and projected to the surface of the scalp where they were detected by EEG. The cortical constraint placed on the BEM was also necessary due to the inability of EEG to reliably resolve electrophysiological activity in subcortical regions of the brain. Source localization (procedure outlined in **Figure 1E**) was conducted independently for every epoch in each dataset in order to prevent regularization parameters and current activations from being computed by concatenating temporally non-contiguous epochs. The solution to the minimum-norm Lagrange multiplier problem for cLORETA is the Tikhonov regularized pseudoinverse solution for cortical current density expressed as follows:

$$J = T\Phi\tag{1}$$

$$T = \left(K^T K + \lambda^2 H^T H\right)^{-1} K^T \tag{2}$$

where J is the current density on the cortex, Φ is the channel space data, **T** is the optimal transformation matrix from channel space to current density space, **K** is the discrete electrical lead field matrix, and **H** is the 2D spatial Laplacian of the cortical head model, introduced to spatially regularize the localized current on the head model grid points. λ is a regularization parameter initialized on a log-space, using generalized cross validation (GCV), and subsequently inferred using an iterative EM formulation derived from the Bayesian structure of the original current localization problem (Trujillo-Barreto et al., 2004). The cLORETA implementation in BCILAB (Delorme et al., 2010) was used for source localization.

Since the regularization parameter λ was recomputed for every epoch, a new transformation matrix (**T**) was also computed, leading to dynamic increases and decreases of spatial regularization in the final current solution as needed. This prevented both over-regularization in cases where current density ramped quickly between grid points, and underregularization when channel space, and therefore current density magnitudes, were lower. A representative signal was extracted from each ROI by computing the signal power of the current localized at each gridpoint within a given ROI and selecting the gridpoint that exhibited the maximum power.

Based on this max-power selection heuristic, representative current signals were generated for the ROIs—a process which was applied to each of the 10 subjects individually, greatly simplifying network-level comparisons across test subjects since the nodes of the cortical network (the ROIs) were localized in the same location every time. This approach differs from traditional ROI determination methods involving ICs because the latter rely heavily on the common spatial localization of ICA outcomes, which is not guaranteed for multiple separate runs of ICA on different data, and, thus, may lead to inconsistencies for group level connectivity analysis (Thompson et al., 2011).

#### Connectivity Analysis–Sift

Application of the maximum power heuristic on the results of cLORETA created a single time series that represented cortical current dynamics for each ROI selected through clustering. These time series representing ROIs (ROI signals) were used to construct linear (MVAR) models from which spectro-temporal causal interactions between ROIs were computed using the Source Information Flow Toolbox (SIFT) (procedure outlined in **Figure 1F**) (Delorme et al., 2011). Specifically, for each participant, the ROI signals were preprocessed by employing detrending and normalization across time and ensemble. The MVAR models were constructed by applying the Vieira-Morf lattice algorithm to the ROI signals, using a sliding window with a length of 550 ms and a window step size of 10 ms, leveraging the window size to decrease the number of the estimated MVAR coefficients (ratio of data points to parameters > 10:1)—a necessary provision arising from the lack of regularization in the Vieira-Morf algorithm. The optimal order for the MVAR model was determined automatically by averaging the optimal orders estimated by six different information criteria: Akaike Information Criterion, corrected Akaike Information Criterion, Hannan-Quinn Criterion, Schwarz-Bayes Criterion, Rissanen Prediction Error, and logarithm of Akaike's Final Prediction Error (Rissanen, 1978; Tu and Xu, 2012). Each of these information criteria was computed for the range of possible model orders between 1 and 30, and the optimal model order for each participant was automatically determined by identifying the minimum value for each of the six generated curves, computing the average, and rounding down to the nearest integer. The model order for the Lift and Reach/Saccade conditions were determined independently. After model order determination, tests of stability and consistency were performed for each MVAR model in order to validate the ability of the model to capture the dynamics of the ROI signal. The autocorrelation function and three different portmanteau tests were used to determine the whiteness of the model residuals. The percent consistency of the model was evaluated for every sliding window within the epoch, indicating, as a percentage, the relative amount of correlational structure within the data that had been captured by the model. A stability test was conducted on every window within the epoch the MVAR coefficient matrix was evaluated.

Pairwise causal relationships between ROIs were computed from the validated MVAR models by employing the (SdDTF) (Korzeniewska et al., 2008). The SdDTF connectivity η 2 ij f , t at a given spectro-temporal point f , t between ROI-i and ROI-j is computed as follows:

$$\eta\_{j \to i}^2 \left( f, t \right) = \frac{\left| H\_{ij} \left( f, t \right) \right| \left| P\_{ij} \left( f, t \right) \right|}{\sum\_f \sum\_{kl} \left| H\_{kl} \left( f, t \right) \right|^2 \left| P\_{kl} \left( f, t \right) \right|^2} \tag{3}$$

where Hmn (f , t) is the MVAR transfer matrix between network nodes (ROIs) m and n, and Pmn (f , t) is the partial coherence of the signals for nodes m and n. The SdDTF connectivity metric is similar to conditional spectral Granger Causality, but only reports direct causal relationship between signals, removing spurious or indirect relationships caused by upstream network topology, particularly common inputs. Unlike Granger Causality, which utilizes AR model residuals to determine causal influence, SdDTF utilizes the AR transfer functions themselves (see Equation 3) allowing for straightforward generalization of the metric to multivariate systems, in this case a network with more than two cortical nodes. To control for the potentially large and unregularized AR transfer functions, DTF normalization is conducted across incoming and outgoing signals, and across frequencies. This normalized term is weighted by the partial coherence in order to incorporate the direct relationships between source nodes, thus mitigating the phantom source problem present in Granger causal metrics and yielding the full form of the direct Directed Transfer Function. This metric is further developed by evaluating it over a short sliding window across the epoch, an additional feature which serves to promote local stationarity of the ROI source signals and observe changes in connectivity dynamics between signals as the epoch progresses. The SdDTF measure was computed separately across epochs and frequencies of interest in this way, yielding a spectrotemporal asymmetrical causal connectivity matrix of SdDTFs with causal source nodes (FROM) defining the columns and sink nodes (TO) defining the rows of the connectivity matrix. Note that the diagonal elements of this matrix do not have relevant physical meaning for the purposes of this network analysis, and were, therefore, omitted.

A permutation-bootstrap approach was employed to determine the statistical significance of each SdDTF value within the asymmetric connectivity matrix (outlined in **Figure 1G**). The permutations were created by randomly mixing ROIcurrent density epochs from the control and test conditions with equal representation, generating an epoch ensemble on which connectivity analysis was conducted. The resulting ROI connectivity matrices from each permutation were accumulated to create a surrogate distribution. The values of each spectrotemporal position in this distribution were, then, sorted in ascending order, allowing for the associated p-value of the tested connectivity value to be computed in a non-parametric manner through its place in the ranked connectivity values at that spectro-temporal position of the surrogate distribution. This process was applied to both the Lift and Reach/Saccade conditions, where the null hypothesis (H0) was the lack of significant increase or decrease in connectivity between the evaluated condition and the case of randomly mixed ROIs. The alternative hypothesis (H1) was that the evaluated condition individually exhibited connectivity structure that was significantly different from the structure acquired by randomly intermixing the two conditions. The p-value for each spectrotemporal point was determined by using an averaged surrogate distribution computed across all 10 participants. A surrogate distribution for each of the 10 participants was computed by conducting connectivity analysis on 300 permutations of the original data, yielding 300 SdDTF values for each spectrotemporal position. The average surrogate distribution was computed by applying spectro-temporal position-wise averages across all 10 surrogate distributions for each of the 300 values and, subsequently, sorting in ascending order the resulting average values at each spectro-temporal position. With 300 permutations computed for each subject, there was sufficient pvalue resolution to establish which time-frequency points passed the α = 0.05 threshold. This threshold was automatically adjusted since the statistical significance of each spectro-temporal point was evaluated independently by applying the False Detection Rate (FDR) control method as a means of accounting for the multiple comparisons problem (Benjamini and Hochberg, 1995).

The difference between statistically significant values of average connectivity matrices corresponding to the Lift condition and Reach/Saccade condition was computed in order to determine the statistically significant changes in the cortical network between the two conditions.

### RESULTS

### IC Clustering Outcome

The ROIs selected for source localization based on the cluster centroids were: bilateral (left/right) inferior occipital (Occ), superior parietal (Par), and precentral gyrus motor (Mot) cortices, as well as, anterior cingulate cortex (ACC), supplementary motor area (SMA), and precuneus, whose left and right hemisphere grid point regions were combined and treated as a single anatomical region during source localization and ROI collapse.

**Figure 3** depicts the clustering outcome for the reach/saccade condition, showing the dipole cluster and average scalp projection for the eight clusters from which ROIs were determined. The 3D dipole maps plotted on a standard MNI MRI volume demonstrate the spatial tightness of the clusters selected as candidates for ROIs, while the average composite topographical maps indicate that there is consistency in both the orientation and di-polarity of the ICs in each cluster. All of the selected clusters feature a relatively tight positive spike in electric potential over the cortical area where the equivalent dipoles are placed; suggesting that the collection of dipoles within the cluster originate from similarly oriented pyramidal cells in gray matter, generating electrical dipoles whose resultant field produces the positive potential present in the scalp maps. Were the IC clusters not as spatially tight, or comprised of dipoles oriented in different directions, the average scalp maps would have electric potential much lower in magnitude and dispersed over a much wider area of the brain. The spatial tightness of the dipoles and the scalp maps suggests that a robust clustering outcome was achieved, and the associated ROIs were reliably utilized with this data in subsequent analysis.

The determination of ROIs from ICs is not necessarily a one-to-one mapping. A single cluster can provide evidence for inclusion of multiple ROIs, as is the case with cluster 5, for example, whose dipole locations span across both hemispheres. So, when identified with the Talairach atlas, the ROIs fall in bilateral regions of the Occipital Cortex (Inferior gyri) that are too far from each other to be grouped across the medial longitudinal fissure, as was the case with all of the other midline ROIs.

All nine ROIs selected for the cortical network are displayed in **Figure 4**, with color-coded grid-points associated with each ROI as follows: midline ROIs are colored Red (ACC, SMA, Precuneus), and bilaterally (left and right) situated ROIs away from the midline (Mot, Par, and Occ) are dark blue, light blue, and yellow respectively.

## MVAR Model Construction and Validation

Model order selection was conducted with an average MVAR parameter-to-datapoint ratio of 0.096 for the Lift condition and 0.034 for the Reach/Saccade condition. The window size was selected so that, in all cases, this ratio never exceeded 0.1, as is recommended for a well-posed MVAR model computed via the Vieira-Morf Lattice Algorithm. For the Lift condition, the average model order across all subjects was 7.4 with a standard deviation of 1.35, and for the Reach/Saccade condition, the model order

average was 7.2 with a standard deviation of 1.03. **Figure 5** shows the output of the model order selection criteria pertaining to the MVAR model order of the Lift condition for a single participant. The top graph shows the average (across epochs and sliding windows) value of each criterion over model orders from 1 to 30. The histograms that follow show the order distribution using each of the criteria.

The MVAR model validation results averaged across all subjects and both conditions are shown in **Figure 6**. The average model performance was highly representative of the individual

performance of each subject's MVAR model since, in all cases, the MVAR model was stable across the entirety of the epoch, as indicated by the natural logarithm of the largest eigenvalue being negative (**Figure 6** Stability Panel). The percent consistency for both Lift and Reach/Saccade conditions was above 90% and the standard deviation in the percent consistency was <1.0% (0.877% for Lift and 0.712% for Reach/Saccade), indicating that the high consistency of the MVAR model with respect to the data was maintained for the duration of the epoch (all 43 sliding windows within each epoch). None of the computed MVAR models passed the residual-whiteness tests at any point during the epoch, suggesting that there was statistically significant correlational structure exhibited by the MVAR model residuals, and that the MVAR modeling did not succeed in capturing the entirety of the temporal dynamics exhibited by the network nodes during the epoch.

#### Spectro-Temporal Connectivity

The connectivity computed between the ROIs in the defined cortical network is displayed in **Figure 7** for the lift condition and **Figure 8** for the reach/saccade condition, with statistically significant spectro-temporal connectivity (p < 0.05) thresholded for each matrix respectively.

The difference between the p-value thresholded connectivity matrices for Lift and Reach/Saccade, shown in **Figures 7**, **8** respectively, was computed and shown in **Figure 9**. Since the lift (control) condition was subtracted from the reach/saccade (test) condition, the warmer (brighter red) regions of the matrix indicate that connectivity was stronger during the reach/saccade and the cooler (brighter blue) regions indicate that connectivity was stronger during the lift.

The majority of lift-dominant causal interactions originate from the supplementary motor area, parietal areas and the Precuneus. During Lift, the SMA exhibits notable outflow to the Left Parietal Area (L Par) and the Precuneus. The SMA drives activity in the L Par over a constant frequency range (8–15 Hz) during both planning and execution (**Figure 9A**). However, a dichotomy between planning and execution exists in the connection form the SMA to the Precuneus, with sporadic low-frequency connectivity (3–14 Hz) present during planning that narrows in frequency band (8–14 Hz) and becomes more consistent during execution (**Figure 9B**). From L Par, the only significant outflow of information occurs to the Precuneus, dominating lower frequencies (3–14 Hz) during part of the early planning phase and the entirety of the execution phase of the lift (**Figure 9C**). From R Par, information flows to the ACC across low frequencies (3–14 Hz) during the planning phase and for some of the execution phase, though the causal connection is limited to a narrow frequency band (10–14 Hz) by the end of the epoch (**Figure 9D**). Similar connectivity is present with the SMA, although in the case of information flow from R Par to SMA, the connectivity does not subside as late in the execution phase as in the case with the ACC (**Figure 9E**). The R Par exhibits another connectivity motif, providing strong causal influence on both the Left Parietal Area and the Precuneus during the planning phase of the lift over a low frequency band (3–10 Hz), and with

relatively higher power compared to its other causal influences (**Figure 9F**). The Precuneus exhibits the greatest amount of causal outflow during lift, acting as an information source for the majority of the other ROIs in the network. The Precuneus has relatively lower-power (lighter blue shading) influence on both precentral gyri, during both planning and execution, with the band of active connectivity narrowing when transitioning from planning to execution (**Figure 9G**). The Precuneus sends information with a similar spatiotemporal structure to the ACC, but with slightly greater causal power, hence the brighter shading of the connectivity during the planning period (**Figure 9H**). The influence the Precuneus has on the SMA is stronger still (even darker shade of blue), and now occupies the entire duration of the epoch and maintains signal bandwidth (3–14 Hz) (**Figure 9H**). The Precuneus influences both parietal areas very strongly during lift, with R Par receiving information only during the planning phase across a low frequency band (3–10 Hz) (**Figure 9I**), and L Par receiving information across a slightly larger band (3–14 Hz) and for both planning and execution. During execution, however, the connectivity is limited to the higher end of that frequency band (10–14 Hz) (**Figure 9J**).

Within reach/saccade-dominant activity regions, which are defined by the red spectro-temporal regions of **Figure 9**, greater connectivity is present in the frontal regions of the cortex. Specifically, both the left and right precentral gyri exhibit stronger connectivity to the left parietal composite area during the planning phase of the epoch across lower frequencies (3–14 Hz), with this connectivity decreasing in power relative to lift during the execution phase (**Figure 9K**). Additionally, the anterior cingulate cortex (ACC) exhibits increased outflow of information across lower frequencies (3–14 Hz) to both precentral gyri again during the planning phase (**Figure 9L**). The ACC also appears to feed information to the left and right visual cortices across a narrower frequency band (8–14 Hz) during the planning phase (**Figure 9M**), with the bandwidth of information flow encompassing lower frequencies from ACC to R Occ as the epoch progresses from planning to execution. For the reach/saccade condition, the SMA exhibits a greater amount of causal influence on the ACC than any other cortical node, with a strong 14–40 Hz communication band present during the planning phase that dips down to lower frequencies (3–8 Hz) as the subject transitions from planning to execution (**Figure 9N**). The ACC reciprocates some of this activity to the SMA, with higher frequency (15–40 Hz) communication during the planning phase that dies down during the execution phase (**Figure 9N**). The SMA also exhibits some information outflow to the visual cortices, with the L Occ flow subsiding by the end of the planning period, and the R Occ flow remaining constant in both frequency band and relative amplitude over the entire epoch (**Figure 9O**). The right parietal composite area (R Par) exhibits causal influence on the Precuneus in higher frequencies (20–25 Hz) during the planning phase, and significantly stronger lower frequency (3–13 Hz) influence during reach/saccade execution in comparison to the lift execution (**Figure 9F**). The

last of the significant communication motif that is strong during the reach/saccade conditions originates from the Precuneus. The influence of Prec on L/R Occ is most prevalent during the execution phase (**Figure 9P**). In both cases, however, the precuneus predominantly drives lower frequency activity (3–14 Hz), with a significant decrease in connectivity magnitude above 14 Hz. There is also evidence of temporally localized higher power bursts of connectivity across low frequency from the Precuneus to the right parietal area, beginning early in the planning period in the 14–20 Hz band and terminating at the end of the epoch after drifting down to 8–13 Hz (**Figure 9I**).

### DISCUSSION

#### Neurophysiological Connectivity

A cortical network associated with the Lift and the Reach/Saccade conditions was identified across 10 subjects, involving distinct ROIs, and the dynamic causal interactions between pairs of ROIs were computed for each condition. The cortical network, the ROIs, and the dynamic causal connections between pairs of ROIs exhibited physiological plausibility consistent with previous results reported in the literature, although some important differences were noted as outlined further below. The bilateral Inferior Occipital Gyri coupled with adjacent Occipital regions, known to form the Visual ROIs, represent primary processing of the visual stimuli present in the reaching task (Blumberg and Kreiman, 2010) . The various motor regions, including bilateral Precentral Gyri (housing the Left and Right Motor Cortices) and the bilaterally grouped Supplementary Motor Area, have been implicated in the planning and execution of the physical movement, e.g., in this case of the participant's arm to reach the presented target as well as the saccade of the eyes to focus on the target (Johnson et al., 1996; Nachev et al., 2008; Wang et al., 2010). Parietal activity is largely implicated by existing literature in guiding action directed at visual targets (Blangero et al., 2009; Filimon, 2010). The Anterior Cingulate Cortex, a grouped ROI consisting of the Left and Right cortices merged across the longitudinal fissure, is thought to be involved in higher level decision making and, of particular interest in this case, the translation of intentions to actions (Kennerley et al., 2006). The Precuneus, also bilaterally grouped across the longitudinal fissure, has been shown to be directly involved in both visual and non-visual reaching tasks in both macaque monkeys and humans, being responsible for the processing of both visual and proprioceptive sensory information during the task (Filimon et al., 2009).

Though entirely possible to treat the three grouped ROIs (SMA, Prec, and ACC) as distinct, the close proximity of the bilateral ROIs to each other within the longitudinal fissure and the spatially-regularized source localization can lead to significant overlap of posterior current density between the ROIs, generating spurious causal interhemispheric coupling, which is avoided by considering those three ROIs to be trans-hemispheric.

Inspection of the connectivity during the lift and reach/saccade condition, shown in **Figure 9**, reveals that the cortical network exhibits relatively lower strength, but

widespread effective connectivity during the reach/saccade condition, and more focal, higher-strength connectivity during the lift condition. In the case of the reach/saccade condition, the ACC exemplifies this widespread, low-power connectivity profile, driving activity in both Left and Right Precentral Gyri during the planning period, and driving activity in both visual cortices for both planning and execution of the reach. The connection from the ACC to L/R Mot (**Figure 9L**) suggests a causal influence of a decision making center in the frontal cortex onto regions directly involved in the planning and execution of a reach and saccade, a connectivity motif that has been observed before using other neuroimaging modalities (Asemi et al., 2015). Furthermore, this causal influence is shown only to exceed lift-connectivity magnitude during the planning period, indicated by the lack of ACC to L/R Motor interaction after 0.5 s, suggesting that after the movement is planned, execution no longer requires direct influence of the ACC onto the Motor cortices. On the other hand, the ACC to L/R Visual connections (**Figure 9M**) remain elevated for the entirety of the epoch during the reach/saccade condition, indicating a possible involvement of the cognitive-control and error detection functions of the ACC (Stevens et al., 2011) in visual processing during execution. The task-positive connections formed by the ACC exhibit higher frequency connectivity (**Figures 9 L–N**), with the lower end of these connectivity bands lying in the α (8–13 Hz) wave range, which is expected considering the motor-relevance of the α-band to indicate inactivation of motor regions (Brinkman et al., 2014), and the higher end of these connectivity bands stretching into the 30–50 Hz range, consistent with the generally higher-frequency content of cognitive cortical processes (Buzsáki and Silva, 2012).

However, this does not seem to be the case for the information flow from the SMA to the ACC (**Figure 9N**), wherein a high α/low β band connection (10–20 Hz) present during the planning period migrates down to lower frequencies (3–8 Hz), occupying the θ band by the end of execution. A possible explanation for this exception stems from the suggested involvement of the SMA in the real-time control of actions during execution, by which it might provide feedback to the ACC regarding the accuracy of the conducted reach (Nachev et al., 2008). The connections from L and R Mot to L Par (**Figure 9K**) exhibit a connectivity profile similar to the ACC → L/R Mot interaction, with lowerpower, reach/saccade-dominant connectivity during planning that is absent during execution. While both the composite parietal regions and the precentral gyri have been implicated in visuo-spatial processing and reaching to identified targets, the analysis indicates that their interaction is limited during the planning phase of a reach, and is not required in order for the reach to be executed after planning. A strong, sharply divided connection is also observed from the Right Composite Parietal Region to the Precuneus (**Figure 9F**)—a connection which is facilitated during lift-planning and reach/saccade-execution, and suppressed during reach/saccade-planning (for low frequencies) and lift-execution. The connectivity profile suggests that during reach/saccade planning, R Par predominantly communicates with the Precuneus through higher frequencies (14–30 Hz), which is of great significance to the Precuneus since it normally communicates with other nodes in the Default Mode Network through this frequency band (Neuner et al., 2014). The transition to lower frequency connectivity (3–14 Hz) during Reach/Saccade execution may simply be a feature of task demands, since similar communication was seen in **Figures 9H,I,M,N,P**, suggesting a generalized tendency for visuo-spatial and motor related cortical regions to communicate through lower frequency bands when the reach/saccade is being executed. The reach/saccadedominant dynamics of the connection from the Precuneus to R Par (**Figure 9I**) also exhibit characteristics similar to the other motor regions, where the high-frequency planning connectivity gives way to lower frequency execution connectivity. However, during the epoch, the frequency range of the connectivity varies periodically, dipping to lower frequencies and, subsequently, returning to higher frequencies (**Figure 9I**, red profile) at the rate of 5 Hz. This behavior suggests that the Precuneus dynamically alters the frequencies over which it communicates with R Par in a periodic manner. The persistence of this behavior through both planning and execution suggests that the behavior is associated with generalized motor-reach information processing, though the decrease (**Figure 9I**, red profile) at the transition between planning and execution at 0.5 s supports the involvement of this periodicity in the motor task being executed by the subject.

Connectivity present in the cortical network is summarized in **Figure 10**, where the Reach/Saccade condition is shown in a 3D (BrainMovie3D) rendering of a cortical volume (Delorme et al., 2011). The ROIs are represented by labeled spheres of varying size and color, with the magnitude of causal information inflow indicated by color and the magnitude of causal information out-flow indicated by size. Larger size spheres indicate larger information out-flow and warmer colored spheres (brighter red) indicate larger information in-flow. Connections between ROIs are represented by tapered cylinders with varying color and diameter. Cylinder tapering indicates the direction of information flow, while cylinder width indicates flow magnitude, and cylinder color represents the dominant frequency of communication, with cooler colors for lower frequencies and warmer colors for higher frequencies. Causal information flow is integrated across all computed frequencies (2–50 Hz), and is shown here (**Figure 10**) for t = 250 ms, that is, during the planning period. The 3D rendering demonstrates the dominant role of the Precuneus as an information source, indicated by its large size, the higher power/frequency coupling of frontal regions (ACC, SMA, L/R Mot), indicated by warmer colored and larger connecting cylinders, and the higher frequency communication present in the network during reach planning.

#### Framework Limitations

The above discussion of neurophysiologically relevant findings serves to demonstrate the nature of the conclusions that can be drawn through the application of this algorithmic framework. Application of this method yields findings that are both consistent with existing literature, and reveal connectivity information in the cortex with an immense degree of detail, allowing for the documentation of source-space spectrotemporal information flow. There are, however, limitations with the current approach which must be considered.

The source orientation assumptions applied to the cortical BEM for the source space held that all cortical dipoles were oriented normally to the surface of the BEM, a canonical assumption that has been previously employed for BEM construction (Fuchs et al., 2002), and is motivated by the perpendicular orientation of cortical pyramidal cells to the cortical surface. Applying this assumption simplifies solving the inverse problem by reducing the number of variables to be estimated to Ng (number of BEM grid-points) from 3 ∗ Ng , as is the case when accounting for the x, y, and z components of the current at each grid-point. However, this assumption greatly restricts the form of the computed current density, and does not allow for any data-driven flexibility with the computed dipoles. Utilizing the position-free BEM and associated Lead Field Matrix would allow for much greater flexibility, allowing for cortical dipole orientation to be computed either jointly with the current density, or after source localization using the current density, and thus identifying orientations through observed data instead of through a prior model.

Another potential drawback of the current source-localization algorithm is the Laplacian regularization utilized by cLORETA. The smoothing features of the discrete Laplacian when used as a Tikhonov matrix cause the estimated current density to be blurred over large regions of the cortical surface. This blurring can cause current density from a single source to spill over into multiple adjacent ROIs and potentially skew causal interactions between these ROIs in a manner that is not corrected by using the SdDTF connectivity metric. Preliminary investigations using this data suggest that this spill-over does not have a significant effect on causal analysis in this proposed framework, but further investigation would be prudent to ensure that the computed connectivity dynamics are not skewed by the choice of cLORETA as the source localization algorithm.

The validity of the presented connectivity dynamics is supported by the statistical analysis of the SdDTF connectivity matrix and the validation conducted on the fitted MVAR models. Though a large amount of the computed connectivity was statistically significant, the residuals of the computed MVAR models failed the whiteness tests because not all correlational structure of the localized current was captured by the models, suggesting the existence of non-linear dynamics. These nonlinear dynamics are readily visible in the model prediction error and suggest that non-linear AR modeling is required to effectively capture these dynamics in the MVAR coefficients. However, all models exhibited a large Percent Consistency in the order of 90%, indicating that a very large percentage of the correlational structure of the localized current in each node was captured by the MVAR. The high percent consistency of the model indicates that these nonlinearities constitute a small portion of the localized current, and the computed MVAR models are able to account for the vast majority of the observed dynamics.

### Extensions to Real-Time BCI Applications

In spite of the potential limitations discussed in the previous section, the proposed algorithmic framework has great potential as a means for extracting connectivity features in real time that can be used for brain state identification and BCI control. While in its current form, the framework is designed for offline processing, Mullen et. al have demonstrated in Mullen et al. (2015) that elements of this framework can be applied in a real-time BCI setting and, with high accuracy, decode brain states using the spectro-temporal connectivity features computed through the use of this framework. The linear nature of the inverse problem formulated through cLORETA allows for rapid localization of current density in the cortex by application of the closed-form solution given in equation 2. Furthermore, the lattice algorithm used for MVAR coefficient estimation, while well suited for off-line analysis where all data is available for fitting, can be replaced by alternative estimation algorithms which take advantage of the convexity of MVAR model-fitting to compute robust models quickly and with relatively few data-points. In addition, the development of online recursive ICA (Hsu et al., 2014) allows for robust online artifact rejection, source identification, and lends to the overall feasibility of this approach in real-time BCI applications. There are other elements of the pipeline, such as the data-driven identification of ROI's for connectivity analysis, which must be adapted to account for the lack of all data present, or meet the computational demands in real-time analysis. Nevertheless, there is great potential use for this framework in both a clinical and research settings for real-time BCI control.

### CONCLUSION

The presented algorithmic and computational framework successfully identified a cortical network, across several subjects, associated with a cognitive task, using EEG signals recorded from individuals performing the task with the aid of a BCI. It also computed statistically significant causal dynamic connectivity profiles among the ROIs comprising the nodes of the identified cortical network. The modular nature of the presented framework allows for flexibility and scalability. The algorithm computing the cortical network can easily scale up to include a greater number of participants without compromising the reliability or accuracy. The ease of incorporating large amounts of data comes as a consequence of utilizing a distributed source localization method with a BEM-defined source space. The flexible-yet-fixed source space geometry allows for warping head models to align with a realistic 3D representation of each participant's cortex by way of the digitized electrode positions, while retaining the identity of anatomical ROIs and their associated gray matter grid points in the head model. Should the user deem it necessary, the modular nature of this framework allows for substitutions of preprocessing algorithms, source localization methods, and connectivity analysis tools as is seen fit. As long as a source-space based localization algorithm is utilized to maintain consistent anatomical regions across subjects, the user is free to create a customized pipeline that best suits their needs, and is not strictly limited to usage of the algorithmic pipeline with the proposed software and algorithms. The proposed analytical approach allows for extraction of information rich connectivity features which account for spectro-temporal connectivity dynamics during the performance of a cognitive task. This feature extraction is made possible by coupling the segmentation-MVAR modeling of cortical signals with the SdDTF—a connectivity metric that is evaluated independently for different frequencies, leading to the asymmetric spectro-temporal connectivity dynamics displayed

### REFERENCES

in **Figures 7**–**9**. The presented framework also includes a nonparametric statistical test to assess the statistical significance of the extracted features—the spectro-temporal connectivity dynamics. The physiological plausibility of the cortical network and the causal dynamic features characterizing interaction among the nodes of the network computed by employing this framework strengthens the framework's applicability in reliably capturing complex brain functionality, which is required by applications, such as diagnostics and BCI.

### ETHICS STATEMENT

This study was carried out in accordance with the recommendations of the Human Subjects Institutional Review Board of the University of California, San Diego with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the Human Subjects Institutional Review Board of the University of California, San Diego.

### AUTHOR CONTRIBUTIONS

HC introduced improvements to the pipeline, conducted reprocessing and analysis of data, and wrote the manuscript. TM assisted with data processing on the first iteration of the pipeline. HP and GC contributed in an advisory manner with strategy and direction of the paper and the underlying computational methods. HP also directed the original project in which the data was collected, and GC assisting with editing the manuscript. JI (senior author) conducted processing and analysis on the first iteration of the pipeline, contributed in an advisory manner more closely than GC and HP, and assisted with editing the manuscript.

### FUNDING

This work was supported in part by the National Science Foundation (NSF) under Grant ENG-1137279 (EFRI M3C), Grant SMA-1041755 and Grant BCS-1460885, and in part by the Office of Naval Research (ONR) Multidisciplinary University Initiative (MURI) under Award No. N00014-10-1-0072.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.