The goals of this study were to develop quantitative techniques to predict the VTA as a function of arbitrary stimulation parameters and electrode geometry, and to use those techniques to customize the electrode design to a target nucleus. Converging theoretical (

McIntyre *et al* 2004b) and experimental (

Hashimoto *et al* 2003) results suggest that DBS generates an excitatory effect on axons surrounding the electrode. While correlations between axonal activation and the therapeutic mechanisms of DBS remain controversial, one leading hypothesis is that high frequency stimulation results in an override of the underlying pathological neural activity patterns (

Montgomery and Baker 2000,

Hashimoto *et al* 2003,

Grill *et al* 2004,

McIntyre *et al* 2004a). Therefore, we concentrated our analysis on the characterization of axonal activation with DBS.

Neural stimulation was estimated with an integrated model that combined FEM electric field solutions with multi-compartment cable models of myelinated axons. The electric field generated during monopolar stimulation by the DBS electrode was calculated from the Poisson equation with a Fourier FEM solver to determine time- and space-dependent voltage within the tissue medium (

Butson and McIntyre 2005). The voltage waveforms were subsequently interpolated onto cable model axons distributed around the electrode, and threshold values of the stimulation voltage necessary for action potential generation were calculated. The second spatial difference of the extracellular potential distribution (Δ

^{2}*V*_{e}/Δ

*x*^{2}) at the site of action potential initiation was determined at the stimulation threshold for each axon (Δ

*x* = internodal spacing of 0.5 mm). The Δ

^{2}*V*_{e}/Δ

*x*^{2} solutions were used to create a VTA prediction scheme as a function of stimulation parameters.

Finite element model

Axisymmetric FEM models of DBS electrodes with approximately 13 000 nodes were constructed in FEMLAB v3.1 (Comsol Inc, Burlington, MA) () as previously described (

Butson and McIntyre 2005). The tissue medium was modeled as homogeneous and isotropic with conductivity

*σ* = 0.3 S m

^{−1}. The axisymmetric volume conductor measured 100 mm tall by 50 mm wide. Electrode contact dimensions ranged from 0.25 mm to 2.5 mm in diameter and from 0.5 mm to 3.8 mm in height. These dimensions were centered around the Medtronic 3387/3389 quadripolar DBS electrode contact dimensions (1.5 mm height, 1.27 mm diameter) (Medtronic Inc, Minneapolis, MN). The electrode lead was modeled as an electrical insulator with the exception of the contact area, which was used for voltage-controlled stimulation. Stimulation voltage was specified at the contact surface, and the Poisson equation was solved for the voltage within the tissue medium using the Fourier FEM solver (

Butson and McIntyre 2005). While the Poisson equation alone provides a spatial voltage solution, it does not account for the time dependence of the stimulus waveform or the capacitance of the electrode–tissue interface. The Fourier FEM solver overcomes these limitations by using a complex stiffness matrix that represents the capacitance of the electrode–tissue interface, and by solving the Poisson equation at multiple frequencies to reconstruct the time-dependent stimulus waveform. The resulting solutions are both time- and space-dependent, and incorporate the effects of electrode capacitance on the stimulus waveform delivered to the tissue during voltage-controlled stimulation (

Butson and McIntyre 2005).

Neural stimulation prediction

Field-axon simulations were conducted using Fourier FEM DBS electrode models coupled to 5.7

*μ*m diameter myelinated axon models (

McIntyre *et al* 2002,

Butson and McIntyre 2005). A collection of 119 model axons were distributed in a 17 × 7 matrix oriented perpendicular to the electrode shaft (). This orientation of axons was used to identify the spatial extent of activation in the vertical and horizontal directions relative to the electrode shaft (localization of activation in axons oriented parallel to the shaft would be ambiguous in the vertical direction). The axons were placed 1 mm to 4 mm lateral to the electrode and +4 mm above to −4 mm below the center of the electrode contact. Each model axon included 21 nodes of Ranvier with 0.5 mm internodal spacing.

The time-dependent potential distribution generated in the tissue medium from the Fourier FEM solution was interpolated onto the length of each axon, and the time-dependent transmembrane potential variations induced by the stimulation were calculated in NEURON v5.7 (

Hines and Carnevale 1997). Threshold stimulus amplitudes were defined that generated action potentials in a one-to-one ratio with a stimulus frequency of 130 Hz, representative of typical clinical DBS parameter settings (

Volkmann *et al* 2002). The threshold values were used to create 2D contours to define the boundary of activation as a function of the stimulus amplitude. These contours were swept around the axis of the electrode to determine the VTA.

