Keywords

1 Introduction

Deep brain stimulation (DBS) is an effective surgical treatment for movement disorders such as Parkinson’s disease (PD). The efficacy of DBS therapy is dependent on the accuracy of surgical targeting and the optimization of stimulation parameters. For the former, complicating the challenge of accurate surgical targeting is brain shift, which occurs from an alteration to intracranial mechanical equilibrium due to the introduction of a burr hole and opening of the dura. Brain shift compromises the alignment between intraoperative patient anatomy and the preoperative magnetic resonance (MR) image volume of the patient, which is used for surgical planning and targeting. Using interventional magnetic resonance (iMR) imaging, Ivan et al. found shift ranging from 0.0 to 10.1 mm throughout the brain with the greatest shift observed in the frontal lobe; 9% of the patients had shift over 2 mm in deep brain structures and 20% over 1 mm [1]. This misalignment between the preoperative presentation and intraoperative anatomy of the patient is particularly problematic for DBS therapy as its targets are often relatively small. Ivan et al. suggested misplacement of the target by as little as 2 mm can introduce inadequate treatment [1]. Current clinical approaches to address brain shift are microelectrode recording (MER) in awake procedures and iMR in asleep procedures. A biomechanical model-based approach that leverages sparse intraoperative measurements, if accurate, could potentially offer another avenue of brain shift correction that could either complement or eliminate the need for MER, or reduce the workflow and economic burdens posed by iMR. Comprehensive efforts of brain shift compensation in DBS surgery via model-based approach are quite limited, especially with respect to validation studies using clinical patient data and possible intraoperative deployment. Early works (e.g. Hamzé et al.) have largely focused on the feasibility of forward solution or preliminary validation with simulated data [2]. In recent work by Luo et al. and followed independently by Li et al., both examined model-based approaches with clinical data: with the former, iMR data was available for comparison, and in the latter only two subsurface targets using preop- and postop- CTs were analyzed; both however, were limited to one patient [3, 4]. The motivation of this work is to refine a model-based approach for brain shift correction that presents minimal disruption to existing workflow for possible intraoperative deployment, and examine model prediction using high fidelity clinical iMR imaging data.

Our emphasis on the feasibility of intraoperative deployment also stems from the desire to further explore the second aspect of a successful DBS therapy, i.e. optimization of programming parameters. To facilitate the optimization of modulation parameters in DBS treatment, previously groups, most notably McIntyre et al., have constructed bioelectric finite element (FE) models to estimate voltage or electric field distribution, which can be used to predict volume of tissue activation (VTA) and its overlap with anatomical structures for surgical planning and target selection [5]. Similarly, we have pursued a bioelectric FE modeling approach based on bioelectric conductive physics; however, unlike previous studies that focus on reconstructing and optimizing stimulation parameters postoperatively, this work aims to establish a framework that bridges the aforementioned biomechanical model with patient-specific bioelectric model with the objective that better surgical navigation and targeting will be facilitated by adjusting for brain shift, and long-term post-implant electrode positions will be enhanced by estimating shifts during intraoperative electrode evaluation. By accounting for both biomechanical and bioelectric effects, the therapeutic and functional impact of DBS therapy can be better planned and predicted. In particular, the impact of shift on neuronal pathway recruitment can be predicted and adjusted without the need of MER assistance or iMR guidance. In the work presented here, we are examining the feasibility of this integrated multi-physics framework of patient-specific biomechanical and bioelectric models and investigating the differences in VTA and neuronal pathway recruitment with and without considering brain shift.

2 Methodology

2.1 Data

Two patients undergoing iMR-guided DBS procedures were studied where preoperative MR and iMR after bilateral DBS implantations were acquired with patient consent and Institutional Review Board approval; details of surgical and imaging procedures can be found in [1]. Examples of acquired data are shown in Fig. 1, noting significant asymmetric shift in frontal lobe, and crosshair indicating the deformed lateral ventricle.

Fig. 1.
figure 1

(a) preopMR, (b) iMR and (c) model updated MR. Better feature agreement exhibited by iMR and model updated MR vs. preopMR, particularly in the frontal lobe. Crosshair indicates agreement at lateral ventricle between model updated MR and iMR vs. preopMR

2.2 Biomechanical Model-Based Approach for Brain Shift Estimation

The developed model-based shift compensation approach relies on the construction of a deformation atlas, which is a pre-computed collection of intraoperative shift possibilities. Subsequently an inverse problem approach driven by sparse intraoperative measurements is employed to provide volumetric shift estimation.

