Computational Modelling of Deep Brain Stimulation for Parkinson's Disease: A Critical Review

We review the use of numerical and computational models to explore deep brain stimulation for Parkinson’s disease (DBS PD). It is a review for the modeler and those interested in PD DBS modelling methods and their value. The main model categories of active fiber, mean field, driving force, and volume of tissue activated are described as well as many modelling techniques. We give the basic requirements for a DBS PD model and current theories of DBS mechanism of action, PD etiology, and movement selection. The emphasis


Introduction
In neuromodulation as a primary treatment for neurological disorders, a secondary goal has been to uncover mechanisms of action (MoA) of treatments to optimize neuromodulation efficacy [16]. Discovered in 1987, deep brain stimulation for Parkinson's disease (DBS PD), a largely reversible procedure with few side effects, was hypothesized to be an electrical 'lesion,' and later began replacing irreversible, surgical ablative procedures (lesions via thalamotomy, pallidotomy) and avoiding their side effects [17]. DBS PD has since proven efficacious to alleviate movement disorders in tens of thousands of implants and opened the door to over 20 neuromodulation treatments and brain targets for neurological disorders, notably for pharmaco-resistant patients, i.e. those for whom dopamine agonists is ineffective (see Table 1 in Udupa & Chen [18]). Attempts to identify the underlying mechanism of DBS as an electrical 'lesion' for movement disorders revealed conceptual naiveté [17,19], and were augmented by hypotheses that DBS suppressed pathological rhythms or entrained restorative rhythms [18]. Another mechanistic hypothesis is removal of an 'information block' [20,21]. However, identifying the mechanism of action of DBS has proven difficult [22][23][24] and modelling studies cite this fact as part of their motivation. The causes of DBS side effects and the means to eliminate them likewise have not been fully uncovered [25]. Estimating parameters accurately is difficult in many biological settings [4], which can be offset by computational modelling [8,26]. One approach to address large error bars is to model parameter sweeps to see behaviour across the plausible or reported ranges of a parameter and to perform sensitivity analyses to determine which parameters affect a given result.
A desirable goal of models is being able to capture only the phenomena of interest, with the possible limitation of missing causative components. This Occam's Razor approach has evolved from experiences with complex models whose development time increases exponentially with complexity, and, especially in neural circuitry with recurrent connections, where results are difficult to decipher.
Efforts to uncover the mechanism of PD etiology now incorporate biochemistry, structural biology, molecular biology, and genetics, though to date no fully MoA-based treatment has emerged (e.g., L-dopa was originally tried for the wrong reasons) [27]. Computer modelling has assumed an increasingly important role due to many advantages such as the ability to explore phenomena inexpensively, quickly, in a large parameter space, and by capturing phenomena that would be difficult or impossible with animal or human subjects [28][29][30]. Numerical and computational models formalize a theory along with empirical findings incorporated into the model. Models are a test platform to explore hypotheses about the underlying MoA and predict the effects of new neuromodulation protocols on side effects and efficacy [31].
Accordingly, we describe numerical methodologies that are innovative and valuable additions to the modeler's toolbox. We include efforts toward closed-loop DBS (aka 'adaptive' or 'aDBS') due to its potential benefit to patientsand because aDBS imposes a rigid set of constraints on models that can aid successful modelling.

Literature Search Methods
We did a PubMed, Europe PMC, and Harvard Library Hollis system search with combinations of keywords 'deep brain stimulation, Parkinson's disease, movement disorder, circuitopathy, neuromodulation, model, computational model, and numerical model'. The inclusive criteria were: 1) The model is an exemplar of important methods or techniques, and 2) The model is transparent and repeatable. The latter criterion did not exclude complex models or difficult mathematics but did require clarity of exposition such that methods could be understood by others and reproduced in programming code. We did not emphasize the significance of modelers' results per se.
However, in an extensive Supplement we compile qualitative and quantitative empirical findings for model calibrations or validations and salient results to reproduce or compare to modelers' own results. The article attempts to help modelers assess tools and paradigms that hold the promise of uncovering PD DBS' MoA.
To date, DBS targeting the posterior subthalamic area and thalamic intermediate nucleus are done primarily for tremor vs. the primary PD targets of STN or GPi, and as such the former targets were not included in our search [32].

Brief Etiology of Parkinson's Disease
Since DBS PD is primarily used to treat PD movement disorders (MD), we focus on MD-relevant etiology here. PD, first described by James Parkinson in 1817, is an age-related neurodegenerative disease in which motor symptoms appear to be related to loss of dopamine-expressing neurons in the substantia nigra pars compacta (SNc, 'dark substance') and to a lesser extent in the ventral tegmental area (VTA). Lesioning the SNc in monkeys results in dopamine loss in the striatum and creates a parkinsonian-like state [27]. It is estimated from a rat model that ~80% of dopaminergic neurons must be lost for PD symptoms to appear, indicating a large amount of redundancy or signalto-noise ratio exists in the BG circuitry; 80% is ~400,000 of the ~500,000 dopaminergic neurons in human SNc/VTA [33]. From low levels in earlier decades, PD prevalence jumps to ~15% for patients over 60 years of age and increases 15-20% in each subsequent decade, constituting 1% overall prevalence over age 60 [34].
Neuropathological changes are observed at the tissue and cellular levels: neuronal loss in the ventrolateral cell group in the substantia nigra that projects to the posterolateral putamen, with accumulation of a-synuclein (Lewy Bodies) in those cells. Dopaminergic loss is thought to be caused by accumulation of pathological alpha-synuclein or mitochondrial dysfunction, or both [27,35].
The basal ganglia nuclei are strongly coupled, and effects of dopamine depletion on the striatum neural receptors are surprisingly complex and subtle [8]. Dopamine depletion correlates with increases in beta band power spectra (13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30) in the BG-cortical loop (see Supplement Table S4), which is also predictive of bradykinesia, and although beta band is a channel used for motor cortex to send instructions for weak-to-moderate contractions to muscle [36], a causal relationship between increased beta power and PD MDs has not been established [30]. The difficulty of sorting out cause and effect in PD etiology suggests elucidation via modelling, yet many studies seem to assume a causal relationship [37]. At the neuronal and circuit systems levels, attempts to correlate observable symptoms of PD with progressive and graded dopaminergic neuron loss are worthwhile modelling goals, as early diagnosis of PD is critical to improving pharmaceutical, electroceutical, or neuromodulation treatment efficacy as well as testing DBS PD MoA hypotheses [27].
An excellent short history of PD etiology is given by Goedert et al. [27]. Etiology of PD at the genetic signaling pathway level is summarized in the Supplement.

