Novel visualization of the spatiotemporal relationship between ictal spiking and LFP supports the involvement of mid-range excitatory circuits during human focal seizures

Calculation of the st-SCA in MEA recordings

To apply the st-SCA to MEA recordings, we must account for the irregular timing and location of spiking activity across the array. A simplified analogy of this approach is to visualize spiking activity as stones being tossed into water. Consider throwing a single stone into water and analyzing the consequent effects by observing the resulting water ripples. We can simulate multiple sources by dropping identical stones from the same height but at different times and locations across the horizontal plane of the water surface, resulting in a complex landscape. To determine the contribution of a single stone to this landscape, we can take a field of view centered around individual stones. According to Eq. 3, averaging across all stones gives us the spatial pattern of activity associated with each stone. If we also include the time interval around each dropped stone, we obtain the stone’s characteristic spatiotemporal perturbation.

To apply this to the analysis of MEA recordings, spikes are detected for each channel in the MEA (Fig. 1B, left column), and the low frequency LFP associated with each spike is determined (Fig. 1B, middle column). This LFP is then spatially translated such that the associated spike position (x, y) is at the origin of a new axes (ξ, ψ) (Fig 1B, right column). This spike detection and LFP translation process is then applied to all channels. Averaging the results across all channels results in a field of view of the spike-associated LFP that is (1) centered around individual spikes and (2) approximately four times larger than the area of the MEA (Fig. 1B, right panel). This field is then calculated for time points τ to result in the st-SCA.

Application of st-SCA to clinical recordings allows visualization of spatiotemporal dynamics of human focal seizures

We applied the st-SCA method to analyze microelectrode recordings of focal seizures in patients undergoing epilepsy surgery evaluation. These recordings were obtained from 96-channel, 4×4mm MEA during and around ictal activity (Schevon et al., 2012; Truccolo et al., 2011). Filtering was used to extract local multi-unit neural firing activity and the associated low frequency component of the LFP of the surrounding network (Methods). A total of 19 seizures from seven patients was used for the study (Table S1). Both ictal and interictal recordings were evaluated, where interictal was defined as being at least two hours away from any known ictal activity. Of the seven patients analyzed for this study, five patients had arrays implanted in the recruited seizure territory (Patients 1-5), and two patients had arrays implanted in unrecruited seizure territory but within the clinically-determined seizure onset zone (Patients 6-7) (Schevon et al., 2012). As previously described, recruited seizure territory is defined as local tissue invasion by the seizure wavefront. Unrecruited territory is located outside of the recruited territory but can show rhythmic EEG activity due to local synaptic activity (Merricks et al., 2021; Schevon et al., 2012).

The representative STAs calculated across ictal and interictal states for two recruited territory recordings and one unrecruited territory recording are depicted in Fig. 2. The black lines represent the STAs, the red lines represent the associated noise estimates, and the vertical dotted lines indicate t = 0, i.e., the timing of the spike trigger. For both recruited and unrecruited territory recordings, the amplitudes of the ictal STAs (Fig. 2C, E, G) were larger than the corresponding the interictal STAs (Fig. 2D, F, H).

• Download figure
• Open in new tab

Figure 2.

Spike-triggered averages (STAs) in recruited and unrecruited cortical territories during ictal and interictal phases show different patterns. The black traces are the signals, and the red traces represent the associated noise estimates. Vertical stippled lines represent the zero of the time-axis.

A—B) Example signal trace of average ictal and interictal LFP activity across MEA channels.

C—F) The STA in the recruited territories show an evolution towards a characteristic negative peak, with or without surrounding oscillations, during the ictal phase. The ictal phase amplitudes are also much higher than those of the interictal phase.

G—H) The STA in the unrecruited territory show much smaller, non-zero ictal signals than the ones in the recruited territories. In contrast to the recruited territories, the small but dominant ictal peak polarity is positive in this case.

The insets in panels D, F, G, H show that unrecruited ictal and across all interictal phases show smaller amplitude, albeit non-zero STAs.

The amplitude for the unrecruited ictal STA (Fig. 2G), however, was much smaller than the recruited ictal STAs (Fig. 2C, E).

