Towards an optimised deep brain stimulation using a large-scale computational network and realistic volume conductor model

Objective. Constructing a theoretical framework to improve deep brain stimulation (DBS) based on the neuronal spatiotemporal patterns of the stimulation-affected areas constitutes a primary target. Approach. We develop a large-scale biophysical network, paired with a realistic volume conductor model, to estimate theoretically efficacious stimulation protocols. Based on previously published anatomically defined structural connectivity, a biophysical basal ganglia-thalamo-cortical neuronal network is constructed using Hodgkin–Huxley dynamics. We define a new biomarker describing the thalamic spatiotemporal activity as a ratio of spiking vs. burst firing. The per cent activation of the different pathways is adapted in the simulation to minimise the differences of the biomarker with respect to its value under healthy conditions. Main results. This neuronal network reproduces spatiotemporal patterns that emerge in Parkinson’s disease. Simulations of the fibre per cent activation for the defined biomarker propose desensitisation of pallido-thalamic synaptic efficacy, induced by high-frequency signals, as one possible crucial mechanism for DBS action. Based on this activation, we define both an optimal electrode position and stimulation protocol using pathway activation modelling. Significance. A key advantage of this research is that it combines different approaches, i.e. the spatiotemporal pattern with the electric field and axonal response modelling, to compute the optimal DBS protocol. By correlating the inherent network dynamics with the activation of white matter fibres, we obtain new insights into the DBS therapeutic action.


Introduction
Abnormal functioning of the basal ganglia (BG)thalamo-cortical circuit plays a central role in several movement disorders.One main characteristic of the altered circuitry in Parkinson's disease (PD) and dystonia is a change in oscillatory patterns of the cortico-BG network [20,59,86].In PD, abnormally pronounced β-activity  emerges in the entire BG network.This is thought to underlie the severe motor symptoms of PD, such as hypokinesia [16,43,76].The synchronous elevated β BG rhythm should, in principle, also affect thalamic activity by sending strong inhibitory signals via GPi (globus pallidus interna) and SNr (substantia nigra pars reticulata) efferents to motor nuclei of the thalamus, i.e. ventral anterior and ventrolateral nuclei (VA and VL) [27].Studies in MPTP-treated monkeys and optogenetic mouse models on activation of inhibitory GPi-thalamic projections show that GPi stimulation results in rebound burst discharges in the thalamic nuclei [44,46] after the initial inhibition [49].The enhanced bursting activity due to hyperpolarisation of thalamic neurons comes from an increased inhibitory BG output tone and a Ttype calcium channel-dependent mechanism [1,46].In summary, this would predict a slowing, malfunctioning of thalamic oscillatory activity [27].
Deep brain stimulation (DBS) of the BG was shown to be an efficient treatment for movement disorders [19,80], but its therapeutic mechanism is still not fully understood.The main hypotheses put forward so far are (a) suppression of neuronal activity of the target area, (b) neuroprotective and plasticity effects (including structural plasticity), and (c) enhanced neuronal activity that results in firing pattern alterations, in particular, disruption of hypersynchronised β-band oscillations [15,32,43,46,64].Indeed, recordings in animal models [15,46,85] and observations from computational network models [28,63,64,67] suggest that DBS in the nucleus subthalamicus (STN) results in a more periodic and regular firing at higher frequencies in the BGthalamic network.
Certain progress in predicting therapeutic DBS effects was achieved using the computational concept of the volume of tissue activated [11,12].Although it became a standard approach when estimating DBS outcomes with computational models, it does not consider the anatomy of brain fibre tracts and, importantly, does not take into account network effects.A qualitative advance in stimulation modelling is an estimation of DBS-induced activation of neuronal pathways and propagation of activity in the respective networks.For various neurological diseases, computational and clinical studies show a correlation between axonal firing in specific pathways and the patient condition [10,29,35,42].Nevertheless, it is not well understood how the pathways' activation affects the neural activity of the respective networks.Prediction of DBS outcome based on network models is not novel, but the induced stimulation current is usually simulated as an increased abstract activation [63,64,69].In this paper, we couple pathway activation modelling (PAM) with a computational network of the BGthalamo-cortical circuit in order to optimise DBS treatment for Parkinson's motor symptoms.
Firstly, we use structural connectivity obtained in previous multi-modal imaging data analysis [13,50,62].Secondly, we simulate neuronal populations in different areas (GPe/GPi, STN, etc), placing and connecting neurons using the aforementioned structural connectivity.Each node-neuron is modelled with a variation of Hodgkin-Huxley's current-balance equations [36,64,73].Then, integrating this high dimensional nonlinear system, we simulate spatio-temporal patterns consistent with healthy (normal) and Parkinsonian states.We define a biomarker for thalamic spatio-temporal activity to quantify the pathological thalamic behaviour in the PD state.Using the network, we optimise pathway per cent activation (i.e. per cent of activated fibres) based on the deviation of the biomarker from the healthy state.After the optimal input is determined, it is assigned as the goal function for the electric field optimisation via adjustment of the electrode position and the stimulating current.The ultimate aim of the study is to obtain a DBS protocol that provides the therapeutic pattern encompassing the entire BGthalamo-cortical network.

Data sources
To describe the structural connectivity of the network, data from different studies on human brain anatomy were utilised.The BG nuclei and their substructures were taken from the DISTAL Atlas [13,22], and the Melbourne Subcortex Atlas [74] was used to delineate substructures of the thalamus, while the relevant cortical regions were selected using the Brainnetome Atlas parcellation [23].Fibre tracts classified to pathways in the vicinity of the STN were taken from [62], and projections of the VA nucleus to motor-cortical regions, required to complete the BGthalamo-cortical network, were extracted from the structural group connectome of 90 PPMI PD-patients GQI [50], post-processed in [22].All data were represented in the Montreal Neurological Institute space.