Theory of Movement Selection by Basal Ganglia and PD Etiology
Hypotheses on the MoA in PD MD center on the complex interaction within the BG, and between the BG, thalamus, and cortex ( Figure 1). Modelers focus on three main signaling pathways that have been proposed, called the direct, indirect, and hyperdirect. The original circuitry metric of healthy vs. dysfunctional states was neuronal firing rate ('rate model'), which has been increasingly augmented by resting state functional connectivity derived from MRI techniques (rsfMRI) [18].

Figure 1
Schematic of the basal ganglia (green), motor thalamus (yellow), motor cortex (blue), and brainstem and spinal cord (brown). There are two input nuclei, the striatum and the subthalamic nucleus (STN), comprised of: 1) Direct pathway: red dotted line from D1 striatal neurons to the globus pallidus internal nucleus (GPi); 2) Indirect pathway: red dashed line from D2 striatal neurons to the globus pallidus external nucleus (GPe); 3) Hyperdirect pathway: blue long-dotted line from cortex to STN. The basal ganglia output is GPi to motor thalamus to motor cortex, or GPe to brainstem to spinal cord. Deep brain stimulation (yellow lightning zigzags) is applied to BG nuclei STN or GPi.

Direct and Indirect Pathways
Synthesizing evidence from anatomy, physiology, and pharmacology, Albin et al. and DeLong et al. proposed a parallel processing action-selection model incorporating a direct and an indirect pathway comprising positive and negative feedback loops between the BG, thalamus and motor cortex, later empirically supported by Nambu et al. [12,13,38]. Functionally, the indirect, polysynaptic pathway from the BG striatum to the GPe and STN is hypothesized to globally inhibit a somatotopically-organized repertoire of motor programs triggered by GPi and SNr. This inhibition continues until the direct, monosynaptic pathway from the striatum inhibits the specific, somatotopic inhibitory group in the GPi/SNr targeting the desired motor program, thereby releasing it from inhibitory control and selecting it for action. The indirect pathway is a deadman's 'brake' signal (i.e. are always on brake) and the direct pathway is a 'go' signal [39]. Or in computer science terms, these circuits are a rapid-response multiplexer for selecting from a repertoire of action sequences. Focal injections of the GABA-A antagonist bicuculline (a chemical lesion) in the striatum results in loss of specificity in BG activity and twitches in specific muscle groups, demonstrating the somatotopic mapping [40].
The targets of the BG output nuclei GPi/STN are the motor regions of the thalamus, which targets the motor cortex, which sends motor program requests to the BG. The BG then processes and relays them to the brain stem and the spinal cord, where they are parsed into an ancient repertoire of motor programs and sent to the relevant muscles via efferent connections from the ventral cord [6-8, 13, 37, 41].
The largest component of the BG, the striatum (10 7 neurons), receives excitatory glutamatergic input from the motor cortex (the primary motor cortex, supplementary motor area, and premotor cortex) and is the main input relay center to the rest of the BG [42]. Via its medium spiny projection neurons (MSNs) that comprise 95% of it, the striatum projects inhibitory dopaminergic input to SNc and VTA [11]. The cortex requests specific motor programs (volitional movement control) via the striatum's fast spiking interneuron excitatory inputs (FSI, ~5% of cells), and the BG regulates the proper selection of each desired program before sending it to brainstem, which processes them and sends them to the spinal cord (automated movement control) [43]. The striatum has two groups of GABAergic MSNs, D1 and D2, characterized by two different dopaminergic inhibitory receptors [8]. Striatal D1 cells project monosynaptically to the excitatory BG output nuclei GPi and SNr, both GABAergic, inhibiting them in the direct pathway, while the D2 cells project polysynaptically to BG output nuclei via inhibition of GPe and STN in the indirect pathway [7,8,13,44]. Dopaminergic input from the SNc increases sensitivity to cortical excitatory input to D1 cells and decreases sensitivity of cortical input to D2 cells. Interestingly, MSNs can receive excitatory glutamatergic cortex input on their dendritic heads and inhibitory dopaminergic SNc/VTA input on their dendritic spines, permitting complex modulation [39]. Similarly, dopamine increases FSI activity by concurrently depolarizing the cell via its D5 receptor and reducing GABA input via the D2 receptor, and other neurotransmitters are involved [43]. The striatum's FSIs play an important role in regulating MSN projection activity via recurrent and feedforward inhibition, and MSNs receive inhibition as well from GPe Type A neurons [40].

