SciELO - Scientific Electronic Library Online

 
 número84Supply chain design using a modified IWD algorithmBiocompatibility of bismuth silicate coatings deposited on 316L stainless steel by sol-gel process índice de autoresíndice de assuntospesquisa de artigos
Home Pagelista alfabética de periódicos  

Serviços Personalizados

Journal

Artigo

Indicadores

Links relacionados

  • Em processo de indexaçãoCitado por Google
  • Não possue artigos similaresSimilares em SciELO
  • Em processo de indexaçãoSimilares em Google

Compartilhar


Revista Facultad de Ingeniería Universidad de Antioquia

versão impressa ISSN 0120-6230

Rev.fac.ing.univ. Antioquia  no.84 Medellín jul./set. 2017

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

Articles

Accelerating the computation of the volume of tissue activated during deep brain stimulation using Gaussian processes

Aceleración del cálculo del volumen de tejido activo durante estimulación cerebral profunda usando procesos gaussianos

Iván De La Pava Panche1 

Viviana Gómez-Orozco1  * 

Mauricio Alexander Álvarez-López2 

Óscar Alberto Henao-Gallo1 

Genaro Daza-Santacoloma3 

Álvaro Ángel Orozco-Gutiérrez1 

1Grupo de Automática, Facultad de Ingenierías, Universidad Tecnológica de Pereira. Carrera 27 # 10-02 Barrio Álamos. A. A. 97. Pereira, Colombia.

2Machine Learning group, Department of Computer Science, University of Shefield. Shefield, UK.

3Instituto de Epilepsia y Parkinson del Eje Cafetero - Neurocentro. Carrera 9 # 25-20 Torre A Consultorio 418. C. P. 660002. Pereira, Colombia.


ABSTRACT

The volume of tissue activated (VTA) is a well-established approach to model the direct effect of deep brain stimulation (DBS) on neural tissue. Previous studies have pointed to its potential clinical applications. However, the elevated computational runtime required to estimate the VTA with standard techniques used in biological neural modeling limits its suitability for practical use. The goal of this study was to develop a novel methodology to reduce the computation time of VTA estimation. To that end, we built a Gaussian process emulator. It combines multicompartment axon models coupled to the stimulating electric field with a Gaussian process classifier (GPC), following the premise that computing the VTA from a field of axons is in essence a binary classification problem. We achieved a considerable reduction in the average time required to estimate the VTA, under both ideal isotropic and realistic anisotropic brain tissue conductive conditions, limiting the loss of accuracy and overcoming other drawbacks entailed by alternative methods.

Keywords: Deep brain stimulation; volume of tissue activated; multicompartment axon model; emulation; Gaussian process classification

RESUMEN

El volumen de tejido activo (VTA) es un enfoque bien establecido para modelar los efectos directos de la estimulación cerebral profunda en el tejido neuronal. Estudios previos han señalado sus posibles aplicaciones clínicas. Sin embargo, el elevado costo computacional requerido para estimar el VTA con las técnicas estándar utilizadas en el modelado neuronal biológico limita su usabilidad a nivel práctico. El objetivo de este estudio fue desarrollar una metodología novedosa para reducir el tiempo de cálculo en la estimación del VTA. Con ese fin, se construyó un emulador basado en procesos gaussianos. Este combina modelos axonales de múltiples compartimientos acoplados al campo de estimulación eléctrica con un sistema de clasificación basado en procesos gaussianos, siguiendo la premisa de que calcular el VTA a partir de un campo axonal es en esencia un problema de clasificación binaria. Se logró una reducción considerable en el tiempo promedio requerido para estimar el VTA, tanto bajo condiciones de conductividad isotrópica idealizada como bajo condiciones realistas de conductividad anisotrópica, limitando la perdida de precisión y superando otros inconvenientes presentes en métodos alternativos.

Palabras clave: Estimulación cerebral profunda; volumen de tejido activo; modelo axonal de múltiples compartimientos; emulación; clasificación de procesos gaussianos

1. Introduction

