Virtual deep brain stimulation: Multiscale co-simulation of a spiking basal ganglia model and a whole-brain mean-field model with The Virtual Brain

Deep brain stimulation (DBS) has been successfully applied in various neurodegenerative diseases as an effective symptomatic treatment. However, its mechanisms of action within the brain network are still poorly understood. Many virtual DBS models analyze a subnetwork around the basal ganglia and its dynamics as a spiking network with their details validated by experimental data. However, connectomic evidence shows widespread effects of DBS affecting many different cortical and subcortical areas. From a clinical perspective, various effects of DBS besides the motoric impact have been demonstrated. The neuroinformatics platform The Virtual Brain (TVB) offers a modeling framework allowing us to virtually perform stimulation, including DBS, and forecast the outcome from a dynamic systems perspective prior to invasive surgery with DBS lead placement. For an accurate prediction of the effects of DBS, we implement a detailed spiking model of the basal ganglia, which we combine with TVB via our previously developed co-simulation environment. This multiscale co-simulation approach builds on the extensive previous literature of spiking models of the basal ganglia while simultaneously offering a whole-brain perspective on widespread effects of the stimulation going beyond the motor circuit. In the first demonstration of our model, we show that virtual DBS can move the firing rates of a Parkinson’s disease patient’s thalamus - basal ganglia network towards the healthy regime while, at the same time, altering the activity in distributed cortical regions with a pronounced effect in frontal regions. Thus, we provide proof of concept for virtual DBS in a co-simulation environment with TVB. The developed modeling approach has the potential to optimize DBS lead placement and configuration and forecast the success of DBS treatment for individual patients. Highlights - We implement and validate a co-simulation approach of a spiking network model for subcortical regions in and around the basal ganglia and interface it with mean-field network models for each cortical region. - Our simulations are based on a normative connectome including detailed tracts between the cortex and the basal ganglia regions combined with subject-specific optimized weights for a healthy control and a patient with Parkinson’s disease. - We provide proof of concept by demonstrating that the implemented model shows biologically plausible dynamics during resting state including decreased thalamic activity in the virtual patient and during virtual deep brain stimulation including normalized thalamic activity and distributed altered cortical activity predominantly in frontal regions. - The presented co-simulation model can be used to tailor deep brain stimulation for individual patients.


Highlights
-We implement and validate a co-simulation approach of a spiking network model for subcortical regions in and around the basal ganglia and interface it with mean-field network models for each cortical region. -Our simulations are based on a normative connectome including detailed tracts between the cortex and the basal ganglia regions combined with subject-specific optimized weights for a healthy control and a patient with Parkinson's disease. -We provide proof of concept by demonstrating that the implemented model shows biologically plausible dynamics during resting state including decreased thalamic activity in the virtual patient and during virtual deep brain stimulation including normalized thalamic activity and distributed altered cortical activity predominantly in frontal regions. -The presented co-simulation model can be used to tailor deep brain stimulation for individual patients.

Abstract
Deep brain stimulation (DBS) has been successfully applied in various neurodegenerative diseases as an effective symptomatic treatment. However, its mechanisms of action within the brain network are still poorly understood. Many virtual DBS models analyze a subnetwork around the basal ganglia and its dynamics as a spiking network with their details validated by experimental data. However, connectomic evidence shows widespread effects of DBS affecting many different cortical and subcortical areas. From a clinical perspective, various effects of DBS besides the motoric impact have been demonstrated. The neuroinformatics platform The Virtual Brain (TVB) offers a modeling framework allowing us to virtually perform stimulation, including DBS, and forecast the outcome from a dynamic systems perspective prior to invasive surgery with DBS lead placement. For an accurate prediction of the effects of DBS, we implement a detailed spiking model of the basal ganglia, which we combine with TVB via our previously developed co-simulation environment. This

