Adaptive Estimation of the Neural Activation Extent in Computational Volume Conductor Models of Deep Brain Stimulation

TL;DR

This paper introduces an adaptive, open-source Python method for efficiently estimating neural activation extents in deep brain stimulation models, accounting for patient-specific tissue properties and reducing computational costs.

Contribution

It presents a novel adaptive scheme embedded in an open-source environment for estimating neural activation in DBS models, incorporating tissue properties and axon parameters.

Findings

Neural activation extent varies with dielectric tissue properties.

Threshold-based estimates significantly reduce computational time.

Patient-specific tissue properties influence activation extents.

Abstract

Objective: The aim of this study is to propose an adaptive scheme embedded into an open-source environment for the estimation of the neural activation extent during deep brain stimulation and to investigate the feasibility of approximating the neural activation extent by thresholds of the field solution. Methods: Open-source solutions for solving the field equation in volume conductor models of deep brain stimulation and computing the neural activation are embedded into a Python package to estimate the neural activation dependent on the dielectric tissue properties and axon parameters by employing a spatially adaptive scheme. Feasibility of the approximation of the neural activation extent by field thresholds is investigated to further reduce the computational expense. Results: The varying extents of neural activation for different patient-specific dielectric properties were estimated with the adaptive scheme. The results revealed the strong influence of the dielectric properties of the encapsulation layer in the acute and chronic phase after surgery. The computational time required to determine the neural activation extent in each studied model case was substantially reduced. Conclusion: The neural activation extent is altered by patient-specific parameters. Threshold values of the electric potential and electric field norm facilitate a computationally efficient method to estimate the neural activation extent. Significance: The presented adaptive scheme is able to robustly determine neural activation extents and field threshold estimates for varying dielectric tissue properties and axon diameters while reducing substantially the computational expense.

Figures

Tables

Table 1. TABLE I: Conductivity of the volume conductor compartments for the different study cases. The values were determined by using the tissue parameters from [ 17 ] .

Study Case	Conductivity [${\mathrm{Sm}}^{-1}$]
	encapsulation	brain	subthalamic
	layer	tissue	nucleus (STN)
Model 1 (homogeneous)	0.064	0.064	0.064
Model 2 (acute phase)	2.000	0.064	0.064
Model 3 (chronic phase)	0.032	0.064	0.064
Model 4 (with STN)	0.032	0.064	0.103

Table 2. TABLE II: Stimulation amplitudes and normalized activation thresholds for approximating the neural activation extent in a distance of 2   mm times 2 mm 2\text{\,}\mathrm{m}\mathrm{m} to 4   mm times 4 mm 4\text{\,}\mathrm{m}\mathrm{m} radial to the active electrode contact.

Study Case	Axon	Amplitude	Maximum
	diameter	range	normalized threshold
	in $\mu m$	in -$\mathrm{mA}$	for $\phi$	for ${\Vert E\Vert }_2$
Model 1	2.0	[1.2, 8.3]	3.37	1.70
	3.0	[0.7, 4.1]	2.91	1.49
	5.7	[0.4, 2.2]	2.64	1.30
	7.3	[0.3, 1.3]	2.18	1.11
	8.7	[0.2, 0.9]	2.13	1.02
	10.0	[0.2, 0.8]	1.96	0.96
Model 2	2.0	[2.6, 13.0]	2.98	1.74
	3.0	[1.4, 6.3]	2.69	1.57
	5.7	[0.9, 3.4]	2.29	1.35
	7.3	[0.6, 2.0]	2.02	1.19
	8.7	[0,4, 1.4]	2.01	1.11
	10.0	[0.4, 1.1]	1.72	1.04
Model 3	2.0	[1.2, 8.1]	3.30	1.67
	3.0	[0.6, 4.0]	3.14	1.54
	5.7	[0.4, 2.2]	2.63	1.29
	7.3	[0.3, 1.3]	2.18	1.11
	8.7	[0.2, 0.9]	2.12	1.01
	10.0	[0.2, 0.8]	1.95	0.96
Model 4	2.0	[1.8, 6.0]	2.11	1.47
	3.0	[0.9, 3.3]	2.18	1.67
	5.7	[0.6, 2.0]	1.97	1.29
	7.3	[0.4, 1.2]	1.79	1.18
	8.7	[0.3, 0.9]	1.73	1.11
	10.0	[0.2, 0.7]	1.83	1.00

Table 3. TABLE III: Number of required axon models, speed up and computation time for the different study cases (models). The number of axons and the computation times are determined for the adaptive scheme, the non-adaptive scheme, and for approximating the neuronal activation extent (VTA) by field threshold values for axons with 5.7   µ ​ m times 5.7 µ m 5.7\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m} fiber diameter. The computation time is measured on a 12 × 12\times 2.4   GHz times 2.4 GHz 2.4\text{\,}\mathrm{G}\mathrm{H}\mathrm{z} , 48   GB times 48 GB 48\text{\,}\mathrm{G}\mathrm{B} workstation. † rotational symmetry was used for computing the neural activation extent. ∗ Values determined by mapping number of axons and computation time to rectangular grids of same extent.

Study Case	Number of Axons			Speed up	Computation time
	VTA	VTA	VTA	VTA	field	VTA	VTA	VTA
	non-adaptive∗	adaptive	field threshold	adaptive	solution	non-adaptive∗	adaptive	field threshold
Model 1†	$152$	$94$	$8$	$54\%$	$4.9\min$	$36.7\min$	$23.8\min$	$3.5\min$
Model 2†	$288$	$173$	$8$	$59\%$	$4.7\min$	$71.3\min$	$44.8\min$	$3.4\min$
Model 3†	$136$	$94$	$8$	$38\%$	$4.8\min$	$35.8\min$	$25.9\min$	$3.4\min$
Model 4	$4,896$	$2,813$	$288$	$66\%$	$5.0\min$	$1,155.8$$\min$	$696.3\min$	$67.0\min$

Equations

$$\nabla \cdot \left[\sigma (\mathbf{r})\nabla (\phi (\mathbf{r}))\right]=0,\mathbf{r}\in \Omega$$

$${\sigma }_c(\omega ,\mathbf{r}):=\sigma (\omega ,\mathbf{r})+\mathsf{j}\omega {ϵ}_0{ϵ}_r(\omega ,\mathbf{r})$$

$${\sigma }_c(\omega )={ϵ}_{\infty }+\frac{{\sigma }_{\mathsf{ion}}}{\mathsf{j}\omega {ϵ}_0}+\sum _{i=1}^4\frac{\Delta {ϵ}_i}{1+(\mathsf{j}\omega {\tau }_i{)}^{1-{\alpha }_i}}$$

$$P=\underset{\Omega }{\int }⟨\sigma (\mathbf{r})\nabla \phi (\mathbf{r}),\nabla \phi (\mathbf{r})⟩\mathsf{d}x$$