Deep brain stimulation (DBS) is a surgical technique used mainly to treat movement disorders, such as Parkinson’s disease, essential tremor, and dystonia, in patients whose symptoms cannot be appropriately controlled with drugs. It involves the application of high frequency electric pulses to a target region in the basal ganglia, the thalamus, or other subcortical structures, via an implanted electrode; this stimulation leads to symptom improvement. Although DBS is widely practiced, its mechanisms of action are still not completely clear, and most of the insight gained in recent years has come from computer simulations [1-5]. In this context, the volume of tissue activated (VTA), defined as the spatial spread of direct neural activation in response to electrical stimulation, is a measure that allows for a computational assessment of the impact of DBS [6-9].

The gold standard for VTA estimation is to couple the electric potential generated by the DBS electrode with a model of multicompartment axons arranged in a field around the electrode shaft [6]. Axons that, as a result of the stimulation, fire action potentials in a one-to-one ratio with the DBS pulses are considered active, and their spatial distribution defines the VTA. The downside of this approach is its elevated computational cost, when considered in the context of its potential clinical applications [10-12]. This obstacle has led to simplified methods that minimize the use of such models by exploiting the relationship between the spatial location of the axons and the electrical stimulation. Such relationship is captured in the form of activation threshold curves, fitted from data obtained from multicompartment axon models, that are then used to estimate the VTA [10-13]. However, these curves do not reproduce the results of multicompartment models accurately and cannot be applied successfully when multiple contacts of the DBS electrode are active [12]. A possible solution to these problems is proposed by Chaturvedi et al. [12]. It consists in training an artificial neuronal network, with the DBS stimulation parameters as inputs and the elliptic profiles defined by the active axons as outputs. Once trained, the neural network can estimate the VTA for any combination of stimulation parameters. However, this methodology only works under the assumption of isotropic tissue conductivity, which affects the size and shape of the VTA [7].

In this paper, we propose an alternative approach to reduce the computational runtime of VTA estimation and overcome the limitations of the methods described above. Active axons fire an action potential for each stimulation pulse delivered by the stimulator. Inactive axons fail to respond in such a way. It follows that active and inactive axons can be thought of as part of two different classes. Furthermore, all axons are independent from one another. Under these premises, determining the volume of tissue activated from a field of axons is equivalent to a binary classification problem. The output of the standard VTA model can be then represented statistically and emulated. The emulation of computer simulations circumvents the problems posed by complex and computationally intensive models by building statistical representations of them used to address the issue under study without additional runs of the original simulation. The main approach to developing emulators uses Gaussian processes [14] [15]. Our aim was then to train a Gaussian process classifier (GPC) to determine whether an axon at a given position in space was active due to DBS, and by doing so, to estimate the VTA. This work is an extended version of a study presented in the 7th Iberian Conference on Pattern Recognition and Image Analysis [16]. The present version contains a more detailed theoretical framework, and a larger experimental setup that includes the estimation of the VTA when realistic anisotropic brain tissue conductivity conditions are considered.

2. Materials and methods

2.1. Electric propagation in the brain tissue

The calculation of the electric potential generated by the DBS electrode was carried out in the finite element method (FEM) software COMSOL Multiphysics 4.2. A simplified 3D model of a clinical electrode (Medtronic DBS 3389 electrode) positioned in the middle of a conductive extracellular medium was built. The simplified electrode model consisted of four conductive contacts (4x106 Sm-1) 1.27 mm in width and 1.5 mm in height separated by insulating bands (1x10-10 Sm-1) 0.5 mm in height, and of an insulating semicircular tip with radius 0.635 mm (Figure 1(a)) [17]. The brain tissue was modeled as a sphere of diameter 10 cm. First, a bulk tissue conductivity of 0.3 Sm-1 was used to account for the isotropic assumption. Then, diffusion tensor based conductivities were used to represent the tissue anisotropy of the basal ganglia region. The diffusion tensors were estimated from the DTI30 dataset, with the RESTORE (Robust Estimation of Tensors by Outlier Rejection) algorithm [18], and then linearly transformed to conductivity tensors [19]. In both cases, a representation of a 0.5 mm encapsulation layer around the electrode was included, and its conductivity was set to three different values (0.680 Sm-1, 0.128 Sm-1, 0.066 Sm-1) to represent low (~500 Ω), medium (~900 Ω), and high (~1,500 Ω) impedance conditions [12-20].