Patients with the MEA located in recruited seizure territory showed STAs with different morphologies (Fig. S2A-E), but all had a dominant negative peak around t = 0. Consistent with previous findings, we found that the STA for Patients 1-3 resembled a sinc-function with a peak embedded in a weak oscillatory component (Fig. 2C, 3A, S2A-C) (Eissa et al., 2018). In contrast, the STA for Patients 4 and 5 did not resemble a sinc-function as Patient 4 showed a dominant peak embedded in a strong oscillation (Fig. S4D) while Patient 5 showed no oscillatory component (Fig. 2E, Fig. 3B, S4E). The STAs for Patients 6 and 7 with the MEA in unrecruited territory were weak with a smaller amplitude deflection around to t = 0 (Fig. 2G, Fig. S4F, G).

• Download figure
• Open in new tab

Figure 3.

Properties of the ictal spatiotemporal spike-centered average (st-SCA) function during a focal seizure show different patterns for two representative patients.

A—B) The temporal average is calculated by averaging the st-SCA over all spatial contributions (±3.6mm).

C—D) A 3D view (azimuth = 0°, elevation = 70°) of the 3D st-SCA summed over time τ =±35ms. The center (ξ, ψ = 0,0) is indicated by the red dot. The two concentric circles are drawn to indicate that the center is surrounded by two rings. Note the apparent radial symmetry of the st-SCA pattern in (C). Grayscale is in μV units.

We then evaluated the relationship between spiking activity and the LFP in the spatiotemporal domain by computing the st-SCA over the entire MEA (ξ, ψ) and times τ = ±1ms (Fig. 3C, D). This 2ms interval averaged across to yield a 2D spatial topography. In the ictal phase for Patients 1-3, we observed a centrally located trough surrounded by a pair of rings with apparent radial symmetry (Fig. 3C, Fig. S2A). The distance between the center and the region indicated by the inner circle was ∼1.5mm, and the distance between center and the region indicated by the outer circle was ∼2.5mm (Fig. 3C). In contrast, the ictal phase for Patients 4 and 5 showed a deep well of stronger negative activity (e.g., Fig. 3D). The st-SCAs during the interictal phase as well as the results obtained in unrecruited territories showed different patterns with relatively smaller amplitude signals (Fig. S2B, D, E, F). For a list of representative temporal and spatial STAs for each patient, see Fig. S2.

In sum, three of five patients with recordings from recruited territories (Patients 1-3) showed st-SCAs resembling sinc-functions in the temporal domain and st-SCAs resembling donut-shaped rings of activity in the spatial domain. In Patients 4 and 5, the STAs were characterized by a non-sinc morphology, and the st-SCAs showed deep and diffuse wells of negative activity (Fig. 3D, S4D, E).

Note that these observations were not attributable to widespread correlations among MEA electrodes. To demonstrate that the observed st-SCA patterns are representative of the spike-LFP relationship and not the global, macroscale correlations amongst network LFPs, we showed that STAs in unrecruited territory show a large and significant oscillatory component only when triggered by spikes from recruited territories (Fig. S4B), and not when triggered by spikes from unrecruited territories (Fig. S4C). This result is a replication of previous studies (Eissa et al., 2017). Furthermore, randomizing the spike times detected across the MEA resulted in complete destruction of the observed st-SCA patterns, emphasizing the importance of spike timing as the driver for these spatiotemporal patterns (Fig. S5).

Finally, calculation of the st-SCA after applying a spatial filter to decorrelate LFP signals across MEA channels did not qualitatively alter the st-SCA patterns (Methods, Fig. S6).

Theoretical model reveals a unique symmetry that allows for prediction of temporal and spatial components of the st-SCA

The st-SCA calculations presented above were obtained indirectly by an averaging procedure in order to reduce the noise not associated with the spike trigger in the spike-LFP relationship. In this section, we simulate an ideal noise-free, controlled environment where spatial and temporal relationships can be measured directly (Fig. 4).

• Download figure
• Open in new tab

Figure 4.

A mathematical model as represented by the recording of macroelectrode measures the underlying network’s st-SCA.

The electrode covers an area of one-dimensional cortex where we record the effect associated with a single centrally located ictal action potential, δ(r, τ) = δ(0,0), i.e., the macroelectrode measures the underlying network’s temporal component of the st-SCA, st-SCA₄(τ). This measurement can be approximated by an unknown action potential’s spatiotemporal cortical activation function, f(r, τ), integrated over the spatial range [−R, R] covered by the electrode.