Hyperdirect Pathway
The hyperdirect pathway, running from motor cortex layer V pyramidal neurons to STN, so-called by Nambu et al. since it bypasses the rest of the BG, is hypothesized to be the pathway through which ipsilateral motor actions are 'requested' by the cortex to the BG via excitatory, glutamatergic connections [11,12,45,46]. It has also been suggested to suppress inappropriate movement and fine-tune motor circuitry signaling (Vanduffel and Arsenault, 2016) [8,11]. The hyperdirect pathway has a latency of just ~1-3 ms compared with ~18-25 ms for the direct pathway, which runs through the GPi/SNr and thalamus to cortex [18,47]. While it is not known whether PD MD-related beta band peaks originate exclusively within the BG or require the cortex, DBS PD was hypothesized to affect the cortex via antidromic activation from the STN [48]. The computational model of Anderson et al. predicts that in high frequency PD DBS, if 15% + of the hyperdirect pathway is affected, antidromic APs suppress synaptic current, which 'lesions' the cortex-STN circuit, preventing generation or propagation of oscillatory activity [49]. However, recent work calls these conclusions into question (see below) [50].

Other Relevant Pathways
While most investigations of the origins of PD movement disorder (MD) focus exclusively on the basal ganglia (BG) (commonly the pallidal-STN loop), other circuitry, notably the thalamus, corticothalamic loop, and the BG-cortico-thalamic loop, may contribute to pathological BG processing and therefore are all open to test by modelling (see Figure 1 in van Albada et al. [8]). Likewise, it should not be ignored that DBS for non-MD symptoms such as pain may be valuable sources for DBS PD for MD model calibration and validation [51].
Modelers may incorporate compensatory mechanisms that reconfigure large-scale motor circuit regulation. Shine et al. found that integration of activity across the medial frontal, lateral parietal, and anterior temporal cortices was inversely correlated with motor symptom severity [52]. Similarly, a recent study pioneering the combination of several empirical techniques, and therefore possibly revealing new phenomena, shows that models may have to include cortical motor areas, reinforcing the significance of the hyperdirect pathway [53]. Ruppert et al.'s multimodal imaging showed reduced rsFC between a dopamine-depleted striatal putamen and cortical motor areas versus healthy state connectivity.

Current Theory of DBS PD Mechanism of Action at the Circuit Level
Four distinct and conflicting MoA have been proposed for DBS PD at the circuit level [23,24]: 1. Direct inhibition of neural activity 2. Direct excitation of neural activity 3. Information 'lesion', jamming, or blocking 4. Synaptic filtering Different MoAs may underly high versus low frequency DBS PD and in either case more than one MoA may be involved [23].
Direct inhibition of net STN or GPi activity was shown to not occur, although suppressing the excitatory sources of aberrant inhibitory interneurons that pathologically are disabled in inhibiting downstream hyperactivity could lead to depressed thalamic activity and akinesia ( Figure 2) [23].
Direct excitatory stimulation could work similarly, but conversely to direct inhibition, by increasing the activity of inhibitory interneurons to suppress pathological over-active STN or GPi neurons. Either MoA ultimately works by preventing DB-related synchronization or other pathology from propagating throughout the motor control network. Net effect is elevated activity through most of the circuit resulting in pathological synchrony. The difference between the diseased and health states could be measured, e.g. by changes in firing rates in the nuclei or power spectra, as could the degree of restoration of the healthy state via PD DBS. GPe: globus pallidus external. STN: subthalamic nucleus. GPi: globus pallidus internal BSt: brain stem. Thal: thalamus. After Humphrey et al., but with effects of depletion added [37].
A problem with the direct action and lesion hypotheses is that they may imply that normal motor control signalling may be impaired during DBS, yet studies show that sensorimotor signalling occurs in the pallidal centers under DBS to STN or GPi, as well as motor sequence learning and rewardbased reinforcement learning, which is another function of the BG [24]. These issues engender a more refined hypothesis that DBS filters out pathological activity while allowing normal motor control signals to pass.
Thus, the synaptic filtering hypothesis also solves the prima facie paradox of DBS PD stimulating higher activity in STN or GPi yet decreasing oscillatory activity: HF DBS PD over-stimulates cells to the extent that synaptic failure occurs, effectively reducing the frequency seen by axons. The resulting lower frequency acts as a filter that attenuates elevated pathological oscillations [54]. A second hypothesis is low-pass thalamic filtering, i.e. the thalamus passes healthy motor signalling (e.g. ≤β-band 13-35 Hz) while HF DBS PD (e.g. at twice the frequency in the circuit) blocks pathological activity [24,55,56].
A recent study by Johnson et al. indicates that the MoA of DBS is orthodromic, as DBS to GPi produced no antidromic activity in the motor cortex, while DBS within STN produced only acute motor cortex activity over a 4-hour period, although therapeutic benefit was maintained as activity waned [50]. Malekmohammadi et al. found that HF DBS (185 Hz, 90 μs) to GPi reduced low beta band power (13)(14)(15)(16)(17)(18) with no corresponding effect in the motor cortex, while high beta power was suppressed in both GPi and motor cortex [20]. Both studies draw attention away from hypotheses that DBS effects occur via the hyperdirect pathway, while not commenting on the role of the hyperdirect pathway in motor control. This work supports the orthodromic theory of DBS PD that, with SNc dopaminergic neuronal loss, GPe inhibition of STN weakens, producing elevated synchrony, in turn producing PD MD symptoms by blocking normal BG-MC signaling. DBS to dorsolateral STN or GPi reduces BG-thalamo-cortical synchrony and synchrony within the cortical motor areas, as indicated by beta-band power, thus removing the 'information block', restoring normal BG-MC signaling via thalamo-cortical projections [20,21,51]. If further studies reinforce this conclusion, it is a major clarification of DBS PD MoA and would help to focus modelers' efforts, beginning with reproducing the results of Johnson et al. and incorporating the required parameters as a calibration of models including the motor cortex.
The varied mechanisms proposed for DBS PD support Ashkan et al.'s call for changing the term 'deep brain stimulation' to 'deep brain modulation' to more accurately characterize DBS' MoAs [23].
We omit investigations into neurotransmitter MoA other than dopamine depletion as they are, as of now, typically incorporated into the overall inhibitory or excitatory connection strength in circuit-level modelling and not as individual channels. Even in active fiber models where there exist, e.g. a dozen different sodium and potassium channels, a 'lumped' approach is taken to capture the net effect of each ion channel type. However, as the science and practice of electroceuticals (neuromodulation replacing the efficacy of drugs), modelling detail at the neurotransmitter level (e.g. glutamate, gamma aminobutyric acid, calcium, serotonin, noradrenaline) may become important. The interested reader is referred to Stefani et al., McIntyre & Anderson, and Jakobs et al. [56][57][58].

