Anatomical registration of intracranial electrodes. Robust model-based localization and deformable smooth brain-shift compensation methods

Precise electrode localization is important for maximizing the utility of intracranial EEG data. Electrodes are typically localized from post-implantation CT artifacts, but algorithms can fail due to low signal-to-noise ratio, unrelated artifacts, or high-density electrode arrays. Minimizing these errors usually requires time-consuming visual localization and can still result in inaccurate localizations. In addition, surgical implantation of grids and strips typically introduces non-linear brain deformations, which result in anatomical registration errors when post-implantation CT images are fused with the pre-implantation MRI images. Several projection methods are currently available, but they either fail to produce smooth solutions or do not account for brain deformations. To address these shortcomings, we propose two novel algorithms for the anatomical registration of intracranial electrodes that are almost fully automatic and provide highly accurate results. We first present GridFit, an algorithm that simultaneously localizes all contacts in grids, strips, or depth arrays by fitting flexible models to the electrodes’ CT artifacts. We observed localization errors of less than one millimeter (below 8% relative to the inter-electrode distance) and robust performance under the presence of noise, unrelated artifacts, and high-density implants when we ran ~6000 simulated scenarios. Furthermore, we validated the method with real data from 20 intracranial patients. As a second registration step, we introduce CEPA, a brain-shift compensation algorithm that combines orthogonal-based projections, spring-mesh models, and spatial regularization constraints. When tested with real data from 15 patients, anatomical registration errors were smaller than those obtained for well-established alternatives. Additionally, CEPA accounted simultaneously for simple mechanical deformation principles, which is not possible with other available methods. Inter-electrode distances of projected coordinates smoothly changed across neighbor electrodes, while changes in inter-electrode distances linearly increased with projection distance. Moreover, in an additional validation procedure, we found that modeling resting-state high-frequency activity (75–145 Hz ) in five patients further supported our new algorithm. Together, GridFit and CEPA constitute a versatile set of tools for the registration of subdural grid, strip, and depth electrode coordinates that provide highly accurate results even in the most challenging implantation scenarios. The methods presented here are implemented in the iElectrodes open-source toolbox, making their use simple, accessible, and straightforward to integrate with other popular toolboxes used for analyzing electrophysiological data.


Definition of GridFit parameters using simulations
Three rounds of localizing simulated CT artifacts was used to determine the best parameters for the GridFit algorithm. The first round was intended to get an overview of the effects of the parameters, the second round to have a fine-grained description of the parameters, and the third one to assess the method's performance.
After the first round, we noticed that kCorr and kDef had opposing effects. For the second round, we defined all parameters for the first fitting step, i.e., kTrans= 1E-6, kDef = 1E5, kCorr = 1E5, σ = IED, ε = 1 [mm], and defined kTrans = 1E-2, kDef = 1E5, and σ= ¼ IED for the second fitting. We evaluated the algorithm's performance by varying the kCorr parameter (kCorr range 10 to 1 E11, localizing in total ~69000 simulated cases). We ran the second fitting step without constraints to incorporate a larger range of parameters with less computational time (Eq. 7, 8, and 9). We observed that optimal kCorr parameters, i.e., the ones producing a smaller normalized median error (Eq. 10) or higher accuracy, depended mostly on: i) electrode type, ii) the number of electrodes, iii) the IED, and iv) the presence of overlaps (for grids and strips). Smaller variations were observed to changes in noise levels or local curvature. Therefore, we defined a function that automatically computes the optimal kCorr parameter, i.e., kCorr Opt, considering the four variables. Supplementary Tables 2, 3, and 4 show the optimal values obtained for each case. Optimal KCorr value for other geometries than the evaluated ones are obtained by a modified Akima cubic Hermite interpolation (using the interp3 function in Matlab).
Supplementary  In the third round, we evaluated the accuracy of our method given optimal or near-optimal parameters. All constraints were used in the two energy-minimization steps, i.e., the procedure was fully-functional. Between 20 and 30 simulations per condition were evaluated for the optimal (kCorrOpt), and near-optimal parameters (0.1 kCorr Opt and 10 kCorr Opt). The latter was to verify the stability of the results. For simplicity, we only report the localization performance using optimal parameters in the final round. Figure 4 shows the accuracy of the localized simulations, and Figure  S1 shows the normalized maximum localization error (Eq. 11).
To evaluate the effect of IED, Noise, array Type, and the presence of Overlaps on the normalized median localization error dLoc Med (Eq. 10) we performed a multiple linear regression using the following model  Tables 5 and 6 show the N-way ANOVA results for these regression models.

Figure S1. GridFit maximum localization error of simulated CT artifacts
3. CEPA brain-shift compensation supplementary material