Structural connectivity of the BG-thalamo-cortical network
To investigate DBS network effects in PD, the classic circuit model of the BG-thalamo-cortical network was employed (figure 1(A)) [55].It consists of three inputs from motor-cortical regions: the direct pathway that passes through the striatum and continues as a GABAergic projection to the GPi and SNr; the indirect pathway that also traverses the striatum but has GABAergic projections to the GPe (external segment of the globus pallidus), which in turn inhibits the STN, GPi and SNr; and the hyperdirect pathway through which the STN receives direct excitatory input from the cortical areas.The glutamatergic efferents of the STN innervate the GPe, GPi and SNr.GABAergic projections of the two latter nuclei to the VA and the VL regions of the thalamus represent the output of the BG circuit [7].The circuit loop is closed by excitatory input of the thalamus to the motorcortical regions.
The intrastructural connectivity of the network was obtained from fibre tracts that were either already assigned to pathways [62], classified based on their position relative to the involved structures (thalamo-cortical projections), or substituted by direct current modulations (striatal efferents, see below).The classified fibre tracts are presented in figure 1.Furthermore, the fibres from [62] were grouped depending on the subregions of brain structures they connect.In this study, the simulated network (figure 1(B)) does not include the substantia nigra (SN).The dopaminergic projections of the SNc (substantia nigra pars compacta) have a low level of myelination [33].They are hence less excitable (by approximately two orders of magnitude) by extracellular electric fields [72].The SNr is often paired with the GPi in the context of the direct/indirect pathways, although the fibre trajectories related to STN-DBS electrodes are distinctly different.To ensure homogeneity of the BG pathways, only those pathways present in [62] were employed, which did not include the striatonigral and the nigrothalamic projections.

Modelling structural connectivity using complex network theory and pathway classification
The topological structure of the brain network plays a critical role in the emergent neuronal activity and brain functionality.Structural brain connectivity features were modelled by defining neuronal elements as nodes and the interconnections among them as edges of the network.The underlying connectivity can be described and simplified using complex network theory [9,70].
We assimilate the network structure of section 2 to build a realistic computational model that provides insights into the mechanisms behind the brain functionality or dis-functionality (i.e.movement disorders) and predict the neuronal dynamics on multiple scales [37,38].Specifically, we build a directed network using the previously described structural connectivity.The intrastructural connectivity within the areas (STN, GPe GPi, VA, motor cortex (MC)) was defined assuming an internal complex network structure (see section 2.3.2).

Construction of the inter-connectivity using structural data analysis
We express the structural connectivity of section 2.2 using network theory, i.e. the connectivity is written as a graph G = (V, E), where V is the set of nodes and E represents the set of edges.The nodes of the structural network are defined as points in the threedimensional space and correspond to the starting and ending points of a fibre tract.Also, due to limitations of imaging, the network resolution is set at 1 mm 3 .The structure of the network is described with the adjacency (or connectivity) matrix A: if there is a fibre tract starting at position p i = (x i , y i , z i ) and ending at p j = (x j , y j , z j ) then A(i, j) = 1; otherwise A(i, j) = 0.At the given resolution of 1 mm 3 , the network contains 134 STN nodes, 244 and 246 GPe and GPi nodes, respectively, 833 thalamic nodes and 2070 cortical nodes [68].

Internal connectivity using small world topology
In this section, we model the internal connectivity of each region (i.e.how the nodes-neurons are connected within a region) using small-world structures [4,9,70,83].In such small-world complex networks [60,83], each node interacts with its k nearest neighbours; additionally, a few randomly chosen remote connections (with a small probability p) within the area are also formed [83].Small-world structures are commonly used/emerged in computational and clinical neuroscience [4-6, 17, 24, 58, 66, 68].Note that the model only considers the respective motor parts of the different structures and not any other functionally distinct areas.
The GPe/GPi, VA and MC layers are modelled as separate small-world networks.Each node increases the initial number of connections (or the degree of the node) by k = 20 degrees on average.The new link lies at a distance ⩽5 mm (defining the local neighbourhood).We assume that 95% of the neurons are interconnected locally.Nonetheless, the small-world topology [83] also allows remote connections (at a distance >5 mm) with a small probability of 5%.The phenomenological values of k and p in this model are chosen such that the network resembles both random and lattice structure as it appears in many studies in neuroscience [4-6, 17, 24, 58, 66].With the chosen values, the connectivity structure resembles real neuronal connectivity as it is shown, for example, in the work of [17] and [58].
For the STN, we chose a modified sparse small world approach inspired by experimental findings of [2,30].In our model, 80% of STN does not form local STN connections, and the remaining 20% of STN neurons have an average of 25 connections each, while few of these are randomly chosen remote connections [30,68,69].

Modelling the dynamics of the thalamo-cortical BG circuitry
We construct neuronal populations in the different areas (GPe, GPi, STN, etc), placing and connecting neurons using the aforementioned connectivity of the augmented network.Depending on the region, each node-neuron is modelled with a variation of Hodgkin-Huxley current balance equations [36,64,73], i.e. the membrane potential of the ith node is given by: where C is the membrane capacitance and V i represents the membrane potential.The first term on the right-hand side of equation ( 1) represents the summation of ionic currents, the second, the summation of synaptic currents and the third term represents the DBS stimulus (which is simulated for all BG areas by different levels of direct current injection, although far from the stimulation site, distant nuclei are activated by synaptic connections and probably remain little affected by the stimulation electrode) or sensorimotor inputs in the cases of the thalamus and motor areas.Different areas show different properties for the currents in equation (1) The model contains the following brain regions: the nucleus STN, the globus pallidus externa (GPe) GPi, the thalamus (VA) and the MC region.We model the connection between neurons according to the augmented structural connectivity of section 2.3.2.In the supplementary section, a detailed mathematical description is given.