Introduction
Deep brain stimulation (DBS) is a neuromodulation technique that has shown beneficial effects for patients suffering from many different neurological disorders (Horn, 2019;Horn & Fox, 2020). DBS is an essential element in the therapeutic regime for movement disorders like Parkinson's disease (PD) (Deuschl et al., 2006;Vitek et al., 2020), dystonia (Kupsch et al., 2006) and essential tremor (Koller et al., 1997). It provides a treatment option for selected cases of medication-refractory epilepsy (Salanova et al., 2015) and obsessivecompulsive disorder (OCD) (Anderson & Ahmed, 2003;Franzini et al., 2010;Nuttin et al., 2008). For major depression (Mayberg et al., 2005), Tourette's syndrome (Ackermans et al., 2011), Huntington's disease (Gruber et al., 2014) and alcohol addiction (U. J. Müller et al., 2009), DBS has shown first treatment successes and is clinically applied on an experimental basis. Albeit the initial implantation surgery, DBS is a reversible neuromodulation technique, in contrast to a permanent effect after a surgical lesion (Horn & Fox, 2020). Despite the benefits of DBS for many diseases, underlying mechanisms are so far poorly understood. At various scales of the brain, attempts have been made to model the outcome of DBS, from single-neuron to whole-brain models (Humphries et al., 2018). However, a multiscale model to bridge these different scales in a single DBS model has yet to be developed.
The most extensive research for DBS has been performed in movement disorders, which share pathology of the interactions between basal ganglia (BG), thalamus and cortex (Plotkin & Goldberg, 2019). The BG are anatomically defined by the striatum and the pallidum, which can be further separated in globus pallidus internus (GPi) and externus (GPe). Functionally, the regions of the subthalamic nucleus (STN) and the substantia nigra, whose degeneration is a key factor in the pathogenesis of PD (Damier et al., 1999;Fearnley & Lees, 1991), are often included in the BG because of their strong interactions with it (Albin et al., 1989). In the following, the term BG refers to "basal ganglia and related nuclei" (Lanciego et al., 2012) according to the widely used understanding as a functional unit of the extrapyramidal system (Heimer, 1983).
The hypothesis that PD patients often suffer from a decreased activity level in the thalamic region causing the motor function to be impaired, resulting in bradykinesia or akinesia, has a long history (DeLong, 1990;Humphries et al., 2018;Jahanshahi et al., 2015). This decreased activity in the thalamus is probably caused by pathological hyperactivity of the globus pallidus as a failure symptom of the dopaminergic system (Dostrovsky et al., 2002), a theory first formulated by the classical rate model of the BG (Albin et al., 1989). The clinically most relevant stimulation targets for PD are the GPi and STN (Horn & Fox, 2020). It is a common approach to model the neurons of these key regions for DBS as a network, employing mathematical descriptions of neuronal behavior and interactions (Yu et al., 2020). Previous studies found alterations in several important pathways through the BG (direct, indirect and hyperdirect pathway, Figure 1) (Nambu et al., 2002) when simulating DBS. Models of a single neuron and its axon cables propose that STN-DBS causes GPi to fire at a regular frequency (Rubin et al., 2012). An extensive amount of previous literature exists modeling the connection from STN to GPe, the striatal microcircuit and different subparts of the cortico-basal-ganglia-thalamo-cortical loop as spiking networks (Yu et al., 2020). These subnetwork studies often include validation of the model details (parameters as well as results) with experimental data and suggest that STN-DBS changes the efferences of the BG to the thalamus by suppressing the burst firing of the GPi (Guo et al., 2008;Rubin & Terman, 2004). A recent application of subnetwork models shows that STN stimulation causes short-term depression of its own activity (Rosenbaum et al., 2014). This short-term depression theory was validated with empirical data from rodents and primates, where STN-DBS eliminated beta-band oscillations in the GPi in Parkinsonian primate brains (Moran et al., 2012). Another class of models including electrical fields and volume of tissue information proposes that STN stimulation causes heterogeneous effects for different neurons depending on their distance to the electrode Humphries & Gurney, 2012;. These volumetric models explain the heterogeneous effects on the firing rates of GPi neurons, observed experimentally in primates (Hahn et al., 2008;Hashimoto et al., 2003).
Most previously established models are based on a priori assumptions about dynamic changes in PD, i.e., assuming differences between PD and healthy subjects with regard to their functional connectivity strengths or their activity levels of striatal projection neurons (Humphries et al., 2018). Though these assumptions are well justified by empirical findings, they critically influence the model outcomes. In contrast, Hamker and colleagues proposed a data-driven spiking model of the BG (Baladron et al., 2019;Maith et al., 2020), that is a generic BG model has been fit to the individual subject data by optimizing its parameters such that features of the simulated activity correlated with the same features of the measurements. Recently, Maith et al. (2020) fitted this BG model for 20 PD patients after DBS implantation and 15 healthy controls with individual resting-state functional magnetic resonance imaging (fMRI) data. However, the whole cortex was so far modeled as a single spiking network node, lacking a whole-brain perspective. This computational model was implemented with the software Artificial Neural Network architect (ANNarchy), used for spike and rate coding of neuronal populations, as well as a combination of both in a single network (Vitay et al., 2015). Network models in ANNarchy are defined through equations written in "natural language". ANNarchy has been used to implement models of the BG pathways (Baladron et al., 2019;Gönner et al., 2020;Maith et al., 2020;Villagrasa et al., 2018), spatial attention and vision (Bergelt & Hamker, 2019;Jamalian et al., 2017;Larisch et al., 2020) and learning and memory (Gönner et al., 2017;J. Müller et al., 2018;Schmid et al., 2019).
The single-neuron and subnetwork models of the BG successfully suggest underlying mechanisms for the improvement of PD hypokinesia symptoms during DBS. However, they are not sufficient in describing the multitude of other effects that DBS potentially has on PD patients, e.g., rigidity, tremor and cognitive or behavioral changes . Therefore, extending local DBS effects of the cortex-BG-thalamus loop towards a largescale network should be the next goal in understanding DBS effects.
Previous studies explored mean-field approaches simulating the whole-brain perspective for virtual DBS (Saenger et al., 2017;van Hartevelt et al., 2014). Mean-field models make use of a physical simplification to enable simulating the average or so-called mean-field behavior of large populations. Simulating the whole brain with mean-field modeling has shown that DBS brought the patients' dynamical regime closer to a healthy one (Saenger et al., 2017;van Hartevelt et al., 2014). Specifically, van Hartevelt and colleagues (2014) showed that STN-DBS has widespread structural and functional effects after long-term use analyzing diffusion tensor imaging (DTI) data of a single PD patient. However, they needed to exclude the STN from their analysis as controls were missing MRI data of this region. Saenger et al. (2017) analyzed the fMRI data of 10 PD patients under both conditions DBS switched on and off and performed virtual DBS for different candidate regions, demonstrating their effects for the whole-brain dynamics. This first application of testing DBS effects on the whole-brain dynamics was performed on the group level, while an extension towards the individual level is required before clinical application.
With respect to whole-brain mean-field simulations, The Virtual Brain (TVB, thevirtualbrain.org) (Ritter et al., 2013;Sanz Leon et al., 2013) offers a neuroinformatics platform to simulate the effects of a virtual DBS. This in silico computation of the whole-brain effects of DBS requires only the MRI data of an individual patient as an input. Simulated brain activity with TVB reproduces empirical phenomena accurately over different modalities . Applying TVB in combination with simulated stimulation has shown the connection between different stimulation targets and functional resting-state networks based on normative surface-based human brain data . Spiegler and colleagues recently reproduced this finding of activating functional resting-state networks through focal stimulation for the mouse brain, where the results were in line with experimental data of optogenetic stimulation (Spiegler et al., 2020). Transcranial direct current stimulation simulated with TVB based on a normative connectome (Kunze et al., 2016) resembled empirical electroencephalography (EEG) findings. However, virtual DBS has not yet been investigated with TVB.
The different computational studies demonstrating the effects of PD and/or DBS on the BG network, from single-neuron studies to whole-brain networks, exemplify the multiscale nature of this research field (Humphries et al., 2018). So far, the whole-brain DBS modeling literature stands isolated from the extensive literature on spiking neural networks of the BG. Only region-wise properties have been compared. None of the dynamical insights from the spiking network literature have been incorporated into the mean-field modeling approaches of DBS. Therefore, in this study, we aim to demonstrate the framework for a multiscale cosimulation approach of virtual DBS. Our goal is to bridge the microscale of single neurons towards the recorded whole-brain signals in one simulation framework, which permits a holistic and comprehensive integration of existing findings. To run whole-brain mean-field simulations and additionally simulate any region's fine-scale neuronal dynamics, including spikes generated by inhibitory and excitatory neurons inside the region, we can use the recently developed TVB-multiscale co-simulation toolbox . TVBmultiscale extends TVB to perform multiscale co-simulations, whereby most of the nodes are simulated with TVB as mean-field models, and a few selected nodes are modeled as spiking networks by another suitable simulator.
In this study, we combine the detailed spiking network model by Maith et al. (2020) for the BG with mean-field simulations in TVB for all cortical regions. We interface the spiking network software ANNarchy with TVB to build the TVB-ANNarchy co-simulation framework . As an underlying connection between BG and cortical regions, we utilize a recently published normative connectivity atlas of these tracts (Petersen et al., 2019) and combine it with individually -that is subject-specific -fitted probabilities and weights from Maith et al. (2020) for the connections among the BG regions. As a first proof of concept, we simulate resting-state conditions for an exemplary control and PD patient network and perform virtual DBS targeting STN and GPi in the patient network. Next, we validate our model by comparing the effects of virtual DBS against results from literature. Our study addresses the following limitations of previous whole-brain DBS modeling studies: 1) We incorporate a previously validated spiking network model of the subnetwork of the BG within our whole-brain modeling. 2) We use an underlying (normative) connectome, which includes the STN, and combine it with individually fitted connectivity data to create an individual patient and control multiscale network. In this way, we offer a computational model that holds the potential to be easily translated towards the individual patient level and used as a 'sandbox' model before future DBS surgeries.

