SciELO - Scientific Electronic Library Online

 
 issue95Efficiency of the removal of microcystin-LR by UV-radiation and hydrogen peroxideDesign of a load carriage system oriented to reduce acceleration forces when carrying a backpack author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

  • On index processCited by Google
  • Have no similar articlesSimilars in SciELO
  • On index processSimilars in Google

Share


Revista Facultad de Ingeniería Universidad de Antioquia

Print version ISSN 0120-6230On-line version ISSN 2422-2844

Rev.fac.ing.univ. Antioquia  no.95 Medellín Apr./June 2020

https://doi.org/10.17533/udea.redin.20190729 

Artículo original

A machine learning approach to support deep brain stimulation programming

Un acercamiento basado en aprendizaje de máquina para apoyar la programación de la estimulación cerebral profunda

Viviana Gómez-Orozco1  * 

Iván de la Pava1 

Andrés Álvarez-Meza2 

Mauricio A. Álvarez3 

Álvaro Orozco-Gutiérrez1 

1Grupo de Investigación Automática, Facultad de Ingeniería, Universidad Tecnológica de Pereira. Carrera 27 # 10-02. C.P. 660003. Barrio Álamos, Pereira, Risaralda - Colombia.

2Grupo de Procesamiento y Reconocimiento de Señales, Universidad Nacional de Colombia. Carrera 27 # 64-60 Sede Manizales. C. P. 170004. Manizales, Caldas - Colombia.

3 Department of Computer Science, University of Sheffield. Western Bank Sheffield S10 2TN UK. C. P. S10 2TN UK. Sheffield, United Kingdom.


ABSTRACT

Adjusting the stimulation parameters is a challenge in deep brain stimulation (DBS) therapy due to the vast number of different configurations available. As a result, systems based on the visualization of the volume of tissue activated (VTA) produced by a particular stimulation setting have been developed. However, the medical specialist still has to search, by trial and error, for a DBS set-up that generates the desired VTA. Therefore, our goal is developing a DBS parameter tuning strategy for current clinical devices that allows defining a target VTA under biophysically viable constraints. We propose a machine learning approach that allows estimating the DBS parameter values for a given VTA, which comprises two main stages: i) A K-nearest neighbors-based deformation to define a target VTA preserving biophysically viable constraints. ii) A parameter estimation stage that consists of a data projection using metric learning to highlight relevant VTA properties, and a regression/ Classification algorithm to estimate the DBS parameters that generate the target VTA. Our methodology allows setting a biophysically compliant target VTA and accurately predicts the required configuration of stimulation parameters. Also, the performance of our approach is stable for both isotropic and anisotropic tissue conductivities. Furthermore, the computational time of the trained system is acceptable for real-world implementations.

Keywords: Volume of tissue activated; kernel-based learning; anisotropy

RESUMEN

Ajustar los parámetros de estimulación es un desafío en la estimulación cerebral profunda (DBS), debido a la gran cantidad de configuraciones disponibles. Como resultado, se han desarrollado sistemas basados en la visualización del volumen de tejido activado (VTA) producido por una configuración de estimulación particular. Sin embargo, el especialista todavía tiene que buscar, mediante ensayo y error, una configuración DBS que genere el VTA deseado. Por lo tanto, nuestro objetivo es desarrollar una estrategia de ajuste de los parámetros de DBS para los dispositivos clínicos actuales que permita definir un VTA objetivo bajo restricciones biofísicamente viables. Proponemos un enfoque de aprendizaje de máquina que permite estimar los valores de los parámetros de DBS para un VTA dado, que consta de dos etapas principales: i) Una deformación basada en K-vecinos más cercanos para definir un VTA objetivo sujeto a restricciones biofísicas. ii) Una etapa de estimación de parámetros que consiste en una proyección de datos para resaltar las propiedades relevantes del VTA, y un algoritmo de regresión/Clasificación para estimar los parámetros DBS necesarios para generar el VTA objetivo. Nuestra metodología permite establecer un VTA objetivo compatible biofísicamente y predice con precisión la configuración requerida de los parámetros de estimulación. Además, el rendimiento de nuestro enfoque es estable tanto para conductividades del tejido isotrópicas como anisotrópicas. Además, el tiempo de cómputo del sistema entrenado es aceptable para implementaciones en el mundo real.

Palabras Clave: Volumen de tejido activado; aprendizaje basado en kernels; anisotropía

1. Introduction

Deep brain stimulation (DBS) is a neurological therapy that consists of sending electrical pulses of particular amplitude, pulse width, and frequency to a target region of the brain through a neurostimulator [1, 2]. DBS is conducted mainly for the treatment of movement disorders, such as Parkinson´s disease, dystonia, and essential tremor [3, 4]. Nowadays, its field of action has expanded to the treatment of various neurological disorders [5-8], despite persistent questions about its mechanisms of action [9, 10]. From a clinical perspective, one of the significant post-surgical challenges for therapy success is the DBS device programming. This process consists of tuning the stimulation parameters (amplitude, pulse width, frequency and electrode contact state) in a way that maximizes the therapeutic benefits while minimizing adverse side effects [11]. Programming the device is inherently complex due to the high number of possible stimulation configurations, without mentioning other patient-dependent variables, such as electrode placement and local brain conductivity, which makes it a long and tedious process that depends heavily on the medical specialist´s experience [12-15].

In the last decade, several applications have been developed to assist in the programming of DBS devices based on the visualization of the spatial extent of direct neuronal activation in response to external electrical stimulation (Volume of Tissue Activated (VTA)) [16-20]. The VTA visualization with accurate 3D reconstructions of the stimulation electrode and its surrounding brain structures allows relating the applied electric pulse train to the therapeutic effects. So, areas of the structures of interest that are affected by the DBS can be extracted based on the VTA [20-22]. Accordingly, systems for VTA visualization allow the medical specialist to vary the stimulator parameters while predicting the VTA, and thereby to gain an idea of the required configuration to achieve the desired therapeutic outcome [15, 23].

The computation of the VTA is grounded on electrophysiological models of neural activity that estimate the brain tissue response to the distribution of electric potential generated by the stimulator [16, 24, 25]. The complexity of such models varies widely depending on the level of detail in both the electrical and the neural characterization of the electrode, and the brain tissue, i.e., homogeneous, heterogeneous, (isotropic or anisotropic brain tissue conductivities); neuronal or axonal models with idealized or tractography-based morphologies and trajectories, among others [25-29]. However, even the simplest model carries a high computational time, so to allow their integration into practical applications, some alternative simulation strategies have been developed [30-32].

