Uncertainty Quantification of Oscillation Suppression during DBS in a Coupled Finite Element and Network Model

TL;DR

This study combines a network model and finite element brain model to quantify how uncertainties in tissue electrical properties affect oscillation suppression during deep brain stimulation, revealing tissue variability impacts DBS efficacy.

Contribution

It introduces a coupled finite element and network model with a stochastic approach to quantify uncertainty in DBS effects on neural oscillations.

Findings

Electrical tissue property variations significantly affect oscillation suppression.

Uncertainty in encapsulation tissue properties has minor influence on DBS efficacy.

The model provides a computationally efficient way to assess DBS variability.

Abstract

Models of the cortico-basal ganglia network and volume conductor models of the brain can provide insight into the mechanisms of action of deep brain stimulation (DBS). In this study, the coupling of a network model, under parkinsonian conditions, to the extracellular field distribution obtained from a three dimensional finite element model of a rodent's brain during DBS is presented. This coupled model is used to investigate the influence of uncertainty in the electrical properties of brain tissue and encapsulation tissue, formed around the electrode after implantation, on the suppression of oscillatory neural activity during DBS. The resulting uncertainty in this effect of DBS on the network activity is quantified using a computationally efficient and non-intrusive stochastic approach based on the generalized Polynomial Chaos. The results suggest that variations in the electrical properties of brain tissue may have a substantial influence on the level of suppression of oscillatory activity during DBS. Applying a global sensitivity analysis on the suppression of the simulated oscillatory activity showed that the influence of uncertainty in the electrical properties of the encapsulation tissue had only a minor influence, in agreement with previous experimental and computational studies investigating the mechanisms of current-controlled DBS in the literature.

Figures

Tables

Table 1. TABLE I: Lower and upper boundary, mean value, and relative deviation of the uniformly distributed conductivities of brain tissue and encapsulation tissue.

Tissue type	Conductivity $\sigma$ $\left[\mathrm{mS}/m\right]$
	Min	Mean	Max	Rel. Deviation
Brain tissue	$64.0$	$77.4$	$90.8$	$16.5\%$
Encapsulation tissue	$32.0$	$38.7$	$45.4$	$16.5\%$

Table 2. TABLE II: Mean value, lower quantile (2.5 percentile), and upper quantile (97.5 percentile) of the uncertain bulk tissue impedance determined from the volume conductor model with uncertain tissue properties.

Stochastic Measure	Impedance [$𝗄\Omega$]	Rel. Deviation
Lower quantile	$15.9$	$15.4\%$
Mean value	$13.8$	$0.0$
Upper quantile	$12.0$	$13.1\%$

Table 3. TABLE III: Sensitivity of the bulk tissue impedance and extracellular potential with respect to the uncertainty in brain σ 𝖻𝗋𝖺𝗂𝗇 subscript 𝜎 𝖻𝗋𝖺𝗂𝗇 \sigma_{\mathsf{brain}} and encapsulation tissue σ 𝖾𝗇𝖼𝖺𝗉 subscript 𝜎 𝖾𝗇𝖼𝖺𝗉 \sigma_{\mathsf{encap}} as well as their interaction.

Quantities of Interest	Sensitivity [$\%$]
	${\sigma }_{\mathrm{𝖻𝗋𝖺𝗂𝗇}}$	${\sigma }_{\mathrm{𝖾𝗇𝖼𝖺𝗉}}$	$({\sigma }_{\mathrm{𝖻𝗋𝖺𝗂𝗇}},{\sigma }_{\mathrm{𝖾𝗇𝖼𝖺𝗉}})$
Bulk Tissue Impedance	$91.1$	$8.9$	$0.0$
Extracellular Potential	$99.9$	$0.1$	$0.0$

Equations

$$\nabla \kappa (\mathbf{r})\nabla \phi (\mathbf{r})=0$$

Full text

I Introduction

Deep brain stimulation (DBS) is a neurosurgical therapy used to treat symptoms of neurodegenerative disorders including Parkinson’s disease (PD). Despite its approval for the treatment of PD and its widely acknowledged success, the mechanisms of action of DBS remain uncertain. Electrophysiological studies in patients and in animal models have shown that PD is associated with excessively synchronized oscillatory neural activity in the cortico-basal ganglia network accompanied by an increased power in the beta frequency range ($15\text{\,}\mathrm{H}\mathrm{z}$-$35\text{\,}\mathrm{H}\mathrm{z}$) [1, 2, 3, 4].