Effects of Dopamine Depletion in Neuron-and Circuit-Level Models
From the movement disorder circuit-level theories described above, it is easy to see how, at least prima facie in PD, the loss of the inhibitory, dopaminergic input from SNc would result in dysregulation of the striatum, the largest BG component and dominant relay to the rest of the movement-controlling BG. To reveal specific MoAs, modelers must, however, take a deeper dive into the range of predicted effects [37].
Following  Table 1 and Table 2) [38,41].  At the circuit level, loss of dopamine results in elevated indirect pathway activity via the STN and weakened hyperdirect pathway activity. Together these results are theorized to prevent normal suppression of the indirect pathway by the direct and hyperdirect pathways, and thereby hinder proper movement program selection and modulation [46]. We detail the cascade of circuit effects from dopamine depletion in D1 and D2 in Figure 2 and Table 1. Figure 2 is representative of a simple schematic to illustrate an experiment to be modelled, where metrics such as firing rates in the nuclei or power spectra can be used to measure the difference between the healthy and diseased states and the degree of healthy state restoration via neuromodulation. It can be seen that the net effect of dopamine depletion at the circuit level is elevated excitatory activity. From this we hypothesize that a background signal-to-noise ratio that is exceeded by motor-program selection commands in the healthy state is replaced by a noisier state characterized by pathological synchrony, and that DBS PD interrupts the synchrony, attenuating noise back toward healthy state levels.
Further calibration data for a dopamine-depleted circuit is provided in Table 2.

Basic Requirements for PDS DB Modelling Studies
Augmenting a similar list from Holt & Nethoff [31], a conceptually complete DBS PD model must include the following: 1. Biomarkers to measure calibration and validation of healthy and diseased states 2. Well-defined descriptions of healthy and diseased states at appropriate systems levels 3. One or more neuromodulation protocols hypothesized to ameliorate a given pathology 4. Hypotheses of how each neuromodulation protocol affects neural activity and restores the healthy state from the diseased state 5. Possible in silico, in vitro, or in vivo experiments to verify or falsify the hypotheses

Model Calibration and Validation
An extensive set of criteria that can be used for calibration and validation of models is given in the Supplement. Any number of the given criteria can be used to calibrate a model, and then running the model to see if it manifests any of the remaining criteria can constitute a validation. In general, biological models tend to be under-validated. Model robustness benefits from incorporating calibration and validation data that are in some sense orthogonal to each other.

Types of DBS PD Computational and Numerical Models
Neural circuitry computational models typically involve a way to represent processing of inputs to neurons from other neurons and output from neurons to other neurons in a software program, some having more detail in such processing, some less, some more biologically inspired, some lessso. Modelling can capture detail at several biological systems levels, which affect methods of calibration, validation, as well as model complexity.

Top-Down Organizing Principles
Holt & Netoff surveyed computational and numerical models of epilepsy as well as DBS PD and classify DBS PD into two categories with different goals: (1) abstract models that replicate behaviours but are not clinically predictive, and (2) models that are clinically predictive, incorporating and solving for patient-specific parameters [31,59]. These are idealized categories, and most modelling is a mix of the two categories and one way or another embody a 'bench-tobedside' paradigm. A purely theoretical model typically must still incorporate clinical parameters to be useful and a purely clinical model whose aim is to improve clinical parameters typically still needs to be distilled down to an abstract model.
In the last two decades the top-down method of imaging correlated functional brain activity to create a high-level resting state functional connectivity connectome (rsFC) of the healthy and diseased states has become increasingly effective in analyzing many neurological disorders [60][61][62]. In rsFC the biomarker dimensions provide a way to describe healthy and diseased states and a 'distance' between diseased and healthy states. Testable hypotheses can be framed in terms of the degree to which neuromodulation restores the healthy state along one or more rsFC dimensions. The rsFC paradigm is likely to provide the over-arching framework into which DBS PD MoA and efficacious parameters will be organized in the future.

