DBS electrode localization and rotational orientation detection using SQUID-based magnetoencephalography

The cylindrical head phantom used in this study is depicted in figure 1(a). It was slightly changed and adapted for this work compared to our earlier studies in [10, 11]. The cylinder body with dimensions comparable with a human head (diameter of 150 mm, height of 250 mm) is made of acrylic glass and is, therefore, neither electrically conductive nor magnetic. It was filled with an isotonic fluid (NaCl 0.9%) to mimic the electric conductivity of a human brain. We used a DBS system by Boston Scientific consisting of a current-controlled DBS neurostimulator (Vercise^™ PC) and a directional electrode (Versice Cartesia^™) as used in the clinical routine. The DBS system was integrated and fixed into the phantom. The DBS electrode can be moved along a circle in the xy-plane by rotating the adjustment wheel and was placed at five different positions; each with a different distance 'a' to the phantom's origin (40, 30, 20, 10, 0 mm). a = 0 mm represents the deepest position at the origin of the phantom, whereas a = 40 mm represents the position nearest to the phantom surface (phantom radius is 75 mm). The vertical position of the electrode remained fixed during these measurements. The electrode is connected to the neurostimulator via an extension wire and connector. In most previous works, the magnetic field around the DBS electrodes was severely contaminated by the presence of artifacts related to ferromagnetic extension wires. These wires were made of stainless steel and generated artifacts in the MEG signals [16–18]. We used non-ferromagnetic extension wires made of platinum-iridium to avoid such artifact [19]. Moreover, we used non-magnetic material for all other components used in the phantom. The screws on top of the electrode holder were made of titanium, while the adjustment wheel and the holders for the electrode and neurostimulator casing were made of plastic. Unintentionally, the ferrule under the titanium screw was made of ferromagnetic material and caused artifact at the corresponding sensor points. The structure of the directional DBS electrode is shown in figure 1(b). It comprised eight individually controlled platinum-iridium contacts (C1–C8), in which the two middle contact levels were segmented into three contacts each spanning 120^∘ of the circumference. Any combination of these contacts can be activated to steer the stimulation current direction. The outer jacket is made of polyurethane. The dimensions of the electrode are as follows: diameter of 1.3 mm, contact length of 1.5 mm, and a vertical spacing between adjacent contacts of 0.5 mm.

Zoom In Zoom Out Reset image size

Before starting the actual MEG recording, several preparatory steps were performed. Four head position indicator (HPI) coils were attached to the surface of the head phantom and a round template with three landmarks (referred as nasion, right pre-auricular, and left pre-auricular) was fixed on the cylinder cap. Then, a point digitizer (Fastrak, Polhemus) was used to localize anatomical fiducial points (nasion, left and right pre-auricular) and HPI coils. Furthermore, 150 additional points around the phantom were localized to improve the accuracy of the spanned coordinate system of the phantom. Figure 2(a) shows the phantom during the digitization process; the glasses represent the reference, and one of the two receivers used for digitization. The digitized locations of the anatomical landmarks in [x, y, z] directions are as follows:

• Nasion: [0, 74.9 mm, 0]
• Right pre-auricular: [74.7 mm, 0, 0]
• Left pre-auricular: [−75.1 mm, 0, 0]

Zoom In Zoom Out Reset image size

The cap of the cylindrical phantom lied exactly on the xy-plane, the longitudinal axis of the cylinder on the z-axis, and the origin of the coordinate system was in the center of the cylinder's circular cap. The radius of the cylindrical phantom was 75 mm, and the digitization process was performed with an accuracy of ±0.3 mm (given by the error of right pre-auricular). The locations of HPI coils, landmarks, additional points, and the position and orientation of the electrode were all expressed in the phantom coordinate frame. The phantom placement in the MEG scanner can be seen in figure 2(b). The MEG sensors were in the MEG coordinate system. Before each recording, MEG sensor positions were automatically determined in phantom coordinates. The relationship between MEG and phantom coordinates was determined by applying current with a particular frequency to the HPI coils and determining their positions with respect to the MEG sensor array. This step was performed by the MEG system itself, which already had the positions of the HPI coils in phantom coordinates. Figure 2(c) illustrates the positioning of the MEG helmet in the phantom coordinate system from side view and shows all digitized points. The coordinates of all five electrode positions were therefore known in advance, and only the position of the MEG sensor array with respect to the phantom was to be determined by the MEG system before the first measurement. Figure 2(d) shows the top view of the phantom with respect to the MEG sensor array and illustrates the orientation of each segmented electrode contact (C2–C7) at electrode position a = 40 mm.