Here, the sum of all measurements across the MEA can be represented by a single macroelectrode, which records cortical activity when a spike occurs at time zero and at the center of the electrode. This is a measurement of the noise-free spatiotemporal local field function f(r, τ) (space (r) and time (τ)) associated with the central spike represented by a delta function, δ(r, τ). Because the potential of cortical generators attenuates sharply with distance, we consider contributions from activity in areas not directly under the macroscopic electrode to be negligible. Under this assumption, the electrode’s signal can be approximated by summing the contributions over the neocortical area under the electrode, and we find the following expression for STA(τ): with rect(−R, R) representing a rectangular window bounded by [−R, R]. Similarly, if we compute the spatial component of the SCA, st-SCA_s(r), by integration over a fixed time epoch, [−T, T], around the seizure onset, we get:

Because the action potential can be represented by a unit impulse, STA(τ) and st-SCA_s(r) are equivalent to a cortical unit impulse response (UIR). The UIR can be used to link the spike train to the associated LFP (Supplementary Information). In most cases the function f(r, τ) cannot be simply derived from measuring the STA(τ) and st-SCA_s(r), but one could apply an ad hoc solution so that Eq. 4 and 5 produce the correct characterization:

In special cases, however, f(r, τ) may be mathematically derived. Note that we determined earlier that in three of five patients with recordings in the recruited territory, the STA(τ) resembles a sinc-function:

This represents a unique scenario in which the f(r, τ) and st-SCA_s(r) may be predicted from the characterization of STA(τ) alone. The sinc-function is the UIR of an ideal filter, and a previous study has used this property to show that an ictal spike train passed through an ideal filter recreates the seizure’s LFP (Eissa et al., 2018). Since the sinc-function is defined as the Fourier transform of a rectangular function, the relationship between time and space conveniently parallels a time-frequency Fourier-transform-pair. Accordingly, we find:

Substitution of this finding into Eq. 5 enables us to find the st-SCA_s(r), which represents the spatial postsynaptic effects of an action potential effective over a fixed time epoch around the seizure onset ([−T, T]):

Thus, we find that in the special case where both STA(τ) and st-SCA_s(r) are described by sinc-functions, the temporal features of the spike-associated LFP can predict the spatial features, and vice versa.

Quantification of the peak-to-peak distance of spatial patterns supports the involvement of mid-range excitatory connections

Since the model showed that a sinc-function in the temporal domain predicts a sinc-function in the spatial domain, we aimed to more quantitatively describe the donut-shaped activity observed in Fig. 3C. Taking advantage of the radial symmetry observed in the st-SCA, we converted the Cartesian coordinates (ξ, ψ) into polar coordinates (r, θ) and focused on the spatial relationship with respect to r (Fig. 5A). This enabled us to depict the st-SCA in two dimensions, (r, τ) (Fig. 5B). A detail of that relationship is depicted in Fig. 5C, and the summed values across this two-dimensional detail are plotted along its margins. These summed values are the two components as a function of space and time (r and τ). Note that the bottom graph in Fig. 5C represents the central trough (τ = ± 35ms) of the function shown in Fig. 3A. As anticipated by the outcome in Eq. 9 we observed a spatial component (Fig. 5C) that shows a central trough with smaller amplitude side lobes---a pattern consistent with the shape of a sinc-function. Note that the resolution and range of the spatial component (r = ±3.6mm) is limited by the size of the MEA (Methods, Fig. 1B). Consistent with the donut-shaped rings observed in Fig. 3C, the peaks of the function shown in the side panel of Fig. 5C were separated by ∼2.5mm (blue arrows, Fig. 5C).

• Download figure
• Open in new tab

Figure 5.

The method to compute the 2D spatiotemporal spike-centered average, st-SCA(r, τ) allows for a detailed visualization of spatial and temporal components for Patient 1

A) The Cartesian coordinates (ξ, ψ) from the 3D st-SCA are converted into polar coordinates (r, θ), resulting in a 2D st-SCA.

B) A color representation of st-SCA(r, τ). The top trace, the temporal component of the st-SCA_t(τ), is obtained by the sum of st-SCA(r, τ) over r (same as the signal in Fig. 3A). Amplitude and color scale are in μV.

