A computational model‐based analysis of basal ganglia pathway changes in Parkinson’s disease inferred from resting‐state fMRI

Previous computational model‐based approaches for understanding the dynamic changes related to Parkinson's disease made particular assumptions about Parkinson's disease‐related activity changes or specified dopamine‐dependent activation or learning rules. Inspired by recent model‐based analysis of resting‐state fMRI, we have taken a data‐driven approach. We fit the free parameters of a spiking neuro‐computational model to match correlations of blood oxygen level‐dependent signals between different basal ganglia nuclei and obtain subject‐specific neuro‐computational models of two subject groups: Parkinson patients and matched controls. When comparing mean firing rates at rest and connectivity strengths between the control and Parkinsonian model groups, several significant differences were found that are consistent with previous experimental observations. We discuss the implications of our approach and compare its results also with the popular “rate model” of the basal ganglia. Our study suggests that a model‐based analysis of imaging data from healthy and Parkinsonian subjects is a promising approach for the future to better understand Parkinson‐related changes in the basal ganglia and corresponding treatments.


| INTRODUCTION
Neuroimaging methods have contributed greatly to the understanding of the human brain. Functional magnetic resonance imaging (fMRI) measures neuronal activation indirectly through changes in blood flow and oxygen saturation. In resting-state fMRI (rs-fMRI; Biswal, Zerrin Yetkin, Haughton, & Hyde, 1995;He et al., 2009), the spontaneous fluctuations of the blood oxygen level-dependent (BOLD) signal are measured while the subjects remain in a constant resting state. Then correlations of BOLD signals between different brain regions can be determined, which are called functional connectivity in the field of rs-fMRI (Friston, 2011). Although there is no direct link between these correlations of BOLD signals and the underlying neuronal connectivity or functionality, several studies showed that the spontaneous fluctuations of the BOLD signals have a neurophysiological basis (Biswal et al., 1995;Brookes et al., 2011;He, Snyder, Zempel, Smyth, & Raichle, 2008). Computational models have greatly complemented the analysis of functional connectivity (for a review see Popovych, Manos, Hoffstaedter, & Eickhoff 2019). Detailed biologically inspired models of the brain can be adjusted to reproduce the observable neuroimaging data and then can be examined in detail on various scales. Thus, this model-based analysis allows to extract information, such as connectivity strengths and spiking activity, that cannot be achieved from the imaging data alone. Several computational studies have demonstrated the potential of this approach (Cabral et al., 2013;Deco & Jirsa, 2012;Deco et al., 2013;Schirner, McIntosh, Jirsa, Deco, & Ritter, 2018;Schmidt et al., 2018;Van Hartevelt et al., 2014).
We here use these techniques to reveal more insight into functional connectivity changes induced by Parkinson's disease (PD), which is a neurological disorder characterized by numerous motor, but also non-motor symptoms (Jankovic, 2008). It is directly or indirectly resulting from the loss of dopaminergic cells in the substantia nigra pars compacta (SNc), which leads to a reduced level of dopamine (DA) in its target structures such as in the striatum and pallidum, main nuclei of the basal ganglia (Bernheimer, Birkmayer, Hornykiewicz, Jellinger, & Seitelberger, 1973), which in turn likely affects processing in the cortex-basal ganglia loops (for reviews Galvan, Devergnas, & Wichmann 2015;Nambu, Tachibana, & Chiken, 2015)]). The findings about the changes in the basal ganglia due to PD have been integrated into conceptual models, such as the influential rate model of the basal ganglia (Albin, Young, & Penney, 1989;DeLong, 1990). As post hoc, conceptual models have less predictive power, neuro-computational models at different levels of detail have been developed to reveal further insight into differences between healthy and Parkinsonian states of processing within the basal ganglia (for reviews see Humphries, Obeso, & Dreyer, 2018;Maia & Frank, 2011;Rubin, 2017;Schroll & Hamker, 2016).
In those models, a PD state has been either introduced by rules that affect the neural activity of striatal projection neurons dependent on DA or DA-dependent learning rules or by directly assuming changes in functional connectivity strengths. While those assumptions are generally well motivated by findings, the particular validity of each assumption is difficult to prove. Thus, the predictions made by these neuro-computational models with respect to how PD changes processing in the basal ganglia critically depend on assumptions built into the model about the Parkinsonian state.
Here, we present a completely different approach to explore PD-related changes in basal ganglia processing which is more data-driven. We fit a spiking neuro-computational model of the basal ganglia to experimental functional connectivity data from PD patients and control subjects. The data from PD patients have been previously obtained by Horn et al. (2019). Instead of including specific assumptions on DA deficiency in the basal ganglia in our model, we only assume that PD leads to altered connectivity in the basal ganglia. The connectivity parameters of the different connections of our spiking neuro-computational model are adjusted individually to match each PD patient's and control subject's functional connectivity data comprising the sensory-motor parts of the cortex, subthalamic nucleus (STN), striatum, internal and external globus pallidus (GPi, GPe) and thalamus providing us an individual model for each PD patient and control subject. While parameter fitting took place at the level of correlations in BOLD signals, we bridge scales and obtain individual spiking resting-state models of the cortexbasal ganglia motor loop. We analyze the data-induced computational differences between the control and patient group and compare them with the rate model of the basal ganglia (DeLong, 1990) and recent data.