In [33], the authors proposed the use of artificial neural networks to estimate the VTA given a combination of stimulation parameters. After training the system with DBS parameters as inputs and elliptic profiles representing the VTA as outputs, it achieved a lower computational time than that of standard models [16]. Likewise, our group developed a hierarchical K-nearest neighbor-based approach that accomplishes the same goal without simplifying the shape of the VTA [34]. Nonetheless, in spite of the considerable benefits offered by the previously described VTA computation and visualization systems, they still depend on a heuristic search (trial and error) of the stimulation parameters necessary to generate a target (desired) VTA. This limitation promises to become increasingly evident because the designs for new generations of DBS devices offer an exponentially higher number of possible combinations of parameters [28, 35-37].

In light of the above, the idea of developing an application that allows estimating the stimulation parameters from a target VTA arises. In such a system, the medical specialist should only define the shape of neuronal activation that he/she deems necessary to improve the patient’s symptomatology. Then, the system would return the set of stimulation parameters that generates the VTA closest to the desired one, following the restrictions imposed by the assumed electrophysiological model.

Unlike the problem of estimating the VTA from the stimulation parameters, the problem in the reverse direction has been less studied [38]. In [39] the authors proposed a particle swarm optimization approach to program DBS devices, assuming each set of stimulation parameters as a particle. Their approach aims to maximize axonal activation in a region of interest, to minimize axonal activation in an area of avoidance, and to minimize power consumption. Their algorithms are based on an optimization function that relies on a smooth activation to compute the axonal response, this is less costly than full electrophysiological models of neural activity and produces accurate predictions for monopolar stimulation. However, activation function predictions are not accurate when dealing with bipolar stimulation [33]. This limitation is not relevant when working with new experimental DBS devices, such as the one used by Peña et al., because their complexity allows obtaining any desired VTA shape using cathodic stimulation only. But it is relevant when working with current DBS clinical devices, such as the Medtronic’s ACTIVA series, because bipolar stimulation is sometimes a necessary strategy to achieve a positive therapeutic outcome [13], given the limitations of the stimulator.

In this work, we propose a machine learning strategy that allows obtaining a combination of DBS parameters for current clinical devices that generates a target VTA. Our approach considers physical constraints of the stimulation system and the brain tissue to find an appropriate DBS tuning. Broadly, the introduced strategy comprises the following stages: i) A K-nearest neighbors-based deformation to define a target VTA preserving biophysically viable constraints. In this regard, the VTA is represented as elliptic profiles, following the approach introduced in [33]. ii) A parameter estimation stage that consists of a data projection based on Centered Kernel Alignment (CKA) [40, 41], to highlight relevant VTA properties regarding DBS parameter variations, and a regression/classification algorithm using kernel-based machines to estimate the required stimulation parameters that generate the target VTA. With the aim of testing the proposal’s flexibility, computational simulations are conducted under different scenarios: low, medium, and high brain tissue impedances values, as well as, isotropic and anisotropic brain tissue conductivity conditions [21, 25]. For each simulation scenario provided, a 10-fold cross-validation strategy is employed. More precisely, a quantitative error value is computed between the VTA simulated using the gold standard algorithm [16], and the VTA belonging to the DBS configuration fixed by our machine learning approach. Our methodology successfully estimated DBS parameter combinations for the target VTAs under the different conditions tested. Concerning the computational time, our strategy requires less than 1.5 seconds in a computer of average specifications (see section 3.3) to find a viable DBS set-up from a target VTA, which would allow its inclusion in practical applications.

2. Materials and methods

With the aim to support DBS programming tasks, we introduced a machine learning-based approach that estimates the required stimulation parameter values from a given set of ellipses representing a target VTA and a fixed stimulation pulsed width, defined in function of the therapeutic requirements identified by the medical specialist. In particular, our approach estimates the stimulation amplitude and the activation state and polarity of the electrode contacts. In what follows, the proposal’s main stages (see Figure 1] are described in detailed.

Figure 1 The main sketch of the introduced machine learning approach to support DBS programming. (a) Inputs. The desired pulse width value [µs] and an ellipse-based representation of the target VTA. (b) Deformation and parameter estimation. A KNN-based deformation is applied to obtain a biophysically compliant and binary coded VTA. Next, a CKA-based projection is carried out to reveal relevant features regarding DBS parameter variations. Then, a SVM regressor and a SVM classifier are trained to compute the DBS parameters. (c) Outputs. Amplitude value A[V ] and electrode contact states c 0, c 1, c 2, c 3 ϵ {0, -1, 1}: The achieved DBS set-up is evaluated by estimating its corresponding VTA  

2.1 A binary coding for the VTA

The gold standard for VTA estimation involves, as a first stage, computing the electric potential generated by the stimulator for a specific set of stimulation parameters [16]. Each set of Q stimulation parameters relevant for VTA estimation can be coded as a scalar , the pulse width in µs, and a vector y = [A, c0, c1, c2, c3], where A is the stimulation amplitude in V, and cr, with r = 0, 1, 2, 3, stands for the activation state and polarity of the r-th electrode contact (-1: cathode, 0: inactive, 1: anode). The electric potential distribution is then used to excite multicompartment axon models arranged in a field in the vicinity of the DBS electrode shaft, as shown in Figure 2(a). Axons that fire an action potential as a response to each stimulation pulse are considered active and their spatial distribution defines the VTA [16]. Thus, a VTA computed from a field of P axons can be coded as a binary vector x ϵ {0; 1} P through axon response concatenation, where xp = 1 if the p-th axon is activated by the stimulation, and xp = 0 otherwise (p ϵ {1, …, P}). Following this representation a dataset of N VTAs and their corresponding combinations of stimulation parameters can be stored as matrices X ϵ NxP and Y ϵ NxQ-1, and a vector ϵ N, where each row vector yi ϵ Q-1 (i ϵ {1, …, N}) in Y, and each element in , hold the stimulation parameters employed to compute the i-th VTA xi ϵ {0, 1}P.

Figure 2 VTA representation. (a) Representation of an axonal field built in the vicinity of a DBS electrode (Medtronics electrode 3389). The VTA is defined by the spatial distribution of the axons activated by the delivered stimulus (green dots). (b) Ellipse-based characterization of the VTA. The input to the deformation stage is a VTA represented through a set of characteristics of two ellipses: their centers [o1, o2], and their axes lengths [a1, a2] and [b1, b2]  

2.2 K-nearest neighbors-based VTA deformation

In this section, we introduce a strategy to obtain a target VTA that closely resembles a volume of activation defined by a specialist. It also follows the restrictions imposed by the assumed properties of the brain tissue and by the characteristics of the stimulation device, i.e. it is constrained by the possible shapes and locations that a VTA can have, based on a given biophysical model. The proposed solution consists of a machine learning technique that uses a characterization of the VTA through ellipses and approximation by nearest neighbors to estimate the target or deformed VTA.

