Clinical deep brain stimulation strategies for orientation-selective pathway activation

Objective. This study investigated stimulation strategies to increase the selectivity of activating axonal pathways within the brain based on their orientations relative to clinical deep brain stimulation (DBS) lead implants. Approach. Previous work has shown how varying electrode shape and controlling the primary electric field direction through preclinical electrode arrays can produce orientation-selective axonal stimulation. Here, we significantly extend those results using computational models to evaluate the degree to which clinical DBS leads can direct stimulus-induced electric fields and generate orientation-selective activation of fiber pathways in the brain. Orientation-selective pulse paradigms were evaluated in conceptual models and in patient-specific models of subthalamic nucleus (STN)-DBS for treating Parkinson’s disease. Main results. Single-contact monopolar or two-contact bipolar stimulation through clinical DBS leads with cylindrical electrodes primarily activated axons orientated parallel to the lead. Conversely, multi-contact monopolar stimulation with a cathode-leading pulse waveform selectively activated axons perpendicular to the DBS lead. Clinical DBS leads with segmented rows of electrodes and a single current source provided additional angular resolution for activating axons oriented 0°, ±22.5°, ±45°, ±67.5°, or 90° relative to the lead shaft. Employing multiple independent current sources to deliver unequal amounts of current through these leads further increased the angular resolution of activation relative to the lead shaft. The patient-specific models indicated that multi-contact cathode configurations, which are rarely used in clinical practice, could increase activation of the hyperdirect pathway collaterals projecting into STN (a putative therapeutic target), while minimizing direct activation of the corticospinal tract of internal capsule, which can elicit sensorimotor side-effects when stimulated. Significance. When combined with patient-specific tissue anisotropy and patient-specific anatomical morphologies of neural pathways responsible for therapy and side effects, orientation-selective DBS approaches show potential to significantly improve clinical outcomes of DBS therapy for a range of existing and investigational clinical indications.


Introduction
Numerous individuals living with medication-refractory neurological and neuropsychiatric disorders have benefited from deep brain stimulation therapy (DBS), which is a treatment involving continuous or intermittent electrical stimulation through one or more leads of electrodes chronically implanted within the brain [1]. Deep brain targets of DBS therapy are embedded within functionally complex networks of axonal afferents, efferents, and fibers of passage that when stimulated can alleviate symptoms. However, stimulation within these networks also can induce adverse side effects [2,3], with occurrence dependent upon the lead position and programmed stimulation settings. There is a strong clinical need to develop patient-specific approaches that can more selectively activate axonal pathways within a DBS target to not only limit stimulation-induced side effects but also increase the overall effect size of the therapy.
As an example of the importance of pathway-selective activation, DBS therapy for Parkinson's disease (PD) often involves implanting a lead of electrodes in the dorsolateral subthalamic nucleus (STN) in order to target axonal pathways coursing approximately perpendicular to the DBS lead (e.g. STN efferents and the so-called hyperdirect pathway afferents). The STN, however, is encapsulated on its lateral edge by the corticospinal tract of internal capsule (IC), which traverses approximately parallel to the DBS lead. If the stimulation intensity required to alleviate parkinsonian motor signs also activates this IC tract, the result is typically a muted therapy with the emergence of involuntary muscle contractions and/ or paresthesias [4]. DBS therapy can also require activation of multiple pathways with different axonal orientations relative to the lead of electrodes. For example, investigational stimulation of the subcallosal cingulate white matter for treatmentresistant depression is thought to require activation of three white matter tracts with different axonal orientations relative to each other: the forceps minor, the uncinate fasciculus, and the cingulum bundle for a sustained therapeutic effect [5].
Previous studies have shown that electrode geometry (i.e. circumference, height, and perimeter of the electrode) can impact stimulus-induced action potential thresholds within axons of varying orientations relative to a monopolar electrode [6,7]. Notably, taller cylindrical electrodes (7.5-8.5 mm) provided greater selectivity for axons oriented perpendicular to the lead and shorter cylindrical electrodes (2-2.5 mm) were more effective at activating axons coursing parallel to the lead [7]. Past studies have shown how maximizing the second spatial derivate of an electric field is most effective for axonal activation [8,9]. More recently, we applied this concept to achieve orientation-selective activation by controlling the primary direction of the electric field in a plane using a threeelectrode bundle [10]. Computational modeling and functional magnetic resonance imaging (fMRI) analysis of corpus callosum stimulation in rodents demonstrated that aligning the electric field to the axonal trajectories of the corpus callosum resulted computationally in the strongest degree of axonal activation and experimentally in the largest fMRI-based cortical response. Further, aligning the primary direction of the electric field tangential to the corpus callosum resulted in the smallest response both computationally and experimentally. Here, we extend these results through computational models, showing how to increase the selectivity of pathway activation by orienting electric fields along and around clinical DBS leads.