| Data acquisition
The rs-fMRI data used here are from the PD patients in deep brain stimulation OFF condition (DBS OFF) and control subjects studied in resting state by Horn et al. (2019). Particularly, we look at (for each subject) the functional connectivity of the basal ganglia, thus the correlations between the BOLD time series obtained in the following brain areas of the motor brain: cortex, STN, striatum, GPe, GPi and thalamus. As these motor brain areas are symmetrical in the left and right hemisphere, our data consist of two functional connectivity matrices per subject, one related to these six brain areas on the left side and the other to the same brain areas on the right side of the brain.
In the study carried out by Horn et al. (2019), PD patients (20 in total, four of them were women) were scanned (with a repetition time of 2,690 ms) first during DBS ON and after an interval of 5-15 min with the impulse generator switched off (DBS OFF). Patients were 63 ± 6.6 (M ± SD) years old when DBS surgery took place and they were scanned 30 ± 21 months after surgery with their usual medication ON (Levodopa). The rs-fMRI time series were detrended; noise related to motion, cerebrospinal fluid, and white-matter were removed; and a bandpass filter between 0.009 and 0.08 Hz was applied. For details about medication, patients' inclusion criteria, scanner, preprocessing of BOLD time series, voxel-wise parcellation and links to open sources please see Horn et al. (2019).
Controls (15 in total, four were female) were 59.5 ± 11.9 years old, and their rs-fMRI data were collected within Parkinson's Disease Progression Marker Initiative (PPMI) database (ppmi-info.org; Marek et al., 2011). These data were processed identically to the patients' data of the same study.

| Network description
We used a spiking neuro-computational model of the basal ganglia based on our previous work (Baladron, Nambu, & Hamker, 2019), but determined the connectivity parameters by an optimization procedure. The model is composed of several populations of spiking neurons, each representing a different nucleus of the basal ganglia, thalamus or cortex. The membrane potential of each cell is computed using the Izhikevich (2004) neuron model: where V is the membrane potential, U the recovery variable, g AMPA and g GABA the synaptic conductances, E AMPA and E GABA the corresponding reversal potentials, C the membrane capacity, I the baseline input current and µ the noise of the synaptic input current. Further, n 2 , n 1 , n 0 , a, b, c and d are parameters that depend on the neuron type of the corresponding population. The fixed parameters for each population are shown in Table 1. When the membrane potential V reaches 30 mV (40 mV for striatal neurons), the action potential is considered being triggered, the membrane potential V is reset to the value c and the recovery variable U is increased by d.
The noise µ of the synaptic input current fluctuates between different random fixed points introduced by changing the value SN of each neuron after every second of simulation. The values of SN are drawn from a normal distribution with a mean value (M SN ) of 0. The standard deviation SD SN of this distribution depends on the corresponding population (Table 1). The cortex is considered having higher levels of noise in order to drive the network, and therefore, instead of a mean M SN of 0 the mean value M SN is drawn from a uniform distribution ranging from −5 to 5. The striatum is modeled by two populations of 200 inhibitory spiny projection neurons. One of these populations represents the direct spiny neurons (dSN) which initiate the direct pathway through their projections to the GPi while the other population represents the indirect spiny neurons (iSN) which initiate the indirect pathway through their projections to the GPe (Figure 1). The parameters for these two neuron types were taken from the model of .  tuned and extended a medium spiny neuron model of Izhikevich (2007; for details see Humphries, Lepora, Wood, & Gurney, 2009) to replicate data of a complex multi-compartment model (Moyer, Wolf, & Finkel, 2007) which matches in vitro wholecell recordings from medium spiny neurons within the nucleus accumbens of rats (Wolf et al., 2005).
The GPi, GPe, STN and the thalamus are modeled each by an additional population of 200 spiking neurons. For the GPi, GPe and STN populations, the parameters were taken from Thibeault and Srinivasa (2013), who got their parameters by recreating the basal ganglia model for action selection (1) of Humphries, Stewart, and Gurney (2006) with Izhikevich spiking neuron models. The model of Humphries et al. (2006) matches in vivo single-cell recordings and tonic firing rates from the STN and the globus pallidus of rats. Both parameter sets were also used in Baladron et al. (2019). For the thalamic population, the parameters of the phasic bursting model proposed by Izhikevich (2004) were used. The values of all cell parameters are given in Table 1. Following our previous modeling approach Baladron et al. (2019), the GPe projects to the GPi and has reciprocal connections with the STN. The GPi combines excitatory input from the STN with the inhibitory input from both the dSN and the GPe and provides a constant inhibition to the thalamus. Further, the thalamus provides feedback to all striatal cells, thus closing a basal ganglia-thalamus loop (Hunnicutt et al., 2016;Parent & Parent, 2004;Sadikot, Parent, & Francois, 1992;Sidibé, Bevan, Bolam, & Smith, 1997;Smith et al., 2014). Additionally, local inhibitory connections were included in the GPi, GPe, dSN and iSN population.
The cortex is modeled by two populations: one excitatory and one inhibitory. The excitatory population is composed of 600 neurons while the inhibitory of 150, following the observed proportion between cortical projection cells and interneurons. For both populations, the parameters of the tonic spiking model proposed by Izhikevich (2004) were used. The populations form an excitatory-inhibitory loop, with connection probabilities and weights set by the optimization process. The excitatory cortical cells project to all striatal cells and to the STN.
Synaptic contacts of cells between nuclei were defined stochastically, given a connection probability determined by the optimization procedure. We use conductance-based synapses in our neuron models. Therefore, the synaptic current of a cell depends on the reversal potential (E SYN ) and the synaptic conductance (g SYN ) of the individual synapses. Excitatory connections are modeled by AMPA synapses and inhibitory connections by GABA synapses. A positive AMPA conductance drives the membrane potential to the value of E AMPA = 0 mV and the GABA conductance drives the membrane potential to the value of E GABA = −90 mV. The conductances of these synapses decay exponentially, with time constants τ AMPA = 10 ms and τ GABA = 10 ms (see Equation 1). The rise of the synaptic conductance after a presynaptic action potential is modeled as an instantaneous increase by the weight value of the connection.
The neuro-computational model has been implemented using the ANNarchy neural simulator (Vitay, Dinkelbach, & Hamker, 2015;version 4.6.9b). To solve the differential equations of the model, the Euler method is used with a time step of 0.1 ms.