$$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial \phi }{\partial r}\right)+\frac{1}{r^2\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial \phi }{\partial \theta }\right)=0$$

$$\frac{1}{\sin \theta }\frac{d}{d\theta }\left(\sin \theta \frac{d\Theta (\theta )}{d\theta }\right)+\lambda \Theta (\theta )$$

$$\frac{d}{dr}\left(r^2\frac{dR(r)}{r}\right)-\lambda R(r)$$

$${\Theta }_n(\theta )=D_n{P}_n(\cos \theta ),P_n(x)=\frac{1}{2^{n}n!}\frac{d^n}{d{x}^n}\left((x^2-1{)}^n\right)$$

$$R_n(r)=C_{1,n}{r}^n+C_{2,n}{r}^{-(n+1)}$$

$$\phi =\{\begin{array}{ll}\sum _{n=0}^{\infty }{A}_n{r}^n{P}_n(x)& ,0\le r\le R_i\\ \sum _{n=0}^{\infty }\left(B_n{r}^n+C_n{r}^{-(n+1)}\right)P_n(x)& ,R_i<r\le R_e\\ -E_{0}rx+\sum _{n=0}^{\infty }{D}_n{r}^{-(n+1)}{P}_n(x)& ,R_e<r\end{array}$$

$${\phi }_1-{\phi }_2$$

$$\mathbf{n}({\mathbf{J}}_1-{\mathbf{J}}_2)$$

$$\left(\begin{array}{cccc}{R}_i^3& -R_i^3& -1& 0\\ {\sigma }_{i,i}& -{\sigma }_{s,i}& 2{\sigma }_s& 0\\ 0& R_e^3& 1& -1\\ 0& {\sigma }_{s,e}& -2{\sigma }_s& 2{\sigma }_e\end{array}\right)\cdot \left(\begin{array}{c}{A}_1\\ B_1\\ C_1\\ D_1\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ -E_{0,e}\\ -{\sigma }_e{E}_{0,e}\end{array}\right)$$

$$\phi (r,\theta )=\{\begin{array}{ll}{A}_{1}r\cos \theta & ,0\le r\le R_i\\ \left(B_{1}r+C_1{r}^{-2}\right)\cos \theta & ,R_i<r\le R_e\\ \left(D_1{r}^{-2}-E_{0}r\right)\cos \theta & ,R_e<r\text{ .}\end{array}$$

Full text

Deep brain stimulation (DBS) is a widely employed effective procedure to treat symptoms of motor disorders such as Parkinson’s disease (PD), essential tremor and dystonia [1, 2] and consists in the implantation of an electrode lead into deep brain target areas. A common target in DBS is the subthalamic nucleus (STN), which constitutes a preferred target for the treatment of PD. The STN consists of different functional zones, which are classified into limbic, associative, and sensorimotor zones [3], from which electrical stimulation of the sensorimotor zone is mostly associated with the relief of motor symptoms of PD [4]. Due to patient-specific parameters, such as brain structure anatomy, dielectric tissue properties, electrode location, and severity of symptoms, the adjustment and optimization of stimulation parameters during and after surgery can be rather time consuming and connected to additional costs. Computational models provide a possibility to estimate the stimulation impact by determining activated areas in the deep brain based on the given patient-specific parameters. The extent of neural activation, or the volume of tissue activated (VTA), is a common computational modeling approach to estimate the size of the activated tissue during DBS and has been applied in various computational studies in this area including homogeneous [5], rotationally symmetric [6], heterogeneous [7, 8], and anisotropic volume conductor models [9] of the human brain and deep brain target areas. In general, the approach is based on positioning a number of models of mammalian nerve fibers (axon models) in a grid located in a plane perpendicular to the electrode lead. For each axon model the computational goal of finding the minimum stimulation amplitude required to activate the axon is solved. From the resulting threshold values at the grid points, a threshold isoline for a given stimulation amplitude is determined. This procedure is repeated for multiple planes rotated around the electrode lead. In case of a rotationally symmetric field model, it is sufficient to compute the threshold isoline in one plane and revolve the solution around the electrode lead. The resulting threshold isolines then provide the measure for computing the VTA. The drawback of the method with respect to computational ressources and adaptivity is that the location of the axon nodes to include a range of desired stimulation amplitudes has to be available prior to the simulation, which involves several pre-simulation runs, which often are carried out manually. The field solution in the target area is commonly computed by creating a volume conductor model of the DBS electrode and the surrounding tissue. For solving the governing equations, which are typically the stationary current field or electro-quasistatic equation for DBS applications [7, 10], often commercial software solutions, such as COMSOL Multiphysics® (http://www.comsol.com) are used [5, 7, 8, 11, 12, 13], while, to our knowledge, no studies employing open-source solutions on the field model have been published in scientific journals yet. Regarding the coupling of the neuronal activation in axon models and the extracellular field distribution, a Python package with the purpose to compute the local field potentials for a given axon distribution of defined activation was presented in [14]. To date, no open-source solution to model the field distribution during DBS and to estimate the resulting neural activation extent exists.

The computation of the VTA is a computationally demanding task. In a previous computational study investigating the relationship between the neural activation and the field solution, an approximation of the extent of neural activation by threshold values of the field solution and their derivates have been investigated [13]. The results suggested that electric field norm thresholds are a good estimator for the extent of neural activation. In such a case only a small neural activation extent has to be computed to get the initial field norm threshold estimate, from which the neural activation extents for varying stimulation amplitudes can be derived. Especially for large diameter axon fibers, the electric field norm constituted a good estimate. The relationship between electric field norm and neural activation was determined by positioning axons normal to the electrode lead along a line originating at the active electrode contact center in a homogeneous volume conductor model for DBS. The threshold value of the electric field norm was equivalent to its value at the maximum neural activation distance for a given stimulation amplitude. The assumption of this approach is that the shape of the neural activation extent is spherical, which is true for a homogeneous volume conductor model and a spherical or point source. For heterogeneous volume conductor models and DBS electrode geometries, the shape deviates from the spherical shape.

The goal of this study is to investigate the feasibility of approximating the neural activation extent by the field solution of a heterogeneous volume conductor model for DBS incorporating a DBS electrode, encapsulation layer, brain tissue, and the STN as target area. To automate the task of determining the neural activation extent and further reduce the computational demand, an algorithm is proposed which determines adaptively the location of axons being activated within a defined stimulation amplitude range or distance range. Besides dropping the need for manually determining the number of axons in the target area for a given stimulation amplitude range, the approach reduces the computational expense by omitting the computation of activation for axons which are located outside the activation volume. The model pipeline for the computation of the field solution and the neural activation is implemented in a Python package and embedding open-source tools for the model generation, meshing, and solving. The field solution and the neural activation are validated using analytical models as well as reference data published in literature. The Python package is designed modular, which allows to interchange field as well as neuron models and to adjust model parameters accordingly. The Python package, as well as the code to replicate the data and figures of this study are made available open-source111https://bitbucket.org/ChrSchmidt83/fanpy/get/fanpy-1.2.zip.

I Methods

I-A Model geometry

Following the approach in [15], the model geometry consists of a DBS electrode model located in a bounding box comprising the different tissue compartments. The geometry of the DBS electrode represents a Medtronic lead model 3387 (Medtronic, Inc., Minneapolis, MN). To account for the inflammatory response of the body tissue to the electrode implant, an encapsulation layer with a thickness of $0.2\text{\,}\mathrm{m}\mathrm{m}$ was incorporated around the electrode body. The bounding box size was determined by an edge length of $100\text{\,}\mathrm{m}\mathrm{m}$. The geometry model of the STN based on a functional zones atlas [16] was generated by creating a surface model out of the right STN threshold maps of the atlas using the open-source software platform 3D Slicer (https://www.slicer.org/, version 4.2.1). The whole geometry model was generated with the open-source software SALOME (http://www.salome-platform.org/, version 7.8.0). After creating and merging of the different compartments, the STN surface model was converted to a solid and positioned in the geometry with the second electrode contact of the DBS electrode located in the sensorimotor zone (Fig. 1).

I-B Manual mesh refinement and subdomain generation

The meshing of the computational domain is carried out with the open-source mesh generator Gmsh (http://gmsh.info, version 2.10.1). Therefore, the model geometry is exported from SALOME in the brep format and loaded into gmsh. Since the geometry contains entities of varying scales (Fig. 1), the mesh was manually refined at the surfaces of the entities by specifying a characteristic length of the finite elements. Based on values for the manual mesh refinement for DBS volume conductor models of the human brain [9], which provided a sufficient refinement of the computational domain, the characteristic length was set to $0.1\text{\,}\mathrm{m}\mathrm{m}$ for the electrode lead, electrode contacts, and encapsulation layer, to $0.2\text{\,}\mathrm{m}\mathrm{m}$ for the STN, and $5.0\text{\,}\mathrm{m}\mathrm{m}$ for the bounding box. Additionally, a refined cubical mesh region with an edge length of $20\text{\,}\mathrm{m}\mathrm{m}$ and a characteristic length of $0.5\text{\,}\mathrm{m}\mathrm{m}$ was defined in the target area. The final mesh contained approximately $240{,}000$ vertices and $1.5$ million cells (Fig. 2). In order to assign material properties and boundary conditions to the model compartments and surfaces, the subdomains of the geometry model were assigned and grouped to physical volumes and surfaces. The information on the manual mesh refinement and the defined physical volumes and surfaces is stored in a geo file, from which automatically the mesh can be generated.

I-C Incorporating the model into FEniCS

The mesh generated with Gmsh is converted into a format readable by the open-source simulation software FEniCS (https://fenicsproject.org, version 2016.2.0) by using the command line tool dolfin-convert. To faciliate the interchange of different model geometries and meshes, information on the model’s subdomains, boundaries, and default boundary conditions and material properties are defined in an xml file. The FieldModel module of the designed Python package loads the definitions in the model xml file, checks the information for consistency, generates the mesh and applies the material properties and boundary conditions. While the conductivity of the electrode lead (insulation) and the electrode contacts (platinum-iridium) are kept constant with a value of $1\text{\cdot}{10}^{-7}\text{\,}\mathrm{S}\mathrm{m}^{-1}$ for the electrode lead, and $1\text{\cdot}{10}^{7}\text{\,}\mathrm{S}\mathrm{m}^{-1}$ for the electrode contacts, the conductivity values for the encapsulation layer, brain tissue, and STN varied for the different study cases. The field equation is determined by a stationary current field problem, which is commonly applied in various computational modeling studies for DBS [13, 15]:

[TABLE]

with the conductivity $\sigma$, the electric potential $\varphi$ and the computational domain $\Omega$. If the capacitive and dispersive properties of human tissue are taken into account, the field equation (1) has to be reformulated for complex materials with

[TABLE]

with the imaginary unit $\mathsf{j}$, the angular frequency $\omega$, the electric field constant $\epsilon_{0}$, and the relative electric permittivity $\epsilon_{\mathsf{r}}$, which placed in equation (1) resembles a quasistatic field problem. The field equations are solved with the finite element method by formulating the variational problem within FEniCS. Dirichlet boundary conditions are applied to the surface of the second electrode contact, located within the STN (Fig. 1), with a unit value of $1.0\text{\,}\mathrm{V}$, and the exterior boundary of the bounding box with a value of $0.0\text{\,}\mathrm{V}$, serving as ground. The resulting linear system of equations was solved for quadratic nodal basis functions using the generalized minimal residual method with a relative tolerance of $1\cdot 10^{-6}$ and an absolute tolerance of $1\cdot 10^{-7}$, employing an algebraic multigrid preconditioner in case of the stationary current field problem. It was ensured that the deviation in the electric potential and electric field norm in the prescribed activation distances (see section I-H) as well as the impedance of the model was below $1\text{\,}\mathrm{\char 37\relax}$ if cubic ansatz functions were employed. Plots of the field distribution were visualized with the open-source visualization application Paraview (http://www.paraview.org/, version=5.0.1).

I-D Electrical properties of human tissue

The conductivity and relative permittivity of biological tissue show a frequency dependence which can be described by different dispersion regions and parametrized by assembled Cole-Cole equations [17] representing a complex conductivity

[TABLE]

with the static ionic conductivity $\sigma_{\mathsf{ion}}$, the relaxation time constants $\tau_{i}$, the dispersion constant $\alpha_{i}\in[0,1]$, and the relative permittivity at high frequency $\epsilon_{\infty}$ as well as the difference of the low and high frequency relative permittivity $\Delta\epsilon_{i}$. The tissue model parameters were taken for white matter, grey matter, and cerebrospinal fluid from [17]. The electrical tissue properties of the encapsulation layer vary over time from an acute phase immediately after surgery to a chronic phase after some weeks due to cell growth in the layer [18]. The acute phase was modeled by the dielectric properties of cerebrospinal fluid [7], while the chronic phase was modeled by dividing the values for white matter by a factor of 2 [15].

I-E Voltage-controlled and current-controlled stimulation

The time-dependent electrical potential in the target area for a given stimulation signal was determined by forming the outer product of the field solution for a unit voltage set to the active second electrode contact (Fig. 1) and the voltage- or current-controlled stimulation signal. The stimulation signals commonly applied in DBS therapy in humans consist of a monophasic square-wave signal with pulse durations in the range of $60\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$ - $100\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$ and a repetition frequency in the range of $130\text{\,}\mathrm{H}\mathrm{z}$ - $150\text{\,}\mathrm{H}\mathrm{z}$ [15, 13, 7, 19]. To avoid charge accumulation in the tissue, the monophasic stimulation pulse is often followed by a reversed charge-balancing pulse of substantially smaller amplitude compared to the active stimulation pulse. Considering that the activation of a neuron is mainly influenced by the amplitude of the stimulation pulse [7], the reversed charge-balancing pulse is not considered in this study. While the time-dependent electrical potential in the target area for voltage-controlled stimulation is provided by the outer product of the field solution for a unit voltage and the voltage-controlled stimulation signal, the time-dependent electrical potential for current-controlled stimulation requires an additional scaling of the field solution by the electrode impedance, which is computed by dividing the square of the unit voltage by the electric power $P$ of the field model.

[TABLE]

with the inner product $\langle,\rangle$, corresponding to a scaling of the unit voltage at the active electrode contact boundary condition by a factor equal to a unit current flowing through its surface.

I-F Neuronal activation model

The neural activation model is based on a myelinated axon cable model, which includes $21$ nodes of Ranvier, paranodal and internodal segmenets as well as the myelin sheath [20], following the assupmtion that activation occurs along the axon [12]. Based on the model parameters given in [20], the model is parametrized with respect to the fiber diameter comprising nine distinct diameters between $5.7\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ and $16.0\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, which were extented by the parameters for $2.0\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ and $3.0\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ fiber diameter taken from [21]. The model is implemented in the open-source simulation environment NEURON (https://www.neuron.yale.edu/neuron/, version=7.4)222https://senselab.med.yale.edu/modeldb/showModel.cshtml?model=3810 with the time-dependent electrical potential for an applied stimulation signal at the location of each axon compartment applied to its extracellular mechanism node neglecting any axon contribution to the extracellular field distribution [15]. The axon activation is determined by solving the linear system of differential equations resulting from the membrane dynamics with the backward Euler method for a time step of $5\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$ [15]. An axon was considered to be activated when the inner potential at the exterior nodes of Ranvier of the model obtained a threshold value of larger than $0\text{\,}\mathrm{m}\mathrm{V}$, representing a generated spike as result of the stimulation pulse, in a $1$-to-$1$ ratio for $10$ delivered stimulation pulses.

I-G Adaptive estimation of the neural activation extent

The proposed algorithm determines adaptively the required location of axons in the target area for a given range of stimulation amplitudes, which ensures that the extent of neural activation for any stimulation amplitude within the given range can be computed from the determined axon locations. The required axon locations are determined in cutting planes, which are located around the electrode lead. First, a seed point, located at a distance of $0.85\text{\,}\mathrm{m}\mathrm{m}$ to the active electrode contact’s center, which corresponds to the extent of the electrode with the encapsulation thickness, is placed in such a plane. Next, the algorithm determines whether the axon, which is positioned perpendicular and centered to the electrode lead at the seed point location, is activated for the minimal given stimulation amplitude. If an activation was recorded, the algorithm continues by placing further axons around the activated axon, with their center node location positioned radial ($\varrho+\Delta s$) and parallel ($z+\Delta s,z-\Delta s$) to the electrode lead with a step size $\Delta s$, and determining their activation. If the axon was not activated, the algorithm stops for this location. The procedure is continued until an inactivated hull of axons is determined, which is achieved if in each line radial to the electrode the axon furthest away to the electrode is not activated by the stimulation. This inactivated hull is then used as seed points for determining the axon locations for the maximum stimulation amplitude within the given range by applying the same algorithm. Finally, the interior points, which are located inside the activated hull of axons for the minimum given stimulation amplitude are removed from the set of locations, resulting in a shell of axon locations for the given stimulation amplitude range. Determining the minimally required stimulation amplitude to elicit an action potential in the axons at these locations represents a root-finding problem, which is solved using the bisection method (binary search method) with a tolerance of $1\cdot 10^{-6}$ (Fig. 3). The algorithm is carried out for $n_{\alpha}=\lceil 360/\Delta\alpha\rceil$ planes around the electrode with the rotational degree step size $\Delta\alpha$. The subsequent computation of the neural activation for each axon model introduced by the algorithm as well as finding the minimally required stimulation amplitude to activate the axon model are carried out by the model pipeline in parallel with worker threads adding and withdrawing axon points to a pool (Python Queue) of axons. In order to prevent the placement of new axon points at the same location to the pool, a new axon point is only added if an axon point with the same location was not added before by another worker thread.

I-H Approximating the neural activation by field thresholds

The proposed adaptive algorithm is used to assess the feasibility of approximating the neural activation extent for various stimulation amplitudes by threshold values of the electrical potential and the electric field norm. A current-controlled stimulation is applied, which is gaining increased interest in clinical application due to the reduced side effects and sensitivity to inter-individual variabilities [22, 23]. Different values for the dielectric properties of the volume conductor model compartments are employed to investigate the approximation quality for different cases of tissue heterogeneity in the target area. The dielectric properties are obtained from equation (3) using the parameters for white matter, grey matter, and cerebrospinal fluid from [17] at a frequency of $2\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$, which constitutes a good approximation of the dispersive nature of the tissue properties for common DBS signals [8]. The corresponding conductivity values are approximately $0.064\text{\,}\mathrm{S}\mathrm{m}^{-1}$ for white matter, $0.103\text{\,}\mathrm{S}\mathrm{m}^{-1}$ for grey matter, and $2.000\text{\,}\mathrm{S}\mathrm{m}^{-1}$ for cerebrospinal fluid. The cases comprise a homogeneous model (Model 1) with the conductivity of the encapsulation layer, the brain tissue, and STN set to the value for white matter, a model with a high conductive encapsulation layer (Model 2) set to the value of cerebrospinal fluid [7] representing the acute phase, as well as with a low conductive encapsulation layer (Model 3), set to half the value of white matter [12] representing the chronic phase, and a heterogeneous model (Model 4) with a low conductive encapsulation layer, brain tissue set to the value of white matter, and the STN set to the value of grey matter. Since the varying dielectric tissue properties result in a variation of the conduction in the volume conductor model, the extent of neural activation is dependent on the tissue properties [11]. Therefore, DBS with the same stimulation amplitude results in different sizes of the VTA depending on the tissue properties in the models. In order to ensure a stimulation of the target region between a distance of $2.0\text{\,}\mathrm{m}\mathrm{m}$ (activation of sensorimotor functional zone) and $4.0\text{\,}\mathrm{m}\mathrm{m}$ (activation of larger parts of the STN) from the electrode center, the required stimulation amplitude to activate a homogeneous volume within this minimum and maximum distance was computed for each model in advance. The determinted stimulation amplitude range is then prescribed to the adaptive algorithm to constrain the extent of the estimated neural activation. Except for the latter model, the neural activation extent is computed exploiting rotational symmetry by computing the extent in a reference plane and revolving the solution around the electrode. A spatial step size of $\Delta s=0.5\,\mathsf{mm}$ and a rotational step size of \Delta\alpha=10\,$$\mathrm{\SIUnitSymbolDegree}. The volumes of the neural activation extent and the corresponding extents determined by the threshold values of the electric potential and the electric field norm were computed by using the Qhull library (http://www.qhull.org/) implemented in SciPy (https://www.scipy.org/, version 0.17.0).

A previous computational study investigated the feasibility of approximating the neural activation extent determined by the coupling of the field distribution with axon models by using constant field threshold values, depending on the stimulation protocol and the axon diameter [13]. The used axon models were based on a multi-compartment mathematical model employing the cable equation, while in this study a double cable axon model is applied [20]. The results showed that iso-volumes for threshold values determined using the electric field norm allowed for a close approximation of the neural activation extent determined from the coupled field-axon models. The deviation between these iso-volumes and the neural activation extents further decreased with increasing fiber diameter of the axons. To determine these threshold values, axon models were used to compute the maximum distance along a line radial to the active electrode center, for which an axon model placed at this distance is still activated by the stimulation. The electric potential or the electric field norm at this distance provided the field threshold value for computing the neural activation extent. This procedure was repeated for a set of stimulation amplitudes in a given range, resulting in averaged field threshold values over the given stimulation amplitude range. The neural activation volumes resulting from this approach and using the electric field norm to determine the field threshold values are denoted as VTAE,const in this study. For evaluating the feasibility of this approach in [13], the threshold values for a given stimulation amplitude range were normalized with respect to the threshold value for the minimum stimulation amplitude.

In addition to this methodology, this study incorporates an approach, which does not average the field threshold values in the given stimulation amplitude range, but computes the threshold-distance relationship along the line radial to the active electrode center, from which the maximum distance to activate an axon for a given stimulation amplitude can be computed in post-processing. In addition, the proposed methodology accounts also for rotationally asymmetric field distributions by computing the threshold-distance relationship in each plane around the stimulation electrode lead and averaging the resulting field threshold values. The corresponding iso-volumes are denoted as VTAφ if electric potential threshold values and VTAE if electric field norm threshold values were applied. To investigate the dependence of the resulting approximation on the fiber diameter as reported in [13], the proposed study cases are carried out for axon fibers with diameters of $2.0\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, $3.0\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, $5.7\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, $7.3\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, $8.7\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, an $10.0\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$. Following the approach in [13], the thresholds are determined by the mean value of the electric potential and the electric field norm at the activation distance radial to the center of the active electrode contact for a given stimulation amplitude within the prescribed range in each plane.

II Results

II-A Validation of the field solution

In order to validate the field solution obtained by the implemented pipeline, a simplified multi-compartment volume conductor model for which the analytical solution is known is used. The analytical example model represents the problem of determining the electric potential distribution in a homogeneous electric field, which is locally disturbed by a conducting layered sphere (Fig. 4). Following the modeling and simulation pipeline described in the methods section, a 3D volume conductor model was composed with a radius for the inner sphere of $0.5\text{\,}\mathrm{c}\mathrm{m}$ and $1.5\text{\,}\mathrm{c}\mathrm{m}$ for the outer sphere. The bounding box around the layered sphere is determined by an edge length of $10\text{\,}\mathrm{c}\mathrm{m}$. A homogeneous electric field of $10\text{\,}\mathrm{V}\mathrm{m}^{-1}$ was generated by applying Dirichlet boundary conditions of $0.5\text{\,}\mathrm{V}$ and $-0.5\text{\,}\mathrm{V}$ on the opposing faces of the bounding box in the $yz-$plane. The conductivity was set to $2.0\text{\,}\mathrm{S}\mathrm{m}^{-1}$ for the inner sphere, $0.1\text{\,}\mathrm{S}\mathrm{m}^{-1}$ for the outer sphere, and $1.0\text{\,}\mathrm{S}\mathrm{m}^{-1}$ for the bounding box. For the quasistatic field problem, relative permittivities of $120$ for the inner sphere, $2\cdot 10^{6}$ for the outer sphere, and $80$ for the bounding box, for a frequency of $35\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ were applied. The mesh was manually refined at the spheres and bounding box surfaces, resulting in a total number of approximately $1.1$ million cells. The analytic solution is determined by solving the field equations (1) analytically using a separation approach, which is explained in detail in the appendix V. The relative deviation of the solution, determined as the norm of the respective field quantity between the computational model and the analytical model along a cut line through the layered sphere was below $7.0\cdot 10^{-3}$ and $2.9\cdot 10^{-2}$ for the electric potential and electric field norm of the stationary field equation and below $1.1\cdot 10^{-2}$ and $2.1\cdot 10^{-2}$ for the norm of the complex-valued electric potential and the electric field norm of the quasistatic field equation (Fig. 5).

II-B Validation of the neural activation solution

To assess the validity of the implementation of the simulation pipeline to determine the neural activation with the myelinated axon cable model [20], the volume conductor model for DBS used in [15] is adapted in order to compare the distance-threshold relation for a fiber diameter of $5.7\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$. The volume conductor model used in the mentioned study comprises the same DBS electrode as in this study as well as an encapsulation layer. A voltage-controlled stimulation signal with a frequency of $150\text{\,}\mathrm{H}\mathrm{z}$ and a pulse duration of $100\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$ is applied. To match the dielectric tissue properties, the conductivity of brain tissue and of the STN was set to $0.3\text{\,}\mathrm{S}\mathrm{m}^{-1}$ and the conductivity of the encapsulation layer to $0.15\text{\,}\mathrm{S}\mathrm{m}^{-1}$. The relative deviation between the threshold-distance relationship was determined to a value below $5.5\text{\,}\mathrm{\char 37\relax}$ by computing the norm between the thresholds determined with the implemented simulation pipeline and the data extracted from [15, Fig. 2], which presents a good agreement of the computed threshold-distance relationship with the data from the original model. Varying the fiber diameter in the considered range described in section I-H resulted in the expected characteristic decrease in the threshold-distance relationships with increasing fiber diameter (Fig. 6).

II-C Approximating the neural activation by field thresholds

The adaptive scheme proposed in I-H was used to determine the neural activation extent for four model setups with varying electrical properties: A homogeneous model (Model 1), a high encapsulation conductivity model (Model 2), a low encapsulation conductivity model (Model 3), and a model including the STN as target area (Model 4). The resulting stimulation amplitude ranges, required to activate a region between $2\text{\,}\mathrm{m}\mathrm{m}$ and $4\text{\,}\mathrm{m}\mathrm{m}$ distance to the active electrode contact, and the corresponding neural activation extents showed a variation for the different models, which results from their varying electrical tissue properties (Fig. 7). This influence is also noticeable in the determined activation thresholds for the electric potential with a threshold value of $-0.25\text{\,}\mathrm{V}$ for the minimum stimulation extent with a normalized increase of $2.64$ for the homogeneous model and a threshold value of $-0.32\text{\,}\mathrm{V}$ with a normalized increase of $1.97$ for the model including the STN. The electric field norm threshold value for the minimum stimulation extent varied between $-129\text{\,}\mathrm{V}\mathrm{m}^{-1}$ and $-141\text{\,}\mathrm{V}\mathrm{m}^{-1}$, but showed a similar normalized increase between $1.29$ and $1.35$ for the models. The varying electrical tissue properties of the models and the resulting varying neural activation extents correspond to changes in the electrode impedance with the high encapsulation conductivity model (Model 2) showing a substantially smaller impedance and also a substantially larger neural activation extent than the other models (Fig. 8).

In order to assess the quality of approximating the neural activation extent by threshold values of the electric potential and the electric field norm, the volumes of the neural activation extent and the extents resulting from the determined threshold values were compared (Fig. 9). While the extents for the neural activation and for the threshold values were in good agreement for the homogeneous model (Model 1), the low encapsulation conductivity model (Model 3), and the model including the STN as target area (Model 4) with average volume deviations of 4.0\text{\,}\mathrm{\char 37\relax}$$\,\pm\,$$1.8\text{\,}\mathrm{\char 37\relax}, the extents for the threshold values understimated the neural activation extent for the high encapsulation conductivity model (Model 2) with average volume deviations of 14.1\text{\,}\mathrm{\char 37\relax}$$\,\pm\,$$2.8\text{\,}\mathrm{\char 37\relax}, as also noticeable in Figure 7B for the minimum and maximum stimulation amplitude. Besides the approach using the threshold-distance relationship to determine the field threshold values, the approach using the mean value of the determined normalized field threshold values as carried out in [13] was used in this study to approximate the neural activation extent by threshold values of the electric field norm. Comparing both approaches for all models by computing the relative deviation between the approximated activation volumes and the neural activation extents, the approach using the threshold-distance relationship ($\mathsf{VTA}_{\mathsf{E}}$) showed smaller average volume deviations of 8.3\pm 3.9$$\text{\,}\mathrm{\char 37\relax} (Min: $3.2\text{\,}\mathrm{\char 37\relax}$, Max: $15.6\text{\,}\mathrm{\char 37\relax}$) compared to 8.8\pm 6.5$$\text{\,}\mathrm{\char 37\relax} (Min: $1.0\text{\,}\mathrm{\char 37\relax}$, Max: $22.2\text{\,}\mathrm{\char 37\relax}$) for the approach using the average value of the normalized field thresholds ($\mathsf{VTA}_{\mathsf{E,const}}$).

Changing the fiber diameter from $5.7\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ to smaller and larger diameters resulted in a variation of the determined stimulation amplitude ranges. The required amplitudes to activate the prescribed distance of $2\text{\,}\mathrm{m}\mathrm{m}$ to $4\text{\,}\mathrm{m}\mathrm{m}$ decreased for increasing fiber diameter, which is also noticeable in the threshold-distance relationship in Figure 6B. The largest increases in the normalized activation thresholds were found for small fiber diameters with values up to $3.37$ for the electric potential and $1.74$ for the electric field norm. With increasing fiber diameter a smaller increase of the normalized activation thresholds for the electric potential and the electric field norm is noticeable, which corresponds to the results reported in [13]. The slope for the normalized activation thresholds for the electric field norm tends to a value of $1$, which allows for approximation of the neural activation extent by using only one threshold value of the electric field norm at an arbitrary stimulation amplitude.

The proposed scheme for the adaptive estimation of the neural activation extent during DBS required approximately $31\text{\,}\mathrm{\char 37\relax}$ to $43\text{\,}\mathrm{\char 37\relax}$ less axons than a non-adaptive approach, where a prescribed number of axons is positioned in a rectangular grid in each rotational plane (Table III). While the estimation whether an axon is activated or not activated by the field for a given stimulation amplitude was carried out by one run of the neuron model, the estimation of the minimum stimulation amplitude for one axon required approximately $23.2$ runs for a tolerance of $1\cdot 10^{-6}$. Considering the number of required optimization runs, a speed up between $38\text{\,}\mathrm{\char 37\relax}$ and $66\text{\,}\mathrm{\char 37\relax}$ is obtained by the adaptive scheme, which reduced the computation time for the rotationally symmetric neural activation computations to $24$ - $45$ minutes saving between $10$ - $27$ minutes on a 12\times$$2.4\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}, $48\text{\,}\mathrm{G}\mathrm{B}$ workstation. The longest computation time was required for the high encapsulation conductivity model (Model 2), since it employed a wide spread of the neural activation extent along the electrode lead due to the highly conductive encapsulation layer. For the heterogeneous case, the adaptive algorithm reduced the computation time from $19$ hours $41$ minutes to $11$ hours $36$ minutes. When the neural activation extent is estimated by the threshold values of the electric potential and the electric field norm, the number of required axons is substantially reduced from determining the activation thresholds at several locations in each plane around the electrode to determining the threshold-distance relationship along a line radial to the active electrode contact. This threshold-distance relationship computation required for a maximum distance of $4\text{\,}\mathrm{m}\mathrm{m}$ radial to the active electrode contact only the computation of the activation threshold of eight axons, requiring approximately 186 runs of the neuron model and a computation time of below four minutes, which is less than the computation time for the field solution, which required approximately five minutes for all models. In case of the heterogeneous model (Model 4), where no rotational symmetry could be exploited, $288$ axons and a computation time of approximately one hour seven minutes was required.

III Discussion

The present paper proposes an adaptive scheme to estimate the neural activation extent during DBS, which is embedded into a Python package using the open-source solutions FEniCS and NEURON. The field solution of the volume conductor model computed with FEniCS as well as the threshold-distance relationship of the axon model computed with NEURON were compared with analytical solutions as well as reference data from literature and show a good agreement with deviations below $2.9\text{\,}\mathrm{\char 37\relax}$ for the field solution and $5.5\text{\,}\mathrm{\char 37\relax}$ for the threshold-distance relationship, respectively (Fig. 5 and Fig. 6). The field model is able to compute stationary current fields (purely resistive material properties) as well as electro-quasistatic fields (complex material properties including conductivity and relative permittivity) for heterogeneous and rotationally asymmetric tissue distributions. The support for incorporating complex material properties further allows for computing the time-dependent field solution dependent on the dispersive electrical properties of biological tissue for any applied voltage- or current-controlled DBS signal using the Fourier Finite element method [8]. The implementation of the field and neuron parts in one Python package made it possible to adaptively estimate the neural activation extent based on the computed field solution and stimulation signal. Instead of solving the field problem and exporting the time-dependent electric potential at the nodes of several axons located at prescribed positions around the electrode lead to determine the minimum stimulation amplitude to elicit an action potential for each axon [7, 11], the adaptive scheme positioned axons only in those regions, where a neural stimulation by the given field solution and stimulation amplitude range would occur. With that, the adaptive scheme requires no pre-knowledge on the neural activation extent for the given volume conductor model and model parameters and, in addition, requires substantially less computational resources and time.

The adaptive scheme was applied to estimate the neural activation extent for model cases with varying tissue properties and axon diameters. The tissue properties were chosen to represent different post-operative stages during DBS as well as a homogeneous and rotationally asymmetric case, where the STN was explicitly modeled as target area. Compared to a non-adaptive approach under the assumption that the minimally required grid size for the axon positions is already known, a speed up of up to $66\text{\,}\mathrm{\char 37\relax}$ was achieved by applying the adaptive scheme.

Nevertheless, even with the adaptive scheme the total computation time for determining the neural activation extent can still be substantially larger than for determining the field solution. To investigate possibilities to further reduce the computational expense, a field threshold approach using the relationship between the field solution and the neural activation suggested in [13] was applied by computing the neural activation along a line radial to the center of the active electrode contact and using the resulting threshold-distance relationship to determine threshold values of the electric potential and electric field norm. In [13], the iso-volumes for these threshold values constituted a good estimate of the neural activation in homogeneous and rotationally symmetric volume conductor models for DBS. In [13], a single-cable axon model was used compared to a double-cable axon model from [20] used in this study. They show a different threshold-distance relationship for the same fiber diameters, with a $3.0\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ single-cable axon model correlating with a $5.7\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ double-cable axon model [13, Fig. 6]. The different axon models and different applied stimulation amplitude ranges impede a direct comparison of the normalized field threshold values in this study with the data from [13]. Nevertheless, larger normalized field thresholds were observered in both studies, when using the electric potential compared to the electric field norm for the field threshold computation (Tab. II and [13, Table I]). In addition, in both studies, the maximum normalized field thresholds decreased for increasing fiber diameters tending to a value close to $1$ when using the electric field norm ($1.06$ for a $7.5\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ axon in [13] and $0.96$ for a $10.0\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ for Model 1, see Table II). Furthermore, the results suggest that also for heterogeneous and rotationally asymmetric field distributions with substantially varying electrode impedances (Fig. 8), threshold values and corresponding iso-volumes of the electric potential and the electric field norm generally constitute a good approximation of the neural activation extent (Fig. 9). However, the results revealed as well that in case of a highly conductive encapsulation layer, as in the phase directly after the surgery, the deviation between the extent approximated by field thresholds and the neural activation extent becomes larger (Fig. 7 and Fig. 9). This increased extent is a result of the increased conductivity of the encapsulation layer, which is spatially connected to the active electrode contact, leading to a higher electric field strength and, with that, an increased activity along the electrode. For this case, the neural activation outside the target area is underestimated, which could have possible implications for determining the impact of unwanted side effects. Nevertheless, the activation in the target area is still approximated with a good quality by the determined field thresholds (Fig. 9). The underestimation of the neural activation for larger distances away from the active electrode by the field threshold approach may be due in part to an approximation artifact. While the response of the axon to a stimulus depends on the electric potential distribution along the whole axon, especially on the second derivative of it [24], the field threshold approximation approach determines the field value for a given stimulation amplitude only at the center node of the axon. For a cathodic stimulation pulse and an axon positioned centered to the active stimulation electrode in an isotropic homogeneous medium, the maximum depolarization occurs at its center node, which means that the determined field threshold is connected to the value of the second derivative of the electric potential along the axon. For larger distances away from the electrode, the electric potential and its second derivative along the axon attenuates. Since the field thresholds were determined for distances between $2\text{\,}\mathrm{m}\mathrm{m}$ and $4\text{\,}\mathrm{m}\mathrm{m}$ from the electrode center, this might present an explanation for the strong correlation between iso-volumes of field threshold values and neural activation extents in the target area, but not for larger distances away from the electrode.

Regarding the given model parameters, the required stimulation distance or the prescribed stimulation amplitude range and whether rotational symmetry can be exploited, the field threshold can be determined by computing the solution of a few hundred to thousand runs of the axon model, which is substantially less than using only spatially distributed axon models in the non-adaptive and adaptive-approach, which requires $10^{4}$ to $10^{5}$ runs of the axon model. The number of required axon model runs directly depends on the neural activation extent. For instance, a larger stimulation amplitude results in a larger neural activation extent, requiring more axon models run to determine its outer shape. Besides the reduction of the computational time achieved by the proposed adaptive and field threshold approximation approach, the computation time could also be decreased by employing a larger time step and smaller number of stimulation pulses. In this study, a time-step of $5\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$ based on [20] and a number of $10\text{\,}$ pulses were used in order to achieve the 1-to-1 ratio in the firing of the axon models with the DBS pulse train. In any case, the accuracy and convergence of the results have to be carefully checked, when larger time steps or shorter pulse trains are used.

The proposed approach reduces the computation time for a model with a rotationally asymmetric field distribution from more than 11 hours to about one hour on a common workstation. Therefore, we belief that this approach has the potential to take the computation of the neural activation extent closer to a real-world application in clinical practice, where (computation) time is an important constraint.

In [13], the neural activation extent is approximated by a constant field threshold for a stimulation amplitude, which is computed from the normalized activation thresholds for the corresponding field quantity, such as the electric potential and the electric field norm, by determining a mean value from the linear fit of the field thresholds for increasing stimulation amplitudes. Field threshold values determined with this approach are used in several computational studies to estimate the neural activation extent during DBS [25, 26, 27]. The determined normalized activation thresholds determined in this study with growth factors between $29\text{\,}\mathrm{\char 37\relax}$ to $35\text{\,}\mathrm{\char 37\relax}$ for a $5.7\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ and $47\text{\,}\mathrm{\char 37\relax}$ to $74\text{\,}\mathrm{\char 37\relax}$ for a $2.0\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ axon suggest that this approach leads to an over- and underestimation of the neural activation extent with volume deviations of up to $24\text{\,}\mathrm{\char 37\relax}$ by using a constant electric field norm threshold for all stimulation amplitudes (Table II). Similar to the results in [13], the deviations decreased for increasing axon diameter (Table II). Considering that generally smaller axon diameters, such as $5.7\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ and below, are used to estimate the neural activation extent during DBS [5, 7, 11, 21], this approach might lead to substantial deviations in the estimated extents when a constant field threshold is used for its approximation. The proposed approach to determine for each model the threshold-distance relationship along a line radial to the active electrode contact to determine the corresponding field threshold for the given stimulation protocol was able to estimate the neural activation extents with deviations below $7.6\text{\,}\mathrm{\char 37\relax}$ using electric field norm threshold values and below $3.2\text{\,}\mathrm{\char 37\relax}$ using electric potential threshold values for the corresponding stimulation amplitudes. The model case representing the acute post-operative phase by a highly conductive encapsulation layer constituted the exception to the approximation quality with deviations of $11.9\text{\,}\mathrm{\char 37\relax}$ and $17.5\text{\,}\mathrm{\char 37\relax}$, respectively, which can be accounted to the spread of the neural activation extent along the electrode lead (Fig. 7). Nevertheless, in contrast to the constant field threshold approach, the suggested threshold-distance field threshold approach ensures that the activation distance in the target region and with that the neural activation extent in the target region equals the extent computed with solely using the axon models. Therefore, the results suggest that this approach is feasible to estimate the neural activation extent dependent on varying dielectric tissue properties, especially for chronic post-operative phases with a low conductive encapsulation layer, and of varying axon diameters while reducing substantially the computational expense.

The axon models to determine the neural activation extent are distributed along normal trajectories in planes around the stimulation electrode, which is a common positioning used for the estimation of the neural activation [5, 7, 11]. This positioning scheme was used in this study in order to compare the simulation results with data from [13]. Considering the anatomy of the target nuclei for DBS, additional knowledge on the orientation and topology of the axons in the target area could provide a more target-specific and realisitic estimation of the neural activation. For instance, the used positioning of the axon models around the electrode is a major simplification compared to the more non-uniform orientation of axonal fibers in and around the STN [28]. The consideration of mutually varying axon fiber orientations and geometries would require a modification of the scheme to estimate the neural activation extent to determine an activation ratio or percentage in the target area rather than a closed volume.

The described process of finding the inactivated hull in each plane around the electrode lead accounts for a monopolar electrode configuration, where one axon seed point is placed in front of the active stimulation electrode in each plane. In case of a multipolar electrode setup, axon seed points have to be placed in front of each active stimulation electrode in order to ensure that a closed inactivated hull is determined for the given electrode configuration. The advantage of a multipolar electrode setup is that the stimulation amplitude for each active stimulation electrode can be adjusted to achieve the optimal neural activation extent. The adaptive algorithm to determine the inactivated hull for all possible combinations of stimulation amplitudes of a multipolar electrode setup would have to be adjusted to include the minimum and maximum stimulation amplitude scenario. Even for multipolar electrode configurations, this would allow to profit from the speed up of the proposed adaptive algorithm in comparison to the non-adaptive algorithm. Regarding the approximation of the neural activation extent by field threshold values, this is not generally given: For a monopolar electrode configuration with one active stimulation electrode, the field threshold approach determines from one threshold-distance relationship one distance value and one corresponding field threshold value for a given stimulation amplitude. This is carried out in each plane around the stimulation electrode lead. In case of a multiplolar electrode configuration, the field distribution varies for the given stimulation amplitude at each active stimulation electrode, which results in interaction effects modifying the threshold-distance relationship for each active stimulation electrode. Therefore, for each given stimulation amplitude configuration of a multipolar electrode configuration, the threshold-distance relationships have to be recomputed for each active stimulation electrode, which reduces the computational speed up. Nevertheless, such an approach would be still substantially faster than recomputing the whole neural activation extent using axon models by the non-adaptive approach. Besides these required modifications on the adaptive neural activation algorithm, the implementation of multipolar electrode configuration support requires also the consideration of varying stimulation amplitudes for the active electrode contacts, which will be the focus of future releases of the FanPy Python package.

The volume conductor model used in this study is embedded into a Python package, which accounts for stationary as well as electro-quasistatic field problems and, therefore, is able to compute the time-dependent field solution employing for resistive as well as dispersive dielectric properties of biological tissue for voltage-controlled and current-controlled stimulation signals. To date, current-controlled stimulation is the favoured stimulation protocol due to the reduced sensitivity of the stimulation impact regarding the dielectric tissue properties and effects of the electrode-tissue interface compared to voltage-controlled sitmulation [22]. For the estimation of the time-dependent field solution during voltage-controlled stimulation, the results of previous studies point out the necessity to incorporate the dielectric effects at the electrode-tissue interface, which show a dispersive as well as non-linear behaviour with respect to the intensity of the stimulation signal [29, 30, 31]. Voltage-controlled stimulation was used in this study to validate the threshold-distance relationship of the implemented axon model [20] neglecting the dielectric effects of the electrode-tissue-interface in order to adapt the model used to generate the reference data [15]. The currently embedded volume conductor model is already able to account for isotropic heterogeneous and dispersive dielectric tissue properties. Besides tissue heterogeneity, the anisotropic dielectric properties of brain tissue can have a substantial influence on the estimation of the neural activation extent and the prediction of side effects during DBS [32]. Therefore, it is planned to incorporate the support for anisotropic conductivity tensors into the field model for future studies.

IV Conclusion

In this study, an adaptive scheme to estimate the neural activation extent during DBS is presented. The computation of the field solution as well as the coupling to axon models and the adaptive computation of their response to the stimulation signal is embedded into an open-source Python package, which was used to estimate the neural activation for varying axon diameters and electrical tissue properties rendering different post-operative stages and target area properties. The determined neural activation extents were used to assess the feasibility of their approximation by field threshold values. By using the threshold-distance relationship for determining the field thresholds and corresponding iso-volumes dependent on the stimulation amplitude, a close approximation of the determined neural activation extents could be achieved, while substantially reducing the computational expense.

V APPENDIX

The field equation (1) for a conducting layered sphere with radius $R_{i}$ of the inner sphere and $R_{e}$ of the outer sphere in an external homogeneous electric field as illustrated in Figure 4 can be formulated using spherical coordinates and a rotational symmetry with respect to the azimuthal angle $\phi$. Within an homogeneous isotropic medium, the field equation for the given problem has the form

[TABLE]

with the radial distance $r$ and the polar angle $\theta$. Since $r$ and $\theta$ can be varied independently, a separation of the potential $\varphi=R(r)\Theta(\theta)$ results in the angular and radial equations

[TABLE]

with the separation constant $\lambda$, where the solution to the angular equation is given by

[TABLE]

using the substitution $x=\cos(\theta)$ and $\lambda=n(n+1)$ with a constant $D_{n}$ and the Legendre polynomial $P_{n}(x)$ for $n=0,1,2,\ldots$ [33]. The radial equation represents an Euler-Cauchy differential equation, which solution is given by

[TABLE]

with the constants $C_{1,n}$ and $C_{2,n}$. Applying the separation equation, the general solution of the potential $\varphi:=\varphi(r,\theta)$ is then given by

[TABLE]

using the substitution $x=\cos(\theta)$ with constants $A_{n},B_{n},C_{n},D_{n}$ for an external homogeneous electric field in $z$-direction $-E_{0}z=-E_{0}r\cos(\theta)$ and a vanishing influence of the conducting layered sphere for $r\rightarrow\infty$. Applying the continuity conditions

[TABLE]

for the electrical potential $\varphi$ and the current density $\boldsymbol{J}$ at the boundaries of the layered sphere and exploting the orthnormality of the Legendre polynomials results in a system of linear equations which has only for $n=1$ a non-trivial solution and non-zero right hand side

[TABLE]

with $\sigma_{i,i}=\sigma_{i}R_{i}^{3}$, $\sigma_{s,i}=\sigma_{s}R_{i}^{3}$, $\sigma_{s,e}=\sigma_{s}R_{e}^{3}$, $E_{0,e}=E_{0}R_{e}^{3}$. The solution of (13) is then given by

[TABLE]