Ellipse based characterization of the VTA

In this stage, we use an ellipse based characterization to represent and manipulate the VTA, following the approach first explored in [33]. The VTA is represented by a set of characteristics that correspond to the centers ([o1; o2]), and the axes lengths ([a1; a2] and [b1; b2]) of two ellipses, as shown in Figure 2(b). The centers and axes lengths of the ellipses are measured with respect to a reference frame with origin in the electrode tip. When the VTA can be characterised by a single ellipse, e.g. only one contact is active, then o1 = o2, a1 = a2, and b1 = b2. So, the input to the deformation algorithm is a vector z = [o1; o2; a1; a2; b1; b2] that contains the geometrical parameters of two ellipses. The use of such simple geometric forms is intended to allow an intuitive manipulation of the VTA by a medical specialist familiar with common VTA graphical representations [15, 17].

VTA deformation method

Given an input vector z* that contains the geometrical parameters of a pair of ellipses representing a target VTA and a dataset of N VTAs X ϵ NxP; a new VTA x* ϵ {0, 1]P can be approximated through a weighted average (Equation (1)], as follows:

where Ωx is a set containing the K nearest neighbors of z* in Z = f(X) and λk ϵ + (k ϵ {1, …, K}). The function f (∙) characterises every binary VTA vector x in X in terms of its associated ellipse based representation z. The neighborhood Ωx (Equation (2)] is defined as:

where zK is the K-th neighbor of z* in Z according to a Gaussian kernel kZ(∙,∙) (Equation (3)]:

where notation de(∙,∙) stands for the Euclidean distance and σZ ϵ +. The kernel bandwidth was fixed as σZ = med(de(zi, zj)), where med(∙) stands for the median operator. The Gaussian kernel was selected because of its representation capacity.

To asses the relative importance of zk, a softmax gating function is used to estimate λk (Equation (4)]:

Where Ωz is the set containing the corresponding elements of Ωx in the ellipse-based representation space. Finally, a thresholding procedure (µ ϵ [0, 1]) is applied to estimate the new binary VTA x* (Equation (5)] as:

Figure 3 outlines the main steps required to implement the above described method. This approach allows estimating a target VTA whose characteristics match those of the VTAs in a precomputed dataset, both in terms of their mathematical representation and biophysical behavior. The deformed VTA and the precomputed dataset can then serve as inputs to the algorithms for parameter estimation.

Figure 3 Outline of the proposed VTA deformation method. Given a pair of ellipses representing a desired VTA z*, a VTA dataset X is mapped to the ellipse based representation space through the mapping function f (∙) to obtain the feature space Z. Next, the Gaussian similarity between z* and every element of Z is computed to determine which ellipses in Z are the closest to the new pair of ellipses. Finally, a new binary VTA x* is obtained through a K-nearest neighbors approximation, that is, a weighted average of the VTAs in X corresponding to the K ellipses in Z closest to z*  

2.3 Parameter estimation

We use an approach based on kernel methods to estimate the configuration of stimulation parameters that generates the VTA closest to the target volume. Figure 4 outlines our approach. It employs Centered Kernel Alignment (CKA) to map the target VTA to a lower dimensional space where the VTA characteristics that better allow predicting the parameters that could have generated the target volume are highlighted. Then, it employs regression and classification Support Vector Machines (SVMs) algorithms to find the stimulation parameters, following the restrictions imposed by the DBS system.

Figure 4 Main scheme of the DBS parameter estimation approach based on kernel methods. Given a low-dimensional VTA obtained by CKA, the ellipse based representation of said VTA, and a pulse width value ꞷ, the stimulation parameters for the target volume are computed through five SVMs: one regression for the amplitude, and four independent classifiers for the state of each electrode contact  

CKA-based data projection

CKA allows quantifying the similarity between two sample spaces, such as the VTAs and the DBS parameters, by comparing two characterizing kernel functions [40]. Provided a VTA dataset X ϵ NxP , and the corresponding DBS parameters Y ϵ NxQ-1 and ϵ N, we define two Gaussian kernel matrices KX ϵ NxN and KY′ ϵ NxN. The first matrix holds elements = k X(x i , x j ) with x i , x j ϵ X, while the second matrix has elements = k Y (y I , y j ) with y i , y j ϵ Y , Y = {Y ,} ϵ NxQ. The CKA alignment can be estimated as (Equation (6)]:

where stands for the centered kernel matrix K obtained as = , where = I - 1 1/ N is the empirical centering matrix, I ϵ NxN is the identity matrix, 1 ϵ N is an all-ones vector, and = denotes the matrix-based Frobenius norm. For the Gaussian kernel k Y′ the pairwise comparison between samples y n and y m is performed using the Euclidean distance operator d e (∙,∙), while for k X the Mahalanobis distance (Equation (7)] is selected. Namely, the distance function of k X is fixed as:

where Γ ϵ PxM, with M < P, is a linear projection matrix, and ΓΓ⊤ is the corresponding inverse covariance matrix. Then, the projection matrix Γ is computed by solving the following optimization problem (Equation (8)]:

where the logarithm function is used for mathematical convenience. After estimating we use it to linearly map x* ϵ {0, 1} P to a lower-dimensional space (Equation (9)]:

This mapping facilitates the estimation of the DBS parameters and decreases computational time of the estimation process.

Parameter estimation algorithms

The mapped VTA ϵ M is then concatenated with the information provided by the ellipse-based characterization of the VTA z* and the pulse width value determined by the specialist to form the input vector v* = [z*, , ] ϵ L, with L < P. The value is added to the input vector to avoid the ill-conditioning of the parameter estimation problem, that is, the fact that different sets of stimulation parameters can lead to the same VTA. Such constraints on have also been included in other parameter estimation approaches [39]. Next, the proposed framework uses two kernel machines to estimate the appropriate amplitude and electrode contacts configuration (maximum two active contacts at a time) for v*. A SVM regressor is used to model the amplitude; while the state of each electrode contact is modeled by an SVM classifier, through a one versus the rest scheme [42].