Macroscopic biomarker of pathological thalamic activity in the Parkinsonian state
The highly synchronised activity (characteristic for PD) coincides with enhanced bursting activity of the BG and thalamus [75,86].Activity in the β-band range is found in both healthy and Parkinsonian states, and it is thought to be associated with the control of ongoing motor activities [59].In the Parkinsonian state, this activity is significantly enhanced [59] [43]; however, the reason for this enhancement is not fully understood [21,59].
The increased BG inhibitory output reduces the excitability of postsynaptic thalamic neurons [44].As a result of the increased inhibitory drive on thalamic motor nuclei, thalamic neurons fire in low-frequency burst discharges [27,46].This pathological activity refers to rebound bursting [27,46,49], as thalamic neurons rhythmically recover from hyperpolarisation, involving T-type calcium channel activation [27,46], and likely also HCN channels [45].This altered behaviour is related to Parkinsonian changes, such as tremors and hypokinesia, but also arousal states in rats, monkeys, and humans [46,49].In order to quantify the thalamic activity, we define a macroscopic variable that takes into account precisely these two main characteristics of the disturbed abnormal thalamic behaviour, i.e. a decreased neuronal activity in conjunction with an increased bursting behaviour: where X represents the stimulated state (i.e.healthy (normal)/Parkinsonian/DBS), meaning that the parameters of the model are tuned in accordance to each simulated state, see also section 3.1, BR is the average number of the bursting blocks (in the thalamic area) and SP represents the average number of spikes (in the thalamus).The parameter a ∈ [0, 1] defines a convex combination and controls the behaviour of the macroscopic function f a , which summarises thalamic activity by counting the relative amount of bursting and spiking activity.
Experimentally and computationally, it is challenging to identify bursts because bursting sets a parameter (mean inter-spike interval) according to intrinsic properties of the detected burst spike trains [14].Here, we follow the simplifying method of [14] to detect bursting blocks.When the mean value of several successive inter-spike intervals is not larger than 15 ms, this behaviour is defined as burst firing.
Another macroscopic variable that is used in this paper is the mean voltage activity V of neurons in an area (e.g.STN); specifically, we define: The mean voltage activity V (related to the local field potential) is utilised for the characterisation of synchronised rhythm (through the Fourier spectrum) under different states (healthy, Parkinsonian or DBS).

Pathway per cent activation of DBS using the network model
The effectiveness of the DBS treatment depends on (a) which specific pathways are activated and (b) on the level of activation [35,40].In the network simulations, we use a one-compartment soma neuron model for all areas.Thus, the activation of an edge (i.e. an axon) due to DBS is modelled as a supra-threshold DBS stimulus on the neuron node, which is the ending point of the axon (i.e. the postsynaptic neuron).
We define the per cent activation A x→y from nucleus x to y as: where N A,x→y is the number of activated edges from x → y and N x→y is the total number of edges (axons) from x to y nucleus.In this study, we analyse the following per cent activations, A 1 = A MC→STN , A 2 = A STN→GPe , A 3 = A STN→GPi , and finally A 4 = A GPi→VA .We developed an algorithm using the Metropolis optimisation [47] for randomly chosen activated nodes (described in the supplementary material) consistent with per cent activation A x→y .
For each set of the parameters A = (A 1 , A 2 , A 3 , A 4 ), the model produces a different neuronal activity whose effectiveness is evaluated from the objective function.The objective function should capture the main impact of DBS, i.e. to increase the thalamic neuronal activity and simultaneously reduce thalamic bursting activity (see also section 2.5) Substituting the form of f, i.e. equation ( 2), we obtain: The values of the model parameters A i , i = 1, . .., 4, were estimated numerically by minimising objective function i.e. equation (6).
The minimisation problem was solved using the MATLAB function fmincont, a gradient-based method for nonlinear problems with constraints for the function class C 1 [52].The step size tolerance was set to tol(X) = 0.001 and the function tolerance was set to tolF = 1.The maximum number of iterations was set to 100.The optimal firing rates now are referred to as A optimised .They will constitute the target vector of equation ( 7) of section 2.8 in order to compute corresponding DBS protocols.

Estimating realistic pathway per cent activation
While the optimised per cent activation, suggested by the network, is defined on abstract neurons, it is important to consider spatial constraints of the activation extent related to the mutual proximity of various neuronal circuits, including regions of avoidance, e.g. the internal capsule implicated in motor sideeffects for STN-DBS [77].Note that the external stimulation alters the transmembrane potential on both axons and somata with dendritic trees, but the former has been shown to be more excitable by external stimuli [61].Therefore, here we quantify activation by Here the electrode is shown at the position that was optimal for the normalisation of the VA-thalamic activity.Note that in order to reduce the computational effort, the axons do not cover the whole length of the pathways but only segments exposed to the high electric field.(B1), (B2) Activation status of the network pathways and other projections (corticofugal and associative pallidosubthalamic), respectively.By solving the cable equation for the computed distribution, responses of passive axons to the stimulation are quantified.If the axon elicits an action potential, it is considered activated.Note that some axons traverse the electrode and/or the encapsulation layer; those are considered 'damaged' and excluded from the simulation.(C) Electrode tip placement is optimized within a 4 mm bounding box centred at the motor aspect of the STN.
the axonal response to DBS.In figure 1, projections that are expected to be reinforced by stimulation are shown in thick lines, and only on these projections were axon models allocated.It is important to note that the DBS-induced extracellular field might affect not only pathways projecting from/to the STN but also passing fibres, such as the pallidothalamic tract.Additionally, we monitored activation in the corticofugal tract to control for current protocols that can cause capsular side effects and the associative portion of the pallidosubthalamic projections, which stimulation might affect cognitive function.
To probe an initiation of action potential in the neuronal tissue, a double cable model of the myelinated axon was employed [54].This widely used model of the mammalian nerve fibre describes a detailed morphology of internodal segments and complex ion channel dynamics.Importantly, no tonic firing was modelled at this stage so that action potentials could be elicited only due to the DBS-induced field.The initial membrane potential was set to −80 mV, and fibre diameters of 5.7 µm were chosen, except for the hyperdirect pathway, where the collateral diameter was set to 3.0 µm.It should be noted that the majority of fibres in the BG are thinner [26,51]; however, larger diameters were also reported [79].An array of previously published works in the computational field uses fibre diameters of, or close to, 5.7 µm [31,35,42,54], and consequently we chose to adapt this value.To reduce computational costs, axon lengths for long projections were limited (e.g.16 mm for the hyperdirect pathway) since the electric field magnitude further from the STN is negligible.Seeding of an axon was always initiated at the fibre tract's segment closest to stimulation contacts, and in total, the algorithm allocated 6735 axon models (see figure 2).
To compute the spatial distribution of the induced extracellular field, a volume conductor model was developed following the approach described in [10].In brief, the heterogeneity of conductivity was obtained from patient-specific MRI segmentations, while tissue anisotropy was assessed based on the normative diffusion data [87], warped into patient (native) space using routines implemented in Lead-DBS toolbox [39].All calculations were performed within the simulation platform OSS-DBS [11], specifically developed for iterative DBS modelling that encompasses geometry and mesh generation, electric field computations and quantification of the induced neuronal activation.The accuracy of the field modelling was ensured by an adaptive mesh refinement algorithm of OSS-DBS that assessed the convergence of the electric potential on the axonal compartments.