The model also integrated the voltage drop at the electrode-tissue interface. The boundary conditions were the same used by Yousif et al. [8]. The stimulation was modeled by imposing Dirichlet boundary conditions on the active contacts of the electrode. Dirichlet boundary conditions were also used to ground the boundaries of the extracellular medium. Zero current flow conditions were imposed on the surfaces of the inactive contacts and insulating components of the electrode. Finally, Poisson’s equation (Eq. (1)) was solved to obtain the voltage propagation in the brain tissue.

Figure 1 (a) DBS electrode model. The electrode contacts were numbered from 0 to 3. (b) Representation of the orientation of the axons around the electrode shaft 

(1)

where is the conductivity of tissue medium, and is the extracellular potential.

2.2. Axonal distribution and VTA estimation

Multicompartment neuron models allow for the computation of temporal and spatial changes in a neuron’s transmembrane potential. They can capture the effects of an external electrical stimulus, by solving the cable equation, Eq. (2), for each compartment in each neuron of the model:

(2)

Where is the transmembrane potential at the compartment, is the specific membrane capacitance, is the axial resistance per unit length, d is the compartment diameter, is the second spatial differences of the transmembrane potential is the second spatial differences of the extracellular potential , and is the total ionic current flowing through the membrane at a given moment in time. The term is the compartment length [21-23].

In the VTA modeling literature, multicompartment neuron models are often restricted to axonal models, since axons have been shown to drive the basic neural response to external electrical stimulation. In addition to this behavior, the fundamental biophysics of how axonal response to external electrical stimulation works is independent of neuron type [24-26]. The axonal model used in this work corresponds to the multicompartment myelinated axon model detailed in [27]. Briefly, each axon includes 21 nodes of Ranvier, 2 myelin attachment segments, 2 paranode main segments, and 3 internode segments between each node. The nodes are modeled by a parallel combination of the membrane capacitance with nonlinear conductances (fast Na+, persistent Na+, and slow K+) and a linear leakage conductance. The paranodal and internodal compartments include two concentric layers, each including a linear conductance in parallel with the membrane capacitance, to represent the myelin sheath and the axolemma.

To estimate the volume of tissue activated, a population of 8,112 straight axons of diameter 5.7μm was built around the electrode shaft. The axonal fibers were oriented in four different directions, perpendicular to the axis of the electrode (Figure 1(b)), and with a distance between axons of 0.5 mm in both the vertical and horizontal directions. This axonal distribution aimed to capture the spatial spread of activation in response to the stimulus in several directions around the electrode, considering that under anisotropic conductivity conditions the electric potential would change as a function of the direction of propagation [28].

To simulate the electrical stimulation of the neuronal tissue, the values of the electric potential were linearly interpolated from the nodes of the FEM solution onto each of the axonal model compartments. An axon was considered active if it fired in a one-to-one ratio with a DBS stimulation pulse. The VTA was defined as the volume enclosed by the locations of the central nodes of Ranvier of the activated axons [28]. The simulations of the axonal response to the electric stimulation were implemented in NEURON 7.3 configured as a Python module [29] [30].

The methodology described here and in section 2.1 corresponds to the gold standard for VTA estimation. A full implementation of it was used to compute volumes of tissue activated, for several combinations of stimulation parameter settings (see Table 1], that would serve as reference data sets to evaluate the performance of the proposed emulator. The base parameters were varied, one at a time, in a range of possible values for the three active contact configurations considered: one cathode (monopolar), cathode-cathode (bipolar {-}), and anode-cathode (bipolar {±}). Every set of stimulation parameters and encapsulation tissue impedances was used under both isotropic and anisotropic conductivity conditions.