Materials and methods
This study consisted of both conceptual and patient-specific computational multi-compartment cable models of axons [11] in the context of clinical DBS leads from three medical device manufacturers (figure 1). The conceptual modeling study involved placing axons in multiple radial orientations around, along, and distal to each DBS lead. The patient-specific modeling study investigated the ability to more selectively stimulate the hyperdirect pathway (HDP) coursing into the STN while minimizing direct activation of the corticospinal tract of IC. In both cases, charge-balanced, pseudo-monophasic rectangular pulses were used consistent with the pulse waveform shapes from clinical neurostimulators.

Conceptual modeling framework
For the conceptual modeling study, three commercial DBS lead designs were incorporated into a finite element model (FEM). These included: the Medtronic 3389 (figure 1(a)), consisting of four cylindrical contacts; the Abbott 6172 ( figure  1(b)), consisting of a cylindrical contact on either side of two rows each with three segmented contacts; and the Boston Scientific 2202 (figure 1(c)), consisting of a tip electrode, two rows each with three segmented contacts, and one cylindrical contact.
The FEMs of each lead design were constructed in COMSOL Multiphysics v5.2a based on manufacturer specifications and surrounded by a 0.25 mm thick isotropic encapsulation layer [12,13], and spherical homogeneous bulk tissue (70 mm-diameter) representing brain tissue with conductivity and permittivity parameters for both the bulk tissue and encapsulation layer calculated from the frequency dependent Cole-Cole dispersion function [14]. The functions were evaluated according to a median frequency of 3.049 kHz from the typical charge-balanced, pseudo-monophasic rectangular pulse waveform used in clinical current-controlled pulse generators. This yielded conductivity measures for bulk tissue (grey matter, σ = 0.106 S m −1 , ɛ r = 65,898), and encapsulation layer (white matter, σ = 0.065 S m −1 , ɛ r = 29 790). The FEMs were meshed using Delaunay Tessellation and variable resolution tetrahedral elements with highest resolution near the electrode (encapsulation layer = 0.08-0.25 mm) and lower resolution further from the lead (bulk tissue = 0.17-7.23 mm). Mesh refinement studies were completed to confirm a higher resolution mesh (~1.5 million elements) resulted in <.3% change in interpolated voltages within 20 mm of the lead and <1% change over the full tissue model. Final meshes consisted of 279 000-281 000 elements depending on the lead design complexity. In each model, the outer surface of the bulk tissue sphere was set as ground via a Dirichlet boundary condition of zero volts and the insulated lead shaft was assigned a Neumann boundary condition of zero flux. Current-controlled stimulation was incorporated in the model using a normal current density, which was calculated by dividing the maximum stimulation pulse amplitude (1 mA) by the geometric surface area of the electrode. The static FEM solved Poisson's equation to calculate the tissue voltage distribution for each electrode configuration. Superposition was used for stimulation paradigms employing multiple electrodes, which was confirmed to be accurate by comparing the interpolated FEM output when running bipolar and multi-cathode solutions versus using superposition of interpolated solutions from each electrode.
Multi-compartment cable models of myelinated axons (40 mm length) were radially distributed adjacent to each of the DBS leads and shifted along, around, and away from the lead as shown in figures 2(a) and (b). Axons were modeled with a 2 µm diameter and each axon was populated with compartments representing nodes of Ranvier, myelin attachment segments, paranodal main segments, and internodal segments connected by an axial resistance. Axon compartment properties were consistent with the multi-compartment axon cable models developed previously [15] and available on ModelDB. Extracellular tissue voltages corresponding to each membrane compartment were interpolated from the FEM solutions and were used to scale the applied pulse waveform (90 µs cathodic phase, 400 µs interphase, 3 ms charge-balancing anodic  phase). To account for the capacitive effects of the electrodetissue interface, the pulse waveform shape was taken from a previous recording of a stimulation artifact measured adjacent to a DBS lead in vivo in the context of current-controlled stimulation, resulting in a filtered waveform with the highfrequency components removed by the electrode-tissue interface [16]. The time-varying extracellular voltage profiles were applied concurrently into each axonal compartment within the NEURON v7.3 programming environment [17]. This was performed by perturbing the axonal membrane voltage using the e_extracellular mechanism. Action potential counters were attached to each node and simple thresholding was used to detect stimulus driven spike activity. An axon was considered 'activated' if an action potential was detected within 5 ms following the stimulation pulse. A binary search algorithm was employed to determine the stimulus amplitude threshold (±0.025 mA) to activate a given axon for a given stimulation setting.
Axonal orientation selectivity for each threshold was quanti fied by fitting a sigmoidal curve to the data (figure 2(c), equation (1)). The selectivity angle was calculated as the axon angle with the minimum activation threshold, which implies that at some stimulation amplitude orientation-selective stimulation can be achieved. The selectivity gain, which is a measure of the window size between orientation-selective activation and omni-directional activation, was quantified as the sigmoid gain (a) divided by the offset from zero (d), which served to normalize the gain (equation (2)). Notably, a gain of 0.5 could come from a = 1.5 and d = 3 mA or a = 3 and d = 6 mA, for example. Results, therefore, include both the normalized gain and the minimum threshold. FEM simulations included runs with either single source (Medtronic 3389 and Abbott 6172) or multi-channel independent currentcontrolled (MICC) sources (Boston Scientific 2202), which enabled unequal distributions of current through multiple electrodes. , and a whole brain diffusion-weighted imaging (DWI) scan (1.5 mm isotropic, 54 directions, b-value = 1500 s mm −2 collected with anteriorposterior and posterior-anterior phase encoding directions). Processing of the DWI was performed using FSL [18] and included eddy current correction using FDT [19], estimation of diffusion parameters using BEDPOSTX [20], and fitting of the diffusion tensor model at each voxel using DTIFIT [20]. All imaging data sets were co-registered to diffusion space. The patient-specific STN model surfaces were obtained by manual segmentation of the 7 T T2 coronal images by three experts. A post-operative CT image (0.6 mm isotropic-Siemens Biograph) was acquired in order to determine final electrode location. The CT image was acquired one month after the surgery in order to enable any potential brain shift accumulated during the surgery to resolve. Finally, in order to visualize each of the electrode locations, the patient's STN and the white matter tracts together, we co-registered with all the images. Because it is very difficult to register the CT image with the 7 T T2 slab image, we chose to use a T1 image from the 1.5 T images as a reference. We first co-registered the 7 T T2 and the post-operative CT to the 1.5 T T1 and then we co-registered the 1.5 T T1 to the 7 T DWI. This enabled all the images to be aligned to the diffusion space.
The imaging data were used to create patient-specific FEMs. The geometry of the model included the DBS lead implant, a 0.25 mm thick encapsulation layer surrounding the lead, and a smoothed surface reconstruction of the head volume. Surface reconstructions of the brain from each subject were derived from the T1-weighted imaging using manual brain extraction and smoothed in Avizo v9.3 (FEI). Brain and head volume reconstructions, were imported using COMSOL Multiphysics CAD import module. The DBS lead was manually constructed within COMSOL and oriented in accordance with the final lead location as determined from postoperative CT imaging [21,22].
Inhomogeneous and anisotropic material conductivities were then assigned according to patient-specific normalized volume diffusion tensor maps [23,24]. The anisotropic conductivity (σ i ) maps were calculated by scaling the diffusion tensor eigenvalues (d i ) for each voxel by the ratio of the geometric volume of the isotropic conductivity tensor to the geometric volume of the diffusion tensor. This equation is given by: in which σ iso was the isotropic conductivity of white matter (σ = 0.065 S m −1 , ɛ r = 29 790) or grey matter (σ = 0.106 S m −1 , ɛ r = 65 898) as defined by the frequency-dependent Cole-Cole dispersion functions at a frequency of 3049 Hz [14,25], and d 1 , d 2 , and d 3 were the diffusion tensor primary, secondary, and tertiary eigenvalues. The FEM was again meshed using tetrahedral elements and Delaunay triangulation with finer resolution near the lead (0.07-0.2 mm) and up to 20 mm at the scalp, final meshes consisted of 717 000-802 000 quadratic tetrahedral elements with an average element quality 0.66 for all patients. Text files containing the brain tissue conductivity tensor maps and relative permittivity values were imported into COMSOL as spatially dependent variables. The tissue properties were assigned by interpolating the tissue conductivity tensor map and relative permittivity onto the finite element mesh using the nearest neighbors function. The material properties for the DBS lead encapsulation (σ = 0.065 S m −1 , ɛ r = 29 790), cortical spinal fluid (σ = 2 S m −1 , ɛ r = 109), and bulk head tissue (σ = 0.23 S m −1 , ɛ r = 1) were calculated from the frequency dependent Cole-Cole dispersion function and values previously identified [13,14,25]. In each model, the head surface and electrode shaft were assigned Neumann boundary condition of zero flux and the base of the head volume was defined as ground via a Dirichlet boundary condition of zero volts. Current-controlled stimulation was again applied using normal current densities.
Multi-compartment axon models of the corticospinal tract within the IC were segmented through probabilistic tractography performed using FSL [18]. The probabilistic tractography approach was guided by seed mask and waypoint masks, which were identified within each subject's T2-and T1-weighted images. A waypoint mask was created in the primary motor cortex of the hemisphere ipsilateral to the DBS lead, and a seed mask was generated within the ipsilateral IC, slightly posterior and ventral of the STN in diffusion space. These masks were transformed into diffusion coordinates using FLIRT [19,26]. The transformed masks and the output of BEDPOSTX were used by PROBTRACKX [27] to create probabilistic maps of the motor region of the IC. Thresholded probabilistic maps were used to generate a digital surface reconstruction that enabled building a population of multicompartment axon models of the IC. The HDP was modeled as horizontal extensions of the IC feeding into the STN as shown from histological tracing [28].
Each IC fiber tract reconstruction was pseudo-randomly populated with 100 multi-compartment axon models and the HDP pathway was pseudo-randomly populated with 200 axons. Axons intersecting the DBS lead or the encapsulation layer were excluded resulting in 94-98 axons in the IC and 86-157 axons in the HDP. The populations of multicompartment axons were constrained by cross-sectional contours from the digital surface reconstruction of the fiber tracts. Axons were modeled with a myelinated diameter of 2 µm for the IC and 1 µm for the HDP pathway in alignment with previous histological studies [29][30][31]. Parameters for 1 µm axons were adapted from [32,44] (node length = 1 µm, MYSA length = 3 µm, FLUT length = 6 µm, STIN length = 43.7 µm, axon diameter = 0.8 µm, node diameter = 0.7, number of myelin lamelle = 20) [15,32]. Stimulus pulses were again applied to these multi-compartment axons as described in section 2.1 with thresholds for activation of each axon found for monopolar, bipolar, and multi-cathode stimulation paradigms.