Patient preoperative MR was used to construct a FE mesh incorporating internal structures such as falx, tentorium and brain stem [6]. Displacement and pressure boundary conditions were prescribed as: (1) brain surface above a preset level is given stress free (freely deforming); (2) brain stem is fixed in displacement; (3) the rest of brain surface and dural septa are given slip (tangential movement allowed yet no normal motion) conditions; (4) nodes above a fluid drainage level have a defined pressure reference value and below have a Neumann condition, i.e. no drainage allowed [6, 7]. Moreover, with the observation of ventricular shape change and hemispheric asymmetric deformation, additional boundary conditions were given to the ventricle. Specifically, the lateral ventricle was modeled as a void and further divided into four segments spatially. For the ipsilateral ventricle associated with the largest asymmetric shift, the ventricle was separated into two segments. Each segmental portion was allowed to assume a type 1 pressure condition with 3 different possible nonzero pressure levels considered, where the introduced pressure gradients reflected ranges within ~7.5 mmHg. Given the combinations available, this provided a total of 9 pressure configurations for consideration in our solution distribution. With respect to the contralateral ventricle, while also having two segments, both were given type 1 pressure level of zero. This pressure treatment described the apparent presence of a pressure gradient due to pneumocephalus. Combining this mode of deformation and the modeled effect of asymmetric CSF drainage, two important factors considered to contribute to asymmetric shift in DBS are compensated for. To account for variability in surgery, different CSF drainage levels and modestly varied head configurations relative to gravity were created assuming the patient was in supine position [3]. Finally, with model-driving conditions reflecting the aforementioned configurations defined, a biphasic biomechanical model was used to resolve the volumetric displacement field for each configuration to form an atlas of deformation [6]. To drive the inverse problem, homologous surface and subsurface points were designated on preoperative and iMR images and served as the source for a least squared error fitting process (n = 27 points for case 1 and n = 31 for case 2). It should be noted here that the rationale for the use of such sparse data (rather than the whole iMR data set) is in anticipation of an approach that would not require iMR, i.e. transcranial or burr hole ultrasound. The optimization involves a linear combination of the deformation atlas solutions evaluated at the sparse data points for fit [3, 6]. Once optimized, the same weighting is used with the entire atlas to provide a volumetric brain shift estimation. As the purpose of this study is to examine the model’s ability to recover subsurface shift, the deformed position of the subsurface points was used to quantitatively evaluate performance. This model predicted volumetric displacement field was also used to update the preoperative MR. This model updated MR was then compared to preoperative MR and iMR for qualitative assessment of model performance as well.

2.3 Bioelectric Model for Volume of Tissue Activation Estimation

Once the updated patient MR is obtained from the biomechanical component of the integrated framework, a bioelectric FE model is constructed via Poisson’s equation:

$$ \nabla \cdot \left( { - \sigma \nabla V_{e} } \right) = I $$
(1)

where Ve is the electrical potential (Volts), σ is the conductivity tensor (Siemens/meter), and I is the injected current from an electrical source.

With the insertion path of the electrode leads visible on the iMR image, shown in Fig. 1(b), the trajectory of the electrode leads was determined by identifying two points along the aforementioned path in order to establish a vector. Furthermore, with the most distal end of the insertion visible on iMR assumed to be tip of the electrode lead and a known DBS model (Medtronic 3389) geometry, electrode leads can be reconstructed and incorporated into the patient-specific bioelectric model with mesh refinement around the electrode contacts. An illustration of the bioelectric model with bilateral electrode leads is shown in Fig. 2. In this preliminary study, a monopolar voltage stimulation was tested where one contact (contact 0) was set to −3 V while the outer brain surface was set to the ground. While patient diffusion tensor imaging (DTI) data was not available in this study, an DTI atlas (HCP1065 at 1 mm available via FSL) in MNI space was mapped to patient space via preservation of principal directions through registration of T1 weighted images in patient and MNI spaces via Advanced Normalization Tools (ANTs) [8, 9]. With simulated patient DTI data, the diffusion tensor of each model FE mesh element was computed via interpolation of the diffusion data from the 8 neighboring voxels with respect to the element centroid and then conductivity was determined from diffusion data via a linear relationship and assumed over the element domain [10], thus providing a heterogenous and anisotropic medium. In addition, the effect of tissue encapsulation surrounding the electrode contacts was considered by designating elements within 0.5 mm of the contacts with a conductivity of 0.1 S/m [11]. The electric potential solution was solved via Eq. (1) and the VTA was determined by a threshold of −0.5 V.

Fig. 2.
figure 2

(a) and (b) Bioelectric model built via model updated MR accounting for shift with reconstructed electrode leads. (c) and (d) VTAs (yellow) superimposed onto model updated MR (Color figure online)

2.4 Integrated Framework—Impact of Shift on Neuronal Recruitment