| BOLD signal computation
The BOLD signal is computed using the Balloon model as described by Friston, Mechelli, Turner, and Price (2000); however, as shown in Table 2, some parameters, including the time constant of the BOLD signal, have been selected from Friston, Harrison, and Penny (2003). Previous studies motivated that the neural input (synaptic activity) is a better predictor of the BOLD signal than the output signal (Logothetis, Pauls, Augath, Trinath, & Oeltermann, 2001;Mathiesen, Caesar, Akgören, & Lauritzen, 1998;Mathiesen, Caesar, & Lauritzen, 2000). Therefore, and similar to recent studies (Schmidt et al., 2018), our BOLD signal computation relies on the synaptic activity ( Figure 2). The procedure used has been included into ANNarchy. A new type of monitor was added which can be attached to neuron populations in order to obtain a simulated BOLD signal. First, synaptic activity is computed at the single neuron level: With its connections, the model forms the direct (cortex-dSN-GPi), indirect (cortex-iSN-GPe-GPi and cortex-iSN-GPe-STN-GPi) and hyperdirect (cortex-STN-GPi) pathway. For each connection, the synaptic contacts of cells between the presynaptic and postsynaptic population were defined stochastically, given a connection probability. A connection probability of 1 means that every neuron of the postpopulation is connected with all the neurons of the pre-population. Each connection has its own weight value, which is set as the weight for all the synaptic contacts of the corresponding connection. The connection probability and weights of all connections are fitted as described in the main text. dSN-direct striatal spiny projection neurons, iSN-indirect striatal spiny projection neurons, STNsubthalamic nucleus, GPe-external globus pallidus, GPi-internal globus pallidus 2282 | MAITH eT Al.
Each neuron has its own synaptic activity s j . It decays exponentially, with the time constant τ syn = 1 ms. Whenever a presynaptic action potential reaches a neuron (no matter if excitatory or inhibitory), its synaptic activity s j is increased by 1 divided by the number of all afferent synaptic contacts n aff,j of this neuron. The theoretical maximum increase in s j in one time step is therefore 1 (if presynaptic action potentials would occur at all synaptic contacts simultaneously). The calculation of the synaptic activity of a given region S R is described by Equation 3. It is calculated as the sum of all mean synaptic activities of the neuron populations belonging to the region p ∈ R with population sizes N p . Only the striatum region (dSN and iSN) and the cortex region (excitatory and inhibitory) consist of more than one neuron population. Finally, before feeding the synaptic activity S R into the Balloon model to produce the BOLD signal of the region, a noise value R (randomly drawn every second from a uniform distribution between 0 and 0.05) is added to the synaptic activity S R in order to model the noise affecting the BOLD signal.

| Fitting procedure
For each subject, a functional connectivity for the left and for the right hemisphere of the brain was available. We fitted our model to the functional connectivity of each hemisphere of a subject independently resulting in 15 controls for the left side and 15 more for the right side and 20 patients for the left side and 20 more for the right side.
To fit the model to the functional connectivity of a subject (patient/control, left/right), we used the Bayesian Adaptive Direct Search (BADS) algorithm (Acerbi & Ji, 2017). This