Conceptual models of the Medtronic 3389 lead
Orientation-selective activation was investigated for the fourcontact Medtronic 3389 lead (C0-C3, figure 3(a)) using combinations of anodes and cathodes. Axons were radially distributed adjacent to the DBS lead with axon angle represented by sigma (σ). The threshold for each axon was determined for two pulse paradigms with the goal of comparing electric fields lines parallel and perpendicular to the DBS lead. In the first case, cathodes were applied to contacts C1 and C2, whereas in the second case, a bipolar configuration was applied with contact C1 as a cathode and contact C2 as an anode ( figure 3(b)). Using a total stimulus threshold of 2.0 mA, for example, the two cathode configuration resulted in orientation selectivity for axons perpendicular to the DBS lead (σ = −15° to 15°) whereas the bipolar configuration resulted in selectivity for axon oriented parallel to the lead (σ = 60° to 120°) (figure 3(c)). The thresholds for activation of all radially distributed axons are shown in figure 3(d).  4(a)). Activation threshold curves for each axon orientation and location distribution as well as for each DBS configuration were fit with a sigmoid curve to calculate the angle of selectivity (σ) and gain (a/d). Single contact monopolar ( figure 4(b)) and dual-contact bipolar (figures 4(c)-(e)) configurations consistently showed selectivity for axons oriented parallel to the DBS lead with average selectivity gains (a/d) of The regions of perpendicular selectivity appeared between the two active contacts and moved away from the lead as the distance between the cathodes increased (figures 4(g) and (h)). The gain value was lowest near the active electrode(s) and increased for more distant axonal elements. This can be explained by the results of figure 4(b), which showed that, for example, at a radius 0.5 mm from the lead, the activation threshold for parallel axons was constant while the threshold for perpendicular axons greatly increased moving along the z-axis, thereby increasing the overall gain value.