Table 1 Realistic stimulation parameter settings used in this study (Medtronic ACTIVA-RC stimulator) 

2.3. Gaussian process emulation of the VTA

Gaussian process classification

A Gaussian process (GP) is defined as a probability distribution over functions such that the values of said functions evaluated at an arbitrary set of points jointly have a Gaussian distribution, that is, . The mean vector and covariance matrix completely specify the GP, and are obtained from associated mean and covariance (kernel) functions, and respectively, that may additionally depend on a set of hyperparameters θ [31-34].

In general, the idea behind binary GP classification works as follows: given a set of inputs and their associated labels , with , the probability is represented by the latent function that describes the data mapped to the unit interval through a sigmoid function . If is symmetric then , and the joint likelihood where , can be expressed as Eq. (3)

(3)

Imposing a Gaussian prior over and assuming a zero mean function , the posterior probability (Eq. (4)) becomes:

(4)

where is the observed data. Then, the class probabilities for a new data pair (Eq. (5)) can be obtained marginalizing

(5)

and computing the expectation as shown in Eq. (6)

(6)

In a Bayesian framework, GPs offer the advantage of analytical tractability: a Gaussian prior (the assumption made before observing the data) and a Gaussian likelihood (how probable the data set is given the assumption) will result in a Gaussian posterior (the uncertainty in the assumption after the data is observed), and even though the discrete nature of the classification problem gives rise to a non-Gaussian likelihood (Eq. (4)), a posterior distribution can still be approximated by methods such as Laplace approximation or Expectation Propagation. In this study, we use the Laplace approximation because of its lower computational cost [32].

Laplace’s method approximates the posterior doing a second order Taylor expansion of around the maximum of the unnormalized posterior as shown in Eq. (7)

(7)

Where and . The hyperparameters θ of the covariance function are selected so that they maximize the evidence or marginal likelihood . From the results of Laplace’s method the logarithm of this quantity can be approximated as Eq. (8)

(8)

Procedure

A random sample of axons was taken from the total axonal population. The axons were sampled from a uniform distribution with respect to the Euclidean distance between their central nodes of Ranvier and the center of the active contact. When more than one contact was active, the above procedure was repeated, sampling, for each of them, a number of axons equal to divided by the number of active contacts. Next, a multicompartment simulation was executed to determine which of the sampled neural fibers were active during the stimulation. The information provided by the multicompartment simulation was converted to a set of labeled data , where is a vector of characteristics for each axon or data point distance between the central nodes of Ranvier and the axis of the DBS electrode. distance between the central nodes of Ranvier and the plane containing the electrode tip and that is perpendicular to the electrode axis. value of the electric potential (V) at each central node of Ranvier), and is its corresponding label (active axon:, inactive axon:). These data were used to train a Gaussian process-based classifier. All inputs to the classification algorithm (pyGPs library for Python 2.7) are depicted in Figure 2; they were selected so that the relationship between activation, the position of each axon with respect to the electrode shaft, and the electric potential was represented in the classification outcome.

Figure 2 Inputs to the classification algorithm 

The classifier was trained with 70% of the data obtained from the multicompartment axon model, assuming a logistic likelihood function prescribing a general purpose kernel (squared exponential covariance function with automatic relevance determination), and using Laplace’s approximation to the posterior. The hyperparameters of the kernel were selected minimizing the Laplace’s approximation of the negative log marginal likelihood. The inference method was chosen due to computational cost considerations [34].

Since a GP classifier assigns class probabilities instead of labels, a cut off probability was determined with a receiver-operating characteristic (ROC) analysis performed on all available data from the multicompartment model. Then, the VTA was estimated by predicting which of the 8,112 central nodes of Ranvier of the entire axonal population, and therefore the axonal fibers they belong to, would be activated by the applied stimulus. Finally, the spreads of activation predicted by the classifier, in 30 independent simulation runs, for each of the stimulation parameters shown in Table 1 were compared against those obtained with activation threshold curves customized to our data [10-13], in terms of the prediction error (Eq. (9)).

(9)