For almost two decades, computational models of DBS have been developed to investigate its mechanisms of action. Among these, biological models of the cortico-basal ganglia network, based on the excitatory and inhibitory connections between the different deep brain nuclei, have been used to investigate the effect of square-wave DBS signals on the network dynamics [5, 6, 7]. In these network models, the DBS signal was applied as an intracellular current to the neurons of the target nuclei, typically the subthalamic nucleus (STN) or globus pallidus interna (GPi) for PD DBS. This modeling approach does not capture the reduction in amplitude of the extracellular potential with increasing distance from the electrode or the effect of the extracellular field distribution along the neuron axons. The distribution of the electric field is determined by the electrode geometry, the electrode-tissue-interface, and the heterogeneity of the brain tissue as well as its electrical properties.

Computational studies based on volume conductor models of the brain have shown that neural activation in the proximity of the stimulation electrode depends strongly on the field distribution caused by the DBS signal and, therefore, on the aforementioned factors [8, 9, 10]. The electrical properties of biological tissue such as brain tissue and encapsulation tissue are subject to uncertainty due to wide variations in the values reported in the literature, most of which is based upon measurement ex vivo [11]. There is considerable variation among reported conductivity and relative permittivity values at lower frequencies in particular. This is due largely to the limitations of available techniques for measuring dielectric properties of biological tissues within the low frequency range, including electrode polarization [11].

Simulation results suggest that the uncertainty in the electrical properties of brain tissue can have a substantial influence on the field distribution and stimulation amplitude required to activate the target nuclei [12]. The impact that the volume conductor effects and material property uncertainty consequently have on the network dynamics and the level of suppression of oscillatory neural activity for a given stimulus, however, is not clear.

To investigate the influence of parametric uncertainties in the electrical properties of brain tissue on the suppression of beta activity in the cortico-basal ganglia network, it is necessary to quantify the uncertainty in the time-dependent extracellular potential distribution in the brain as well as in the firing patterns of the neuronal network. These kind of stochastic problems can be solved with standard brute-force approaches such as Monte-Carlo simulation This approach, however, requires a large number of subsequent evaluations of the computational model. Alternatively, computationally efficient approaches such as that based on the generalized Polynomial Chaos (gPC) technique may be used. This approach substantially reduces the number of model evaluations required by generating a surrogate model of the quantities of interest that may be more easily evaluated. The reduction of the number of required model evaluations for generating this surrogate model is carried out by a multi-variate interpolation, which requires the model to be evaluated only for specific samples of the brain tissue electrical properties to compute the multi-variate interpolation polynomials.

The aim of this study was to explore how uncertainty in the electrical properties of brain tissue can effect predicted levels of suppression of pathological neural oscillations in the cortico-basal ganglia network during DBS. To address this, a coupled model of volume conduction within a rodent’s brain and neural activity within the cortico-basal ganglia network during DBS was developed. The gPC method was then applied to quantify the influence of uncertainty in the electrical properties of brain tissue and encapsulation tissue. The network model is based on a model of the cortico-basal ganglia, which was previously used to investigate the relative effects of antidromic and orthodromic activation of cortico-STN afferents [13]. To enable coupling between the field distribution in the rodent’s brain and the network model, the nodes of the cortico-STN afferent axons were first identified to allow application of the corresponding extracellular potential. The coupled model was then expanded to a stochastic model with the electrical properties of the brain and encapsulation tissue modeled as random variables. This expanded model was used to investigate the influence of varying tissue electrical properties on the suppression of beta-band oscillatory activity in the network. A preliminary version of this study has been reported in [14].

II Methods

II-A Cortico-Basal Ganglia Network Model