Conceptual models of the Abbott 6172 lead
Orientation-selective activation was further investigated for the Abbot 6172 segmented DBS array, which is composed of two cylindrical contacts separated by two rows of three segmented contacts in a 1-3-3-1 configuration ( figure 5(a)). Activation thresholds were found for axons radially distributed 1 mm from the DBS lead centered in the middle of four segmented contacts -that is, two contacts in row 2 and two contacts in row 3. When multi-cathode and bipolar configurations were applied to the segmented contacts by grouping contacts C2-C3 and contacts C5-C6 ( figure 5(b)), similar orientation selectivity to that found for the Medtronic 3389 DBS lead was observed ( figure 3(c)). The multi-cathode configuration showed orientation selectivity for axons perpendicular to the lead (σ = −22.5° to 22.5°), whereas the bipolar configuration showed selectivity for axons parallel to the lead (σ = 60° to 120°) ( figure 5(c)). Additional orientation selectivity could be achieved by applying multi-cathode and multianode configurations to independent segmented contacts in rows 2 and 3 ( figure 5(d)). For instance, a multi-cathode configuration through two diagonally separated segmented contacts resulted in selectivity orthogonal to the active contacts (σ = 127.5° to 165°) while a bipolar configuration through  DBS lead with radially distributed axons along (z) and around (θ) the lead. In the orientation activation plots, dot size represents minimum threshold, dot color represents the gain or degree of selectivity, and arrow direction indicates direction of selectivity. Configurations are shown for (b) four cathodes including contacts C2, C3, C5, C6; (c) two cathodes with contacts C2, C3; (d) two cathodes with contacts C3, C5; (e) bipolar with contacts C2, C3 as cathodes and contacts C5, C6 as anodes; (f) bipolar with contact C3 as a cathode and C2 as an anode; (g) bipolar with contact C3 as a cathode and C5 as an anode. the same contacts resulted in selectivity for axons parallel (σ = 45° to 67.5°) ( figure 5(e)). The thresholds for each configuration are shown in figure 5(f).
The segmented contacts on the Abbott 6172 DBS lead enabled additional orientation-selective activation for axons located between the active contacts. This selectivity was characterized by shifting the radially distributed axons in figures 5(c) and (e) along (z) and around (θ) the DBS lead with all axons centered at r = 1 mm ( figure 6(a)). Multi-cathode and bipolar configurations were then applied through the segmented contacts to characterize how the orientation selectivity changes around and along the lead. Selectivity for multicathode configurations (figures 6(b)-(d)) showed selectivity for axons perpendicular to the cathodes. However, axons further along the lead shank from the two cathodes showed selectivity for axons parallel to the lead (figure 6(d)), which was also found for the cylindrical electrodes (figure 4). Bipolar configurations (figures 6(e)-(g)) showed selectivity for axons aligned with the primary direction of the electric field, where that primary direction was along the lead in figure 6(e), perpend icular to the lead in figure 6(f), and diagonal to the lead in figure 6(g).