C) Detail of the central part of Panel B. The left margin shows the resulting wave from summation over time, generating the spatial component of st-SCA. The green arrows on the left indicate the distance (∼2.5mm) between the peaks seen in this function. The bottom margin depicts the resulting wave from summation over space, generating the temporal component of the st-SCA.

Patients

Seven patients with pharmacoresistant focal epilepsy underwent chronic intracranial EEG studies to help identify the epileptogenic zone for subsequent removal. Patients 1, 4, 6, and 7 were recruited at Columbia University Medical Center, and Patients 2, 3, and 5 were recruited from Massachusetts General Hospital/Brigham and Women’s Hospitals (Table S1). Procedures were approved by the Internal Review Board committees at Columbia University Medical Center, The University of Chicago Comer Children’s Hospital, and Massachusetts General Hospital/Brigham and Women’s Hospitals. The patients’ surgeries and treatment plans were not directed by or altered as a result of these studies.

Signal acquisition and pre-processing

A 96-channel, 4 × 4mm MEA (Utah array; Blackrock Microsystems) was implanted along with subdural electrodes (ECoG) with the goal of recording from seizure onset sites. Additional details of study enrollment and surgical procedures have been previously published (Schevon et al., 2012; Truccolo et al., 2014). Signals from the MEA were acquired continuously at a sample rate of 30 kHz per channel (0.3-7500Hz bandpass, 16-bit precision, range ±8 mV). The reference was epidural. Up to three seizures from each patient were selected for detailed analysis to avoid biasing the dataset from the patients from whom many seizures were recorded. Seizure recordings were categorized as recruited or unrecruited territory using previously described methods (Schevon et al., 2012). Channels and time periods with excessive artifact or low signal-to-noise ratio were excluded.

Unit activity was identified using filtered 0.3-3kHz signals with spikes defined as deflections ≥4 standard deviations below the mean. The low frequency component of the local field potential (LFP) activity across the array was created by averaging the artifact-free LFP activity from all micro-electrode signals filtered 2-50Hz. The averaged LFP procedure has been shown to generate signals that are representative of and comparable to nearby electrocorticography signals (Eissa et al., 2017; Eissa et al., 2018).

Spatiotemporal spike-centered average (st-SCA) calculations and signal analysis

All signal processing and statistical analyses were performed in MATLAB (MATLAB, Natick, MA, USA). The st-SCA was determined as follows. First, as described in (Eissa et al., 2017; Eissa et al., 2018), we detected the spikes in the multi-unit activity. Next, we collected the spatiotemporal data around each spike and translated the time and position of all associated LFPs relative to each spike’s time and position (Fig. 1). Finally, we summed all translated data and computed the average at each time and position by dividing the sum by its number of contributions. Note that this position-dependent average is necessary because not every position receives the same number of contributions during the translation of the LFP’s axes.

Evidence of radial symmetry of the st-SCA (Fig. 3C) allowed conversion from Cartesian coordinates (ξ, ψ) coordinates to polar coordinates (r, θ). By ignoring the minor deviations from radial symmetry, we focused on the spatial component of the st-SCA with respect to r (Fig. 5A), which enabled us to depict the spatiotemporal properties in two dimensions (Fig. 5). Furthermore, if we compute the sum across space, we obtain purely the temporal component of the st-SCA, which is equivalent to the STA. Similarly, summation over time τ generates the spatial component of the st-SCA. With these results, we can assess to what extent our model of the ictal network, a linear time-invariant (LTI) system with unit impulse response C(τ) ∝ sinc(r, τ), fits the data. A stepwise description of this method can be found in the Supporting Information.

For calculations involving the spatial filtering of LFP signals, we applied the spatial whitening process as described in Hyvärinen et al. (2001) and Telenczuk et al. (2017). As previously published, a signal is spatially filtered by matrix multiplication with a whitening matrix W, where W is the inverse square root of the signal’s covariance matrix, C: where E is a matrix of eigenvectors of C_signal,and D is a diagonal matrix with inverse square roots of eigenvalues λ_ion its diagonal, such that and D_ij (Telenczuk et al., 2017). In this study, the signals being transformed were the MEA channel signals bandpass filtered at 2-50Hz.