Optimized programming algorithm for cylindrical and directional deep brain stimulation electrodes

Objective. Deep brain stimulation (DBS) is a growing treatment option for movement and psychiatric disorders. As DBS technology moves toward directional leads with increased numbers of smaller electrode contacts, trial-and-error methods of manual DBS programming are becoming too time-consuming for clinical feasibility. We propose an algorithm to automate DBS programming in near real-time for a wide range of DBS lead designs. Approach. Magnetic resonance imaging and diffusion tensor imaging are used to build finite element models that include anisotropic conductivity. The algorithm maximizes activation of target tissue and utilizes the Hessian matrix of the electric potential to approximate activation of neurons in all directions. We demonstrate our algorithm’s ability in an example programming case that targets the subthalamic nucleus (STN) for the treatment of Parkinson’s disease for three lead designs: the Medtronic 3389 (four cylindrical contacts), the direct STNAcute (two cylindrical contacts, six directional contacts), and the Medtronic-Sapiens lead (40 directional contacts). Main results. The optimization algorithm returns patient-specific contact configurations in near real-time—less than 10 s for even the most complex leads. When the lead was placed centrally in the target STN, the directional leads were able to activate over 50% of the region, whereas the Medtronic 3389 could activate only 40%. When the lead was placed 2 mm lateral to the target, the directional leads performed as well as they did in the central position, but the Medtronic 3389 activated only 2.9% of the STN. Significance. This DBS programming algorithm can be applied to cylindrical electrodes as well as novel directional leads that are too complex with modern technology to be manually programmed. This algorithm may reduce clinical programming time and encourage the use of directional leads, since they activate a larger volume of the target area than cylindrical electrodes in central and off-target lead placements.