Conceptual models of the Boston Scientific 2202 lead
Models with the Boston Scientific 2202 segmented DBS lead included capabilities of MICC stimulation, which is a feature of the 16-channel Boston Scientific Vercise System. The ability to apply unequal amounts of current through the segmented contacts provided opportunities to investigate shifting the primary direction of the electric field (PDEF) from σ = 0° perpendicular to the lead to σ = 90° parallel to the lead. Similar to the Abbott 6172 lead models, axons were arranged 1 mm from the lead surface and shifted along (z) and around (θ) the lead ( figure 7(a)). The primary direction of the electric field was shifted from parallel to the lead (figure 7(b)) to perpendicular to the lead (figure 7(f)) in increments of 22.5°. Figure 7(b) shows that when the primary direction of the electric field was parallel to the lead (−90°) the selectivity was also parallel to the lead. As the primary direction of the electric field was shifted to −67.5° in figure 7(c), the angle of selectivity also shifted along the line between the diagonal contacts 3 and 5. This trend continued through figures 7(d) and (e) until the primary direction of the electric field was perpendicular to the lead in figure 7(f) where selectivity was perpendicular to the lead.

Patient-specific models of the Medtronic 3389 DBS lead
Computational models were also constructed using preoperative MR and post-operative CT imaging data from three STN-DBS patients who were implanted with the Medtronic 3389 DBS lead ( figure 8(a)). For each patient, orientationselective stimulation strategies were investigated to limit activation of the motor territory of the IC and enhance activation of the HDP, which courses from the IC to the STN. The rationale of this choice being that the HDP has been considered to be a therapeutic target pathway underlying STN-DBS therapy [33][34][35][36][37][38], whereas direct IC activation is thought to underlie involuntary motor contraction side effects of DBS [39,40]. Activation thresholds for axons in both pathways were quantified and compared in terms of the percent of the HDP that could be activated at the lowest total current intensity before activation of any axon within the IC pathway. Stimulation configurations included single contact monopolar stimulation, multi-cathode stimulation, and bipolar stimulation. In all cases, the contact(s) closest to the HDP were used.
For patient 1, while monopolar stimulation using contact C3 activated 93.0% of the HDP prior to activation of any modeled IC axons, multi-cathode stimulation was able to achieve 98.1% activation of the HDP ( figure 8(b)). Patient 2 also showed the highest HDP activation with a multi-cathode configuration (68.6% activation), whereas a single monopolar configuration and bipolar configuration activated 62.8% and 66.3% of the axons within the modeled HDP respectively ( figure 8(c)). While configurations for patient 3 showed lower activation levels of the HDP compared to those found in patients 1 and 2 ( figure 8(d)), the multi-cathode paradigm again showed the largest effect size (21.7% HDP activation) compared to 14.2% activation using a single monopolar cathode. It should be noted that while the multi-cathode stimulation configuration resulted in a larger electric field comp onent tangential to the active contacts and more consistent with the orientation of the HDP, these multi-cathode configurations often used higher total current.