Given X, Y and , we define V = [Z, XΓ, ] ϵ NxL, with XΓ = XΓ ϵ NxM. From Y we extract a vector ψR ϵ N and a matrix ΨC ϵ {-1, 0, 1}Nx4 holding the amplitude values and the contacts configuration, respectively. Next, we build two functions: the first function fR: M( estimates the amplitude parameter ϵ in ψ R as: = f R (v i ) = φ R (v i )w+b, where w ϵ MR, b ϵ , and φR: M(MR. Then, a support vector regression optimization problem (Equation (10)] can be defined as follows:

where γR ϵ +is a regularization parameter, ui = , and = (u-ϵ)2 if u i ϵ, otherwise, = 0 (ϵ ϵ +).

Regarding the second function, fC: M( {-1, 0, 1}, which allows computing the DBS contact configuration vector ϵ {-1, 0, 1}4 in Ψ C , we built a soft margin support vector classifier over V to compute the r-th contact value as (Equation (11)]:

where ϵ is the weight of training sample j for the r-th classifier and a r ϵ is a bias term [42].

2.4 Experiments and validation scheme

To evaluate the performance of the proposed methodology, we built six VTAs datasets, each one generated by 1, 000 pseudo-randomly selected combinations of commonly used stimulation parameters for the Medtronic ACTIVA-RC stimulator [13, 43]. We limited the DBS stimulation parameter space as follows: a maximum of two electrode contacts were considered simultaneously active, i.e. monopolar (one active contact) and bipolar (two active contacts) configurations. The electric potentials of each dataset were generated for 200 monopolar and 800 bipolar stimulation configurations, with a set frequency of 130 Hz. Finally, the amplitude was restricted to the range A ϵ [0.5, 5.5] V in steps of 0.5 V, and the pulse-width to ϵ [60, 450] µs in steps of 30 µs.

To simulate the electric potential generated by each combination of stimulation parameters in the datasets, we modeled the brain tissue as a cube of length 10 cm, and positioned a 3D model of a clinical electrode (Medtronic DBS 3389 electrode) [27] in the middle of it. For half of the datasets we assumed an isotropic bulk tissue conductivity of 0.2 Sm-1. For the other half we used diffusion tensor-based conductivities to model anisotropic conditions in the basal ganglia region [21, 44, 45]. We also included a representation of a 0.5 mm encapsulation layer around the electrode and set its conductivity to three different values (0.680 Sm -1 , 0.128 Sm -1 , 0.066 Sm -1 ) to represent low (~500 Ω), medium (~ 1. 000 Ω), and high (~ 1, 500 Ω) impedance conditions [25, 33]. Each of the three datasets built per bulk brain tissue conductivity condition corresponds to a different conductivity value for the encapsulation layer. The model also integrated a voltage drop of 48% at the electrode-tissue interface. We used the same boundary conditions as in [19] to model the active contacts of the electrode and the ground. Afterwards, we used the finite element method (FEM ) [46] software COMSOL Multiphysics 4.3 to obtain the electric potential distribution in the brain tissue model for a particular combination of stimulation parameters by solving the Poisson’s equation.

Then, to simulate the neural response to the stimulation for the isotropic case, the electric potential was interpolated onto each section of a field of 4, 144 multi-compartment axon models of diameter 5.7 µm, which were oriented perpendicularly to the axis of the electrode, and positioned into a grid of width 9 mm, height 27.75 mm, and inter-axonal space 0.25 mm in both the vertical and the horizontal directions [47]. For the anisotropic case, the field consisted of 8, 112 axon models oriented in four different directions perpendicular to electrode axis, with an inter-axonal spacing of 0.5 mm. The response of the axons to the stimulating potentials defined each VTA: axons that fired an action potential per stimulation pulse were considered active and their positions in space shaped the VTA [16]. All the axonal response simulations were implemented in NEURON 7.3 configured as a Python module [48].

Finally, to asses the performance of the proposed methodology we characterized each VTA in the datasets in terms of their ellipse based representation. These ellipses and their associated pulse widths were used as inputs to the machine learning algorithms. The CKA optimization problem (equation (8)] was solved through iterative gradient descent optimization. The support vector regression optimization problem (equation (10)] and each classifier (equation (11)] were solved by quadratic optimization from the well-known SVM dual formulation [42]. The parameters of these algorithms were selected through nested cross-validation. Then, the returned DBS parameters were used to estimate new VTAs, through the methodology described in [34]. Lastly, the new VTAs were directly compared with the corresponding reference VTAs in the datasets through the error metric (Equation (12)] introduced in [33]:

where FP and FN are the false positives and false negatives with respect to the reference VTA, and AA is the reference number of active axons. We chose this evaluation approach because the final aim of our methodology is to estimate a set of parameters that would generate a VTA as close as possible to the target volume. Also, we chose a simplified VTA estimation method for the last part of the validation scheme, instead of the gold standard, because any practical application would require a fast computation of the VTA generated by the estimated parameters (that may differ slightly from the defined target) and the gold standard is too computationally intensive for that purpose. It is important to note that any other fast VTA computation method of similar performance, such as the one in [33], could have been used as part of the validation scheme. Additionally, we carried out tests to evaluate the performance of the machine learning algorithms in the different stages of the proposed methodology, e.g. the precision in the estimation of the stimulation parameters by the SVMs. All experiments were implemented in Matlab R2015a following a nested cross-validation scheme with 10 repetitions, where 70%of the samples were used as a training set and the remaining 30 % as a test set. According to the nested scheme, the free parameters of the algorithms were tuned using 80 % - 20 % 5-fold cross-validation over the training set for each of the main repetitions. Finally, in order to study if there are statistically significant differences in the performance of the proposed method regarding the impedance and conductivity conditions tested, we used the Kruskal-Wallis test over the error distributions obtained from the cross-validation scheme. If the null hypothesis for equal medians is rejected, we performed multiple comparison tests using the Tukey-Kramer scheme [49]. All significance levels were measured at 5%.

3. Results

3.1 VTA deformation

The first stage of the proposed methodology resides in the VTA deformation from a set of ellipses and their associated parameters. In turn, the free parameters of the algorithm, the number of neighbors K and the threshold µ, were tuned through a 10-fold grid search in the ranges K = [1, 99] and µ = [0, 1] aiming for the combination of parameters that generated the smallest error. That is, the smallest error between the deformed VTAs, obtained from the ellipse based characterization of the datasets, and the reference VTAs. Figure 5 shows the average error surface obtained in function of K and µ for the VTA deformation algorithm in the isotropic [Figure 5(a)) and anisotropic [Figure 5(b)) conditions, for an impedance of 1, 000 Ω. For both conductivity conditions it was found that the combination of parameters that generated the smallest error were K = 3 y µ = 0.4. The behavior of these error surfaces with respect to the number of neighbors K indicates smooth variations among VTAs, since small neighborhoods yield the smallest errors (the same behavior was observed for the remaining impedance values).

Figure 5 Average error in the deformation of the VTA from a set of ellipses. (a) Isotropic conditions and (b) anisotropic conditions, for an impedance of 1000 Ω. The error variance remained relatively constant, at around 0.16 for the isotropic case and 0.25 for the anisotropic case, for the values of K and µ analysed  

Once the parameters of the algorithm were tuned, its performance was evaluated for each impedance value and conductivity condition. Figure 6 shows the errors obtained. Overall, the performance of the deformation algorithm is not significantly affected by the impedance value, except for low impedance in the anisotropic condition. Namely, for the anisotropic condition the pairwise statistical comparison of the obtained error distributions for low-medium, low-high, and medium-high impedances yields p-values of 0.042, 4e-4, and 0.166, respectively. While for the isotropic condition an analogous analysis yields p-values of 0.695, 0.329, and 0.065. On the other hand, a notorious difference between the error levels in the isotropic and anisotropic conditions is observed, with the latter doubling those of the former. However, even for that condition the error remains low. These low error levels indicate that the proposed deformation method allows for a satisfactory reconstruction based on a binary representation from the parameters of a set of ellipses [33], following the restrictions imposed by the tissue and by the stimulation system. Also, the method achieves the latter by estimating a new VTA from a set of precomputed VTAs through a weighted average. Thus, the new VTA inherits the properties of the neighboring VTAs, under the assumption of smooth variations among the theoretically possible shapes that the VTA can hold.

Figure 6 Box-whisker plots of the errors obtained for the ellipse based deformation method, discriminated by impedance value for the isotropic (bottom row) and anisotropic (top row) conditions  

3.2 Parameter estimation

Five independent SVMs were trained, one machine for regression (amplitude) and four machines for classification (contacts state), to estimate the parameters of stimulation that allow generating the target VTA. As evaluation metrics we used the relative error for the regression machine and the classification accuracy for the classification machines, where the relative error was defined as the ratio of the absolute error of the estimated value of A to the true value of A. Figure 7 shows the relative error in the estimation of the stimulation amplitude A, and the estimation accuracy for contacts c 0 to c 3, discriminated by impedance condition and brain tissue conductivity, isotropic (7(a) to 7(c)) and anisotropic (7(d) to 7(f)). The relative error in the estimation of the amplitude for the isotropic condition oscillates between 0.15 and 0.17. For the state of the electrode contacts the estimation accuracy is higher than 90%. For the anisotropic condition, relative errors of about 0.25 for the amplitude and accuracies of about 90% for the state of each of the electrode contacts are obtained. Except for the low impedance value where the lower bound of the accuracy is about 80%. In comparison with the isotropic condition, these relative errors are higher and the accuracies are evidently lower. For both conditions there is also an upward trend in the mean values of the accuracy as the value of the impedance of the encapsulation tissue increases, with an accompanying reduction in the dispersion of the accuracy data for the estimation of the activation state of each contact.

Figure 7 Relative error and estimation accuracies obtained with the SVMs for the stimulation amplitude A, and the state of the electrode contacts c 0 to c 3, respectively, discriminated by impedance value and brain tissue conductivity condition: isotropic (top row) and anisotropic (bottom row) 

Figures 8(a) to 8(c) show the number of support vectors, expressed as a percentage of the total number of elements of the training set, used by the classification machines trained to estimate the state of each of the contacts, for each impedance condition under study in the isotropic condition. Figures 8(d) to 8(f) show the support vectors used in the anisotropic condition. Overall, the small percentage of support vectors points to the fact that the classifiers are not overfitting the data. Besides, a downward trend is observed in the number of vectors used by the machines in the training as the impedance increases, which indicates that discerning the configuration of each contact for low impedance values was more challenging for the classification SVMs, since the percentage of support vectors tends to increase as the underlying function or decision boundary becomes more complex [42]. This behaviour, and in general the lower performance for low impedance conditions, is due to the fact that the volumes under those conditions tend to cover the spatial regions adjacent to the activated contacts, including areas corresponding to non-active contacts as shown in Figure 9(a), making it difficult to identify the activation state and polarity of each of the four electrodes. In contrast, for high impedance values the extent of tissue activated by the electrical stimulation is reduced, which leads to clearly separated volumes for each contact, as shown in Figures 9(b) and 9(c). On the other hand, the largest errors obtained for the anisotropic tissue condition can be attributed to the fact that in this case the activation volume presents relatively complex morphologies in comparison to the symmetric shapes of the volumes computed for the isotropic assumption, as shown in Figure 9(d).

Figure 8 Number of support vectors (as a percentage of the total number of elements in the training set) selected in the classification task for each of the electrode contacts c 0 to c 3, discriminated by impedance value and brain tissue conductivity condition: isotropic (top row) and anisotropic (bottom row)  

Figure 9 VTAs obtained for a single combination of stimulation parameters (A=2 V, =90µs, c 0 = 0, c 1 = -1, c 2 = 0, c 3 = 0) for different impedance values, (a) 500 Ω, (b) 1000 Ω, and (c) 1500 Ω under isotropic conditions; and (d) 1000 Ω under anisotropic conditions  

3.3 Testing the proposed machine learning approach

In order to evaluate the proposed methodology in terms of its ability to generate a set of parameters that result in a VTA close to the target one, the error between the VTAs computed from the estimated parameters and the reference VTAs was obtained through equation (12). The estimation of the VTAs from the parameters was made using a hierarchical K-nearest neighbor (HKNN) approach that consists of a data-driven approximation of an expected squared error functional to obtain a local weighted VTA averaging [34]. Figure 10 shows the final reconstruction error for the isotropic and anisotropic cases, respectively. The error distributions, involving information of both monopolar and bipolar configurations, are in ranges of values similar to those reported in the literature for machine learning based VTA estimation strategies that solve only the traditional VTA estimation problem [33, 34]. Also, in congruence with SVM performance a slight downward trend is observed in the median error as the impedance value increases, and higher errors appear when the conductivity of the tissue is assumed to be anisotropic (with respect to the isotropic condition) due to the greater complexity in the form of these activation volumes; however, these errors stayed below those obtained with common techniques for the solution of the direct problem, such as the use of activation functions [33]. Table 1 presents the pairwise statistical comparisons (p-values) of the final reconstruction error distributions, for each impedance and conductivity condition tested. For all impedances there are statistically significant differences in the performance of the proposed methodology between the isotropic and anisotropic conditions. Within the isotropic condition there are not statistically significant differences in the system’s performance as a function of the impedance tested, while for the anisotropic case the system generates errors that are statistically larger for the low-impedance scenario as compared to the medium and high cases.

Table 1 Pairwise statistical comparisons (p-values) of the error distributions obtained with the validation scheme, for each impedance and conductivity condition tested. The p-values highlighted in blue represent statistically significant differences in the performance of the proposed methodology for the involved conditions, at a significance level of 5% 

Figure 10 Box-whisker plots of the errors obtained with the validation scheme, discriminated by impedance condition for the isotropic (bottom row) and anisotropic (top row) conditions  

Finally, our group has developed an application that implements the proposed methodology for manipulating the VTA and obtaining the corresponding estimation parameters, together with a 3D visualization system of brain structures and the implanted stimulation electrode reconstructed from magnetic resonance and computed tomography images. Figure 11 shows the brain structures of interest in DBS for Parkinson’s disease, such as the subthalamic nucleus and the thalamus, in addition to the target VTAs and the VTAs obtained from the parameters estimated for one [Figure 11(a)) and two active contacts [Figure 11(b)). This implementation in C++ and Python 2.7, with all machine learning algorithms previously trained, allowed (once the desired VTA was defined) to execute tasks of visualization and deformation of the VTA, and an estimation of the DBS parameters in an mean time of less than 1.5 s in a computer of average specifications (Dell optiplex 990 with an Intel Core i7-2600 processor and 8 GB of RAM). This points to the feasibility of including the proposed methodology in practical applications.

Figure 11 The proposed methodology was implemented as part of a 3D visualization system of brain structures and the VTA. Here we show representations of the subthalamic nucleus, the thalamus, the DBS electrode, target VTAs (red) and VTAs computed from the parameters estimated for such targets (cyan) in the case where there is (a) a single active contact and (b) two active contacts with opposing polarities (anode-cathode)  

4. Discussion

One of the main challenges for the success of DBS therapy is to find a set of appropriate stimulation parameters that will result in a therapeutic benefit for the patient [11]. For a long time, this process of adjustment was done exclusively by trial and error. However, in the last decade computer models, and software packages based on them, have emerged to aid in this task [15, 17, 20]. Although these tools represent an important advance, their philosophy is still based on the manual search of the appropriate parameters. In view of the above, we proposed a methodology that allows to automatically find a configuration of stimulation parameters, for current clinical DBS devices, from a target VTA determined by the medical specialist.

The proposed methodology allows to first obtain a target VTA using a set of ellipses and their associated characteristics, by means of KNN approximation. From this deformed VTA and a given pulse width, the remaining stimulation parameters are obtained through classification and regression algorithms (SVMs) in conjunction with dimensionality reduction techniques (CKA). The deformation of the VTA by means of KNN allows translating and morphologically manipulating the target volume, not in an arbitrary way but following the restrictions imposed by the characteristics of the brain tissue and by the stimulation device. The need for this stage arises because the target VTA must conform with physical constrains that a random set of ellipses does not necessarily satisfy. While this may not be overly evident under an isotropic tissue conductivity assumption, when the VTA has nearly idealized ellipsoidal shapes, it becomes obvious under anisotropic conditions when a set of ellipses can only loosely represent the VTA. The use of CKA prior to the estimation of the stimulation parameters from the deformed VTA allows for the reduction of the data dimensionality while encoding its relevant information in order to decrease the computational time of the subsequent stage. Lastly, the generalization capability of the SVMs allows to handle non-linear changes in the VTA shape, in addition to the conversion between the dimensionally reduced VTAs and the values of the amplitude or the space of the polarities of the electrode contacts. On the whole, the proposed methodology generates errors with respect to the gold standard (once the VTA is re-calculated from the estimated parameters) in ranges comparable to those reported in the literature for other machine learning strategies that seek to solve only the traditional or direct problem of VTA estimation [33, 34].

4.1 Study limitations

Despite the positive results obtained we must emphasize that they are based on analyses performed for reduced ranges of stimulation parameters and for a single stimulation system (Medtronics ACTIVA-RC stimulator and 3,389 electrode). The selected ranges are based on the recommendations of clinical programming guidelines, especially concerning the amplitude of stimulation and the number of contacts simultaneously active [13]. The DBS device was selected because nowadays this family of stimulators is still the most widely implanted system [50]. Even so, the proposed methodology can be generalized to other parameter ranges as long as the corresponding datasets are available; the same observation applies to more complex DBS devices (e.g. devises based on current stimulation and with a higher number of contacts). Nonetheless, the generation of this new databases would entailed a large number of upfront simulations. Another limitation of this study is the fact that one of the stimulation parameters, the pulse width, was included as one of the system inputs. This was done because a large number of pulse width-amplitude combinations generate extremely similar VTAs, making the problem of interest very ill-conditioned. Taking into account that in clinical practice the pulse width tends to be set to a fixed value to then vary the other parameters of stimulation [13], we decided to establish this magnitude as system input to minimize the impact of the aforementioned hurdle. In the future, in order to include the pulse width in the parameters to be estimated, we propose the use of multiple output models. In these models, the machine learning systems would not deal with the stimulation parameters independently but simultaneously, this in turn should allow them to capture relations between the parameters that would otherwise go unnoticed. Additionally, restrictions regarding power efficiency [23, 39] could be included, this would give additional information to the system when selecting a suitable combination of stimulation parameters. Besides, a Bayesian optimization approach could be used to simultaneously tune all the free parameters of the proposed methodology, which could also allow capturing hidden patterns in the VTA data and to improve the overall performance of the system.

Finally, the models used to calculate the VTAs of the datasets are highly idealized [16]. Although the inclusion of the diffusion tensors in the calculation of the electric potential aims to improve that shortcoming their use along with axonal models of ordered, straight trajectories is not ideal [21, 25]. Furthermore, the approach used to obtain the diffusion tensors may lead to underestimated extents of stimulation [29]. However, as argued by [33], systems based on simple isotropic models have been shown to be useful in clinical practice [15, 23].

5. Conclusions

We have presented a simulation methodology based on machine learning techniques to estimate a DBS stimulation parameter setting for current clinical devices. The definition of a target VTA is accomplished through a deformation stage based on a KNN approximation that ensures biophysical compliance. The parameter estimation stage proper employs kernel methods, namely CKA and SVMs, to independently obtain the values of each of the DBS parameters under consideration. Since our methodology is based on data driven techniques it can be used for VTA data obtained under different assumptions, e.g. isotropic or anisotropic tissue conductivity conditions, and could be easily adapted to other DBS devices. Moreover, the trained algorithms have a low computational time. Also, an implementation of this methodology has been included as part of a VTA visualization and manipulation software. As future work, we will include the pulse width in the parameters to be estimated by using multiple output models and by taking into account power efficiency criteria, and we will extend our approach to experimental DBS electrode models.

6. Acknowledgements

This work was supported by project 1101-807-63808 titled “Caracterización morfológica de estructuras cerebrales por técnicas de imagen para el tratamiento mediante implantación quirúrgica de neuroestimuladores en la enfermedad de Parkinson” funded by Colciencias. Authors Viviana Gómez Orozco and Iván De La Pava Panche were supported by the program “Doctorado Nacional en Empresa - Convoctoria 758 de 2016”, also funded by Colciencias.

References

[1] D. Cunha and et al, “Toward sophisticated basal ganglia neuromodulation: Review on basal ganglia deep brain stimulation,” Neuroscience & Biobehavioral Reviews, vol. 58, pp. 186-210, Nov. 2015. [ Links ]

[2] C. C. McIntyre and R. W. Anderson, “Deep brain stimulation mechanisms: the control of network activity via neurochemistry modulation,” Journal of Neurochemistry, vol. 139, no. S1, pp. 338-345, 2016. [ Links ]

[3] K. L. Collins, E. M. Lehmann, and P. G. Patil, “Deep brain stimulation for movement disorders,” Neurobiology of Disease, vol. 38, no. 3, pp. 338-345, 2010. [ Links ]

[4] M. S. Okun, “Deep-brain stimulation for Parkinson’s disease,” New England Journal of Medicine, vol. 367, no. 16, pp. 1529-1538, 2012. [ Links ]

[5] X. Chen, Y. Xiong, G. Xu, and X. Liu, “Deep brain stimulation,” Interventional Neurology, vol. 1, no. 3-4, pp. 200-212, 2013. [ Links ]

[6] T. E. Schlaepfer, B. H. Bewernick, S. Kayser, B. Mädler, and V. A. Coenen, “Rapid effects of deep brain stimulation for treatment-resistant major depression,” Biological Psychiatry, vol. 73, no. 12, pp. 1204-1212, 2013. [ Links ]

[7] N. G. Laxpati, W. S. Kasoff, and R. E. Gross, “Deep brain stimulation for the treatment of epilepsy: circuits, targets, and trials,” Neurotherapeutics, vol. 11, no. 3, pp. 508-526, 2014. [ Links ]

[8] D. R. Cleary, A. Ozpinar, A. M. Raslan, and A. L. Ko, “Deep brain stimulation for psychiatric disorders: where we are now,” Neurosurgical Focus, vol. 38, no. 6, p. E2, 2015. [ Links ]

[9] A. Veerakumar and O. Berton, “Cellular mechanisms of Deep brain stimulation: activity-dependent focal circuit reprogramming?” Current Opinion in Behavioral Sciences, vol. 4, pp. 48-55, 2015. [ Links ]

[10] A. A. Kühn and J. Volkmann, “Innovations in deep brain stimulation methodology,” Movement Disorders, vol. 32, no. 1, pp. 11-19, 2017. [ Links ]

[11] S. Miocinovic, S. Somayajula, S. Chitnis, and J. L. Vitek, “History, applications, and mechanisms of deep brain stimulation,” JAMA neurology, vol. 70, no. 2, pp. 163-171, 2013. [ Links ]

[12] K. Hunka, O. Suchowersky, S. Wood, L. Derwent, and Z. H. Kiss, “Nursing time to program and assess deep brain stimulators in movement disorder patients,” Journal of Neuroscience Nursing, vol. 37, no. 4, p. 204, 2005. [ Links ]

[13] J. Volkmann , E. Moro, and R. Pahwa, “Basic algorithms for the programming of deep brain stimulation in Parkinson’s disease,” Movement Disorders , vol. 21, no. S4, pp. S284-S289, 2006. [ Links ]

[14] M. Astrom and et al, “Method for patient-specific finite element modeling and simulation of deep brain stimulation,” Medical & Biological Engineering & Computing, vol. 47, no. 1, pp. 21-28, 2009. [ Links ]

[15] M. H. Pourfar and et al, “Model-based deep brain stimulation programming for Parkinson’s disease: the GUIDE pilot study,” Stereotactic and Functional Neurosurgery, vol. 93, no. 4, pp. 231-239, 2015. [ Links ]

[16] C. R. Butson andC. C. McIntyre , “Tissue and electrode capacitance reduce neural activation volumes during deep brain stimulation,” Clinical Neurophysiology, vol. 116, no. 10, October 2005. [Online]. Available: https://doi.org/10.1016/j.clinph.2005.06.023Links ]

[17] C. Butson, A. Noecker, C. Maks, andC. C. McIntyre , “StimExplorer: deep brain stimulation parameter selection software system,” Acta Neurochirurgica. Supplement, vol. 97, no. pt2, pp. 569-74, 2007. [ Links ]

[18] N. Yousif, R. Borisyuk, N. Pavese, D. Nandi, and P. Bain, “Spatiotemporal visualization of deep brain stimulation-induced effects in the subthalamic nucleus,” European Journal of Neuroscience, vol. 36, no. 2, July 2012. [Online]. Available: https://doi.org/10.1111/j.1460-9568.2012.08086.xLinks ]

[19] N. Yousif and et al, “Evaluating the impact of the deep brain stimulation induced electric field on subthalamic neurons: A computational modelling study,” Journal of Neuroscience Methods, vol. 188, no. 1, April 30 2010. [Online]. Available: https://doi.org/10. 1016/j.jneumeth.2010.01.026Links ]

[20] P. M. Lauro and et al, “DBSproc: an open source process for DBS electrode localization and tractographic analysis,” Human Brain Mapping, vol. 37, no. 1, pp. 422-433, 2016. [ Links ]

[21] C. R. Butson , S. E. Cooper, J. M. Henderson, andC. C. McIntyre , “Patient-specific analysis of the volume of tissue activated during deep brain stimulation,” Neuroimage, vol. 34, no. 2, January 15 2007. [Online]. Available: https://doi.org/10.1016/j.neuroimage.2006.09.034Links ]

[22] T. A. Dembek and et al, “Probabilistic mapping of deep brain stimulation effects in essential tremor,” NeuroImage: Clinical, vol. 13, pp. 164-173, 2017. [ Links ]

[23] A. M. Frankemolle and et al, “Reversing cognitive-motor impairments in Parkinson’s disease patients using a computational modelling approach to deep brain stimulation programming,” Neuroimage, vol. 133, no. 3, March 2010. [Online]. Available: https://doi.org/10.1093/brain/awp315Links ]

[24] C. R. Butson andC. C. McIntyre , “Current steering to control the volume of tissue activated during deep brain stimulation,” Brain Stimulation, vol. 1, no. 1, pp. 7-15, 2008. [ Links ]

[25] A. Chaturvedi, C. R. Butson , S. F. Lempka, S. E. Cooper , andC. C. McIntyre , “Patient-specific models of deep brain stimulation: influence of field model complexity on neural activation predictions,” Brain Stimulation , vol. 3, no. 2, pp. 65-77, 2010. [ Links ]

[26] G. Walckiers, B. Fuchs, J. P. Thiran, J. R. Mosig, and C. Pollo, “Influence of the implanted pulse generator as reference electrode in finite element model of monopolar deep brain stimulation,” Journal of Neuroscience Methods , vol. 186, no. 1, pp. 90-96, 2010. [ Links ]

[27] T. C. Zhang and W. M. Grill, “Modeling deep brain stimulation: point source approximation versus realistic representation of the electrode,” Journal of Neural Engineering, vol. 7, no. 6, December 2010. [Online]. Available: https://doi.org/10.1088/1741-2560/7/6/066009Links ]

[28] K. J. Van Dijk and et al, “A novel lead design enables selective Deep brain stimulation of neural populations in the subthalamic region,” Journal of Neural Engineering, vol. 12, no. 4, p. 046003, 2015. [ Links ]

[29] B. Howell and C. C. McIntyre, “Analyzing the tradeoff between electrical complexity and accuracy in patient-specific computational models of deep brain stimulation,” Journal of Neural Engineering , vol. 13, no. 3, p. 036023, 2016. [ Links ]

[30] C. R. Butson andC. C. McIntyre , “Role of electrode design on the volume of tissue activated during deep brain stimulation,” Journal of Neural Engineering , vol. 3, no. 1, March 2006. [Online]. Available: https://doi.org/10.1088/1741-2560/3/1/001Links ]

[31] E. Peterson, O. Izad, and D. J. Tyler, “Predicting myelinated axon activation using spatial characteristics of the extracellular field,” Journal of Neural Engineering , vol. 8, no. 4, p. 046030, 2011. [ Links ]

[32] M. Astrom, E. Diczfalusy, H. Martens, and K. Wardell, “Relationship between neural activation and electric field distribution during deep brain stimulation,” IEEE Trans. Biomed. Eng., vol. 62, no. 2, February 2015. [Online]. Available: https://doi.org/10.1109/TBME.2014.2363494. [ Links ]

[33] A. Chaturvedi, J. L. Luján, andC. C. McIntyre , “Artificial neural network based characterization of the volume of tissue activated during deep brain stimulation,” Journal of Neural Engineering , vol. 10, no. 5, October 2013. [Online]. Available: https://doi.org/10.1088/1741-2560/10/5/056023Links ]

[34] I. De La Pava and et al, “A hierarchical K-nearest neighbor approach for volume of tissue activated estimation,” in Iberoamerican Congress on Pattern Recognition, Lima, Perú, 2016, pp. 125-133. [ Links ]

[35] M. F. Contarino and et al, “Directional steering: A novel approach to deep brain stimulation,” Neurology, vol. 83, no. 13, pp. 1163-1169, 2014. [ Links ]

[36] M. Hariz, “Deep brain stimulation: new techniques,” Parkinsonism & Related Disorders, vol. 20, no. 13, pp. S192-S196, 2014. [ Links ]

[37] C. Pollo and et al, “Directional deep brain stimulation: an intraoperative double-blind pilot study,” Brain, vol. 137, no. pt7, pp. 2015-2026, 2014. [ Links ]

[38] V. Gómez and et al, “A kernel-based approach for DBS parameter estimation,” in Iberoamerican Congress on Pattern Recognition, Lima, Perú, 2016, pp. 158-166. [ Links ]

[39] E. Peña, S. Zhang, S. Deyo, Y. Xiao, and M. D. Johnson, “Particle swarm optimization for programming deep brain stimulation arrays,” Journal of Neural Engineering , vol. 14, no. 1, p. 016014, 2017. [ Links ]

[40] C. Cortes, M. Mohri, and A. Rostamizadeh, “Algorithms for learning kernels based on centered alignment,” Journal of Machine Learning Research, vol. 13, no. 1, pp. 795-828, 2012. [ Links ]

[41] D. Cárdenas, D. Collazos, A. Álvarez, and G. Castellanos, “Algorithms for learning kernels based on centered alignment,” Engineering Applications of Artificial Intelligence, vol. 68, pp. 10-17, 2018. [ Links ]

[42] B. Schölkopf and A. J. Smola, Learning with kernels: support vector machines, regularization, optimization, and beyond. Cambridge, USA: MIT press, 2001. [ Links ]

[43] A. Machado and et al, “Deep brain stimulation for Parkinson’s disease: surgical technique and perioperative management, Movement Disorders, vol. 21, no. S14, pp. S247-S258, 2006. [ Links ]

[44] L. C. Chang, D. K. Jones, and C. Pierpaoli, “ARESTORE: robust estimation of tensors by outlier rejection,” Magnetic Resonance in Medicine, vol. 53, no. 5, May 2005. [Online]. Available: https://doi.org/10.1002/mrm.20426Links ]

[45] D. S. Tuch, V. J. Wedeen, A. M. Dale, J. S. George, and J. W. Belliveau, “Conductivity tensor mapping of the human brain using diffusion tensor MRI,” in Proceedings of the National Academy of Sciences, USA, 2001, pp. 11 697-11 701. [ Links ]

[46] O. C. Zienkiewicz, R. L. Taylor, and J. Zhu, Finite Element Method: Its Basis and Fundamentals. Elsevier, Incorporated, 2013. [ Links ]

[47] C. C. McIntyre , A. G. Richardson, and W. M. Grill, “Modeling the excitability of mammalian nerve fibers: influence of afterpotentials on the recovery cycle,” Journal of Neurophysiology, vol. 87, no. 2, pp. 995-1006, 2002. [ Links ]

[48] M. L. Hines, A. P. Davison, and E. Muller, “NEURON and Python,” Frontiers in Neuroinformatics, January 28 2009. [Online]. Available: https://doi.org/10.3389/neuro.11.001.2009. [ Links ]

[49] J. Pizarro, E. Guerrero, and P. L. Galindo, “Multiple comparison procedures applied to model selection,” Neurocomputing, vol. 48, no. 1-4, pp. 155-173, 2002. [ Links ]

[50] M. Okun and P. Zeilman, “Parkinson’s disease: Guide to Deep brain stimulation therapy. Hagerstown, MD: National Parkinson Foundation,” 2017. [ Links ]

Received: February 04, 2019; Accepted: July 05, 2019

* Correspondig author: Viviana Gómez Orozco e-mail: vigomez@utp.edu.co

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License