1 Introduction

Deep brain stimulation (DBS) is a well-established treatment for Parkinson’s disease, dystonia and essential tremor [1]. To investigate the therapeutic effects and undesirable side effects, it is necessary to estimate the neural response of brain tissue during DBS, known as volume of tissue activated (VTA) [2]. Currently, DBS devices operate in an open loop, where the stimulation parameters are configured in the surgical planning. For this reason, patients are stimulated with fixed values of amplitude, frequency and other stimulation parameters previously set by neurosurgeon [3]. Therefore, the selection of DBS-parameters performed by the neurosurgeon is an empirical and complex process that requires great skill and clinical experience, and for some cases the therapeutic results are not optimal [4, 5].

Nowadays, the methodology for programming the DBS device is based on changing the amplitude and frequency until some symptoms disappear. However, this procedure may cause side effects and a rapid deterioration of device batteries [6, 7]. To avoid side effects, it is necessary to make a good adjustment of DBS-parameters and determine the VTA previously defined by the neurosurgeon [8]. Some procedures have been proposed in literature to deal with this issue. For example, in [9], a methodology is proposed to clarify the mechanism of DBS by using computer simulations and finite element method (FEM). This first attempt allowed to relate the VTA with therapeutic results in PD patients. Nevertheless, it did not resolve the inverse problem: to find optimal parameters from a desired VTA. Then, the authors of [10] developed an interesting approach to find optimal parameters of DBS through detailed models of electrode-tissue interface, achieving interesting clinical outcomes and proving that optimal parameters can be found by analyzing the VTA. On the other hand, due to technological advances, it is possible to collate information on electrode positions with clinical effects, providing a visual representation of the electrical field related to the stimulated nucleus. It allows to incorporate patient-specific computer models to customize DBS parameter settings for each patient. Any examples include frameworks closed-source such as: the Boston Scientific Guide DBS system, Medtronic Optimise system and Renishaw NeuroInspire system [5, 11], and frameworks open-source such as: DBSproc [12]. However, these systems do not identify with acceptable accuracy the parameters of neuromodulation in the treatment of Parkinson’s disease. For this reason, the VTA analysis and optimal parameters of DBS are research topics with considerable gaps, because there are not significant studies regarding the inverse problem of the DBS.

The proposed methodology is performed from initial settings of the DBS, coupled with medical images obtained in the surgical procedure and brain atlas established in a previous work within the Research Group in Automatics of the Universidad Tecnológica Pereira, called NEURONAV [1]. The main objective is to estimate configuration parameters in the DBS device during postoperative therapy of PD patients. To do this, we developed a graphical model with deformable registration of medical images, where the neurosurgeon can set the desired VTA to achieve good therapeutic results and minimal side effects in patients with PD. Then, the proposed system based on a support vector machine (SVM) defines a configuration of the DBS device that adjusts the desired electrical propagation. This proposal focuses on the solution of the inverse problem, which consists of estimating and setting the stimulation parameters for DBS electrodes (voltage/current, electrode impedance, frequency, pulse widths, and contacts), given a previous VTA set by the neurosurgeon.

2 Materials and Methods

2.1 Toolboxes

To build the graphic module, we use open source tools written in c++. To achieve an interactive tool where the neurosurgeon can identify, visualize and manipulate the VTA, we employed the visualization tool kit VTKFootnote 1. To build a graphical user interface compatible with VTK and 3D graphics support via OpenGL, we used FLTKFootnote 2 toolkit that incorporates a GLUT emulation. Finally, to simulate the geometric variations of the VTA, we use the GetFEMFootnote 3 tool, which offers a framework for solving potentially coupled systems of linear and nonlinear partial differential equations with the finite element method (FEM). GetFEM is interconnected with languages such as python and allows to solve problems in 1D, 2D and 3D.

2.2 Surface Reconstruction

The 3D representation of VTA is done with tetrahedral mesh elements using the Delaunay method. This technique connects polygons to plot the surface of the VTA previously estimated with a GetFEM. In the bipolar configuration, it is necessary to take into account the electrode dimensions to decide which are the points corresponding to the active contacts that interact with the VTA surface. Figure 1 illustrates the pipeline of surface reconstruction inside the graphical module. The data objects (source) represent and provide access to data, and the process objects (filter) operate over data objects. Here, the connections are made using the methods setInput()/getOutput(). Once the model is built, the execution of the pipeline must be carefully controlled. In this regard, VTK uses a process of implicit distributed update.

Fig. 1.
figure 1

Pipeline for the reconstruction of a VTA previously estimated with GetFEM.

Full size image

2.3 Geometric Variations