Discussion
This computational study provides a detailed framework for understanding how electrode configurations on clinical DBS systems can be interpreted in the context of increasing the selectivity of axonal pathway activation within the brain. The models showed that leads with only cylindrical electrodes (e.g. Medtronic 3389 lead) have the capability of more selectively activating axons oriented parallel or perpendicular to the lead depending on the type of bipolar or monopolar configuration applied. DBS leads with segmented contacts enable steering the lowest activation orientation diagonal to the lead (e.g. Abbott 6172 and Boston Scientific 2202 leads). Further, with MICC stimulation, one can effectively steer the electric field to align with any direction of axons surrounding a segmented DBS lead. These approaches are poised to increase the therapeutic windows for DBS therapies.

Mechanism of orientation-selective stimulation
The models indicated that the driving force underlying the orientation-selective activation results was the spatial orientation of an axon relative to the direction of the imposed electric field generated by the DBS electrode configuration and stimulation waveform parameters. Transmembrane currents induced by extracellular electrical stimulation are thought to stem from the second spatial difference in the extracellular tissue electric potential along axonal membrane compartments. Thus, to maximize the depolarizing transmembrane currents within an axon, one would want to maximize the second spatial difference in the extracellular electric potential along the axon [8,9,41], which effectively means aligning the electric field with the primary direction of the axonal trajectory. This concept was recently modeled and empirically validated in rodents with fMRI [10]. A similar concept has been discussed in the literature in terms of axonal bends that create local 'hot spots' of strong depolarizing membrane currents as shown in the globus pallidus (GPi) for DBS applications [42] and in the cortex for transcranial magnetic stimulation applications [43].

Cylindrical contacts-comparisons between single and multi-contact monopolar stimulation
For the cylindrical contacts (1.5 mm height) used in clinical DBS leads, simulations with monopolar stimulation through one cylindrical contact showed highest selectively for activating axons parallel to the lead shank. The height of the active electrode contact is an important design component as has been shown previously [7]. Cylindrical electrodes with shorter heights are more likely to activate axons parallel to the lead, whereas longer heights more readily activate axons tangential to the lead. Given the results from our study, this is likely due to the difference between electric field lines having a radial distribution for short contacts ( figure 9(a)) and more of a perpendicular distribution away from longer contacts in the case of monopolar stimulation in which the return electrode is distant from the DBS lead. We extended this concept with clinical DBS leads by grouping electrodes together in two cathode configurations. This effectively created an electrode height of 3.5 mm such that electric field lines were primarily perpendicular to the electrodes as shown in figure 9(b). In contrast, using bipolar stimulus configurations, the electric field lines were oriented parallel to the lead both between the active contacts as well as further away from the active contacts, leading to the selectivity for axons parallel to the lead (figure 9(c)).
In the case of the Medtronic 3389 DBS lead, dual cathode configurations indeed resulted in orientation-selective activation of axons oriented perpendicular to the DBS lead. However, upon moving the active cathodes apart (e.g. C0 and C2, or C0 and C3), virtual anodes were generated such that the electric field lines became more radial in orientation around the active contacts, which not only muted the perpendicular orientation selectivity at each active contact but also shifted this orientation selectivity radially farther from the DBS lead. It is also important to note that for the single contact monopolar stimulation, the gain value for the parallel orientation was weak immediately adjacent to the active contact suggesting that the 1.5 mm height electrodes used clinically conveniently activate a wide range of axon orientations. However, farther along the lead shank, the orientation selectivity gain becomes much stronger for axons parallel to the lead, presumably by the presence of virtual anodes in the tissue [44,45] and the distance effect [46]. That is, for axons centered farther from the active contact, only axons oriented parallel to the lead will course through the high electric field region immediately adjacent to the active electrode(s).