Computing optimal electrode position using the per cent activation of the neuronal network model
The optimised pathway per cent activation defined the goal function for the electric field optimisation where A i is the per cent activation in the ith pathway.
Here, we quantified the difference between activation patterns with the Canberra distance [48] that was previously shown to be effective in a correlation model for the PD motor symptoms' improvement [10].Furthermore, we rejected solutions that yielded per cent activation above 5% in the corticofugal tracts, which are implicated in motor contractions.The following DBS parameters were set as the optimisation variables: the coordinates of the implantation site, the tilt of the electrode in the sagittal and coronal planes, and the delivered currents through the electrode contacts.The coordinates of the implantation were bounded by a 4 mm box centred on the STN, and the sagittal and coronal tilts were varied between 29.2 • -45 • and 10 • -30 • , respectively [34].The DBS signal was defined as a 130 Hz 60 µs train pulse.Its amplitude was limited by 2.5 mA for concentric and 2.0 mA for segmented contacts.This not only ensured that the charge density per phase lies below the threshold of 0.03 mC cm −2 [3,53] but also restricted the optimisation algorithm to more energyefficient protocols.
In this study, an axon was considered to have been activated if the discrete event of an action potential initiation was observed in response to the altered extracellular field.Therefore, the standard derivative-based methods are less suitable for optimisation.Furthermore, the goal function and the distribution of the projections implicated a non-convex problem with multiple local minima.After preliminary testing, the generalised simulated annealing [84], implemented in the Python package Scipy [81], was chosen as the optimisation algorithm.The optimisation of the current was performed utilising the linearity of Laplace's equation.It allows for avoidance of the recalculation of the electric field for any new combination of currents/potentials for the given electrode position.The algorithm works in the following order for an electrode with N contacts: The computationally intensive operations, such as the iterative solution of the electric field problem and the cable equation for thousands of axons were efficiently parallelized on a 24-Core 3.8 GHz workstation with 256 GB of memory.The electrode placement optimizer converged after ≈34 days of continuous computations.

Simulating normal state
Parameters were tuned to simulate healthy conditions during the initiation of motor activity (as described in the supplementary material).Figure 3 depicts the neuronal activity of the nuclei under normal conditions.The STN, GPe, and GPi show sparse firing with low levels of correlation of neuronal activity between STN and the pallidal nuclei (except for one rather synchronised firing period, which intrinsically emerged in the network at approx.650 ms, followed by an approx.50 ms pause in the GPe, and to a lesser extent in the GPi).1: Choose a new electrode position with the outer loop optimiser 2: For frequencies in the power spectrum of the DBS rectangular pulse, compute spatial distributions of the electric potential for each contact, setting this contact to 1 V, others to floating potentials and grounding the exterior of the computational domain 3: Apply Inverse Fourier Transformation (IFT) to obtain N solutions ϕ 1V (r, t); 4: Find the optimal currents by scaling the solutions with vector S and superpositioning (i.e.solve the optimisation problem for S): • If available, choose S from the previous iterations as the initial guess • With the scaled potential ϕ scaled (r, t) on the axonal compartments, compute the pathway activation and evaluate equation ( 7); 5: Return the result to the outer loop optimiser Generally, STN neurons form an almost random pattern and fire irregularly.Neuronal activity usually consists of spiking (B) and rare bursting activity (see figure 3(B)).The Fourier analysis of the macroscopic mean voltage V STN (see equation ( 3)) shows a wide, nearly stable spectrum, with a small peak at around 20 Hz and a drop from 500 Hz onwards (expected to reflect intercalated firing of different neurons, with ever smaller inter-spike intervals becoming increasingly unlikely; see also supplementary material).
The GPe-GPi activity is plotted in figures 3(b) and (c).Like the STN activity, it is also characterised by sparse firing, while individual neurons can show both bursting and spiking activity (b).The frequency of the mean voltage V for both GPe and GPi shows peaks in low β-band, underlying weak rhythmicity at approx.15 Hz.
The resulting time-dependent activity of neurons in the VA thalamus and MC is shown in figures 3(d) and (e), respectively.In contrast to STN and GP, both thalamic and motor areas show a dense spiking activity in γ-band with a peak at 48 Hz, with a broader range of γ-band (between 45 and 60 Hz) for the MC (see also supplementary material).

Simulating Parkinsonian state
In PD, the degeneration of nigrostriatal dopaminergic neurons leads to a loss of dopaminergic innervation in the striatum.The resulting reduction of D1/D2 receptor-mediated activity affects direct/indirect pathway functionality.In the direct pathway, reduced D1-receptor activation results in an increase of the GPi neuronal activity as a consequence of disinhibition, which in turn, results in higher levels of inhibitory activity in projections to the thalamus.In the indirect pathway, the reduction of suppressive D2-mediated receptor activation leads to an inhibition of GPe, thus enhancing STN activity (note that at this stage, we do not consider feedback activation of the STN from the GPe, which actually would lead to unphysiological delayed STN inhibition).The overactive STN will enhance neuronal activity in the GPi even more-which again leads to even more pronounced thalamic inhibition.
Consistent with disturbed SNc-striatal pathway activation, in the model, we increase the conductances g STN , g STN→GPe , g STN→GPi , while decreasing the internal GPe conductance.Simultaneously, we increase the inhibition from GPi to the thalamus by increasing the synaptic conductance g GPiVA (see supplementary material, section 3). Figure 4 shows the network dynamics under Parkinsonian conditions.The main emergent characteristic of the BG network activity is the synchronous activity, see figures 4(A) and (B), particularly in STN and GP, peaking at 24 Hz.The neurons, in turn, often fire in bursts (B).
The increased GPi activity is projected directly through the network to the VA and VL regions of the thalamus.An immediate consequence of this projection is that the ventral thalamic neurons change from regular spiking (healthy behaviour, at approx.47 Hz) to bursting activity (at approx.8 Hz).Due to the complex network structure of the thalamus, the ventral thalamic bursting activity, in turn, is projected to the VA thalamus through remote inter-thalamic connections [65].For motor-cortical activity, this results in a slight slowing, reaching a peak now at 37 Hz, rather than 47 Hz under healthy conditions.
The time-dependent activity of all neurons in the thalamus and MC is shown in figures 4(d) and (e), respectively.As already described, both thalamic and motor areas show an altered power spectrum compared to the healthy conditions: The power spectrum of thalamic activity peaks at low frequency, i.e. 8 Hz (low α), resulting mainly from the silencing of the thalamus and the emergence of low-frequency bursting.The spectral analysis of the MC shows a low γ to high β band activity, with the main component shifted to 37 Hz (from 47 Hz in the healthy state); see the supplementary material.