Mean Field, Driving Force, and Volume of Tissue Activated Models
Gunalan et al. 2018 sought to compare the accuracy of three categories of DBS PD models: 1) field-cable, which includes active fiber models coupled with FEM of the electric field, and which is most computationally complex and time-consuming, but also most accurate, 2) driving force, which is in the middle of the complexity spectrum, and 3) volume of tissue activated (VTA), which is least complex, and as they found, least accurate [77].
Neural mean field theory is a high-level, computationally efficient mathematical framework typically used for modelling large neural circuits. The term 'mean field' derives from the use of spatial and temporal averages in discrete, coarse-grained volumes and time intervals to describe distances between neural populations at spatio-temporal position and numbers of APs in a given time, e.g. firing rate Qa (where a is for average) at position r and time t: Qa(r,t) [78]. The dramatic reduction in degrees of freedom results in computational efficiency and transparency [79]. Rise and fall sides of the AP can be described by exponential, sigmoid, or step functions. A potential shortcoming is the use of just two neuron types, excitatory or inhibitory. Mean field theory's framework matches that of rsFMRI and describing the difference between healthy, diseased, and neuromodulated states [15]. We see significant advantages in marrying neural field models with more detailed bottom-up approaches that incorporate, e.g., different neural types as do Kumaravelu et al. [29,80].
VTA models typically incorporate several simplifying assumptions: 1) isotropic tissue conductivity, 2) straight axons oriented perpendicular to the applied electric field, 3) finite volumes of tissue proximate to the electrode, all for parameters of fiber diameter, specific electrode array activation, pulse width, and amplitude [77,81,82]. Note that assumptions 1 & 2 can be loosened for greater accuracy by using anisotropic conductivities in multiple tissues and axons of any orientation to the field. Also, while VTA commonly requires an additional assumption that AP frequency match stimulation frequency ('1:1 response'), we suspect this requirement is spurious and misleading, based on other neuromodulation paradigms (e.g. spinal cord stimulation and vagus nerve stimulation), in which efficacy can be achieved without 1:1 response.
The intermediate complexity category, driving force models, significantly reduces complexity of field-cable models by using a metric for activating an axonal AP [77,[83][84][85][86][87]. The inherent limitations of the simplified approach can be ameliorated by incorporating distance-weighted potential contributions to the excitation of a given node of Ranvier (2 adjacent nodes in Rattay's 'activating function' and an unlimited number in Warman's 'driving function'), or calibration to exogenous passive or active models [85,[88][89][90][91]. The accuracy gain of ~20% under certain circumstances may not be necessary for most model goals or as a first pass model.
Ultimately, DBS PD model types are strategies to efficiently bridge several biological systems levels (multi-scale modelling) by incorporating no more detail from each than is needed to capture the phenomena of interest, thus, models must be judged in the context of their goals and no method is 'best' [29]. Our categorization of the models reviewed here is given below.

van Albada et al.
A notable exemplar of modelling DBS PD is the paradigm of van Albada et al., which captured the direct, indirect, and hyperdirect pathways, 5 types of progressive dopamine loss, and sought to identify differences between healthy and PD states [8,92]. Their study is also an early example of using functional connectivity to distinguish the healthy and diseased states and map the model to known anatomical regions. The work exemplifies a balance of the qualities of 1) exemplary knowledge of PD and literature review for hypotheses to test and calibrations, 2) sophisticated numerical methods, and 3) a rich set of phenomena captured at the appropriate systems level, in this case using mean field theory for efficient simplification of a complicated system. Most modelling studies are under-weighted in one or two of those categories. Concomitantly, the downside of their framework is that it is on the complex side of the modelling spectrum and involves sophisticated mathematics, which limits access and reproduction of results to a subset of modelers. A reasonable use of the van Albada et al. work would be for modelers to create a simpler way to capture a salient subset of their paradigm, e.g. progressive dopamine loss, and compare results with those of the more complex model.

Gunalan & McIntyre
Gunalan & McIntyre sought to model DBS PD recruitment of the hyperdirect pathway, as suggested by Li et al., who noted that DBS targeting STN will generate axonal APs in both directions along the axon, i.e. orthodromically to STN's BG target, GPi, or antidromically back to the cortex, and cited supporting evidence of antidromic effects [48,93] (but see conflicting evidence in [50]). The Gunalyn-McIntyre model is quite sophisticated in a number of respects, such as the detailed MRI data underlying its design and calibration, the use of an active fiber nerve model that is, in itself, quite sophisticated [64], and the design of the neural circuitry. Notably, Gunalan & McIntyre conclude that the evoked potential (EP) ~3 ms EP1 and ~5 ms EP2 signals are correlated with antidromic activation of the hyperdirect pathway, and specifically with the 10 and 5.7 µm fibers in their model. While prima facie this seems significant, their model only included the three diameter axons, and these seem too large for the cortico-BG circuitry, which tracing studies specifically of the hyperdirect pathway show to be 2-3 µm in diameter (and the article describes them as 'large caliber' in the context of cortical axons) [46]. A further criticism is that their nerve model was calibrated to certain motor fiber behaviour such as threshold and AP shape, and the thresholds may be unreliable below the typical large diameter of motor axons. While these observations suggest possible flaws in the Gunalyn-McIntyre methodology, the model is sophisticated, original, and worth re-building with accurate calibration parameters or re-using its components in a different, and perhaps simpler, model.