Where FP and FN are the false positives and false negatives with respect to the reference data set, and AA is the reference number of active axons.

3. Results and discussion

In this work, we explored a variation to the standard approach to VTA estimation: the VTA is defined by the axons activated for a given set of stimulation parameters (see Figure 3). These axons constitute a subset of the axonal field that satisfies certain conditions that differentiate its elements from the rest of the axonal population; active and inactive axons can be viewed as belonging to two different classes.

Figure 3 (a) Locations of the central nodes of Ranvier (dots) of the axons activated by a bipolar stimulation pulse. (b) The VTA is the volume defined by the spatial distribution of the active axons 

Based on this premise, we aimed to develop a classification system that could reduce the computation time required for VTA estimation, without the shortcomings of alternative methods like activation threshold curves. Figure 4 summarizes our approach. Instead of a large model with enough axons to describe accurately the region of interest around the electrode shaft, we sampled a comparatively small number of axons from the total axonal population, and with this reduced model we executed a full multicompartment model simulation in order to obtain labeled data to train a GP classifier (Figure 4 (a)). Then, we used the classifier to predict which elements of the entire axonal population would be active (Figure 4 (b)), obtaining an estimation of the VTA (Figure 4 (c)).

Figure 4 Scheme of the proposed Gaussian process emulation of the VTA 

The reduction in the computational runtime comes from the use of fewer axons in the multicompartment axon model. Figure 5 shows how the simulation runtime grows linearly with the number of sampled axons (n). These average values were obtained for VTAs estimated for 20 randomly selected combinations of stimulation parameters from Table 1, assuming either isotropic or anisotropic conductivity conditions, and holding the pulse width at 90 μs. Besides, for the same data set, the prediction errors (Eq. (9)) stabilize after the number of axons reaches 250. was chosen to evaluate the performance of the proposed method because it is past that point and falls in the range between 5% and 10% of the total axonal population. The average time to estimate the VTA for all trials and stimulation parameter settings, sampling 500 axons, was 267.7 s. This figure represents a reduction of 92.3% compared with an average time of 3,467.9 s required to solve the model described in section 2.2 [6] (models running in a Dell OptiPlex 990 with an Intel Core i7-2,600 processor, and 8 GB RAM.).

Figure 5 Box-whisker plots of the prediction errors and average simulation runtimes (mean±SD) as functions of the number of axons used in the classification procedure 

Activation threshold curves can achieve simulation runtimes of under 1 s, that is, if they are not fitted every time the stimulation parameters vary. Otherwise, it is necessary to run full multicompartment axon models to obtain new data to fit them. Whatever the case, they are prone to error [12]. Figure 6 shows the total errors generated by the classifier and by threshold curves with respect to reference data sets for the monopolar configuration and several combinations of stimulation parameters. The Gaussian process emulator easily outperformed the threshold curves, independently of whether isotropic or anisotropic conductivities were assumed, although with larger uncertainties in the latter case. Also, unlike the curves, the proposed emulator allows for the estimation of the VTA when more than one of the electrode contacts are active.

Figure 6 Prediction errors for the monopolar configuration: (a) amplitude, (b) pulse width, and (c) impedance variations 

Figure 7 shows four different volumes of tissue activated for two different bipolar stimulation parameter settings, but under different conductivity conditions: (a) anode-cathode configuration (bipolar{±}) under isotropic tissue conductivity, (b) anode-cathode configuration under anisotropic tissue conductivity, (c) cathode-cathode configuration (bipolar{-}) under isotropic tissue conductivity, (d) cathode-cathode configuration under anisotropic tissue conductivity. This shows the ability of the emulator not only to work with multi-contact stimulation, but also to represent the effects that tissue conductivity has on the VTA. The tissue conductivity affects the VTA via changes in the propagation of the stimulus waveform and consequently on axonal activation. The volumes of tissue activated when isotropic conductivity is assumed exhibit radial symmetry with respect to the electrode axis, as expected. Such symmetry disappears for the anisotropic case. In terms of the prediction error, the two bipolar configurations studied, follow distinct trends (see Figure 8). While the error for the cathode-cathode configuration remains low, stable, and with similar differences between the isotropic and anisotropic scenarios across the range of stimulation parameter settings simulated, the error for the anode-cathode configuration is higher and presents large variations between tissue conductivity conditions. Figure 9 shows these tendencies more clearly. More than the effect of anisotropy, it is having two spatially separate regions of activation what affects the most the performance of our method. Median errors for the anode-cathode configuration, which produces independent volumes of tissue activated for each active contact, double those of the other two configurations. However, even in this worst case scenario, the GP emulator matches the performance of activation threshold curves for the simpler monopolar case.

