1 Introduction

The deep brain stimulation (DBS) is a well established treatment method for late-stage Parkinson’s disease (PD), essential tremor, dystonia and other movement disorders. It consists of surgical placement of a permanent stimulation electrode into subcortical structures and chronic electrical stimulation using a stimulator implanted commonly in the chest cavity. In order to achieve a high level of symptom suppression with low side-effects, a highly accurate positioning of the stimulation contact is necessary, yet challenging. In case of the DBS for PD, which is the main focus of this study, the most common target – the subthalamic nucleus (STN) – is around 10 mm long along its longest axis and has a relatively low contrast in pre-operative MRI scans (see Fig. 1). Moreover, the optimal stimulation target is even smaller and lies in the dorsolareral motory subregion of the STN.

A typical implantation procedure starts with a pre-operative MRI scanning, which is used for target nucleus localization and surgical planning. In order to mitigate brain shift and other inaccuracies, occurring during the surgery, most surgical teams then employ intra-operative microelectrode recording (MER) of electrophysiological activity in the vicinity of the planned target, using typically up to five parallel microelectrodes. In clinical practice, the individual MER signals are evaluated manually by a neurologist and the target nucleus is identified based on a characteristic firing pattern.

Over previous years, researchers have suggested several automatic classification methods for the MER signals, based most commonly on signal power and spectral properties of the MER, some of which got recently included into clinical software tools for microexploration [1, 2]. Despite the apparent benefits these methods may have for implantation efficiency, they provide no spatial mapping of the electrophysiological findings or explicit MER localization within the nucleus, necessary for both clinical and research applications.

In our recent study [3], we presented a probabilistic model, which allows mapping of an anatomical STN atlas to the recorded multi-electrode MER directly and thus provides MER classification and localization at the same time. However, due to the inherent anisotropy and low spatial distribution of the MER (we used a common “Ben-gun” setting with 5 parallel MER trajectories, spaced 2 mm apart in a cruciform configuration, with signals recorded at steps of 0.5 mm), the MER data provide accurate information about size of the STN along the axis of the electrodes around the planned target but provide substantially less information about the shape of the STN in other anatomical directions.

In this paper, we investigate the possibilities of fusion of our previous model with additional information obtained from the pre-operative MRI imagery, by combining atlas rotation and scaling based on pre-operative MRI landmarks with additional position refinement and brain shift estimation using the MER data. We validate the properties of the extended models on a set of 27 multi-electrode trajectories from 15 PD patients. As both approaches are not without limitations, we also outline the possibility to perform a complete fusion of MRI and MER data to estimate the patient-specific stn shape, as well as the brain shift directly at the same time.

Fig. 1.
figure 1

Illustration of STN size and contrast in an axial slice of pre-operative T2-weighted MRI image: STN contour (red, right hemisphere) and characteristic hypointensity (green circle, left hemisphere) (Color figure online)

Full size image

2 Methods

2.1 Common Definitions

Throughout this text, we use transformation vector \(\varvec{r}\) with 9 degrees of freedom to transform a 3D surface-based atlas into patient-specific coordinates.

$$\begin{aligned} \varvec{r} = [ \varvec{t}, \varvec{s}, \varvec{\gamma }], \end{aligned}$$
(1)

where \(\varvec{t}=[t_x, t_y, t_z]\) is the translation, \(\varvec{s}=[s_x, s_y, s_z]\) scaling and \(\varvec{\gamma }=[\gamma _x, \gamma _y, \gamma _z]\) rotation along/around the x (medial \(\rightarrow \) lateral), y (posterior \(\rightarrow \) anterior) and z (ventral \(\rightarrow \) dorsal) axis. We use the 3D STN atlas from [4] in a form of standard 3D triangular mesh but any surface-based STN atlas can be used as well. As a reference, we use a set of 12 characteristic STN landmark points (plus the anterior and posterior commissure: AC and PC) as in [5], which were identified by an experienced neurologist on the atlas, as well as on pre-operative MRI data of each patient.

The MER recordings are represented in the feature vector \(\varvec{x} = \{x_1,...,x_N\}\), recorded at corresponding spatial locations \(\varvec{L} = \{\varvec{l_1},...,\varvec{l_N}\}\). The vector \(\varvec{x}\) consists of a single feature, the normalized signal root-mean-square of the whole MER signal (NRMS) as in [6]. For the purposes of validation, we use manual annotation of each MER signal as STN or non-STN, done by an experienced neurologist during the surgery.

2.2 Imaging-Only Method (allPoints)

As a reference, we use a method based solely on pre-operative MRI data and STN landmark points annotated therein, the allPoints. This method uses 12 characteristic landmarks on the STN boundaries, defined previously in [5] and coordinates of the anterior and posterior commissure. The method then finds a full 9-DOF transformation to minimize the least-square distance between the characteristic points on the atlas and in given patient’s manually annotated MRI data.

2.3 Basic Electrophysiology-Only Model (nrmsCon)

