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.

The following steps are performed for modeling subdural grids and strips:

is minimized by varying the electrode coordinates to the final position x_i, where E_t is the translation energy, E_d is the deformation energy, and K is a constant value. The energy minimization is constrained to the electrode coordinates x_i being closely located over the SCE surface, i.e., ∀i, ||x_i−s_i|| ≤ ϵ, where s_i is the closest point in the SCE surface to electrode i, and ε is a tolerance distance.

E_t (Figure 2B) represents the energy needed to translate the electrode coordinates from the original position of each electrode xi0 to the final position x_i on the SCE surface:

where N is the number of electrodes.

E_d (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 d_ij is the distance between electrodes, dij0 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). E_d is typically interpreted as the deformation of “springs” connecting the electrodes (Dykstra et al., 2012).

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.

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:

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

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:

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.

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 N(0,σcorr2) distribution with a kernel provided by Britten et al. (2004). Uncorrelated noise values were drawn from a N(0,σuncorr2) distribution. Therefore, the total noise variance was σtotal2=σcorr2+σuncorr2. 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 σtotal2 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.

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.

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.

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

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

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]).

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

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 = β₀ + β₁ 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.010, β₁ = 0.493, R² = 0.413, F(1,92) = 64.409, p = 3.187e-12).

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

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

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.

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.

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.

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.

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 toolbox¹ (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.

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.