The electrical stimulation comprised of biphasic rectangular pulses with a stimulus pulse phase and a passive charge-balancing phase. The parameters of stimulation used are mentioned in table 1. MEG data was collected with an Elekta Neuromag VectorView® MEG scanner at the Universitätsklinikum Düsseldorf (UKD) in Düsseldorf. The scanner comprises of 306 individual channels, corresponding to 102 magnetometers and 204 gradiometers. It can be noted that only the magnetometer data was used in this study. The noise level according to the datasheet is about 3 fT $({\sqrt{\mathrm{Hz}}})^{-1}$ for both magnetometers and gradiometers. The acquisition parameters of the MEG scanner were set as follows: ***The sampling rate was set to its maximum of 5 kHz. The low-pass filter was set to the highest possible value of 1660 Hz and the cut-off frequency of the high-pass filter was set to direct current (no high-pass filtering) in order to obtain measurements of the DBS signal with the maximum allowed acquisition bandwidth. The duration of each measurement was three minutes. A total of 38 measurements with bipolar electrode configurations were taken (see table 1), with different activated electrode contacts (1–8: contact numbers, first number: anode, second number: cathode), stimulation parameters (amplitude, pulse width, frequency), and electrode positions within the phantom. Before each recording of the first 29 measurements, the phantom was removed from the MEG scanner, the position of the electrode was changed by rotating the adjustment wheel and the phantom was placed back again. It was not possible to place the phantom at the exact same position in the MEG every time, so the estimation of the phantom position in the scanner had to be performed each time by the MEG system, as described in detail in the previous section. A direct comparison between the localization and rotational orientation results across these 29 measurements can therefore not be drawn. However, in the 30th measurement, the electrode position was set once to a = 40 mm and was not changed until the end. Here, the obtained results could be compared with each other.

For the sake of clarity, the table is divided into five groups. In the first group, the bipolar electrode configuration with the maximum possible stimulation parameters was set. In the second group, measurements were performed with the activation of contacts at different electrode heights and positions. The third and fourth group provide directional stimulation with the activation of the tip of the electrode (contact C1) against an individual segmented contact (C2–C4). Notice that the phantom was not moved during the measurements in the fourth group and afterwards. In the fifth group, the stimulation amplitude was reduced to a value as used in clinical routine, i.e. 3 or 1.5 mA. The good quality of recorded data was ensured by an empty room measurement prior to the start of the experiment and by visual inspection of about 1 min of MEG recording before each measurement. Electrode localization was performed for measurements with non-directional stimulation (groups 1, 2, and 5), and rotational orientation detection was performed for measurements with directional stimulation (groups 3 and 4).

The processing steps of the MEG measurements are depicted in figure 3. For each measurement, we obtained 102 time signals, one from each MEG sensor. The measured data was imported into MATLAB® (Version R2018a) and preprocessed using the FieldTrip toolbox [20]. Each time signal was 180 s long and high-pass filtered (6th order Butterworth) with a cutoff frequency of 60 Hz without signal loss since the signal contains higher frequencies (stimulation frequency f₀ and its harmonics) as shown in table 1. Each signal was divided into short segments of length equal to the inverse of stimulation frequency, e.g. for 170 Hz stimulation $\mathrm{T}_s = 1/f_s = 5.88$ ms (first five measurements in table 1) and $\mathrm{T}_s = 7,69$ ms for 130 Hz of stimulation (remaining measurements in table 1). These segments were then averaged which improved the signal-to-noise ratio (SNR) by $\sqrt{N}$, where $N = 180~\mathrm{s} \cdot f_0$. The maximum value from this averaged time segment is then taken for each MEG sensor, resulting in a total of 102 values representing the maximum magnetic field $\mathrm{B}_{\mathrm{meas}, i}$ measured by each sensor at location i. Since these were phantom measurements, the measured MEG signals were not contaminated by biological artifacts such as cardiac muscle, skeletal muscles, or eye movements. Thus, the measurements with a phantom represent an ideal case. The ambient noise solely consisted of the power line interference, which only affected non-important frequencies (50 Hz and harmonics). The only magnetic noise in our MEG recordings arose from a non-biological source, i.e. from the ferromagnetic ferrule, which was placed under the screw, that over time became magnetized and generated an artificial magnetic field. This artifact was modeled in the electromagnetic simulations using the known position and dimensions of the ferrule and was subtracted from the measured data. The results are given in section 3.2.