Simulating DBS using optimal per cent activation of the network
In this section, we present the optimisation results for the network interaction dynamics based on the minimisation of the biomarker of the thalamic activity defined in section 2.6 with equation (6).For this, the network structure and the model parameters reflect the Parkinsonian state.The optimal values of the per cent activation found are: A MC→STN = 0.878 (from MC to STN), A STN→GPe = 0.464, A STN→GPi = 0.562,  A GPi→Tha = 0.052, see also figure 5. DBS induces dramatic alternations in firing dynamics in all three BG nuclei.In the STN, neurons fire tonically, and overall they show a significant peak in the Fourier spectrum at 43 Hz, see figure 6(a)-(III), with further dampened peaks at harmonic frequencies, as reported by [78].Looking more specifically at individual neurons (all identical in properties, but heterogeneously connected), four classes of response behaviour can be differentiated in this nucleus: About 34% follow the stimulation at 130 Hz, showing additionally subthreshold depolarization (that might be attributed to either intrasubthalamic connections or GPe inhibition).Another 30% do not follow this rhythm but fire action potentials at 43 Hz.About 22% of the neurons remain unaffected by the stimulation and just fire single action potentials during the entire 300 ms period.The remaining cells stay utterly silent at hyperpolarised potentials (see section 5 and figure 1 of the supplementary material).
As a result of this activity in STN, GPe and GPi regions, neurons also show fast tonic activity again with the first component at ≈43 Hz, corresponding to one-third of the neurons (STN driving these cells).However, an even higher peak appears at ≈130 Hz, in resonance with the DBS and corresponding to fastfiring cells of the STN (another third of the neurons), presumably driving these cells.
As a consequence of STN DBS, thalamic activity is reset to fire tonically, albeit not at the 47 Hz rhythm of the healthy state, but at a slightly faster pace  3).For comparison reasons, the power spectrum of the healthy conditions is also depicted with blue colour.Under DBS, the activity in STN, GPe and GPi remains highly synchronised, particularly between GPe and GPi, at even higher peak frequencies (43 Hz).The representative neurons in the thalamus and cortex, i.e.II (d) and (e), become tonic again, at higher frequencies (58 Hz), in the γ range, and thus assume an activity pattern similar to the healthy condition again.(58 Hz).In turn, neurons in the MC revert to firing faster, close to the second peak of the healthy state (58 Hz).So overall, the activity is approaching healthy conditions again.As outlined above, the mean activity V in the thalamus peaks at 58 Hz, with a wider activation range in the interval  Hz, which shows a broad γ-band activation close to a previous healthy condition.Analogous behaviour is presented in MC with a characteristic γ-band activation and the highest peak at 58 Hz.These findings underscore that during DBS, in all nuclei, there is a bandwidth shift to the γ range, overwriting the pathological β synchronisation of Parkinsonian state, see figure 4.

Matching the optimal per cent activation in the volume conductor model
The question addressed in this study is whether an electrode position and current protocols can be found that will approximate optimised per cent activation A optimised predicted by our network model.While it may be possible to disrupt pathological activity by optimising stimulation patterns in the network model, it might still be problematic to translate this into feasible stimulation settings if the pattern is not reproducible with standard DBS systems.
The electrode placement and current optimisation allowed us to reach nearly the same per cent activation (see figure 9) using the directional lead (Boston Scientific Vercise TM , Marlborough, MA, USA), except for the hyperdirect pathway, the activation of which was limited by the detrimental capsular stimulation.

Discussion
This study proposes a biomarker-based optimisation paradigm for DBS treatment by combining complex spatio-temporal patterns emerging in large-scale networks with modelling of pathway activation.Rather than optimising a volume of tissue activated, the combined model extends this approach and aims to optimise realistic spatio-temporal discharge patterns in the motor nuclei of the thalamus.This directly associates the network's structural connectivity and the inherent dynamics to trajectories of white matter fibres, whose activation can be quantified by PAM.It allowed us to evaluate different electrode positions and stimulation protocols with respect to pathological bursting vs. spiking thalamic activity, serving as biomarkers.
The large-scale network model produces realistic spatio-temporal patterns with enhanced β activity (24 Hz, see figure 4) in the BG during the simulation of the Parkinsonian state.This pathological emergent behaviour is in good agreement with experimental findings in animal models of PD and with increased burst firing in patients [16,43,75,76].Specifically [76], analysed eight patients with PD and DBS and showed that pathological symptoms of bradykinesia, dyskinesia and tremor were accompanied by STN multi-unit activity showing a β-peak (≈20 Hz ± 6 Hz) [76].
Additionally, in our model, not only the frequency spectrum of the BG activity is shifted in the Parkinsonian state, but also that of the thalamic activity: the latter shows a slowing of activity, with a shift in the spectrum, from γ-to slow α-band (8 Hz; see figure 4).These emergent spatio-temporal dynamics of the thalamus are characterised by decreased thalamic tonic spiking and increased bursting activity.
These results are confirmed by several studies in patients and human and animal models [27,44,56].Both the reduced thalamic activity (i.e.reduced firing rate) [27,56], and, more specifically and importantly, an increase of power in the 3-13 Hz band (mean at 8 Hz, i.e. in the low α-band) with the reduction of γband spiking activity in VA neurons of the thalamus was shown in Parkinsonian nonhuman primates [44].In the current study, simulating the neuronal network with Parkinsonian conditions, the number of bursting periods of VA thalamic neurons increased as a result of the enhanced BG inputs (see figure 8(B)).An increased number of bursts was also reported in the study of Kim et al [46]: using an optogenetic mice model in which specific BG outputs were selectively perturbed, multiple Parkinson-like symptoms could be evoked.Specifically, Kim et al [46], using light trains of 20 Hz in GPI (imitating Parkinsonian β rhythm), induced high-frequency muscle and lowfrequency tremor activity.The authors showed that these symptoms depended on the number of thalamic neurons that exhibit bursting behaviour.
The optimal per cent activation of fibres under DBS (computed by our method) re-establishes higher frequency dynamics in the BG network.Thus, the activity peak in STN, GPe and GPi shifts from β-to γ-band, i.e. from 24 to 43 Hz (see figure 6).Similarly, the thalamus shows an activity shift from 8 Hz (α) to 58 Hz (γ), and in the MC, this results in a shift from 37 to 58 Hz (i.e. from low to higher γ).DBS simulations at optimal activation rate showed that a proportion of STN neurons depict unusually wide action potentials during DBS (33% of STN neurons).We assume that these neurons experience unphysiologically strong activation during DBS.Although there are no measurements of action potentials in vivo under these conditions, our model predicts that under such conditions, calcium currents are activated, which are strong enough to result in underlying [Ca2+]-spikes.