F I G U R E 2
Simplified scheme of how the BOLD signal in a particular region is computed from the neural activity targeting this region. We take here the example of the striatum region. First, the rate activity in the cortex, dSN, iSN and thalamus targeting the cells in the striatum (top dSN, bottom iSN) is shown. For simplicity, we only plot 6 neurons of the dSN and iSN (12 small circles). How the synaptic activity s j changes at the single neuron level is described in the main text. For each population of the region (here dSN and iSN), the mean of the synaptic activities s j is calculated. Then, the means of the populations are added to obtain the synaptic activity S R of the region (here the striatum). Finally, the obtained synaptic activity in the striatum, with the addition of noise ( R ), is sent to the Balloon model which computes the BOLD signal for the striatum (see the plot with the BOLD signal over time). The plots are obtained from simulating the model for 20 s and recording spikes, synaptic activity and bold signal 10 times per second. dSN-direct striatal spiny projection neurons, iSN indirect striatal spiny projection neurons is a general model-fitting tool and a free MATLAB package. For the fitting of each functional connectivity, we carried out 20 optimization processes with BADS. Each optimization process consisted of multiple optimization steps.
In an optimization step, BADS first selected the values of the 19 weights and probabilities. Then, a 250 s fitting simulation with the model was carried out, in which for each second, BOLD signal values of the cortex, STN, striatum, GPi, GPe and thalamus were stored, resulting in six BOLD time series. Then, the Pearson correlation matrix between these six signals was computed, which is the functional connectivity of the model. Finally, the loss value was obtained as defined in Equations 4 and 5: where the loss is the Frobenius norm of the difference (D) between the functional connectivity of the model (FC M ) and the functional connectivity of the subject (FC s ).
In each optimization process, BADS runs multiple optimization steps, each with different values for the 38 connectivity parameters, until the loss value converges to a minimum. In each of the 20 optimization processes, different initial values were randomly drawn for the 38 parameters. The probabilities and weights were initialized from a random uniform distribution ranging from 0.1 to 0.3 and from 0.005 to 0.01, respectively. The parameter values were limited in the range between 0.1 and 0.5 for the probabilities and between 0.005 and 0.015 for the weights. Finally, for each subject, we selected the 38 connectivity parameters of the optimization process that resulted in the lowest loss value.

| Rest period simulations
We compared the obtained individual models between the different groups (PD patients and controls) with respect to the mean firing rates at rest. Therefore, additional simulations were performed with the models (see Table 3 for an overview).
To obtain the mean firing rates at rest, a 15-s resting period was simulated for each model (fitted to a specific subject). The same seed value was set for all simulations. Thus, the simulations differed only in the fitted connection parameters. The cortex input current in these simulations was set to I = 7. In addition, the random mean cortical firing rate fluctuations have been deactivated by setting the mean value of the normal distribution from which the noise values SN are drawn to a constant value of 0, as with the other populations. Thus, ensuring that the cortex population has an average firing rate around 10-15 Hz, which is common for neurons of the motor cortex at rest (Velliste et al., 2014).

| Sensitivity analysis
There are many predefined neuron parameters (Table 1) in our model. Although most of them are well motivated by previous studies (Humphries, Lepora, et al., 2009;Thibeault & Srinivasa, 2013), it is unlikely that these parameters are exactly like those of the individual human subjects. Variations in the predefined parameters may affect the values of the fitted connectivity parameters. Therefore, we conducted a sensitivity analysis in which we varied all the predefined non-zero parameters of Table 1 to estimate how sensitive our results are to variations. Although the parameters of the neuron types of the excitatory (cortex) and inhibitory (cortexI) cortical neurons are given separately in Table 1, we used the same parameter values for cortexI as for cortex, except for I and SD, which are always zero for cortexI. As it was not possible for us to investigate how the predefined parameters would affect the fitted connectivity parameters due to the long duration of the fitting procedure (over one month for the 70 fitted models), we only investigated the sensitivity of the loss of the already fitted models. For each fitted model, we varied each of the 59 parameters from −5% to +5% in 0.3125% increments of its value specified in Table 1 while leaving the other 58 parameters constant. The simulation procedure was the same as in the  Table 3) and differed only in parameter initialization. The connectivity parameters (weights and probabilities) were not randomly initialized instead the parameter values obtained by the previous optimization procedure for each fitted model were used. Further in each simulation, one predefined parameter deviated from its standard value as described above.

| Fitting results
Models fitted to the left and right hemispheres of each subject are treated as individual models for the control and patient group in the analyses, providing us with 30 control and 40 patient models. Examples of experimental and their corresponding fitted functional connectivities of two control and two patient subjects are shown in Figure 3. The mean loss of all the fitted functional connectivities is M = 1.10 with a standard deviation of SD = 0.20. Further, we calculated the pairwise correlation (Pearson correlation coefficient) between each fitted and the corresponding experimental functional connectivity. The mean correlation of all the fitted functional connectivities is M = 0.89 with a standard deviation of SD = 0.04. To verify if the fitted models can be clustered into different groups of models (control/patient), we tested if the control models fit significantly better to the experimental control data than to the experimental patient data and vice versa for the patient models. Thus, for each group (controls/patients) we calculated the losses between the functional connectivities of the models and the mean experimental functional connectivity (MEFC) of both groups, according to Equation 4. Then for each group of models, we compared the losses obtained from the control MEFC with the losses obtained from the patient MEFC by performing an one-tailed t test for dependent samples. The significance of the two tests was determined with a significance level of α = 0.05 and further Bonferroni correction. As shown in Table 4, the functional connectivities of the control models have a significantly lower loss to the control MEFC than to the patient MEFC and vice versa for the patient models, meaning that our fitting processes led to two different groups of models.