Spiking network model for the basal ganglia
The spiking network model and its dynamics (including parameters) were taken from a previous publication   (Figure 1). Eight neuronal populations were included, each with different properties. The cortex consisted of 600 excitatory neurons coupled with 150 inhibitory neurons (possessing a self-inhibitory connection). From the excitatory population of the cortex, spikes were transmitted to the STN, as well as the striatum. The striatum was modeled with two different inhibitory neuronal populations, the direct (dSN) and the indirect (iSN) striatal spiny projection neurons, each with a selfinhibitory connection. The GPe and GPi were each represented by inhibitory neurons with a self-inhibitory connection. The thalamus was also modeled as a spiking network node. . We implemented this model inside our TVB-ANNarchy framework with the underlying previously optimized connection weights and probabilities for the data of one control and one PD patient (taken from ). The direct pathway is shown here as the path from the excitatory cortical neurons over the direct striatal projection neurons to the GPi. Similarly, the indirect pathway goes from the cortex, over the indirect spinal projection neurons and the GPe towards the GPi. The third pathway through the BG is the cortex-STN-GPi pathway, which is also called the hyperdirect pathway. CxExcit: excitatory population of the cortex; CxInh: inhibitory population of the cortex; GPi: internal globus pallidus; GPe: external globus pallidus; STN: subthalamic nucleus; dSN: striatum, direct striatal spiny projection neurons; iSN: striatum, indirect striatal spiny projection neurons; Thal: thalamus; excit.: excitatory; inhib.: inhibitory; conn.: connection.

Figure 1: Structure of the basal ganglia spiking model. Previously published detailed basal ganglia (BG) model by
Each spiking network population was modeled by an enhanced version of the Izhikevich model (Izhikevich, 2004;Maith et al., 2020). For details of this previously published model, we refer to the Supplementary Material. Maith et al. (2020) optimized the connection probabilities and weights between the nodes to fit empirical fMRI blood-oxygen-leveldependent (BOLD) signal correlation data for each individual and each hemisphere separately. We used this optimized data from one of the controls and one of the patients (left hemisphere only). We selected these subjects as representatives of their groups because their regional firing rates were close to the respective mean values.