Lempka Howell McIntyre
Integrated circuit pulse generators (IPG) are used to produce waveform shapes intended by DBS protocol designers, but what waveforms are actually produced? An influential paper by Geddes illuminated this question from the standpoint of voltage-driven vs. current-driven IPGs, showing that current-driven IPGs preserve rectangular waves with high-fidelity while voltage-driven IPGs do not, and consequently manufacturers converted some protocols from voltage-driven to currentdriven [94]. Lempka et al. analysed actual IPG fidelity specifically for protocols and waveforms used in DBS, simulating IPG circuits with simplified electronically-equivalent model circuits calibrated to bench IPG output recordings [95]. They found significant differences for long bipolar pulses, and average errors were ~3-24% for constant-voltage stimulation and ~2-11% for constant-current stimulation. Lempka et al. note the importance of using accurate waveforms in computational modelling DBS to predict efficacy, e.g. via metrics such as VTA. Another factor motivating manufacturers to adopt constant-current stimulation was its ability to preserve the strength and shape of the electric field at the target dorsal column as electrically-resistive scar tissue forms under the electrodes (fibrosis) [96]. Uneven fibrosis not only can weaken the field, but counter-intuitively, can strengthen the field and distort the intended distribution, such as when scar forms under a 'guarding anode' that is restricting the penetration of cathodic pulses [90]. These theoretical predictions were borne out in a study showing greater patient satisfaction and pain relief with constant-current vs. constant-voltage stimulation [97].

Lindahl et al.
Noting that previous computational and animal models had failed to uncover the global effects of numerous local effects in the BG of dopamine loss, Lindahl et al. attempted to model the BG in a larger, coarse-grained network context, calibrated with known local dopamine loss effects [40].
While coarse-grained (7 neural groups), Lindahl et al.'s spiking network model is nonetheless quite complex (80,000 neurons, 11 neuron parameters, 5 neuron types, ~100 synaptic parameters, ~10 6 connections, and various other parameters) and reproducing their results, re-purposing, or modifying the model would be a significant undertaking. Thus, the paper is a rich source of parameter values and calibrations for modelers.

Izhikevich
An insightful description and comparison of spiking neuron (aka 'integrate-and-fire') models, their benefits, computational efficiency, usage in network models, and tuning phenomenologically (i.e. to firing characteristics vs. to empirical neuron parameter values) is given by Izhikevich [98]. These models derive from the general theory of nonlinear dynamical systems, especially bifurcation theory, the essence of which is described. The term 'hybrid' refers to the two phases of simulation implemented with a piecewise algorithm: 1) computationally detailed, continuous dynamics of the prelude to neuron threshold and the upstroke of spike generation, and 2) computationally simple, discontinuous and linear 'hard reset' to initial conditions after a spike peak voltage defined by the modeler. A clever hack is used to transition time-stepping from a relatively coarse-grained step (typically 1 ms) leading up to a spike and a finer-grained timestep in the vicinity of the spike to capture the rapid spike dynamics. However, a trade-off occurs between the computational efficiency of using large time steps and numerical instability triggered by large synaptic currents. Instead of using computationally-intensive numerical methods (e.g. backward or 'implicit' Euler, higher-order Runge-Kutta) to solve the differential equation for membrane potential, Izhikevich sacrifices accuracy in the down-spike shape by using a linear equation.

Thibeault et al.
A considerably simpler and more manageable spiking neural network model to compare with Lindahl et al. is that of Thibeault & Srinivasa, containing 4 neuron types and groups, 5 neuron parameters, 5 synaptic parameters, and 12 sets of connections modulated by just 2 parameters [99]. Notably, the Thibeault & Srinivasa study is an exemplar of reducing model complexity in their critical selection of four previous hybrid models of the BG and their simplification of the models. They provide succinct overviews of a selection of earlier models and the hybrid spiking model comparing their different topologies. They analyzed the action-selection theory of BG and the PD diseased state but fall short of being able to offer improvements on DBS protocols.
While not using the term 'functional connectivity, ' Thibeault & Srinivas, in their usage of correlated firing rates in healthy vs. PD states, implicitly align with the functional connectivity paradigm. Although the model is considerably simpler than others (e.g. Lindahl et al.), Thibeault & Srinivas nonetheless largely succeeded in calibrating their BG neural firing patterns to empirical data, comparing their results to those of earlier models, e.g. Humphries et al. [41]. They argue that some BG firing patterns require enough complexity to model individual neurons as well as network-or population-level effects. This issue in general will remain central to ongoing discussions of realistic and efficient connectome modelling [29]. Thibeault & Srinivas' theme is that the use of simpler models may reveal mechanisms of action that are difficult to parse out in more complex models. Besides their Izhikevich hybrid neuron, they point to the adaptive-exponential-integrate-and-fire model [100], and secondarily, rate-and population-based models. They stress the unreliability of using hybrid neurons in small networks vs. large-scale ones and recommend biologically more realistic conductance-based models for simulating individual neurons or small networks.

Arle et al.
Along those lines, the hybrid UNCuS universal network simulator captures both individual neuron and network-level phenomena efficiently [101]. The model can simulate individual neurons with any characteristic current-voltage curve and networks containing them, including ion channel conductances, threshold accommodation, and connectivity strengths set globally, or locally along the dendritic tree, as well as in the density of source-target connections. A lumped nonlinear conductance parameter permits user-defined cell model expansion by using the parameter to tune to a given neuron's characteristics, store the neuron in the repertoire of model neurons, and call it as needed at the population level. Notably, incorporating computationally efficient features from MacGregor, Rall, and McNeal, their framework scales efficiently to model small-to-medium sized networks such as DBS PD [101,102] or as large as the human spinal cord containing 10 6 neurons [87,[103][104][105]. In contrast to models on the more complex side of the spectrum such as van Albada [8], the Arle et al. model is accessible and can be programmed by anyone familiar with numerical solutions of differential equations. The principal calibration for modelling PD DBS was background BG firing rates in humans and the key result was DBS 'resetting' firing rate variance in PD and predicting efficacy at stimulation frequencies > 130 Hz in agreement with clinical reports.