Segmented leads
For a DBS lead with segmented contacts, similar orientation selectivity for axons oriented parallel and perpendicular to the lead could be achieved as with cylindrical electrodes by simply grouping the segmented contacts. Further selectivity beyond 0/90° relative to the lead could be achieved by applying anodes and cathodes using a single source to the segmented contacts either around the lead (figures 6(c) and (f)) or at an approximately ±45° angle relative to the lead (figures 6(d) and (g)). Additionally, by splitting a single source into two cathodes with a single anode, stimulation could be applied at ±22.5° and ±67.5°. It should be noted that for these segmented contacts, the effect was seen close to the lead as shown for r = 1 mm, but similar to the results shown in figure 4 for the Medtronic lead, the orientation selectivity was muted farther from the DBS lead. Another feature of DBS leads with segmented rows was enhanced selectivity for nearby axons oriented perpendicular to the lead using anodes and cathodes in the same row (figures 6(f) and 7(f)). This configuration showed greater orientation-selective gains than the multi-cathode paradigm using columnar cylindrical electrodes (figure 4(f)).

Comparisons between monopolar, bipolar, and MICC stimulation
The inclusion of local anodes can be used to reduce current spread to neighboring regions beyond the target volume of interest around a DBS lead. This study showed that bipolar configurations could be leveraged to generate orientation selective modulation of axonal fibers, especially in cases in which an axon projects parallel to the electric field and between the active contacts. This selectivity gain was reduced directly over the electrodes suggesting that axonal trajectories and branching patterns at the mm to sub-mm scale near DBS leads may be critical to optimizing spatiotemporal stimulation patterns to affect the desired axonal pathway modulation. While bipolar configurations exhibited higher orientation selectivity gains than did monopolar configurations, the angles exhibiting orientation selectivity in both single source cases were limited by the number, location, size, and spacing of the electrode contacts on the DBS lead. The models indicated that using a single source with two or more anodes or cathodes could improve the angular resolution of the orientation selectivity; however, the amount of current delivered through each contact in such cases would depend on the relative impedance between cathodes and between anodes. One approach to circumvent this issue is the use of multiple independent current sources. The models showed how MICC configurations that incorporated cathodes and anodes with unequal splitting of current enabled flexible steering of the electric field across a much broader range of axon orientations. Further, using MICC with segmented DBS leads, one effectively has control to achieve orientation selectivity for any axon angle between −90 to 90° relative to the lead.

Clinical significance
This study suggests an alternative approach to consider when programming DBS systems, with particular relevance to optimizing therapeutic windows for providing effective treatment without eliciting stimulation-induced side effects. This approach is based on maximizing the alignment of the induced electric field with the target pathway(s) of interest, while minimizing alignment with axonal pathways implicated in side-effects of DBS. For example, in a target such as the STN, the HDP is approximately perpendicular to typical DBS lead implant trajectories while main IC fibers are approximately parallel to the lead. DBS leads with electrodes segmented around the lead (such as the Abbott 6672 or Boston Scientific 2202 leads) could be used to better align the primary direction of the electric field to achieve robust HDP activation. As shown in this study, when attempting to activate the HDP with a four-contact lead with cylindrical electrodes, the primary direction of the electric field can be aligned indirectly with using multi-cathode configurations. This was shown in retrospective analysis of imaging data from three STN-DBS patients in which multi-cathode configurations consistently activated the HDP to a larger degree before activating any axon within the IC. Such orientation-selective methods may be useful in cases in which one wishes to widen the therapeutic window between therapy and side effect thresholds. To our knowledge a comparison between bipolar and multi-cathode configurations has not yet been fully characterized leveraging prospective subject-specific computational models.
Orientation-selective DBS programming methods may also be useful to probe how modulating individual pathways more selectively in and around a target region affects therapeutic outcomes and stimulus thresholds for eliciting side effects. One pathway implicated in the physiological mechanisms of DBS therapy for PD is direct modulation of STN efferents to the pallidum [47]. In this case, the axonal trajectories are similar in orientation to those of the HDP, suggesting that the results from the three STN-DBS patients could extrapolate to modulation of the STN output as well. However, there are other pathways adjacent to the STN that have different orientations from those of the HDP, IC, and STN efferent axons and that may be of therapeutic relevance to treatment of PD. These include the pallidofugal pathway within the thalamic fasciculus [48], the cerebellothalamic pathway [49], and fibers of passage coursing to and from brainstem nuclei.
There are a number of other clinical indications and DBS targets in which orientation-selective approaches could be useful. For instance, DBS targeting the internal GPi for PD [50] and dystonia [51] involves implanting a lead along an oblique trajectory such that the lead is roughly tangential to afferents to GPi and efferents extending from GPi, and roughly parallel to the IC. Targeting GPi-related pathways [21,42], one might consider using multi-cathode stimulation, which has been noted in the literature [52]. Another example in which orientation-selective stimulation approaches would be useful is the ventral intermediate nucleus of thalamus for treating Essential Tremor [53], where the cerebellothalamic and corticothalamic/reticular nucleus fibers are oriented approximately tangential to one another, with the former fibers coursing roughly parallel to the lead [54][55][56]. Further, if one considers that several groups are now advocating targeting the posterior subthalamic area for Essential Tremor, one might consider using a DBS lead design with multiple rows of segmented contacts, especially the distal row [22,57] in order to align the electric field to target the superior cerebellar peduncle more selectively, but also to avoid strong electric field orientations for more distant fibers that can induce side effects including paresthesias when stimulated.
The orientation-selective approach could also be useful for DBS targets embedded in white matter with multiple crossing fiber pathways. Retrospective modeling studies of responders to DBS for treatment-resistant major depression, for example, indicated that three pathways (forceps minor, uncinate fasciculus, and the cingulum bundle) as well as subcortical nuclei very likely need to be modulated in order to achieve a therapeutic response [5,58]. These pathways have different axonal orientations, with the cingulum bundle roughly orthogonal to the forceps minor and uncinate fasciculus. Orientationselective activation in this case may improve activation of each of these pathways without modulating other pathways. Additionally, the orientation-selective approach could be leveraged to more precisely identify how modulation of each of these pathways contributes to the reported therapeutic outcomes of the responder populations using an intrasubject experimental design.

