- 1Department of Psychology, University of Oslo, Oslo, Norway
- 2RITMO Centre for Interdisciplinary Studies in Rhythm, Time and Motion, University of Oslo, Oslo, Norway
- 3Department of Neurosurgery, Oslo University Hospital, Oslo, Norway
- 4Department of Neuropsychology, Helgeland Hospital, Mosjøen, Norway
- 5Department of Psychology, Helen Wills Neuroscience Institute, University of California, Berkeley, Berkeley, CA, United States
Introduction: Intracranial electrodes are implanted in patients with drug-resistant epilepsy as part of their pre-surgical evaluation. This allows the investigation of normal and pathological brain functions with excellent spatial and temporal resolution. The spatial resolution relies on methods that precisely localize the implanted electrodes in the cerebral cortex, which is critical for drawing valid inferences about the anatomical localization of brain function. Multiple methods have been developed to localize the electrodes, mainly relying on pre-implantation MRI and post-implantation computer tomography (CT) images. However, they are hard to validate because there is no ground truth data to test them and there is no standard approach to systematically quantify their performance. In other words, their validation lacks standardization. Our work aimed to model intracranial electrode arrays and simulate realistic implantation scenarios, thereby providing localization algorithms with new ways to evaluate and optimize their performance.
Results: We implemented novel methods to model the coordinates of implanted grids, strips, and depth electrodes, as well as the CT artifacts produced by these. We successfully modeled realistic implantation scenarios, including different sizes, inter-electrode distances, and brain areas. In total, ∼3,300 grids and strips were fitted over the brain surface, and ∼850 depth electrode arrays penetrating the cortical tissue were modeled. Realistic CT artifacts were simulated at the electrode locations under 12 different noise levels. Altogether, ∼50,000 thresholded CT artifact arrays were simulated in these scenarios, and validated with real data from 17 patients regarding the coordinates’ spatial deformation, and the CT artifacts’ shape, intensity distribution, and noise level. Finally, we provide an example of how the simulation platform is used to characterize the performance of two cluster-based localization methods.
Conclusion: We successfully developed the first platform to model implanted intracranial grids, strips, and depth electrodes and realistically simulate thresholded CT artifacts and their noise. These methods provide a basis for developing more complex models, while simulations allow systematic evaluation of the performance of electrode localization techniques. The methods described in this article, and the results obtained from the simulations, are freely available via open repositories. A graphical user interface implementation is also accessible via the open-source iElectrodes toolbox.
Introduction
Intracranial subdural grids and depth electrodes are implanted in patients with drug-resistant epilepsy as part of their pre-surgical evaluation. Electrophysiological and neuroanatomical data are used to delineate the seizure onset zone and functional areas that will guide respective surgery (Rosenow and Lüders, 2001). Intracranial electroencephalography (iEEG) recordings provide insights into human brain electrophysiology and functional mapping with unparalleled spatial and temporal resolution, offering both clinical and research applications. Knowing the exact location of electrodes in relation to the individual cortical or subcortical anatomy is a prerequisite for a complete understanding of the electrophysiological data; leading to a precise resection of the epileptic foci and the anatomical localization of specific brain functions (Lachaux et al., 2003; Parvizi and Kastner, 2017; Frauscher et al., 2018; Stolk et al., 2018).
One of the most common approaches to define electrode coordinates is to localize their artifacts in post-implantation computer tomography (CT) images co-registered to pre-implantation magnetic resonance imaging (MRI) scans. Typically, electrode CT artifacts are extracted by thresholding CT images, and then electrode coordinates are manually or semi-automatically computed from these (Blenkmann et al., 2017; Stolk et al., 2018; Li et al., 2019; Centracchio et al., 2021; Rockhill et al., 2022). In both cases, errors associated with the procedure are rarely quantified. It is well known that other CT artifacts than the ones of interest (e.g., connection cables or clips), adjacent electrodes, overlapping grids or strips, and noise or low image resolution make precise localization problematic (Brang et al., 2016; LaPlante et al., 2016; Narizzano et al., 2017).
Other localization approaches are based on the use of post-implantation MRI and CT images (LaViolette et al., 2011b; Hinds et al., 2018), clinical neuronavigational data (Gupta et al., 2014), or electrophysiological data (Branco et al., 2018b). However, these data might not be available in some clinical setups.
Over the last few years, the spatial resolution of grids and depth electrodes has reached inter-electrode distances (IED) of 2–3 mm, and will likely improve even more (Gupta et al., 2014; Chang, 2015; Martin et al., 2018). Given their size and the low signal-to-noise ratio (SNR) of images, high-density (HD) arrays are harder to localize with currently available methods [see novel attempts in Hamilton et al. (2017), Narizzano et al. (2017), Branco et al. (2018a), and Erhardt et al. (2020)].
Although methods to localize CT artifacts and co-localize them to pre-implantation MRI are common, there is no reliable gold standard to quantify their precision and robustness against noise. Methods are usually validated with different datasets, where noise is not quantified. This lack of standardization in the validation makes the comparison of methods imprecise and leaves unanswered questions: what is the localization error associated with a particular SNR? Which localization method performs better for a particular case? These are important factors undermining the validity of claims regarding specific brain areas related to a function or pathology.
To aid in solving these problems, we propose a new platform to model realistic intracranial electrode implantation scenarios and their respective CT artifacts, enabling systematic quantification of errors in electrode localization algorithms.
In the next sections, we will introduce methods for fitting subdural grid and strip electrode models onto the smooth hull surface, and depth electrode array models targeting subcortical sites. We then describe the modeling of thresholded CT artifacts and their noise. Thereafter, we report the simulation results and their validation with real data. Moreover, we provide two examples: (i) performance evaluation of localization methods, and (ii) modeling implantations on a patient’s native anatomy. Finally, we discuss the results in the context of existing electrode localization methods, the limitations of our approach, and future challenges.
The methods and results presented in this article are publicly available (except for individual patient data).
Materials and methods
The proposed platform allows the systematic evaluation of intracranial electrode localization algorithms. The workflow consists of several steps. Concisely, they can be grouped into pre-processing of images to obtain a Smooth Cortical Envelope (SCE), the simulation of electrode coordinates, the simulation of CT artifacts, and the performance evaluation. Figure 1 shows the main steps, and the following subsections describe them in detail.
Figure 1. Workflow for simulating electrode CT artifacts and evaluating localization algorithms. The main steps of the workflow are shown in the diagram. First, a Smooth Cortical Envelope (SCE) is obtained from an individual or average atlas image. Electrode arrays are modeled according to design parameters and fitted over the SCE (grids) or placed within the brain tissue (depths) to obtain the electrode coordinates. Then, CT artifacts are simulated at each electrode coordinate, considering the electrode type and spatial orientation. Subsequently, noise is added to simulate different realistic scenarios. CT artifacts are then fed into the localization algorithm under evaluation. Finally, localized coordinates are compared with the ground truth electrode coordinates. Knowing the true electrode coordinates allows a correct estimation of the localization error and its sensitivity to different noise levels. Moreover, the results can be easily replicated. Orange parallelograms represent data, and purple squares represent processes.
Unless otherwise specified, the methods were implemented in Matlab R2019 (The MathWorks Inc., USA) and iElectrodes toolbox for Matlab (Blenkmann et al., 2017). The most relevant functions needed to perform the simulations are mentioned along the method’s subsections.
In general, median values were used to describe the central tendencies of the variables of interest since most distributions were not normal. Their corresponding CIs were obtained by bootstrapping.
Smooth cortical envelope surface and seed points
The implantation of real grids and strips onto different cortical areas introduces variations in how these arrays bend to follow the brain curvature. To reproduce realistic scenarios, we simulated grids implanted over 3D cortical surfaces. High-resolution pial surfaces extracted from the MNI atlas were used as a starting point (or individual structural MRI images for simulations in native space; see subsection “Example II: Interactive modeling of implanted electrodes in patient’s native brain space” for details). In our processing pipeline, we used Freesurfer software (Dale et al., 1999), and individual cortical parcellation images were obtained using the Destrieux atlas (Destrieux et al., 2010), but other atlases or software could also be used. We then computed an SCE surface for each cerebral hemisphere by enclosing the corresponding pial surface with a 30 mm radius sphere. A mesh smoothing was applied to remove small local protuberances (low-pass spatial filter, 100 iterations, alpha weight = 0.5, Iso2Mesh toolbox, Fang and Boas, 2009). Supplementary Figure 1 shows an example of the pial surface and SCE computed for an individual patient’s brain. The procedure is implemented in the buildSCE.m function.
We visually selected 57 seed points over the SCE surface and used them as reference points to model grids, strips, and depth electrode arrays. We computed the local curvature of a smoothed version of the SCE (Rusinkiewicz, 2004). Seed points were visually selected on regions where the local curvature was relatively “Low,” “Medium,” or “High” within the range of curvature values. We modeled different electrode arrays depending on the location and the local curvature of the SCE surface.
Modeling grid and strip electrode coordinates
The following steps are performed for modeling subdural grids and strips:
(1) A 2D flat model of the array is generated and placed tangentially to the SCE in a given “seed” point (e.g., the purple grid in Figure 2D) using the function tangentialGrid.m. For strip cases, the number of rows is set to three to avoid unrealistic geometrical deformations (e.g., undulating shapes).
Figure 2. Modeling grid and strip electrode coordinates over the SCE surface. Electrode arrays are initially modeled as 2D (grids) or 1D (strips) objects, and then fitted to the SCE surface by minimizing the translation (Et) and deformation (Ed) energies (Equation 1). (A) First, diagonal, and second neighboring connections are used to compute the deformation energy Ed, accounting for normal, shear, and bending, respectively. (B) The translation energy Eti associated with electrode i is proportional to the distance between the initial (x0i) and the final location (xi) over the SCE surface. (C) The deformation energy Edij associated with neighboring electrodes i and j is proportional to the amount of deformation between the original (d0ij) and final distance (dij). Equation details are explained in Subsection “Modeling grid and strip electrode coordinates”. (D) Example of a 2D grid model (4 × 8, 10 mm IED) before (pink) and after (green) being fitted to the SCE surface on the left frontal lobe. (E) Zoom-in detail showing the normal vectors at each electrode location. Normal vectors are subsequently used to orient the simulated CT artifacts. IED, inter-electrode distance; SCE, smooth cortical envelope.
(2) The 2D model is fitted over the SCE surface using an energy minimization algorithm (Figure 2D, Dykstra et al., 2012; Trotta et al., 2018) implemented in the function projection2mesh.m. The energy function
is minimized by varying the electrode coordinates to the final position xi, where Et is the translation energy, Ed is the deformation energy, and K is a constant value. The energy minimization is constrained to the electrode coordinates xi being closely located over the SCE surface, i.e., ∀i, ||xi−si|| ≤ ϵ, where si is the closest point in the SCE surface to electrode i, and ε is a tolerance distance.
Et (Figure 2B) represents the energy needed to translate the electrode coordinates from the original position of each electrode to the final position xi on the SCE surface:
where N is the number of electrodes.
Ed (Figure 2C) represents the energy required to change the inter-electrode euclidean distance between neighbors from the original to the final location. It is defined as
where dij is the distance between electrodes, is the initial distance between contacts, and αij takes values of 1 or 0 if electrodes i and j are neighbors or not, respectively. To control for normal, bending, and shear deformations, first, second, and diagonal (grids only) neighbors are considered, as shown in Figure 2A (Trotta et al., 2018). Ed is typically interpreted as the deformation of “springs” connecting the electrodes (Dykstra et al., 2012).
(3) A normal vector at each electrode coordinate is computed. We select the coordinates from each electrode of interest and the nearest neighbors and compute a principal component analysis (PCA) of these points. The normal vector is the one associated with the smallest component (Figure 2E).
(4) For strips, only the middle row coordinates are kept from the three rows, while the lateral ones are discarded.
We simulated electrode arrays of multiple dimensions representing models frequently available in the market by different manufacturers. These cover various combinations of size and IED, including grids of 2 × 4, 4 × 4, 4 × 8, 8 × 8, 8 × 16, and 16 × 16 contacts, strips of 1 × 4, 1 × 6, and 1 × 8 contacts. Arrays were simulated at the 57 seed points in the standardized MNI SCE surface with different local curvatures (Supplementary Figure 2).
Grid and strip implantation scenarios were simulated several times at each seed point, rotating the grids at angles multiple of 30 degrees around the center. The location and the local curvature of the SCE surface surrounding the seed points determined which arrays were simulated. Large grids were simulated over Low curvature regions, whereas medium and small size grids were simulated over Medium and High curvature areas. For example, over the lateral fronto-temporo-parietal cortex, it is realistic to simulate an 8 × 8, 10 mm IED grid, but unrealistic (in usual clinical settings) to simulate the implantation of such a big grid over the frontal pole.
For the fitting of grid and strips onto the SCE surface, coefficient K was set to 1,000, and the tolerance distance ε = 0.1 mm for all arrays, except for the 16 × 16 grid cases where K = 100 and ε = 0.5 mm (Equation 1).
Deformations are needed to fit a grid or strip (originally plane objects) over the brain surface (a curved object). When localizing grids, we typically observe that these deformations are small. Therefore, we measured the IED after fitting grids and strips to the SCE surfaces and discarded simulations if one or more neighboring electrodes had IED variations over 5% from the original values. A simple example of the simulation is implemented in the simulateArraySCE.m function.
Modeling depth electrode coordinates
In clinical practice, the implantation of depth electrodes is defined by two points, and therefore a unique trajectory connecting them. The points are typically defined as an “entry point” on the cortical surface and a “target point” at the deepest brain location reached. We will adopt this nomenclature throughout this paper.
It is common practice to use trajectories orthogonal to the skull surface to avoid sliding of the drill and minimize bone damage during surgery. Therefore, given an entry point on the SCE surface, we define the array trajectory vector orthogonal to this surface. Two strategies are applied to distribute the electrodes along the trajectory uniformly:
(1) If the array is shorter than the amount of tissue intersected by the trajectory, the target point is defined as a random point within the trajectory while keeping all contacts inside the brain.
(2) If the array is longer than the intersected tissue, the target point is defined as the deepest electrode point in the trajectory within the brain tissue (enclosed by the SCE surface). Contacts outside the brain tissue are removed.
As an optional step, a curvature deformation can be applied within a random orientation 2D plane containing the principal axes of the array. Using the symmetric Lanczos window, an arc-shaped electrode array is obtained. A simple example of the simulations is implemented in simulateDepthArray.m function.
We selected 33 entry points over the MNI SCE from the previous set of 57 seed points, restricted to locations where the implantation of depth electrodes were realistic. At each point, we simulated the implantation of depth electrodes as combinations of 4, 8, 10, 12, 15, and 18 contacts; 3, 5, and 10 mm of IED; and linear or curved deformation (maximum deformation of 1% of the total length).
Modeling overlapping grids and strips coordinates
Overlapping grids or strips are obstacles for methods that aim to detect intracranial electrodes automatically (e.g., Figure 3F). We developed a method to realistically simulate these cases (Figure 3E), implemented in the function simulateOverlaps.m. Briefly, it consists of:
Figure 3. Modeling overlapping grids. Example of an overlapping grid in an 8 × 8, 10 mm, IED case. (A) Original electrode coordinates (green) were transformed into the principal component space, where a surface (Sfit) was fitted to the coordinates. (B) In a 2D space defined by PC1 and PC2, the overlapping grid (magenta) was defined with a random orientation. (C) The overlapping grid was translated (in the arrow direction) until the desired number of overlapping electrodes was achieved. (D) Coordinates were back-projected to the original 3D space surface. (E) Voxel artifacts were simulated in the new overlapping grid electrode coordinates, overlaying the SCE surface. Voxel artifacts’ orientations were interpolated from close-by electrodes of the original array. (F) Example of real overlapping grids (5 × 4, 10 mm IED and 4 × 12, 6 mm IED) overlaying the left SCE surface. The inset shows the corresponding thresholded CT artifact voxels. Plots (B,C) show the PC1 axes direction inverted for illustrative purposes. IED, inter-electrode distance; SCE, smooth cortical envelope; A, anterior; S, superior.
(1) A surface is fitted to the background array (Figure 3A).
(2) The coordinates are projected to a 2D principal component space (PCA) where the overlapping array is initially positioned (Figure 3B).
(3) The overlapping array is displaced in a random direction until the desired amount of overlap is achieved (Figure 3C).
(4) The coordinates are back-projected to the initial surface (Figure 3D).
A more detailed explanation of the procedure can be found in Supplementary Section 3.
Grids and strips were simulated with a 10% overlap and no overlap. Overlapping grids were defined with a random orientation and translated in a random direction. In principle, grids and strips of any size and IED can be overlapped. However, for the sake of simplicity, we used arrays with the same dimensions.
Modeling electrode computer tomography artifacts and their noise
Grid and depth arrays have a flat disc and cylinder-shaped metallic contacts, respectively. Their corresponding thresholded CT artifacts usually show ellipsoidal shapes, and their voxels’ intensity typically decreases with their distance to the center of the electrode (Blenkmann et al., 2017; Stolk et al., 2018). Therefore, we modeled the CT artifacts representing electrodes as ellipsoidal-shaped clusters of voxels of varying intensity. In a parsimonious approach, we defined the intensity as a linear function that declined with the scaled radial distance r of a given location (x, y, z) to the ellipsoid center as:
where
and a, b, and c are the semi-axes of the ellipsoid model. Semi-axes length parameters were defined depending on electrode type, their inter-electrode distance, following manufacturers’ specified dimensions, and visually adjusted to match real CT artifacts’ observed size and shape. Specifically, we used a = 2.2, b = 2.2, and c = 1.5 for 10 and 5 mm IED grids, a = 1.1, b = 1.1, and c = 1 for 3 mm IED grids, a = 1.25, b = 1.25, and c = 1.75 for 10 and 5 mm depth electrode arrays, and a = 1.1, b = 1.1, and c = 1.5 for 3 mm IED depth electrode arrays.
The function simulateVoxels.m was used to simulate the clusters of voxels for each CT artifact. First, voxel coordinates were generated by using a 0.5 mm resolution 3D square lattice grid of 10 mm × 10 mm × 10 mm around each electrode to represent the effect of discrete sampling in CT images. All sampling lattice grids had a common random origin and orientation mimicking the random location and orientation of the patient’s head in the image reference system. Then, we centered the reference coordinate system at the electrode’s coordinate. We rotated the reference system to align its z-axis to the electrode’s specific normal vector, i.e., the surface’s orthogonal vector for grids and strips (Figure 2E) or the trajectory vector for depths. Finally, we computed the intensity of each voxel using Equation 4.
The localization of intracranial electrodes is typically sensitive to the signal-to-noise ratio of the processed CT images. To model the noise in the CT artifacts, we added spatially correlated and uncorrelated noise to each voxel’s intensity using the function addnoise2electrodes.m (Britten et al., 2004; Gravel et al., 2004). Spatially correlated noise was computed by convolving uncorrelated noise from a distribution with a kernel provided by Britten et al. (2004). Uncorrelated noise values were drawn from a distribution. Therefore, the total noise variance was . The two noise sources contributed differently to the deformations observed in the CT artifacts.
Finally, voxels with intensity values below zero were removed after adding noise.
We simulated the CT artifacts at each electrode location and manipulated the noise levels affecting these. Twelve noise levels were simulated for each scenario, where was logarithmically distributed between 0.2 and 2.2 (from now on denoted as noise levels 1–12). The ratio between correlated and uncorrelated noise variances was set to 20:1, a heuristically determined ratio obtained by observing the shape of CT artifacts. The individual impact of these two noise sources was clear at high noise levels. Higher ratios of uncorrelated noise left many “disconnected” voxels surrounding the CT artifacts, whereas higher ratios of correlated noise introduced unrealistically deformed shapes.
Validation
Patients
To validate our modeling platform, we compared our simulations with real data. We analyzed MRI and CT images from 17 adult patients with drug-resistant epilepsy who underwent iEEG recording as part of the pre-surgical evaluation for respective surgery.
Three patients were implanted with only subdural grids, seven with subdural grids and depth electrodes, and seven with only depth electrodes. IED ranged from 3 to 10 mm for both grids and depth arrays. Patients implanted with depth electrode arrays had on average 69 depth contacts (range 8–172 contacts per patient, total of 1,068 contacts, DIXI Medical, France, Ad-Tech Medical Instrument Corporation, USA, or PMT Corp., USA). Patients implanted with subdural grids had on average 146 subdural contacts (range 29–306 contacts per patient, total of 1,656 contacts, Ad-Tech Medical Instrument Corporation, USA, or PMT Corp., USA).
Electrode localization in real cases
We followed a routine procedure to localize intracranial electrodes (Stolk et al., 2018; Blenkmann et al., 2019). Pre-implantation T1-weighted MRI images were processed using the FreeSurfer standard pipeline (Dale et al., 1999), where individual brain segmentation images, pial surfaces, and cortical parcellation images (Destrieux atlas) were obtained (Destrieux et al., 2010). Post-implantation CT images were co-registered to the pre-implantation MRI using SPM 12 software (Studholme et al., 1999). Realigned MRI and CT images were resampled to 0.5 mm × 0.5 mm × 0.5 mm resolution. CT images were thresholded to visualize the clusters of high-intensity voxels (also known as CT artifacts). Threshold values were visually defined to identify clusters of voxels representing individual electrodes following the procedure described in Blenkmann et al. (2017).
Each cluster of high-intensity voxels was extracted using the k-medoids clustering algorithm (as implemented in iElectrodes toolbox). Then, the weighted center coordinates of each cluster were computed. In noisy situations, the center of each cluster was visually identified instead of using automatic clustering.
Validation of electrode coordinates
To validate the modeling of electrode array coordinates, we visually inspected the resulting coordinates. Grid and strip coordinates were expected to follow the SCE surface, while only minimally deforming the IED. Moreover, we described deformations in the distance between their first, second, and diagonal neighbors. To integrate the impact of these deformations in the framework of electrode localizations, we qualitatively compared these deformations with those obtained from manually localized real arrays of 10 mm IED. We anticipated simulation deformations to be smaller than the localized ones, since the latter includes array deformations and localization errors.
In addition, we evaluated the deformations of depth arrays. First, array lengths were normalized to allow their comparison. Then, we defined the main axis connecting the first and last contact. Finally, we measured the shortest distance between each contact and the axis.
Validation of computer tomography artifacts
To validate our model of electrode CT artifacts, we compared the characteristics of these synthetic artifacts with those of real artifacts. Noise characterization was particularly relevant since it can affect the accuracy of localization algorithms. We applied an information theory measure, entropy, to delineate the noise level of simulated and real data. With this aim, we first computed the principal axis of each electrode. For grids, we computed the orthogonal direction at each electrode given its closest neighbors, and for depth electrodes, we computed the principal axis as the direction connecting the first and last electrodes of the array. We aligned the principal axes to the z-axis and positioned their centers at the coordinate system’s origin. Then, the 2-dimensional intensity-radius histograms were computed, showing the bivariate distribution of voxels in terms of intensity (on the y-axis) and radial distance to the center (on the x-axis). Histograms were computed in 0.25 mm bins between 0 and 5 for the radial distance, and 0.05 bin size for the normalized intensity values between 0 and 1, while excluding intensity values below percentile 2.5 and above percentile 97.5. Afterward, we computed the pairwise 2D cross-correlation matrices between all histograms within an array. A cross-correlation matrix indicates the degree of similarity between two intensity-radius histograms as a function of the displacement of one relative to the other. Sharp cross-correlations are indicative of similar histograms, whereas smother (more distributed) cross-correlations indicate less consistency between the histograms. In our case, we expected the noise to differently affect each CT artifact within an array. Therefore, this should be reflected in smother cross-correlations of their corresponding 2D histograms when noise levels are high, and sharp and peakier cross-correlations when noise levels are low. Finally, we estimated the array’s noise by measuring the entropy of these cross-correlation images and computing their average. Entropy values represent the amount of disorder present in the cross-correlation matrices. Entropy has been previously used to estimate the focus of images, by evaluating the uniformity of the 2D spectral decomposition of images (Bove, 1993; Kristan and Pernu, 2004). Similarly, in our case, entropy values will increase as the cross-correlation images get close to uniformly distributed (e.g., smoother images) and decrease as they get farther from uniform (e.g., peakier images). Therefore, low entropy values indicate more consistent, i.e., less noisy, artifacts within an array, and vice versa. Entropy values were sensible to the size of the electrodes. Therefore, entropy values were normalized by subtracting each electrode size group’s mean value.
To evaluate the similarity between simulation and real data arrays, we computed the average correlation across their electrode histograms (normalized by the mean correlation within electrodes of the real array). If our model is able to capture the effect of noise in real data, real electrodes with high noise levels should correlate better with high noise level simulations. To test this hypothesis, we used a linear model to evaluate the relationship between the estimated noise level in real electrodes and the noise level of the highest correlating (i.e., best fitting) simulation.
Grids with less than twelve electrodes were discarded. These grids have a high proportion of edge-located electrodes which can lead to less accurate estimations of the orthogonal vector and affect the following analysis.
The spatial organization of electrodes in strips precludes the computation of orthogonal vectors and, therefore strips were not included in the current analysis. However, strip electrodes are usually the same size as grid electrodes, and we assume that the conclusions drawn from the analysis of grids also apply to strips.
Example I: Evaluation of cluster-based localization methods
We demonstrate the applicability of our modeling platform by evaluating the performance of two localization approaches implemented in v1.000 and v1.010 of the iElectrodes toolbox. These are based on the k-means and k-medoids clustering algorithms, respectively (Hartigan and Wong, 1979; Kaufman and Rousseeuw, 1990). Before running the localization scripts, artifacts were thresholded between 0.25 and 0.6 depending on the noise level.
Localized coordinates were contrasted with simulated coordinates to evaluate their accuracy under several noise levels. The methods were compared with the Wilcoxon signed-rank tests, and effect sizes were approximated from their estimated z-scores (Pallant, 2007).
The function localizeArray.m and dependencies provide an exemplary code for testing localization algorithms. The script can be easily modified for comparison of other localization algorithms.
Example II: Interactive modeling of implanted electrodes in patient’s native brain space
To demonstrate the algorithm’s usefulness in more realistic scenarios, we applied our methods in a patient’s native brain space. We processed pre-implantation MRI and post-implantation CT images from an adult patient with drug-resistant epilepsy, a potential candidate for respective surgery, as described in Section “Patients.” The SCE surfaces for each hemisphere were computed following the steps described in Section “Smooth cortical envelope surface and seed points.”
We simulated grids, strips, and depth electrodes over multiple center points (seeds) and orientations using the interactive Graphical User Interface (GUI) provided by the iElectrodes toolbox (Blenkmann et al., 2017). For grids and strips, the 2D flat model arrays were first manually translated and rotated over the 3D SCE surface (Figure 2D). Once the desired location was reached, they were fitted to the SCE surface. Depth electrodes were defined by visually setting the target and entry points in the 2D views.
Results
In the following subsections, we will present the simulation results, their validation with real data, and examples.
Simulation and validation of electrode array coordinates
A total of 3,646 scenarios for grids and strips were simulated, and in 3,321 instances the arrays were successfully placed over the SCE surface. In 9% of the cases, the IED deformations between at least one pair of contacts exceeded the 5% threshold and were discarded.
For depth electrodes, the procedure resulted in 858 simulation scenarios for depth electrodes within the MNI brain.
We first validated our results by visual observation. Grid and strip electrode coordinates closely followed the SCE surface and showed an almost uniform spatial distribution. Changes in the IED were difficult to detect. For example, Figures 4A,B,D,E show simulated and real localized grid coordinates. Another important aspect was the simulation of overlaps, comparable to those observed in real situations (Figures 3E,F). Besides, the curved trajectories of simulated depth arrays mimicked those of real cases (Figures 4C,F).
Figure 4. Examples of real and simulated electrode coordinates. Observe that in real and simulated cases, the grid coordinates follow the brain curvature and that depth electrode trajectories are slightly bent. The top row shows real electrode coordinates from an 8 × 8, 3 mm IED array over the lateral temporal cortex (A), a 4 × 4, 10 mm IED array over the frontal pole (B), and an 8 contacts, 5 mm IED, depth electrode array into the temporal lobe (C). The bottom row shows simulated electrode arrays from an 8 × 8, 10 mm IED grid over the lateral fronto-temporo-parietal cortices (D), an 8 × 8, 3 mm IED high-density array over the frontal pole (E), and a 10 contacts depth electrode array into the temporal lobe (F). SCE and pial surfaces are depicted semi-transparently in the top plots to enhance visualization.
The implantation of subdural electrode arrays often introduces a brain shift, and this causes post-implantation CT images to capture the electrode artifacts “buried” under the pial or SCE surface reconstructed from pre-implantation MRIs. This deformation of the brain surface is not modeled in our current implementation.
To evaluate the quality of the simulated grid and strip coordinates, we quantified the deformations introduced by projecting them onto the SCE. We measured the distance between 1st, diagonal, and 2nd neighboring contacts, normalized in each case by the IED. The median deformations were 0.174, 0.169, and 0.243% of the IED, respectively. The median deformation increased for higher values of IED and SCE curvature (Figures 5A,B). As an upper bound reference, we measured the IED deformations from CT localized real grids (13 arrays, 484 electrodes). These deformations were one order of magnitude larger than the simulated ones (Supplementary Figure 3A).
Figure 5. Deformation between contacts for simulated grids and strips. Bar plots show that the median deformation between contacts increases with the IED (A), and the SCE curvature (B). As a reference, the observed deformations in real cases were an order of magnitude larger (Supplementary Figure 3A). Error bars denote 95% CI of the median. SCE, smooth cortical envelope; IED, inter-electrode distance.
To better quantify the effect of curvature and the number of electrodes on the simulations, we measured the distance between the fitted electrodes and the SCE. Although these values were constrained in the fitting procedure, some variability was observed. The electrodes’ median distance to the SCE decreased with the number of contacts and increased with the SCE curvature (Supplementary Figure 4). Overall, the median distance of electrodes to the SCE was 0.012 mm.
To validate the bending model of depth arrays, we measured deformations in real cases. Across real arrays, the maximum deformation (i.e., the maximum normalized distance to the main axis) was 0.82%, 95% CI [0.72, 1.01], comparable to the 1% maximum deformation implemented in our approach. Supplementary Figure 3B shows the bending profile over distance and the Lanczos function used in the simulations. Interestingly, the profile indicates that the maximum deformation is not centered, but closer to the deepest contact (normalized distance to target 0.37, 95% CI [0.33, 0.43]).
Simulation and validation of electrode computer tomography artifacts
We simulated the CT artifacts at each electrode location and manipulated the noise levels affecting these. By using the normal orientation vector, the procedure generated grid and strip CT artifacts aligned with the brain’s surface curvature. Similarly, depth electrode artifacts were aligned with the arrays’ principal axes. Altogether, we produced ∼40,000 simulations of grids and strips and ∼10,000 depth electrode CT artifacts.
To validate our procedure, we first visually inspected the CT artifacts in a subset of 100 randomly selected simulations. Grid and depth arrays have flat disc and cylinder-shaped metallic contacts, respectively, whereas their CT artifacts usually have ellipsoidal shapes following these geometries (Figures 6, 7). Importantly, similar shapes were achieved by the simulations (Figures 6, 7).
Figure 6. Examples of real and simulated electrode CT artifacts. Images on the top row (A,D,G,J) show electrode artifacts on CT images superimposed on co-registered MRIs from real cases. Plots on the middle row (B,E,H,K) show thresholded CT voxels corresponding to electrodes in relation to the SCE surface for the same cases as above. Plots on the bottom row (C,F,I,L) show simulations similar to the real examples above. Observe the similarity between the real and simulated CT artifacts. The examples selected have different noise levels to reflect some of the variability observed. The real grid case in (A,B), and the correspondingly simulated grid (C) have relatively low noise, whereas the real case depicted in (D,E), and simulation (F) show a higher noise level. The case depicted in (E) also shows a cable passing over the grid. Both real and simulated depth electrodes depicted in (G–I) have low noise levels, whereas cases shown in (J–L) have medium noise levels. Array details: (A,B) 4 × 8 grid, 10 mm IED, implanted over the temporal cortex. (D,E) 10 × 25 grid, 3 mm IED, implanted over the temporal cortex. (G,H) 12 electrodes depth array, 3.5 mm IED, implanted into the superior frontal cortex. (J,K) 8 electrodes depth array, 5 mm IED, implanted into the temporal cortex. (C) 8 × 8 grid, 10 mm IED, over the fronto-temporo-parietal cortex. (F) 4 × 8 grid, 3 mm IED, over the temporo-parietal cortex. (I) 18 electrodes depth array, 3 mm IED, implanted into the occipital cortex. (L) 10 electrodes depth array, 5 mm IED, implanted into the frontal cortex. For illustration purposes, the SCE surfaces are semi-transparent. IED, inter-electrode distance; SCE, smooth cortical envelope. A, anterior; P, posterior; R, right; S, superior.
Figure 7. Examples of real and simulated electrode CT artifacts and their 2D intensity-radius histograms. The three top rows show examples of real CT artifacts obtained from a 10 mm IED grid (A), a 3 mm IED high-density grid (B), and a 3 mm IED depth electrode array (C). The three bottom rows (D–F) show simulated cases with noise levels that best fit the real ones above (present the highest correlation). The left column plots show the spatial distribution of the voxels’ artifacts. The middle column shows the same artifacts after all individual electrodes were recentered at the origin (0,0,0), and their principal axes (orthogonally to the cortical surface for grids) or main axes (connecting outer electrodes for depths) aligned to the z-axis. For visualization purposes, only half of the electrode voxels are shown. The right column plots show their respective intensity vs. radius histogram. For illustrative purposes (C,F) middle and right column plots were obtained from 10 depth electrodes in the corresponding patient, or 10 simulated depth arrays, and not only from the individual arrays shown on the left column.
Even under the effect of brain shift, we observed grids and strips electrodes in the CT images to be approximately parallel to the SCE obtained from pre-implantation MRIs (Figures 6B,E). Accordingly, simulated CT artifacts followed the SCE curvature (Figures 6C,F).
Moreover, noise levels gradually affected the shape of real artifacts. In general, larger noise levels hampered the visual identification of individual artifacts (e.g., Figures 6E,F). Importantly, similar effects were obtained in the simulations. Other noise sources, like cables or metallic clips, could also affect the visual identification, but were not simulated.
To further characterize electrode artifacts and their noise, we computed 2D histograms to capture the bivariate distribution of voxels in terms of intensity and radial distance to the center. In the real cases (15 grids across 11 patients, IED 3 and 10 mm; and 77 depth electrodes from 14 patients, IED 3 and 5 mm), we observed a tendency for intensity to decrease with radius (e.g., Figures 7A–C, right column), a pattern that was mimicked by the simulations (e.g., Figures 7D–F, right column). These 2D histograms allowed us to estimate the noise level that affected the intensity distribution of voxels in space. For this purpose, we computed the mean entropy of the 2D cross-correlation of histograms as a proxy for the noise level in each array. Supplementary Figure 5 shows examples of low and high noise arrays and their corresponding average 2D cross-correlation images.
As expected, incrementing the noise level in the simulations produced a systematic increase in their average entropy (Supplementary Figures 6C,D, 7C,D). We also observed that the entropy measure not only captured the noise level, but was also affected by the size of the CT artifacts, being higher for bigger artifacts in real and simulated cases. Therefore, entropy values were normalized by subtracting each electrode size group’s mean value to account for the artifact size’s effect on the entropy measure.
Entropy was also computed for real CT artifacts, allowing us to compare real and simulated CT artifacts in terms of noise content (Supplementary Figures 6A,B, 7A,B).
To validate the modeled noise against the noise present on real CT artifacts, we correlated the 2D histograms of simulated artifacts with those of real artifacts (Supplementary Figure 8). Accordingly, we obtained a “best fit” simulation for each real case. As expected, low-noise real CT artifacts correlated better with low-noise simulations and vice versa.
To obtain an overall perspective of this relationship, we fitted a linear regression model between the estimated noise level (i.e., entropy) of the best fit simulations and the estimated noise level of the real data artifacts. Consequently, the model real_entropy_normalized = β0 + β1 simulation_entropy_normalized was used to assess their relationship. We observed a significant linear correlation between the estimated noise levels of the real CT artifacts and the estimated noise level of their best-fitting simulation counterparts (Figure 8; β0 = −0.010, β1 = 0.493, R2 = 0.413, F(1,92) = 64.409, p = 3.187e-12).
Figure 8. Correlation between real and simulated noise. The figure illustrates that real data noise linearly correlates with simulation noise. The entropy of the cross-correlation between 2D histograms was used as a proxy for noise levels. Individual marks represent a “best fit” simulation – real data pair. The best fit simulation was determined by the highest correlation between their corresponding 2D intensity-radius histograms (see Supplementary Figure 8). Entropy values were normalized (mean subtracted) to account for the effect of artifact size.
Example I: Evaluation of cluster-based localization methods
To show the usefulness of our simulation framework, we evaluated two electrode localization algorithms in ∼3,300 grid and strip scenarios under 12 different noise level realizations each (∼40,000 simulations in total). Coordinates were localized using k-means and k-medoids methods and contrasted with the ground truth simulated coordinates to evaluate their accuracy (Supplementary Figure 9 shows two exemplary results).
As expected, the median and maximum localization errors were affected by noise levels (Figure 9). K-medoids median localization error was below 15% of the IED when considering all noise levels and below 5% for medium and low noise levels. Instead, K-means showed larger median and maximum errors in all noise levels for IEDs of 3 and 5 mm arrays, and for the three highest noise levels in 10 mm IED arrays (Wilcoxon signed-rank test, p < 0.05, FDR corrected).
Figure 9. Localization error is affected by the noise level. Localization errors from k-means and k-medoids clustering algorithms. Localized coordinates were contrasted with the ground truth simulated coordinates to obtain localization errors. The median (A) and the maximum (B) localization error increase with the simulated noise level. Grids and strips of 3 and 5 mm IEDs are more likely to be affected by noise than arrays of 10 mm IED. The top images show the effect size of the difference between k-means and k-medoids methods. Non-significant differences are depicted with semitransparent colors.
Large effects (r > 0.5, Cohen, 1992) in the median localization error were observed across all noise levels in 5 mm IED grids. Instead, 3 and 10 mm IED arrays showed large effects for the largest three and the largest noise level, respectively (Figure 9A, top). Similar effects were observed for the maximum localization error (Figure 9B, top).
Example II: Interactive simulations in patient’s native brain space
To further demonstrate the usefulness of our platform, we tested the algorithms on a patient’s native anatomical brain space (i.e., unnormalized space). We successfully simulated grids, strips, and depth electrodes using the dedicated controls in the GUI. The interactive interface allowed us to define their locations and orientations precisely. The 2D and 3D views provided a clear interpretation of the array locations in relation to the brain’s pial surface and subcortical structures. Moreover, observing the subject-specific parcellation atlas aided in defining their anatomical location.
Figure 10 shows multiple electrode arrays simulated over the pial surface of a single patient. A high-density grid (5 mm IED) was simulated over the right fronto-temporo-parietal region, and its fifth row was particularly aligned with the superior temporal gyrus. Smaller grids were simulated over the fronto-parietal cortex, and strips were simulated to cover the anterior lateral frontal and orbitofrontal cortices, and the middle and inferior temporal gyrus. Note that grid electrodes were not forced to contact the pial- but the SCE surface.
Figure 10. Interactive electrode simulations on a patient’s native brain anatomy. Example of electrodes simulated on a native space brain anatomy using the interactive GUI provided in iElectrodes. The controls (magenta box) are used for displacing and rotating 2D flat grid or strip models over the SCE surface (similar to Figure 2D). The 3D view is particularly useful for defining the location of the arrays in relation to the lateral pial surfaces, whereas the 2D views help to identify their relation to sulci and deep structures. When the desired location and orientation are achieved, arrays are fitted to the SCE surface. For depth electrodes, target and entry points are manually defined in the 2D views. The 2D slices show the Destrieux atlas (color coded), the T1-weighted MRI image, and the simulated coordinates (green ‘x’), making the anatomical identification of the arrays easier. Note that electrodes were fitted to the SCE surface, but the pial brain surface is shown for illustrative purposes. SCE, smooth cortical envelope. GUI, graphical user interface.
Discussion
Summary of results
Several automatic or semi-automatic electrode localization algorithms have been developed in recent years (Arnulfo et al., 2015; Hamilton et al., 2017; Branco et al., 2018a,b; Granados et al., 2018; Centracchio et al., 2021). However, given the lack of ground truth data, the validation of these methods lacks standardization, imposing a major limitation on one of the most important benefits of analyzing brain activity from intracranial recordings, namely its excellent spatial resolution. Therefore, how can we compare the precision of different localization methods validated using different, usually small, datasets? Moreover, how do we know that these datasets represent important features (e.g., noise) of the data we are analyzing in a particular study? These questions are not easily answered, but we expect that simulations can provide a closer approximation to the ground truth in the attempt to answer them.
We provide the first modeling platform of intracranial electrodes to establish a test bench for electrode localization algorithms. Our simulations covered a wide range of realistic scenarios that can be useful for testing localization algorithms. To achieve our aim, we developed novel methods to model intracranial EEG electrode coordinates and the CT artifacts typically produced by these. We mimicked realistic scenarios by simulating implanted arrays of various geometrical dimensions in a wide range of cortical curvatures. We simulated grid and strip coordinates by fitting models to the SCE surface. In addition, we modeled depth electrodes between anatomical target and entry points and precisely defined their bending profile.
Then, we simulated the CT artifacts at the electrodes’ coordinates. The signal intensity distribution over space, the shape of the artifacts, and the orientation of each cluster of high-intensity voxels were carefully modeled. Moreover, different noise levels were simulated and validated with real data.
An extensive database with ∼50,000 cases is freely available for this purpose (see Data availability statement). We focused on providing control over those situations where localization algorithms might fail or encounter difficulties, such as high-density arrays (grids and depth electrodes), high noise levels, overlapping grids, or highly curved implants.
Finally, we showed the impact of noise on the accuracy of two electrode localization algorithms.
Modeling implantation coordinates
From a visual perspective, simulated coordinates in grids, strips, and depth electrode arrays resemble those of real data. Grid and strips naturally followed the SCE curvature with almost uniformly distributed coordinates, whereas depth electrodes depicted realistically curved trajectories. Moreover, we modeled overlapping CT artifacts for grids and strips, controlling the size, orientation, and the number of overlapping electrodes. These features add realism to the models not shown before.
To objectively evaluate the accuracy of grid and strip implantation models, we measured their deformation (i.e., changes in the inter-electrode distance). Fitting grids or strips over the SCE surface introduce deformations since plane- or line-shaped models must bend to fit the curved brain envelope surface. Importantly, deformation values were relatively small (<0.3% of the IED), indicating that the resulting models are reliable. As an upper bound reference for our method, we measured the IED deformations in real arrays. These deformations were one order of magnitude larger, indicating that our method did not introduce unrealistically high deformations. In our approach, we discarded simulations where the individual IED deformations exceeded 5% of the IED, which might be overly conservative.
As shown in Figure 5A, grid deformations increased with larger IEDs. We theorized that bigger grids had to deform more to cover a larger extent of the curved SCE surface, in contrast to smaller grids that suffer less deformation in the reduced area they overlay. Similarly, higher levels of curvature produced larger deformation values (Figure 5B).
Parameter K (Equation 1) controls the deformation introduced in the fitting procedure. We defined K = 1,000, which allowed very small deformations (Trotta et al., 2018). However, for large high-density grids (IED of 3 mm), we had to reduce the value to K = 100 since the algorithm did not reach feasible solutions in most scenarios. We suggest using the highest K value that reaches a feasible solution, therefore introducing the smallest deformations.
Grid and strip electrodes fitting to the SCE were constrained by ε, a tolerance distance between these two. Overall, the resulting distances between electrodes and SCE were negligible compared to the IEDs, indicating successful algorithm performance.
Depth electrode arrays, on the other hand, bend during implantation, mainly due to brain shift.
Our simulations successfully resulted in deformation values equivalent to those observed in real data (Supplementary Figure 3A) and those reported in clinical practice (∼1–3 mm at the tip; Cardinale et al., 2013; Vakharia et al., 2017; Granados et al., 2018). Interestingly, the bending pattern was not symmetrical as expected, and further research is needed to understand its behavior.
Our approach provides precise bending control by applying an arc-shaped function that departs from the original straight linear model. We used a symmetric Lanczos window function for this purpose. However, the use of other arc-shaped functions is straightforward.
Our approach does not consider deviations from the planned implantation trajectory, which can introduce other errors (Vakharia et al., 2017; Granados et al., 2018).
Besides the simulations on MNI space, we also showed the feasibility of simulating on individual (real) subject space. This adds to building even more realistic simulations that could help map expected localization errors at the individual level.
Simulation of computer tomography artifacts, noise, and overlapping electrodes
CT artifacts corresponding to grids and strip electrodes were placed over the SCE surface, keeping the electrode’s principal axes orthogonal to the surface. Meanwhile, the ones corresponding to depth electrode arrays penetrating the brain were aligned to the arrays’ main axes (Figure 6). From a visual perspective, simulated artifacts echo those observed in real data.
Previously, models of individual electrode CT artifacts were done as uniform intensity cylinders without considering the details examined in the current study (Brang et al., 2016). We defined ellipsoid-shaped artifacts, with the intensity changing as a function of the radius. The electrode models produced a histogram profile resembling real artifacts, i.e., an intensity decrease with increasing radius (Figure 7).
The introduction of noise played a substantial role in producing realistic CT artifacts and was captured by the intensity-radius histograms (Figures 6, 7). Adding correlated noise introduced realistic shape deformations and was visually perceived to affect the center of mass of the electrodes. This feature provides a wider range of difficult scenarios and makes localization more challenging.
A classical approach for characterizing noise in medical images is the subtraction of a smoothed mean image. However, this approach is inappropriate for intensity transitions like those present in intracranial electrode artifacts (Gravel et al., 2004). To circumvent this limitation, we used an information theory approach inspired by techniques quantifying image focus (Bove, 1993; Kristan and Pernu, 2004). In this respect, the measure’s sensitivity to simulated noise was an important consistency check (Supplementary Figures 6, 7), which allowed us to show that noise levels in simulated artifacts mimic those observed in real data (Figure 8). Note that our model explained ∼40% of the real noise variance, indicating their high similarity.
The relevance of intracranial electrode models for localization algorithms
Over the last decade, several approaches have been proposed to localize intracranial electrodes based on CT and MRI images. The majority used post-implantation CT and pre-implantation MRI images.
The detection of CT artifacts has typically been a manual process (Princich et al., 2013), but has recently been approached by semiautomatic techniques such as clustering voxels of high intensity (Taimouri et al., 2014; Blenkmann et al., 2015, 2017; Brang et al., 2016; Qin et al., 2017; Branco et al., 2018a; Granados et al., 2018), shape analysis (Centracchio et al., 2021), or the interpolation of coordinates given entry and target points in depth electrodes (Arnulfo et al., 2015; Li et al., 2020). Moreover, several approaches have been integrated in novel processing pipelines (LaPlante et al., 2016; Blenkmann et al., 2017; Groppe et al., 2017; Stolk et al., 2018; Li et al., 2020; Davis et al., 2021; Rockhill et al., 2022), providing users several alternatives and even handling group studies (Deman et al., 2018).
Noise and overlapping electrodes are two well-known difficulties for localization algorithms, which hamper the success of automatic methods. For example, Brang et al. (2016) excluded overlapping electrodes from the analysis, given the resulting difficulties, while others treated such cases manually (Taimouri et al., 2014; LaPlante et al., 2016; Branco et al., 2018a; Centracchio et al., 2021). In the same vein, Narizzano et al. (2017) observed errors in their depth electrode estimations associated with other electrodes in close proximity, requiring manual intervention from the user. Moreover, noise signals could be mistakenly detected as electrodes (LaPlante et al., 2016), whereas CT image resolution affects localization accuracy (Brang et al., 2016). Apart from the studies above, the effect of SNR on electrode localization algorithms’ accuracy was rarely discussed, most likely due to the lack of standardized measures to quantify the noise level. The framework proposed here provides a controlled simulation of noise levels and overlapping electrodes, allowing performance evaluation of localization algorithms.
Recently, Centracchio et al. (2021) proposed a novel approach to detect electrode artifacts using a Gaussian support vector machine. Classification accuracy was very high in the analyzed datasets. However, these were restricted to a small number of geometries and electrode sizes. Similar results were obtained by applying deep learning to depth SEEG and DBS electrodes (Vlasov et al., 2021). Approaches like these could strongly benefit from large databases, such as the ones offered here, making their results generalizable across a wider range of geometries.
The spatial resolution of grids and depth electrodes has increased over the last years (Chang, 2015; Erhardt et al., 2020), and high-density arrays are more informative than low-density ones in both cognitive (Gupta et al., 2014; Jiang et al., 2018) and clinical research (Stead et al., 2010). High-density electrodes require additional spatial precision and can be an obstacle for many of the frequently used localization algorithms (Hamilton et al., 2017; Narizzano et al., 2017; Branco et al., 2018b). Simulations can be a reliable platform for developing novel localization techniques for high-density electrode arrays.
Importantly, our models introduced deformations on the order of 1/1,000 of the inter-electrode distance, and distances between electrodes and the SCE on the order of 1/100 mm, which ensure precise modeling of the grids and strips. These errors and deformations were negligible compared to those observed with previous localization (∼0.2–0.6 mm; Blenkmann et al., 2017; Narizzano et al., 2017) and brain-shift correction algorithms (∼2–3 mm; Brang et al., 2016; Branco et al., 2018a,b; Trotta et al., 2018).
As a proof of concept, we evaluated two electrode localization algorithms, k-means and k-medoids. Knowing the ground truth location allowed us to characterize how noise levels affected their accuracy, and to have strong evidence to conclude that k-medoids perform better than k-means in most situations. In a similar way, the current framework has been used to develop a novel localization algorithm called GridFit, which is embedded in the iElectrodes toolbox1 (Blenkmann et al., 2017; more details about the algorithm will be provided in a forthcoming publication). The simulation platform was used to define the optimal set of parameters needed to precisely localize electrodes under controlled levels of noise, overlap, and curvature. The simulated scenarios helped to identify the sensitivity of the parameters to these variables independently. The same achievement would have been impossible using real data, given the scarce nature of the data and the lack of control of the variables of interest. It is our intent that others find the simulation platform useful to develop novel localization tools.
Assumptions, limitations, and future directions
Although we provide a substantial number of scenarios and a good starting point to model implanted electrodes, there are some limitations in the current models.
First, we took a simple approach to the spatial intensity distribution of CT artifacts. The simplistic assumption allowed us to build realistic models of large arrays of electrodes.
However, more sophisticated approaches could be implemented, considering the x-ray interaction with metallic electrodes (i.e., beam hardening, scatter effects, and Poisson noise), and the corresponding alteration produced in the image reconstructions (Boas and Fleischmann, 2011; Katsura et al., 2018). Developments in this direction could pave the way to model the artifacts produced by microwires at the tip of depth electrodes (e.g., Behnke–Fried electrodes, Ad-Tech Medical) and the design of novel localization algorithms for this specific and unsolved problem.
Moreover, alterations in the electrode shapes might be caused by beam hardening (Boas and Fleischmann, 2011; Centracchio et al., 2021), and in some situations, streak artifacts can be prominent in CT images, but these are typically removed by thresholding (LaViolette et al., 2011b). Accordingly, we modeled electrode artifacts after thresholding, as most algorithms are applied to this type of data. However, there might be some cases where thresholding is a fundamental step to assess (e.g., Davis et al., 2021). CT image augmentation and more complex models need to be developed for such instances.
Second, the implantation of intracranial grids and strips is a procedure that results in brain tissue deformation. Deformations of 10 mm or more can occur on the brain surface around the electrodes or in deeper brain structures due to cerebrospinal fluid loss in the ventricles (Studholme et al., 2001; LaViolette et al., 2011a). Implantation of depth electrodes might also produce brain deformations, but to a much lower extent, with a smaller amount of cerebrospinal fluid loss, if any (Elias et al., 2007). Modeling brain deformations is a challenging problem, where multiple variables have to be considered, such as the size and location of the skull opening, the amount of cerebrospinal fluid loss, and the swelling of soft tissue, among other factors (Studholme et al., 2001). Given the complexity of the problem, we assumed non-deformed brains in our simulations, which precludes their use to evaluate brain-shift correction algorithms. The use of non-linear finite element methods can be a successful way to model these more complex brain deformations (Wittek et al., 2007).
Finally, we should mention that the volume covered by the target coordinates in the depth electrode simulations is far from covering the full brain volume. Similarly, we did not cover the complete lateral cortex with grids and strips. However, this is possible in the current framework by carefully selecting the entry and target coordinates, seed points, electrode sizes, and geometries.
Conclusion
Intracranial EEG recordings allow us to study brain function with excellent spatial resolution and rely on precisely localizing the implanted electrodes. Here, we presented the first platform to model electrode coordinates and CT artifacts of implanted grids, strips, and depth electrodes.
Implanted electrodes under realistic scenarios were successfully modeled with high accuracy, resembling real cases. These methods enable the systematic and quantitative evaluation of electrode localization strategies, contributing to the development of more accurate techniques. The platform could be a starting point for more sophisticated models, such as brain tissue deformations or microwires.
The modeling methods and simulation results are freely available to the research community via open repositories. Moreover, a graphical user interface implementation is also available via the open-source iElectrodes toolbox.
Data availability statement
The scripts and simulated datasets presented in this study can be found in online repositories. The name of the repository and accession number can be found below: Center for Open Science (COS) Open Science Framework (OSF), https://osf.io/9fsm3/, doi: 10.17605/OSF.IO/9FSM3. A graphical interface for the interactive simulation of electrode coordinates is available within the iElectrodes toolbox at https://sourceforge.net/projects/ielectrodes/. The patients’ datasets analyzed in this study are not publicly available due to our ethical approval conditions that do not permit public archiving of anonymized study data.
Ethics statement
The studies involving human participants were reviewed and approved by the Regional Committees for Medical and Health Research Ethics, Region North Norway (REK 2015/175), and the Human Subjects Committees at UCSF, UC Irvine, and UC Berkeley. The patients provided their written informed consent to participate in this study.
Author contributions
AOB, A-KS, and TE designed this study. AOB coded and tested the software tools, performed the analyses, and wrote the manuscript. TE, A-KS, PL, JI, and RK provided the data and validated the results. All authors revised the manuscript, read, and approved the final manuscript.
Funding
This work was partly supported by the Research Council of Norway, project numbers 240389 and 314925, through its Centres of Excellence scheme, project number 262762, and Intpart, project number 274996, NINDS Grant R37NS21135, NIMH CONTE Center P50MH109429, and Brain Initiative U01-NS108916.
Acknowledgments
We are grateful to the patients for kindly participating in our study. We thank FRONT Neurolab/RITMO members for their rich discussions and the reviewers for their useful comments and suggestions that improved the quality of this manuscript. This study has previously appeared online at https://www.biorxiv.org/content/10.1101/2021.04.28.441828v1.
Conflict of interest
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.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fninf.2022.788685/full#supplementary-material
Abbreviations
IED, inter-electrode distance; SCE, smooth cortical envelope; SNR, signal to noise ratio.
Footnotes
References
Arnulfo, G., Narizzano, M., Cardinale, F., Fato, M. M., and Palva, J. M. (2015). Automatic segmentation of deep intracerebral electrodes in computed tomography scans. BMC Bioinformatics 16:99. doi: 10.1186/s12859-015-0511-6
Blenkmann, A., Phillips, H., Princich, J. P., and Kochen, S. (2015). “Grid and depth intracranial electrodes localization in a normalized space using MRI and CT images,” in VI Latin American Congress on Biomedical Engineering CLAIB 2014, Paraná, Argentina 29, 30 & 31 October 2014, A. Braidot and A. Hadad (eds) (Cham: Springer) Vol. 64, 413–416. doi: 10.1007/978-3-319-13117-7_106
Blenkmann, H. N., Princich, J. P., Rowe, J. B., Bekinschtein, T. A., et al. (2017). iElectrodes: A Comprehensive Open-Source Toolbox for Depth and Subdural Grid Electrode Localization. Front. Neuroinformatics 11:14. doi: 10.3389/fninf.2017.00014
Blenkmann, O., Collavini, S., Lubell, J., Llorens, A., Funderud, I., et al. (2019). Auditory deviance detection in the human insula: An intracranial EEG study. Cortex 121, 189–200. doi: 10.1016/j.cortex.2019.09.002
Boas, F. E., and Fleischmann, D. (2011). Evaluation of two iterative techniques for reducing metal artifacts in computed tomography. Radiology 259, 894–902.
Bove, V. M. (1993). Entropy-based depth from focus. J. Optical Soc. Am. 10:561. doi: 10.1364/josaa.10.000561
Branco, M. P., Leibbrand, M., Vansteensel, M. J., Freudenburg, Z. V., and Ramsey, N. F. (2018b). GridLoc: An automatic and unsupervised localization method for high-density ECoG grids. NeuroImage 179, 225–234. doi: 10.1016/j.neuroimage.2018.06.050
Branco, M. P., Gaglianese, A., Glen, D. R., Hermes, D., Saad, Z. S., Petridou, N., et al. (2018a). ALICE: A tool for automatic localization of intracranial electrodes for clinical and high-density grids. J. Neurosci. Methods 301, 43–51. doi: 10.1016/j.jneumeth.2017.10.022
Brang, D., Dai, Z., Zheng, W., and Towle, V. L. (2016). Registering imaged ECoG electrodes to human cortex: A geometry-based technique. J. Neurosci. Methods 273, 64–73. doi: 10.1016/j.jneumeth.2016.08.007
Britten, A. J., Crotty, M., Kiremidjian, H., Grundy, A., and Adam, E. J. (2004). The addition of computer simulated noise to investigate radiation dose and image quality in images with spatial correlation of statistical noise: An example application to X-ray CT of the brain. Br. J. Radiol. 77, 323–328. doi: 10.1259/bjr/78576048
Cardinale, F., Cossu, M., Castana, L., Casaceli, G., Schiariti, M. P., Miserocchi, A., et al. (2013). Stereoelectroencephalography: Surgical methodology, safety, and stereotactic application accuracy in 500 procedures. Neurosurgery 72, 353–366. doi: 10.1227/NEU.0b013e31827d1161
Centracchio, J., Sarno, A., Esposito, D., Andreozzi, E., Pavone, L., Di Gennaro, G., et al. (2021). Efficient automated localization of ECoG electrodes in CT images via shape analysis. Int. J. Comput. Assist. Radiol. Surgery 16: 543–554. doi: 10.1007/s11548-021-02325-0
Chang, E. F. (2015). Towards Large-Scale, Human-Based, Mesoscopic Neurotechnologies. Neuron 86, 68–78. doi: 10.1016/j.neuron.2015.03.037
Dale, A. M., Fischl, B., and Sereno, M. I. (1999). Cortical surface-based analysis. I. Segmentation and surface reconstruction. NeuroImage 9, 179–194. doi: 10.1006/nimg.1998.0395
Davis, T. S., Caston, R. M., Philip, B., Charlebois, C. M., Anderson, D. N., Weaver, K. E., et al. (2021). LeGUI: A Fast and Accurate Graphical User Interface for Automated Detection and Anatomical Localization of Intracranial Electrodes. Front. Neurosci. 15:769872. doi: 10.3389/fnins.2021.769872
Deman, P., Bhattacharjee, M., Tadel, F., Job, A. S., Rivière, D., Cointepas, Y., et al. (2018). IntrAnat Electrodes: A Free Database and Visualization Software for Intracranial Electroencephalographic Data Processed for Case and Group Studies. Front. Neuroinformatics 12:40. doi: 10.3389/fninf.2018.00040
Destrieux, C., Fischl, B., Dale, A., and Halgren, E. (2010). Automatic parcellation of human cortical gyri and sulci using standard anatomical nomenclature. NeuroImage 53, 1–15. doi: 10.1016/j.neuroimage.2010.06.010
Dykstra, A. R., Chan, A. M., Quinn, B. T., Zepeda, R., Keller, C. J., Cormier, J., et al. (2012). Individualized localization and cortical surface-based registration of intracranial electrodes. NeuroImage 59, 3563–3570.
Elias, W. J., Fu, K.-M., and Frysinger, R. C. (2007). Cortical and subcortical brain shift during stereotactic procedures. J. Neurosurg. 107, 983–988.
Erhardt, J. B., Lottner, T., Pasluosta, C. F., Gessner, I., Mathur, S., Schuettler, M., et al. (2020). Fabrication and validation of reference structures for the localization of subdural standard- and micro-electrodes in MRI. J. Neural Eng. 17:046044.
Fang, Q., and Boas, D. A. (2009). “Tetrahedral mesh generation from volumetric binary and grayscale images,” in 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro, (Piscataway, NJ: IEEE), 1142–1145.
Frauscher, B., von Ellenrieder, N., Zelmann, R., Doležalová, I., Minotti, L., Olivier, A., et al. (2018). Atlas of the normal intracranial electroencephalogram: neurophysiological awake activity in different cortical areas. Brain 141, 1130–1144. doi: 10.1093/brain/awy035
Granados, A., Vakharia, V., Rodionov, R., Schweiger, M., Vos, S. B., O’Keeffe, A. G., et al. (2018). Automatic segmentation of stereoelectroencephalography (SEEG) electrodes post-implantation considering bending. Int. J. Comput. Assist. Radiol. Surg. 13, 935–946. doi: 10.1007/s11548-018-1740-8
Gravel, P., Beaudoin, G., and DeGuise, J. A. (2004). A Method for Modeling Noise in Medical Images. IEEE Trans. Med. Imaging 23, 1221–1232. doi: 10.1109/TMI.2004.832656
Groppe, D. M., Bickel, S., Dykstra, A. R., Wang, X., Mégevand, P., Mercier, M. R., et al. (2017). iELVis: An open source MATLAB toolbox for localizing and visualizing human intracranial electrode data. J. Neurosci. Methods 281, 40–48. doi: 10.1016/j.jneumeth.2017.01.022
Gupta, D., Hill, N. J., Adamo, M. A., Ritaccio, A., and Schalk, G. (2014). Localizing ECoG electrodes on the cortical anatomy without post-implantation imaging. NeuroImage 6, 64–76. doi: 10.1016/j.nicl.2014.07.015
Hamilton, L. S., Chang, D. L., Lee, M. B., and Chang, E. F. (2017). Semi-automated Anatomical Labeling and Inter-subject Warping of High-Density Intracranial Recording Electrodes in Electrocorticography. Front Neuroinform. 11:62. doi: 10.3389/fninf.2017.00062
Hartigan, J. A., and Wong, M. A. (1979). Algorithm AS 136: A K-Means Clustering Algorithm. Appl. Statist. 28:100. doi: 10.2307/2346830
Hinds, W. A., Misra, A., Sperling, M. R., Sharan, A., Tracy, J. I., and Moxon, K. A. (2018). Enhanced co-registration methods to improve intracranial electrode contact localization. NeuroImage 20, 398–406. doi: 10.1016/j.nicl.2018.07.026
Jiang, T., Liu, S., Pellizzer, G., Aydoseli, A., Karamursel, S., Sabanci, P. A., et al. (2018). Characterization of Hand Clenching in Human Sensorimotor Cortex Using High-, and Ultra-High Frequency Band Modulations of Electrocorticogram. Front. Neurosci. 12:110. doi: 10.3389/fnins.2018.00110
Katsura, M., Sato, J., Akahane, M., Kunimatsu, A., and Abe, O. (2018). Current and novel techniques for metal artifact reduction at CT: Practical guide for radiologists. Radiographics 38, 450–461. doi: 10.1148/rg.2018170102
Kaufman, L., and Rousseeuw, P. J. (1990). Partitioning Around Medoids (Program PAM). In Finding Groups in Data: An Introduction to Cluster Analysis. Hoboken, NJ: John Wiley & Sons, Ltd, doi: 10.1002/9780470316801.ch2
Kristan, M., and Pernu, F. (2004). “Entropy based measure of camera focus,” in Proceedings of the thirteenth electrotechnical and computer science conference ERKI (Portorož: IEEE), 179–182.
Lachaux, J. P., Rudrauf, D., and Kahane, P. (2003). Intracranial EEG and human brain mapping. J. Physiol. Paris 97, 613–628. doi: 10.1016/j.jphysparis.2004.01.018
LaPlante, R. A., Tang, W., Peled, N., Vallejo, D. I., Borzello, M., Dougherty, D. D., et al. (2016). The interactive electrode localization utility: software for automatic sorting and labeling of intracranial subdural electrodes. Int. J. Comput. Assist. Radiol. Surg. 12:1829–1837. doi: 10.1007/s11548-016-1504-2
LaViolette, P. S., Rand, S. D., Raghavan, M., Ellingson, B. M., Schmainda, K. M., and Mueller, W. (2011b). Three-dimensional visualization of subdural electrodes for presurgical planning. Neurosurgery 68, 152–161.
LaViolette, P. S., Rand, S. D., Ellingson, B. M., Raghavan, M., Lew, S. M., Schmainda, K. M., et al. (2011a). 3D visualization of subdural electrode shift as measured at craniotomy reopening. Epilepsy Res. 94, 102–109.
Li, G., Jiang, S., Chen, C., Brunner, P., Wu, Z., Schalk, G., et al. (2019). iEEGview: An open-source multifunction GUI-based Matlab toolbox for localization and visualization of human intracranial electrodes. J. Neural Eng. 17:016016. doi: 10.1088/1741-2552/ab51a5
Li, G., Jiang, S., Chen, C., Brunner, P., Wu, Z., Schalk, G., et al. (2020). IEEGview: An open-source multifunction GUI-based Matlab toolbox for localization and visualization of human intracranial electrodes. J. Neural Eng. 17:016016.
Martin, S., Iturrate, I., Millán, J., del, R., Knight, R. T., and Pasley, B. N. (2018). Decoding inner speech using electrocorticography: Progress and challenges toward a speech prosthesis. Front. Neurosci. 12:422. doi: 10.3389/fnins.2018.00422
Narizzano, M., Arnulfo, G., Ricci, S., Toselli, B., Tisdall, M., Canessa, A., et al. (2017). SEEG assistant: a 3DSlicer extension to support epilepsy surgery. BMC Bioinformatics 18:124. doi: 10.1186/s12859-017-1545-8
Pallant, J. (2007). SPSS survival manual: a step by step guide to data analysis using SPSS for Windows (version 15). 335.
Parvizi, J., and Kastner, S. (2017). Human intracranial EEG: Promises and Limitations. Nat. Neurosci. 21, 474–483. doi: 10.1038/s41593-018-0108-2
Princich, J. P., Wassermann, D., Latini, F., Oddo, S., Blenkmann, A., Seifer, G., et al. (2013). Rapid and efficient localization of depth electrodes and cortical labeling using free and open source medical software in epilepsy surgery candidates. Front. Neurosci. 7:260. doi: 10.3389/fnins.2013.00260
Qin, C., Tan, Z., Pan, Y., Li, Y., Wang, L., Ren, L., et al. (2017). Automatic and Precise Localization and Cortical Labeling of Subdural and Depth Intracranial Electrodes. Front. Neuroinform. 11:10. doi: 10.3389/fninf.2017.00010
Rockhill, A., Larson, E., Stedelin, B., Mantovani, A., Raslan, A., Gramfort, A., et al. (2022). Intracranial Electrode Location and Analysis in MNE-Python. J. Open Sour. Softw. 7:3897. doi: 10.21105/joss.03897
Rosenow, F., and Lüders, H. (2001). Presurgical evaluation of epilepsy. Brain 124, 1683–1700. doi: 10.1093/brain/124.9.1683
Rusinkiewicz, S. (2004). “Estimating curvatures and their derivatives on triangle meshes,” in Proceedings - 2nd International Symposium on 3D Data Processing, Visualization, and Transmission. 3DPVT,(Piscataway, NJ) Vol. 2004, 486–493. doi: 10.1109/TDPVT.2004.1335277
Stead, M., Bower, M., Brinkmann, B. H., Lee, K., Marsh, W. R., Meyer, F. B., et al. (2010). Microseizures and the spatiotemporal scales of human partial epilepsy. Brain 133, 2789–2797. doi: 10.1093/brain/awq190
Stolk, A., Griffin, S., Van Der Meij, R., Dewar, C., Saez, I., Lin, J. J., et al. (2018). Integrated analysis of anatomical and electrophysiological human intracranial data. Nat. Protoc. 13, 1699–1723. doi: 10.1038/s41596-018-0009-6
Studholme, C., Hill, D. L. G., and Hawkes, D. J. (1999). An overlap invariant entropy measure of 3D medical image alignment. Pattern Recognit. 32, 71–86. doi: 10.1016/S0031-3203(98)00091-0
Studholme, C., Novotny, E., Zubal, I. G., and Duncan, J. S. (2001). Estimating tissue deformation between functional images induced by intracranial electrode implantation using anatomical MRI. NeuroImage 13, 561–576. doi: 10.1006/nimg.2000.0692
Taimouri, V., Akhondi-Asl, A., Tomas-Fernandez, X., Peters, J. M., Prabhu, S. P., Poduri, A., et al. (2014). Electrode localization for planning surgical resection of the epileptogenic zone in pediatric epilepsy. Int. J. Comput. Assist. Radiol. Surg. 9, 91–105. doi: 10.1007/s11548-013-0915-6
Trotta, M. S., Cocjin, J., Whitehead, E., Damera, S., Wittig, J. H., Saad, Z. S., et al. (2018). Surface based electrode localization and standardized regions of interest for intracranial EEG. Hum. Brain Mapp. 39, 709–721. doi: 10.1002/hbm.23876
Vakharia, V. N., Sparks, R., O’Keeffe, A. G., Rodionov, R., Miserocchi, A., McEvoy, A., et al. (2017). Accuracy of intracranial electrode placement for stereoencephalography: A systematic review and meta-analysis. Epilepsia 58, 921–932. doi: 10.1111/epi.13713
Vlasov, V., Bofferding, M., Marx, L., Zhang, C., Goncalves, J., Husch, A., et al. (2021). Automated Deep Learning-based Segmentation of Brain, SEEG and DBS Electrodes on CT Images. Informatik aktuell 92–97. doi: 10.1007/978-3-658-33198-6_22
Keywords: SEEG, ECoG, iEEG, intracranial electrodes, depth electrodes, subcortical grids, subdural grids
Citation: Blenkmann AO, Solbakk A-K, Ivanovic J, Larsson PG, Knight RT and Endestad T (2022) Modeling intracranial electrodes. A simulation platform for the evaluation of localization algorithms. Front. Neuroinform. 16:788685. doi: 10.3389/fninf.2022.788685
Received: 03 October 2021; Accepted: 05 September 2022;
Published: 06 October 2022.
Edited by:
Pedro Antonio Valdes-Sosa, University of Electronic Science and Technology of China, ChinaReviewed by:
Torbjørn Vefferstad Ness, Norwegian University of Life Sciences, NorwayAntonio Sarno, University of Naples Federico II, Italy
Thomas Lottner, University Medical Center Freiburg, Germany
Copyright © 2022 Blenkmann, Solbakk, Ivanovic, Larsson, Knight and Endestad. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Alejandro O. Blenkmann, YWJsZW5rbWFubkBnbWFpbC5jb20=