Figure S2. Projection error comparison of Brain-shift compensation methods
Comparison of projection error of Normal-Grid, Normal-SCE, and Springs methods with CEPA. For strip electrodes, the reference location was obtained as the mean between the Normal-SCE and Springs method results. Therefore, the errors for these methods are the same. The normal-Grid method is not available for strips. Projection error statistics for strips and their comparison with CEPA. Z-values were obtained from the Wilcoxon Signed Rank Test. Results indicate that CEPA significantly outperformed the alternative methods. Moreover, the projection of GridFit localized CT artifacts outperformed the Visual one. Note: Normal-Grid method is not available for strips. Therefore, the reference location was obtained as the mean between the Normal-SCE and Springs method results, and the errors and statistical results for these methods were the same.

Comparison of projection roughness of Normal-Grid, Normal-SCE, and Springs methods with CEPA.
All comparisons produced significant statistical differences. Data were obtained from 40 grids and 59 strips.

Comparison of projection roughness for GridFit and visually localized CT artifacts. It can be inferred from the above that visual localization variance in the inter-electrode distance propagates to the back-projected coordinates as reduced roughness. Median projection errors were obtained from 40 grids and 59 strips. Error bars denote 95% CI of the median obtained by bootstrapping. Note that statistical tests (Supp. Tables 9 and 10) were done on paired samples and error bars in the plots do not depict errors associated with those tests.
Supplementary Table 9. Projection roughness statistics for grids.

Normal-Grid Normal-SCE Springs
Normal-Grid

Projection distance vs. local deformation
We studied the local deformations of grids and strips in relationship to the projection distance. More precisely, local deformations were defined as the relationship between the median distance of an electrode ej with its neighbors ek normalized by the IED, i. e., median({||ej -ek||}) where k index neighbor electrodes. The projection distance was computed from the localized CT artifact e0j to the back-projected coordinate ej, i.e., ||ej -e0j||, normalized by the distance to the brain center. The effect of the distance was evaluated with linear mixed-effects models (LME) of the deformations, with fixed effects for distance and random effects for distance (and intercept in Springs and CEPA) grouped by array and patient. Specifically:

Deformation ~ Distance + (Distance | Patient:Array)
for Normal-Grid and Normal-SCE projection methods, and

Deformation ~ 1 + Distance + (1 | Patient:Array) + (Distance | Patient:Array)
for Springs and CEPA projection methods. Model parameters were estimated using the Maximum Likelihood Estimation method. Including an intercept factor in the last two LME models (Springs and CEPA) considers that local deformations could occur even when the projection distance is zero. This is because each projected coordinate can be locally deformed by the complete set of projections. In contrast, the first two methods (Normal-Grid and Normal-SCE) should not produce this effect. Data outliers from each regression were detected using Grubbs's test implemented in Matlab's function "isoutlier". Outliers from all models Springs 38,and CEPA 19) were combined (N = 114) and removed to perform the final data regressions (out of N = 2327). Fewer samples were available for the Normal-Grid method since it is not applicable for strips. Only data obtained from localizing CT artifacts with the GridFit algorithm were analyzed for simplicity.
T-tests were performed to determine the significance of fixed effects, and Likelihood Ratio Tests (LRT) were used to compare the LME models. LRT showed negative results for including correlated random effects in the second model. LRT showed positive results when models were compared against simpler linear regressions with no random effects. Along the same line, the confidence intervals of the random effects excluded zero. It can be observed in Figure 6E that all methods except Springs were sensitive to the projection distance. Normal-Grid showed the strongest effect of distance over deformation, followed by CEPA and Normal-SCE. Moreover, variable levels of dispersion were observed for the different projections as indicated by the adjusted R 2 . Supplementary Table 11 shows modeling details. Supplementary

Alternative reference location
In the main article, we described the localization errors in relation to a mean location reference, i.e., the mean of Normal-SCE, Normal-Grid, and Springs methods' results. Even if this is a parsimonious approach, it might not reflect the best possible back-projection. Here, we explored the possibility that any of the evaluated methods could be more reliable as a reference. Therefore, each method was defined as a reference, and the errors were computed against these locations. The figure below (Fig S7) shows the median errors of using each method as a reference location. Observe that in all situations but one (Springs-ref, CEPA vs. Normal-Grid), the CEPA method statistically outperformed the other evaluated methods. These results suggest that in the worst scenario assessed here, the CEPA method produced the second-best localization alternative, with a median error below 1.3mm in all cases. See statistical details in Supplementary Tables 12  and 13.

Figure S7. Projection error when selecting other brain-shift compensation methods as a reference.
In all situations but one (Springs-ref, CEPA vs. Normal-Grid)