GABAergic synaptic plasticity effects during DBS
The simulation highlights the importance of the depression of GABAergic activity (either by direct silencing or desensitisation) for the efficacy of DBS.For the latter, a postsynaptic loss of sensitivity for the ligand still bound to the receptor is well described for GABAergic synapses [8].In addition, a depletion of the readily releasable transmitter pool as a presynaptic loss of synaptic efficacy is possible.This depletion of the transmitter has been described for hippocampal GABAergic interneurons, and it is dependent on the firing history of the synapse and, as a consequence, on the frequency of applied DBS.Thus, at high frequencies (>100 Hz), the transmitter release is suppressed shortly after the onset of the stimulus [57,82].It is possible that also pallidal neurones show this behaviour.
Our model proposes that precisely this suppression (be it via desensitisation or depletion) is one part of the DBS mechanism: high-frequency DBS induces loss of GABAergic synaptic efficacy, leading to a disruption of the pathological signal transmission through the pallidothalamic network.In figure 7, we show three single thalamic neurons under DBS, but under two different conditions: In the first (blue traces), desensitisation of GABAergic synapses is taken into account, and in the second, not (black traces).These examples demonstrate that when desensitisation of GABAergic pallidal projections is present, thalamic neurons continue to fire tonically.Otherwise, the majority of cells in the motor thalamus actually fall silent figures 7(A) and (B) due to the ongoing inhibition.Only a small fraction escapes it and continues to fire figure 7(C).

Two possible mechanisms of DBS: pre-and postsynaptic GABAergic silencing
The modelling of DBS suggests two possible mechanisms of DBS acting in parallel, as shown in figure 8.In this figure, the critical importance of the level of pallido-thalamic inhibition is demonstrated.Under healthy conditions, this level was simulated as extremely low (close to 0.01 pA), with minimal fluctuations, assuming a locomotive activity of the network.In essence, this releases thalamic neurons into firing tonically, and the firing frequency is finetuned by minimal changes in inhibitory drive figure 8 (A)-(I)−(III).Under Parkinsonian conditions, the GPi becomes much more active (due to disinhibition in the direct pathway and due to facilitation in the indirect pathway), resulting in much higher pallidothalamic inhibition levels (0.1-0.4 pA; baseline and peaks, respectively).This silences thalamic neurons for most of the time and only occasionally allows for burst-like behaviour in figures 8 (A)-(I)−(III), particularly after rhythmic and strong inhibitory episodes (B)-(II).Under these conditions, in all examples, the inhibition level is generally much higher than under healthy conditions, varying between 0.1 pA as basal level, and 0.4-0.6 pA at peak inhibition.As a consequence, thalamic activity is strongly reduced, and only intercalated firing can be observed.While in principle, intermittent bursting of thalamic neurons can be seen in all cases, the higher the rhythmicity and organisation of inhibition become (see example 2 with very rhythmic inhibitory drive), the more pronounced the rebound burst firing in thalamic neurons upon release from inhibition.(C) Under conditions of DBS, two different scenarios lead to restoration of thalamic firing: in example 1, inhibition is low but still present, indicating that pallido-thalamic projections are active, but the response is desensitised postsynaptically.As a result, the thalamic neuron fires tonically at a lower frequency.In example 2, a second mechanism takes over: In this case, inhibition is completely lost, possibly due to presynaptic silencing of GPe (due to GPi projections, for example).Again, this results in tonic firing, but now at a higher frequency.Example 3 shows that DBS is not acting perfectly: Here, the inhibition level is not very low (neither desensitisation nor silencing work perfectly), and as a consequence, the thalamic neuron is mainly silent, as under Parkinsonian conditions.The optimal electrode position was found in the motor aspect of the STN (highlighted in black).All contacts of the directional lead were active with both polarities, but the amplitude on a single contact did not exceed 1.3 mA, and the total current summed up to 0.9 mA.The obtained per cent activation deviated by less than 10% from the network-based optimal rates, except for the hyperdirect pathway, where the targeted rate was 30% higher.The reason for the limited activation is the spread of the current to the corticofugal pathway (not shown here); in the model, its per cent activation was capped at 5% to avoid capsular side effects.GPa and GPm-associative and motor portions of the globus pallidus, respectively.(B) Power spectrum of the VA/Tha activity in three cases, i.e. in a healthy state (black), under conditions of DBS, using thalamic biomarker to optimise network behaviour (red, network), and under conditions of DBS using electric field optimisation around the electrode (blue, electric field).(C) Power spectra of GPi in three cases, i.e. in a healthy state (black), under conditions of DBS using the thalamic biomarker to optimise network behaviour (red, network), and for per cent activation achieved by electric field optimisation (blue, electric field ).The power spectra for the latter case is computed in the network using, however, the values of per cent activation that are depicted in A.
Looking at this example more closely, the burst firing occurs in the wake of an inhibitory event but before the inhibition has entirely ceased, that the inhibitory hyperpolarisation is instrumental in initiating the burst.Under DBS, the high-frequency activation of the STN efferents also drives the GPi to higher activity.Nevertheless, pallido-thalamic inhibition levels drop to levels close to the healthy situation.Why this?As already explained, one mechanism of this reduction of the inhibitory drive is actually due to postsynaptic desensitisation [8] at high-frequency pallidal firing-an example is seen in figures 8(C)-(I).Besides this postsynaptic process, also a presynaptic one appears to occur: a complete silencing of inhibition due to DBS.One possible explanation is that massive STN activation due to DBS also activates the excitatory STN-GPe pathway (or, actually, the pathway is activated directly), and via this route, also the inhibitory GPe-GPi pathway, which would silence at least some of the GPi neurons.Such an example is seen in figures 8(C)-(II).The simulation also allows for a third insight: DBS is not always perfectly effective.In figures 8(C)-(III), the pallido-thalamic inhibition escapes suppression and remains nearly as active as under Parkinsonian conditions (in fact, pallidal afferents appear to fire in the range of [60 − 70] Hz, two times lower than DBS frequency (130 Hz)).As a consequence, the thalamic neuron falls silent, also because the inhibitory firing frequency appears to be too low to allow for postsynaptic desensitisation.