Figure 7 Volumes of tissue activated estimated by the classifier for the two bipolar configurations used, under different tissue conductivity conditions 

Figure 8 Box-whisker plots of the prediction errors for the two bipolar configurations used: (a) amplitude, (b) pulse width, and (c) impedance variations 

Figure 9 Box-whisker plots of the GP classifier prediction errors across all stimulation parameter settings, under isotropic and anisotropic tissue conductivity conditions 

At this point, we must note that the error metric used to assess the performance of the emulator is geared towards classification problems, and it is in that sense appropriate for our case. However, it is only partially bounded and, more importantly, it does not account for the morphology of the VTA. The problem of directly accounting for morphological errors in the estimated VTA, which we think could have more practical meaning, is a pending task and constitutes a future line of work.

A second important point is that our results were obtained applying a single DBS stimulation pulse to the multicompartment axon models. The VTA criterion is clear in that an axon is considered active if it fires an action potential for each DBS pulse, but there is less clarity about the number of pulses that should be applied to verify that condition. Miocinovic et al. [2] report that they did not observe any changes in their results after applying a train of 25 DBS pulses to their experimental setting. Yousif et al. [8] used up to 12 pulses in their experiments. In our case, increasing the number of pulses applied to the reference model, described in section 2.2, did not produce any observable changes in our results. We then set up a smaller axonal field of 738 fibers arranged as in [6] but with a higher axonal density (0.25 mm separation between axons in the horizontal and vertical directions). This new set up revealed small changes in the periphery of the VTA as the number of applied pulses increased (Figure 10). However, the combined effect of increasing the time of simulation, to represent the delivery of a train of pulses, and the number of axons, to obtain the same axonal density in the entire region of interest around the DBS electrode, would result in a disproportional increase in the computational burden compared to the benefits of a more detailed model. Our initial set-up was detailed enough to reproduce successfully experimental results presented elsewhere [6] [8].

Figure 10 Two-dimensional representation of the changes in the extension of the VTA due to the number of applied DBS pulses used for its estimation 

We must also highlight the fact that the number of parameter settings used in this study is only a small fraction of the thousands of possible combinations of stimulation parameters offered by DBS devices. The proposed methodology is intended to work for any of such combinations, but we chose to focus on the parameter ranges advised in clinical guidelines [35], especially in regard to the stimulation amplitude and to the number of active contacts. Our approach also assumes highly idealized axonal orientations and trajectories. Nonetheless, models based on these simplifications have yielded positive results in studies that compared their predictions of the VTA with experimental data from Parkinson’s disease patients implanted with DBS systems [11]. A last important point is that although the computational runtime of our method is much shorter than that of full multicompartment axon models, it is still high for real-time applications. Further reductions in the computation time required for VTA estimation can be achieved by optimizing the use of computational resources, e.g., running a version of the code that parallelized the execution of NEURON [36], once per each of the stimulation parameters shown in Table 1, reduced the average simulation runtime by an additional 22.6%.

5. Conclusions

In this study, we developed a new methodology to reduce the computation time required to estimate the volume of tissue activated (VTA) during deep brain stimulation (DBS). We built a Gaussian process emulator that combined multicompartment axon models coupled to the stimulating electric field with a Gaussian process classifier. Our approach cut down to a tenth the average computational runtime of VTA estimation compared with the gold standard. In addition, it allowed for the estimation of the VTA for non-monopolar stimulation, overcoming a limitation of commonly used alternatives, such as activation threshold curves. It also outperformed these, in terms of their prediction errors, for all the combinations of parameter settings studied where direct comparison was possible. Finally, the emulator also succeeded in estimating the VTA when realistic anisotropic brain tissue conductivities were included in the simulation model.