Multiscale co-simulation of TVB and ANNarchy
Every node in the TVB network represents a brain region and its dynamics are simulated with a mean-field approximation. The nodes are connected with weights and delays (computed from tract lengths given a transmission speed) that can be determined for individual subjects employing DTI. As a mean-field model for the cortical regions, we chose the reduced Wong-Wang-model (Deco et al., 2013), which is often used to replicate fMRI data (Aerts et al., 2018;Klein et al., 2021) (details in the Supplementary Material). For an overview of all variables used in this study, we refer to Supplementary Table 1. In the TVBmultiscale framework , co-simulation is based on the concept of TVB "proxy" nodes that are created inside the spiking network ( Figure 2). TVB "proxy" nodes are either stimulating devices, thereby mimicking TVB cortex node dynamics (i.e., mean-field spiking rates) and coupling to the spiking nodes, or output (e.g., recording) devices, thereby extracting spiking dynamics to be transmitted to TVB. Thus, TVB and the spiking network simulator communicate on the level of neuronal populations' mean-field activities. TVBmultiscale, which is continuously expanding, is freely available on github (github.com/thevirtual-brain/tvb-multiscale) and interfaces TVB with different spiking network simulators (currently Neural Simulation Technology (NEST) (Eppler et al., 2008) and ANNarchy).
Since the previous BG model implementation was fitted with empirical data using ANNarchy , we built an interface between ANNarchy and TVB. We developed python code to incorporate the ANNarchy simulator into TVB-multiscale (details in the Supplementary Material). We validated our implementation of the spiking network by Maith et al. (2020) against the authors' original ANNarchy code by performing short simulations without noise for the two selected subjects (Supplementary Table 3). Each TVB cortex mean-field node (prime notation for nodes modeled only as mean-fields nodes in TVB) couples to a node modeled in ANNarchy (notation without prime for the spiking regions) via the instantaneous spike rate variable , which drives a population of neurons (same size as for the excitatory cortex node of the spiking network by Maith et al. (2020)) generating correlated spike trains , where stands for the spike time of the neuron with index in the population of the "proxy" node and is the Kronecker delta. The generated spikes were weighted by and delayed by based on the TVB connectome and the optimized weights for each subject (see below). For details of the spike trains' generation, we refer to the Supplementary Material.

Figure 2: Implementation of the interface for the multiscale model. (A) TVB to
In the other direction, each node modeled in ANNarchy updates the state of the corresponding TVB mean-field node since it is still represented in the TVB model and couples to TVB nodes . The update utilizes an ANNarchy monitor that records spikes for each TVB time step. The recorded spikes are converted to an instantaneous population mean rate that overwrites an auxiliary TVB state variable, called the input rate . The latter drives a linear integration equation of another auxiliary TVB state variable, named integrated rate , which, in its turn, acts as a smoothing low pass filter to have time series similar to the TVB mean-field ones, where is the time constant of the integration and provides the number of spikes per second. Finally, the integrated rate overwrites the state variable of the TVB model. All the rest of the TVB mean-field nodes follow the equations of the mean-field model described in the Supplementary Material. We simulated two ANNarchy time steps (of ) for every TVB time step ( ).

Underlying connectivity
To connect the TVB nodes and the spiking network simulator, we needed to assign connectivity weights for the paths between the BG regions and the cortex. Acquiring accurate data for those tracts is challenging because structural MRI data inherits many limitations (Jones et al., 2013;Thomas et al., 2014). Recently, Petersen et al. (2019) published a state-of-the-art axonal pathway atlas for the human brain that combines previous results from histological and imaging data literature with expert knowledge of neuroanatomists and brain-imaging scientists who collaborated on defining those tracts applying a holographic visualization technique (Petersen et al., 2019) (details to be found in the Supplementary Material). We used this normative tract data by Petersen et al. (2019) to include a fine-grained parcellation for the BG and the thalamus (based on CIT-168 brain atlas (Pauli et al., 2018)) and detailed data of their pathways to and from the cortical regions because of its current use for clinical DBS planning (Noecker et al., 2021). Whereas  used the motoric parts of the BG regions only , we used the complete BG regions as a first approach. For the cortex, the automated anatomical labeling (AAL) atlas parcellation was applied (Rolls et al., 2015;Tzourio-Mazoyer et al., 2002). The files were transformed to the DBS Intrinsic Template Atlas (DISTAL) space (Ewert et al., 2018) and the number of streamlines between each region pair was counted. This procedure resulted in a whole-brain matrix for the pathways between the cortex and the BG structures.
Some additional preparation steps have been performed on the connectome. As a first demonstration and because Maith et al. (2020) also treated the hemispheres in isolation, we focused on the left hemisphere only. Thus, all regions belonging to the right hemisphere and the vermis have been deleted from the connectome together with all their connections. Additionally, the connections from the inhibitory neuronal populations in the BG (GPe, GPi and striatum) to any cortical regions have been set to zero as they are not biologically plausible from a functional perspective of movement regulation, leaving in this direction only the connections from the thalamus and the STN to the cortex. The resulting connectome included 57 regions (for a list of all included regions: Supplementary Table 4). Its weights were normalized by the sum of all the incoming connection weights of the corresponding region and by the 99 th percentile of all weights.
The previous work of Maith et al. (2020) optimized the connection probabilities and weights among BG regions per individual to best fit the empirical fMRI data. To personalize the normative connectome, we replaced the network among the BG and thalamus regions with the optimized weights computed by Maith et al. (2020) for the control subject and the PD patient, respectively ( Figure 3). This 'hybrid' connectome constituted normative connectome weights among the cortex regions and between cortex and BG but included individually fitted connection weights and probabilities for the spiking network of the BG. For the presentation and for determining the couplings between the two scales, we adjusted the normative weights to be in the same range of values as the optimized connection weights by scaling them with the ratio between the 95 th percentiles of both weight distributions. The global coupling of the TVB mean-field model was set for each subject to , i.e., we are canceling the above normalization for the weights among the TVB nodes (see next section for the exact procedure of selecting the value of 15). The conduction speed was set to , thus, determining the time delay of couplings among all nodes of the multiscale model. The tract lengths among all regions were approximated by the Euclidean distance between their center coordinates (Supplementary Figure 2).  Maith et al. (2020) for the left hemisphere of the analyzed control and patient (upper panels) overrode the within-BG connection weights inside the connectome (red square) based on (Petersen et al., 2019). Each entry in any of the three colored matrices represents the normalized number of streamlines that start in the region marked on the vertical axis and end in the region marked on the horizontal axis. The brain network in the lower panel shows all connections taken from the individually fitted weights and the BG regions in red and the other regions in blue with the connections taken from the normative connectome of the atlas by Petersen et al. (2019)