Zoom In Zoom Out Reset image size

Finite element modeling (FEM) and model-fitting approach were used to determine location and rotational orientation of the DBS electrode in the phantom. Using FEM, the magnetic field was calculated at measurement points i with different electrode positions p and rotations r. The values calculated by the model are represented by $\mathrm{B}_{\mathrm{model}, i}(p,r)$. That is, a single model with a given electrode position p₁ and electrode rotation r₁ resulted in 102 modeled values corresponding to a total of 102 measurement points. Modeled data with different electrode positions was used for electrode localization that was performed for measurements with non-directional stimulation (groups 1, 2, and 5 in table 1). The corresponding contacts of the electrode were activated in the model resulting in a symmetrical flow of current along the length of the electrode. Modeled data with different electrode rotations was used to perform orientation detection using measurements with directional stimulation (groups 3 and 4 in table 1). The bottom contact C1 (the tip) of the electrode was activated against one of the segmented contacts C2, C3, or C4. This led to a current flow that was steered in a direction compared to the non-directional stimulation. The model is described in detail in the next paragraph. A fitting approach was then used that searched for the model that best fit the measured and preprocessed data $\mathrm{B}_{\mathrm{meas}, i}$, which was described in the previous subsection. Therefor, the algorithm calculated the normalized root mean square error (NRMSE) between each modeled and measured data, and determined the model with the smallest error. From the selected model $\mathrm{B}_{\mathrm{model}, i}(p_{\textrm{best}},r_{\textrm{best}})$, $p_\textrm{best}$ was the computed electrode position and $r_\textrm{best}$ the computed electrode orientation. The error in electrode localization was calculated with the Euclidean distance between $p_\textrm{best}$ and $p_\textrm{real}$, the real known position of the electrode in the phantom. The error in determining the electrode orientation was calculated using the angle between $r_\textrm{best}$ and $r_\textrm{real}$, the real known orientation of the electrode in the phantom.

FEM electromagnetic simulations were based on a model created in COMSOL Multiphysics®, as shown in figure 4. For each structure in the model, the corresponding material was taken directly from the built-in material library and value of a uniform electrical conductivity σ, relative permeability µ_r , and relative permittivity _r were assigned. The electrode was placed parallel to the longitudinal axis of the cylinder. The cylinder body and the outer jacket of the electrode were modeled with acrylic plastic ($\sigma \approx 0~\mathrm{S\ m^{-1}}$, µ_r = 1, _r = 4.2), the eight independent electrode contacts with platinum-iridium ($\sigma = 10^6~\mathrm{S\ m^{-1}}$, µ_r = 1, _r = 1), and the content of the phantom with saline solution ($\sigma = 1.6~\mathrm{S\ m^{-1}}$, µ_r = 1, _r = 80). These values were taken from the literature. The phantom including the modeled electrode was surrounded by a cuboid box with a length of 500 mm that is filled with air ($\sigma \approx 0~\mathrm{S\ m^{-1}}$, µ_r = 1, _r = 1). The MEG helmet was modeled as a surface and all measurement points were located on that surface. Tetrahedral mesh of 'Extra Fine' resolution was utilized. To model the bipolar electrode configurations, the 'Terminal' boundary condition of 'Current Type' was used for the cathode and the 'Ground' boundary condition was used for the anode. Simulations were performed using the 'Magnetic and Electric Fields (mef)' interface of the 'AC/DC module' of COMSOL and magnetic flux densities in x, y, and z direction ($\mathrm{B}_x$, $\mathrm{B}_y$, and $\mathrm{B}_z$) on the surface of the MEG helmet were calculated. In this experiment, SQUID sensors (fixed in the helmet) measured the magnetic field perpendicular to the helmet surface. To find the magnetic field $\mathrm{B}_{\mathrm{model}, i}$ perpendicular to the helmet surface in the simulation, we used the dot product $\mathrm{B}_{\mathrm{model}, i} = \vec{\mathrm{B}}_{i}\cdot\vec{\mathrm{s}}(i)$, where $\vec{\mathrm{s}}(i)$ is the unit vector perpendicular to the helmet at the sensor location i and $\mathrm{B}_{i} = \mathrm{B}_x + \mathrm{B}_y + \mathrm{B}_z$.