Study limitations
This computational study contains several limitations that should be considered. First, conceptual models were simplified using a homogeneous, isotropic sphere for the brain tissue that did not take into account the complexity of conductivity and geometric variability in the brain. In addition, axons were modeled as straight lines with no variability in the center node location. Activation is more likely at bends in axons [43], suggesting that the modeled axons likely under-predict activation. There is a growing need to map with finer resolution the local axonal trajectories and branching patterns near DBS targets. These simplifications were included in the models with the purpose of gaining intuition and confirming the theory of aligning the primary direction of the electric field with clinical DBS lead designs. Additional complexity was added with the patient-specific anatomical models, which included the complexity of inhomogeneous, anisotropic tissue and patientspecific axon trajectories. Yet, there remain several unknowns with patient-specific models of DBS. For one, the diameter of axons in each pathway were constant; histological studies have noted that axonal diameters can have a wide distribution within a given pathway [30,31], but these diameter distributions have not been fully characterized for each DBS target and pathways within each target. Additionally, the study focused on myelinated axons, while many targets contain both myelinated and unmyelinated axons. Activation of unmyelinated axons often requires larger stimulus intensities, and thus different orientation selectivity may be observed between unmyelinated and myelinated axons [42]. Future studies are needed to determine how myelination, axon diameter, and electrode-tissue interface complexity affect both orientationselectivity and orientation-selective gain.
These patient-specific models were fairly complex, but simplified the capacitive effects by using a recorded waveform rather than a discrete Fourier transform approach [59,60]. The fine precision, and its clinical implications, heavily depends on the accuracy of the patient-specific model as well as the precise localization of the DBS lead in relationship to the anatomical features. While such imaging capabilities exist in specific academic institutions that specialize in high-field imaging research, the applicability of such technology is still limited for the general community. Finally, as with most computational modeling studies, when incorporating the simplifications and sources of errors in these models, the effect size can be lower than predicted and clinical validation is needed to confirm the novel results and predictions found here.

Conclusions
This study sought to characterize the degree of orientationselective activation that can be achieved using clinically available DBS systems. It was found that DBS leads composed of cylindrical contacts could be used to selectively activate axons parallel and perpendicular to the lead shaft. Segmented leads provided further selectivity for axons approximately 0°, ±22.5°, ±45°, ±67.5°, and 90° relative to the lead shaft using a single current source. Multiple current sources combined with a segmented lead provided selectivity that could be shifted to any angle between −90° and 90°. Patient-specific results showed that employing multiple cathodes among cylindrical contacts provided greater activation of a target oriented perpendicular to a DBS lead while avoiding adjacent tracts coursing parallel to the lead shaft. The results found here are poised for evaluation in both preclinical and clinical studies.