where *k* is the number of frequencies of each set (in this case 2), and *fi* and *fj* are the individual frequencies taken from the union of the sets *x* and *y*. Note that here *x* and *y* refer to sets of two frequencies while in Equation (4) *x* and *y* referred to individual frequencies. This distance expresses the arithmetic mean of all possible pairs in the set resulting from the union of the initial sets *x* and *y*. Note that for *k* = 1, this distance measure is exactly the same distance (Equation 4) as used in the first part of the algorithm. In this second part 3*N/*8 iterations are performed, which is *N/*2 (the number of disjoint sets) times 3/4 (see above).

#### *Second closed-loop in the ACL-algorithm: online auditory feedback of SSVEP magnitudes*

In order to offer additional dynamic information to the subject related to his/her brain activity beyond the SSVEP detection confirmation cue, we provide a continuous online auditory feedback during the trials which represents the distance between the actual state and the pre-defined goal. The feedback signal consists of a 20 possible sinusoids with a range between 100 and 575 Hz which are updated every 0.25 s. The represented distance measure is defined as the difference between the EEG-SSVEP *signal to noise ratio* for the target frequency (Starget *f* ) and the threshold. Once Starget *f* has reached this threshold level, the auditory feedback is muted. Previously, subjects are instructed that their goal is to raise the pitch of the sinusoids as high as possible, and that after possible success their further goal would be trying to keep the sounds muted for 1.75 s; after this silence, the program automatically proceeds to the trial's next step. This kind of continuous auditory feedback aims to help subjects to learn to gain control *in their particular way* over SSVEP magnitudes by attracting their attentional resources to these voluntary attempts to increase self-regulation of their resonating brain states.

Concluding, there are two assisted closed loops in our system: the first one operates over the stimulation frequency set with the aim to directly improve the ITRs of each subject. This closed-loop informs the system about subject and environment specificities. The second one informs the subject about his/her brain activity in relation to the use of the interface and helps him/her to do so faster and more accurately. This closed loop works several times for each step of a trial.

#### **SSVEP DETECTION**

In order to reduce the experiment's complexity in terms of a reductionistic paradigm, we choose a simple SSVEP detection strategy in our study. During the *top* and *prefixed* frequency stimulation, the Starget *f* value is calculated every 0.25 s. If this value exceeds the threshold for 1.75 consecutive seconds, then this SSVEP is defined as "detected." The threshold value was set to 10 which reflects the observed noise flow (see **Figure 5**). To avoid longer waiting periods when the subject is unable to exceed the threshold, a time limit of 4 s is used, after which that step is considered as fault.

During the ACL, to favor SSVEP detection in case that the subject exceeds the threshold and more time than the 1.75 s is needed to be classified as "detected," there is a small modification in this protocol to allow adaptive time extensions. When Starget *f* exceeds the threshold in a given 0.25 s time step, the time limit is increased for another 0.25 s.

#### **EFFICIENCY MEASURES**

After each iteration of the algorithm, both the success rate and time needed are saved. For the *prefixed* and *top* frequencies, *standard Information Transfer Rate* (ITR) is calculated:

$$\text{ITR}(\text{SR}, \ t) = \left(\log\_2(N) + \text{SR} \cdot \log\_2(\text{SR}) + (1 - \text{SR}) \cdot \right)$$

$$\log\_2((1 - \text{SR})/(N - 1)) \cdot \text{Norm}/t \qquad (6)$$

where *N* is the number of targets (*N* = 4 in our case). The value *SR* represents the success rate and *t* is the time taken in minutes. *Norm* is a normalization value set to 960 (60 s times 16 steps in each iteration). Note that if *SR* ≤ 1/*N,* then ITR(*SR, t*) = 0.

In contrast to the conditions *prefixed* and *top,* ITR is measured *several* times during the ACL. Thus, for further a-posteriori analyses these ITR distributions have to be represented by descriptive statistics: for condition *ACL* therefore *M* and *Mdn* of success rates and needed times are used to calculate ITRMean and ITRMedian, completed by maximum ITR (ITRMax).

#### **CONVERGENCE MEASURE**

For a-posteriori analyses, a *convergence measure* for the algorithm in terms of the stimulus frequency exploration was defined: the duration of the 2 freq. search of the algorithm is divided into two parts. For each part, the numbers of explored frequencies are determined and divided by the maximal number of possible frequencies which could be explored (twice the number of iterations). The decrease comparing this measure in the second part vs. in the first part is a sign for how much the frequency exploration is converging. As can be seen in **Table 1**, the number of iterations varies over the subjects. The convergence measure is not reported for the first part because in our sample all subjects had the same maximal value 1, i.e., all possible frequencies were explored. We will use this measure to discuss how the ACL algorithm seems to adapt to subjects' interindividual differences.

#### **STUDY DESIGN**

A three conditions (*ACL, top, prefixed*) balanced *within*-subjects design with three times full permutation of presentation order (ABC, ACB, BAC, BCA, CAB, CBA) and with random assignment of subjects, resulting in *N* = 18 was employed.

#### **BASELINE RESTING STATE EEG MEASURES AS POSSIBLE INTERINDIVIDUAL CORRELATES OF ITR PERFORMANCES**

Aiming to investigate possible correlations between baseline resting state EEG measures and the variables of the experiment, the 30 s baseline EEG (see **Figure 3**) at all eight electrodes reported above were manually cleaned from artifacts with the result of *M* = 20*.*02 s, *SD* = 5*.*54 artifact free epochs.  Then, preprocessed EEG signals were converted into *frequency domain* by a sliding windows FFT transform of 2 s window length (2048 sample points) with 3.906 ms displacement (4 sample points, which correspond to a 256 Hz sample frequency in the resulting frequency domain signals), after linear detrending and treatment by a *Hann-*window function. Obtained FFT coefficients were squared to obtain the power spectrum and then normalized by dividing by 2048 sample points. In order to obtain *absolute PSD*s for the defined EEG frequencies bands of interest, corresponding coefficients were summed: *thetaLow* (3.5–6.5 Hz), *thetaHigh* (6.5–7.5 Hz); *alphaLow* (7.5–9 Hz), *alphaHigh* (9–12.5 Hz); *betaLow* (12.5–18 Hz), *betaMid* (18–24 Hz), *betaHigh* (18–30 Hz); *totalSpectrum* (0.5–70 Hz). In a first step, those absolute frequency domain *PSDs* signals were normalized dividing every sample point by the corresponding one of *totalSpectrum* which resulted in dimensionless ratios. These ratios indicate for every 256 time points per second the *relative* energy contribution of the frequency band of interest to the EEG total energy at this particular moment. In a last step, in order to represent EEG baseline resting state activities in the analyzed artifact free epochs by one single value for every frequency band, means of these normalized signals were computed over all corresponding time points. Thus, finally we obtained the desired baseline resting state EEG measures as *relative mean PSDs* for further correlational analyses*,* single values for every frequency band over all subjects.

Another measure of interindividual EEG variability is the resting state *individual alpha frequency* (IAF), because it has been found to be remarkably stable *within* subjects, but relatively variable *between* subjects (Kondacs and Szabó, 1999). In order to determine IAF in our experiment, coefficients of PSDs corresponding to the frequency band 8–13 Hz at *Oz* were normalized by *totalSpectrum* PSDs and averaged over all sliding windows in the artifact free baseline resting state epochs. In this averaged and normalized power spectrum the *alpha* frequency with the highest PSD was manually measured and defined as IAF (*peak frequency method*).

#### **STATISTICAL ANALYSES**

All statistical analyses were computed using SPSS 17.0 and STATISTICA 6.0. Previously, *Shapiro–Wilk* tests were calculated to check each of the three conditions for normal distribution in the underlying populations. If one or more conditions showed significant departures from normality, *non-parametric* tests were preferred for further analyses: a *Friedman* test was performed as an *omnibus* test to investigate whether the *central tendencies* of one or more conditions differed significantly from the rest. In case of such a significant result, *post hoc pairwise comparisons* were performed in order to find out what conditions exactly differed significantly from each other, based upon comparison of *mean*

*Note: Information transfer rates (ITRs) in bits/min as measures of individual BCI performances under the different experimental conditions and all Mdn values are highlighted in bold for further analyses.*

*N trials refers to the number of iterations in the first part of ACL (using two flicker LEDs).*

*Convergence measure first half is not reported in the table because all subjects had the same value 1.*

*SNR SSVEPs in Scanning phase are means over all used 20 flicker frequencies.*

*rank differences* using as significance criteria the *critical rank differences* proposed by the more progressive approach of Conover (1980) vs. the more conservative of Schaich and Hamerle (1984).

In order to quantify the *effect sizes* of those *post hoc* pairwise comparisons which resulted in significant differences, we used the probability of superiority of dependent scores, *PSdep*, recommended by Grissom and Kim (2012) and developed in Grissom (1994). It expresses the probability that in a randomly sampled matched pair the value from the condition containing the higher scores is indeed larger than that from the one containing lower scores. *PSdep* is calculated by dividing the number of *positive* differences between the condition *containing the higher scores* minus the condition *containing the lower scores* by the total number of matched pairs. For classifying *PSdep* into small, middle and large effect sizes based upon the standards of Cohen (1988), the *cut-off* values reported by Grissom (1994) are used: *small* 0.56, *medium* 0.64, and *large* 0.71. The same author offers a table to directly convert *PS* into equivalent *Cohen*'*s* . Thus, as effect size measures both *PSdep* and *Cohen*'*s* are reported with standards *small* = 0*.*20, *medium* = 0*.*50 and *large* = 0*.*80 (Cohen, 1988).

In order to check whether significant differences over all six possible permutations of the presentation order might be found, a *mixed-design repeated measures ANOVA* was computed with *stimulation condition* as repeated *within-*subjects factor with three levels (i) ACL algorithm represented as ITR*Median*, (ii) prefixed and (iii) top and *presentation order* as *between-*subjects factor with the six possible permutations as levels (ABC, ACB, BAC etc.). Previously, *Levene*'*s* tests were performed in order to check for homogeneities of error variance. Moreover, the assumption of *sphericity* of the covariance matrix was verified previously by a *Mauchly*'*s sphericity test* in order to assure that the *F* ratios match an *F* distribution. If there was a significant departure from sphericity*, Greenhouse-Geisser* estimates were used to correct degrees of freedom which results in fractions instead of usual integers. Although data may not follow a normal distribution, *ANOVA* has been demonstrated to be relatively robust against moderate deviations from normality (see e.g., Khan and Rayner, 2003). Univariate analyses were used to examine whether there is a significant *between*-subjects main effect of *presentation order* and further if there is a significant *interaction* effect between *presentation order* × *stimulation condition*. Analyses were repeated representing condition (i) ACL algorithm also as ITRMean vs. ITRMax.

For the investigation of linear correlational relationships, *Spearman's rank order correlation coefficient Rho* was additionally used apart from the common *Pearson product-moment correlation coefficient r* due to its relative robustness firstly against outliers, but also against other than linear, but still monotonic relationships and against departures from normality or homoscedasticity. Whenever relevant influence of outliers was suspected, *Spearman's rank correlation coefficient Rho* was preferred.

A-priori *statistical test power* analyses with the program *G*∗*Power 3* (Faul et al., 2007) show that *Pearson* correlation significance tests in the employed sample size of *N* = 18 and with standard significance level α = 0*.*05 have test powers (1 − β*)* ≥ 0*.*80 as recommend by Cohen (1988), when they have effect sizes in the underlying population ρ ≥ 0*.*60, as compared to *H*0 : ρ = 0*.*00. For ρ = 0*.*50 test power is (1 − β*)* ≥ 0*.*60, for ρ = 0*.*40 *(*1 − β*)* = 0*.*40 and for ρ = 0*.*30 *(*1 − β*)* ≈ 0*.*20. Thus, although the employed sample size *N* = 18 is relatively small, hypothesis testing of *Pearson* correlations with full recommended strictness is definitely possible at the level of assumed large effect sizes.

#### **RESULTS**

**Table 1** reports the data for all *N* = 18 subjects under the three experimental conditions, representing (i) *ACL algorithm* as ITRMean, ITRMedian and ITRMax. Inferential statistical hypotheses testing that (i) outperformed the other two flicker stimulation conditions is reported below.

**Figure 6** shows the SSVEP frequency-response curves in our experiments. For all subjects, the 20 flicker frequencies in the scanning phase were presented in the same order: 23, 37, 30, 31, 36, 22, 29, 33, 39, 24, 35, 21, 25, 27, 32, 34, 28, 20, 26, and 38 Hz. Sequential randomness of this order is confirmed with *Z* = −0*.*230 and *p*exact = 0*.*828 (*Wald–Wolfowitz* runs test after *Mdn* split dichotomization). Our findings that in the 20–39 Hz range, lower flicker frequencies *over all subjects* (**Figure 6A**) evoke higher SSVEP magnitudes are in line with other studies which reported a global maximum SSVEP amplitude around 10 Hz with additional local maxima around 20, 40, and 80 Hz (Regan, 1989; Herrmann, 2001; Bayram et al., 2011). In our sample, we found that SSVEP frequency-response curves differed remarkably *between subjects* (**Figure 6B**) probably due to trait and state variabilities which justifies that they are determined in our experiment in the scanning phase for every subject *individually.*

Analyzing **Figure 6C**, higher frequencies ≥30 Hz lead to higher correlations; no relevant differences can be seen comparing the three experimental conditions. Interestingly, following e.g., Zschocke and Hansen (2012), 30 Hz is the upper boundary of *beta* activity observable in scalp EEGs by conventional amplifiers.

### **SIGNIFICANT AND LARGE IMPROVEMENT OF SSVEP-BCI EFFICIENCY BY THE NOVEL ACL ALGORITHM**

Analyzing the differences in the *central tendencies* between the three experimental conditions (i) *ACL algorithm* (ii) *prefixed* (iii) *top* we represented condition (i) based upon three different descriptive statistics, (a) ITRMean, (b) ITRMedian, (d) ITRMax (see section "Materials and Methods" and **Table 1**). Applying nonparametric inferential statistics we found a very significant and very large superiority of condition (i) *ACL algorithm* over the other two (ii) and (iii) which is independent of its three types of representation (a), (b), and (c), while there is no significant difference between (ii) and (iii). The used statistical methods and measures for the following results are found in section "Statistical Analyses."

(a) A *Friedman* omnibus test comparing the ITRs between the three experimental conditions (i) *ACL algorithm* **represented as ITRMean**, (ii) *prefixed* and (iii) *top* shows a significant overall difference with χ2*(*2*)* = 10*.*116, *p* = 0*.*006.

*Post-hoc pairwise comparisons* based upon *critical mean rank differences* 0.82 (Schaich and Hamerle, 1984) vs. 0.58 (Conover, 1980) indicate that ITRs are significantly higher in (i) *ACL algorithm* as compared to (ii) *prefixed* (*mean rank difference* = 1.03, very large effect size *PSdep* = 0*.*83, = 1*.*37) and also as compared to (iii) *top* (*mean rank difference* = 0.64, large effect size *PSdep* = 0*.*72, = 0*.*83). Comparison of (ii) *prefixed* with (iii) *top* results in a non-significant difference (*mean rank difference* = 0.39).

(b) A *Friedman* omnibus test comparing the ITRs between the three experimental conditions (i) *ACL algorithm* **represented as ITRMedian**, (ii) *prefixed* and (iii) *top* shows a significant overall difference with χ2*(*2*)* = 9*.*262, *p* = 0*.*01.

*Post-hoc pairwise comparisons* based upon *critical mean rank differences* 0.82 (Schaich and Hamerle, 1984) vs. 0.57 (Conover, 1980) indicate that ITRs are significantly higher in (i) *ACL algorithm* as compared to (ii) *prefixed* (*mean rank difference* = 0.94, very large effect size *PSdep* = 0*.*81, = 1*.*25) and also as compared to (iii) *top* (*mean rank difference* = 0.64, very large effect size *PSdep* = 0*.*76, = 1*.*21) applying the less conservative criterion of (Conover, 1980). Comparison of (ii) *prefixed* with (iii) *top* results in a non-significant difference (*mean rank difference* = 0.31).

(c) A *Friedman* omnibus test comparing the ITRs between the three experimental conditions (i) *ACL algorithm* **represented as ITRMax**, (ii) *prefixed* and (iii) *top* shows a significant overall difference with χ2*(*2*)* = 22*.*986, *p* = 0*.*00001.

*Post-hoc pairwise comparisons* based upon *critical mean rank differences* 0.82 (Schaich and Hamerle, 1984) vs. 0.41 (Conover, 1980) indicate that ITRs are significantly higher in (i) *ACL algorithm* as compared to (ii) *prefixed* (*mean rank difference* = 1.47, extremely large effect size *PSdep* = 0*.*94,  = 2*.*25) and also as compared to (iii) *top* (*mean rank difference* = 1.19, extremely large effect size *PSdep* = 0*.*94, = 2*.*25). Comparison of (ii) *prefixed* with (iii) *top* results in a non-significant difference (*mean rank difference* = 0.28).

## **THE ACL ALGORITHM SEEMS TO ADAPT TO SUBJECTS' INTERINDIVIDUAL DIFFERENCES**

*NTrials* in condition (i) *ACL algorithm* using two flicker LEDs (see **Table 1**) is deterministically given by 3/4 of the total number of the SSVEP-*SNR* responses under the 20 flicker frequencies in the scanning phase of the experiment which had exceeded the defined threshold value of 10 (*suitable* frequencies), see *ACL Algorithm* of section "Materials and Methods." Thus, in order to make the investigation of possible interindividual associations between the SSVEP-*SNR* magnitudes with the convergence measure second half (see section "Materials and Methods") relatively independent from *NTrials*, all subjects with *NTrials <* 25th percentile (8*.*75 ≈ 9) were excluded, # subject 6, 8, 9, 12, 13, and 16. The resulting rest of *N* = 12 subjects showed a relatively small variability with range of *NTrials* between 11 and 15. The measure SSVEP-*SNR* mean magnitudes in the scanning phase of the experiment (a) over all flicker frequencies from 20 to 39 Hz was split into two measures, one for (b) *lower* frequencies from 20 to 29 Hz and the other for (c) *higher* frequencies from 30 to 39 Hz. In this subsample, convergence measure second half shows large and highly significant correlations with (a) of *r* = 0*.*839, *p* = 0*.*001, with (b) of *r* = 0*.*843, *p* = 0*.*001 and with (c) of *r* = 0*.*763, *p* = 0*.*004. Checking these relationships against the remaining variability of *NTrials* and age as controlled third variables in *partial correlation* analyses, indeed no changes are observed; those found relationships can be considered as linearly independent from *NTrials* and age. Hence, these findings show that the convergence of the *ACL* algorithm highly depends on the subjects' *trait* ability to generate higher SSVEP-*SNR* magnitudes, with no relevant differences observed between *lower* vs. *higher* flicker frequencies: focusing on a subsample with a more or less constant number of *suitable* frequencies, the ACL algorithm explored the more distinct frequencies in those subjects who displayed the *larger* SSVEP-*SNR* magnitudes in the scanning phase of the experiment.

In conclusion, these findings imply that the ACL algorithm shows a distinct exploration behavior for different subjects and thus indeed is able to adapt to subjects' interindividual differences. Whether this adaptation is the *cause* for the ACL algorithm's outperformance of (ii) *top* and (iii) *prefixed* cannot be examined in depth with the employed experimental design and has to be investigated in further studies.

#### **BASELINE RESTING STATE EEG MEASURES AS CORRELATES OF INTERINDIVIDUAL DIFFERENCES**

Searching for significant and relevant associations between interindiviudal variabilities of ITR performances under the three experimental conditions vs. of baseline resting state EEG relative mean PSDs in all computed frequency bands at all eight used electrodes, effects were only found in *thetaHigh* (6.5–7.5 Hz) and *betaMid* (18–24 Hz). In all the other bands nothing could be observed.

Whereas *Pearson* correlations showed no relationships between the resting state relative mean *thetaHigh* PSDs at *Oz* vs. ITRs in conditions (iii) *prefixed* (*r* = 0*.*034, *p* = 0*.*894) and (ii) *top* (*r* = 0*.*196, *p* = 0*.*436), a significant positive correlation with condition (i) *ACL algorithm* was found (*r* = 0*.*467, *p* = 0*.*048) representing the performance as ITRMedian. Searching for similar relationships in the other seven used electrodes, no associations were observed; these effects exclusively occur at *Oz* in our sample*.* Following the effects size classifications of Cohen (1988), this correlation is to be considered as *moderate. Partial correlation* analyses confirmed that this correlation is linearly independent against age and all means of SSVEP-*SNRs* in the previous scanning phase of the experiment over (a) *all* 20 flicker frequencies, (b) also over the *lower* frequencies 20–29 Hz and (c) also over the *higher* frequencies 30–39 Hz.

At least in the examined sample, interindividual variability in relative mean *thetaHigh PSD* at *Oz* seems to differentiate between *ACL algorithm* and the other two conditions: the larger the observed relative mean PSDs among subjects in the baseline resting state are, the better will be their later SSVEP-BCI performance exclusively under the use of *ACL algorithm.*

At first sight, analyzing baseline resting state relative mean *betaMid* PSDs, an exclusive relationship with only the ITRs in condition (iii) *top* was found for *F3* (*r* = 0*.*484, *p* = 0*.*042), although its neighbor electrodes also showed relationships not very far away from significance, probably due to small sample

*distances* **in brackets calculated in a linear regression analysis with the**

size: *F4* with *r* = 0*.*425, *p* = 0*.*117 and *Fz* with *r* = 0*.*410, *p* = 0*.*091. All the other used electrodes showed no associations. After further graphic inspection of relevant scatterplots and *Box-Whisker-Plots,* a possible negative relationship between baseline resting state relative mean *betaMid* PSDs at *Oz* and ITRMean in condition (i) *ACL algorithm* was suspected, hidden by outliers. *Box-Whisker-Plots* suggested case 15 and 11 as outliers, so for further analysis *Mahalanobis distances* were computed in a linear regression analysis with the ITRsMean of condition (i) *ACL algorithm* as criterion variable and baseline resting state relative mean *betaMid* PSDs at *Oz* as predictor variable. The inspection of *Mahalanobis distances* and the scatterplot (see **Figure 7**) suggest that subject 15 and 11 might be considered as outliers. Excluding them changes the correlation from *r* = −0*.*262, *p* = 0*.*294 to significant *r* = −0*.*530, *p* = 0*.*042. *Partial correlation* analyses confirmed that this correlation is linearly independent against age and all means of SSVEP-*SNRs* in the previous scanning phase of the experiment (a), (b), and (c) mentioned above.

Interestingly, excluding case 15 and 11, baseline resting state relative mean PSDs *betaMid* vs. *thetaHigh* both at *Oz* show an almost significant correlation over the subjects with *r* = −0*.*482 and *p* = 0*.*059, probably due to the small sample size, which is stable against the third variables age and all SSVEP-*SNRs* in the previous scanning phase of the experiment (a), (b), and (c), mentioned above.

*r* = −0*.*262, *p* = 0*.*294 to significant *r* = −0*.*510, *p* = 0*.*043.

In conclusion, baseline resting state relative mean *betaMid* PSDs seem to predict ITR performances under (i) *ACL algorithm* vs. (iii) *top* in an opposed fashion depending on the electrodes: the *lower* baseline resting state relative mean *betaMid* PSDs are at *Oz*, the *higher* will be the ITRs under condition (i); and the *higher* baseline resting state relative mean *betaMid* PSDs are at frontal electrodes (*F3, Fz, F4)* the *higher* will be the ITRs under condition (iii). In addition to these findings in *betaMid*, the higher the baseline resting state relative mean *thetaHigh* PSDs at *Oz* are, the higher will be the ITRs exclusively under condition (i).

Returning to the above described subsample of *N* = 12 obtained by exclusion of all subjects with *NTrials <* 25th percentile (8*.*75 ≈ 9), an interesting observation was found: IAF shows differentiating relationships with ITR performances: a significant correlation of *r* = 0*.*577, *p* = 0*.*0496 was only found with ITRs under (i) *ACL algorithm* (see scatterplot **Figure 8**), but neither under (ii) *top* with *r* = 0*.*394, *p* = 0*.*205 nor under (iii) *prefixed r* = 0*.*283, *p* = 0*.*373. The higher subjects' IAF are in the subsample, the better will be their ITR performance exclusively under the ACL algorithm. *Partial correlation* analyses confirmed that this association is linearly independent against age. Repeating this analysis for the entire sample of *N* = 18 *no* significant correlations between individual alpha frequency (IAF) and ITR performances under the three experimental conditions become apparent (i) with *r* = 0*.*282, *p* = 0*.*257, (ii) *r* = 0*.*198, *p* = 0*.*432 and (iii) *r* = 0*.*243, *p* = 0*.*332. These findings imply that subjects with *low* ITRs in all three conditions might represent another population as compared to the rest. Further studies may try to replicate these findings and identify dimensions which discriminate between these possible two different populations. Moreover, these findings could be relevant for the understanding

**as continuous line, 95% confidence regression bands as dotted lines).** A significant *Pearson* correlation with *r* = 0*.*577, *p* = 0*.*0496 was found in the remaining subsample of *N* = 12 (blue points), removing subjects with *NTrials <* 25th percentile (8*.*75 ≈ 9) (red points), while over the entire sample of *N* = 18 the correlation is hidden with *r* = 0*.*282, *p* = 0*.*257 (all points). This relationship seems to exist exclusively for condition (i) *ACL algorithm*: the higher subjects' IAF are in this subsample, the better will be their ITRMean performance exclusively under (i). *Partial correlation* analyses confirmed that this association is linearly independent against age.

of the so-called *BCI illiteracy* phenomenon (Blankertz et al., 2010; Vidaurre and Blankertz, 2010; Volosyak et al., 2011), see section "Discussion."

Inspired by the findings of Koch et al. (2008) who found correlations of IAF with both magnitudes of *visually evoked potential*s (VEPs) and also with cortical oxygenation measured by *nearinfrared spectroscopy* (NIRS), *Spearman* rank order correlations were computed between IAF and means of SSVEP-*SNR* magnitudes in the scanning phase of the experiment (a) over *all* 20 used flicker frequencies 20–39 Hz, (b) over the *lower* frequencies 20–29 Hz and (c) over the *higher* frequencies 30–39 Hz in the described subsample of *N* = 12. Although not fully reaching significance level, probably due to the relatively small sample size, an interesting pattern was found: IAF vs. (a) with *rho* = 0*.*561, *p* = 0*.*058, IAF vs. (b) with *rho* = 0*.*183, *p* = 0*.*568 and IAF vs. (c) *rho* = 0*.*557, *p* = 0*.*060. Although not fully significant, probably due to the small sample size, interindividual differences in SSVEP-*SNR* magnitudes under the employed *higher* flicker frequencies seem to show a tendency of positive association to higher IAFs while this relationship might not exist for the stimulation with the *lower* frequencies (or if so, it may presumably be lower). These findings motivated the re-analysis of the found relationship in **Figure 8** by *partial correlations* whether it would be linearly independent against SSVEP-*SNR* magnitudes in the scanning phase of the experiment (a), (b) and (c) as described above. While (a) and (b) showed no relevant influence on this relationship, controlling for (c) resulted in a reduction from former *Pearson r* = 0*.*577, *p* = 0*.*0496 to *r* = 0*.*396, *p* = 0*.*228. Hence, these findings imply that IAF and (c) the magnitude of SSVEP responses to only the employed *higher* flicker frequencies share remarkably amounts of common interindividual variability while explaining variability of ITRMean under the ACL algorithm.

#### **EFFECTS OF THE PERMUTATION OF PRESENTATION ORDER**

Investigating possible effects of the permutation of presentation order, a *mixed-design repeated measures ANOVA* was computed with *stimulation condition* as repeated *within-*subjects factor with three levels (i) ACL algorithm represented as ITRMedian, (ii) prefixed and (iii) top and *presentation order* as *between-*subjects factor with the six possible permutations as levels (ABC, ACB, BAC etc.). *Levene*'*s* tests showed homogeneities of error variances. There was no significant *between*-subjects main effect of *presentation order* with *F(*5*,* 18*)* = 2*.*26, *p* = 0*.*115, η2 p = 0*.*485. Because *Mauchly*'*s sphericity* test indicated a significant departure from the assumption of *sphericity* with χ2*(*2*)* = 6*.*54, *p* = 0*.*038, *Greenhouse-Geisser* estimates were used to correct degrees of freedom (ε = 0*.*691). There was no significant *interaction* between *presentation order* × *stimulation condition* with *F(*10*,* 18*)* = 0*.*67, *p* = 0*.*738, η2 p = 0*.*219. *ANOVA* analyses were repeated also for condition (i) ACL algorithm represented as ITRMean and ITRMax which resulted in similar findings. In conclusion, neither significant main effects nor significant interactions could be found over all six possible permutations of presentation order. Hence, the found effects in the *central tendencies* reported above with regard to all ITR performances can be considered as independent from possible presentation order effects.

## **DISCUSSION**

Although electrophysiology-based closed-loop interactions with biological nervous systems have been used since the 1940s, modern computers and online software control techniques allow a wide variety of novel activity dependent protocols in neuroscience research and related applications. Current BCI bring up a number of problems related to relatively long previous training times and still relatively low efficiencies (ITRs). This calls for novel techniques which can also address context and subject specificities, e.g., adaptive detection of SSVEPs (e.g., Krauledat et al., 2008).

In this paper we described an *assisted closed-loop protocol* which enhances BCI efficiency, as compared to classic BCI protocols, by providing both the subject and the system with online information which helps them to reach the BCI goal in their interaction. We used a reductionistic paradigm to constrain the inherent complexity of closed-loop exploration: four simultaneous frequencies, a basic SSVEP detection strategy and a relatively simple task to be accomplished by the user. More complex BCI systems might further benefit from the described approach. Our paradigm calls for many possible improvements, ranging from advanced SSVEP detection algorithms, stimuli which inform the user more effectively, up to a more adaptive online control of the interface itself by measuring and exploring additional dimensions (multimodality).

The literature on SSVEP-BCIs does not report general recommendations for the selection of the properties of the visual stimuli (Wu et al., 2008; Zhu et al., 2010), although it is known that the SSVEP magnitudes depend on extrinsic and intrinsic dimensions (Ding et al., 2006; Lopez-Gordo et al., 2011). Our study shows that a closed-loop subject-specific selection of the stimulation frequencies together with the closed-loop auditory feedback lead to increased BCI ITR performance which outperformed the employed control conditions.

Although *assisted closed-loop protocols* seem to enhance BCI efficiency, their use is limited by the additional time needed for the exploration process. In the protocol discussed in this paper, the average time to perform the experiment was around half an hour, flicker frequency selection took most of this time. Due to time restrictions, the parameter space can never be explored completely, so BCI efficiency improvement might remain suboptimal. Thus, there is some unknown *trade-off* between improvement and time needed, which should be explored in further studies. Furthermore, the question how replicable the found flicker frequencies are in the same subjects over multiple follow-up time points could be explored. Probably, observing this stability over time (e.g., *test-retest reliability*) may help to discover important trait vs. state dimensions related to variability of BCI performance. Another limitation due to the SSVEP physiology is that the time window for the auditory feedback is relatively short, so subjects have to establish control over the BCI goal in the range of a few seconds. This implies possible interactions with subjects' traits and states related to cognitive *processing speed* and dimensions of learning abilities.

ACL algorithms offer new possibilities as compared to traditional open-loop paradigms, but require additional decisions and new perspectives for their design and analysis, e.g., with regard to online measurement of actual states and performance, parameter search responding to the particular dynamic behavior of the system, properties of the feedback stimuli, *actuation laws*, etc. However, our findings imply that this additional effort can improve BCI efficiency and contribute to reveal dynamics of the nervous system which would remain hidden under traditional paradigms. Because our analyses showed that EEG resting state measures can predict assisted closed-loop SSVEP-BCI performance, our novel approach seems to flexibly adapt/interact with interindividual cerebral variabilities. Although found in the context of a *sensory motor rhythms* (SMRs) based BCI, other recent work also demonstrated that EEG resting state measures can be relevant predictors of BCI performance (Blankertz et al., 2010). In this emerging field, it could be fruitful to identify possible EEG resting state measures which can differentiate/predict between BCI performances based on biosignals originating from distinct physiological mechanisms: SSVEPs, P300*,* SMRs, *slow cortical potentials* (SCPs)*, electrocorticogram* (ECoG), *magnetoencephalography* (MEG), NIRS *or blood-oxygen-level-dependent* (BOLD). Apart from these biosignals reflecting *brain* activity, *peripheral* psychophysiological measures have been investigated in the context of BCIs, especially as performance predictors, such as *parasympathic/vagal* parameters of resting state *heart rate variability* (HRV) (Kaufmann et al., 2011).

Our proposed approach of new adaptive-interactive paradigms might offer innovative ways how to address the problem of the so-called *BCI illiteracy,* i.e., the incapacity of some subjects to achieve control of BCIs (Blankertz et al., 2010; Vidaurre and Blankertz, 2010; Volosyak et al., 2011). It might be fruitful to explore the possible different impact of ACL algorithms in BCIs based on the mentioned distinct physiological mechanisms, especially with regard to their specific BCI illiteracies.

As mentioned in section "Baseline Resting State EEG Measures as Possible Interindividual Correlates of ITR Performances," the IAF is a measure of *interindividual* EEG variability because it is remarkably stable *within* subjects, but relatively variable *between* subjects (Kondacs and Szabó, 1999). IAF seems to be highly heritable, e.g., Posthuma et al. (2001) found in a study comparing mono- vs. dizygotic twins, analyzing a large representative sample of healthy Dutch adults (*N* = 688), that 71–83% of total IAF variance could be ascribed to genetic variances. Thus, IAF may be considered as an *endophenotype* following the definition of Gottesman and Gould (2003). Klimesch (1997) found in a sample of age matched subjects that the IAF of good working memory performers is about 1 Hz higher vs. that of bad performers. Jin et al. (2006) found that IAF is positively correlated with conflict reaction time. Severity of *Alzheimer*'*s* disease is positively related to the extent of typical IAF slowing in this pathology (Rodriguez et al., 1999). On the neurophysiological level, Steriade et al. (1990) reported that IAF depends on membrane properties of the thalamic neurons which project to the cortex, implying thalamo-cortical feedback loops as one of the important generators of alpha activity (Lopes da Silva, 1991). Mayer et al. (2007) successfully modeled the synchronization of locally coupled bistable thalamic oscillators as controlled by the influence of corticothalamic projections, probably responsible for widespread spindle oscillations in the thalamus. Given these findings, IAF might be understood as a positive correlate of thalamo-cortical information processing speed. With regard of possible correlations of IAF with SSVEP magnitudes, Koch et al. (2008) found interesting correlations of IAF with both magnitudes of VEPs and cortical oxygenation measured by NIRS. Concluding, IAF seems to open new insights into the understanding of the neural circuits underlying BCI performance and thus should be considered as a promising predictor for further studies.

In this study, only eight EEG electrodes were used to investigate EEG resting state measures as performance predictors, but further works might use more electrodes of the 10–20 system to allow a-posteriori offline analyses of *scalp maps* and the use of *source localization* techniques, e.g., *LORETA* (for a review see Grech et al., 2008). Findings of research concerning the cerebral *resting-state networks* call for further studies which use simultaneous EEG/fMRI recordings (for reviews see e.g., Fox and Raichle,

## **REFERENCES**

*Electrochemistry in Biology and Medicine,* ed T. Shedlovsky (New York, NY: Wiley), 121–140.

2007; Van den Heuvel and Hulshoff-Pol, 2010; for typical studies see e.g., Damoiseaux et al., 2006; Van den Heuvel et al., 2009; Yuan et al., 2012).

Opening the scope to other uses, the demonstrated advantage of our adaptive-interactive BCI protocol can be expanded conceptually, e.g., to innovative applications such as diagnostic/therapeutic tools in clinical contexts: exploring the subjectspecific dynamical trajectory of machine-subject interaction could extract information which otherwise would remain undiscovered. Thus, far beyond an engineering focus, the proposed approach might be employed as a new paradigm for basic neuroscientific and biomedical research.

## **ACKNOWLEDGMENTS**

We thank Víctor Bonilla for technical help. This work was supported by UAM CEMU 2012-004, MINECO TIN2012-30883 and TIN-2010-19607.

Federation. *Electroencephalogr. Clin. Neurophysiol.* 10, 370–375.

performance of independent binary SSVEP-BCIs. *Clin. Neurophysiol.* 122, 128–133.

of the brain at rest–exploring EEG microstates as electrophysiological signatures of BOLD resting state networks. *Neuroimage* 60, 2062–2072.

**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.

*Received: 31 October 2012; paper pending published: 22 December 2012; accepted: 06 February 2013; published online: 25 February 2013. Citation: Fernandez-Vargas J, Pfaff HU, Rodríguez FB and Varona P (2013) Assisted closed-loop optimization of SSVEP-BCI efficiency. Front. Neural Circuits 7:27. doi: 10.3389/fncir. 2013.00027*

*Copyright © 2013 Fernandez-Vargas, Pfaff, Rodríguez and Varona. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any thirdparty graphics etc.*

## Coupling BCI and cortical stimulation for brain-state-dependent stimulation: methods for spectral estimation in the presence of stimulation after-effects

#### *Armin Walter 1 \*, Ander R. Murguialday2,3, Wolfgang Rosenstiel 1, Niels Birbaumer 2,4 and Martin Bogdan1,5*

*1 Department of Computer Engineering, Wilhelm-Schickard-Institute, Eberhard Karls Universität Tübingen, Tübingen, Germany*

*2 Institute of Medical Psychology and Behavioural Neurobiology, University Hospital Tübingen, Tübingen, Germany*

*3 Health Technologies Department, TECNALIA, San Sebastian, Spain*

*4 Ospedale San Camillo, IRCCS, Venice, Italy*

*5 Department of Computer Engineering, University of Leipzig, Leipzig, Germany*

#### *Edited by:*

*Ahmed El Hady, Max Planck Institute for Dynamics and Self Organization, Germany*

#### *Reviewed by:*

*Yang Dan, University of California, USA*

*Pratik Y. Chhatbar, Medical University of South Carolina, USA*

#### *\*Correspondence:*

*Armin Walter, Department of Computer Engineering, Wilhelm-Schickard-Institute, Eberhard Karls Universität Tübingen, Sand 14, 72076 Tübingen, Germany. e-mail: armin.walter@ uni-tuebingen.de*

Brain-state-dependent stimulation (BSDS) combines brain-computer interfaces (BCIs) and cortical stimulation into one paradigm that allows the online decoding for example of movement intention from brain signals while simultaneously applying stimulation. If the BCI decoding is performed by spectral features, stimulation after-effects such as artefacts and evoked activity present a challenge for a successful implementation of BSDS because they can impair the detection of targeted brain states. Therefore, efficient and robust methods are needed to minimize the influence of the stimulation-induced effects on spectral estimation without violating the real-time constraints of the BCI. In this work, we compared four methods for spectral estimation with autoregressive (AR) models in the presence of pulsed cortical stimulation. Using combined EEG-TMS (electroencephalography-transcranial magnetic stimulation) as well as combined electrocorticography (ECoG) and epidural electrical stimulation, three patients performed a motor task using a sensorimotor-rhythm BCI. Three stimulation paradigms were varied between sessions: (1) no stimulation, (2) single stimulation pulses applied independently (open-loop), or (3) coupled to the BCI output (closed-loop) such that stimulation was given only while an intention to move was detected using neural data. We found that removing the stimulation after-effects by linear interpolation can introduce a bias in the estimation of the spectral power of the sensorimotor rhythm, leading to an overestimation of decoding performance in the closed-loop setting. We propose the use of the Burg algorithm for segmented data to deal with stimulation after-effects. This work shows that the combination of BCIs controlled with spectral features and cortical stimulation in a closed-loop fashion is possible when the influence of stimulation after-effects on spectral estimation is minimized.

**Keywords: brain-computer interfaces, cortical stimulation, spectral estimation, brain-state-dependent stimulation, autoregressive models**

#### **1. INTRODUCTION**

Cortical stimulation is being used to study cortical function, e.g., (Matsumoto et al., 2007). In clinical settings, it is employed for surgical planning (Lefaucheur and de Andrade, 2009) and therapy (Tsubokawa et al., 1991). Furthermore, preliminary studies on the use of cortical stimulation for stroke rehabilitation which used stimulation together with physiotherapy in order to modulate cortical excitability have been conducted (Brown et al., 2008; Levy et al., 2008). Taking the current brain activity of the patient into account when selecting stimulation parameters has been proposed as a possible improvement (Plow et al., 2009). Such an activity-dependent stimulation paradigm has been used by Jackson et al. (2006), who were able to show that cortical microstimulation associated in time with brain activity during a motor task can induce neural reorganization lasting for several days after stimulation in primates.

The effects of transcranial magnetic stimulation (TMS) as well depend on brain states of the stimulated person (Mitchell et al., 2007). Recently, Bergmann et al. (2012) applied TMS coupled to electroencephalography (EEG) to investigate the dependency of stimulation effects on the phase of slow EEG oscillations during sleep. In general, such activity-dependent or brain-state-dependent stimulation (BSDS) paradigms allow to investigate cortical networks at specific activation levels, making BSDS a potentially useful tool in cognitive neuroscience (Jensen et al., 2011) as well as in clinical studies improving consistency of the stimulation effects (Plow et al., 2009).

For effective BSDS, reliable decoding of the brain-state from the ongoing brain activity is necessary. Over the last decades in the field of brain-computer interfaces (BCIs) several different strategies were investigated (Birbaumer et al., 1999; Birbaumer and Cohen, 2007). Especially in the case of movement-related brain states during active or imagined limb movements, spectral power has been shown to be useful for their decoding. In particular, event-related (de-) synchronization of sensorimotor rhythms is an informative measure for discriminating movement and resting states (Wolpaw et al., 2002). Therefore, if one wants to combine BSDS with a movement task, one has to minimize the interference of the stimulation on the estimation of the spectral features to detect the brain-state properly.

The stimulation effects involve problems with spectral estimation caused by the stimulation artefact and the evoked neural activity. A stimulation pulse evokes an artefact in the signal (**Figure 1A**) with an amplitude in the range of several hundred millivolts or even volts, thus often exceeding the dynamic range of the amplifier (Veniero et al., 2009). In the vicinity of stimulation, evoked potentials are recorded (**Figure 1B**) which can reach amplitudes of several hundred microvolts (Matsumoto et al., 2007). Thus, if an analyzed window contains a stimulation pulse, the estimation of the spectrum of this window is difficult, because it is not stationary. This is evident in **Figure 1C**, showing that each stimulation pulse results in strong jumps in the estimated spectral power. Waiting long enough after the pulse is one solution. This approach results in non-continuous brain-state decoding with waiting periods after a stimulus of at least several hundred milliseconds. It dictates a longer inter-stimulus interval (ISI), because a robust

components are covered. **(C)** Time-course of the power at 12 Hz of the signal displayed in **(A)**, resulting from a time-frequency analysis with auto-regressive models (order 16) when a window of 500 ms is shifted in 40 ms steps over the data. Hence, a single stimulation pulse distorts the spectrum for the next 500 ms because it remains in the data window. **(D)** Zoom on the region of the last stimulation pulse. Power at 12 Hz without stimulus processing (solid line) and when the gap is defined as in **(B)** and either MEMgap (dashed line) or AR modeling with order 16 (dashed-and-dotted line) are applied to deal with it.

estimate of the brain-state is needed before the next pulse can be applied.

If small ISIs and/or continuous decoding of the brain-state is necessary, methods that enable spectral estimation of data containing stimulation after-effects are mandatory. One potential solution for this, which has been used mainly in offline studies (no BSDS), is to separate the stimulation effects from the signal, as for example in Litvak et al. (2007). This places restrictions on the recording setup, such as the need for an amplifier with high dynamic range to cover the entire amplitude of the artefact and it is unclear whether such a procedure can be performed online without resulting in residual artefacts which would still lead to distortions of the spectrum. We present in this paper another solution suitable for online BSDS: we ignore the short segment of data dominated by the after-effects of stimulation when estimating the spectrum, leaving us with the challenge to estimate the spectrum when portions of the data are missing from a continuous data flow. We term such an excluded data segment a *gap*. In online experiments, using either signals synchronized with the stimulator or a peak detection algorithm, one can mark a sample before the stimulation pulse as the beginning of the gap. The number of following samples marked as belonging to the gap (i.e., the *gap size*) should be chosen in advance such that the gap, ideally, encloses just the stimulation artefact, and the largest evoked components. The dashed line and the dashed-and-dotted line in **Figure 1D** show the results of two approaches introduced in this work to extract the spectral power when the artefacts are masked by the gap shown in **Figure 1B**. They are much closer to the power before and after the stimulus, compared to the power without any processing of the stimulus (solid line).

In this paper we compare different online brain-state decoding methods on their suitability to perform spectral estimation using autoregressive (AR) models on data containing stimulation pulses and gaps. We consider here only stimulation paradigms with pulsed stimuli and restrict ourselves to data acquired with EEG or electrocorticography (ECoG) and stimulation performed using TMS or epidural electrodes. First, we introduce the methods for spectral estimation in the presence of gaps and investigate

the effects of parameter estimates such as AR model order and gap size on the resulting spectrum. We present results from a simulation study in which gaps are artificially inserted into a BCI data set recorded without stimulation. We then show the different results of the algorithms on short data segments of two BCI training experiments, one with simultaneous TMS and one with simultaneous epidural electrical stimulation to illustrate the effects of cortical stimulation on spectral estimation and the results of correcting stimulation after-effects. Finally, we investigate the separability of intended hand movement and rest for different experimental paradigms (no stimulation, open-loop, or closed-loop stimulation) using non-invasive and invasive data during BCI experiments in three chronic stroke patients.

## **2. METHODS**

#### **2.1. PARTICIPANTS**

Data was recorded from three chronic stroke patients (**Table 1**) suffering from paresis of the left hand. None of the patients was able to produce voluntary finger movements with the left hand. All procedures were approved by the local ethics committee of the medical faculty of the university hospital in Tübingen. Each stroke patient was implanted with 16 epidural platinum iridium disk electrodes (Resume II, Medtronic, Fridley, USA) with a contact diameter of 4 mm placed over the ipsilesional S1, M1, and pre-motor cortex on four strips with an inter-electrode center-tocenter distance of 10 mm. They were arranged in a 4 × 4 grid-like pattern (**Figure 2**). During pre-surgical evaluation, all subjects completed the task described below with combined EEG-TMS (*non-invasive case*) and repeated the same task after the surgery using electrical epidural stimulation and recordings from the implanted electrodes (*invasive case*). The BCI and stimulation experiments were conducted during a period of 4 weeks following the implantation.

#### **2.2. TASK**

The patient was facing a 19- monitor. The left upper limb of the patient was fixed using two straps, one at the forearm and one around the wrist and magnets fixed the fingertips to the actuators

**post-surgical CT for the three patients.** From left to right: P1–P3.

of a mechatronic hand orthosis (Tyromotion Amadeo HTS, Graz, Austria). This device was controlled by a BCI and moved the fingers of the paralyzed hand between an opened and a closed position. The range of the movement was adjusted in each session (Ramos-Murguialday et al., 2012) because it was limited by the spasticity of the patient. Each trial of the task consisted of three phases: preparation (2 s), feedback (6 s), and rest (8 s). During preparation, the subject received an auditory cue ("Left Hand") but was instructed to wait with the execution until the next auditory command ("Go!") was given at the start of the feedback phase. During the feedback phase starting with a closed position of the left hand, the patient had to try to open the left hand until the end of the feedback phase. At that point, another auditory cue ("Relax!") was given. During the rest period, the left hand of the patient was returned to its original closed position (2–3 s) and the patient was instructed to relax. An experimental session was divided into a 4–16 runs, each of these consisting of 11 trials. Runs with clear non-stimulation-related artefacts (e.g., amplifier saturation) on the analyzed channels were excluded from further analysis, resulting in a minimum of three runs per session for analysis and an average of 8*.*7 ± 4*.*3.

#### **2.3. ELECTROPHYSIOLOGICAL RECORDING**

Both EEG and ECoG were recorded with monopolar 32 channel amplifiers (BrainAmp MR plus, BrainProducts, Munich, Germany) with a sampling rate of 1000 Hz. The data was acquired in a packet-wise fashion, where the recording computer received every 40 ms one packet of data consisting of 40 samples per channel. The same behavior was modeled in our simulations of an online BCI. A high-pass filter with a cutoff frequency at 0.16 Hz and a low-pass filter with a cutoff frequency at 1000 Hz were applied. We recorded 32 channels of EEG in the standard 10–10 system, referenced to FCz, using circular Ag-AgCl electrodes. ECoG data was referenced to an electrode at the medio-frontal corner of the electrode grid over pre-motor cortex. Signal acquisition, signal processing and control of the orthosis and (if present in the experiment) the TMS or electrical stimulator were performed using the general-purpose BCI framework BCI2000 (http://www*.*bci2000*.*org) (Schalk et al., 2004) extended with custom-developed features for the control of these devices.

#### **2.4. STIMULATION**

We applied stimulation in the non-invasive case over the hotspot for *extensor digitorum communis* (EDC) activity, identified by a standard mapping paradigm (Wassermann et al., 2008). TMS was applied with a figure-of-eight coil (NeXstim, Helsinki, Finland) with single biphasic pulses (sinusoidal coil current, positive phase first, pulse width 280μs) and an intensity of 110% of the resting motor threshold. The ISI of successive pulses was set to 3 s.

For epidural electrical stimulation we used single biphasic anodal square-wave pulses with a length of 500μs. Stimulation intensity was selected individually per patient and session and chosen to reliably evoke MEPs on the paretic upper limb of the patient. The minimum ISI was set to 2 s in most experiments except when stimulation was applied coupled with the BCI output. In this case, a minimum ISI of 500 ms was chosen. The pulses were applied using a constant current stimulator (STG4008, Multichannel systems, Reutlingen, Germany) with the anode as the epidural electrode that evoked the strongest MEPs on the left upper limb and the cathode being a 50 × 90 mm adhesive electrode placed on the left clavicle of the patient. The current source of the stimulator was switched off 2 s after the last stimulation pulse if no other pulse was triggered before due to a software error, leading to a small but visible step in the recorded signal (**Figure 1B**).

#### **2.5. AUTOREGRESSIVE (AR) MODELS**

A popular choice for spectral estimation in BCI research is to use an AR model for which the coefficients are estimated with the maximum entropy method (Krusienski et al., 2006; McFarland and Wolpaw, 2008). An AR model can be viewed as a linear predictor of the signal samples *x(tk)*, defined as:

$$\mathfrak{x}(t\_k) = \sum\_{i=1}^{p} c\_i \mathfrak{x}(t\_{k-i}) + e$$

where *p* is the order of the model and *e* a sample of a white noise process. If one uses a continuous window of length *N* with *N p* consisting of samples *x(t*0*)* to *x(tN*−1*)*, one could solve the following equations with a least-squares procedure to get the coefficients *ci*:

$$\mathbf{x}(t\_k) = \sum\_{i=1}^{p} c\_i \mathbf{x}(t\_{k-i}) \tag{1}$$

$$\mathbf{x}(t\_{k-p}) = \sum\_{i=1}^{p} c\_i \mathbf{x}(t\_{k-p+i}) \text{ for all } k = p, \dots, N-1$$

However, the resulting coefficients do not guarantee a stable AR model (de Waele and Broersen, 2000). Burg proposed a recursive algorithm for the solution of this system that provides stable models with less variance compared to least squares solutions and the Yule-Walker algorithm (Kay and Marple, 1981; de Waele and Broersen, 2000). The Burg algorithm computes the AR coefficients in *p* steps by evaluating in the *i*-th step the residuals of forward and backward prediction of the samples using the coefficients obtained in the *(i* − 1*)*-th step. It is described in Appendix 1.1 in more detail. Spectral estimation with AR models is briefly introduced in Appendix 1.2.

The Burg algorithm requires that the input data is sampled continuously without gaps, a condition which is shared by most of the other algorithms for AR model estimation. Therefore, we need to either fill or remove the gaps before applying one of these algorithms to our data or modify the AR model estimation algorithms to be usable for data with gaps.

#### **2.6. SPECTRAL ESTIMATION IN THE PRESENCE OF GAPS**

This section contains a short description of the different algorithms we compare in this paper that deal with the pre-processing of data containing gaps for spectral estimation with AR models. The input for these algorithms are a segment of data and a vector that contains for each sample in the segment either a 1 (sample belongs to a gap, it has to be excluded from spectral estimation) or a 0 (sample is "clean").

Four methods for dealing with gaps in the data are described below: (1) *linear interpolation*, (2) *AR modeling* which fill the gap with generated data, (3) the *joining of data segments* that removes the gap, and (4) a *modified Burg algorithm for segmented data*. After application of the methods (1)–(3), the standard Burg algorithm is used to estimate the AR model and the spectrum.

#### *2.6.1. Linear interpolation*

We can bridge gaps in the data by linear interpolation between the last sample before and the first sample after the gap:

$$\hat{\mathbf{x}}(t\_{\mathcal{S}+k}) = \mathbf{x}(t\_{\mathcal{S}-1}) + \frac{k+1}{l+1} \cdot \left(\mathbf{x}(t\_{\mathcal{S}+l}) - \mathbf{x}(t\_{\mathcal{S}-1})\right), \ 0 \le k \le l-1 \tag{2}$$

where *x* are the signal samples recorded at times *ti*, *l* is the length of the gap in samples and *tg*−1 is the index of the last sample before the gap.

While this might work for offline analysis of a data set, in the case of online analysis during a BCI experiment, in which data is received in a sample- or packet-wise system, one might have not yet received the first clean sample after the gap when trying to produce an estimate for *x(tg* + *k)* within the gap. We used a simple approach to solve this problem which consists of filling the gap with the value of the last sample before the gap (*x*ˆ*(tg* + *k)* = *x(tg* − 1*)*) as long as we have not received the packet containing the end of the gap and using linear interpolation for the rest of the gap otherwise. We term this approach *on-line compatible linear interpolation*.

#### *2.6.2. AR modeling*

As a somewhat more sophisticated technique compared to linear interpolation, we generated data from an AR model to fill the gap. For this we used the coefficients *ci* of the AR model estimated for the data window directly before the gap to predict the missing samples *x*ˆ:

$$\hat{\mathbf{x}}(t\_{\xi+k}) = \sum\_{i=1}^{p} c\_i \mathbf{x}'(t\_{\xi+k-i}) + \sigma \cdot e(t\_{\xi+k}), \ 0 \le k \le l-1,\tag{3}$$

$$\mathbf{x}'(t\_{\xi+j}) = \begin{cases} \mathbf{x}(t\_{\xi+j}) & \text{if } j < 0\\ \hat{\mathbf{x}}(t\_{\xi+j}) & \text{otherwise} \end{cases}$$

*x* can refer to either actually recorded samples before the gap or estimated samples by the AR procedure. σ is the standard deviation of the white noise component in the estimated AR model and *e(t)* one value of a white noise process. While this approach has the property to generate data for the gap consistent with the previously measured data, one might prefer to use a mixture of AR modeling and linear interpolation for the online case. This would avoid jumps in the data when merging generated data within the gap with new samples acquired after the gap. These jumps occur for all AR model orders we have tested in our simulations (see Appendix 1.4 for details). We have used this combination here by performing AR extrapolation when information about the first sample after the gap was not available and using linear interpolation otherwise. The signal was received in packets with a length of 40 ms and for each packet, one of three actions were taken: (1) if a packet contained the start and the end of a gap, then linear interpolation was used to fill the gap. (2) If it contained only the start or if the whole packet was part of the gap, then the AR model was used as a linear predictor to fill the gap. (3) If it contained only the end of the gap, then the last sample of the last packet and the first sample after the gap were connected by linear interpolation.

#### *2.6.3. Joining two segments*

If one chooses to ignore the information of the gap altogether when estimating the model, one might consider simply joining the two segments around the gap, therefore sacrificing information about the timing in the vicinity of the gap. In practice, this means that we update the data window only with those samples from a newly acquired data packet that do not belong to a gap. In order to keep the window size for spectral estimation constant, this has the consequence that older samples are used to compute the spectrum with this method compared to the other algorithms.

## *2.6.4. Burg algorithm for segments (MEMgap)*

For standard algorithms that compute the AR coefficients, the samples within the data window need to be continuous. We can make the least-squares estimation of the AR coefficients compatible with data containing gaps by eliminating all equations from (Equation 1) that contain samples from within a gap and then solving the rest of the equations for the coefficients *ci*. As the Burg algorithm (see Appendix 1.1) yields more stable AR models than the least-squares estimation, we modified it to work with gaps based on the Burg algorithm for segmented data proposed in de Waele and Broersen (2000). This was achieved by limiting the computation of forward and backward prediction errors in each step of the algorithm to those samples that are far enough away from a gap. In the remainder of this paper, this algorithm is called MEMgap (Maximum Entropy Method for data with gaps) for brevity. A detailed description of the algorithm is given in Appendix 1.3.

#### **2.7. SIMULATIONS ON CLEAN DATA**

To study empirically the influence of gaps on the estimated spectrum, we performed simulations on 12 data sets that were recorded without stimulation by artificially inserting gaps, then applying the methods described above to estimate the spectrum. The results of the different methods were compared with a reference time-frequency analysis obtained when using the original data set without gaps. Each data set has a length of 182 s. These data sets, each containing 11 trials, were recorded with ECoG in patient P1 in one experimental session. For clarity reasons, we restrict ourselves to one channel (an electrode over right M1). For spectral computation we kept the length of the window constant at 500 ms and the update rate at 25 Hz = 40 samples. We estimated the power at frequencies between 5 and 99 Hz in 2 Hz increments and varied for each method the gap size (0–100 ms in steps of 5 ms) and the model order (values: 16, 32, and 64). We computed the normalized bias, root mean squared error (RMSE) and variance (var) of the stimulus processing algorithms as follows:

$$\text{bias}(f) = \frac{\frac{1}{n} \sum\_{i} (P(f, i) - P\_0(f, i))}{\overline{P\_0(f)}}$$

$$\text{RMSE}(f) = \frac{\sqrt{\frac{1}{n} \sum\_{i} (P(f, i) - P\_0(f, i))^2}}{\overline{P\_0(f)}}$$

$$\text{var}(f) = \text{Var}\left(\frac{P(f) - P\_0(f)}{\overline{P\_0(f)}}\right)$$

*P(f,i)* is the spectral power of data window *i* for frequency bin *f* , *P*0*(f, i)* is the power of the original data window without gaps and *P*0*(f)* is the average power of the full original recording without gaps for frequency bin *f* . *n* is the number of data windows that are affected by gaps (i.e., data windows where *P(f,i)* − *P*0*(f,i)* is not zero). var*(f)* is the variance of the difference between the power values of the original data and the power values of the data with gaps for all data windows affected by gaps and frequency bin *f* , divided by the average power for frequency bin *f* in the data set without gaps. For example, a normalized bias of −0*.*1 means that the estimated power after application of the stimulus processing algorithm is on average 10% smaller than the power of the original data set if a gap is present.

The statistical evaluation of the spectral bias results in the simulations was performed as follows: we obtained the bias for each data set, resulting in 12 values, and performed a nonparametric Wilcoxon signed-rank test for zero median. If the *p*-value for this test was below 0.01, we regarded the bias as significant.

## **3. RESULTS**

First we show results of simulated gaps on data without stimulation to assess the influence of gaps and the stimulation-processing algorithms on the estimated spectrum. Then we illustrate the influence of real single TMS and epidural stimulation pulses on the spectrum if they are left untreated and how the methods of this paper deal with their after-effects. Finally, we apply the algorithms to data sets of BCI experiments with open-loop or closed-loop stimulation and investigate the effect of each method on the discrimination between the brain states during intended movements and rest.

#### **3.1. GAP SIZE**

**Figures 3A**–**C** show the influence of gap sizes between 5 and 100 ms on the error in spectral estimation for three particular frequencies (9, 21, and 81 Hz) and a model order of 32. We find that the RMSE increases with the gap size for all methods. This happens, because the information of the samples that are excluded by the gap is missing for the AR estimation, leading to a greater deviation from the AR coefficients without gaps for increasing gap size. The linear interpolation methods exhibit a negative bias and the AR-prediction shows positive bias (**Figures 3D**–**F** and **5**). The negative bias of the linear interpolation methods occurs because a section of the data window is reduced to a straight line which has a power of almost 0 for higher frequencies, leading to a decrease in the estimated power for these frequencies. This effect increases with greater gaps. AR modeling can lead to jumps in the data, because the extrapolated signal from the start of the gap is not necessarily connected to the actual recorded signal at the end of the gap. Such jumps result in higher estimated power across all frequencies and thus a positive bias. For longer gaps this bias increases because the potential deviation from the true values after the gap (the jumps) becomes larger. The mixture of linear interpolation and AR prediction is in general closer to 0 than the other two, but the sign of its bias depends on data packet size and gap size. The joining and MEMgap algorithms exhibit a bias close to zero, but the RMSE is smaller for MEMgap than for joining. The variance (**Figures 3G**–**I**) also scales with gap size but there are strong differences between the methods visible, with MEMgap and the linear interpolation methods having the lowest variance.

#### **3.2. MODEL ORDER**

Variations of the model order have the largest effect on the AR modeling and the MEMgap algorithm. While AR modeling exhibits a significant positive bias at 21 Hz for gaps longer than 60 ms at a model order of 16 (**Figure 4D**), it is not significantly biased for a model order of 32 and 64 (**Figures 4E**,**F**). As shown in section 3.3 and **Figure 5**, this is due to the frequency-dependency of the bias for AR modeling which has a global minimum around 20 Hz for model order 32 and 64. For MEMgap we find no significant bias for all model orders (**Figures 4D**–**F**) and that the absolute error of the power estimation, captured by the normalized RMSE, as well as the variance, increases rapidly with increasing model orders (**Figures 4A**–**C**,**G**–**I**). This is probably due to the lower number of samples fully available for AR estimation with MEMgap compared to the standard Burg algorithm: for MEMgap, forward or backward prediction errors can not be calculated for up to 2*p* samples around each gap, where *p* is the model order. Higher values of *p* only increase this difference, leaving MEMgap with less and less samples for AR estimation, thus probably leading to greater errors. In general, MEMgap has the lowest RMSE for orders 16 and 32 and gaps longer than 30 ms and the lowest RMSE of all methods with a bias close to 0 at an order of 64.

#### **3.3. FREQUENCY**

In **Figure 3**, we show the results for low and high frequencies with 9 and 81 Hz as parts of the μ and high γ bands, respectively, in

addition to the "intermediate" frequency of 21 Hz as part of the β-band. For 81 Hz, the linear interpolation methods already show a significant negative bias for gaps of 5 ms, whereas for 9 Hz this only becomes significant for gaps greater than 35 ms. This is easily understandable considering that one cycle of a 9 Hz oscillation lasts for more than 100 ms, therefore linear interpolation over a gap of 10–20 ms would be fairly consistent with the real shape of the undisturbed signal. The bias of MEMgap is not significant for any frequency (**Figures 3D**–**F**). The joining method on the other hand exhibits a negative bias for 9 Hz and gaps smaller than 40 ms and a significant positive bias for 81 Hz. For 21 Hz, The bias is significant only for gaps smaller than 10 ms. In terms of RMSE and variance (**Figures 3A**–**C**,**G**–**I**), MEMgap always displays the lowest values for gap sizes greater than 50 ms.

The results in **Figures 3D**–**F**, especially for AR modeling and joining, suggest that the bias might be frequency-dependent. In **Figure 5**, the bias is shown relative to the frequency bin for model orders of 16, 32, and 64 for a gap size of 100 ms where it should be most pronounced. We find that for the joining method, the bias is negative, although non-significant, for frequencies lower than 25 Hz and positive otherwise (significant for most frequencies *>*60 Hz). For AR modeling, the bias is in general positive (significantly for all frequencies for model order 16 and above 55 Hz for 32) and increases with frequency, has a minimum around 20 Hz for a model order of 32 and 64 and is also increased for lower frequencies. For linear interpolation, there is a bias close to −0.2, indicating a reduction in power of about 20%, for frequencies higher than 20 Hz. This value can be explained with the fact that 20% (100 ms of a 500 ms window) of the data had to be filled by linear interpolation which removes the high-frequency content. MEMgap exhibits no significant bias across all frequencies and model orders, except for very low frequencies and high model orders where all methods show a positive bias. Although the bias for the mixture of AR modeling and interpolation is also not significant for most frequencies above 10 Hz and higher model orders, this is due to the interaction between gap size and packet size for this particular method. As seen in **Figures 3D** or **F**, for example, a gap size of 80 would lead to a positive bias.

#### **3.4. APPLICATION ON DATA WITH STIMULATION**

In our experiments, we received the data in packets with a length of 40 ms. This leads to the jumps seen in the bias relative to the gap size for the combination of AR modeling and linear interpolation (e.g., **Figure 3F**, magenta line) as either linear interpolation of AR modeling dominate the outcome. The packet length might be different for other recordings, so we excluded this method from the rest of the experiments, as the conclusions would be very specific for our setup. Further experiments are needed to investigate the influence of this specific parameter. As the simulation results of the two versions of linear interpolation did not differ much, we restricted ourselves to four of the six methods for the remainder of the paper: online-compatible linear interpolation, AR modeling, joining, and MEMgap.

A model order of 16 was chosen for spectral estimation and AR extrapolation as this section is mostly illustrative in nature and serves the purpose to study whether the results from the simulations are transferable to data with actual stimulation. In terms of the estimation bias, we found the clearest effects for a model order of 16: a negative bias for linear interpolation and a positive bias for AR modeling. The latter bias was not present for higher model orders around the studied frequency bin of 21 Hz.

**Figure 6** illustrates the effect of epidural stimulation on the recorded ECoG activity, the evoked activity after stimulation and their influence on spectral estimation for one representative stimulation pulse. **Figure 6A** shows the raw trace of data with a single pulse of electrical epidural stimulation occurring at time point 0. **Figure 6B** displays a zoom on the first 100 ms after the pulse. The stimulation artefact itself is contained within the first 10 ms after the pulse. After that, one can find evoked activity with its peak occurring 13 ms after the pulse and an amplitude of 240μV. This is much higher than the short-term amplitude fluctuations found in our ECoG data without stimulation.

**Figure 6C** demonstrates the importance of adjusting the length of the gap to the actual stimulation effects on the signal. Applying a gap of 10 ms to the data might be enough to cover the stimulation artefact itself, but the spectrum then still shows a clear positive bias due to the influence of the evoked activity. The results for short gap sizes are very similar to those without any gap. Only gaps greater than 20 ms cover the extent of the artefact and the initial evoked activity, leading to power values that are similar to those obtained for data windows without the stimulation event (windows 16–26). There is no clear difference in the outcome of the gaps greater than 20 ms.

differences between the power at a gap size of 0 and 50 and above. Window numbers correspond to the brackets shown in **(A)**, where the first one is 1, the second one (shifted by 40 ms) is 2 and so on. The computation of windows −1, 0, and 15–26 used data that is outside the margins of **(A)**. **(D)** Comparison of linear interpolation (red), AR modeling (green), joining (blue), and MEMgap (gray) with gap sizes of 10 (dashed) and 50 (solid) applied on the data in **(A)**. The solid black line with circles in **(C)** and **(D)** shows the result of spectral estimation without processing of the stimulation after-effect (gap size = 0).

In **Figure 6D**, linear interpolation, joining the data segments, MEMgap, and AR modeling are compared when applied to the stimulation event both for short and long gaps. All methods perform poorly for a gap size of 10 ms, but there are differences for 50 ms. Applying joining and AR modeling results in higher power values than linear interpolation and MEMgap with a clear difference in estimated power between the windows with and without stimulation. Assuming perfect exclusion of all stimulus-related effects, we expect that the power does not differ strongly between e.g., window 15 (which includes a small portion of the gap) and window 16 (without the gap), therefore the result of MEMgap and linear interpolation is more realistic than the output of the other methods. At least for the AR modeling method, the increased estimate of the power compared to, for example, linear interpolation is consistent with the positive bias shown in **Figure 3B**. A reason for the positive bias of the joining method for this example data set might be that drifts of the signal after a stimulation on epidural electrodes are common. If we take a data segment with post-stimulus drifts, exclude the gap and join the data before and after the gap into one window, it will contain a sharp discontinuity and have a comparatively high spectral power. With linear interpolation, the discontinuity will be less severe and have a smaller impact on the signal power. For MEMgap it does not play a role as data before and after the gap is always separated during estimation of the AR coefficients.

Stimulation artefacts and evoked activity are found for combined EEG and TMS in a similar way as for stimulation over implanted electrodes with the strength of the evoked activity depending on the distance to the stimulation site. We illustrate this in **Figure 7** with the result of a TMS pulse on the activity recorded on a distant EEG channel. There is no strong evoked potential visible after the stimulation, therefore, as is evident in **Figures 7C**,**D**, a short window of 10 ms is already sufficient to cover the artefact and to produce an estimation of spectral power that is similar in value compared to that resulting for data windows long after the stimulation when using either linear interpolation, joining or MEMgap to correct for the gap.

#### **3.5. INFLUENCE ON DECODING PERFORMANCE**

The stimulation-processing algorithm can bias the estimated spectrum, or will at least produce deviations from the original spectral power without gaps. This poses the question, how strongly these errors influence the actual brain-state decoding during a BCI experiment. For example, if the bias of linear interpolation toward underestimation of the signal power directly influences, how well we can differentiate data packets obtained during a movement from those recorded during rest, then this algorithm is not suitable for BSDS because it might induce a bias in the subject's performance in an online experiment.

To investigate this, we used data sets with different stimulation paradigms and recording methods (EEG and ECoG) to assess the influence of the algorithms and gap size on the decoding abilities of a BCI system. The patients always performed the

same cued attempted hand movements but we varied the stimulation paradigm between no stimulation, stimulation with fixed ISI and stimulation coupled to the output of the BCI (i.e., BSDS). In the last paradigm, the stimulation pulses were only applied while the BCI detected an intention to move from modulations of the power in the β-band and therefore moved the orthosis. If stimulation was used, stimulation artefacts were identified online with a peak detector if the voltage of two consecutive samples differed by more than 1 mV. The start of the gap was set 2 ms before this artefact and the gap size was adjusted for each patient and session depending on the length of the evoked activity as determined by several test stimuli applied before the start of the session. This resulted typically in a gap length between 30 and 70 ms. Stimulus processing was performed during recording with the online-compatible linear interpolation method.

In the offline analysis, we applied the four methods: joining, linear interpolation, AR modeling, and MEMgap on these data sets and varied the gap size between 0, 10, 50, and 100. We simulated the two different stimulation conditions on the data without stimulation by varying, in which phases of the trial gaps are placed: in the uncoupled condition the whole trial was valid, so the placement of gaps was independent of the activity and brainstate of the patient. For the coupled condition only time points within the movement phase were used as gaps, thus simulating a BSDS paradigm. In both cases the ISI was fixed at 2 s.

After applying the respective stimulation processing algorithm, we computed the spectral power between 16 and 22 Hz on channels located over the right motor cortex. For EEG measurements we used FC4, C4, and CP4 as defined by the 10–10 system (Society, 2006), whereas for ECoG measurements the electrodes were selected individually per patient based on the results of a screening session. We used a window size of 500 ms and a model order of 16 for spectral estimation. These were the same parameters, channels and frequencies that had been used during the online feedback experiments in which the data was recorded. Furthermore, our simulations showed a positive estimation bias for AR modeling at 16–22 Hz only for a model order of 16, not for 32 or 64. Thus, we only used an order of 16 for the simulations on data without stimulation. In order to investigate, whether higher model orders have a substantial effect on the processing of real stimuli, we used model orders of 16, 32, and 64 on the data with open-loop and closed-loop stimulation. For each run (consisting of 11 trials), we calculated the area under the ROC curve (AUC) for the sum of the logarithm of the power values within each data window in the movement phase versus those in the rest phase. We used this as a measure of the separability of these phases on a single-packet level. Taken together from all three patients, we analyzed 87 runs of EEG recordings without stimulation, 24 runs with uncoupled EEG-TMS, 131 ECoG runs without stimuli, 51 runs of ECoG with uncoupled, and 82 runs of ECoG with coupled stimulation. For each recording and stimulation condition, algorithm and gap size, this resulted in a distribution of AUC scores, one per run.

The conditions without stimulation allowed us to test for the bias and absolute error introduced by the gaps into the AUC scores. Thus, we computed the pair-wise differences between the AUC scores of a gap size of 0 and those of all combinations of algorithms and gap size for these conditions. Using Kruskal– Wallis tests, Bonferroni-corrected for multiple comparisons, we tested which algorithm leads to the smallest absolute differences in AUC scores. We also applied Wilcoxon signed rank tests to assess, whether the median of the differences deviates significantly from 0, indicating a systematic bias in the AUC scores. As there is no "true" reference distribution of the AUC scores possible for data with stimulation, we used Bonferroni-corrected non-parametric Friedman tests which account for possible effects of using the same sessions in all conditions to test whether gap sizes greater than 0 lead to different AUC scores compared to a gap size of 0 and to test whether there is a difference between the algorithms at a certain gap size.

**Figure 8** shows data without any actual stimulation, so ideally the difference in AUC scores between a gap size of 0 and of 100 should be zero for all runs. In **Figures 8A**,**B**, stimuli were simulated throughout the trial, thus independent of the task or the output of the BCI. Session-wise comparison of the AUC values with Friedman tests for each gap size show significant differences between MEMgap and the other algorithms only for long gaps. There is a slight decrease in the average AUC value for all algorithms for a gap size of 100 compared to 0 and MEMgap and joining yield significantly smaller absolute differences in AUC scores compared to AR modeling and interpolation (*p <* 0*.*000001). In **Figures 8C**,**D**, the stimuli were simulated only throughout the movement phase. We find a significant decrease in AUC values for AR modeling and a significant increase for linear interpolation. This means that linear interpolation artificially "improves" the decoding power. As shown in the simulation studies, linear interpolation of large gaps leads to a decrease of the power between 16 and 22 Hz (negative bias) which increases the event-related desynchronization effect of sensorimotor rhythms during attempted movements (Wolpaw et al., 2002). The MEMgap method shows a significantly smaller deviation of the AUC values at gap sizes of 50 and 100 from the AUC values without gaps than all the other methods (*p <* 0*.*000001). In contrast to the other algorithms, the median of the AUC differences after MEMgap never differs significantly from 0 for the BSDS condition, except for the ECoG data set with a gap size of 100 (*p* = 0*.*005). The other algorithms differ significantly from MEMgap in almost all cases of the coupled conditions.

Patients in the data sets shown in **Figures 9A**,**B** were stimulated independent of the task. We show only a model order of 16, because the results for an order of 32 and 64 are very similar. It is evident from gap sizes of 0 and 10 that untreated stimulation after-effects are detrimental for decoding. Online decoding will be more successful if enough samples are excluded after a stimulus (in these examples: a gap size of 50 ms seems to work well, although this varies between patients). Using Friedman tests for session-wise comparison of the AUC scores, significant differences of the algorithms are found, although the mean absolute differences are very small (≤0*.*01). In case of the uncoupled ECoG condition, AUC scores with untreated stimulation after-effects are significantly lower than AUC scores for gap sizes of 50 and 100, independent of the applied algorithm (*p <* 0*.*01). This effect is due to residual stimulation after-effects for small or absent gaps that lead to a very high power of data windows that contain

electrical or magnetic pulses. In particular such data windows in the movement phase will be classified incorrectly. If the strong after-effects are removed by longer gaps, the classifier is more likely to produce a correct result which is reflected in the increased AUC score for gaps of at least 50 ms.

Finally, in **Figure 9C**, stimulation was given only during the movement phase. The average AUC value for a gap size of 0 is smaller than 0.5, indicating a higher power during movement than during rest, as opposed to the expected event-related desynchronization. This is due to the task-dependent existence of the stimulation effects: the large stimulation after-effects that occur only during the movement phase lead to a very high spectral power of this phase. Thus, the spectral power of the movement phase is very well separable from the power of the rest phase for a gap size of 0. For a gap size of 10, there is a large variability in the AUC scores. This is because for one of the three patients, a gap of 10 ms was not sufficient to cover all artefact-related jumps in the recording, resulting in AUC scores lower than 0.5. If the aftereffects are dealt with by using a gap size of 50, the relationship between the power during rest and feedback reverses and resembles the expected ERD/ERS pattern. For a gap size of 100, we find in **Figure 9C** that the largest average AUC value is reached for linear interpolation and the smallest one for AR modeling, both differing significantly from the AUC values for MEMgap (*p <* 0*.*000001). This relationship is found for all tested model orders, where joining and AR modeling are on average worse than MEMgap by more than 0*.*02 and 0*.*05, respectively, while linear interpolation yields higher scores by at least 0*.*01. This is consistent with the simulation results in **Figures 8C**,**D** indicating an artificial over- and under-estimation of class separability by these methods. It supports the hypothesis that MEMgap is probably best suited to deal with large gaps in the data, especially for BSDS, because based on the simulation studies the deviation from the true AUC value is significantly smaller than for the other methods.

## **4. DISCUSSION**

One challenge when trying to combine online brain-state decoding from spectral data and direct cortical stimulation is that the after-effects of stimulation such as artefacts (Taylor et al., 2008) or evoked activity (Matsumoto et al., 2007) can have a much higher amplitude than the background brain signals. Therefore, estimation of the brain-state from a segment of data has been unreliable,

if such stimulation after-effects are contained in this segment. This leaves us with three options: we can (1) use only data segments for decoding that are free of any after-effects, (2) attempt to separate stimulation after-effects from background brain activity, e.g., by fitting a template of the expected shape of the effects to the recording, or (3) isolate the portions of the data segment that are "contaminated" by stimulation effects and use only the "clean" parts for decoding.

level for a purely random classifier (Fawcett, 2006). **(A)** Average AUC values for the separation of movement and rest from experiments with epidural

In earlier studies combining TMS and EEG without BSDS (i.e., without the necessity to perform real-time brain-statedecoding from the EEG), options (1) and (2) have been used. In such studies, either a fixed length window around the stimulus was removed offline (Fuggetta et al., 2006), a decomposition into artefact-free and contaminated data was attempted in postprocessing (Litvak et al., 2007; Morbidi et al., 2007; Erez et al., 2010) or a sample-and-hold circuit was used during recording to fix the amplifier output at a constant level during the pulse (Ilmoniemi et al., 1997). The latter method is especially helpful for amplifiers that recover from TMS pulses only after a delay of several hundred milliseconds (Ilmoniemi and Kiciˇ c, 2010 ´ ), although some current amplifiers are able to keep this delay lower than 10 ms (Veniero et al., 2009). The drawback of the sample-and-hold approach is that information on the brain signal directly after the pulse is invariably lost and that the signal contains gaps.

Option (1) has also been used by Bergmann et al. (2012) in their study on EEG-guided TMS, making a waiting period of several seconds between stimulation pulses necessary. If the brainstate is decoded from spectral features and for example 500 ms of data is needed to estimate these features robustly, one has to wait for 500 ms plus the expected duration of the stimulation after-effects for making the first estimate of the brain-state after a stimulation pulse. This duration is therefore also the absolute minimum ISI in this scenario. Removal of the after-effects by template subtraction is only possible, if several constraints are met: the full amplitude range of the stimulation effects has to be within the dynamic range of the amplifier, as portions of the data in which the amplifier is in saturation can not be recovered with this method, resulting in the necessity to correct for gaps in the signal as in option (3). If the recorded effects are not sufficiently stable, attempting to remove them will lead to residuals in the signal. Like the original after-effects, these residuals can have a detrimental effect on the quality of the estimated spectrum and, thus, the decoding process. The employed removal algorithms need to be suitable for an online BCI, so they need to work on a single-trial level and therefore should not be too computationally demanding.

for this algorithm and gap size differ significantly from MEMgap (*p <* 0*.*05,

Friedman test, Bonferroni-corrected).

We have chosen approach (3) for this work, the deliberate introduction of gaps into the signal covering the strongest aftereffects of stimulation and correcting for these gaps during spectral estimation. This allows continuous decoding without influence of the stimulation after-effects, as long as the duration of the aftereffects is estimated properly. To apply this approach, methods are needed that do not depend on continuous data segments for brain-state decoding and can deal with gaps in the data.

In our experiments, we analyzed the spectral power in the μ- and/or β-band to detect the patient's intention to move the paralyzed hand. We compared different approaches (linear interpolation, AR modeling, joining of data segments, and the Burg algorithm adapted for segmented data) on their ability to estimate the spectrum with gaps in the data. To this end, we used an ECoG BCI training data set and analyzed the normalized RMSE, bias and variance of the difference between the estimated spectrum with and without gaps. The RMSE increased with the gap size, although the slope of the error increase was smaller for MEMgap and joining than for algorithms that fill the gap with artificial data (linear interpolation and AR modeling). We found a clear systematic negative bias for linear interpolation and a systematic positive bias for AR modeling. We studied the frequency range between 16 and 22 Hz in most detail, where the bias of AR modeling was only apparent for a model order of 16, but a clear bias of AR modeling can be found for other frequencies at higher model orders, making this method also potentially unreliable. The joining method produces a bias close to 0 around a frequency of 20 Hz, but can lead to a positive bias for higher frequencies, whereas the MEMgap method always results in a bias close to 0. For gaps smaller than 40 ms, linear interpolation typically has the smallest absolute deviation from the true power values while MEMgap outperforms the other methods for longer gaps.

As our simulations show, the RMSE for linear interpolation is smaller than for MEMgap, thus at first glance making linear interpolation superior to MEMgap for large model orders and/or small gaps. However, in the context of a continuous BCI decoding for BSDS, the negative bias exhibited by the linear interpolation methods will bias the output of the BCI in favor of ERD, thus distorting the real performance of the participant to some extent if stimulation is coupled to the detected brain-state. We therefore think that MEMgap is most suited for BSDS as it is superior or at least equal to the other methods in terms of RMSE and variance, does not introduce a systematic bias and outperforms the other methods in minimizing the stimulation after-effects in our BCI paradigm.

Whether this approach of identifying and ignoring the segments of data dominated by stimulation after-effects is feasible in any given experimental setting depends on the duration of stimulation-evoked potentials after the pulse. As we showed here in the simulation studies, if the strongest evoked activity is contained within the first 50–100 ms after the pulse, then a decoding approach using MEMgap is feasible. If no strong evoked activity is observed, e.g., in the case of a remote recording location as illustrated in **Figure 7**, then a short gap of 10 ms covering the

## **REFERENCES**

Kotchoubey, B., Kübler, A., et al. (1999). A spelling device for the paralysed. *Nature* 398, 297–298.

stimulus artefact together with linear interpolation or MEMgap would be sufficient. In Ferreri et al. (2011), evoked EEG activity following single pulse TMS was found for up to 300 ms after the pulse with amplitude fluctuations of less than 20μV for late components. Although we do not expect that such small potentials would have a large impact on the estimated spectrum, especially compared to the stimulation artefact itself or early evoked activity, for every experiment of BSDS with continuous decoding, the size and shape of the evoked activity should be carefully studied to get a proper estimate of the duration of strong after-effects. As was shown by Casarotto et al. (2010), the after-effects depend on a number of parameters, such as stimulation intensity, location and (in the case of TMS) coil orientation. If gaps longer than the 100 ms tested here are necessary to cover all stimulation-related activity, one should either wait long enough until all effects have ceased before making the next brain-state decoding attempt, or increase the size of the data window on which the spectrum is estimated to ensure that it contains enough clean samples to compute a valid estimate.

In conclusion, we have shown that the application of cortical stimulation coupled to the output of an online brain-state decoder based on spectral features is feasible as long as the employed algorithms remove both the stimulation artefact and large early components of evoked activity and allow spectral estimation on non-continuous data. Especially if closed-loop BSDS is used, algorithms that do not introduce a strong bias into the estimated spectrum such as MEMgap are to be preferred over biased methods like linear interpolation to ensure a reliable decoding of the brain-state. In general, the methods investigated here are not restricted to applications with cortical stimulation but can be employed whenever spectral estimation has to be performed on non-continuous data sets with missing blocks of samples.

#### **ACKNOWLEDGMENTS**

This work was funded by the European Research Council (ERC) grant 227632 "BCCI—A Bidirectional Cortical Communication Interface" and the Reinhardt-Koselleck Award of the Deutsche Forschungsgemeinschaft (DFG). We thank Alireza Gharabaghi, Maria Teresa Leão, Georgios Naros, Alexandros Viziotis, and Martin Spüler for insightful discussions and their help in data acquisition and the two reviewers for helpful comments on the manuscript.

*ONE* 5:e10281. doi: 10.1371/journal.pone.0010281

G., et al. (2011). Human brain connectivity during single and paired pulse transcranial magnetic stimulation. *Neuroimage* 54, 90–102.

Neuronal responses to magnetic stimulation reveal cortical reactivity and connectivity. *Neuroreport* 8, 3537–3540.

hemiparetic stroke: a multicenter feasibility study of safety and efficacy. *J. Neurosurg.* 108, 707–714.

*Handbook of Transcranial Stimulation (Oxford Handbooks)*. New York, NY: Oxford University Press.

Wolpaw, J. R., Birbaumer, N., McFarland, D. J., Pfurtscheller, G., and Vaughan, T. M. (2002). Brain-computer interfaces for communication and control. *Clin. Neurophysiol.* 113, 767–791.

**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.

*Received: 06 August 2012; accepted: 29 October 2012; published online: 16 November 2012.*

*Citation: Walter A, Murguialday AR, Rosenstiel W, Birbaumer N and Bogdan M (2012) Coupling BCI and cortical stimulation for brain-state-dependent stimulation: methods for spectral estimation in the presence of stimulation aftereffects. Front. Neural Circuits 6:87. doi: 10.3389/fncir.2012.00087*

*Copyright © 2012 Walter, Murguialday, Rosenstiel, Birbaumer and Bogdan. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

## **1. APPENDIX**

#### **1.1. THE BURG ALGORITHM**

The Burg algorithm is used for the estimation of the coefficients *ci* of an autoregressive (AR) model

$$\varkappa(t\_p) = \sum\_{i=1}^{p} c\_i \varkappa(t\_p - i) + e$$

with order *p* for samples *x(tk)*, 0 ≤ *k < N*,*e* a sample from a white noise sequence. The algorithm needs *p* recursive steps and in each step *j*, the coefficients *cj,i* for an autoregressive model of order *j* are computed by the following procedure:

An initial estimation of the power of the white noise component in the AR model is obtained by

$$P\_0 = \frac{1}{N} \sum\_{k=0}^{N-1} |\mathbf{x}(t\_k)|^2.$$

Each new coefficient *ci,i* is computed by minimizing the forward and backward prediction errors

$$f\_{\mathfrak{P},k} = \mathfrak{x}(\mathfrak{z}\_k) - \sum\_{i=1}^{\mathfrak{p}} c\_{\mathfrak{p},i} \mathfrak{x}(\mathfrak{z}\_{k-i}) \qquad \qquad \text{with } k = \mathfrak{p}, \dots, N - 1$$

$$b\_{p,k} = \mathbf{x}(t\_{k-p}) - \sum\_{i=1}^{p} c\_{p,i} \mathbf{x}(t\_{k-p+i}) \quad \text{with } k = p, \dots, N-1$$

with the formula

$$c\_{i,i} = \frac{-2\sum\_{k \in I\_i} f\_{i-1,k} \cdot b\_{i-1,k-2}}{\sum\_{k \in I\_i} \left( |f\_{i-1,k}|^2 + |b\_{i-1,k-1}|^2 \right)}, \quad I\_i = \{i+1, \dots, N-1\}. \tag{A1}$$

Each previously computed coefficient *ci* − 1*,k* is then adjusted by

$$c\_{i,k} = c\_{i-1,k} + c\_{i,i} \cdot c\_{i-1,i-k}$$

and we update the power estimation to

$$P\_i = (1 - \left|c\_{i,i}\right|^2) \cdot P\_{i-1}$$

and the forward and backward prediction errors:

$$f\_{i,k} = f\_{i-1,k} + c\_{i,i} \cdot b\_{i-1,k-1}$$

$$b\_{i,k} = b\_{i-1,k-1} + c\_{i,i} \cdot f\_{i-1,k} \cdot b\_{i,k}$$

After *p* steps, this results in the AR coefficients *ci* = *ci,p*, *i* = 1*,..., p*.

#### **1.2. ESTIMATING THE SPECTRUM FROM AN AR MODEL**

An AR model can be interpreted as an all-pole infinite-impulseresponse filter with order *p* and coefficients *ci* which is applied to a white noise process with a power of *Pp* (Pardey et al., 1996). Thus, after finding the *p* autoregressive coefficients *ci*, one can estimate the spectrum by evaluating the transfer function *H(z)* =

$$\sqrt{P\_{\mathcal{P}}} \left(1 - \sum\_{k=1}^{\mathcal{P}} c\_k z^k \right)^{-1} \text{ of the filter to find power values}$$

$$P(\omega) = \frac{P\_p}{|1 - \sum\_{k=1}^{p} c\_k e^{-jk\alpha}|^2}$$

at (normalized) frequencies ω.

#### **1.3. THE MEMgap ALGORITHM**

If one assumes that a sequence *g* of length *N* exists and that *g(n)* = 1 only if the corresponding sample *x(tn)* is part of a gap in the data and 0 otherwise then we just have to make sure that none of the samples *x(tn)* with *g(n)* = 1 influence the estimation of the model coefficients. The Burg algorithm computes the AR coefficients for order *p* in *p* steps, yielding in the *i*-th step the coefficients of an AR model with order *i*. If we use in the *i*-th step only those samples fully for computation of the AR coefficients that are at least *i* + 1 time steps away from a sample with *g(n)* = 1, we achieve the desired effect. To be more precise, the coefficients are computed by evaluating forward and backward prediction errors (see Appendix 1.1). In the MEMgap algorithm, forward prediction errors are only computed for samples that are at least *i* + 1 time steps after a gap, backward prediction errors only for those at least *i* + 1 time steps before a gap. Formally, this is done by modifying *Ii* in Equation (A1) to the set

$$I\_i = \{k \mid \mathbf{g}(k) = 0 \land i < k < N \land \}$$

$$(k - n < 0 \lor k - n > i) \,\forall n \text{ with } \mathbf{g}(n) = 1\}.$$

This set can also be computed iteratively in each step of the Burg algorithm as *Ii* = *Ii*−1 ∩ *I*- *i*−1, where *I*- *i*−1 is the set *Ii*−1 with each entry incremented by 1 and *I*0 = {*k* | *g(k)* = 0 ∧ 0 *< k < N*}. This resembles a "forbidden zone" that initially contains only the gaps but grows in each step of the algorithm by one sample. The estimation of the white noise power *P*0 has to be calculated only with samples outside of gaps: *P*0 = 1 |*I*0| *k*∈*I*0 |*x(tk)*| 2. The rest of the algorithm works as the standard Burg algorithm described in Appendix 1.1.

Obviously, this method can only work if the continuous segments between gaps are long enough. Therefore, there needs to be at least one segment of samples with a length that is at least equal to the model order *p* in order to make an estimation of the spectrum with this method. In practice, it is preferable if the number of samples in such a segment is several times higher than the model order in order to reduce the bias and variance of the estimator. If there are *ns* segments of clean data with a length greater than *p* and the total number of samples in these segments equals *M*, then only *M* − 2*nsp* forward and backward prediction errors are available in the *p*-th step of the MEMgap algorithm, although all *M* samples are evaluated to compute these errors. This means that even if the total number of samples within gaps might be the same, one can expect that the variance of the spectral estimation will be smaller if there are only a few large gaps in the data compared to having many small gaps because less samples contribute fully in the second case. According to de Waele and Broersen (2000), the same holds for the estimation bias which is inversely proportional to the number of available samples.

#### **1.4. AR MODELING INDUCES JUMPS IN THE SIGNAL**

The approach to fill the gap with samples that were extrapolated by an AR model fitted to the data before the gap can be problematic, if the extrapolation diverges strongly from the measured data. When actually measured samples are added to the data buffer after the gap, there can be a large amplitude difference between the last (extrapolated) sample within the gap and the first measured sample after the gap (**Figure A1A**). Such a "jump" in the signal results in high power across all frequencies, thus distorting the spectrum. To assess the influence of the model order and the gap size on the jumps, we ran simulations on the data from the ECoG recordings without stimulation used in **Figure 8A**: in total, 11266 stimuli were simulated on 397 minutes of data, the gap size was varied between 0 and 100 ms in steps of 5 ms and model orders 16, 32, and 64 were tested.

We applied AR modeling to deal with the gaps and measured the jump as the absolute voltage difference between the last sample generated by the AR model at the end of the gap and the first sample after the gap. The results are shown in **Figure A1B**.

and 600 ms. An AR model of order 32 is estimated from the 500 ms before the gap and applied as a linear predictor to generate 100 samples to fill the gap (red). The original ECoG samples within For comparison, the average absolute voltage difference between neighboring samples of the original ECoG data without any gaps or stimuli is 1*.*68 ± 0*.*43μV (mean ± std).

We found that for all model orders, the average height of the jump increases sharply up to a gap size of 10, yielding 12*.*35 ± 5*.*0μV for order 16, 11*.*82 ± 4*.*43μV for 32 and 11*.*29 ± 4*.*37μV for 64. For further increasing gap size, the average jump height increases more slowly for higher model orders than for lower ones. For a gap size of 100 ms we find average jump heights of 28*.*08 ± 14*.*94μV for order 16, 23*.*45 ± 12*.*51μV for 32 and 20*.*23 ± 11*.*94μV for 64. Thus, while the jump height at the end of the gap is significantly smaller for a model order of 64 compared to 32 and 16 (gap size = 100, paired Wilcoxon signed rank tests, both *p <* 10−17), it is still vastly higher than the average sample-tosample difference for ongoing ECoG activity. Therefore, we cannot conclude that higher model orders prevent jumps after the gap.

Furthermore, if the gaps are used to cover the effects of real stimulation, there has to be a jump at the end if the gap is filled by AR modeling and evoked activity is present. AR modeling attempts to extrapolate the pre-stimulus signal which almost certainly differs in its time course from the stimulation-evoked activity, therefore extrapolation can not work perfectly, regardless of the model order.

gap is the jump height. **(B)** Average jump height after filling the gap with AR modeling for model orders 16 (black), 32 (red), and 64

(green).

## Erratum: Coupling BCI and cortical stimulation for brain-state-dependent stimulation: methods for spectral estimation in the presence of stimulation after-effects

*Armin Walter1 \*, Ander Ramos Murguialday2,3, Martin Spüler 1, Georgios Naros 4,5, Maria Teresa Leão4,5, Alireza Gharabaghi 4,5, Wolfgang Rosenstiel 1, Niels Birbaumer 2,6 and Martin Bogdan1,7*

*1 Department of Computer Engineering, Wilhelm-Schickard-Institute, Eberhard Karls University Tübingen, Tübingen, Germany*

*3 Health Technologies Department, TECNALIA, San Sebastian, Spain*

*4 Department of Neurosurgery, University Hospital Tübingen, Tübingen, Germany*

*5 Department of Neuroprosthetics, Centre for Integrative Neuroscience, Eberhard Karls University Tübingen, Tübingen, Germany*

*6 IRCCS, Ospedale San Camillo, Venice, Italy*

*7 Department of Computer Engineering, University of Leipzig, Leipzig, Germany*

*\*Correspondence<EMAIL_ADDRESS>
#### *Edited by:*

*Ahmed El Hady, Max Planck Institute for Dynamics and Self Organization, Germany*

#### *Reviewed by:*

*Ahmed El Hady, Max Planck Institute for Dynamics and Self Organization, Germany*

#### **A commentary on**

#### **Coupling BCI and cortical stimulation for brain-state-dependent stimulation: methods for spectral estimation in the presence of stimulation after-effects**

*by Walter, A., Murguialday, A. R., Rosenstiel, W., Birbaumer, N., and Bogdan, M. (2012). Front. Neural Circuits 6:87. doi: 10.3389/fncir.2012.00087*

The article "Coupling BCI and cortical stimulation for brain-state-dependent stimulation: methods for spectral estimation in the presence of stimulation after-effects" (Walter et al., 2012) incorrectly omitted four authors. The list of authors should have been: Armin Walter, Ander Ramos Murguialday, Martin Spüler, Georgios Naros, Maria Teresa Leão, Alireza Gharabaghi, Wolfgang Rosenstiel, Niels Birbaumer, Martin Bogdan.

All further references to this paper should cite the corrected author list as displayed above.

Furthermore, the text of the "Acknowledgments" section needs to be changed to: This work was funded by the European Research Council (ERC) grant 227632 "BCCI—A Bidirectional Cortical Communication Interface" and the Reinhardt-Koselleck Award of the Deutsche Forschungsgemeinschaft (DFG). We thank Alexandros Viziotis for his work with the patients and the two reviewers for helpful comments on the manuscript.

#### **REFERENCES**

Walter, A., Murguialday, A. R., Rosenstiel, W., Birbaumer, N., and Bogdan, M. (2012). Coupling BCI and cortical stimulation for brain-state-dependent stimulation: methods for spectral estimation in the presence of stimulation after-effects. *Front. Neural Circuits* 6:87. doi: 10.3389/fncir.2012. 00087

*Received: 09 January 2013; accepted: 09 January 2013; published online: 11 February 2013.*

*Citation: Walter A, Murguialday AR, Spüler M, Naros G, Leão MT, Gharabaghi A, Rosenstiel W, Birbaumer N and Bogdan M (2013) Erratum: Coupling BCI and cortical stimulation for brain-state-dependent stimulation: methods for spectral estimation in the presence of stimulation after-effects. Front. Neural Circuits 7:6. doi: 10.3389/fncir.2013.00006*

*Copyright © 2013 Walter, Murguialday, Spüler, Naros, Leão, Gharabaghi, Rosenstiel, Birbaumer and Bogdan. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

*2 Institute of Medical Psychology and Behavioural Neurobiology, University Hospital Tübingen, Tübingen, Germany*