Fitting the co-simulation model to individual dynamics
We implemented the previous BG model by Maith et al. (2020) inside our TVB-ANNarchy framework (Figure 4). For the multiscale model ("TVB-cortex model"), we replaced the spiking node "cortex" with the whole brain connectomic model in TVB ( Figure 5). However, the input from the multitude of the TVB mean-field nodes leads to different driving dynamics of the spiking network than in Maith et al. (2020). We aimed for TVB driving dynamics that would exhibit (a) a mean firing rate across all TVB nodes of similar to motor cortex neurons at rest (Velliste et al., 2014), where the range of rate values across TVB regions is determined by the structural connectome; (b) low amplitude random fluctuations of the rate around the equilibrium point of the above mean rate, resembling the rate dynamics of the cortex node in Maith et al. (2020) (c) some correlation among the neurons' spiking, which in Maith et al. (2020) is due to the internal connectivity of the spiking cortex node populations. We set the operation point of the TVB mean-field network by progressively increasing global coupling until an equilibrium point was reached (for ) with a mean firing rate across the whole TVB brain approaching from below via a few "trial and error" simulations. After the equilibrium point was approximated, we increased the additive white noise to a standard deviation of 10 -4 allowing small fluctuations around the equilibrium point without changing the pattern of nodes with higher firing rates (Supplementary Figure 3 displays a characteristic TVB time series during co-simulation).

Figure 4: Schematic overview of our study design. The basal ganglia (BG) model of one control and one
Parkinson's disease patient were taken from the previous study . Next, we implemented the previous model inside our TVB-ANNarchy framework, not yet activating TVB, the socalled "spiking-cortex model". We confirmed that this implementation reaches similar firing rates as the one from the previous study (step "confirm"). As a second step, we replaced the spiking-cortex node with mean-field simulations using TVB. To stay in the range of the previously confirmed firing rates for the BG regions, we fine-tuned the connection weights from TVB to ANNarchy for the TVBcortex models of the control and the patient. So far, all of the described modeling steps were taken in resting-state conditions. As a third step, we stimulated the STN and the GPi as two frequently targeted regions virtually (virtual DBS) and analyzed the effects for the BG spiking network as well as for the cortical regions. We analyzed whether virtual DBS could bring the patient's brain dynamics closer to the healthy one. Whenever there is a brain next to the model (even when it is grayed out), the simulation took place inside the TVB-ANNarchy environment.
For the multiscale TVB-cortex model, the three connections from cortex to STN, dSN and iSN were substituted by the respective set of connections from each of the corresponding TVB nodes ( Figure 5). For scaling these connections, we created a "spiking-cortex model" by substituting the cortex node of the network from Maith et al. (2020) with an ANNarchy spike generator identical to the one used as TVB "proxy" nodes (Supplementary Material). The spiking-cortex model acted as the "bridge" between the noisy TVB cortex driving the multiscale model and the Izhikevich population spiking cortex of Maith et al. (2020). With this model, we performed resting-state simulations for both subjects. Then, we tuned -again via a few "trial and error" co-simulations -three interface factors , , scaling all TVB connections to STN, dSN and iSN, respectively, to approximate the mean population rates of the spiking-cortex model ( Figure 5, Supplementary Table 1). These interface factors multiply the TVB connectome weights resulting in the interface weights ( Figure 2A). This step was also taken to ensure that dSN and iSN can receive different input while the striatum connectivity is equally strong for them. This way, the sum across all TVB nodes is the resulting total weight of the cortex input to the BG spiking populations, which then has an effect close to the one of the optimized weights in Maith et al. (2020). The final mean firing rates were approximated within for all spiking populations except for the thalamus of the control network, which was within (Supplementary Table 5, Supplementary Figure 4).  For the results of the TVB-cortex simulations, we simulated each condition 10 times, randomly selecting initial conditions for the TVB state from a normal distribution with mean equal to the initial conditions used originally for fitting the resting-state simulations and standard deviation 0.1 (Supplementary Material). The simulation length for all of our simulations was . After each simulation, we computed the mean firing rate over the last .

