Theoretical Optimization of Stimulation Strategies for a Directionally Segmented Deep Brain Stimulation Electrode Array

Abstract

Programming deep brain stimulation (DBS) systems currently involves a clinician manually sweeping through a range of stimulus parameter settings to identify the setting that delivers the most robust therapy for a patient. With the advent of DBS arrays with a higher number and density of electrodes, this trial and error process becomes unmanageable in a clinical setting. This study developed a computationally efficient, model-based algorithm to estimate an electrode configuration that will most strongly activate tissue within a volume of interest. The cerebellar-receiving area of motor thalamus, the target for treating essential tremor with DBS, was rendered from imaging data and discretized into grid points aligned in approximate afferent and efferent axonal pathway orientations. A finite-element model (FEM) was constructed to simulate the volumetric tissue voltage during DBS. We leveraged the principle of voltage superposition to formulate a convex optimization-based approach to maximize activating function (AF) values at each grid point (via three different criteria), hence increasing the overall probability of action potential initiation and neuronal entrainment within the target volume. For both efferent and afferent pathways, this approach achieved global optima within several seconds. The optimal electrode configuration and resulting AF values differed across each optimization criteria and between axonal orientations. This approach only required a set of FEM simulations equal to the number of DBS array electrodes, and could readily accommodate anisotropic-inhomogeneous tissue conductances or other axonal orientations. The algorithm provides an efficient, flexible determination of optimal electrode configurations for programming DBS arrays.

Fig. 1

Three-dimensional thalamic nuclei reconstructions were generated from (A) NHP susceptibility-weighted imaging and (B) warped brain atlas overlays. (C) Sagittal and (D) coronal view of the reconstructed VPLo and VPLc with the implanted DBSA.

Fig. 2

Procedural flowchart for this study. Three different optimization criteria (MD, QP, LP) were considered separately.

Fig. 3

Discretization of thalamic volumes. The top row (A,B) shows the discretization process for the efferent direction while the bottom row (C,D) shows the process in the afferent direction. Both the VPLo and VPLc are shown in coronal orientation. In either case, the left image shows grid points arranged in serial layers spanning either volumes. The red arrow indicates the orientation of the neuronal processes. The right image shows only those grid points that fall within the volume (red points). The internodal distance between successive layers of grid points are shown. (A) and (C): VPLo. (B) and (D): VPLc.

Fig. 4

Algorithm-generated electrode configurations for the thalamic efferent pathway approximations. The left, middle, and right columns show outcomes from the MD, QP and LP optimization criteria, respectively. Active contacts (> 1 μA) in each case are shown in red (top row) with the precise amount of current calculated by the algorithm shown in indexed colors (second row). Axial views (third row) and oblique views (fourth row) of VPLo and VPLc are shown in the context of the DBSA with active contacts shown in red.

Fig. 5

Algorithm-generated electrode configurations for the thalamic afferent pathway approximations. Labeling is identical to that described for Figure 4.

Fig. 6

Comparison of Max Curve (red) to solutions obtained by MD, QP, and LP for efferent (top) and afferent (bottom) data. Grid points in the region of interest are sorted based on their maximum achievable values. For each optimization criterion, the actual value at each grid point is plotted underneath its maximum possible value. Therefore, the closeness of the grid points to the Max Curve is a measure of optimization performance. (A–C) and (G–I): AF values at all grid points are presented in units of V/mm². (D–F) and (J–L): Grid points with AF values less than −0.01 V/mm² were omitted in these plots. The remaining AF values were made positive by adding 0.01 and the natural logarithm of the resulting values was computed.

Fig. 7

AF values resulting from DBSA stimulation using algorithm-generated electrode configurations. (A–C): Efferent data, axial view of the AF values resulting from stimulation configurations generated using the MD, QP, and LP criteria, respectively. (D–F): Coronal view of the AF values shown in (ac). (G–I): Afferent data, axial view of the AF values resulting from stimulation configurations generated using the MD, QP, and LP criteria, respectively. (J–L): Coronal view of the AF values shown in (j–l). For visualization purposes, all AF values greater than −0.01V/mm² were made positive by adding 0.01 and the natural logarithm of the resulting values was computed. The logarithm values were used as indexed colors. The color bar in this figure ranges from −6 to −3 in logarithm values. Points with values outside of this range were directly assigned the values of −6 or −3. Points with AF values less than −0.01V/mm² were assigned logarithm values of −6. Refer to Figure 1(D) for borders between VPLo and VPLc in the coronal view. Refer to Figure 4 or 5 for the border in the axial view.

Fig. 8

Performance comparison of MD, QP, and LP to one million random electrode configurations (gray normalized histogram), in efferent and afferent data. The three measures considered (Max Deviation, Sum of Squares, Sum) each indicate deviation from maximum possible activating function values. Thus, lower values correspond to better performance.

Fig. 9

Comparison of runtime and sampling robustness for MD, QP, and LP. The original number of grid points for efferent data was 27,173. For each level of sampling (½, ¼, etc.) a random subset of the original points were selected. (A, D, G) Average height of active electrodes with respect to multiple sampled subsets. Standard deviation of height progressively increases as smaller subsets of original points are used. (B, E, H) Average angular direction of active electrodes. Similarly to average height, standard deviation increases toward random distribution. (C, F, I) Runtime and logarithm of runtime with respect to sampling.