With model estimated VTA in the intraoperative configuration determined and a model predicted brain shift displacement field obtained, active nodes forming the VTA could be mapped to preoperative space using the inverse of the model displacement field. This process provided estimates of surgical target shift; it also allowed for a comparison of VTAs with and without the consideration of brain shift. With active nodes forming the VTAs determined with and without shift consideration, these active nodes could be transformed to a standard MNI space. The transformation was achieved by registration of T1 weighted images in patient and MNI spaces (ICBM152) via ANTs [8]. Once VTAs were mapped to the standard space, a Human Connectome Project (HCP) dMRI population-averaged template (HCP1021 at 1 mm) available via DSI Studio (http://dsi-studio.labsolver.org) was used to perform tractography [12]. As an experiment to understand the shift impact on neuronal pathway recruitment, this transformation enabled the examination of the differences in pathway recruitment due to VTA differences introduced by shift. Tractography was performed in DSI Studio using the default quantitative anisotropy (qa) threshold, angular threshold of 60°, step size 0.5 mm, smoothing parameter 0.2, length 20.0–200.0 mm with a 500,000 seed count [12]. The resulting fiber tracts were compared between the VTAs with and without shift considerations, qualitatively via visual examination, and quantitatively via number, length and volume coverage of recruited tracts.

3 Results

3.1 Brain Shift Compensation Performance

Qualitative assessment of model shift correction was performed by comparing preoperative MR, iMR and model updated MR on corresponding slices. An example is shown in Fig. 1. Model updated MR exhibits better feature agreements with only sparse data (particularly in the frontal lobe) with the iMR data as compared to preoperative MR. The crosshair in Fig. 1(a)–(c) indicate better shift recovery at the lateral ventricle by the model updated MR image. Also of note, the inverse problem that produced the model prediction and subsequent patient MR update were computed in <1 min.

With respect to a quantitative assessment, designated corresponding subsurface features near the trajectory of electrode leads were examined by comparing this intraoperative measurement to its model predicted counterpart. For each patient, 16 features were examined for a total of 32 measurements in gauging model performance. The model reduced shift induced misalignment from 6.71 ± 1.89 to 1.87 ± 0.64 mm (~72.1% correction) in case 1, and from 8.64 ± 1.42 to 2.89 ± 0.92 mm (~66.47%) in case 2. Overall the model reduced misalignment from 7.67 ± 1.92 to 2.38 ± 0.94 mm (~68.94%).

3.2 VTA Estimation

The bioelectric model with reconstructed electrode leads are shown in Fig. 2(a) and (b), here noting model predicted asymmetric shift. With contact 0 (red in Fig. 2(b)) set to −3 V for both implants, the VTAs are obtained by thresholding at −0.5 V. The estimated VTAs are shown in yellow and superimposed onto the model updated MR in Fig. 2(c)–(d) for case 1 and case 2, respectively.

3.3 Neuronal Pathway Recruitment with and without Shift Consideration

The mean computed shift experienced by VTAs (surgical targets) across 2 cases was 1.13 mm. VTA shift was 0.20 and 1.19 mm for case 1 for left and right implants respectively, and 2.60 and 0.52 mm for case 2 for left and right implants. These magnitudes are consistent with previous reports in the literature [1]. Quantitative difference of neuronal pathways resulted from VTAs with and without shift considerations is summarized in Table 1.

Table 1. Comparison of number of recruited tracts, tract length and tract volume for VTA distributions with and without shift considerations (left implant/right implant)

Figure 3(e) provides a 3D view of the recruited pathways for both cases without shift consideration. For qualitative visual comparison, recruited pathways of the target with significant shift (>1 mm) in each case based on VTAs are illustrated in Fig. 3, where results with shift consideration are shown in the top panels and without shift consideration in the bottom panels (axial, sagittal, and coronal views are shown in Fig. 3(a), (b), and (c), respectively). Figure 3(d) further provides a zoomed coronal view highlighting the difference in the extent of recruited pathways with and without shift consideration at target with significant shift (1.19 mm for case 1 and 2.60 mm for case 2).

Fig. 3.
figure 3

Compare neuronal pathway recruitment for VTAs with shift consideration (top) vs. without (bottom) in axial (a), sagittal (b) and coronal view (c). Notable difference in extent in (d). (e) 3D view of neuronal pathways with bilateral implants without shift consideration

4 Discussion and Conclusions

In this work, we presented an integrated framework of patient specific biomechanical and bioelectric models for DBS burr hole procedures. A model-based brain shift correction approach was developed that can provide updated patient MR image volumes and demonstrated good agreement in comparison to iMR data. The correction methodology presents minimal disruption to existing workflow and enables the potential intraoperative deployment of this system in aiding surgical navigation, targeting and direct visualization. Furthermore, by coupling this biomechanical model to a patient specific bioelectric model, VTAs could be estimated and neuronal pathway recruitment could be predicted, which has implications regarding functional outcome and therapy quality. Interestingly, differences in neuronal pathway recruitment are readily observable with and without shift consideration, both qualitatively (different extents in Fig. 3) and quantitatively (e.g. tract number in Table 1). This potentially highlights the need for intraoperative shift correction with the complement of a tool capable of providing bioelectric information to assist surgical targeting and functional therapy. Regarding limitations of the work, the methods must be implemented on a larger population with likely efforts using MER for validation of neuronal recruitment impact. For the biomechanical model, better understanding of the interplay of pneumocephalus and ventricular effects are needed. For the bioelectric model, increased sophistication, e.g. accounting for the frequency-dependency of the stimulation, and enhanced VTA modeling, is desired. Nevertheless, the work herein is provocative in that even small deformations are shown to induce considerable functional change in neuronal pathway recruitment.