An alternative approach to the computationally intensive field-neuron simulations described above is the use of an activating function-based technique (

Rattay 1986). The second difference of the extracellular potential distribution along a neural process (Δ

^{2}*V*_{e}/Δ

*x*^{2}) provides a quantitative estimate of the polarization of the neuron in response to an applied electric field. One attractive feature of this method is that for a given stimulus pulse width, Δ

^{2}*V*_{e}/Δ

*x*^{2} threshold values are relatively constant across a wide range of electrode designs (see the results section, ). However, Δ

^{2}*V*_{e}/Δ

*x*^{2} alone is not an accurate predictor of neural activation (

Warman *et al* 1992,

Zierhofer 2001,

Moffitt *et al* 2004). As a first approximation, a Δ

^{2}*V*_{e}/Δ

*x*^{2} threshold value can be selected and the spread of stimulation can be defined by the volume of tissue encompassed by that criteria. However, this overestimates activation near the electrode and underestimates activation distant to the electrode. An alternative is to create a Δ

^{2}*V*_{e}/Δ

*x*^{2} threshold curve fit as a function of electrode to axon distance (). This technique has the advantage of providing accurate threshold values across all electrode–axon distances, but a separate curve must be generated for each pulse width. Further, this method fails to accommodate different electrode designs and breaks down for electrodes with a small diameter/height aspect ratio ().

A third approach that addresses the limitations mentioned above is to determine Δ^{2}*V*_{e}/Δ*x*^{2} threshold values as a function of pulse width and voltage. Specifically, Δ^{2}*V*_{e}/Δ*x*^{2} threshold values are recorded, and these values are expressed as a function of cathodic voltage (V) times pulse width (PW, *μ*s) (). This expression allows two stimulation parameters to be condensed into a single number for prediction of thresholds. Further, threshold values recorded this way were found to be valid for a wide range of electrode designs and stimulation parameters. These values can then be used to create 2D spatial contours that are swept around the *z*-axis to define the VTA volume (). For purposes of volume calculations, it is often convenient to describe the VTA contours with analytical functions. To do so, each contour is described by an ellipse:

where *x*_{0}, *y*_{0} is the center of the ellipse, and *a* and *b* are the semimajor and semiminor axes, respectively (assuming *b < a*). The semimajor and semiminor coefficients are calculated from the following: *a* = distance of threshold value from electrode contact along *x*-axis; *b* = maximum *y* value of 2D threshold contour. Under the model conditions used in this study, the electrode contact is centered on the origin and the center of each ellipse is *x*_{0} = *a*, *y*_{0} = 0. With this method, Δ^{2}*V*_{e}/Δ*x*^{2} threshold values and VTA volumes can be predicted for a wide range of electrode designs and stimulation parameters.

In this study, we solved the complete field-axon model for each electrode design shown in to build the Δ

^{2}*V*_{e}/Δ

*x*^{2} threshold relationships. Our motivation for this exercise was twofold. First, we wanted to provide analytical functions that could be used by other DBS investigators, who may not have the expertise to develop full field-axon models, thereby allowing relatively accurate predictions for their own studies. Second, we believe that generalized Δ

^{2}*V*_{e}/Δ

*x*^{2} threshold relationships, developed under idealized conditions, can be useful for rapid evaluation of a given electrode design and/or stimulation parameter settings in a more complicated environment. For example, we have recently developed DBS models derived from human diffusion tensor magnetic resonance imaging data that can be customized to individual patients (

Butson *et al* 2004,

2005b,

McIntyre *et al* 2004c). However, the prospect of running full field-axon simulations to make VTA predictions for pre-operative electrode targeting or post-operative stimulation parameter selection is highly unrealistic given the thousands of permutations that may be necessary to find an optimal electrode location and/or stimulation parameter setting. Our Δ

^{2}*V*_{e}/Δ

*x*^{2} threshold relationships provide a computationally efficient VTA prediction technique to apply DBS modeling to clinical research. Further, we have used our Δ

^{2}*V*_{e}/Δ

*x*^{2} threshold relationships in this study to evaluate the VTA generated by electrode contacts designs with heights and diameters that fall between electrode geometries for which we have full field-axon data ( and ). In turn, we are able to more quickly and efficiently search the parameter space for electrode design optimization applications.