This model forms the basis for the extended models below and has been presented in our previous paper [3]. Simply put, the model shifts, scales and rotates and the 3D atlas around the MER recording sites in a way that the high NRMS values are encapsulated in the STN atlas volume and the low NRMS values are excluded (owing to the higher neuron density and thus higher NRMS values inside of the STN).

In more formal terms, model assumes different distribution of NRMS values observed inside and outside of the STN (emission probabilities, modeled using separate log-normal distributions), and fuzzy boundaries of the STN atlas (membership probabilities modeled using a logistic function). These parameters form together the parameter vector \(\varvec{\varTheta }\), which is estimated from training data. In order to fit the atlas to MER recordings of a particular patient (NRMS values \(\varvec{x}\) measured at locations \(\varvec{L}\)), the model finds parameters \(\varvec{r}^*\), which maximize the likelihood, defined as:

$$\begin{aligned} \varvec{r}^*\ = \underset{\varvec{r}}{arg \max } ~ \mathcal {L}(\varvec{r}|\{\varvec{x},\varvec{L}\},\varvec{\varTheta }) = \underset{\varvec{r}}{arg \max }~p(\{\varvec{x},\varvec{L}\}|\varvec{r},\varvec{\varTheta }) \end{aligned}$$
(2)

where the probability of a single observation \(\{x_i,\varvec{l_i}\}\) being in state s is given by the product of the emission probability and membership probabilities

$$\begin{aligned} p(\{x_i, \varvec{l_i} \in s \}|\varvec{r},\varvec{\varTheta }) = p(x_i |\varvec{l_i} \in s,\varvec{r},\varvec{\varTheta }) \cdot p(\varvec{l_i} \in s |\varvec{r},\varvec{\varTheta }) \end{aligned}$$
(3)

The joint probability for a single observation is then computed as a summation over both states possible states (INside and OUTside the STN):

$$\begin{aligned} p(\{x_i,\varvec{l_i}\}|\varvec{r},\varvec{\varTheta }) = p(\{x_i, \varvec{l_i} \in IN \}|\varvec{r},\varvec{\varTheta }) + p(\{x_i, \varvec{l_i} \in OUT \}|\varvec{r},\varvec{\varTheta }) \end{aligned}$$
(4)

To compute the joint probability of the whole observation sequence of N MER, we naïvely assume conditional independence given model parameters and compute the joint probability as:

$$\begin{aligned} p(\{\varvec{x},\varvec{L}\}|\varvec{r},\varvec{\varTheta }) = \prod _{i=1}^{N} p(\{x_i,\varvec{l_i}\}|\varvec{r},\varvec{\varTheta }) \end{aligned}$$
(5)

The maximum shift is constrained to \(\pm 5\) mm in any direction, maximum scaling to \(\pm 25\%\) in each direction and rotation maximum \(\pm 15^\circ \) around each axis, the model is thus abbreviated nrmsCon. For more details on model structure and fitting, please refer to the aforementioned publication [3] or the thesis [7].

2.4 The Proposed Combined Model (nrmsBrainShift)

We introduce the following way to fuse the electrophysiology-based model with prior information about STN size and rotation, contained in the pre-operative MRI landmarks: In the first step, the allPoints landmark-based transformation is used to compute atlas scaling and rotation. Subsequently, the MER-based model described above is used to estimate the translation parameters \(\varvec{t}\), and thus to estimate the brain shift. The model is also capable of additional modification of the scaling and rotation parameters, which are regularized. The probability density function for all observations from Eq. (5), is modified as follows:

$$\begin{aligned} p(\{\varvec{x},\varvec{L}\}|\varvec{r},\hat{\varvec{\varTheta }}) = \prod _{i=1}^{N} p(\{x_i,\varvec{l_i}\}|\varvec{r},\hat{\varvec{\varTheta }}) \cdot \prod _{m\in \{\varvec{s},\varvec{\gamma }\}} p(r_m|\hat{\varvec{\varTheta }}) ^w, \end{aligned}$$
(6)

where the additional term \(p(r_m|\hat{\varvec{\varTheta }})\) penalizes deviation from the initial allPoints scaling and rotation, using likelihood of the normal distribution (\(p(r_m|\hat{\varvec{\varTheta }}) = 1/{\sqrt{2\pi \sigma _m^2}} \cdot exp(-\frac{(x-r_m)^2}{2\sigma _m^2})\)), centered at the initial value of given parameter \(r_m\), with standard deviation \(\sigma _m\) estimated from the training data and stored in the extended parameter vector \(\hat{\varvec{\varTheta }}\). The exponent w represents a weighting parameter, which can be used to set the trade-off between MER-based (\(w\rightarrow 0\)) and MRI-based (\(w \rightarrow \infty \)) fitting. We evaluated the results for \(w \in \{0,0.01,0.025,0.05,0.1,0.25,0.5,1\}\).

2.5 Performance Evaluation