Zoom In Zoom Out Reset image size

In the model, the electrode was successively moved in the phantom with a step size of 1 mm (in all x, y, and z directions) and the magnetic field was calculated at each position. This data was used to analyze measurements with non-directional stimulation. The fitting approach found the electrode position which resulted in minimum NRMSE with the measured data. To analyze directional stimulation measurements, the electrode was successively rotated along its own axis with an angle step of 1^∘ and the produced magnetic field was found for each step. Again, the fitting approach found the orientation that minimized the NRMSE. To improve detection accuracy, contiguous measurements could be evaluated together using a-priori information of electrode geometry. When determining electrode position, measurements with different levels of activated contacts along the electrode (lower, middle, upper level stimulation) were used, and the distance between adjacent levels was known to be 2 mm. When determining electrode rotation, measurements with different directional segments were used, and the angle between adjacent segments was known to be 120^∘. In contrast to the analysis of a single measurement, where the fitting approach determined only a single model which gave minimum NRMSE, when analyzing contiguous measurements (i.e. three measurements with stimulation electrode at different level or direction), the fitting algorithm determined three contiguous models (i.e. three related models with stimulation level separated by 2 mm or direction separated by 120^∘) whose calculated values minimized the sum of the individual NRMSE with the corresponding measured values. The results based on consideration of multiple measurements are presented in section 3.5.

The magnetic flux densities obtained from the time signals are summarized in table 2. The maximum amplitudes from all sensors from all measurements in groups 2 and 5 were considered. From the data shown in the table, we can conclude the following two points:

• Decreasing the stimulation amplitude leads to a proportional decrease in the magnetic field.
• The expected mean magnetic fields from MEG measurements in DBS patients with bipolar stimulation amplitude of 1.5–3 mA will be in average in the range of 1–2 pT.

The magnetic field distribution for measurement #11 (non-directional) can be seen in figure 5(a). Since we had data for only 102 points across the surface, we used interpolation to get a surface plot. This figure represents a typical measured field distribution for bipolar recordings. It shows the magnetic field that is generated by the current flowing between the activated contacts (see figure 5(b)) and the undesired field that is created due to the magnetization of the ferromagnetic ferrule (see figure 5(c)). Our phantom was made of non-magnetic materials (titanium screws and plastic parts). As already mentioned in section 2.4, this ferrule described the only magnetic noise which produced an artificial magnetic field. This artifact was observed in each measurement and was removed by modeling (in electromagnetic simulations) and by subtracting it from the measured data. Figure 6 compares the electromagnetic simulation results of a non-directional (left column) and directional (right column) stimulation. Figures 6(a) and (b) displays the top and side view of the electrode with the current density generated due to the active contacts. In the non-directional case, the current flows symmetrically along the electrode and in the directional case, it is asymmetrical. The asymmetry in the current flow also causes an asymmetry in the magnetic field distribution (see right side of figure 6(c)) which determines the rotational orientation of the electrode. We did not observe an effect of the connector, excess electrode wire, and neurostimulator on the measured magnetic field. This may be because the return current in the connector and wire cancels the field generated by the current flowing in the other direction (bipolar configuration). The current flowing within the neurostimulator may not be symmetrical and could have generated some magnetic field, but this is difficult to isolate and measure.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Before the results of electrode localization and orientation determination are presented in the following sections, this section investigates the overall measurement accuracy since it describes the spatial resolution of the entire system. The system measurement accuracy depended on the inaccuracy of the following measuring steps:

• The phantom was constructed in-house and had an accuracy in sub-millimeter range (±0.1 mm).
• The digitizer device showed an accuracy of ±0.3 mm (see section 2.2).
• The position of the HPI coils on the phantom surface, localized internally by the MEG system, provided the accuracy of the initial phantom position determination. The distances between the fitted coil positions and the digitized points are compared in figure 7, which shows the HPI fitting errors for all performed measurements. The average error was ±1.6 mm for the pair 1–3, 1.4 mm for the pair 1–4, and 0.6 mm for the pair 3–4. The pairs with coil number 2 had larger fitting errors and were therefore not considered by the MEG system.

Zoom In Zoom Out Reset image size