Optimised electrode position closely matches the clinically efficient placement for motor symptoms
The optimisation results for the current and electrode placement demonstrated that a conventional DBS system could induce activation patterns similar to those proposed by the network to restitute VAthalamic activity close to healthy conditions.A complete restitution is not reached since (a) DBS induces artificial activity that does not match physiological conditions and (b) the anatomical constraints limit stimulation (e.g.due to the proximity of the corticofugal tract).The optimised electrode position closely matched the clinically efficient placement in the motor aspect of the STN [18].It should be emphasised that the computational network operates with abstract activations and can generate physiologically implausible solutions, which need to be further evaluated with the biophysical model (see figure 9(A)).In this case, the suggested hyperdirect activation level was high.As was shown in [41], activation of the HDP is also highly correlated with activation of the corticofugal tract, which is associated with motor contractions, dysarthria and other side effects.This was taken into account for the current optimisation in our biophysical model, where the HDP activation reached 57.4%, comparable to the values reported in [41].The per cent activations computed in the biophysical model were again tested in the network model and created comparable VA dynamics (see figure 9(B)).
Furthermore, the optimised currents used both polarities, resulting in a strong local electric field, with only a 0.9 mA flow to the remote grounding, i.e. the casing (see figure 9).Importantly, one should bear in mind that various modelling approximations and uncertainties affect the electric field computations.Hence, the provided values should be rather used as the initial setting, further adjusted depending on electrophysiological recordings and patient response.

Model's limitations
The present study constitutes a computational approximation of the BG-thalamo-cortical network with the following assumptions and limitations.The connectivity structure was based on normative atlases; however, the methods can be adapted using individual patient tractography, replacing the network adjacency matrix A of section 2.3.1.The synaptic coupling was tuned according to the theory of direct indirect pathways' activation in Parkinson's to produce beta-band oscillatory activity within the BG as a pathophysiological marker of PD.Further, intrastructural connectivity in the nuclei was assumed to take the form of small-world complex structures.This novel approach in BG modelling is reasonably justified in previous modelling and experimental publications [4-6, 17, 24, 58, 66].As a limitation of the model, the exact structure of the connectivity on this microscopic level is not known.Hence, it remains to be clarified in the future how this can be analogously modelled.In the current study, the striatal input to GPe and GPi was simplified as homogeneous constant currents to all neurones in the nuclei but with different values between GPe and GPi, as well as between the healthy and Parkinsonian conditions.
The approach for the modelling and quantification of pathway activation from the network model in this study differs from [25].For the current approach, as explained in the methods section 2.6, neurons in the nucleus (STN, GPE etc) are modelled as singlesoma-neuron; however, in the computation of pathways activation, we used a multicompartment axon model.Based on experimental findings [27,44,46,49], in the current study, we introduced new biomarkers (bursting and spiking behaviour of thalamic neurons, inspired by the experiment of [46]), thus focusing on restoring thalamic functional, dynamic behaviour.Other possible options for biomarkers could be envisaged, for example, the local field activity of the MC or BG, which would focus on β-band suppression [25,43,63,71].Finding optimal biomarkers to guide DBS, with the clinical scores as the outcome measure, would be the gold standard in future studies.
From a future modelling perspective, one important topic will be the study of dynamic network variations with respect to DBS parameters (targeting, pulse width, frequency, shape, etc).Functional brain partitioning, which can be reached by combining model dynamics and network structure, may offer a computational approach and a systematic method that connects the clinical score with the biomarkers used for DBS optimisation.

Figure 1 .
Figure 1.Basal ganglia-thalamo-cortical circuits.(A) Circuit model comprising all main connections.MC/PMC-motor and premotor cortical regions, respectively.(B) Simulated reduced circuit model with highlighted (bold arrows) connections affected by STN-DBS.The synaptic connections between Striatum and Globus Pallidus pars externa/interna (GPe and GPi) were modelled by different constant currents.(C) Structural connectivity of the simulated network is based on the pathway atlas of the human motor network constructed from multimodal data, including diffusion, histological and structural MRI data, fused to a virtual 3D rendering [62].(C1) Projections from GPi to the thalamus (VA nucleus) are shown in red, connections between the nucleus subthalamicus (STN) and GPe are shown in green, and projections from STN to GPi are shown in violet.(C2) Connections between the motor cortex (MC/PMC) and the thalamus projections were obtained by classifying fibre tracts from [50] and are shown in orange.Projections from the motor cortex to STN (hyperdirect pathway) are shown in blue.Nuclei are shown in the following colours: GPe in light grey, GPi in dark grey, STN in dark orange, thalamus (VA nucleus) in yellow and motor cortex (M1) in light grey.Adapted from [68].CC BY 4.0.

Figure 2 .
Figure 2. Pathway activation modelling.(A) Extracellular electric potential distribution on computational axons is computed in the volume conductor model for a given current protocol.Here the electrode is shown at the position that was optimal for the normalisation of the VA-thalamic activity.Note that in order to reduce the computational effort, the axons do not cover the whole length of the pathways but only segments exposed to the high electric field.(B1), (B2) Activation status of the network pathways and other projections (corticofugal and associative pallidosubthalamic), respectively.By solving the cable equation for the computed distribution, responses of passive axons to the stimulation are quantified.If the axon elicits an action potential, it is considered activated.Note that some axons traverse the electrode and/or the encapsulation layer; those are considered 'damaged' and excluded from the simulation.(C) Electrode tip placement is optimized within a 4 mm bounding box centred at the motor aspect of the STN.