Introduction
Deep brain stimulation (DBS) has become widely accepted as a surgical intervention for movement disorders such as Parkinson's disease (Rodriguez-Oroz et al 2001, essential tremor (Benabid et al 1996) and dystonia (Vidailhet et al 2005(Vidailhet et al , 2007, and is currently being explored for a growing number of psychiatric disorders, such as Tourette Syndrome (Porta et al 2009), obsessive-compulsive disorder (OCD) (Greenberg et al 2010) and treatment-resistant depression (Mayberg et al 2005). DBS is now being used for a relatively large number of applications, but clinical DBS technology remains largely unchanged from its initial design, which was created nearly four decades ago (Brice and Mclellan 1980). In recent years, numerous academic and industry groups have explored more complex lead designs with greater numbers of contacts and increased field shaping ability (figure 1) (Buhlmann et al 2011, Vitek and Starr 2013, Contarino et al 2014, Pollo et al 2014, Boston Scientific Corportation 2015, van Dijkand et al 2015, Willsie and Dorval 2015a, 2015b, St. Jude Medical 2015. As the DBS lead designs are trending toward increased complexity, there is an associated need to update manually-based DBS programming methods to accommodate new technology. Modern DBS programming consists largely of manual, trial-and-error methods in which clinicians adjust settings based on patient responses (Rizzone et al 2001, Moro et al 2002, Volkmann et al 2002. Manual programming poses a substantial time burden on the clinician while potentially subjecting the patient to discomfort, and is unlikely to yield optimized parameter settings. Implants for the overwhelming majority of DBS patients involve leads with four evenly-spaced cylindrical contacts, which, despite being the simplest of lead designs, require approximately 18-32 h of clinical programming in the first year post-surgery (Hunka et al 2005). Programming time varies depending on the target and the programming center, and the number of adjustment visits that a patient must make to the programming center can range widely, from four to 17 visits in the first year (Ondo and Bronte-Stewart 2005). There are over 25 000 possible combinations of programming parameters such as pulse width, frequency, and voltage, and 65 combinations of contact selection for a standard four-contact lead such as the Medtronic 3389 (Kuncel and Grill 2004). The number of theoretically possible configurations for a lead is exponential in the number of contacts, so that with the emergence of increasingly complex lead designs, manual programming will soon be infeasible as the number of possible parameter combinations increases by many orders of magnitude.
A number of factors contribute to positive outcomes for DBS, including careful patient selection, accurate lead placement in the target area, and the effective programming of DBS (Walter and Vitek 2004, Machado et al 2006, Volkmann et al 2006. However, once a patient has been given implants, the only aspect that can be modified to improve clinical outcomes is the programming of the DBS stimulator. Programming becomes especially important if the DBS lead does not precisely hit its target because stimulation that spreads away from the intended target region may lead to psychiatric and motor side effects (Volkmann et al 2002, Krack et al 2003, Appleby et al 2007, Rossi et al 2016. Even with state-of-the-art surgical techniques, it is typical to see a deviation from the target area of 2 mm in lead placement (Mobin et al 1999, Patel et al 2002, Burchiel et al 2013. Additionally, the brain itself can shift by 2 mm to 4 mm during surgery, especially after the first lead in bilateral DBS surgery has been implanted (Winkler et al 2005, Khan et al 2008, Hunsche et al 2009, Walter et al 2009, Pallavaram et al 2010. In the case of DBS for Parkinson's disease (PD) patients, for both subthalamic nucleus (STN) and globus pallidus pars interna (GPi) targets, if the distance from the target is greater than 3 mm, the electrode does not reach any portion of the area physiologically defined as optimal. Unfortunately, such placement errors have been reported to occur in over 40% of implantations (Guridi et al 2000, Okun et al 2005, Ellis et al 2008, Rolston et al 2016.
Cylindrical contacts produce nearly spherical regions of axonal stimulation, which are distorted based on the level of anisotropy in the surrounding tissue. Regardless, such cylindrical electrodes are ill-equipped to steer stimulation back into a target structure that was missed during surgical implantation. The prevalence of lead placement deviation has motivated the creation of directional DBS electrodes. With a greater number of smaller contacts, directional leads can generate asymmetric stimulation profiles and steer stimulation back toward a missed target. However, the increasing complexity of directional leads compounds the clinical difficulty and time burden associated with their programming.
The objective of this paper is to create an automated method for optimal parameter selection that can compute patient-specific contact settings and configurations in near real-time. The algorithm has the ability to target and avoid nuclei and/or fiber targets while limiting charge density to safe levels (<30 µC cm −2 ) as mandated by the U.S. Food and Drug Administration (FDA) (McCreery et al 1990, Shannon 1992. The algorithm can be used for a variety of lead designs and computes contact configurations using pre-computed finite element method (FEM) solutions of the electric potential distribution, taking into account a patient's lead position, as well as the anisotropic tissue conductivities measured from diffusion tensor imaging (DTI). Quick programming results are possible through use of the activating function approximation, based on the second spatial difference of voltage along an axon, to estimate neuronal activation (McNeal 1976, Rattay 1986, Warman et al 1992, McIntyre et al 2004. To extend the activating function approximation to 3D space, we used the Hessian of the electric potential: a matrix of all secondorder partial derivatives with respect to space. The methods developed in this paper can be used to target either (i) nuclear regions in the brain or (ii) fiber tracts. For both objectives, we develop a family of constraints that, when enforced, simultaneously avoid other regions or fiber tracts that might not be desirable to stimulate, based on disease state.
As DBS programming remains a manual process, sifting through the parameter space is already a difficult and timeconsuming task for leads with cylindrical electrodes and even more complex for directional leads. Our approach to DBS parameter optimization can accommodate new directional leads and compute contact configurations for these complex lead designs in near real-time, thus lowering the barrier of programming complexity for novel directional leads. A major goal of our automated, patient-specific DBS programming algorithm is to reduce programming time for all lead designs and to individualize DBS therapy.

Materials and methods
We tailored the optimization algorithm to an example patient using MRI to locate relevant neural targets and DTI to generate anisotropic tissue conductivities and tractography. In a typical DBS implant patient, MRI data and postoperative CT are used to determine lead location, but for this study, we used the atlas brain created from patient data as defined in Wakana et al (2004) as a model patient to demonstrate the capabilities of the optimization algorithm.

Finite element method (FEM)
We used the finite element method (FEM) in SCIRun 4.7 (SCI Institute, University of Utah, Salt Lake City, UT) to solve the bioelectric field problem for three lead designs: the Medtronic 3389, direct STNAcute, and the Medtronic-Sapiens, which have 4, 8 and 40 contacts, respectively (figure 1). We chose these three lead geometries to explore algorithm performance because the Medtronic 3389 is commonly used clinically, and the direct STNAcute and Medtronic-Sapiens designs have been implanted in human patients (Contarino et al 2014, Pollo et al 2014. For each lead design, a tetrahedral mesh was generated inside a 100 × 100 × 100 mm cube, with the lead in the middle of the cube. The mesh was generated in SCIRun using the interface module to TetGen (Si 2015), in which the surfaces of the leads and node locations with 0.1 mm resolution on a 20 × 20 × 20 mm regular grid were pre-defined. The resulting tetrahedral FEM mesh consisted of ~ 9 million nodes and ~52 million elements. At a resolution of 52 million elements, we saw convergence of both the solution of the electric potential and the related second difference. When we compared the electric potential linearly solved over the mesh of 52 million elements to a cubically solved mesh of 2.5 million elements, the electric potential had a relative error of 0.01%. Most of the nodes were distributed densely around the electrode boundaries. Isotropic conductivities were used for each lead geometry, with electrode contacts set to σ = 1 × 10 6 S m −1 and the shaft set to σ = 1 × 10 −10 S m −1 (Wei and Grill 2005, Miocinovic et al 2006, Zhang and Grill 2010. Anisotropic conductivities for the surrounding tissue were derived from DTI. We solved the Poisson equation (equation (1)) to calculate electric potential (V e ) at every node in the 3D space under Dirichlet boundary conditions given conductivity tensors (σ) and a single source at the center of each contact (i) (Butson and McIntyre 2005, Butson et al 2007a.
Here, Ω C is the volume of the FEM mesh and Γ Neu is its outer boundary. For voltage-controlled stimulation, V e,0 is the stimulation voltage, Γ Dir = { x 0 }, i.e. the position of the voltage source inside the active contact, and i = 0. For current-controlled stimulation, Γ Dir is an empty set and i = I 0 δ x0 is a point current source of strength I 0 at position x 0 . We solved the forward problem for each individual contact at −1 V to form a set of bioelectric field solutions for each lead. We then took advantage of the system's linearity to approximate all possible combinations of active contacts through the principle of superposition.

Conductivity tensors
We used DTI to improve the accuracy of our FEM model by incorporating anisotropic tissue conductivities. The use of anisotropic conductivities contributes to the patient specificity of the model. As a result, we obtain a nonsymmetric voltage distribution. Additionally, ignoring heterogeneity can reduce a model's predictive ability of neural activation (Hyde et al 2012, Howell andMcIntyre 2016). We used the volume-normalized approach to process the tensor data following equation (2) The volume of each anisotropic conductivity tensor, D, was derived from DTI data and was scaled to the equivalent volume of an isotropic conductivity tensor for the isotropic conductivity, d iso , to obtain anisotropic conductivity tensors, σ aniso . Therefore, D was scaled by the ratio of the d iso to the cubic root of the product of the eigenvalues, d k . The chosen isotropic conductivity was 0.2 S m −1 , a value commonly used for brain tissue conductivity (Ranck 1963, Wei and Grill 2005, Zhang and Grill 2010. T1-weighted MRI sagittal and axial slices shown with the segmented STN (green) and thalamus (yellow). (B) A tensor field generated from DTI MRI demonstrates anisotropic flow through brain tissue; yellow coloring indicates high anisotropy while black coloring indicates isotropy. (C) The FEM model incorporates heterogeneous tissue conductivities from DTI MRI to capture anisotropic voltage spread through neighboring brain tissue. (D) Tracts from the IC, containing cortico-motor axons, were generated from a seed region placed laterally from the STN. A total of 708 IC trajectories were generated from diffusion tensors shown in panel B and had a minimum fractional anisotropy value of 0.15. (E) Multi-compartment cable models were solved for each IC trajectory to quantify activation by optimized settings produced by the algorithm. The stimulus waveform consists of a 90 µs cathodic pulse, followed by a 100 µs delay, then a 900 µs anode pulse at 10% of the cathodic voltage. An example spike train can be seen for 130 Hz stimulation at −1 V.

Patient imaging data
2.3.1. Segmentation. To apply the optimization algorithm on a patient-specific level, the target area must be precisely defined for that patient. The target area was segmented from the patient atlas MRI (figure 2) (Wakana et al 2004). The chosen target for this study was the STN, a common target for DBS treatment of Parkinson's disease (Limousin et al 1998).
Using the STN as the neural target makes for an ideal example since it shares a boundary with the internal capsule (IC), a bundle of white matter tracts that are associated with negative side effects if stimulated (Ashby et al 1999, Tamma et al 2002, Krack et al 2003.

Tractography.
We generated deterministic tractography from DTI tensor data taken from Wakana et al (2004) and used a deterministic streamline algorithm in 3D Slicer (Fedorov et al 2012) to determine fiber trajectories of the IC. The 708 trajectories were generated from a seed region of size 5 mm × 2 mm × 1 mm lateral from the STN, and diffusion tensors had a fractional anisotropy value of 0.15 or greater (Chaturvedi et al 2010, Makris et al 2016). Additionally, we used a stopping track curvature of 0.4, an integration step of 0.5 mm, and restricted the minimum and maximum path lengths to 15 mm and 400 mm, respectively. We would also like to point out that our representation of the IC is only a small portion of the whole IC, which is a large white matter structure in which critical sensory and motor pathways travel to and from the cerebral cortex. However, generating the solutions of such cable models is time-consuming and does not align with our goal of creating a real-time optimization algorithm. Conveniently, extracellular electrical stimulation can be approximated by the second spatial derivative of the electric potential along the nodes of Ranvier of an axon (Rattay 1986(Rattay , 1999. The activating function approximates the depolarizing influence of stimulation along an axon, and neural activation can be estimated from the second difference threshold corresponding to the start of axon firing. Using the second difference threshold to approximate activation may not fully capture the behavior of neurons, and could result in either overestimation or underestimation of neuron activation, but we attempted to limit the error in our estimates of neural activation by verifying the activating function threshold with cable models. We determined the activating function (A.F.) by calculating the second spatial difference of extracellular potentials, V e,n , at the nth node of Ranvier. The second difference was taken across neighboring nodes of Ranvier, at n − 1 and n + 1, with an internodal spacing of Δx (equation (3)). We used a scaled version of the activating function (S.A.F.) that removes the denominator, Δx 2 , to generalize neural activation in our model of optimization (equation (4)) (McIntyre et al 2004, Martens et al 2011. The activating function is able to approximate activation only along a 1D fiber. To approximate neural activation in a 3D space, we use the Hessian matrix (equation (5)) of second spatial partial derivatives of the electric potential. The second derivative of the function V e at location x in the direction of the unit vector u(x) can be written using the inner product, as in equation (6).
The Hessian of the electric potential data was computed on a 20 × 20 × 20 mm cube around the lead with an 0.1 mm resolution. We included the grid points on which we calculated the Hessian as nodes in the 3D mesh, so that the electric potential could be solved directly at those nodes in the FEM forward solution. While the Hessian was calculated at all node points on the grid, the optimization algorithm excluded Hessian values on nodes that were located within the lead body. The use of the Hessian matrix enables us to generalize the activating function to predict activation along not only fiber tracts but also boundaries of nuclei-a novel application of the activating function that we discuss further in section 2.6.1. We also note that because the Hessian is a linear derivative of electric potential, the principle of superposition can be applied to the Hessian as well, to approximate neural activation for all combinations of contact configurations and settings.

Verifying neural activation using NEURON.
To verify our approximation of neural activation using the activating function, we built multicompartment cable models for each IC trajectory in NEURON 7.4. The IC fibers were modeled as 10 µm diameter, myelinated axons with nodes of Ranvier, paranodal, and intermodal sections explicitly defined (McIntyre et al 2002, Chaturvedi et al 2010. Electric potentials calculated from the FEM model were interpolated onto the multicompartment cable models, to determine whether algorithm-determined DBS settings resulted in action potential generation in the IC fibers. Transmembrane currents are implicitly calculated by NEURON as part of the extracellular mechanism. In the cable models, we set frequency and pulse width to 130 Hz and 90 µs to simulate typical DBS settings. The biphasic pulse consisted of a 90 µs cathodic pulse followed by a 100 µs interpulse interval and a charge-balancing 900 µs anodic pulse at 0.1 of the cathodic amplitude (figure 2(E)). In our study we used a 90 µs pulse width, as an example of a typically used setting and for a more conservative approximation of neural activation.

Optimization problem
Our algorithm has the ability to accommodate programming guidelines that include any combination of the following: (1) target nuclei; (2) avoid nuclei; (3) target fiber tracts; (4) avoid fiber tracts, and; (5) limit charge density. The patient-specific model used in this study to simulate programming is an example DBS patient in which we simulate Parkinson's disease, with the nucleus of interest being the STN. Additionally, we include the thalamus as a nucleus of avoidance and the IC as a fiber tract of avoidance. In our example patient STN DBS scenario, we do not have a fiber tract of interest as it is difficult to generate tractography data for basal ganglia structures with the DTI resolutions that can be achieved with a 1.5 T scanner. This might be possible in the future, with higher resolution scans.
For safety purposes, we determined the maximum allowable voltages that would keep charge density below 30 µC cm −2 for each electrode for the frequency of 130 Hz and pulse width of 90 µs (McCreery et al 1990, Shannon 1992. To establish the maximum allowable voltages, we either drew clinical impedance values taken from the literature if available, or calculated the impedance values from the FEM model. To keep the charge density below the maximum 30 µC cm −2 charge per phase, we used a conservative impedance value of 500 Ω for the Medtronic 3389, calculated the direct STNAcute to have an impedance value of 1.5 kΩ, and used the reported impedance value of 2.5 kΩ for the Medtronic-Sapiens (Butson et al 2006, Contarino et al 2014, Pollo et al 2014. It is important to note that these impedance values are at the lower end of clinical impedance values (Satzer et al 2014), so providing more conservative estimates for optimal contact amplitudes and configurations, since a lower impedance value allows a greater spread of the electric potential into the tissue.
We placed leads so that the equivalent center of the most distal contact of the Medtronic 3389 was located at the centroid of the STN. An alternate position, in which the lead was displaced 2 mm in the lateral direction, was chosen as a typical, yet challenging, lead placement in which the lead is positioned problematically close to the IC, a fiber tract of avoidance in our patient model.

Optimization algorithm
In DBS, the goal is to increase therapeutic benefit by stimulating a target region while limiting stimulation spread to other areas that might be responsible for side effects. Such a scenario fits problems that can be solved through linear convex optimization by maximizing an objective function, neural activation in our case, subject to constraints such as limiting stimulation of regions inducing side effects. We use the interior point method of linear convex optimization to determine ideal contact configurations and contact amplitudes based on neural targeting criteria. We use CVX, a package for MATLAB for the defining and solving of such convex optimization problems (Grant and Boyd 2008, 2014).
2.6.1. Definitions. We optimize the selection of active contacts in each lead for n individually controllable contacts. The algorithm has the ability to target a nucleus, a user-defined volume, or fiber tracts by maximizing the activating function over the region. Similarly, activation in an area can be avoided by keeping the second difference value along the boundary of a region below a threshold, α. The threshold, α, which we refer to as the sensitivity parameter, represents the scaled activating function in units of millivolts (equation (4)). We also keep the charge density within the safety limits for each electrode. It is important to impose charge density constraints since smaller electrode designs generally have higher contact impedance values, and certain combinations of settings can result in charge density values that exceed 30 µC cm −2 on the surface of a contact.
To either target or avoid neural regions, we use a set of equations to manipulate the Hessian matrix, H i , to generate voltage values, c i , for each contact, i, out of n individually controllable contacts. To target a nucleus or a user-defined region, Ω, we maximize the targeting objective function d Ω (equation (7)), the sum of the eigenvalues of H i averaged over region Ω. This objective function quantifies activation in all principal directions to accommodate multiple possible fiber orientations within nuclei.

Nucleus/Region Targeting
In order to limit activation of a nucleus or user-defined region, we limit the maximum eigenvalue as defined in the avoidance function, e Ω (equation (8)), below a threshold at each point along the boundary of the nucleus or region. The maximum eigenvalue corresponds to the maximal second spatial derivative in the region of interest, so that limiting this value below the sensitivity parameter, α, corresponds to a reduction in firing throughout the entire region. In order to fully suppress firing in the region, the parameter α must be chosen to correspond to the threshold of the first firing of neurons representative of that region. Selection of the sensitivity parameter is explored in section 3.1.

Nucleus/Region Avoidance
(8) Finally, we interpret that the functions f v,Ω and g v,Ω measure the second difference along a fiber tract in the direction of vector v(x). The vectors, v, are derived from tractography data and must have a magnitude that corresponds to the distance between the flanking nodes of each node of Ranvier for the fiber tract of interest. For example, the distance between flanking nodes for a 10 µm fiber is 2.3 mm (McIntyre et al 2002). In order to target a fiber tract, one must maximize the objective function f v,Ω (equation (9)), and to avoid a given fiber tract, one must limit the avoidance function g v,Ω (equation (10)) below the sensitivity parameter. It's important to note that in equations (9) and (10), x is a node point on the 20 × 20 × 20 mm grid in the fiber tract region, Ω. The fiber tract region is represented as a vector field, in which directionality vectors from the fiber tracts were interpolated onto the grid nodes.

Fiber Targeting
Fiber Avoidance (10) If firing occurs at one point along an axon, the entire axon is seen as activated because action potentials will rapidly spread from the point of action potential generation. In order to properly constrain such activation in a nucleus or along a fiber tract (equations (8) and (10)), the algorithm must check that the second difference is below the sensitivity parameter α at every point x in the region of interest, Ω. If there is no charge present inside the region of interest, it is sufficient to check only around the outer boundary, to reduce computation time. Node-by-node checking for activation increases computation time, but is necessary to avoid regions that may be responsible for side effects when stimulated.

STN DBS optimization paradigm.
The algorithm follows the paradigm in which targeting involves maximization of either equation (7) or (9) and region avoidance involves keeping the second difference values in equations (8) and (10) below a threshold. The optimization algorithm is able to consider any combination of the above equations, including multiple instances of the same equations for different regions or fiber tracts. To demonstrate a potential programming case for the optimization model, we used the patient-specific example of STN DBS in which we employed equation (7) for STN activation, equation (8) for thalamus avoidance, equation (10) for IC avoidance, and a 30 µC cm −2 charge density constraint, using the following set of equations: where Ω STN is the physical region the STN occupies for every x in the IC region, Ω IC .
(11) We chose these targeting and avoidance constraints to demonstrate the capabilities of the optimization algorithm; we do not make claims as to which regions are optimal or suboptimal for therapeutic stimulation. More studies will have to be conducted to explore optimal stimulation targets, and we hope that our algorithm can eventually be used as a tool in such studies, to precisely stimulate or avoid neural targets.

Results
We ran our algorithm for three lead designs: the Medtronic 3389, the direct STNAcute, and the Medtronic-Sapiens (figure 1). The atlas brain (Wakana et al 2004) serves as the patient-specific model in which we targeted the STN, avoided the thalamus and IC, and limited the charge density as an example programming application of the algorithm. We assessed algorithm performance in two DBS lead positions, one in which the DBS lead was centrally placed in the STN, and another in which the lead was placed 2 mm laterally (figure 3(A)). We included the 2 mm lateral placement case to simulate a challenging lead placement adjacent to the IC. The distance of 2 mm was chosen to represent typical lead positioning deviations that might result in well-characterized side effects through stimulation of the IC (Ashby et al 1999, Tamma et al 2002, Krack et al 2003.

Sensitivity threshold parameter
When not considering the many possible combinations of targeting and avoidance criteria for nuclei and fiber tracts, the optimization algorithm has two free numerical parameters: the maximum programmable contact voltage and the sensitivity parameter, α. The maximum allowable contact voltage is determined by voltage values that meet the criterion that charge density per phase should be below 30 µC cm −2 . For this study, those values were calculated to be 10 V for the Medtronic 3389, 5 V for the directional contacts of the direct STNAcute, and 3.4 V for the Medtronic-Sapiens. The maximum allowable voltage decreases in the directional leads because smaller-sized contacts have larger impedance values and, as a result, voltages must be reduced appropriately to stay within the charge density constraints.
The second free parameter in the algorithm is the sensitivity α, which represents a user-input value corresponding to the first firing of neuron fibers. This value is analogous to the activating function threshold value, which can vary based on fiber diameter, distance to fiber, and the orientation of the fiber paths (McIntyre et al 2002). A variety of threshold values could apply to a patient, based on the neurons of interest, but based on values in the literature we estimate this threshold to be within a range of 5-30 mV (Rattay 1986, McIntyre et al 2004, Martens et al 2011. For the STN DBS case at hand, we determined a reasonable sensitivity parameter based on the second difference value, which corresponded with the initiation of IC firing using the Medtronic 3389 for both the STN center and STN lateral positions ( figure 3(B)). When the second difference was greater than 20 mV, the IC capsule fibers began to fire. The second difference profile with the directional electrodes is more complex since the field can be highly asymmetric along a fiber tract. The impact of complex contact configurations on the second difference has not been well studied for directional electrodes.
We explored how the sensitivity parameter altered the output configuration of the algorithm for the Medtronic-Sapiens lead ( figure 3(D)). As the sensitivity parameter was lowered, either fewer contacts were turned on or the same contacts were turned on at lower amplitudes. As the sensitivity input parameter was adjusted, the size of the stimulation field changed, but the directional shaping of the field was conserved. It can be seen that in the centrally placed case, contacts facing all directions on the Medtronic-Sapiens were utilized. However, in the lateral case, the contacts (colored purple, blue, and green) that face the STN were preferentially activated over the orange contacts that face the IC.

Optimal configurations for cylindrical and directional electrodes
We ran the optimization algorithm for the STN DBS patient model in the center and lateral positions. For this patient DBS model, we maximized activation of the STN by maximizing equation (7), which corresponds to the sum of the eigenvalues within the STN region. Due to the convexity of the optimization problem, results from each optimization run represent unique contact configurations and contact amplitudes based on the target and avoidance criteria chosen. The resulting contact configurations and amplitudes are consistent across every re-run of the optimization algorithm.
In our optimization paradigm, maximization of STN activation is constrained by keeping the activating function on the IC and the thalamus below the sensitivity threshold, α. In the centrally placed DBS lead condition, figure 4(A) shows the optimal contact configurations and amplitudes in 2D illustrations. For the center placement, since contacts 0 and 1 of the Medtronic 3389 are located inside the STN, the optim ization algorithm determined they are most effective in stimulating the STN and set contacts 0 and 1 to approximately −1 V each. In the directional leads, the lateral contacts that face the IC are minimally used. For the 2 mm lateral placement, very little of the STN is activated with the Medtronic 3389; all four contacts turn on at minimal voltages since they cannot be programmed at greater amplitudes without activating the IC. Considering the lateral position for both the direct STNAcute and the Medtronic-Sapiens, usage of the lateral contacts is avoided, thus steering the stimulation spread back to the STN. Activation volume was calculated by measuring the volume that fell within the 20 mV isosurface of the maximum Fewer contacts are turned on at lower amplitudes if the sensitivity parameter is lower, thus generating more conservative contact settings. (ii) In the lateral case, the lateral contacts, which are facing the IC (colored orange) are kept off or low to avoid IC stimulation. The sensitivity parameter changes the contact amplitudes but the configuration pattern used to steer stimulation away from the IC is maintained across all parameter values. eigenvalue of the Hessian matrix and excluded the lead shaft; activation volume comparisons across leads are shown in figure 4(B). For the central STN position, the Medtronic 3389 was able to activate 39.9% of the STN, whereas the directional leads, the direct STNAcute and Medtronic-Sapiens, were able to activate 62.9% and 57.1%, respectively. When the lead was placed into the lateral position, 2 mm closer to the IC, the Medtronic 3389 was unable to activate a reasonable amount of STN-only 2.9%-without also activating the fibers we isolated from the IC. In contrast, the directional leads were able to activate a substantial amount of the STN: 39.5% for the direct STNAcute and 56.8% for the Medtronic-Sapiens. For the Medtronic-Sapiens, the amount of STN activated is indistinguishable from the center position, demonstrating that judicious contact selection can compensate for lead placement deviation.
We evaluated the algorithm-determined settings by ramping the amplitudes from 0% to 300% and measuring levels of STN activation and IC activation in the NEURON models, as shown in figure 4(C). As the values suggested by the algorithm were scaled below 100%, the amount of STN activation decreased toward zero. As the amplitudes were increased past 100%, the amount of STN activation increased at the expense of additionally activating IC fibers we isolated with tractography. Additionally, when the amplitudes were increased past the suggested values, charge density limits were exceeded in the directional electrode cases. Through this analysis, we determined that the setting produced by the algorithm finds a balance between maximizing target activation, complying with charge density limits, and limiting activation of undesired regions that might induce side effects.
While the results throughout the paper are for voltagecontrolled stimulation, the algorithm is independent of the methods used in the FEM solution calculation, and principally independent of using voltage or current control. However, showing both voltage and current values may be valuable since new directional leads primarily use current-controlled stimulation. The algorithm-determined, optimal voltage settings are listed in table 1 along with the corresponding multichannel, independent current control results for each contact. The direct STNAcute activated more STN than the Medtronic-Sapiens, potentially due to the favorable orientation of the directional contacts; only one column of directional contacts faces the IC. For the lateral placement, the Medtronic-Sapiens electrode was able to activate a similar STN volume as in the centrally placed lead case while the other leads experienced reduced STN activation, especially for the Medtronic 3389 at 2.9%. (C) Contacts were ramped in amplitude to verify that the algorithm produces optimal settings given the constraints used. Lowering the amplitudes results in less STN activation, but increasing the amplitudes past optimal elicits IC activation. Data points represented by an 'x' notation mean the contact configuration exceeds safe charge density limits.
The relationship between the voltage and current control results is dependent of the impedance of each electrode contact. A lower impedance value, as used in this study, will result in a larger current draw across the electrode. A low impedance model would more likely represent the acute period immediately after implantation, while a higher impedance model would represent stimulation observed in the weeks or months afterward.

Targeting constraints
At minimum, one targeting constraint must be used-for example, targeting a nucleus of interest-and any additional constraints thereafter will be used for increased targeting specificity. In this section, we investigate how the use of targeting constraints alters the results of the optimization algorithm.
The constraints we used for this example patient were: targeting the STN, limiting charge density, avoiding the IC, and avoiding the thalamus. Figure 5 visualizes and quantifies the influence of various targeting constraints on contact configurations, amplitudes, and the activation metrics. It is worth noting that when charge density is not controlled, contacts will be turned up to the maximum allowable voltage. In the absence of avoidance constraints, knowing the target region is sufficient for the selection of the contacts that are most useful in stimulating the target region. The addition of avoidance criteria increases specificity by reducing activation of certain areas outside the STN. In cases in which avoidance criteria were added, the contact configurations were adjusted to have fewer contacts enabled and/or voltage amplitudes on select contacts were reduced.
Minor changes to the contact configurations can result in drastic changes to field shaping, so smart selection of contacts can avoid certain regions altogether while still maintaining high levels of STN activation. For example, in the center placement case, when charge density is limited and there are no other avoidance constraints, STN activation is 68%, but IC activation is high at 81%. By adding the IC avoidance criteria, the activation volume is reduced by only 10%, but IC activation is reduced to 0%. Similarly, in the case of the lateral placement, activation of the STN is at 60% and IC activation is at 27%, but once the targeting constraints are added, 57% of the STN is still activated whereas no IC fibers are activated.

Computational efficiency
We used computational time measured across ten runs of the algorithm to quantify its efficiency. The optimization algorithm runs in near real-time, ranging from 0.1 to 10 s based on lead complexity and the number of nodes that define targeting regions. The algorithm was run on a personal computer with a 4.0 GHz quad-core Intel Core i7, 32 GB of RAM, and 8 GB of 1867 MHz DDR3 memory. As the number of contacts increases, the computational time increases linearly, as seen Table 1. Algorithm-determined voltage control and current control optimal settings in units of volts and milliamps, respectively. Each pair of values represents voltage control and current control amplitudes for contacts, as shown and arranged in figure 4(A). C0-C3 indicates contacts labeled from the distal end to the proximal end in the Medtronic 3389 lead. R0-R3 and R0-R9 indicate rows of contacts for the direct STNAcute and Medtronic-Sapiens leads, respectively. In the direct STNAcute, R0 and R1 each have three columns, which represent the directional electrodes, while R2 and R3 have only one column, which represents the cylindrical contacts. Impedance values for the electrodes are approximately 500 Ω for the cylindrical electrodes, 1.5 kΩ for the direct STNAcute directional electrodes, and 2.5 kΩ for the Medtronic-Sapiens electrodes.
Voltage source (V)/Current source (mA)   figure 6. Additionally, each avoidance constraint adds complexity to the algorithm proportional to the number of nodes over a region where the activating function must be checked to be below the sensitivity parameter. A longer computation time is needed for thalamus avoidance than for IC avoidance because the thalamus segmentation we used has about 44 000 nodes compared to fewer than 7000 nodes for the IC. This computation time can be reduced by decreasing the number of points that belong to a region of avoidance while still maintaining the shape of the region. Computation time was approximately 0.1-0.2 s without any avoidance constraints. In most cases computation time was below 1 s, with the exception of the Medtronic-Sapiens case in which the thalamus was used as a constraint. Across all conditions tested, the algorithm computed an optimal contact configuration quickly and without the need for highcomputing resources, demonstrating the feasibility of running this algorithm on a personal computing device.

STN Center Placement
However, it is important to note that preliminary work is necessary to build the FEM model, including the generation of MRI, CT, and DTI scans for the patient as well as the processing required to segment nuclei, generate anisotropic conductivities, identify lead placement, generate tractography, and build the finite element mesh. The construction of a patientspecific finite element model takes several hours and manual interventions may be required. While many of the necessary steps can be performed automatically (e.g. image registration and FEM computations) or semi-automatically (image segmentation, tractography, mesh generation), supervision and manual corrections are still necessary. Thus, it is important to keep in mind that for the algorithm to return patient-specific contact settings, a patient-specific model must first be built, and the time that must be taken to build the model will vary, depending on the complexity of the model. Our hope is that generic models can be developed that can offer quick results for most patients, together with more complex models for challenging programming scenarios.

Discussion
The success of DBS as a surgical therapy, and the emergence of novel DBS lead designs with increased numbers of contacts, call for improved programming techniques that do not rely on manual programming. Our optimization model of DBS programming was applied to an STN DBS example patient and three lead designs: the Medtronic 3389 and two directional leads, the direct STNAcute and the Medtronic-Sapiens. We achieved our goal of creating an automated procedure for selecting contacts and their amplitudes so as to maximize activation of a target while obeying a variety of user-defined targeting constraints.

Targeting using directional electrodes versus cylindrical electrodes
An interesting result from our exploration of the optimization algorithm was the difference in targeting ability by each of the leads in the example scenario of targeting the STN for Parkinson's disease. Even when the lead was placed centrally within the STN, both directional leads activated more STN volume (60%) than the Medtronic 3389 (40%). Thus, even for properly located leads, the directional leads activated an STN volume 50% greater than that activated by the Medtronic 3389. This increase in activation can be attributed to the fact When avoiding the IC tracts we generated, the laterally facing contacts are not used. (B) The contacts that face away from the IC are preferentially activated in the 2 mm lateral placement case. When charge density limits are considered, there is only a 3% reduction in STN activation when bringing the number of modeled IC fibers activated from 27% to 0%. that, even with the portion of the IC we modeled being located 3 mm away, directional leads can still be used to avoid IC activation in the case of central lead placement. The ability of directional electrodes to steer stimulation is largely due to the selection of sources on the portion of the lead that faces away from non-target regions. The presence of the shaft creates an insulating barrier that limits electric potential spreading to areas behind the lead shaft.
A more generally accepted benefit of directional leads is apparent when the electrodes are not centrally placed in the target. We used an example case where the lead was moved 2 mm lateral of the STN center, and thus closer to the IC-a region associated with side effects when stimulated (Ashby et al 1999, Tamma et al 2002, Krack et al 2003. The Medtronic 3389 lead was unable to directionally steer stimulation and could thus activate a mere 2.9% of STN without exciting regions inducing side effects. In contrast, the direct STNAcute and the Medtronic-Sapiens achieved 39.9% and 56.8% activation of the STN, respectively. Thus, the directional leads could activate as much target tissue when they were 2 mm off-center as the cylindrical contacts could activate when they were precisely on target. The use of directional leads might help mitigate negative outcomes resulting from imperfect lead placement, if the contacts could be programmed appropriately.
The amount of STN activation will vary based on the radial orientation of the contacts in directional leads. Each directional electrode has a radial resolution, so changes to the radial orientation of the contacts may either slightly improve or hinder targeting. The directional contacts on the direct STNAcute have a radial resolution of 60 degrees; the staggered arrangement of contacts on the Medtronic-Sapiens gives the electrodes 22.5 degrees of resolution (Martens et al 2011, Pollo et al 2014, Rossi et al 2016. For the center placement, the direct STNAcute activated more STN volume because it was placed radially to maximize STN activation, since only one column of directional contacts faced the IC. A 60-degree rotation of the direct STNAcute might decrease the amount of STN activated and bring the amount of STN activation down to the level of the Medtronic-Sapiens for the center placement case.

Approaches to DBS programming
Historically, DBS programming has been performed primarily through manual approaches that involve trial-and-error testing of settings. Clinicians who carry out programming operate using established protocols, and program within a range of recommended settings that have resulted in positive clinical outcomes in the past. For novel disorders, programming might rely more on straightforward programming strategies to fully explore the options of pulse width, frequency, and voltage combinations. For more established DBS treatments, such as in STN DBS for Parkinson's disease, the recommended protocols begin with commonly successful pulse width and frequency values, with voltage initially being the only free parameter that a clinician manipulates. Clinicians step up the voltage and observe whether or not there is a therapeutic benefit. Depending on the patient response, the clinician can further increase the voltage to maximize therapeutic benefit, or decrease voltage if the patient shows side effects. If the patient experiences only limited clinical benefit, the pulse width and frequency can then be modified in the effort to improve clinical outcomes (Volkmann et al 2006). The process of trial-and-error programming is timeconsuming, ranging from 18 to 32 h across numerous postoperative visits in the first year (Hunka et al 2005). Another challenge to trial-and-error programming is that, once a setting is changed, the patient's symptoms do not change instantly (often referred to as the wash-in time). Depending on the patient, tremor, bradykinesia, and rigidity can take seconds or several minutes to stabilize after a DBS setting is switched (Krack et al 2002, Temperli et al 2003. This leads to great difficulty in assessing whether clinician-determined settings are indeed effective in the treating symptoms of during the appointment window, and often requires multiple programming appointments. The purpose of automated DBS programming algorithms is to provide a starting point for clinicians, in order to reduce the time burden and improve the accuracy of DBS programming. Such algorithms are based on the computational modeling of target structures and optimizing the spread of activation to target regions. Previous work has highlighted the use of the volume of tissue activated (VTA) as predictions of activation with NEURON models (Butson et al 2007a). The goal of the VTA is to visualize an activation volume for a set of DBS stimulation parameters. The VTA is built around NEURON models, which consider the nonlinear firing dynamics and temporal influences of pulse duration and frequency, to characterize and visualize activation. However, using a NEURON model approach does not allow for fast optimization because Figure 6. Summary of computation times required for all lead designs, given different targeting constraints-typically 0.1-10 s. Running the algorithm for the Medtronic-Sapiens is the slowest because it has the largest number of contacts to optimize. The thalamus avoidance criteria substantially increase computation time because there are 44 000 thalamus nodes at which the algorithm must limit second difference values below the sensitivity threshold. the individual multicompartment cable models are slow to compute.
A number of software programs have been created to visually guide clinicians to reach therapeutic settings based on the representation of VTA spread in target and non-target areas (Frankemolle et al 2010, Butson et al 2013, Pourfar et al 2015. In order to visualize the VTA quickly, a database of all possible VTA shapes must be created for different amplitudes, pulse widths, and impedance values. To work around the computational burden, groups have approximated the VTA with generic spheres known as a volume of activated tissue (VAT), in which the radius of the sphere is determined by a set of fit parameters (Mädler and Coenen 2012). The VAT, and now the VTA, are used in Lead-DBS software to visualize stimulation spread in DBS (Horn and Kühn 2015). Another technique has been developed to reduce computation time in estimating the VTA, using the radius information from pre-computed VTAs as the training set in an artificial neural network (Chaturvedi et al 2013). The use of the artificial neural network technique has been implemented in the StimVision clinical software tool to facilitate tractography-based DBS targeting (Noecker et al 2017). The algorithms discussed above can quickly estimate activation volumes to help manual programming efforts, but they do not generate optimal contact configurations and amplitudes based on patient-specific data. While the above approaches still rely on the manual determination of optimal settings, in the future these software packages could utilize these VTA approximations to produce optimal settings. One software package, StimExplorer, already provides recommendations of DBS parameters calculated from a pre-computed database of VTAs so that settings produce maximal volume inside a target while limiting stimulation outside the target , Butson et al 2007a. However, this approach relies heavily on the pre-computation of NEURON models for all voltage, pulse width, and impedance combinations.
To avoid pre-computation of NEURON models, the activating function is useful for approximating neuron firing. With this metric, one can use linear convex optimization, which was first introduced by Xiao et al (2016), for neural targeting with directional DBS leads. The methods used in Xiao et al (2016) consider a few approaches in which the difference between the activating function and the maximal activating function for each contact is minimized. While using the linear convex optimization method is a quick way to compute settings for multiple electrodes contacts, it relies on the calculation of the maximal activating function specific to a set of pre-defined fiber pathways. We expanded on this method by using manipulations of the Hessian matrix to approximate activation of nuclei in all directions, beyond a pre-defined neuron pathway. We additionally outlined a method using the Hessian for target avoidance, whether the target is a fiber tract or a nucleus. Other algorithms use particle swarm optimization machine learning to identify contact parameters to maximize target activation, minimize side activation of side effect regions, and minimize power consumption (Peña et al 2017). Unlike with linear optim ization methods, the solution of the particle swarm algorithm is able to minimize and maximize constraints simultaneously and, through the use of a Pareto front, to generate multiple contact configurations with varying levels of prioritization in targeting regions of interest, avoiding regions inducing side effect, and limiting power consumption. This offers flexibility by offering multiple solutions in the same run. However, by using the machine learning approach, re-runs of the same criteria might not converge on the same results and have a run time in the order of several minutes.
Our approach used the second difference as an estimate of activation, so as to avoid the necessity of pre-computation with NEURON models. The only pre-computation that is required is of the voltage FEM solutions for each contact, which can be linearly combined to a multitude of possible parameter combinations. Using the second difference allows our algorithm to include a number of features, such as flexibility with targeting criteria, patient specificity, and the generation and adjustability of results based on patient sensitivity in real-time.

Flexible targeting of nuclei, regions of interest, and fiber tracts
In this study, DBS for Parkinson's disease, in which the target is the STN, was used to demonstrate the algorithm's performance. Targeting structures consisted of segmentations of the STN and thalamus and tractography of the IC from the atlas brain model. However, we emphasize that this algorithm can be applied to other DBS scenarios as well, such as for the globus pallidus pars interna (GPi) which is a target for DBS for Parkinson's disease (Anderson et al 2005), as well other disorders such as essential tremor and dystonia (Benabid et al 1996, Vidailhet et al 2005. While most movement disorders have nuclear targets, the algorithm can also accommodate programming for other disease states with fiber tracts as targets. These include: the fornix bundle for Alzheimer's disease (Laxton et al 2010, Laxton andLozano 2013); the subcallosal cingulate (SCC) white matter, the inferior thalamic peduncle, the medial forebrain bundle, or the anterior limb of the IC for treatment-resistant depression (Jiménez et al 2005, Mayberg et al 2005, Malone et al 2009, Bewernick et al 2010, Schlaepfer et al 2013; and the medial dorsal tegmental tract (DTTM) for traumatic brain injury (Baker et al 2016).
We have enabled our algorithm to allow for any combination of targeting and avoidance constraints, whether the targets are a single nucleus, multiple nuclei, a user-defined region determined by a clinician, or fiber tracts. The flexibility in choosing targets and constraints allows the algorithm to be applied indiscriminately to a variety of neurological indications if the computational model is first constructed for them. If future clinical studies isolate regions of stimulation that correspond to positive outcome measures, we can use those regions as targets instead of nucleus boundaries. Additionally, because the algorithm is quick to compute, an adjustment of targeting parameters may be made in real-time, thus giving a clinician the option to toggle targeting constraints to see which might produce a better result for the patient. Quick targeting adjustments might not only enable better patient treatment, but could also be used as a research tool to investigate the influence of stimulation of specific regions on clinical outcome scores.

Patient specificity
The algorithm is based on a patient-specific FEM model generated from MRI and DTI patient data. Because the response to DBS might vary from patient to patient, the algorithm accommodates the unique qualities of a patient's brain morphology, nuclei size and shape, fiber tracts, tissue anisotropy, and lead placement (Deuschl et al 2003, Tarsy et al 2003, York et al 2009. Incorporating these elements allows for patientspecific DBS settings tailored to optimize neuronal activation. The level of patient specificity can be adjusted based on the patient information available. For example, many patients do not receive DTI scans because scans are expensive and time-consuming. As a result, it might not be possible to generate anisotropic tissue conductivities or tractography data for those patients. One possible solution would be to register an anisotropic atlas brain to the patient brain. However, as the complexity of the model grows, more time must be dedicated to building the patient-specific model. Alternatively, it might be more desirable to simplify the FEM model, for example, by using isotropic tissue conductivities to develop robust yet generic models for quick application to patients. Simplification of the FEM model has been shown to reduce the accuracy of simulated electric potential spread, especially in the case of using isotropic conductivities (Chaturvedi et al 2010, Howell andMcIntyre 2016). Understandably, in simplifying the patient model, there is a trade-off between speed in generating an optimal solution and accuracy of the solution.

Adjustability of optimized contact configurations and amplitudes
Because a number of assumptions underlie any computational model, it is possible that our model does not capture all qualities of patients and their responses to DBS stimulation. Therefore, we included the sensitivity parameter that allows for adjustments to contact configuration and amplitudes while conserving the shape of the stimulation fields. The sensitivity parameter should be viewed as a sliding scale that can be adjusted in real-time during programming. For safety purposes, it should start low and be ramped up as tolerated by the patient. For example, some patients may be more sensitive to stimulation than our model predicts, and side effects could arise earlier than expected; changes in sensitivity to stimulation could be the result of imperfect segmentation of nucleus geometries or poor estimation of activation thresholds of nearby fiber tracts. A decrease in the sensitivity parameter will result in more conservative contact settings, by reducing the amplitudes of the contacts being used and/or reducing the number of active contacts.
The sensitivity parameter is based on the activating function approximation and represents the threshold of firing. This threshold value can be selected based on literature values.
The value varies depending on axon diameter and position, but commonly falls in the range of 5-30 mV (Rattay 1986, McIntyre et al 2004, Martens et al 2011. For the fibers we modeled as part of the IC, and given our stimulation parameters, we found the threshold to be 20 mV for our specific fiber type and orientation. In the future, it would be useful for studies to fully characterize the relationship between activating function and NEURON models, to serve as guidance for the sensitivity parameter selection. Some studies have already created guides to match activating function thresholds with NEURON models for a wide range of fiber diameters (Peterson et al 2011). However, such studies have looked only at axon activation due to stimulation fields produced by cylindrical electrodes in a single orientation. The second difference relationship is more complex when considering directional electrodes and asymmetric fields, and would become even more complex when looking at bipolar stimulation fields. Further experimental studies are needed to characterize the relation between second difference values and axon firing in the context of directional leads.

Real-time optimization results
The algorithm is able to perform in near real-time, taking only a few seconds to generate contact configurations that maximize target activation while obeying the imposed constraints. The time of computation largely depends on the number of contacts and linear elements that make up the nucleus and fiber regions. In the case of nucleus or fiber tract avoidance, the algorithm must conduct a pointwise check at each location to verify that the second difference is below the sensitivity threshold; a second difference value below threshold implies that the neurons in that region were not activated. As the number of points to check for the threshold crossings of the sensitivity parameters increases, the time to compute an optimal solution increases linearly. The algorithm's approach to targeting any given region is to maximize the second difference averaged across that region; the implementation of averaging over a region enables faster computation times.
In the example case we studied, the Medtronic-Sapiens was the slowest to compute in all conditions, due to its greater number of contacts. Additionally, the thalamus segmentation we used had about 44 000 boundary points and thus resulted in longer computation times when the thalamus was considered as a region of avoidance. In the future, the number of points on nucleus regions could be reduced while still maintaining the important features of the thalamus, to improve computation time. The contribution of nuclear regions and fiber tracts on computational time should be taken into consideration when setting up the patient-specific model to retain its near realtime capabilities.

Limitations
The contact configurations and amplitudes are generated based on regions defined in the computational model. With any computational model, a number of assumptions must be made. Errors in the segmentations of nuclei regions or the tractography data can alter the results of the DBS programming algorithm, and the settings may not fully avoid regions when applied to the patient for this reason. To combat such deviations, we added the sensitivity parameter to accommodate the differences that might exist between the computational model and the patient's brain. Additionally, the algorithm is able to produce only amplitude settings and does not optimize frequency or pulse width. Changes in the frequency and pulse width values can be compensated for through modification of the sensitivity parameter. The manner in which changing those DBS settings affects activation sensitivity could be estimated from the literature or simulated with time-dependent multicompartment cable models in NEURON. Ultimately, more experimental research is needed, and we encourage future studies to characterize second difference thresholds for cylindrical and directional electrodes alike.
Another limitation is that the FEM solutions for each contact must be computed in advance of running the programming algorithm. Pre-computing allows for the use of anisotropic conductivities-enabling patient-specific programming-without reducing the efficiency of the algorithm. Pre-computing also allows quick access to spatial voltage data for each contact, and through this strategy, we are able to compute contact configurations in near real-time. The full FEM solutions are precomputed, but only the data in a 20 mm × 20 mm × 20 mm region surrounding the lead is stored. Neural targets that fall outside this area cannot be considered by the algorithm because no stored voltage or Hessian solutions exist for those regions. However, nothing inherent to the algorithm limits us to saving data on a 20 mm × 20 mm × 20 mm grid. So, as computational power grows in the near future, the algorithm could optimize stimulation at the scale of the whole brain, for example, for both bilateral DBS leads at the same time. For the time being, this limitation has minimal impact on programming configurations because it is unlikely that substantial stimulation will reach areas 10 mm away from the electrode.
The Medtronic 3387/3389 leads have approximately 25 000 parameter combinations and 65 active electrode combinations (Kuncel and Grill 2004); it is unreasonable to assume that clinicians are able to sift through the entire parameter space and manually determine optimal contacts. Because the number of parameter combinations exponentially increases with the addition of more contacts, manual programming becomes even more challenging. The goal of the optimization algorithm is not necessarily to produce an absolute setting, but rather to explore the parameter space quickly and provide guidance to programming clinicians for more efficient programming procedures. The settings produced by the algorithm should serve as a starting point for programming clinicians. The clinicians can then make adjustments to algorithm-determined settings, either through modifications on the algorithm parameters or via slight manual manipulations as they see fit. Such a programming structure allows for interactive programming that may expedite DBS programming, allow clinicians to see more DBS patients, enable the use of complex lead designs, and achieve nearoptimal settings tailored to the patient.

Conclusions
We demonstrated that the optimization model for DBS programming we created offers a quick and robust solution to manual programming challenges for commonly used cylindrical electrodes and novel directional electrodes. The algorithm accommodates a range of lead and programming complexity, incorporates various levels of targeting and avoidance constraints, and abides by FDA-imposed charge density limits. Future work will involve applying the algorithm to patients who have been implanted with either cylindrical leads or directional leads to validate algorithm-chosen settings, and determining its impact on clinical DBS programming efficiency.