The sum of all inaccuracies leads to a value of ±2 mm and describes the spatial resolution of the entire system. This value is taken into account in the electrode localization and rotation determination in the next sections. The dominant error is caused by the estimation of the phantom position by the MEG scanner.

In this section, the results of electrode localization and rotational orientation determination from the first three groups in table 1 are presented. Note that after each measurement in these groups, the phantom was taken out of the MEG scanner to reposition the electrode and could not be placed exactly in the same position with respect to the MEG scanner and, thus, the measurements conducted in these groups cannot be compared. Therefore, each measurement had to be analyzed separately.

3.4.1. Localization

Electrode localization was performed for the first 20 measurements in table 1. In these measurements with non-directional stimulation electrode position (in the xy-plane), the level of activated contacts along the z-axis was varied. The real electrode locations and the localization results from the measured magnetic data are depicted in figure 8. It shows the top view of the cylindrical phantom which has a radius of 75 mm. Five blue dots represent the five electrode placements (a = 40, 30, 20, 10, 0 mm) used in this work. The locations of the DBS electrode found using MEG are illustrated by red asterisks. This deviation is similar at all electrode positions which suggests that the localization accuracy is not affected by the proximity of the electrode to the phantom surface and to the sensor array. A similar distribution of localization results was observed for activated contacts at different electrode levels. The localization accuracy (represented by the Euclidean distance between real and computed electrode position) is calculated to be around 3 mm on average.

Zoom In Zoom Out Reset image size

3.4.2. Orientation detection

The rotational orientation of the DBS electrode was found for measurements 21–29. Directional current was applied in the phantom in different directions using bipolar stimulation by activating contact C1 as anode and segmented contacts C2, C3, and C4 one-by-one as cathode. Real orientation of each segmented contact at three different electrode positions in the phantom, as well as the rotation results from the measured data are depicted in figure 9. It displays the top view of the phantom, where blue lines represent real orientations of each segment and red arrows orientations estimated by magnetic field measurements. Computed orientations are located close to the corresponding real orientations. The detection accuracy was calculated to be around 15^∘.

Zoom In Zoom Out Reset image size

In this section, the results of electrode localization and rotational orientation determination from the last two measurement groups in table 1 are presented. The position of the phantom remained the same across these measurements since it was not taken out of the MEG scanner to reposition the electrode. Since the distance between adjacent contacts of the electrode (which is 2 mm) and the angle between adjacent segmented contacts (which is 120^∘) is known, measurements at different stimulation heights and directions can be considered together to improve detection accuracy.

3.5.1. Localization

Electrode localization was performed for the measurements 33–35 and 36–38, where the phantom was not moved during measurements. Figure 10(a) depicts the localization results for measurements 33–35 with colored points. Different colors represent different levels of activated contacts along the electrode, where green refers to the lower, red to the middle, and blue refers to the upper level. The positions of localized sources reveal an average localization error of 2.8 mm, which is closer to the accuracy obtained in the previous section. Since the position of each electrode level was estimated separately, the distance between adjacent points in the figure deviated significantly from 2 mm. If the a-priori known distance is applied as a condition for localization as delineated in figure 10(b), the result improves slightly to the value 2.2 mm. Moreover, it can be seen that localized points are not arranged along the electrode but along a different axis. This error is mainly due to the inaccurate estimation of the phantom position in the MEG scanner. When looking at the figure, it already becomes clear that the localization accuracy could be further improved if the localized points would be slightly shifted in negative z-direction and tilted in the xz-plane. This observation suggests that a slight readjustment of the MEG sensor array, i.e. an improvement of the system spatial resolution, would significantly improve the localization accuracy.

Zoom In Zoom Out Reset image size

3.5.2. Orientation detection

The determination of the rotational orientation of the directional electrode was performed for measurements 30–32. The electrode was directionally stimulated in the phantom by activating the tip of the electrode against all segmented contacts successively. The obtained rotation results are shown in figure 11(a) with red arrows, while the blue lines indicate real orientations of the segments. The accuracy of rotation determination is calculated to an average of 15^∘, which corresponds exactly to the accuracy obtained in the previous section for different electrode positions. The orientation for each segment was estimated separately. If the a-priori known angle is applied as a condition for rotation detection, the result improves slightly to a value of 11^∘, as illustrated in figure 11(b). This error is mainly due to the inaccurate estimation of the phantom position in the MEG scanner. A minimization of this error can be achieved by a more precise measuring system.

Zoom In Zoom Out Reset image size