Figure 3 .
Figure 3. Representation of the whole network dynamics under healthy conditions.Each row represents the indicated area (STN, GPe, GPi, Tha/VA and MC).Under healthy conditions, there is only loosely correlated activity in STN, GPe and GPi, with tonic activity in the thalamus and motor cortex, at frequencies peaking at 47 Hz, i.e. in the γ range.(I) Raster plot representation.Black dots represent activated neurons (i.e.action potentials defined as transients passing V = −15 mV to positive values) against time (in ms) and space (i.e.neuron index in the nuclei).(II) The middle column depicts the time series of one representative neuron in each of the nuclei.(III) The right column shows the power spectrum of the mean voltage activity; see equation (3).The axes in III are on a logarithmic scale.

Figure 4 .
Figure 4. Representation of the network dynamics under Parkinsonian state.Each row represents the indicated area (STN, GPe, GPi, Tha/VA and MC).Under Parkinsonian conditions, there is a marked increase in synchronisation among STN, GPe and GPi, with frequency peaking at 24 Hz, i.e. within the β range.Altogether, this results in a lower thalamic firing frequency (major peak at 8 Hz, and a secondary one at 37 Hz, i.e. in the low β and low γ ranges).(I) Raster plot representation.Black dots depict activated neurons (i.e.action potentials defined as transients passing V = −15 mV to positive values) against time (in ms) and space (i.e.neuron index in the nuclei).(II) The middle column depicts the time series of one representative neuron in each of the nuclei.(III) The right column shows the power spectrum of the mean voltage activity see (3).For comparison reasons, the power spectrum of the healthy conditions is depicted with blue colour.

Figure 5 .
Figure 5. DBS-induced per cent fibres activation, optimised in the network using the biomarker of the VA-thalamic activity.The optimal values of the percentages of activated fibres found are: A MC→STN = 0.878 (from MC to STN), A STN→GPe = 0.464, A STN→GPi = 0.562 and A GPi→Tha = 0.052.

Figure 6 .
Figure 6.Representation of the network dynamics under DBS conditions using the optimal per cent activation values of section 3.3.Each row represents the indicated area (STN, GPe, GPi, Tha/VA and MC).(I) Raster plot representation.Black dots depict activated neurons (i.e.action potentials defined as transients passing V = −15 mV to positive values) against time (in ms) and space (i.e.neuron index in the nuclei).(II) The middle column depicts the time series of one representative neuron in each of the nuclei.(III) The right column shows the power spectrum of the mean voltage activity; see equation (3).For comparison reasons, the power spectrum of the healthy conditions is also depicted with blue colour.Under DBS, the activity in STN, GPe and GPi remains highly synchronised, particularly between GPe and GPi, at even higher peak frequencies(43 Hz).The representative neurons in the thalamus and cortex, i.e.II (d) and (e), become tonic again, at higher frequencies (58 Hz), in the γ range, and thus assume an activity pattern similar to the healthy condition again.

Figure 7 .
Figure 7. Single unit activity of three neurons in the ventral anterior thalamus under STN-DBS.In all cases (A)-(C), the blue traces depict neuronal activity considering GABAergic desensitisation, and the black traces show neuronal activity with non-desensitising GABAergic inhibition.If desensitisation is present, thalamic neurons continue to fire tonically in the γ-band in all instances (A)-(C).If desensitisation is not occurring, most of the neurons are actually inhibited and silent (A) and (B), and only some continue to be active (C).

Figure 8 .
Figure 8.Time series representation of activity (right ordinate, red traces) of three different thalamic neurons in the network (three rows) under normal (A), Parkinsonian (B) and (C) DBS treatment.The total synaptic current (from GPi to THA) is given on the left ordinate.(A) Modelling the healthy state, the three examples show the tonic firing of the thalamic neuron, whose frequency depends on the inhibition level.Although this is generally low, there is still some inhibition present in examples 1 and 3. Still, virtually none in example 2. (B) Simulating the Parkinsonian state, essentially three similar behaviours are observed:Under these conditions, in all examples, the inhibition level is generally much higher than under healthy conditions, varying between 0.1 pA as basal level, and 0.4-0.6 pA at peak inhibition.As a consequence, thalamic activity is strongly reduced, and only intercalated firing can be observed.While in principle, intermittent bursting of thalamic neurons can be seen in all cases, the higher the rhythmicity and organisation of inhibition become (see example 2 with very rhythmic inhibitory drive), the more pronounced the rebound burst firing in thalamic neurons upon release from inhibition.(C) Under conditions of DBS, two different scenarios lead to restoration of thalamic firing: in example 1, inhibition is low but still present, indicating that pallido-thalamic projections are active, but the response is desensitised postsynaptically.As a result, the thalamic neuron fires tonically at a lower frequency.In example 2, a second mechanism takes over: In this case, inhibition is completely lost, possibly due to presynaptic silencing of GPe (due to GPi projections, for example).Again, this results in tonic firing, but now at a higher frequency.Example 3 shows that DBS is not acting perfectly: Here, the inhibition level is not very low (neither desensitisation nor silencing work perfectly), and as a consequence, the thalamic neuron is mainly silent, as under Parkinsonian conditions.

Figure 9 .
Figure 9. Optimised electrode placement for the suppression of thalamic rebounding.(A) The optimal electrode position was found in the motor aspect of the STN (highlighted in black).All contacts of the directional lead were active with both polarities, but the amplitude on a single contact did not exceed 1.3 mA, and the total current summed up to 0.9 mA.The obtained per cent activation deviated by less than 10% from the network-based optimal rates, except for the hyperdirect pathway, where the targeted rate was 30% higher.The reason for the limited activation is the spread of the current to the corticofugal pathway (not shown here); in the model, its per cent activation was capped at 5% to avoid capsular side effects.GPa and GPm-associative and motor portions of the globus pallidus, respectively.(B) Power spectrum of the VA/Tha activity in three cases, i.e. in a healthy state (black), under conditions of DBS, using thalamic biomarker to optimise network behaviour (red, network), and under conditions of DBS using electric field optimisation around the electrode (blue, electric field).(C) Power spectra of GPi in three cases, i.e. in a healthy state (black), under conditions of DBS using the thalamic biomarker to optimise network behaviour (red, network), and for per cent activation achieved by electric field optimisation (blue, electric field ).The power spectra for the latter case is computed in the network using, however, the values of per cent activation that are depicted in A.