The network model is based on a previous cortico-basal ganglia model in which the stimulus was applied to cortico-STN afferent axons using an analytical, purely resistive point source to model the field distribution during DBS [13]. In this network model, the neurons in the STN, and globus pallidus external (GPe) are modeled as single compartment, conductance based Hodgkin-Huxley type models, with threshold crossing models representing the interneurons within the cortex. The STN was based on a parkinsonian rat model which captures the generation of plateau potentials [15]. The cortex was composed of a threshold crossing model representing cortical interneurons [16] and a multi-compartment model of cortical neurons consisting of a soma, an axon initial segment, a myelinated main axon, comprising five nodes of Ranvier and four internodes, and an axon collateral which connects to the STN. The soma was based on the regular spiking neuron model developed by Pospischill et al [17], while the axon initial segment, axon, and collateral were based on a mouse model developed by Foust et al [18]. Each nucleus and neuron type consisted of a population of $30$ neurons. The nuclei were connected by excitatory (glutamergic) and inhibitory (GABA) synapses (Fig. 1). The type and direction of connections between nuclei of the cortico-basal ganglia were based on data available for the network in the diseased state [19], however, it is more difficult to ascertain the exact number of connections between individual neurons in different nuclei. The STN neurons received inhibitory input from the GPe and excitatory input from three randomly chosen cortical neurons. The PD state of the network was simulated by applying a DC current to the cortical neurons, as the STN receives substantial direct cortical input [20]. The model connections and parameters were further designed to match antidromic firing rate results from a parkinsonian rat model [21], which showed an increasing frequency of antidromic spikes as the frequency of stimulation was increased, until reaching a peak at $125\text{\,}\mathrm{H}\mathrm{z}$ stimulation frequency. and increasing the synaptic gains within the network as described in [22], leading to oscillatory activity in the STN neurons within the beta frequency range. The network model was implemented in NEURON 7.1 [23]. A time step of $1\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$ was used and the random seed was reused to ensure same connections between the nuclei and synaptic noises for each computation. The model was initialized for 5 seconds to reach the PD steady-state before the application of DBS. DBS was then modeled by applying the computed extracellular potential for 2 seconds to both the axon collateral and the main axon of the cortical neurons. Pulse trains representing the firing times of all STN neurons were summed and the power spectrum of the resulting composite pulse train was computed. The level of beta activity within the STN was then estimated as the integral of the power spectrum of the composite spike train in the frequency range between $15\text{\,}\mathrm{H}\mathrm{z}$ and $30\text{\,}\mathrm{H}\mathrm{z}$ [24]. The power spectra and level of beta activity in the GPe and cortical neurons were similarly estimated.

II-B Extracellular Field Distribution in the Rodent’s Brain