Implementation of the DBS stimulus
Besides the resting-state co-simulations, we applied a stimulus to our multiscale model (starting at and lasting till the end) inside either GPi or STN as possible target regions ( Figure 5). We simulated the propagation of these stimuli and the whole-brain response to them to provide a first proof of concept of the possibilities of this kind of multiscale modeling. We tested the virtual DBS stimuli within the spiking-cortex network and the TVB-cortex simulation model.
To the GPi, we applied a constant inhibitory current stimulus of an amplitude of aiming at reducing its firing rate and therefore disinhibiting the thalamus. For the other DBS simulations, we applied two realistic stimuli to STN, a monophasic and biphasic pulse-like current because the former is the most commonly implemented stimulus in previous DBS simulation studies (Yu et al., 2020) and the latter is used in clinical practice (Krauss et al., 2021) (Figure 6). The monophasic stimulus is adapted from (Michmizos & Nikita, 2011) and the biphasic stimulus is similar to the one used in (Liu et al., 2020) (details in the Supplementary Material).

Effects of the stimuli on cortical regions
To investigate the effects of the different stimuli on the cortex, we compared our restingstate TVB-cortex simulations with the simulations including the stimuli for the patient network. We also investigated cortical differences between the resting-state condition of the control and the patient. For these comparisons, we calculated the region-wise difference of cortical firing rates between simulations. Firing rates were averaged over the last of a simulation. In the resting-state case of comparing the patient and the control, we subtracted the average firing rates of the control's resting-state simulation from the ones obtained with the patient network. For evaluating the cortical effects of the different stimuli, we subtracted the resting-state average firing rates from the ones of the stimulus-induced time series. In addition, the resulting regional differences were normalized by the mean difference over the obtained regional differences, for each subtraction separately.

Results
For the whole-brain co-simulation model, we visualized the raster plots of the spikingnetwork regions for the four conditions, resting-state, GPi-DBS, STN-DBS with a monophasic and STN-DBS with a biphasic stimulus (Figure 7). Comparing the resting-state firing rates, the largest difference between the patient and the control can be found in the thalamus (Figure 8 and Supplementary Table 4). The stimuli caused the biggest changes in firing rate in the stimulated regions themselves (STN or GPi, respectively) and also in the thalamus (Figure 8 and Supplementary Table 4). The GPi-DBS simulation induced disinhibition of the thalamus from the GPi, allowing the thalamus to fire more than in the resting-state condition. Both STN-DBS simulations, however, also showed increased thalamic activity compared to the resting state but together with an increased firing in the GPi (Figure 8). Compared with the resting-state firing of the control, the patient seems to come closer to the rates of the control in multiple regions of the BG during all DBS scenarios (Figure 8). The common mechanism over all three stimulation protocols was the increase in thalamic activity. Thus, the thalamus firing rate seems to "normalize" toward the healthy regime during virtual DBS.   For the spiking-cortex model, the effects of the stimuli could not be traced further towards the cortical regions since this model is isolated and lacks embedding in the whole-brain network. Comparing the resting-state activities of the cortical regions between patient and control showed an increased average firing rate in the frontal regions and a decreased firing rate in the postcentral gyrus for the patient ( Figure 9A). Regarding the cortical effects of stimulation for the TVB-cortex model, we plotted the differences measured by the average firing rate between the resting-state and each stimulus simulation per cortical region on the template brain in Figure 9B-D. In all three virtual DBS simulations, the largest induced changes among the cortical regions were found in the frontal lobe and additionally in the postcentral gyrus for the monophasic STN stimulus. Altered levels of firing rates induced in the GPi and STN by the stimuli appear to be conveyed towards these cortical regions altering their activity with respect to the resting state. Concerning the specific regions, the middle frontal gyrus and insula were for all three stimuli among the top five regions regarding the most increased firing rates induced by the stimulus. Interestingly, only the STN monophasic stimulus created a slight reduction of firing rate in the supplementary motor area. Figure 9: Effects of the stimuli on cortical regions. We plotted the differences in averaged firing rate over the last 1000ms of the simulation time on the template brain, subtracting the average rate obtained from one condition from the other. In addition, the resulting regional differences were normalized by the mean difference over the obtained regional differences for each subtraction separately.