| Comparison between patient and control models
In order to compare our results with the popular rate model of PD (Albin et al., 1989;DeLong, 1990), we computed the connection strengths between nuclei as a multiplication of weight value and connection probability. We compared the connection strengths of the control group with the connection strengths of the patient group by a two-tailed t test for independent samples. The significance of the resulting 19 tests was determined using the false discovery rate (FDR) method (Benjamini & Hochberg, 1995) with a significance level of α = 0.05. The results are shown in Table 5. Several significant differences in the connection strengths were found. Almost all significant differences have large effects according to the Cohen's d greater than 0.8 everywhere. Most connections were stronger in the patient group. Looking at the connection strengths in the direct, indirect and hyperdirect pathway, all three pathways were stronger in the Parkinsonian models. In the direct pathway, the dSN-GPi connection strength was increased by 69.4% and in the hyperdirect pathway the cortex-STN connection strength was increased by 22.9%. Especially noticeable were the differences in the indirect pathway. The increase in the cortex-iSN connection strength was particularly high at 110.2% and had also by far the largest effect size at 1.75. Further, the GPe-STN connection strength was increased by 81.5% and the GPe-GPi connection strength by 33.0%. To obtain mean firing rates at rest, a 15-s rest period was simulated with each model. The last 10 s of these 15-s rest periods were used to obtain the mean firing rates. To calculate the mean firing rate of a region, the number of spikes, which occurred in the respective population, was divided by the time and population size. We compared the mean firing rates between the control and patient models by a two-tailed t test for independent samples for each region. The significance of the resulting 8 tests was determined using the FDR method (Benjamini & Hochberg, 1995) with a significance level of α = 0.05. The results are shown in Table 6. The relative changes of the mean firing rates and connection strengths from the control group to the patient group are summarized in Figure 4. For the mean firing rates, several significant differences with medium or large effect sizes were found. Mostly, the mean firing rates were lower in the patient group models. Particularly, high reductions in the mean firing rate occurred for GPe (−16.0%, Cohen's d = −1.04) and thalamus (−31.9%, Cohen's d = −1.05). Again, the iSN was especially noticeable. Here we observed by far the largest increase in the mean firing rate with a relative change of 109.8% and an effect size of 1.81.
We further investigated how homogeneous the models of the control group and the patient group were. Therefore, we calculated the coefficient of variation (CV) of the mean firing rates and of the connectivity strengths for each group. The CV value results from dividing the standard deviation by the mean of a sample. The CV values of the control and Note: M c /SD c are the means/standard deviations of the losses between the functional connectivities of the models and the mean experimental functional connectivity of the control group. M p /SD p are the means/standard deviations of the losses between the functional connectivities of the models and the mean experimental functional connectivity of the patient group. An one-tailed t-test for dependent samples was made to compare these two samples. A * indicates significance after Bonferroni correction for a significance level of α = 0.05. Note: M c /SD c are the means/standard deviations of the connection strengths of all models of the control group. M p /SD p are the means/standard deviations of the connection strengths of all models of the patient group. A two-tailed t-test for independent samples was made to compare these two samples. A * indicates significance after FDR correction for a significance level of α = 0.05.  Note: M c /SD c are the means/standard deviations of the mean firing rates at rest of all models of the control group. M p /SD p are the means/standard deviations of the mean firing rates at rest of all models of the patient group. A two-tailed t-test for independent samples was made to compare these two samples. A * indicates significance after FDR correction for a significance level of α = 0.05.

T A B L E 6
The mean resting firing rates

F I G U R E 4
The relative changes of the mean firing rates and connectivity strengths. The numbers below each region name represent the relative change of the mean firing rate from the control group to the patient group. The relative change in connectivity strength is represented for each connection by the line thickness and color. The line color shows the direction of the changes. Dark gray indicate positive changes (increase of connectivity strength), light gray negative changes (decrease of connectivity strength). The line thickness shows the amount of the relative change, the thicker the more the connectivity strength changed. A black border around connections or regions indicates a significant difference in the mean firing rate at rest or mean connectivity strength between the control and patient group T A B L E 7 The coefficients of variation of the mean firing rates or connectivity strengths of the components of the direct pathway for the control group (CV C ) and the Parkinsonian group (CV P )

T A B L E 8
The coefficients of variation of the mean firing rates or connectivity strengths of the components of the indirect pathway for the control group (CV C ) and the Parkinsonian group (CV P ) models of the patient group were more heterogeneous. To investigate the heterogeneity of the different pathways, we calculated the mean CV values for each pathway from all its components (shown in Tables 7, 8 and 9). All pathways were more heterogeneous in the patient group. The indirect pathway had a mean ΔCV of 0.05 and thus the smallest increase of heterogeneity. The direct pathway with a mean CV value of 0.47 was the most heterogeneous pathway in the patient group models. The cortex-iSN connection and the iSN rate were more homogeneous instead of more heterogeneous in the patient group; thus, the increased activity of iSN cells was very common in the Parkinsonian models.

| Sensitivity analysis
In the sensitivity analysis, we obtained the loss values calculated from the differences between the functional connectivities of controls/PD patients and the functional connectivities of the corresponding fitted models (see Equations 4 and 5) for different parameter variations in the fitted models. For each parameter, which was varied, we calculated the relative change in the loss depending on the change of the parameter value. This was done for each of the 70 models. 90% of all loss changes over all models and parameter deviations were smaller than ±16%. They ranged from −43% to +363%. For each parameter, we calculated the standard deviation of the corresponding loss values. The loss values of the two parameters with the highest standard deviation (SD cortex n1 = 0.64, SD thalamus n1 = 0.59) and the two parameters with the lowest standard deviations (SD STN I = 0.05, SD thalamus a = 0.05) are shown in Figure 6. For cortex n1 and thalamus n1, the loss increased on average with increasing parameter deviation, but in some models, there was almost no change or even a reduction of the loss. Positive and negative parameter deviations did not have exactly the same effect but if the loss changed on average, it increased for both positive and negative parameter deviations. Thus, to indicate how sensitive the loss is to each individual parameter, we fitted a linear regression model without an intercept to predict the relative loss changes of all models using the absolute values of the relative parameter deviations. The slope coefficient obtained for each parameter indicated by what percentage the average loss of the models changed due to a 1% parameter deviation.

T A B L E 9
The coefficients of variation of the mean firing rates or connectivity strengths of the components of the hyperdirect pathway for the control group (CV C ) and the Parkinsonian group (CV P )  The slope coefficients for all the parameters are shown in Figure 7. They ranged from 0.12 to 23.61 with a median of 0.19. The median indicated that most slope coefficients were quite low. Regarding the different neuron populations, parameter deviations in the cortex and thalamus population in particular had a strong influence on the loss. This is probably because the cortex drives the entire model and the thalamus has a very similar influence through its feedback connections to the striatum. In terms of the various parameters, n0, n1 and n2, which were specified by (Izhikevich, 2003), clearly had the strongest effect. These are the parameters of the quadratic membrane potential equation, which therefore have a direct influence on the activity of the neurons.

| DISCUSSION
We have used a primarily data-driven approach to generate changes related to PD in a subject-specific basal ganglia computational model, based on the assumption that Parkinsonian and control models differ in their functional connectivity strengths. By tuning the connection strength and synapse density of each connection in our model to match the functional connectivity obtained from the model with the functional connectivity of individual control subjects and PD patients, we obtained individual basal ganglia models of healthy controls and PD patients. This allowed us to compare these two groups of models in terms of individual connection strengths and mean firing rates at rest. Many computational studies already investigated PD with basal ganglia models (reviews, Humphries et al., 2018;Schroll & Hamker, 2016). However, these models were not fit to data of healthy and PD subjects (Holgado, Terry, & Bogacz, 2010;Kumar, Cardanobile, Rotter, & Aertsen, 2011;Lindahl & Kotaleski, 2016;Mandali, Rengaswamy, Chakravarthy, & Moustafa, 2015;McCarthy et al., 2011;Terman, Rubin, Yew, & Wilson, 2002). To compare healthy and PD states in basal ganglia models, the most common approach so far is to implement the modulatory effects of DA into the model, where DA is either an external variable (Guthrie, Myers, & Gluck, 2009;Humphries et al., 2006;Leblois, Boraud, Meissner, Bergman, & Hansel, 2006;Lindahl & Kotaleski, 2016) or implemented within the model by neurons representing the SNc activity (Frank, 2005;Schroll, Vitay, & Hamker, 2014). States of DA then differentially affect the activity of striatal spiny projection neurons or the degree of long-term potentiation or depression. Thus, the predictions made by these neuro-computational models with respect to how PD changes processing in the basal ganglia critically depend on assumptions about the effects of DA in the basal ganglia built into the models. F I G U R E 6 Relative loss changes depending on relative deviations of four example parameters, cortex n1, cortex n2, striatum n1 and striatum n2. Each observation is presented by a small circle. The data were collected over the 70 fitted models; thus, there are 70 observations for each parameter variation. The 33 parameter variations ranging from −5% to +5% in 0.3125% increments. The black line indicates the mean over the 70 models. The plots are quite symmetrical and thus show that the change of the loss depends mainly on the absolute value of the relative parameter deviation. ΔLoss-absolute loss change of the corresponding model, Loss0loss value of the corresponding model with standard parameters from Table 1 | 2289 MAITH eT Al.

| The "rate model" of PD and corresponding neurophysiological findings
Due to the popularity of the so-called rate model of the basal ganglia, particularly the version of DeLong (1990), we compared the individual predictions of this rate model and our results with respect to experimental findings compiled in several reviews (DeLong & Wichmann, 2009;Galvan et al., 2015;Galvan & Wichmann, 2008;Lee & Masmanidis, 2019;McGregor & Nelson, 2019;Obeso & Lanciego, 2011;Obeso et al., 2000;Wichmann & Delong, 2002).
The rate model suggests that DA activates the dSN cells of the direct pathway and inhibits the iSN cells of the indirect pathway. Due to the reduction of DA in PD, iSN cells become more active which leads to a stronger inhibition of GPe, whereupon the STN is disinhibited. The hyperactive STN then causes a stronger stimulation of the GPi. In addition, the DA loss reduces the activity of the dSN cells, whereby the GPi is additionally disinhibited. Convergingly, this leads to a hyperactive GPi, resulting in a strong inhibition of the thalamus, whereby the cortex receives less excitatory input from the thalamus. The inhibited thalamus and the less excited or less responsive cortex are proposed to result in bradykinesia.
Although we have not implemented modulation by DA or DA changes in our models, similarities between the results of our study and the predictions of the rate model are apparent. The most consistent change that occurred in the Parkinsonian models was the stronger cortex-iSN connection and higher iSN firing rate. Despite some variability in experimental observations, more active iSN cells seem to be a very consistent finding in experimental studies. Increased firing rates were measured in PD animal models (Kita & Kita, 2011;Mallet, Ballion, Le Moine, & Gonon, 2006;Sharott, Vinciati, Nakamura, & Magill, 2017), and further findings on synapse density (Day et al., 2006), increased metabolism in GPe (Crossman, Mitchell, & Sambrook, 1985;Mitchell, Cross, Sambrook, & Crossman, 1986;Schwartzman & Alexander, 1985) and increased GABA levels in GPe (Bianchi, Galeffi, Bolam, & Della Corte, 2003;Robertson, Graham, Sambrook, & Crossman, 1991) indicate increased activity of the iSN cells. Only a few studies did not reproduce an increased iSN firing rate in PD (Ketzef et al., 2017;Ryan, Bair-Marshall, & Nelson, 2018). Further, our model predicts a lower firing rate in GPe. Also, these changes were reported in many animal model studies (Boraud, Bezard, Guehl, Bioulac, & Gross, 1998;Filion & Tremblay, 1991;Miller & DeLong, 1987Pan & Walters, 1988;Schwartzman & Alexander, 1985;Soares et al., 2004) but a few reports argue against it, like Bezard, Boraud, Bioulac, and Gross (1999), which observed an increase in STN but no decrease in GPe firing rate.
Our Parkinsonian models showed lower firing rates in the thalamus, as supported by many previous findings (Kammermeier, Pittard, Hamada, & Wichmann, 2016;Molnar, Pilliar, Lozano, & Dostrovsky, 2005;Ni, Gao, Benabid, & Benazzouz, 2000;Schneider & Rothblat, 1996;Vitek, Ashe, & Kaneoke, 1994). However, the findings for the thalamus are quite inconsistent (Galvan et al., 2015) including findings of no rate change (Pessiglione et al., 2005) or even increased rates (Bosch-Bouju, Smither, Hyland, & Parr-Brownlie, 2014). Concurrently, the thalamus in our PD models had the most heterogeneous mean firing rates. In fact, some PD patients' models showed mean thalamus firing rates above 20 Hz, which is higher than in some controls' models. Therefore, our results also suggest that some PD patients may have an increased rather than a decreased thalamus firing rate. F I G U R E 7 Slope coefficients of the linear regression model predicting the relative loss changes using the absolute values of the relative parameter deviations. The data were collected over the 70 fitted models and 33 parameter variations ranging from −5% to +5% in 0.3125% increments for each parameter. The slope coefficient of each parameter is shown by the width of the vertical bars. The background colors indicate the corresponding population of the parameters A further consistency of our results with the predictions of the rate model is the decreased firing rates of the cortex, which were previously demonstrated in studies on monkeys (Pasquereau, DeLong, & Turner, 2016;Pasquereau & Turner, 2011). According to the rate model, the reduced activity of the cortex should result from the excessively inhibited thalamus. This cannot be the case in our model as we closed the loop via the thalamus. As we have drastically simplified the cortex in our model and have neglected cortical inputs like cortex-thalamus loops, it is unlikely that our model plausibly explains why the cortex has a lower activity. Furthermore, the difference in the cortical firing rate is relatively small in our models.
In some points, our findings contradict estimations of the rate model. First, the activity of the dSN in our PD models is higher than the rate model would predict. We do not observe a significant decrease in the PD group models but even the opposite tendency. However, there are some sparse findings in the literature that describe increased dSN firing rates in PD (Kita & Kita, 2011;Mallet et al., 2006;Ryan et al., 2018) albeit they are less frequent than findings of dSN decreased rates. Some studies show that the metabolism in the GPi is increased (Crossman et al., 1985;Schwartzman & Alexander, 1985). This is usually understood to be resulting from increased input from the hyperactive STN. However, the increased metabolism could also be a consequence of increased afferent synaptic activity caused by the dSN cells. In addition, our results show that cortex-dSN connectivity strength, dSN firing rate and dSN-GPi connectivity strength are among the most heterogeneous effects in our Parkinsonian models.
Second, contradicting predictions of the rate model our findings indicate a less active GPi in PD. Although strongly debated, GPi overactivity predicted by the rate model has been partially supported by some studies (Boraud et al., 1998;Filion & Tremblay, 1991;Miller & DeLong, 1987Soares et al., 2004;Wichmann et al., 1999). Findings on increased metabolism in the thalamus also suggest an increased inhibitory input from the GPi (Mitchell et al., 1989;Rolland et al., 2007). However, these findings could also speak for a stronger GPi-thalamus connection, as was the case in our Parkinsonian models. Nevertheless, the findings of previous studies do not indicate that the GPi is less active in PD.
In summary, our data-driven approach illustrates PD changes similar to predictions of the rate model, which are in coherence with many experimental findings. This suggests a certain plausibility of the rate changes in our models, which result exclusively from fitting functional connectivity of patients and healthy controls. However, there are also some inconsistencies in our data with the rate model predictions, especially the lower firing rates in GPi and the lack of increased firing rates in STN.
The data to which our models were fitted were collected during active dopaminergic medication of the PD patients (Horn et al., 2019). This medication may have affected basal ganglia firing rates, which should be considered when interpreting our results. Levy et al. (2001) showed that DA agonist apomorphine administration leads to a generally lower GPi rate and a lower STN rate during movement in PD patients. This could be a reason why we did not find an increased STN or GPi rate in our models. In addition, the DBS electrodes implanted in the patients despite being switched off probably caused artifacts in the BOLD signal, especially in the STN signal, which reduces the validity of our results about the STN.

| Firing rates to explain PD
So far, we have only compared findings, which suggest or show certain neuronal activity changes with our results and the rate model. Of course, the rate model does not only describe the rates themselves, but also tries to explain PD symptoms by changes in the rates. Meanwhile, it has been shown that the predictions the rate model makes about the effects of certain rate changes on motor function cannot be replicated very consistently (DeLong & Wichmann, 2009;Galvan & Wichmann, 2008;Obeso & Lanciego, 2011;Wichmann & Delong, 2002) and other features such as firing pattern with prolonged bursts and pauses are most likely more important as a simple output rate (Hutchison et al., 1994;Levy, Hutchison, Lozano, & Dostrovsky, 2000;Molnar et al., 2005;Vitek et al., 1999;Wichmann & Soares, 2006). Moreover, invasive recordings in humans via DBS electrodes over the last two decades have shown that abnormal oscillatory activity in the cortico-BG network is related to motor symptoms such as increased beta band activity in PD bradykinesia and rigidity (Kühn & Volkmann, 2017).
We, therefore, do not want to draw any direct conclusions about the symptoms in PD with our findings on the rate changes. Rather, our study should be seen as a proof of concept, in that we obtained plausible rate changes, in accordance with many experimental findings, solely from the functional connectivity data of controls and patients.

| Limitations
Our approach is complementary to previous neuro-computational accounts to study the effects of PD in the basal ganglia. Most importantly, it directly relies on data from both, control subjects and PD patients, it requires no specific assumptions of how DA deficiency associated with PD changes processing in the model and it leads to subjectspecific neuro-computational models at the level of spiking activity. However, we shall also mention the limitations of the present design. Although we are able to bridge scales from BOLD data to neural spiking dynamics, one has to be aware that the free model parameters are tuned to match BOLD correlation data. Thus, no direct access to firing rates or local field potentials exists that could further constrain the model. Given BOLD data, the model is already fairly complex, but certainly far away from the true biological detail, as several neuron types and their connectivity are not included. An important prerequisite for the validity of our results is that the differences in functional connectivity between patients and controls were caused by PD. This is not completely fulfilled in our data. Even though the subjects were age-matched, the data were collected in two different laboratories. The patients' data were scanned at 3 T, and the controls' data were scanned at 1.5 T. In addition, the patients received a DBS electrode which may cause artifacts in the BOLD signal even when switched off. Besides the free parameters obtained by the optimization process, our model has a large number of fixed parameters that were not optimized by the optimization process. Our sensitivity analysis shows that some of these parameters have a strong influence on the losses of the models. Thus, changes in these parameters could also cause changes in our results. However, it is important to remember that these parameters should not be arbitrarily different from ours. Most of our parameters were chosen based on previous studies (Humphries, Lepora, et al., 2009;Izhikevich, 2004;Thibeault & Srinivasa, 2013) to mimic the firing properties of the respective neurons based on physiological data. The optimal case would be to have access to the specific parameters for each individual patient and control subject, but this would require single-cell recordings for all the different nuclei for each patient and control subject, which would not be feasible. Thus, we took a compromise and combined general assumptions (e.g., the neuron models/ firing properties) with easily accessible individual data (like MRI) to obtain individual models. Despite these limitations, we found a good match to existing data as outlined above, although present data are also quite variable.
However, not each prediction of the model should be interpreted literally. While the increase in iSN activity is very consistent among the subject-specific models of PD and also observed in recent studies (Kita & Kita, 2011;Mallet et al., 2006;Sharott et al., 2017), the massive increase in iSN firing rate is likely an overestimation. Thus, despite the promising observations, future studies should explore how to add additional constraints to the model or add more data to the optimization procedure. Our approach used in this study could be further extended with more complex models, albeit there is a limitation in the number of parameters being fit. Models that would replicate the neuronal mechanisms relevant for PD symptoms can then also be used to investigate treatment options for PD in more detail. As the method used in this study provides models of both, healthy and Parkinsonian states, it offers a promising environment to systematically search for treatments that bring the Parkinsonian models closer to the healthy ones. In addition, the method provides individual models of patients, a potential that has not been exploited in our present study. With individual models, tailored PD treatments could be investigated, comparable to the "virtual epileptic patient" (Jirsa et al., 2017). This might help to plan and improve treatment efficacy for the heterogeneous PD phenotypes.

| CONCLUSION
The approach to fit functional connectivity data with a spiking neuro-computational basal ganglia model and thus inferring firing rates and connectivity strengths has proven to be quite successful. We have indeed found many similarities with the rate model of the basal ganglia and corresponding recent physiological findings. Our results show quite meaningful changes in firing rates in the basal ganglia, associated with PD. These changes emerged only due to differences in the functional connectivity of PD patients and controls and were thus obtained without implementing effects of DA and DA changes associated with PD in the models. As firing rates alone have become somewhat outdated in terms of explaining PD symptoms, in the future, our approach could be further developed to include effects of synchrony, oscillations and bursts in PD, by using additional data in the fitting procedure. In addition, our study is a confident first step toward obtaining subject-specific models from non-invasive imaging data of PD patients and healthy controls that could be useful for tailored treatments of PD.