The volume conductor model of the rodent’s brain is based on a three-dimensional digital atlas database of the adult C57BL/6J mouse brain. The segmented brain atlas data comprises an average from magnetic resonance microimages of 10 adult males and has an isotropic resolution of $47\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ [25]. The brain surface was parametrized using non-uniform rational B-splines with the software Geomagic Studio®  (http://www.geomagic.com). The finite element model of the rodent’s brain was generated using Comsol Multiphysics™  (http://www.comsol.com) by importing the parametrized brain surface model and a model of a DBS microelectrode (diameter: $125\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, electrode contact height: $320\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ [26]), as illustrated in Fig. 2. The microelectrode was positioned in the STN area ventral and lateral to the thalamus by comparing the image data with the Allen Brain Atlas of the mouse brain (http://atlas.brain-map.org). The electrical ground of the model was set to the area of the cerebellum to model the placement of the ground electrode in the neck of the mouse. A square-wave current-controlled DBS signal with a frequency of $130\text{\,}\mathrm{H}\mathrm{z}$ and a pulse duration of $60\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$ was applied to the electrode contact via a Dirichlet boundary condition [27]. The time-dependent field distribution was computed using the Fourier Finite Element Method [8] by solving the Laplace equation

[TABLE]

where $\kappa(\boldsymbol{r})$ is the electric conductivity and $\varphi(\boldsymbol{r})$ the electric potential at location $\boldsymbol{r}$ in the computational domain $\Omega$. Due to reactions of the body to the implant, scar tissue forms an encapsulation layer around the electrode. The electrical properties of this highly resistive layer as well as brain tissue vary considerably within the literature [9, 28]. To capture these variations, their electrical properties were modeled as uniformly distributed random variables with a standard deviation of $10\text{\,}\mathrm{\char 37\relax}$, which corresponds to a relative deviation with respect of the 97.5 percentile and 2.5 percentile of $16.5\text{\,}\mathrm{\char 37\relax}$ (Table I). The mean value of brain tissue conductivity was chosen according to the value of white matter at a frequency of $2\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ [29], which constitutes a good approximation of the dispersive nature of the electrical properties in volume conductor models of the brain [10]. The mean value of the encapsulation tissue was set to half the mean value of the brain tissue [12].

DBS was applied to the cortical collaterals and axons, following the results of an optogenetic study in 6-OHDA rats that found that theraputic effects of STN DBS could be accounted for by direct selective stimlation of afferent fibers projecting from the cortex to the STN [30]. The collaterals of the cortical neurons branching to the STN were randomly distributed within a distance of $250\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ from the center of the electrode contact using uniformly distributed random variables for their Cartesian coordinates (Fig. 2). The collaterals were orientated perpendicular to the main axis of the electrode, with the longitudinal axes of each collateral centered at the electrode. The axons of the cortical neurons, from which the collaterals branched, were orientated perpendicular to the collaterals and parallel to the electrode. To ensure a sufficient resolution of the finite element mesh with respect to the size of the collaterals and the axons, manual mesh refinement was applied by setting a maximum mesh size of $25\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ in the area of the collaterals, $50\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ in the area of the axons, and $10\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ at the electrode contact surface. The final mesh consisted of approximately $1.2$ million elements resulting in approximately $1.7$ million degrees of freedom using quadratic basis functions. The linear system was solved using generalized minimal residual method with a relative tolerance of $1\cdot 10^{-6}$.

II-C Coupling of the Extracellular Field Distribution to the Network Model

The spatial distribution of the cortical neurons resulted in differences in the extracellular potential along the collaterals and main axons of each neuron. To incorporate the effect of variations in extracellular field distribution on the dynamics of the neural network, the extracellular field was coupled to the nodes of the target axons in the discretized network model. Following the approach of previous computational studies to couple the field distribution from a finite element model to neuron models [8, 31], the time-dependent potential of the finite element model was computed at the location of each segment of the cortical axons and collaterals in the neural network model and was applied as a discrete time-dependent extracellular potential.

In a previous study, the extracellular potential during DBS was assumed to preferentially activate the branching collaterals and the extracellular potential was, therefore, applied solely at the collateral [13]. To capture potential effects of the extracellular potential on both the cortical axons and collaterals, in the present study the time-dependent extracellular potential computed using the volume conductor model was applied to each segment along both collaterals and axons at a total of $10$ equidistant nodes along each collateral and $55$ nodes along each axon.

II-D Uncertainty Quantification Using Polynomial Chaos

During the application of DBS, the firing times of neurons within the cortex and STN will vary with the time-dependent extracellular potential. The beta-band activity of the neural network model will thus depend on the time-dependent extracellular potential distribution computed using the volume conductor model. The field distribution, in turn, depends on the conductivity of the encapsulation tissue and brain tissue. In volume conductor models, the isotropic conductivity, and permittivity if considered, of the biological tissues are typically modeled using a single value. There is, however, considerable uncertainty in the true values of the dielectric properties, due to challenges in the measurement of the electrical properties of biological tissue, especially in the low $\text{\,}\mathrm{H}\mathrm{z}$ and $\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ frequency regime [11]. To capture the uncertainty arising from this variation, in this study, the conductivity of each tissue type was described using a random variable to model its uncertainty. To compute the statistics of the quantities of interest, standard stochastic methods such as Monte-Carlo simulation could be used. To provide a sufficient accuracy in the stochastic measures, Monte-Carlo simulation requires a large number of random samples of the model parameters [32]. For each sample, the whole computational model must be solved from which a corresponding sample of each quantity of interest can be computed. If the deterministic model is computationally expensive, as it is the case in the present study, this approach becomes computationally inefficient.

To reduce the computational effort for determining the stochastic measures of the quantities of interest, a non-intrusive variant of the gPC technqiue is applied. This technique determines a surrogate model of each quantity of interest by expanding them in a series of multi-variate orthogonal polynomials and computing the model solution at predefined nodes in the parameter space, which is spanned in this study by the uncertain electrical properties of brain tissue [33]. These nodes are determined by combining the one-dimensional Clenshaw-Curtis quadrature rule with a combinatorical algorithm to generate a multi-dimensional sparse grid representation of the parameter space (sparse-pseudospectral approach). These model solutions at the grid nodes are used to compute the coefficients of the multi-variate polynomial series expansion with multi-dimensional numerical integration. The number of grid nodes depend on the polynomial degree $p$ of the gPC expansion as well as on the number of model parameters $M$ [34]. In general, gPC is computationally advantageous over Monte-Carlo simulation for a small number of model parameters (e.g. $M=6$ and $p=3$ requires $389$ evaluations of the model). In addition, the series expansions of the quantities of interest provided by the gPC approach allow for a simple global sensitivity analysis by Sobol’ indices. These indices determine the sensitivity of these quantities with respect to each uncertain parameter as well as their stochastic interactions by evaluating their conditional variances [35]. In this relation, the term conditional variance describes the variance of the quantity of interest with respect to each uncertain parameter, the uncertain electrical properties of brain tissue, or any combination of them.

The gPC and sparse-pseudo-spectral approach has previously been successfully applied to other bio-electrical applications including optimizing the stimulation protocol in a multi-electrode transcranial direct current stimulation setup as well as for an electro-stimulative hip-revision system [36, 37].

The method was implemented in Matlab as a function package, which determines the deterministic parameter samples on the chosen sparse grid based on the predefined properties of the random model parameters and the polynomial degree $p$, as described in [36]. The list of parameter samples was used to subsequently initialize the solution process of the volume conductor model with COMSOL and of the network model with NEURON. From the resulting list of solutions, the corresponding sample list of the desired quantities of interest including the extracellular potential and the level of beta suppression was derived, from which the respective series expansion, as mentioned above, were computed. The evaluation of these series expansions of the quantities of interest for any samples in the parameter space require only the subsequent evaluation of the product of basis polynomials, which provides a suitable model for Monte Carlo simulation [34]. Therefore, the determined surrogate models of the quantities of interest allow for the application of Monte Carlo simulation with a large number of random samples to provide sufficiently accurate stochastic measures.

III Results

Under simulated parkinsonian conditions, STN neurons in the model exhibited synchronized activity within the beta frequency range (Fig. 3). To investigate the response of the network model to DBS of increasing stimulation intensity, the level of beta suppression in the STN neurons within the network model was investigated for stimulation amplitudes between $50\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$ and $250\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$. At each stimulation amplitude, the corresponding extracellular field distribution was computed using the mean conductivity values of encapsulation and brain tissue, as given in Table I. The beta-band power of the point process representing the firing times of the STN neurons progressively decreased as the amplitude of the DBS signal was increased (Fig. 4(d)). As the DBS amplitude was increased, an increase in the number of activated collaterals and axons was observed, leading to an increase in the region of directly stimulated tissue. (Fig. 3). The collateral and axonal activation orthodromically activated the STN leading to a gradual entrainment of the STN neurons to the stimulus pulse and subsequent reduction of the level of beta-band oscillations (Fig. 4(a)), which is also noticeable for the GPe neurons (Fig. 4(b)). The cortical neurons and inhibitory interneurons were simultaneously antidromically activated, reducing beta-band activity from the cortex to the STN (Fig. 4(c)). For stimulation amplitudes of $150\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$ and higher, a slight increase in the amplitude of the beta-band high-frequencies between $30\text{\,}\mathrm{H}\mathrm{z}$ and $40\text{\,}\mathrm{H}\mathrm{z}$ can be seen. The increased activity in this frequency range is accompanied by a more dense spiking pattern at $250\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$ compared to that at $150\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$ (Fig. 3). Nevertheless, a clear suppression of most of the power spectrum due to an almost complete entrainment of all STN, GPe, and cortical neurons, is observed for all stimulation amplitudes between $140\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$ and $250\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$ (Fig. 4(d)).

Following evaluation of the network response with the tissue conductivities in the volume conductor model set to their respective mean values, the gPC method with a polynomial order of $p=4$ was used to investigate the influence of uncertainty in the tissue conductivity on the quantities of interest. Using the sparse-pseudo spectral approach and the multi-dimensional implementation of the numerical integration method based on the Clenshaw-Curtis rule, the computation of the corresponding gPC series expansions required the evaluation of the computational model for $29$ deterministic samples of the tissue conductivity for the computation of the uncertain extracellular potential. The uncertainty of the extracellular potential within the proximity of the stimulation electrode was approximately $15\text{\,}\mathrm{\char 37\relax}$ and $18\text{\,}\mathrm{\char 37\relax}$ with respect to the upper quantile (97.5 percentile) and the lower quantile (2.5 percentile), respectively. Exemplarily results are presented for a reference cortical collateral positioned perpendicular to the electrode at a radial distance of $100\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ and for a stimulation amplitude of $250\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$ in Figure 5(a). The slightly larger deviation with respect to the lower quantile compared to the upper quantile arises from the inversely proportional relationship between conductivity and potential distribution for current-controlled stimulation, which results in an asymmetric probability density function and, therefore, unequal differences in the mean value [12]. The level of uncertainty in the time-dependent extracellular potential at a specific location was similar at the other locations of the collateral and axon segments (Fig. 5(b,c)). Comparing the average relative deviation of approximately $16.3\text{\,}\mathrm{\char 37\relax}$ in the extracellular potential to the uncertainty in the tissue conductivities of approximately $16.5\text{\,}\mathrm{\char 37\relax}$ (Table I), an almost 1:1 relation between the level of their uncertainty can be noted. A similar relationship was observed for the uncertainty in the bulk tissue impedance of the volume conductor model (Table II). An increase of the conductivity of brain tissue as well as encapsulation tissue resulted in a decrease in the bulk tissue impedance as well as the extracellular potential distribution, and vice versa. Regarding the accuracy of the stochastic measures in these quantities of the volume conductor model, a doubling in the polynomial order of $p=8$, which required the computation of additional $36$ model evaluations, resulted in a difference of below $9\cdot 10^{-5}\,$\text{,}\mathrm{\char 37\relax}$$ and $0.03\text{\,}\mathrm{\char 37\relax}$ in their mean values and relative deviations, respectively. Since the extracellular potential distribution within the volume conductor model scales linearly with respect to the applied stimulation amplitude due to the linear material properties and the Laplace equation (1), the solution for different stimulation amplitudes was obtained by simply multiplying the reference solution by the corresponding stimulation amplitude.

While this linear property of the volume conductor model allowed for a reuse of the reference solution, which only had to be computed once, the deterministic solution of the network model had to be reevaluated for each stimulation amplitude due to the nonlinearities in the dynamics of the network model and individual neurons [13]. This nonlinearity in the response of the network model with respect to the stimulation amplitude is apparent in the suppression-amplitude relation shown in Figure 4. Recent studies have demonstrated suppression of beta band activity witihn the basal ganglia local field potential (LFP) during DBS, accompanied by a reduction in the parkinsonian motor symptoms of rigidity and bradykensia [38]. The presence of beta band LFP oscillations are interpreted as being indicative of pathological synchronization between neurons and have attracted substantial interest as a potential biomarker of the parkinsonian state [39]. The influence of uncertainty in the tissue electrical properties on the level of suppression of beta frequency neural oscillations in the STN was thus chosen for investigation. The effect on uncertainty in the material properties of the brain and encapsulation tissue on uncertainty in the level of suppression of synchronous beta band activity among the STN neurons is shown in Figure 6. For lower conductivity values, as in the case of the lower quantile values, a decrease in the STN beta-band activity occurred at lower stimulation amplitudes as the potential within the bulk tissue for a given stimulus amplitude was higher. Higher stimulation amplitudes were correspondingly required to elicit comparable levels of beta suppression when higher conductivity values, as in the case of the upper quantile values, were assumed within the model. Therefore, the uncertainty in the tissue conductivities introduced a large uncertain area in the window of effective stimulation spanned by the lower and upper quantiles (Fig. 6(a)).

Stimulation amplitudes below this window were too low to sufficiently modify the membrane potential along the cortical axons and collaterals to have an effect on the spiking pattern in the STN neurons. Conversely, stimulation amplitudes above this window seem to have no additional effect on the spiking pattern. This effect can be explained by investigating the membrane potential along the cortical neurons. These neurons were already entrained by the DBS stimulus and, therefore, a further increase in the stimulation amplitude did not influence the dynamics of the already stimulated cortical axons (Fig. 3). The resulting level of uncertainty was greatest for a stimulation amplitude of $120\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$ with an average value of approximately $260\text{\,}\mathrm{\char 37\relax}$ (Fig. 6). The asymmetry in the relative deviation with respect to the upper and lower quantile of the extracellular potential and bulk tissue impedance (Fig. 5 and Table II) is also noticeable in the uncertainty of the level of beta suppression. In general, the uncertainty in the level of beta suppression in the network model was substantially larger than the predefined uncertainty in the tissue conductivities, which can be ascribed to the nonlinearity in the mathematical formulation of the network model. This nonlinear behavior affects also the accuracy of the stochastic measures with respect to the polynomial order of the gPC. The uncertain beta supression of beta activity in the network model was computed using gPC with a polynomial order of $p=8$, which required the evaluation of $65$ network models for each stimulation amplitude. A doubling of the polynomial order of $p=16$ resulted in a difference of less than $0.1\text{\,}\mathrm{\char 37\relax}$ and $9.1\text{\,}\mathrm{\char 37\relax}$ in the mean value and relative deviations of the level of beta suppression, respectively, evaluated at a stimulation amplitude of $120\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$, which is substantially larger than for the stochastic measures of the volume conductor model. Nevertheless, the accuracy in the stochastic measures of the network model is sufficiently large for an investigation of the overall trend of its uncertainty. Regarding the required number of random samples to determine the stochastic measures of the surrogate models obtained by the gPC method, a number of $100.000\text{\,}\mathrm{r}$andom samples was sufficient to provide an accuracy of approximately $0.1\text{\,}\mathrm{\char 37\relax}$.

The investigation of the sensitivity of the quantities of interest by evaluating the Sobol’ indices revealed that the major influence on the uncertainty in the extracellular potential and the bulk tissue impedance is due to the uncertainty in brain tissue conductivity (Table III). As expected, due to the linearity of the volume conductor model, no interaction effects of the uncertainty in both tissue conductivities is noticeable. Since the dynamics of the network model depend on the extracellular potential along the axons and collaterals of the cortical neurons distributed in the volume conductor model, a similar behavior would be assumed for the sensitivity in the level of beta suppression. Investigating this sensitivity for stimulation amplitudes between $120\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$ and $170\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$, this is confirmed (Fig. 7). Nevertheless, the influence of uncertainty in encapsulation tissue conductivity as well as its interaction with the uncertainty in brain tissue conductivity seem to increase for stimulation amplitudes below and above that range, which is assumed to be caused by dynamic effects related to the random noise in the network model and high sensitivity of the neurons to the extracellular potential when close to the threshold for excitation.

IV Discussion

The coupling of a cortico-basal ganglia network model with the extracellular field distribution in the target area of DBS allows for an investigation of the mechanisms of action of DBS in a simulation model that captures both network and volume conductor effects. Analytic models of the extracellular field distribution add the spatial location and extracellular stimulation of the neurons to the network model, but are unable to consider effects of the geometry of the electrode, the brain anatomy, and the electrode-tissue interface. Therefore, in this study a 3D finite element model of a rodent brain was coupled to a network model to investigate the influence of uncertainty in the conductivity profiles of biological tissue on beta-band suppression in the network model during DBS. The application of the gPC method allowed for a computationally efficient quantification of the uncertainty in the quantities of interest of the volume conductor model and the network model on a common workstation, which would be infeasible if standard stochastic methods such as Monte Carlo simulation were used. The computation time and memory consumption for one realization was approximately 2 hours and $3.5\text{\,}\mathrm{G}\mathrm{b}$ on a $2.4\,$\text{,}\mathrm{G}\mathrm{h}\mathrm{z} core of a common workstation, with over $95\text{\,}\mathrm{\char 37\relax}$ of the computation time utilized by the network model. The computation of the realizations for all tissue conductivity samples and stimulation amplitudes required less than $10$ days on $12\times 2.4\,$\text{\,}\mathrm{G}\mathrm{h}\mathrm{z} cores of the workstation. A Monte Carlo simulation with comparable accuracy, which would have required a large number of several thousand random samples, would take an estimated computational time of several months and a substantially larger amount of storage. Further, the surrogate models determined by the gPC method of the investigated quantities of interest consist of a series expansion of multi-variate polynomials, which are evaluable for a large number of random samples by Monte Carlo simulation in a comparably small amount of time.

The window of stimulation amplitudes, at which the level of beta suppression is sensitive to uncertainty in the electrical properties of brain and encapsulation tissue, ranges from $110\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$ to $185\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$ (Fig. 6) for the assumed uncertainty in the tissue conductivities of the volume conductor model (Table I). The stimulation amplitude required to effectively modify the spiking patterns in the neurons of the network model, and subsequently the level of beta suppression, depends not only on the tissue conductivity, but also on various factors such as the exact location of the neurons in the proximity of the electrode, the properties of the target neuron axons, the electrode geometry, and the modeling of the encapsulation layer. Despite these factors, similar stimulation amplitudes to suppress PD symptoms in experimental studies in rodents have been reported with values between $30\text{\,}\mathrm{a}$nd $100\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$ for a $130\text{\,}\mathrm{H}\mathrm{z}$ DBS stimulus effective [40, 41]. In [41], an effective stimulation radius of $250\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ for a stimulation amplitude of $100\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{A}$ was reported. The axon collaterals in the present study were positioned within a similar radius around the stimulation electrode. The magnitude of the beta-band oscillatory activity within the network was substantially reduced at stimulation amplitudes for which not all STN neurons were directly entrained by the DBS (Fig. 4). This result suggests that the estimation of the “volume of tissue activated” alone to predict neural activation during DBS [8, 9, 12] may also overestimate the required stimulus amplitude to obtain therapeutic effects of DBS, as the propagation of stimulation effects through the connections within the cortico-basal ganglia network is not considered.

The volume conductor model of a rodent’s brain used in this study consists of a DBS electrode, which is surrounded by an encapsulation tissue layer and is placed within the region of the STN, and the bulk brain tissue medium. The model allows for computing the field distribution around the electrode with respect to the influence of the encapsulation and brain tissue. Given that this medium is homogeneous, the extracellular potential around the cylindrical electrode resembles an almost spherical field distribution with increasing distance to the electrode with the exterior boundary of the rodent brain having most likely only a minor influence. Nevertheless, neural elements can be located very close to the encapsulation layer, as it is the case in this study, which requires a sufficiently detailed model of the electrode and the neighboring areas. In addition, the generated 3D rodent brain model allows arbitrary orientations and locations for the electrode and neuronal elements inside the brain to be chosen and the brain structures to be incorporated. In the future, the model will be extended with models of basal ganglia brain structures, which will result in heterogeneous and anisotropic electrical properties of the corresponding tissue types within the model and, therefore, result in a complex, non-spherical field distribution.

The network model used in this study focuses on the STN-GPe loop and cortical inputs to the STN. The application of DBS to the model is based upon the hypothesis that DBS exerts its effects, at least partially, through the activation of STN afferent fibers projecting from the cortex [42]. In this model, DBS activates corticofugal fibre afferent projecting from the cortex to the STN. The resulting activity antidromically modulates the firing patterns of the cortical neurons while also orthodromically activating the single compartment STN neurons. Previous computational studies have simulated the effect of extracellular potential on isolated neurons [43, 44] or used a simplified point source approximation to apply the potential to the network [13]. Though the present study is the first to combine a detailed volume conductor model with a network model, it is likely that DBS may also effect other neurons and network circuits, not included here, either directly or indirectly [45]. In particular, the question as to whether DBS results in excitation or inhibition of STN neurons remains a controversial topic with both excitatory and inhibitory effects reported [46], alongside the possibility that DBS may inhibit activity of the STN soma while simultaneously activating the afferent output of the STN [44]. There is greater agreement regarding the effects of DBS on the afferent output of the STN, with it generally accepted that STN DBS results in stimulation of the STN output time-locked to the stimulus [47]. Previous models of the corticobasal ganglia network have stimulated STN DBS through the intracellular stimulation of single compartment STN neurons [7, 48, 6, 5]. The single compartment representation of STN neurons in the present study approximates the “input-output” characteristics of the STN neurons and does not rule out the possibility that the cell-body of the STN neurons may be inhibited while the nerve axons are simultaneously stimulated. Multi-compartment models of these nuclei, especially of the STN, would increase the physiological accuracy of the model. However, this would also substantially increase the computational time and the use of computational resources. For example, the model used in this study considered in total $1,950$ neural segments along the cortical neurons for extracellular stimulation and required for one realization approximately $3.5\text{\,}\mathrm{G}\mathrm{b}$ memory. Replacing the STN neurons in this study by multi-compartment models, as used for example in [8], would require extracellular stimulation of additional $6,630$ neural segments, which would result in a substantial increase of the required computational resources and computational time. The number of neurons in each nucleus of the network model, while consistent with what has been used in previous computational studies [5, 49], is lower than it is physiologically realistic. The relatively low numbers of neurons, which are highly interconnected throughout the basal ganglia-network, make the model sensitive to changes in electric potential as even small changes in potential can quickly increase the percentage of neurons activated.

The development of realistic volume conductor models coupled to network models of the basal ganglia move towards the prediction of required stimulation amplitudes to obtain therapeutic effects of DBS, based on realistic brain anatomy and electrical properties as well as the electrode-tissue-interface and neural dynamics. For this problem class and computational models, the results of this study are of crucial importance, since they show the substantial influence of brain tissue uncertainty on pathological neural oscillations which are one possible objective function of the corresponding optimization problem. Therefore, an objective function omitting this uncertainty may not sufficiently render patient-specific stimulation profiles, but only represent one possible solution.

Regarding the minor influence of the uncertainty in encapsulation tissue conductivity on the level of beta suppression and the extracellular potential, the results of this study are in agreement with other computational and experimental studies concerning current-controlled DBS [8, 9, 50]. For voltage-controlled stimulation, a more precise knowledge of the encapsulation layer properties in the post-operative state may be crucial for modeling the mechanisms of action of DBS in the rodent’s brain [12].

V Conclusion

Uncertainty in the electrical properties of biological tissues, including grey and white matter and the encapsulation tissue which forms around the electrode, represents a significant limitation of existing models of the electric field in the vicinity of the DBS electrode. The potential impact of this uncertainty on the dynamics of neural firing and pathological neural oscillations within the cortico-basal ganglia network during DBS is unclear. The methodology presented in this paper offers a computationally efficient method with which to quantify these effects by coupling a finite element model of the electric field with a physiologically-based network model, and applying a stochastic approach based on the generalized Polynomial Chaos. The results of this study indicate that modest variations in the electrical conductivity of brain tissue of $13-15$ $\text{\,}\mathrm{\char 37\relax}$, can result in substantial deviations, of up to $300$ $\text{\,}\mathrm{\char 37\relax}$ about the mean, in the level of beta suppression among STN neurons during DBS. The methods proposed offer an efficient method with which to incorporate known uncertainties in parameters of interest in models of DBS.