6. Acknowledgments

This work was supported by the project 1110-744-55860, and authors Iván De La Pava and Viviana Gómez were supported by the program “Doctorado Nacional en Empresa - Convoctoria 758 de 2016”, both funded by Colciencias. We would like to thank Hernan Dario Vargas Cardona for providing us the brain conductivity tensors used in this study.

7. References

1. C. McIntyre, W. Grill, D. Sherman, and N. Thakor, “Cellular effects of deep brain stimulation: model-based analysis of activation and inhibition,” Journal of Neurophysiology, vol. 91, no. 4, pp. 1457-1469, 2004. [ Links ]

2. S. Miocinovic et al., “Computational analysis of subthalamic nucleus and lenticular fasciculus activation during therapeutic deep brain stimulation,” Journal of neurophysiology, vol. 96, no. 3, pp. 1569-1580, 2006. [ Links ]

3. C. McIntyre and P. Hahn, “Network perspectives on the mechanisms of deep brain stimulation,” Neurobiology of Disease, vol. 38, no. 3, pp. 329-337, 2010. [ Links ]

4. R. So, A. Kent, and W. Grill, “Relative contributions of local cell and passing fiber activation and silencing to changes in thalamic fidelity during deep brain stimulation and lesioning: A computational modeling study,” Journal of Computational Neuroscience, vol. 32, no. 3, pp. 499-519, 2012. [ Links ]

5. M. Birdno, W. Tang, J. Dostrovsky, W. Hutchison, and W. Grill , “Response of human thalamic neurons to high-frequency stimulation,” PLoS One, vol. 9, no. 5, pp. 1-10, 2014. [ Links ]

6. C. Butson andC. McIntyre , “Tissue and electrode capacitance reduce neural activation volumes during deep brain stimulation,” Clinical Neurophysiology, vol. 116, no. 10, pp. 2490-2500, 2005. [ Links ]

7. C. Butson , S. Cooper, J. Henderson, and C. McIntyre , “Patient-specific analysis of the volume of tissue activated during deep brain stimulation,” Neuroimage, vol. 34, no. 2, pp. 661-670, 2007. [ Links ]

8. N. Yousif 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, pp. 105-112, 2010. [ Links ]

9. 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, pp. 2252-2259, 2012. [ Links ]

10. M. Åström, E. Diczfalusy, H. Martens, and K. Wårdell, “Relationship between neural activation and electric field distribution during deep brain stimulation,” IEEE Trans. Biomed. Eng., vol. 62, no. 2, pp. 664-672, 2015. [ Links ]

11. A. Frankemolle et al., “Reversing cognitive-motor impairments in Parkinson’s disease patients using a computational modelling approach to deep brain stimulation programming,” Brain, vol. 133, no. 3, pp. 746-761, 2010. [ Links ]

12. A. Chaturvedi, J. Luján, andC. McIntyre , “Artificial neural network based characterization of the volume of tissue activated during deep brain stimulation,” Journal of neural engineering, vol. 10, no. 5, pp. 1-8, 2013. [ Links ]

13. C. Butson andC. McIntyre , “Role of electrode design on the volume of tissue activated during deep brain stimulation,” Journal of neural engineering , vol. 3, no. 1, pp. 1-8, 2006. [ Links ]

14. L. Bastos and A. O’Hagan, “Diagnostics for Gaussian process emulators,” Technometrics, vol. 51, no. 4, pp. 425-438, 2009. [ Links ]

15. A. O’Hagan , “Bayesian analysis of computer code outputs: A tutorial,” Reliability Engineering & System Safety, vol. 91, no. 10-11, pp. 1290-1300, 2006. [ Links ]