Kumaravelu et al.
Kumaravelu et al. offer specific analyses of how DBS works and what protocols are optimal, according to their model [80]. Noting that earlier computational studies replicated primate animal models of PD but were not calibrated to the 6-OHDA rat PD model, they used that calibration and wanted to find what this particular model showed about DBS PD to the STN. Interestingly, Kumaravelu et al. used 7 complex neuron models incorporating up to 11 different ion channel dynamics, and compensated for the lower-level complexity at the network level by using just 10 neurons per nucleus-striatum, STN, GPe, GPi (which in rat is homologous entopeduncular nucleus, EP), thalamus and cortex (containing regular-and fast-spiking neurons). Thus, their emphasis is 180 degrees opposite to rate-or population-based models, and their emphasis on where complexity is needed to capture the key DBS PD phenomena strikingly different than the hybrid spiking models. Despite the bottom-up complexity, the calibration to the diseased state is simple and transparent, involving just 3 neuron parameters capturing 3 effects of striatal dopamine loss: 1) reducing M-type potassium channel conductance, 2) reducing cortico-striatal conductance, and 3) increasing GPe synaptic conductance. The model of Kumaravelu et al., calibrated uniquely to the 6-OHDA lesioned rat model, showed, similar to findings in Arle et al. [102], that 130 Hz DBS PD disrupted a pathological synchrony in STN, GPi, and GPe nuclei. The notable feature of their study is its focus on a specific model and developing a model that is complex at the lower biological systems level offset by minimal complexity at the immediately higher systems level.

Muller et al.
The Muller et al. study is an exemplary, transparent use of mean field modelling to capture largescale network behavior where the underlying methodology is explicated in detail. They compare the pathological power spectra and synchrony changes resulting from DBS PD to the STN and GPi. The stated goal is to provide a framework to systematically test the effectiveness of different DBS protocols that has been validated against empirical and clinical findings. The focus on steady state firing rate perturbations around a mean value demonstrates the importance of tying mean field to those of rsfMRI.

Kang & Lowery
Kang & Lowery hypothesized that DBS PD disrupted pathological synchrony using beta band power as the metric [106] and sought to distinguish orthodromic effects on the STN from antidromic effects on the cortex, and examine the combined effects [106,107]. They present the empirical evidence from human and rat studies for combined effects, and contrast those with empirical and theoretical arguments that antidromic effects may be negligible or minor. Analyzing and separating out orthodromic and antidromic effects in DBS PD via simulating physical or chemical lesions, and sweeping across parameter ranges that are simply impossible outside of numerical methods, are good examples of what can be done relatively quickly and inexpensively with a computational model vs. in vivo studies [107].
Kang & Lowery employed numerous simplifications: three neural centers (cortex, GPe, and STN), a point electrode source, homogeneous conductivity,~100 hundred neurons/nucleus,random vs. topographic, connectivity, a point source electric field, and a volume conductor model [108]. At the neural level, however, they attempted to capture more detailed effects drawing on different model neurons from other labs for the different nuclei [106].
At the cortical neural level, the model followed Rattay [85] and was detailed anatomically, modelling excitable elements of the neural cell body, axon initial segment (AIS), main axon and axon collateral, while omitting the dendritic level, which may be important for the timing of AP reception (cf. Arle et al. Figure 2, Yousif et al., and Rall [101,105,109]). STN and GPe neurons derived from previous models developed specifically for those nuclei. The specialized STN neuron followed the conductance-based model of Otsuka et al., who calibrated to data on how STN is regulated by 'plateau potentials' manifesting during hyperpolarization of particular STN neurons [110]. The STN neuron followed the detailed, specific model of Rubin & Terman [111]. In these cases a simplified model capturing the essential STN neuron current-voltage behaviour would likely suffice, whereas in their previous study STN neuron detail was necessary to capture the effects of dopamine depletion [106]. Last, the computationally efficient Ishikevich spiking model was used for interneurons (see above description [98]).
Kang & Lowery's calibrations was top-down and notably simple. They initiated beta-band oscillations in the cortex, adjusted cortico-striatal synchrony via synaptic gain adjustment (a simple scaling factor) from cortex to STN, between self-inhibitory GPe neurons, and mutually between STN and GPe until beta synchrony manifested in the STN and GPe [107]. They then used computational 'lesions' to isolate the antidromic and orthodromic pathways under DBS and examine the resulting effects in each circuit separately and with both activated.
Kang & Lowery's key results are: 1) Orthodromic DBS entrains STN neurons to the stimulus frequency and suppresses pathological beta activity throughout the network; 2) Antidromic DBS achieves the same effect via a different mechanism, to wit, disrupting cortical activity that otherwise (per their calibration) causes or reinforces beta activity in the BG (STN and GPe), and 3) The effects are synergistic. The antidromic result predicts that the cortex is a primary target of DBS PD as well as STN, providing an alternative hypothesis to activation of STN alone and may illuminate identification of DBS side effects via axons of passage (cf. [25,26,93,[112][113][114]). HF DBS (60 μs pulse width @ 130 Hz) is predicted to reduce STN beta-band power via disruption in both directions and is roughly proportional to frequency between 80 and 130 Hz, while LF DBS would reinforce it via resonance effects. (Note that clinically, DBS PD frequency >130 Hz is more effective than lower frequencies and is always used except to treat dystonia.) DBS at 20 Hz had no effect in either direction.

