Multivariate autoregressive models with exogenous inputs for intracerebral responses to direct electrical stimulation of the human brain

Abstract

A multivariate autoregressive (MVAR) model with exogenous inputs (MVARX) is developed for describing the cortical interactions excited by direct electrical current stimulation of the cortex. Current stimulation is challenging to model because it excites neurons in multiple locations both near and distant to the stimulation site. The approach presented here models these effects using an exogenous input that is passed through a bank of filters, one for each channel. The filtered input and a random input excite a MVAR system describing the interactions between cortical activity at the recording sites. The exogenous input filter coefficients, the autoregressive coefficients, and random input characteristics are estimated from the measured activity due to current stimulation. The effectiveness of the approach is demonstrated using intracranial recordings from three surgical epilepsy patients. We evaluate models for wakefulness and NREM sleep in these patients with two stimulation levels in one patient and two stimulation sites in another resulting in a total of 10 datasets. Excellent agreement between measured and model-predicted evoked responses is obtained across all datasets. Furthermore, one-step prediction is used to show that the model also describes dynamics in pre-stimulus and evoked recordings. We also compare integrated information-a measure of intracortical communication thought to reflect the capacity for consciousness-associated with the network model in wakefulness and sleep. As predicted, higher information integration is found in wakefulness than in sleep for all five cases.

Keywords: MVARX model; cross-validation; evoked response; integrated information; intracerebal EEG.

Figure 1

Recording and stimulation electrode placements for the subjects. Black dots represents recording channels while black “X” represents stimulating channel(s). (A) Subject A, right hemisphere is shown. 1, Inferior frontal opercular; 2, anterior horizontal lateral fissure; 3, middle frontal gyrus; 4, middle frontal sulcus; 5, superior temporal sulcus; 6, inferior frontal sulcus; 7, middle temporal gyrus; 8, middle frontal gyrus; 9, middle temporal gyrus; 10, orbital gyrus; 11, precentral gyrus; 12, superior frontal sulcus; and X, middle frontal gyrus. (B) Subject B, right hemisphere is shown. 1, Inferior frontal gyrus; 2, superior temporal sulcus; 3, posterior lateral fissure; 4, postcentral solcus; 5, superior temporal gyrus; 6, transversal temporal sulcus; 7, superior frontal gyrus; 8, subcentral gyrus; X (L1), precentral gyrus, and X (L2), subcentral sulcus. (C) Subject C, right hemisphere is shown. 1, Precentral gyrus; 2, posterior middle temporal gyrus; 3, inferior parietal lobule; 4, postcentral gyrus; 5, postcentral sulcus; 6, angular gyrus; 7, supramarginal gyrus; 8, anterior middle temporal gyrus; 9, inferior temporal gyrus, and X, superior parietal lobule.

Figure 2

Tukey-windowed median filtering for eliminating volume conduction artifacts. The upper trace depicts an example of raw data (blue solid line) and the Tukey-windowed median filter output (red dashed line). The lower trace depicts the weighting applied to the raw data (blue solid line) and the median filtered data (red solid line) to eliminate the volume conduction artifact.

Figure 3

Schematic diagram of the MVARX model. y_{i, n} denotes the recorded signals at the electrodes while x_n represents current stimulation and w_{i, n} is model error, or equivalently, a random input that generates spontaneous activity. a_{i, j} captures the a priori unknown connectivity between recording sites while [B]_i,: represents the a priori unknown transmission characteristics between the stimulus and recording sites.

Figure 4

Comparison between average CV evoked and average CV model responses of Subject A to two different stimulation strengths in wakefulness. In panels (A) and (B) the black dotted lines indicate the origin while the error bars denote the standard error of the mean. (A) Average CV evoked and average CV model responses of channels 1, 7, 4, and 11 with 1 mA current stimulation. (B) Average CV evoked and average CV model responses of channels 1, 7, 4, and 11 with 5 mA current stimulation. (C) Normalized mean-squared difference in each channel for 1 mA stimulation. (D) Relative root mean-squared energy in each channel for 1 mA stimulation. (E) Normalized mean-squared difference in each channel for 5 mA stimulation. (F) Relative root mean-squared energy in each channel for 5 mA stimulation.

Figure 5

Comparison between average CV evoked and average CV model responses of Subject A to two different stimulation strengths in sleep. In panels (A) and (B) the black dotted lines indicate the origin while the error bars denote the standard error of the mean. (A) Average CV evoked and average CV model responses of channels 1, 7, 4, and 11 with 1 mA current stimulation. (B) Average CV evoked and average CV model responses of channels 1, 7, 4, and 11 with 5 mA current stimulation. (C) Normalized mean-squared difference in each channel for 1 mA stimulation. (D) Relative root mean-squared energy in each channel for 1 mA stimulation. (E) Normalized mean-squared difference in each channel for 5 mA stimulation. (F) Relative root mean-squared energy in each channel for 5 mA stimulation.

Figure 6

Comparison between average CV evoked responses and average CV model responses of Subject B with two different stimulating locations in wakefulness. In panels (A) and (B) the black dotted lines indicate the origin while the error bars denote the standard error of the mean. (A) Average CV evoked and average CV model responses of channels 1, 3, 6, and 8 when the stimulating channel is L1. (B) Average CV evoked and average CV model responses of channels 1, 3, 6, and 8 when the stimulating channel is L2. (C) Normalized mean-squared difference in each channel when the stimulating channel is L1. (D) Relative root mean-squared energy in each channel when the stimulating channel is L2. (E) Normalized mean-squared difference in each channel when the stimulating channel is L1. (F) Relative root mean-squared energy in each channel with the stimulating channel is L2.

Figure 7

Normalized mean-squared response difference [see Equation (24)] in each dataset.

Figure 8

Comparison between recorded signal and one-step prediction of Subject B when the stimulating site is L2. 1.5 s pre-stimulus is shown followed by two-and-a-half epochs of evoked data. The black dotted lines in the panels indicate the origin. The model is estimated from data beginning with the fourth epoch. (A) Wake recorded and predicted signals in channels 1, 3, 6, and 8 ordered from top to bottom. (B) Non-REM sleep recorded and predicted signals in channels 1, 3, 6, and 8 ordered from top to bottom.

Figure 9

Normalized mean-squared one-step prediction error [see Equation (25)] in each dataset.

Figure 10

Exogenous input filters B for each channel as a function of time. The identical colormap is used for each row. (A) Subject A wake, 1 mA. (B) Subject A sleep, 1 mA. (C) Subject A wake, 5 mA. (D) Subject A sleep, 5 mA. (E) Subject B wake, stimulation site L1. (F) Subject B sleep, stimulation site L1. (G) Subject B wake, stimulation site L2. (H) Subject B sleep, stimulation site L2. (I) Subject C wake. (J) Subject C sleep.

Figure 11

Integrated information of Subject A, when the stimulation current is of 5 mA.

Figure 12

(A) Average maximum values of integrated information with error bars indicating one standard error. (B) Average lag at which maximum integrated information is achieved, with error bars indicating one standard error.