16. I. De La Pava et al., “A Gaussian Process Emulator for Estimating the Volume of Tissue Activated During Deep Brain Stimulation,” in 7th Iberian Conference on Pattern Recognition and Image Analysis (IbPRIA), Santiago de Compostela, Spain, 2015, pp. 691-699. [ Links ]

17. T. Zhang and W. Grill, “Modeling deep brain stimulation: Point source approximation versus realistic representation of the electrode,” Journal of neural engineering , vol. 7, no. 6, 066009 pp. 1-11, 2010. [ Links ]

18. L. Chang, D. Jones, and C. Pierpaoli, “Restore: Robust estimation of tensors by outlier rejection,” Magnetic Resonance in Medicine, vol. 53, no. 5, pp. 1088-1095, 2005. [ Links ]

19. D. Tuch, V. Wedeen, A. Dale, J. George and J. Belliveau, “Conductivity tensor mapping of the human brain using diffusion tensor MRI,” Proceedings of the National Academy of Sciences of the USA, vol. 98, no. 20, pp. 11697-11701, 2001. [ Links ]

20. A. Chaturvedi , C. Butson , S. Lempka, S. Cooper , andC. 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-67, 2010. [ Links ]

21. F. Rattay, “Analysis of models for extracellular fiber stimulation,” IEEE Transactions on Biomedical Engineering, vol. 36, no. 7, pp. 676-682, 1989. [ Links ]

22. R. Plonsey and R. Barr, Bioelectricity: A quantitative approach, 3 rd ed. New York, USA: Springer Science & Business Media, 2007. [ Links ]

23. D. Sterratt, B. Graham, A. Gillies, and D. Willshaw, Principles of computational modelling in neuroscience, 1 st ed. Cambridge, UK: Cambridge University Press, 2011. [ Links ]

24. C. McIntyre andW. Grill , “Excitation of central nervous system neurons by nonuniform electric fields,” Biophysical journal, vol. 76, no. 2, pp. 878-888, 1999. [ Links ]

25. D. McNeal, “Analysis of a model for excitation of myelinated nerve,” IEEE Transactions on Biomedical Engineering , vol. 23, no. 4, pp. 329-337, 1976. [ Links ]

26. F. Rattay, “Analysis of models for external stimulation of axons,” IEEE Transactions on Biomedical Engineering , vol. 33, no. 10, pp. 974-977, 1986. [ Links ]

27. C. McIntyre , A. Richardson, andW. 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 ]

28. G. Walckiers, “Bio-electromagnetic model of deep brain stimulation,” Ph.D. dissertation, École polytechnique fédérale de Lausanne, Lausanne, Switzerland, 2009. [ Links ]

29. M. Hines, A. Davison, and E. Muller, “Neuron and python,” Frontiers in Neuroinformatics, vol. 3, no. 1, pp. 1-12, 2009. [ Links ]

30. M. Hines and N. Carnevale, “The neuron simulation environment,” Neural Computation, vol. 9, no. 6, pp. 1179-1209, 1997. [ Links ]

31. C. Bishop, Pattern recognition and machine learning, 1 st ed. New York, USA: Springer, 2006. [ Links ]

32. C. Rasmussen and C. Williams, Gaussian Processes for Machine Learning, 1 st ed. Cambridge, USA: MIT Press, 2006. [ Links ]

33. M. Kuss and C. Rasmussen, “Assessing approximate inference for binary gaussian process classification,” The Journal of Machine Learning Research, vol. 6, pp. 1679-1704, 2005. [ Links ]

34. H. Nickisch and C. Rasmussen, “Approximations for binary gaussian process classification,” Journal of Machine Learning Research, vol. 9, pp. 2035-2078, 2008. [ Links ]

35. 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. 14, pp. 284-289, 2006. [ Links ]

36. M. Hines andN. Carnevale , “Translating network models to parallel hardware in neuron,” Journal of Neuroscience Methods , vol. 169, no. 2, pp. 425-455, 2008. [ Links ]

Received: September 14, 2016; Accepted: August 07, 2017

* Corresponding 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