Discussion
In this study, we introduced a multiscale modeling strategy for the brain network, which allows to model the spiking network dynamics of the BG subnetwork in detail while simultaneously offering a whole-brain perspective of the evolving dynamics. We showed a first proof of concept that this new resulting TVB-multiscale model generates biologically plausible activity in resting state and after virtual DBS. This model has the potential to forecast DBS effects for different locations and different configurations on an individual patient level.
Our presented results show that the DBS stimulus introduced on our patient network causes disinhibition of the thalamus, leading to an increased firing rate during stimulation compared to resting state. Even though our results are in line with the hypotheses formulated by the classical rate model (Albin et al., 1989), conflicting evidence from clinical studies suggests a broader perspective as reduced thalamic activity alone neither explains all symptoms of PD nor all existing therapeutic effects (Eisinger et al., 2019;Marsden & Obeso, 1994;Rodriguez-Oroz et al., 2009). As for the direct effect of the stimulus on the target region, the recent theory of short-term depression states that STN-DBS blocks the transfer of lowfrequency oscillations downstream, e.g., towards GPe and GPi, and brings the thalamic activity back to healthy functioning (Humphries et al., 2018). Other theories exist about the effects of the DBS stimulus being of excitatory, inhibitory or disruptive nature on its neighboring areas and a consensus has yet to be reached in this research field (Chiken & Nambu, 2016). Still, our results show that the thalamic activity was brought back to healthy functioning by DBS, which is in line with the general mechanism of DBS (Humphries et al., 2018).
The increased firing rates of the thalamus during our STN-DBS simulations are not caused by decreased GPi activity, which cannot be explained by the classical direct/indirect pathway model of BG. Empirical evidence supports the observed increased firing rates of GPi during STN stimulation (Reese et al., 2011), which were assumed to overwrite pathological activity patterns. One recent computational modeling study with optogenetic data of rodents has shown that increased GPi activity, when synchronized, is able to drive excitatory thalamic responses despite the inhibitory nature of the connection (Liu et al., unpublished results). The proposed underlying mechanism is that bursts of inhibition from GPi to thalamus can cause hyperpolarization and then post-inhibitory rebound firings of thalamus neurons. Postinhibition spikes or bursts are characteristic behavior of the Izhikevich neuronal model used in this study (Izhikevich, 2004). Taking this unclear mechanism of pacing into account, our model provides computational evidence supporting a more connectomic effect leading to thalamic activation.
There have not been previous studies of multiscale co-simulation of DBS. PD is a multiscale disease (Kerr et al., 2013) with pathological mechanisms at many different scales, from deterioration observed in single neurons up to large-scale brain dynamics. Thus, in the attempt of modeling the broad perspective of potential treatment effects, one should also no longer focus on a single scale. One previous study embedded a spiking network for BG regions inside a neural field model for the cortex (Kerr et al., 2013). However, this previous modeling strategy did not subdivide the cortex mean-field model further into separate regions nor did the authors simulate DBS.
Compared to spiking models that encompass the BG regions only, our presented model can show whole-brain effects of stimulation going beyond the motor cortex. The presented results show an increase in overall activity in cortical regions for all of the three applied stimuli. This result is in line with the theory that PD patients have lower thalamic activity and, thus, a weaker driving activity from the thalamus towards the cortex. Subsequently, the cortex reacts with an increase of activity to the DBS-induced disinhibition of the thalamus. The frontal regions and the insula seem to be most impacted by all three different stimuli, measured by an increase in firing rate. The insula is linked strongly with non-motor symptoms in PD (Christopher et al., 2014) and a previous study reported a BOLD signal increase in the insula during STN-DBS (Kahan et al., 2012). The middle and inferior frontal gyrus also demonstrated one of the biggest shifts between DBS-OFF and DBS-ON condition measuring fMRI (Saenger et al., 2017). Interestingly, the monophasic stimulus applied on STN provoked a slight decrease of activity in the supplementary motor area in our results. This finding is in line with experimental results showing that DBS weakens excessive phaselocking interactions in the motor areas of PD patients (de Hemptinne et al., 2015). Supplementary motor areas, which are located at the transition between primary motor areas and prefrontal cortex, are involved in intentional movement initiation (Goldberg, 1985) and their impaired function is supposed to contribute to PD symptoms (Jacobs et al., 2009). Direct stimulation of supplementary motor areas with transcranial magnetic stimulation leads to improved freezing of gait symptoms in PD (Kim et al., 2018;Shirota et al., 2013), while dopaminergic medication can be related to improved supplementary motor area activation and improved motoric functions (Jenkins et al., 1992;Rascol et al., 1994). STN DBS in PD has been shown in fMRI (Stefurak et al., 2003) and positron emission tomography (Ceballos-Baumann et al., 1999) studies to activate motor as well as premotor areas, concordant with the simulated patterns in this work.
Our spiking network relies on a high level of biological realism with regards to spatiotemporal dynamics. Space refers to the fact that the spiking network receives input from different brain regions of TVB, which is closer to the reality regarding the multitude of different white matter connections between the cortex and the BG (Lenglet et al., 2012). More realistic time modeling implies the specific mean-field model dynamics that are chosen, as opposed to other studies, in which spiking networks are driven by Poisson spike trains, white noise or harmonic oscillations (Humphries et al., 2006;Park et al., 2011;Terman et al., 2002). The former approach of driving these spiking networks with noise seems to be an abstract view of the biologically underlying phenomena (Kerr et al., 2013).
There are still several open challenges in the field of DBS research that a virtual testing environment could potentially address. First, the exact placement of the electrodes seems crucial for the clinical outcome for patients. For PD and OCD, recent studies have shown that the connectivity profile of the brain area encompassing the inserted electrode predicts clinical outcome measures for patients (Baldermann et al., 2019;Horn et al., 2017Horn et al., , 2019Joutsa et al., 2018). This phenomenon was validated for dystonia (Corp et al., 2019;Okromelidze et al., 2020), essential tremor (Al-Fatly et al., 2019) and epilepsy (Middlebrooks et al., 2018). Testing the effects of different placement strategies before surgery could provide simulation-based advice for neurosurgeons. In this first co-simulation approach for DBS, we modeled stimuli targeting the GPi or STN area directly and completely. The clinical reality looks more complex (Krauss et al., 2021). The different sub-areas within the STN, for example, are involved in different pathways, i.e., the sensorimotor, associative and limbic loop. As most DBS systems provide several lead contacts to choose, the precise stimulus location is a common problem in clinical fine-tuning of DBS. With the upcoming of more detailed brain atlases, one could easily extend our used parcellation towards a finer grid inside the BG and model these subparts separately. Here, we presented the scaffold model that can be fine-tuned towards a more realistic model in a straight-forward manner. Second, so far, little individual information is considered for each patient and often the electrodes are placed based on normative data (Horn & Fox, 2020). Fitting an individual TVB model for patients provides a more personalized approach based on individual structural and functional imaging or electrophysiological data. TVB has previously been applied to help with predictions of clinical features for individual patients. Using individual positron emission tomography images, EEG slowing in patients with AD could be inferred from Abeta accumulation with the help of TVB (Stefanovski et al., 2019). Recently, a study has shown that the TVB feature of simulated mean local field potential frequency per brain region significantly improves the classification of individuals as AD patients, mild cognitive impairment patients or healthy controls using machine learning . For epilepsy patients, TVB has successfully been applied to optimize the determination of the resection and epileptic zone per individual before surgery .
A personalized virtual brain including structural data and dynamics based on MRI data is flexible in exploring other neuromodulation techniques with little extra effort. The hypothesis is that neuromodulation techniques can move the brain network dynamics between the diseased and healthy state ( Figure 10). With the current study, we have made a first attempt to "control" brain network dynamics by modeling stimulation in the brain of a PD patient. There is evidence that PD patients could also benefit from other neuromodulation techniques (Brittain & Cagnan, 2018). For example, a first study found that transcranial magnetic stimulation (TMS) of the supplementary motor area helps to improve the motoric symptoms of PD patients (Shirota et al., 2013). With our co-simulation framework, we can potentially analyze the impact of such a stimulation originating on the surface and follow the complete loop of cortico-basal-ganglia-thalamic-cortex connectivity. The flexibility of the presented virtual model could help with finding the best therapy for each individual patient. Figure 10: Schematic overview of "controlling" brain dynamics from one state to another using The Virtual Brain. Different interventions, e.g., deep brain stimulation (DBS), transcranial magnetic stimulation (TMS) or pharmacological interventions, can shift the brain dynamics from one state to another. These neuromodulation techniques hold the potential to alter the brain dynamics from a diseased brain towards a healthy target brain. Using The Virtual Brain, we aim to explore the different pathways leading to healthy functioning in a virtual environment for individual patients.
The applied data-driven model from Maith et al. (2020) does not make use of any prior assumptions regarding the pathological PD activity within the BG network, which stands in contrast to many previous models (Leblois et al., 2006;Lindahl & Hellgren Kotaleski, 2016). Fitting the outcomes of a model with empirical data from patients and controls offers an alternative approach to determining BG and whole-brain model dynamics. With this primarily data-driven approach, Maith et al. (2020) found many similarities of the obtained personal models of individual PD patients with physiological findings of PD, such as lower firing rates in the thalamus.
Our study inherits some limitations. In this proof-of-concept study, we modeled the TVB input that drives the spiking BG network with the reduced Wong-Wang mean-field model (Deco et al., 2013). In an improved version of this model, we could adjust the TVB meanfield dynamics to qualitatively correspond better with the original spiking network of Izhikevich neurons  by taking advantage of existing mean-field approximations of such networks (Nicola & Campbell, 2013;Visser & Van Gils, 2014). Such a choice would allow a more accurate analytical and computational determination of the large-scale brain dynamics (e.g., involved bifurcations) and inform the interface modeling between the two scales accordingly (e.g., in terms of scaling or more complex transformations). Further, alternatives to correlated spike trains' generators for converting the TVB mean-field nodes' rates into spike trains of TVB "proxy" nodes could be more effective in mimicking the Izhikevich spiking cortex node dynamics. All of the above options can be better explored by an upcoming computationally optimized version of the TVBmultiscale toolbox, implementing parallel co-simulation, allowing for a systematic exploration of the parameter space of the multiscale model to better fit individual neuroimaging data. So far, we fitted the virtual co-simulation brains to two individuals, which can be easily extended to larger cohorts with the only necessary data being DTI and either fMRI or electrophysiological data to fit the model dynamics accurately. The well-known characteristic of PD patients to demonstrate hyper-synchronization in the beta band ( ) in the sensorimotor network and the STN (Cruz et al., 2011;Whitmer et al., 2012) is reversed by DBS (Kühn et al., 2008;Wingeier et al., 2006). Our approach did not yet incorporate modeling the electrophysiological signatures of virtual DBS. However, TVB has often been used to monitor EEG-like activity from simulated time series and TVB-multiscale will soon also be equipped for this monitoring. Short-term plasticity probably plays an essential role in DBS effects (Milosevic et al., 2018), which has not yet been implemented in our model. Similarly, long-term plasticity effects due to DBS probably exist in structural and functional networks (van Hartevelt et al., 2014). With the spiking model allowing for an implementation of plasticity rules, we could explore its effects on the whole-brain dynamics with our model in future work. Moreover, our BG network misses the substantia nigra region as a crucial factor influencing PD dynamics and so far, we limited our analyses to a single (left) hemisphere.

Conclusions
In this study, we presented a co-simulation model for the BG as a spiking network together with TVB mean-field simulations for the whole brain. Our results show biologically plausible effects of virtual DBS performed in this multiscale modeling framework, bringing the patient's network dynamics of the BG closer to the healthy regime. The presented model offers a bridge between the different scales affected by DBS in the brain. It has the potential to be used as a 'sandbox' model for individual patients suffering from different neurological disorders prior to surgical interventions. Different strategies for DBS lead placements and configurations can be tested and evaluated. Future work needs to validate this model in larger patient cohorts and establish its link with clinical post-surgery improvement.