Band power modulation through intracranial EEG stimulation and its cross-session consistency

We used data that are publicly available as part of the Restoring Active Memory project (managed by the University of Pennsylvania; http://memory.psych.upenn.edu/RAM). As stated in the project's website 'Informed consent has been obtained from each subject to share their data, and personally identifiable information has been removed to protect subject confidentiality'. The original research protocol for data acquisition was approved by the relevant bodies at the participating institutions. Furthermore, the University Ethics Committee at Newcastle University approved the current project involving the data analysis reported here (Ref: 12721/2018). We extracted data from all patients (n = 87) that underwent at least one stimulation session while performing memory tasks. We excluded eight patients that either had substantial stimulation artefacts in almost all channels or their data were limited (single session with < 18 stimulation trials). Thus, we analysed data from 79 subjects (mean age 35.4, range 16–63, 44 females) from which 36 had at least 2 stimulation sessions with the same stimulation location (totalling 101 pairs of sessions with same stimulation location).

Stimulation was delivered using charge-balanced biphasic rectangular pulses (300 µs pulse width) at 10, 25, 50, 100, or 200 Hz frequency 0.25–3.5 mA amplitude. The duration of the stimulation was 500 ms or 4.6 s, depending on the experiment. Note that, for some subjects the amplitude of stimulation varies from session to session and we analyse it in figure 4 and account for it in figure 6. Stimulation duration does not vary between sessions within a subject. Stimulation frequency varies in 2 out of the overall 79 subjects but these subjects were not used to evaluate cross-session consistency. Sessions with stimulation duration 4.6 s always have 50 Hz stimulation frequency.

To measure stimulation effect, 1-second segments were extracted from the iEEG signals around every stimulation trial; that is, we extracted one segment before (pre) and one after (post) the stimulation event, with a 50 ms buffer between each segment and the event (following [15]). To assess baseline fluctuations, 'pre' and 'post' segments were also extracted from the baseline activity during baseline epochs, with a pre-post interval equal to the one around the stimulation trials of the same session. A baseline epoch was considered to be any inter-stimulus interval which was at least 20 s long and 5 s away from the stimulation itself. Figure 1 shows a schematic of the session timeline and the process of segment extraction. Since the stimulation trials were temporally organised in groups of three in a typical session (i.e. less than 10 s interval between trials in the same group), we extracted baseline pre/post segments from each baseline epoch in groups of three as well (see figure S1 in Supplementary Material (https://stacks.iop.org/JNE/17/054001/mmedia)), such that the number of segments taken around stimulation and the number of baseline segments were approximately equal in each session.

Zoom In Zoom Out Reset image size

The time series of each segment were centred around zero and de-trended. De-trending was achieved by applying linear regression and then removing the least-squares fit from the signal. Any channels with repeated artefacts were excluded (see below). A common average re-referencing was applied to the remaining set of channels. The stimulation channels were excluded from the common average calculation, but the calculated common average was applied to them. The band power of each segment was calculated in five different bands [delta (2–4 Hz), theta (4–8 Hz), alpha (8–12 Hz), beta (12–25 Hz), and gamma (25–55 Hz)] after estimating the power spectral density of the segment using Welch's method (with window length equal to half of the segment length and overlap length equal to a quarter of the segment length). Finally, the band powers were log-transformed. Figure 1 shows a schematic of the preprocessing.

Channels with repeated stimulation artefacts (i.e. voltage deflection) were excluded. A repeated stimulation artefact was detected based on two criteria. Either one of these two criteria was sufficient to indicate a channel with repeated artefacts. First, a strong effect of stimulation on the average (across time) voltage of the first half of the post segment compared to the average (across time) voltage of the second half of the pre segment. The effect was quantified using the t-statistic of a paired t-test. Second, the second half of the average (across trials) post signal had a slow return to the average (across time) voltage value of the pre segments. This was detected by linear regression.

Box plots were used to summarise various distributions in the Results. Central lines indicate median values, while the boxes extend from the 25th to 75th percentile (interquartile range) of the distribution. Whiskers extend to the upper and lower adjacent values, that is, the most extreme values that are not outliers. Outliers are considered to be values that lie more than 1.5×[interquartile range] away from the 25th or 75th percentile.

The effect of stimulation on band power from pre to post was considered as the z-statistic (indicated by U throughout) produced by the Wilcoxon sign rank test (paired non-parametric test; signrank function in MATLAB). A positive U indicates an increase in band power from pre to post, whereas a negative U indicates a decrease from pre to post. To also quantify the baseline fluctuations of band power, the same measure was used on the 'pre/post' pairs of the baseline activity (figure 1).

The overall difference in stimulation effect between two sessions (across all channel/band combinations) in figure 4(B) was quantified by using the absolute t-statistic of a paired t-test on the absolute effect U of the two sessions. We used absolute effects as we wanted to generally assess changes in effect magnitude.

The consistency of stimulation effect was measured for each pair of sessions with the same stimulation location in the same subject. All possible combinations of two sessions were considered, totalling 101 pairs. The consistency was computed by first pairing the effect values of corresponding channel/band combinations between the two sessions. Note that only the intersection of valid channels between the two sessions was considered (a channel can be excluded due to artefacts in one session but not the other). The consistency coefficient was given by the Fisher-transformed zero-centred Pearson's correlation. Considering the effect values of the two sessions as random variables S₁ and S₂, then the consistency coefficient is given by: $r_0 = \mathbf{E}[S_1S_2]/(\hat{\sigma}_1\hat{\sigma}_2)$, where $\mathbf{E}$ denotes expected value and $\hat{\sigma}$ refers to the average deviation from 0 ($\hat{\sigma} = \sqrt{(\sum^n_{i = 1}{s^2_i})/n}$). We use the zero-centred Pearson's correlation to only detect a zero-translated agreement between the random variables, that is, in the form of $S_1 = kS_2$, with 0 intercept and k a non-zero constant.

The consistency curve was used to express the consistency between two sessions by gradually considering fewer pairs of effect values at low effect sizes. Considering a scatter plot of all the effect value pairs, it was computed by gradually increasing the radius of an exclusion circle emanating from (0, 0). The consistency curve at radius = 0 gives the consistency when all points are included in the consistency calculation, whereas the consistency curve at radius = x expresses the consistency as computed after excluding every pair of effect values that lie inside a circle with centre (0, 0) and radius x. The circle was gradually enlarged with a step of 0.2 and the enlargement stopped just before covering 98% of the scattered values. We used this procedure to ensure that we can detect consistency even if only a few channels exhibited consistency, without the consistency being masked by low effect channels.

Each consistency curve is represented by its maximum consistency coefficient. The maximum consistency coefficient is the value on the curve that deviates the most from 0, being positive or negative. Thus, it expresses the strongest correlation or anti-correlation found between the effect values of the two sessions.

To explore which factors determine consistency across all 101 session pairs, we modelled the maximal value of consistency as a linear combination of the following variables:

• session time difference: absolute time difference between the sessions' starting timestamps.
• difference in baseline (band power) means: mean absolute paired difference between the sessions' mean values of band power during baseline (both 'pre' and 'post').
• difference in baseline (band power) standard deviations: mean absolute paired difference between the sessions' standard deviations of band power during baseline (both 'pre' and 'post').
• average max effect: average (between sessions) maximum effect (across all channel/band combinations).
• average min effect: average (between sessions) minimum effect (across all channel/band combinations).
• average stimulation amplitude: average stimulation amplitude between sessions.
• stimulation amplitude difference: difference in stimulation amplitude between sessions.
• stimulation frequency: frequency of stimulation pulse train (always common between examined session pairs).
• depth of the stimulation location: distance of stimulation location (midpoint between anode and cathode) from brain surface.
• task difference: difference in memory tasks performed by the subject during recording; that is, 0 for same and 1 for different tasks between sessions (categorical variable).

The stimulation depth was computed as the Euclidean distance of the anode-cathode midpoint from the subject's brain surface. If that midpoint was found to be outside the provided surface, its depth was set to negative (minus the Euclidean distance).

In order to quantify the explanatory power of all the different independent variables on the consistency we used ANOVA test on the model built by the Multiple Linear Regression Analysis. We built the model and assessed the ANOVA effects 200 times through bootstrapping. We used this bootstrapping approach to check for the robustness of the model. The ANOVA effect, the R², and the Adjusted R² are reported.

The code used for the analysis of the pre-processed data to produce the results is available online https://github.com/cnnp-lab/iEEGstim_consistency.git.

Figure 2 shows the measured stimulation effects across channels and frequency bands for one example session in each of two example subjects 1022 and 1069. These example sessions represent sessions with weak (figure 2, left) and strong (figure 2, right) stimulation effects. As a reference, the upper panels show the 'effect' during baseline, that is, the background fluctuations of band power. The lower panels show the stimulation effect in terms of band power changes, based on multiple pre- and post-stimulation pairs (see example inset panels on the right and figure 1). Notice that, even in the example subject 1069, where some strong stimulation effects are seen, these are restricted to a handful of channels and specific frequency bands. This observation is typical for all the sessions that exhibited a strong effect. Similarly, the example session on the left is a typical example of all the sessions that have a stimulation effect that is indistinguishable from the baseline fluctuations.

Zoom In Zoom Out Reset image size

The lower panels in figure 2 show the spatial layout of the iEEG stimulation and recording channels in the brain, with electrodes colour-coded by their corresponding stimulation effect sizes. Note that a strong stimulation effect, in this case on theta band, is not limited to contacts close to the stimulation site but also affected remote contacts (lower right panel).

In order to assess if the effect of stimulation exceeded baseline fluctuations in general across all 165 sessions and 79 patients, we compared the extrema of the stimulation effect to the extrema of the baseline fluctuations for each frequency band. Figure 3(A) shows the distributions of minima and maxima effect on theta band for baseline and stimulation. These extrema were taken across channels to capture the strongest effect during a session. Generally, it is evident that, even the channel with the strongest stimulation effect does not have a substantially larger effect size compared to the baseline fluctuations. In theta band, only 10.2% of the sessions exhibit a minimum (negative) stimulation effect that exceeds the adjacent value of the baseline minima. Similarly, only 18.1% of the sessions exhibit a maximum (positive) stimulation effect that exceeds the adjacent value of the baseline maxima (see figure 3(A)).

Zoom In Zoom Out Reset image size

The limited stimulation effect across all sessions was also evident when we computed the paired differences in effect between stimulation and baseline conditions. The histograms for the effect minima and maxima in figure 3(B) indicate that, in most sessions, even the most extreme effect sizes do not exceed the band power fluctuations during baseline. However, these distributions are not zero-centred (paired t-test for minima: p = 2.4 · 10⁻⁵, effect size for minima: −0.430; paired t-test for maxima : p = 4.9 · 10⁻⁵, effect size for maxima: 0.437), indicating that across patients and sessions there is a small but significant difference between baseline and stimulation conditions in our dataset. Similar results were found for all frequency bands (see figure S2 in Supplementary Material).

Next, we investigated whether the low stimulation effect size in most sessions can be attributed to the stimulation amplitude of the session. Figure 4(A) shows that there is no correlation between the effect size achieved in the session and the session stimulation amplitude. The distributions of effect sizes for each session, across all channels and frequency bands, are represented by their minima and maxima. Neither of these two measures tend to increase or decrease with the stimulation amplitude (range: 0.25–3.5 mA; see also figure S3 in Supplementary Material for band specific results).

Zoom In Zoom Out Reset image size

Furthermore, we considered all the pairs of stimulation sessions with the same stimulation location in the same subject (101 pairs). We tested whether their difference in effect size is correlated with the difference of stimulation amplitude between the sessions. Figure 4(B) shows that the absolute difference in effect size is not correlated with the absolute difference in stimulation amplitude. Thus, even for the same subject and the same stimulation location, an increase in stimulation amplitude does not necessarily produce a stronger effect.

Similar analysis was performed on the influence of stimulation frequency and duration. We found a significant anti-correlation between the stimulation frequency and the minimum values of effect across all sessions (Pearson's correlation, p = 0.008). The minimum effect becomes more negative with higher stimulation frequency, that is, stronger decrease in band power with higher frequencies of stimulation (specifically in theta, alpha, and beta bands, see figure S7 in Supplementary Material). We also found that the short stimulation duration (0.5 s) can elicit more extreme minimum and maximum effect values (see figure S8 in Supplementary Material).

To investigate the consistency of stimulation effect across sessions, we focused on pairs of sessions in the same subject and stimulation location. Figure 5 shows two examples of these session pairs. The example on the left (subject 1022) does not show positive correlation between the two sessions in terms of stimulation effect in different channels and frequency bands. The example on the right (subject 1069) shows that the patient's two sessions are positively correlated. Notice that this correlation is mainly driven by channels that exhibit a strong positive stimulation effect in the first place.

Zoom In Zoom Out Reset image size

In order to assess the level of consistency in stimulation effect across all 101 session pairs, we computed the consistency curve for each pair. Figure 6(A) shows the consistency curves of the two session pair examples in figure 5 alongside some illustrations on how the curve is computed: a circle of exclusion emanating from (0,0) is gradually enlarged and the consistency coefficient is calculated for varying values of the circle's radius (see Methods). The consistency curve (as a function of the radius) captures the consistency coefficient of the session pair when all effect values are considered (at radius 0) but also while increasingly excluding channel and frequency band combinations with weaker stimulation effect (at higher radii). The exclusion of channel/band combinations with weak effect serves to minimise the influence of the inherently inconsistent band power fluctuations on the consistency calculation. In addition, considering the selective connectivity of brain areas, it is expected that only a subset of channels will respond to a localised stimulation. The two curves shown in figure 6(A) capture the difference between high and low consistency as shown in the scatter plots of figure 5, not only when all values are considered, but also when only the strong effect values are considered. This approach of gradually excluding the weaker stimulation effects (around the level of baseline fluctuations) essentially allows us to capture consistency in the few channels that display a discernible stimulation effect in the first place.

Zoom In Zoom Out Reset image size

Figure 6(B) shows the consistency curves for all 101 session pairs and the confidence interval of consistency coefficients of the baseline periods (blue background). The overall consistency in this dataset is limited: 32.7% of the session pairs have higher consistency than the 97.5th percentile of the baseline 'effect' at radius = 0; 12.9% of the session pairs have higher consistency than the 97.5th percentile of the baseline 'effect' at radius = 3; and only 34.6% of the session pairs have a maximum consistency that is higher than the maximum value of the baseline 'effect' confidence interval. Four examples of maximum consistency coefficients on four of these curves are indicated with brown markers in figure 6(B). We will consider these maximum consistency coefficients as the representative values of the session pair consistency in the following (i.e. highest consistency achieved after exclusion of some not stimulation-related channels).

In figure 6(C), we demonstrate a strong and significant correlation between the average maximum effect and the maximum consistency coefficients (Pearson's $r = 0.536, p = 7.4\times10^{-9}$). Theta and alpha bands contribute more to this correlation (see figure S4 in Supplementary Material). As a comparison, we applied the same procedure to simulated data (normal distribution with mean 0 and standard deviation matching the the sessions' baseline), and the correlation is not present (Pearson's $r = -0.003, p = 0.438$). Essentially, the stronger stimulation effects also tend to be more consistent across sessions.

Finally, we built a multiple linear regression model to explain the maximum consistency coefficients as a linear combination of multiple independent variables including the average maximum effect (R² = 0.406, Adjusted R² = 0.340). The high explanatory power of the average maximum effect is also evident after running ANOVA on the multiple linear regression model, with the results shown in figure 6(D) (distributions produced after 200 bootstrap samples). Other than the strong effect of the average maximum effect on consistency, figure 6(D) shows a fair effect of both the task difference and the difference of baseline mean on consistency (p = 0.008 and p = 0.017, respectively), which are both anti-correlated with the maximum consistency coefficient.