Yousif Purswani et al.
A thesis of the study by Yousif et al. (2010) of DBS PD to STN was that a detailed model of neural anatomy including the dendritic tree, soma, and axon would reveal significant phenomena missed by previous studies [115]. Following Arle et al. [101,102], they coupled a finite element model to predict the electric potential distribution in space and time to their neural model. Notably, they compared the DBS metric, VTA, using their neural anatomical model to the typical axon models. Not unexpectedly, to directly fire APs by stimulating neuron cell bodies required over 10x the stimulus amplitude as did axons, but they also found that neurons are stimulated in high current density regions near the electrodes perhaps in sufficient numbers to affect DBS. Yousif et al. found significant variation in VTA across different DBS frequencies, and noted the importance of capturing the effects of DBS frequency as well as the resulting frequency of APs in targets such as the STN. Their statement that frequencies below 100 Hz worsen symptoms has been proven simplistic by subsequent studies such as those described in the meta-study by Su et al., which found differential effects in the HF and LF DBS ranges [10].

Closed-Loop, Adaptive DBS (aDBS)
We give an overview of the relatively new goal for DBS PD, 'closed-loop' technology, because it may supersede the current non-closed loop paradigm and modelling will play an important part in its development. Closed-loop refers to a technology wherein recordings of effects produced by a neuromodulation stimulus are used as feedback to modify the stimulus in real-time ('adaptive' or 'aDBS'). Besides the advantage of superior protocols tuned to optimize efficacy for different tremor types, a significant potential advantage is greater energy efficiency due to lower duty cycle (ON time to ameliorate MDs), and thus less frequent invasive surgical battery changes for PD patients [22]. Wang et al. make a strong case that 1) MoA understanding is necessary to understand neuromodulation side effect etiology, 2) mechanistic models are necessary to uncovering MoA, and 3) closed-loop protocols developed via the models are necessary to reduce side effects [116].
Closed-loop DBS PD is in an early stage. The closed-loop concept floated for DBS PD is to identify signatures for the several types of PD MDs-'biomarkers' or 'signatures', such as changes in LFPs in the STN or GPi-and then select from a menu of efficacious DBS protocols a response to each signature. However, neither has a reliable signature been identified for any DBS MD, nor has a specific protocol addressing a type of DBS MD yet been developed [22,117]. Modelling may help overcome these two critical and necessary hurdles. Developing protocol responses to specific MD signatures will lead to reducing side effects such as psychiatric symptoms and speech interference [118,119]. Along with improving the accuracy of MD signatures and developing efficacious responses, a third hurdle is to identify the signatures rapidly enough that aDBS algorithms can respond in real-time [118]. However, segmented leads may allow programming to eliminate most side effects without the need of a closed-loop component.
Biometrics based on beta-band power are being replaced by more precise metrics [120] and single signature components, e.g. beta-band power or high-frequency oscillations, have proven insufficient; detecting multiple signature features increases accuracy [121][122][123]. Hirschmann et al. found that training the machine learning algorithms on a single patient was more accurate than across patients, indicating signatures differ significantly between patients [121]. This finding further strengthens the need to use machine learning in individual aDBS devices and for modelers to encompass the range of a given signature in patient populations. Thus, modelers will need to simulate machine learning to keep up with aDBS technology.
Advances in the application of machine learning techniques may be the route to identifying signatures for tremor types, to which modelling techniques can be applied to develop the effective protocols and open the door to aDBS. The first human clinical trial comparing conventional DBS with aDBS has been planned [120]. The study acknowledges the limitations of beta band power as a biomarker for tremor and will focus on its correlation in STN and GPi with bradykinesia and rigidity, will compare STN versus GPi targeting, test the hypothesis that the intended low duty cycle of aDBS will reduce dysarthria (impaired speech) due to over-stimulation in conventional DBS, and test for side effects in executive function.

Conclusion
Computational modelling of DBS PD has extensively explored the main phenomena of interest in PD etiology, such as dopamine depletion and dysfunction of the direct, indirect, and hyper-direct signaling pathways. We have reviewed a representative sample of neural simulation methods that have been developed by modelers. While the work done to date has fallen short of decisive unraveling of PD etiology and DBS MoA, it provides a rich paradigm for selecting tools with which to target specific PD phenomena and DBS effects, suggest hypotheses to test, and incorporate refined empirical parameters. At a meta level, computational modelling of DBS PD is an example of the general need for modelers, theoreticians, and empiricists to couple their work more tightly. Empirical results should be rapidly incorporated into the models and model predictions should be more rapidly tested empirically [30].
Ultimately, the measure of model utility is the ability to predict therapeutic outcomes and improve protocols and devices. Modelers must judiciously weigh the level of model detail needed to capture diseased vs. healthy states, the etiology of the diseased state, and the restoration of the healthy state due to the applied therapy, vs. limitations, such as transparency and computational efficiency [29,99]. Looking to the future of individualized closed-loop PD DBS, portable, energyefficient hardware containing a model of the patient is unlikely to be feasible unless the model is simple enough to run computationally, to understand, and to validate that it runs correctly [99,124].