In order to estimate the out of sample performance of the proposed method and due to the relatively small sample size (in terms of whole patient sets), we employed the leave one subject out (LOSO) procedure. In each iteration we kept one subject’s data (maximum two 5-electrode trajectories for bi-laterally implanted patients) for model fitting and evaluation, while all other data were used to obtain the parameters \(\varvec{\varTheta }\).

To compute performance metrics, we use two approaches:

  1. (i)

    Machine-learning metrics where we count the number of STN MER recordings (according to expert MER labels), correctly encapsulated in the atlas volume at the final position (true positives), or falsely excluded from the atlas volume (false negatives). True negatives and false positives are computed analogously from the non-STN labeled MERs. Standard performance measures are calculated: sensitivity, specificity and accuracy.

  2. (ii)

    Evaluation of transformation parameters, obtained from the tested model, compared to least-squares transformation of the atlas to the STN landmark points in the pre-operative MRI data of given patient (see the allPoints method below). Here, we assume that the pre-operative data provides accurate information about the rotation and scaling of the atlas, but does not provide a good estimate of the translation vector \(\varvec{t}\) due to the non-negligible brain-shift.

3 Results and Discussion

3.1 Collected Data

For validation of the method, we use a dataset from 27 explorations in 15 PD patients with complete 3D information and another 9 explorations from 4 patients without information on spatial recording locations, used in the training phase to estimate \(\varvec{\varTheta }\) (or \(\hat{\varvec{\varTheta }}\)) only. But was excluded from validation. Altogether, we used 175 electrode trajectories and 4538 recorded MER signals from 19 PD patients.

Table 1. Classification results (LOSO validation-set)
Full size table
Fig. 2.
figure 2

Evaluation of the dependency of the proposed nrmsBrainShift method on the weighting parameter sigma, by computing Pearson’s correlation coefficient for each transformation parameter with the reference allPoints method for varying values of the weighting coefficient w. Note that the translation parameters \(t_x\), \(t_y\),\(t_z\) are not penalized and are thus unaffected by the value of w.

Full size image
Fig. 3.
figure 3

Examples of model fit using the proposed nrmsBrainShift method (left), fusing electrophysiology data with MRI landmarks, the MRI-landmark-only allPoints method (center) and the electrophysiology-only nrmsCon (right) on data from patient No. 5 (right STN). The final STN atlas position after fitting is shown in purple, width of the five microelectrode trajectory cylinders denotes the NRMS value, while colors denote manual labels: STN in yellow, non-STN in grey. MER positions inside the resulting model are denoted by black points, planned target by red o. The nrmsBrainShift method provides an anatomically more reasonable fit at the cost of slightly lower accuracy. (Color figure online)

Full size image

3.2 Fitting Results

The classification results on validation data is shown in Table 1. Although the electrophysiology-based nrmsCon achieves much better fit to the electrophysiology data, than the MRI-based allPoints method, there was almost no correlation of the scaling and rotation parameters with the reference allPoints method. As we assume the prevailing source of inaccuracy during surgery to be due to displacement, rather than deformation, we expected the only differences in terms of shift/translation. On the presented dataset, the nrmsCon model achieved the pre-set optimization constraints (min/max scaling and rotation) in 47 cases, which accounted for more than 30% of the fits (and 19 out of 27 trajectories). The method apparently leads to overfitting, providing good classification at the cost of diverging from anatomically reasonable transformation.

In contrast, the newly proposed combined nrmsBrainShift technique achieved only slightly lower classification accuracy, while maintaining reasonable transformation - a fact illustrated also in the Fig. 2.

The impact of the weighting coefficient on divergence from anatomically relevant location is illustrated in Fig. 3. While high values of w lead to a highly constrained fit, where the only changes are in the translation parameters, low values of w lead to more flexible fit to the MER data at the cost of lower veracity of the transformation. Comparing the properties to the Table 1, it is clear that increasing w towards one leads to only a minor drop in classification accuracy, and a more marked drop in sensitivity. Overall, the sensitivity is the most problematic parameter for all methods, which is likely due to the inability of the model to adapt more flexibly to patient-specific STN shapes.

4 Conclusion

While the previously published electrophysiology-only nrmsCon model [3] proved electrophysiology-based fitting feasible, a subsequent detailed investigation revealed strong overfitting with too harsh model transformation, leading to unlikely results.

Fortunately, the proposed model using the pre-operative landmarks to initialize (and potentially constrain) the fitting, achieves comparable accuracy - i.e. ability to correctly contain STN-labeled MER locations - while maintaining anatomically accurate scaling and rotation. The main drawback of the method is thus in the necessity to identify the 12 landmark points in pre-operative data. We believe, that similar probabilistic framework could be used for direct automatic fusion of pre-operative MRI data, which would eliminate the need for the manual landmark labelling and increase the ability of the model to adapt to inter-individual differences in STN shape.

Overall, the fusion of pre-operative MRI data with electrophysiology provides a promising option for increasing accuracy of electrode localization both intra-operatively, as well as during offline evaluation in research studies on DBS mechanisms and STN physiology.