The typical procedure to represent a point in a 3D environment is through cartesian vectors of three elements, \(\mathbf {p}=(x,y,z)\). The geometric transformations are represented in homogeneous coordinates by a \(4\times 4\) transformation matrix denominated \(\mathbf {T}\), such that homogeneous coordinates of a point \(\mathbf {P}=(x,y,z,1)\), becomes geometrically at point \(\mathbf {P'}=(x',y',z',1)\) with the following operation \(\mathbf {P'}= \mathbf {T}\times \mathbf {P}\), known as affine transformation:

To find the area of stimulation expected within the STN region, we apply the affine transformation to the actor object for simulating the neural activation desired by the neurosurgeon. Accordingly, we perform a transformation operation (vtkTransform) over the object (vtkActor). To apply a scaling matrix on the actor object, we center the object in the coordinates (0, 0, 0) of the render window.

2.4 Parameters Estimation with Machine Learning

We train a set of machine learning algorithms using 1000 isotropic VTA models validated by the medical team of The institute of epilepsy and Parkinson of the eje cafetero-Neurocentro, Colombia. Specifically, we train a SVM-regressor for the amplitude and pulse width. For active contacts \((C_0-C_3)\) we employ SVM classifiers. In this context, the labels \(\mathbf {y}\) are the parameter values. While, the numeric flag of VTA (0 refers to non-VTA and 1 is VTA positive) concatenated with the spatial coordinates correspond to the training features \(\mathbf {X}\).

2.5 Experimental Setup

For the graphical module of VTA analysis, we use the Universidad Tecnologica de Pereira (DB-UTP) database. It contains recordings of MRI studies (T1 and T2 sequences with \(1\times 1 \times 1\) mm\(^3\) voxel size) from four patients with Parkinson’s disease. The STN region in MRI studies was labeled for specialists of Neurocentro. We use the DBS electrodes 3387 and 3389 manufactured by MedtronicFootnote 4 (Medtronic, Inc. USA). Both electrodes consist of four contacts with a length of 1.5 mm separated by 1.5 mm (3387) and 0.5 mm (3389). The diameter is 1.27 mm and the contact area has 6 mm\(^2\). Each contact can be used as anode or cathode in bipolar electrode configuration or as cathode in monopolar stimulation setting [13]. The electric field was simulated in the neighboring of STN, with an electric potential in a range of 0.5 V to 10.5 V with steps of 0.5 V, and pulse width of 60 \(\upmu \)s to 450 \(\upmu \)s (steps of 30 \(\upmu \)s). Active contacts are \(C_0, C_1, C_2, C_3\). GetFEM is used to solve specific VTA models for each patient.

Figure 2 shows the flow diagram to run experiments with isotropic models. First, it is performed the data acquisition: medical imaging, VTA points, and the database for training the learning algorithms (to identify the stimulation parameters for a new geometric variation of the VTA). Second, we observe the processes performed to modify interactively the VTA: reconstruction of the surface and affine transformation used for the deformation of the VTA. Finally, we have the output data (desired VTA) and the corresponding parameter values given by the SVMs.

Fig. 2.
figure 2

Flow diagram for identifying DBS parameters given a desired VTA (inverse problem solution).

Full size image

3 Experimental Results and Discussion

3.1 Parameters Estimation

In Table 1, we show accuracy results for the DBS parameters estimated with the SVM. The dataset for training the SVM has 1000 VTA models. We employ a hold-out validation scheme, where \(70\%\) of data is used for training, and \(30\%\) is kept for validation. We repeat the experiment 100 times, with a random selection of training and validation samples. These results show high accuracy to identify the parameters for a desired VTA. Therefore, the inverse solution problem can be solved with machine learning techniques.

Table 1. Table of results of the accuracy in estimating the parameters.
Full size table

In Fig. 3, we observe a simulated monopolar VTA given a geometric variation desired by a neurosurgeon. This procedure allows to control and visualize the neural activation, reducing side effects and minimizing the battery consumption of the implanted pulse generator. The electrode is located within the subthalamic region and the first contact is active and it stimulates a large section of the STN.

Fig. 3.
figure 3

VTA monopolar within the STN for a desired configuration. The blue cylinder is the stimulating electrode (Medtronic-3387) with four contacts (green areas), the VTA corresponds to the red surface and the STN is the cyan colored structure. (Color figure online)

Full size image

3.2 Surface Reconstruction and Geometric Variations

Best computational performance of the module with an affine transformation to the vtkActor in the render window. We can see in Fig. 4 the points cloud which returns the GetFEM software, which constitutes a discrete estimation of VTA. For a simulation neural response generated by DBS in the brain, the Laplace equation describes the mathematical model. For this reason, a partial differential equation must be solved with finite element method. Once, the VTA points (spatially distributed) are estimated with GetFEM, we transform them in a surface depending of the type of stimulation (monopolar or bipolar).

Fig. 4.
figure 4

Monopolar and bipolar VTA, (A) Point cloud, which returns GetFEM. (B) Reconstruction delaunay as from the point cloud. (C) Geometric variations on applying affine transformation vtkActor.

Full size image

4 Conclusions and Future Work

In this work, we presented an interactive module for DBS parameters identification through supervised machine learning. The estimation of these parameters depends on a geometric variation of isotropic VTA models given by a neurosurgeon. This is known as the inverse problem of DBS. This module is integrated to the previously developed software: NEURONAV [1].

Results demonstrated that inverse problem can be successfully solved using machine learning methods. Also, the geometric variations allows to estimate monopolar and bipolar configurations of parameters. In addition, The module allows to visualize the variations of VTA without considerable delays. This property could be used in DBS therapy for controlling therapeutic outcomes and side effects in PD patients.

Three main tasks are established as future works: First, the proposed scheme can be used with anisotropic models studies using diffusion tensor imaging (DTI) of each patient. Second, to apply volume reduction techniques for generating representative models instead of a manual manipulation of the VTA. Third, to use codifications of GPU through the API of CUDA for improving processing times.