Transient Response of Basal Ganglia Network in Healthy and Low-Dopamine State

Abstract The basal ganglia (BG) are crucial for a variety of motor and cognitive functions. Changes induced by persistent low-dopamine (e.g., in Parkinson’s disease; PD) result in aberrant changes in steady-state population activity (β band oscillations) and the transient response of the BG. Typically, a brief cortical stimulation results in a triphasic response in the substantia nigra pars reticulata (SNr; an output of the BG). The properties of the triphasic responses are shaped by dopamine levels. While mechanisms underlying aberrant steady state activity are well studied, it is still unclear which BG interactions are crucial for the aberrant transient responses in the BG. Moreover, it is also unclear whether mechanisms underlying the aberrant changes in steady-state activity and transient response are the same. Here, we used numerical simulations of a network model of BG to identify the key factors that determine the shape of the transient responses. We show that an aberrant transient response of the SNr in the low-dopamine state involves changes in the direct pathway and the recurrent interactions within the globus pallidus externa (GPe) and between GPe and subthalamic nucleus (STN). However, the connections from D2-type spiny projection neurons (D2-SPN) to GPe are most crucial in shaping the transient response and by restoring them to their healthy level, we could restore the shape of transient response even in low-dopamine state. Finally, we show that the changes in BG that result in aberrant transient response are also sufficient to generate pathologic oscillatory activity in the steady state.


Introduction
Parkinson's disease (PD) is a debilitating neurodegenerative brain disease with multiple cognitive and motor symptoms. Etiologically the disease is attributed to the progressive loss of dopaminergic neurons in the substantia nigra pars compacta. Dopamine affects neuronal excitability, synaptic strength and synaptic plasticity. Consistent with this, data from human patients and animal models show that persistent dopamine deficit results in a number of changes in the neuronal activity especially in the basal ganglia (BG). At the level of neuronal activity, in PD, synchronized b band oscillations (15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30) in the globus pallidus externa (GPe) and subthalamic nucleus (STN; Raz et al., 2000;Brown et al., 2001;Mallet et al., 2008;Tinkhauser et al., 2017) emerge along with an increase in spike bursts (Tachibana et al., 2011;Nambu et al., 2015). Recent experimental studies also pointed out the role of GPe subpopulation (arkypallidal: GPe-TA, and prototypical: GPe-TI) in orchestrating the oscillatory activity in the BG subnuclei (de la Crompe et al., 2020). In the striatum, firing rate of D2-type dopamine receptors expressing spiny projection neurons (D2-SPN) is increased whereas firing rate of D1-SPNs is reduced (Mallet et al., 2006;Sharott et al., 2017) in the PD conditions. Moreover, while cortical inputs to D2-SPN are enhanced, inputs to D1-SPN are weakened (Parker et al., 2016;Ketzef et al., 2017;Filipovi c et al., 2019). The aforementioned changes in the activity and structure of the BG are persistent and indicate a change in "operating point" of the BG. But these observations do not provide mechanistic links between behavior deficits of PD and BG activity.
During action-selection or decision-making tasks the BG receives transient inputs (Gage et al., 2010) from different cortical regions. It is therefore important to understand how the response of the BG network to a transient cortical input is altered during PD condition. In a healthy state, transient cortical stimulation elicits predominantly a triphasic response (composed of early excitation, inhibition, and late excitation) in most neurons of the BG output nuclei, i.e., globus pallidus interna (GPi) or substantia nigra pars reticulata (SNr; Chiken and Nambu, 2013;Sano et al., 2013;Ozaki et al., 2017). The triphasic response is consistent with the predictions of a simple feedforward model of the BG involving the so-called direct, indirect and hyperdirect pathway (Albin et al., 1989;Jaeger and Kita, 2011). However, a small fraction of neurons in SNr (Sano and Nambu, 2019) or GPi (Iwamuro et al., 2017) can also show biphasic or monophasic responses. In lowdopamine conditions, the fraction of neurons showing triphasic, biphasic and monophasic responses is changed resulting in an altered population response.
To identify what determines the shape of BG transient responses we used a BG network developed by Lindahl and Kotaleski (2016). We found that, consistent with experimental data (Sano and Nambu, 2019) and the feedforward model of the BG (Albin et al., 1989), in healthy state, the SNr showed triphasic responses for brief cortical inputs at the population level. In the low-dopamine state, with the default settings, the SNr transient response was biphasic. However, by changing the strength of synapses along the direct (D1-SPN !SNr) and indirect pathways (D2-SPN !GPe-TI, and GPe-TI !STN) it was possible to observe the triphasic responses even in low-dopamine state. Interestingly, we found that changes in the transient response properties in PD state involve not only changes in the feed-forward connections (e.g., D1-SPN !SNr) but also recurrent interactions within BG subnuclei, e.g., the recurrent connections within the GPe (GPe-TA $GPe-TI) and between GPe and STN (GPe $STN). Next, we show that by restoring the connection from D2-SPN to GPe (D2-SPN !GPe-TI) to a normal value, even in low-dopamine state we can recover a transient response similar to that observed in healthy/normal state. Thus, the D2-SPN !GPe-TI emerged as the most important descriptor of the aberrant transient response. Interestingly, the same connections can also unleash b band oscillations (Kumar et al., 2011;Mirzaei et al., 2017). That is, the same changes underlie both the emergence of pathologic b band oscillations, and pathologic transient response.

Neuron model
In order to achieve a good trade-off between simulation efficacy and ability to capture the neuronal dynamics, we used two types of neuron models in our BG network. Striatal D1 and D2 type dopamine receptor expressing spiny neurons (D1-SPN and D2-SPN), fast-spiking interneurons (FSIs) and STN neurons were realized using the standard leaky-integrate-fire neuron (LIF) model with conductance-based synapses. The subthreshold dynamics of the membrane potential V x (t) was governed by the Equation 1: where x [ {D1-SPN, D2-SPN, FSI, and STN}. In Equation  1, C x , G x , V rest represent membrane capacitance, leak conductance and resting potentials, respectively. When V x reaches the threshold potential V x th , a spike is elicited and V x is reset to V x rest for refractory duration t ref = 2 ms. I syn (t) models the total synaptic input current received by the neuron (see Fig. 1 for the various sources of inputs to these neurons).
GPe-TA, GPe-TI and SNr neurons were modelled as a LIF neuron with exponential adaptation (AdEx) to capture rebound firing on release from hyperpolarization and spike trigerred adaptation as well (Nakanishi et al., 1987;Cooper and Stanford, 2000;Bugaysen et al., 2010). The subthreshold dynamics of these neurons were defined as: where x [ {GPe-TA, GPe-TI, SNr}. In Equation 2, V x T represents the spike-threshold, D T represents the slope factor, t w is the time constant of the adaptation variable w, and a controls the adaptation term. Given Equation 2, when V x reaches the spike-cutoff potential then a spike is generated and V x , as well as w x are reset at values V rest , w x 1 b, respectively, where b denotes the spike-triggered adaptation. I syn (t) models the total synaptic input current received by the neuron (see Fig. 1 for the various sources of inputs to these neurons).
The parameters for GPe-TA, GPe-TI and SNr neurons are given in the Tables 6,7,9,respectively. Whether the response is shaped by the neuron complexity or network interactions has been highly debated without any clear conclusion (Prinz et al., 2004;Marder and Taylor, 2011;Sahasranamam et al., 2016). Here, we have chosen to use simplified models, so that we can exclusively focus on network interactions. Moreover, it should be noted that while LIF neuron model may appear simple, we can vary its inputs and parameters to fit many diverse input-output firing rate relationships (indeed, that is what we have done here).

Synapse model
Neurons were connected using static conductancebased synapses. Each incoming spike elicited an a function shaped conductance transient, after a fixed delay since following the spike in the presynaptic neurons. The time course of the conductance transient was given as the following: where syn [ {exc, inh} and x 2 D1ÀSPN; D2ÀSPN; FSI; GPeÀTA; GPeÀTI; STN; SNr. In Equation 3, J x syn is the peak of the conductance transient and t x syn is synaptic time constant. Each incoming synaptic current induces a current transient as given by the following: where V x rev is the reversal potential of the synapse for a neuron in population x 2 fD1ÀSPN; D2ÀSPN; FSI; GPeÀTA; GPeÀTI; STN; SNrg. All synaptic parameters are specified in Table 2.

BG network
The BG comprises of striatum, STN, GPe, SNr, and GPi in primates or entopeduncular nucleus (EPN) in rodents ( Fig. 1). Although GPi and SNr are the output nuclei of the BG, but in this work, we only focus on the SNr activity. To model BG, we adapted a previously published model by Lindahl and Kotaleski (2016). However, unlike that model (Lindahl and Kotaleski, 2016), here, we reduced the time complexity of our proposed network by scaling down the size of striatum (D1-SPN, D2-SPN, FSI). Also a few synaptic and neural parameters were adjusted to achieve the network performance in healthy and PD conditions. The main differences between these two models are detailed in the later part of Materials and Methods.
Our reduced model of the BG consisted of 6539 neurons. Number of neurons in each subpopulation, number of connections and synaptic connectivity parameters are provided in Table 1.

Dopamine induced changes in neuron and synapse parameters
To model the effect of dopamine we followed the approach taken by Lindahl and Kotaleski (2016). Dopaminergic effects on SPNs, FSIs, STN, GPe and SNr neurons and their synaptic connections were modelled by modulating parameters such as the resting state potentials (E L ), spike threshold (V th ), and synaptic strengths. The dopamine modulation was modeled using a parameter a dop ranging between 0 (PD condition) and 1 (high dopamine). The normal state was mapped to a dop = 0.8 (a normal ). The effect of the a dop on the neuron and synaptic properties are presented in the subsequent sections.

Dopamine effects on neuron properties
In D1-SPNs, D1 type dopamine receptor activation not only shows a hyperpolarizing effect by increasing potassium inward rectifier (KIR) current, but also induces depolarizing effects on the resting membrane potential (Gruber et al., 2003). We modelled these two contributions by changing the spike threshold and resting membrane potential: th and E D1ÀSPN L were reduced. The parameters b V th and b E L (see Table 3) were chosen based on Humphries et al. (2009). Although Planert et al. (2013) suggested that dopamine concentration modulates the excitability of D2-SPN, in low-dopamine state no significant changes in their excitability were observed. Therefore, following the reasoning given by Lindahl and Kotaleski (2016) in this model we also ignored the effects of dopamine on the D2-SPNs. However, to test whether this assumption affects our results, we simulated dopamine induced change in D2-SPNs properties and measured the transient response (see Table 13). We confirmed that dopamine modulation of D2-SPN has a negligible effect on the transient response.
We modelled the dopaminergic depolarizing effect induced through D1 type receptor activation on the FSIs, by modulating their resting membrane potential: where b E L (see Table 5) was set such that E FSI L at low dopamine level was 5 mV lower than that of the high dopamine level (Bracci et al., 2002).
The dopaminergic depolarizing effects on the GPe neurons (both TA and TI) are manifested as up-regulation of the hyperpolarization-activated cyclic nucleotide-gated (HCN) channels (Chan et al., 2011) which essentially results in a change in the resting membrane potential of the neurons. To mimic this effect, we changed the resting membrane potential of the GPe neurons in the following manner: The values of b E L for both the GPe-TA and GPe-TI neurons (see Tables 6, 7) were set such that the resting state potential of the GPe neurons at low dopamine level was 10 mV lower than that of its value at high dopamine level.
Dopaminergic effect on the SNr neurons (Zhou et al., 2009) was realized by changing their resting membrane potential: where b E L (see Table 9) was taken such that the resting potential at low dopamine level was 5 mV lower than its value at high dopamine level. The scaling factors b i ði 2 fE L ; V th g), for the linear modulation (f ða dop Þ ¼ a dop À a normal ) were tuned for each parameter to match their experimentally reported results (Lindahl and Kotaleski, 2016) in both normal and PD conditions (rodent models).

Dopamine effects on synaptic weights
High dopamine strengthens cortical projection on to D1-SPN and weakens cortical projections on to D2-SPN (Hernández-Echeagaray et al., 2004). The decrease in connectivity both in terms of synaptic strength and number of recurrent connections among SPNs is also attributed to dopamine depletion (Taverna et al., 2008). In addition, dopamine depletion is reported to enhance the strength of GABAergic synapses (Bracci et al., 2002) between FSI-FSI and increases the number of connections between FSI and D2-SPN (Gittis et al., 2011), but not D1-SPN. Within the GPe, dopamine depletion strengthens the GPe $GPe (Miguelez et al., 2012) and GPe !FSI connections (Bracci et al., 2002). In addition to that, it also strengthens the GPe-TA !SPN synapses (Glajch et al., 2016).
Dopamine depletion also strengthens the D2-SPN projections on to GPe neurons through reduced D2-receptor activation (Chuhma et al., 2011). Similarly, reduced dopamine concentration strengthens the STN !GPe synapses  and also responsible for increasing the synaptic efficacy in GPe-TI !STN synapses (Baufreton and Bevan, 2008). Galvan and Wichmann (2008) and Chu et al. (2017) claimed that cortico-STN transmission is reduced because of dopamine loss, but, Shen and Johnson (2000) suggested strengthening of the cortico-STN synapses and enhancement of the responsiveness of cortico-STN-SNr pathway at low dopamine level. Experimental data Sano and Nambu, 2019;Chiken et al., 2021;Wahyu et al., 2021) also reported that the strength of early excitation zone of the transient response in PD condition is either comparable to or much stronger than healthy state. This could be caused either by strengthening of cortico-STN synapses keeping STN !SNr synaptic property unchanged or by weakening of the cortico-STN synapses but increasing the STN to SNr weight. As dopamine receptor D1 and D2 activation induce two opposing effects, i.e., facilitating and depressing, respectively, on STN-SNr EPSC , it is not well understood how dopamine depletion affects the same in PD state (Lindahl and Kotaleski, 2016). Thus, in our model, we did not change the strength of STN !SNr synapse in low dopamine state and modeled the dopamine depletion induced changes in cortico-STN synapse by increasing the synaptic weight (Holgado et al., 2010;Lindahl and Kotaleski, 2016). On the other hand, at low dopamine, the D1-SPN to SNr connection strength was reduced (Chuhma et al., 2011), therefore, I GABA from D1-SPN to SNr was modelled to reflect the same.
Dopaminergic effect on the synaptic strength (g x!y syn ) was modelled as g x!y syn ¼ g x!y syn ð11b x y f Þ, where x, y [ {FSI, D1-SPN, D2-SPN, STN, Cortex, GPe, SNr} and the values of b x y were given in the Table 10.

External inputs
In our network model, all the neuronal populations received uncorrelated excitatory Poisson input spike-train. This input was provided to obtain baseline firing rate in the neurons. For the striatum this input corresponds to the cortical and thalamic inputs. For the STN, this input corresponds to the cortical inputs. For the GPe and SNr neurons, this input may correspond to either the endogenous activity or other inputs from outside the BG. Each neuron in a given population received a different realization of Poisson type spikes with the same parameters. The input rates were tuned both in normal and PD conditions to ensure that the basal firing rates (FR) of different subnuclei were consistent with the in vivo recordings in anaesthetized rats. For example, in normal condition baseline firing rate (in Hz) of D1-SPN and D2-SPN [ [0.01, 2.0] (Miller et al., 2008;Lindahl and Kotaleski, 2016), FSI [ [10, 20] (Gage et al., 2010), STN [ [10,13] (Fujimoto and Kita, 1993;Paz et al., 2005), and SNr [20,35] Benhamou and Cohen, 2014). The baseline activities of GPe-TA (11.8 6 1.1 Hz) and GPe-TI (24.2 6 0.7 Hz) is consistent with experimental data Mallet et al. (2008Mallet et al. ( , 2012. Similarly, in PD condition, frequencies of the background noise (in Hz) were also tuned to achieve range of basal firing rate of D1-SPN [ [0.1, 0.5 (Ruskin et al., 2002;Kita and Kita, 2011;Wahyu et al., 2021) had not observed firing rate changes in PD state. Given this, we kept basal firing rate of SNr the same as it is in normal state.
To characterize the effect of a transient cortical stimulation on the neuronal responses of the SNr, we stimulated striatal and STN neurons with a brief stimulus which amounted to injection of a rate modulated Poisson spiketrain ( Fig. 3, top panel). The fraction of stimulated neurons is specified in corresponding figures in Results.
This input was modelled by using the inhomogeneous_-poisson_generator device in NEST (Gewaltig and Diesmann, 2007). Because the transient stimulation was modelled as injection of Poisson spike-train for a brief period of time, we could control the strength of input stimulation by varying the amplitude of the EPSP generated by the injected spike train. Moreover, this allowed us to modulate the strength of input in relation to dopamine levels (see above, Dopamine effects on synaptic weights, for how dopamine affected synaptic weights). Transient response was measured in both normal and PD conditions.
Main differences between our model and the one by Lindahl and Kotaleski (2016) Here we build on the model by Lindahl and Kotaleski (2016); however, we made a few changes in the neuron and synapse models and changed the number of neurons in some of the BG subnetworks. The main focus of our work is to investigate how the structure of connectivity within and between BG subnetworks. Therefore, it was important to reduce the model complexity when possible.
Unlike their model (Lindahl and Kotaleski, 2016), striatal and STN neurons were modelled as simple LIF neurons without any kind of adaptation and, all the synapses were static as opposed to the dynamic ones. As we have argued later, this simplification has no specific effect on our key results. Striatal SPNs spike at low firing rate in an asynchronous manner, despite their recurrent connectivity and inputs. It is not necessary to model the striatum with %75,000 neurons as was done by Lindahl et al. (2013). When the parameters are appropriately scaled, we can obtain low firing rate asynchronous activity in a network of 4000 neurons. Therefore, we also reduced the size of striatal neuronal population. To this end, we changed the synaptic strengths and a few neuronal model parameters, such that the average synaptic input to GPe and GPi/SNr neurons was identical to the model used by Lindahl and Kotaleski (2016). This ensured that the model had the same repertoire of dynamical states as that of the model by Lindahl and Kotaleski (2016). Finally, to generate triphasic shaped transient responses, we also changed the values of b x y ðf Þ (see Table 10). Besides these changes, we followed the model closely while modeling the effects of dopamine on neuron and synapse parameters.

Limitations of the model
Unlike in the experimental data, in our model, all neurons responded with similar response profile. This is because the model is homogeneous in terms of neuron and synapse properties. It was important to keep the model homogeneous to isolate the various interactions that lead to triphasic or other shapes of transient response. Furthermore, all synapses are static in this model. We note that Lindahl et al. (2013) suggested that synaptic short term plasticity is important for the triphasic response when STN is stimulated. However, as we show in this study, triphasic responses do not require synaptic short-term dynamics. Moreover, to the best of our knowledge, there is no experimental evidence for short-term plasticity to be the cause of triphasic response. Furthermore, given that the short-term dynamics time constants are of the order of 100 ms, the effects of short-term dynamics may not be strong in the b band. It is possible that short-term dynamics of synapses may affect the properties of b band bursts but at least in this study we do not investigate such transient oscillation.
We have only considered effects of changes in the dopamine baseline. Transient responses of BG could also be accompanied by phasic change in DA levels. Such effects have been ignored. We did not explicitly model the effect of low-dopamine on D2-SPNs by changing the neuron and synapse properties. Instead, we mimicked this effect indirectly by increasing their baseline activity in PD condition. Next, we have only modelled the FSIs and ignored other types of interneurons. Only recently a detailed microcircuit has been modelled with numerical simulations (Hjorth et al., 2020). In future studies it may be possible to use a reduced version of that network for BG modeling. Finally, our model does not address the changes in the spatiotemporal dynamics of BG nuclei given cortical stimulation, as connectivity within each subnetwork is independent of spatial distances among the neurons.

Simulation tools
All the simulations were performed using the simulator NEST 2.20 (Jordan et al., 2019). All differential equations were integrated using Runga-Kutta method with a time step of 0.1 ms.

Code accessibility
The code to simulate the BG network is available on GitHub: https://github.com/arvkumar/Basal-Ganglia-Transient-Response. The Readme. _ txt file provides the necessary instructions to run the code. The simulation code was written in Python 3.0 and requires NEST 2.20 (Jordan et al., 2019) to run. The code is also provided as Extended Data 1 file.

Data analysis Transient response analysis
To get better estimate of the transient response we performed 100 trials and recorded the response over 1200 ms. The onset of the transient input (T stimulation ) was randomly chosen between 700 and 900 ms for every trial. Note that, the stimulation point was chosen between 700 and 900 ms to discard the initial transient that appears at the beginning of the simulation. To understand the effect of the transient stimulation on the SNr activities, the neuronal responses of SNr neurons were observed before and after the cortical stimulation point T stimulation . A 350-ms window size was defined around T stimulation to extract responses from each trial. For this, we used a time window of 100 ms before and 250 ms after the stimulation point. T stimulation is marked as 0 in the figures mentioned in the subsequent sections.
The responses were evaluated by constructing peristimulus time histograms (PSTH), using 1-ms rectangular bins for each trial data. The baseline neuronal activity was computed based on the mean and SD of the PSTH data for 100-ms duration before the onset of the stimulation point.
Visual inspection of the experimental measurement of the transient responses shows it has actually four distinct phases (schematically shown in Fig. 2). However, strangely in the literature such response has been referred to as triphasic. Here, to analyze the transient response of the SNr, both in normal and PD conditions, the PSTH were divided into four zones (see Fig. 2) based on the change in firing rate. These zones consisted of two excitatory (EE and LE for early and late, respectively) and two inhibitory (EI and LI) zones. The change, i.e., increase or decrease in firing activity was marked as excitation or inhibition, respectively, if the firing rate was significantly above or below the baseline (p , 0.05, one-tailed Z-test) discharge rate for at least two consecutive time bins (2 ms; Sano et al., 2013). The latency of each zone was measured as the time when the first bin exceeded the baseline. Similarly, the zone terminated when activities during two consecutive bins fell below the significance level. The end time was determined as the time of the last bin exceeding the significance level. The total time duration from the first bin to the last bin (of a significant response) was considered as the duration of each zone. The sum of heights of bins within a particular zone is considered as the area as well as the strength of the zone, whereas the area per unit time (area/time) indicates the average strength of that zone. Thus, we extracted the following features from PSTH for each zone: latency (L), duration (D), absolute area indicating strength (A) of that zone, mean (H m ), and SD (H s ) of bin heights. In addition, we also measured the peak amplitudes (H p ) of each zone (i.e., H p 2 fH EE max ; H LE max ; H EI min ; H LI min g).
. F i s were computed for each zone in different network conditions (F for all four zones). Such networks, for example, were simulated by restoring selected synaptic strengths from PD condition to the normal condition. These can be referred to as TestNetworks. The similarity between such a TestNetwork and networks tuned in healthy condition and PD condition was calculated using Euclidean distance metric: In order to observe the statistical variation of the above features within the normal and PD conditions, a subpopulation of SNr neurons was considered. From the whole population of SNr neurons, a percentage of neurons (N S %) was randomly chosen to represent an observation.
The responses of these subpopulation of SNr neurons were averaged over multiple trials (O S = 100 in the present simulation) and then the features of the four zones were extracted. The choice of the subpopulation of N S % was also varied for every observation and a large number of such observations (O S ) were made. These observations were used to derive the mean and SD of the features for each of the above-mentioned EE, EI, LE, and LI zones. Here, we have considered O S = 100 and N S = 50%.
In experimental studies (Ozaki et al., 2017;Sano and Nambu, 2019), in healthy state, the transient response is often characterized by dividing the response in three zones EE, EI, and LE. In low-dopamine state the experimentally observed biphasic response pattern could consist of any two zones out of the three zones (EE, EI, and LE; Sano and Nambu, 2019;Wahyu et al., 2021). Experimental data also shows that in both normal and PD conditions, LE is followed by an LI zone . Therefore, here we defined four zones to characterize the transient response in the healthy state ( Fig. 2A). Thus, although we have defined four zones, we still refer to it as a triphasic response to be consistent with the terminology used in the literature. In our simulations, biphasic response (Fig. 2B) observed in PD condition consisted of EE and LE, while the EI zone was not observed. We In these cases, the EI was not detectable using statistical test, hence the two excitations (EE and LE) were merged during computation of the parameters. ND denotes that the zone was not detected using significance test (see Materials and Methods, Data analysis). merged the EE and LE zones together for computing the features of the PD-biphasic response.

Global network activity
The oscillatory behavior of population activities was assessed in PD as well as normal conditions. We ran the simulation for longer duration (5 s), without any transient input, to allow the oscillations to set in to their steady state. This experiment was also conducted over 100 trials.
Synchrony in the firing rates of a neuronal population was estimated using Fano factor (FF pop ; Kumar et al., 2011): where E pop and V pop are the mean and variance of the neuronal activity for the same population, respectively. For an uncorrelated ensemble of Poisson processes, FF pop = 1 and when neurons tend to correlate, FF pop . 1.
Here, we binned the neuronal activity using rectangular bins of 3-ms duration. This window size was similar to the one used in previous studies (Mallet et al., 2008;Lindahl and Kotaleski, 2016).
To determine the strength of oscillatory neuronal activities in the b band, we estimated the oscillation index (OI pop ). To this end, we estimated the spectrum of the population activity [S pop (f)]. As we used 3 ms bins to calculate the PSTH, the sampling frequency (F s ) was 333.3 Hz. To estimate the oscillation index, we measured the relative power confined in the b band: The phase relationships of the firing patterns between the two types (GPe-TA and GPe-TI) of GPe nuclei as well as with STN were computed from the PSTH, having bin size of 1 ms. As we were interested in analyzing the pathologic b oscillation, the individual PSTH responses were bandpass filtered between 12 and 30 Hz.
Initially, at every time instance, corresponding to each bin, the instantaneous phase was calculated using the Hilbert transform. Then the differences of the instantaneous phases were obtained between a pair of nuclei, for every 1 ms. Finally, the histogram of the difference in the phase was obtained with 100 bins in the range of 0 to p .

Results
The standard feedforward model of the BG (Albin et al., 1989) predicts that transient cortical stimulation Figure 5. Effect of strength of cortical stimulation on BG transient response shape. To vary the strength of cortical stimulation we varied the fraction of striatal and STN populations that received cortical inputs from 10% to 100%. A, Average transient response (100 trials) in SNr in normal (blue color) and PD state (orange color). Lighter (darker) color-shades indicate smaller (larger) size of stimulated population. Note that in PD condition, even the strongest cortical input failed to elicit a response similar to that seen in healthy state. B, Changes in the duration of the four zones of the transient response in normal state. C, Changes in the area per unit time (area/time) of the four zones of the transient response in normal state. D, Same as in B but for PD state when the network responded with a triphasic response. E, Same as in C but for PD state when the network responded with a triphasic response. Note missing colors in a given bar implies that we could not detect the corresponding zone.
will result in what has been denoted as the triphasic response in the SNr as the stimulus induced activity is propagated over the hyperdirect, direct and indirect pathways. Indeed, many neurons, at least in a healthy state, do show a triphasic response in vivo. However, in both healthy and dopamine-depleted conditions, response pattern of a sizeable fraction of neurons deviates from the triphasic response shape Sano and Nambu, 2019;Chiken et al., 2021;Wahyu et al., 2021) indicating the role of recurrent interactions within and between BG nuclei. To understand how different neurons and network parameters shape the output of SNr when the striatum and STN are transiently stimulated, we used numerical simulations of the BG network with spiking neurons. In the model, we systematically varied the dopamine level and studied how strength of different connections in the BG affects the shape of the transient response in both healthy and PD conditions. Here, we set the dopamine level to 0.8 and 0.0 to tune the model into healthy and PD conditions, respectively (Lindahl and Kotaleski, 2016).  Cortically evoked transient response in SNr Transient stimulation of the neocortex results in a brief excitation followed by inhibition caused by the recurrent inhibition. Therefore, we stimulated the striatal and STN neurons with a rate modulated Poisson process that mimicked the excitation-inhibition pattern cortical stimulus response (see Fig. 3; Materials and Methods). Consistent with the predictions of a feedforward model of the BG and in vivo experimental data, in healthy state SNr neurons responded with a triphasic response consisting of early excitation (because of STN), inhibition (because of the D1-SPN projections), and late excitation (because of indirect pathway), i.e., the EE-EI-LE response (see Fig. 3A). By contrast, in PD condition, SNr neurons responded with a biphasic response (from here on referred to as default PD condition), consisting of a prominent early excitation and late excitation (i.e., EE-LE; see Fig. 3B). Thus, the model suggests that persistent dopamine depletion (see Materials and Methods; Table 10) not only affects the steady-state of the BG network (i.e., b band oscillations)    (Gertler et al., 2008) Reset value for v_m after spike V_th À45 mV (Bahuguna et al., 2015) Spike threshold tau_syn_ex 0.3 ms (Bahuguna et al., 2015) Rise time of excitatory synaptic conductance tau_syn_in 2 ms (Bahuguna et al., 2015) Rise time of inhibitory synaptic conductance E_L À85.4 mV Leak reversal potential E_ex 0 mV Excitatory reversal potential E_in À64 mV Inhibitory reversal potential I_e 0 pA Constant input current C_m 157 pF (Gertler et al., 2008) Membrane capacitance g_L 6.46 nS (Gertler et al., 2008) Leak conductance t_ref 2 ms Duration of refractory period but also impair the transient inhibitory effect of the striatal inputs to the SNr because of weak D1-SPN projections as well as stronger activity along the indirect pathway. In PD conditions, the missing EI phase is attributed to the weaker D1-SPN !SNr connection and the stronger inhibitory influence of the D2-SPN !GPe-TI connection which disinhibited the SNr neurons and subsequently resulted in the prolonged LE phase. In fact, in our simulations, in extreme PD conditions (when a dop = 0) D1-SPN response was almost zero. Experimental data show that even in the PD condition, 15-40% of the SNr neurons respond in a triphasic manner (Sano and Nambu, 2019;Wahyu et al., 2021). In our model, to generate a triphasic response in PD condition (Fig. 3C), we needed to make additional changes other than those brought in by low dopamine. In particular, we decreased D2-SPN !GPe-TI, increased D1-SPN !SNr, and decreased GPe-TI !STN connections (see Table 10 for numerical values). Note that, although these synaptic changes in the network parameters affected the shape of the transient response but did not affect the oscillation and synchrony in both the PD states (OI for PDbiphasic = 0.24 and PD-triphasic = 0.23, FF for PDbiphasic = 18.01 and PD-triphasic = 13.49; see Fig. 8 for more details). In other words, small changes in the synaptic connectivity can affect the transient response without changing the baseline activity in a qualitative manner.
In the above, we ignored that dopamine may also affect the properties of D2-SPNs and thereby may affect the transient response. To check whether dopamine dependent modulation of D2-SPNs may affect the transient response we changed the excitability of D2-SPNs as a function of dopamine levels (Day et al., 2008;Damodaran et al., 2015). However, these changed in D2-SPNs did not affect the shape of the transient response (see Table 13) in either of the two PD states.
Next, we also tested how progressive change in the dopamine levels may affect the shape of the transient response. To this end, we tuned the model in either PD biphasic or PD triphasic state and systematically increased the level of dopamine. We found that the four phases of the transient response change gradually as a function of dopamine level (Fig. 4A,B) in PD-triphasic condition. By contrast, in PD-biphasic state, there appears to be a threshold below which EI is not detectable (Fig. 4C,  D).
While qualitatively in both healthy and PD conditions cortical stimulation evoked an early excitation but in PD condition (both biphasic and triphasic) the duration and amplitude of early excitation were higher than that of in the healthy condition. This was because dopamine depletion amplified the excitation of SNr neurons through the hyperdirect pathway. Moreover, in PD condition when we could generate triphasic response pattern, the duration  and amplitude of early inhibition (i.e., EI) were much smaller than that observed in the healthy condition (early inhibition was completely absent in the biphasic responses). Finally, the late excitation phase (LE) of the triphasic response was longer in PD condition than in the healthy condition. The details of further differences in transient response properties are provided in the Table 11. The trend of the features in normal and PD conditions is consistent with the experimental data (Ozaki et al., 2017;Sano and Nambu, 2019;Chiken et al., 2021;Wahyu et al., 2021). These results suggest that dopamine depletion primarily affected the EI and LE zones. On one hand dopamine depletion reduced excitability of D1-SPN (Gruber et al., 2003) and reduced basal firing in D1-SPN while increasing in firing rate of D2-SPNs. Therefore, the direct pathway was weakened and resulted in reduced EI in the SNr. On the other hand the basal firing rate of GPe-TI was reduced and GPe-TA was increased (Mallet et al., 2008) because of the strengthening of the striato-pallidal pathway. This resulted in prolonged "LE" zone.

STN evoked transient response in the SNr
To separate the contribution of direct and hyperdirect pathways we measured the SNr response when only STN was stimulated (Fig. 3, bottom row, green trace). In healthy state, consistent with experimental data (Maurice et al., 2003) and a previous modeling study (Lindahl et al., 2013), STN stimulation alone generated a triphasic response in the SNr; however, there were notable differences: the EI zone was weaker, LE zone was both weaker and delayed and, LI zone was absent. In this type of stimulation codition, STN to SNr connections shaped the EE zone, and STN$GPe-TI connections shaped the EI zone and the LE zone. Because cortical inputs to the striatum promote inhibition in the GPe-TI and SNr, their removal  To obtain a triphasic response in PD condition, we had to change a few parameters of the network tuned in default PD state (biphasic). These changes are marked in boldface.
made the EI and LE zone much weaker and diminished the LI zone. In PD condition, STN stimulation induced a transient response with only EE zone (essentially, EE and LE zones observed in a healthy state were merged into a single excitatory zone). The magnitude of the EE zone was much higher in the PD triphasic configuration as compared with the normal state because of the stronger hyperdirect pathway. These results confirm that in healthy state the hyperdirect pathway shapes the EE zone, and suggest that the contribution of the striatal activity (direct pathways) is in increasing the magnitude while decreasing the duration of the EI and LE zones.

Effect of the strength of cortical stimulation on the transient response
The aforementioned transient responses were measured by stimulating 50% of the striatal and STN population. Next, we asked whether differences in the shape of the triphasic response observed in PD and healthy conditions could be reduced by stimulating more neurons. To this end, we systematically increased the number of striatal and STN neurons that received cortical stimulation (to mimic the strength of cortical stimulation). To quantify the changes in the shape of the triphasic responses we measured the duration and area per unit time of the four zones in healthy (Fig. 5B,C) and PD condition (Fig. 5D,E).
We found that in both healthy and PD conditions, the amplitude of the four zones are monotonically increased and saturated at a maximum value (Fig. 5A). On the other hand, area per unit time of the excitatory zones monotonically increased. Thus, an increase in the stimulus strength increased both excitation and induced stronger inhibition in healthy state. However, in healthy state a very weak cortical stimulation (10%) failed to elicit a detectable EE response (Fig. 5B,C) but in PD condition (Fig. 5D,E), the same weak stimulation elicited a strong EE response, again indicated the strengthening of the hyperdirect pathway in low-dopamine state.
Overall, these results show that even with the strongest stimulation in PD condition, we could not reproduce the transient response properties observed in healthy state even with the weakest cortical stimulation. This suggests that the differences in the transient response are not simply because of the altered cortico-BG projections but are primarily because of the altered connectivity within the BG.

Effect of change in synaptic connections on cortical evoked transient response in SNr
Above, we demonstrated the transient response pattern for a specific combination of synaptic strengths. The total space of different synaptic parameters is 22-dimensional (Table 2), and therefore, it is not feasible to test the robustness of our results in a systematic manner by varying all the connection parameters. The structure of BG connectivity suggests that the triphasic response pattern is shaped by D1-SPN !SNr (early inhibition), GPe-TA $GPe-TI, STN $GPe-TI, and D2-SPN !GPe-TI (late excitation/inhibition) connectivity. Therefore, we individually varied these six connections and quantified the duration and area per unit time (area/time) of the four zones of the transient response. The minimum and maximum values of each of the synaptic weight (except D1-SPN!SNr synapses) corresponded to their values in L-dopa-induced Early inhibition (EI) *ND Latency (ms) 11.0 6 0, 11.0 6 0 12.0 6 0, 11.97 6 0.17 Duration (ms) 6.06 0.0, 5.99 6 0. Here, the variations in the features were obtained using multiple observations (100 in number) of the simulation output. In each observation, 50% and 20% of SNr neurons were randomly chosen whose results are shown in bold and italic, respectively. The format of the result is mean 6 std. It can be seen that when 20% of neurons are used to compute the statistics, then the SD (std) increases, as expected because of the reduced number of neurons. Such a statistics is useful to study the variations that the population neuron brings into the network. ##the statistics for EE in PD-biphasic are given for the complete excitatory response comprising of both EE and LE. In this case, the EI was not detectable using statistical test, hence the two excitations (EE and LE) were merged during computation of the parameters. Here, the deviation of peak amplitude (H p ) was measurable with respect to baseline (H bas ). *ND denotes that the zone was not detected using significance test (see Materials and Methods, Data analysis). dyskinesia (LID) and PD conditions, respectively (see Materials and Methods; Table 10). For the case of D1-SPN !SNr synapses minimum and maximum values corresponded to PD condition and LID, respectively (see Table 10). To mimic low-dopamine and high-dopamine states, we varied the connection strengths in the following manner. Let's consider that the synaptic weight in normal condition is v and in PD condition is mv (where the scaling factor m was derived from Table  10). To mimic high-dopamine (greenish colors) we reduced v to v/m in three steps. Similarly, to mimic lowdopamine (reddish colors) we increased v to mv in three steps. Thus, the simulations were done for seven different configurations of synaptic strengths including the normal. We found that synaptic weight changes can affect the duration of the LE and LI zones, but not of EE and EI zones (Fig. 6A-F). By contrast, the area per unit time of all four zones was sensitive to synaptic weight changes ( Fig.  6G-L). Our simulation experiments showed that mainly the connection between the D1-SPN !SNr (Fig. 6A,G), the STN !GPe-TI (Fig. 6E,K), and D2-SPN !GPe-TI (Fig.  6B,H) controlled the shape of the three zones namely EI, LE, and LI. To analyze the relative contribution of these  Table 11; however, by pooling synaptic weights corresponding to Figure 6 Normal ( Table 11) PD-biphasic PD-triphasic (same as Table 11 Column 1 and column 3 are taken from Table 11. Here, the variations in the features were obtained using multiple observations (100 in number) of the simulation output. In each observation, 50% and 20% of SNr neurons were randomly chosen whose results are shown in bold and italic, respectively. ##: Statistics for EE in PD-biphasic are given for the complete excitatory response comprising of both EE and LE. In this case, the EI was not detectable using statistical test, hence the two excitations (EE and LE) were merged during computation of the parameters. Here, H p À H bas denotes the deviation of peak amplitude (H p ) with respect to baseline (H bas ). *ND denotes that the zone was not detected using significance test (see Materials and Methods, Data analysis).
connections, we calculated the coefficient of variation (CV) of duration (Fig. 6M) and area/time (Fig. 6N) of the four zones for each of the six connections we changed. Smaller the CV values, smaller is the influence of that connection on the transient response properties. From this analysis the D2-SPN !GPe-TI connection emerged as the most crucial parameter in shaping the transient response in both low and high-dopamine conditions. For extreme values of D2-SPN !GPe-TI connection, area/time of late excitation was very high in low dopamine condition and late inhibition zone was completely absent in high dopamine state. Besides D2-SPN!GPe-TI connection, STN !GPe-TI connection was the second most important parameter as it affected the EI, LE, and in particular, the LI zones. Finally, it was somewhat surprising that the GPe-TI $GPe-TA and GPe-TI !STN connections did not affect the area/ time of any of the zones (Fig. 6I,J,L).

Effect of restoration of dopaminergic synaptic connection on the transient response in SNr
To get further insights into the network mechanisms underlying the aberrant transient response, we asked whether we could restore the healthy state of transient response by restoring specific connections to their healthy levels. In the previous section we showed that the D1-SPN !SNr, STN !GPe-TI, and D2-SPN !GPe-TI (Fig. 6M,N) have the strongest effect on the triphasic response. To further confirm their role in the generation of aberrant transient responses, we restore selected connections along the indirect pathway (i.e., D2-SPN !GPe-TI or STN $GPe-TI loop, or GPe-TA $GPe-TI).
To this end, first, we tuned the BG network in the PD state (Table 10) such that the SNr shows a triphasic response. Then restored the strength of D2-SPN !TI, STN $GPe-TI loop, GPe-TA $GPe-TI, and D1-SPN !SNr one by one. During the restoration of the synaptic connection between a We changed the connection strength from v/m to mv in seven steps (for more information, see Materials and Methods). Cortical input was given to 50% of the striatal and STN populations. A-F, Change in duration of the four zones as the strength of one of the connections was increased (3 values) or decreased (3 values) to mimic low and high dopamine states. The duration of each zone was normalized with respect to duration in the normal condition for that particular zone (beige colored bars). G-L, Same as in panels A-F but for area per unit time for each of the four zones. LI: late inhibition, LE: late excitation, EI: early inhibition, EE: early excitation; $: not detected. M, CV of the duration of the four zones duration as a function of the six synaptic connections (computed from A-F). High value of CV means that changes in the particular connection result in higher changes (increase/decrease) in the duration of a zone. N, Same as in panel M but for area per unit time.
pair of nuclei, (1) the synaptic weights and delays were made equal to normal, (2) the basal firing rates were made similar to the normal by changing the background firing rate or the background current, (3) basal firing of SNr was kept the same as that of the normal condition. To compare the triphasic response in healthy and PD states (with and without restoration of certain synaptic weights) we measured the distance between two network conditions (Eq. 5; see Materials and Methods, Data analysis).
We found that restoring the D2-SPN !GPe-TI alone is sufficient to bring the shape of the transient response close to the one observed in healthy state (Fig.   7A,E). However, by restoring the GPe-TI $GPe-TA or STN $GPe-TI connections, we only restored the shape of EI and LE zone but not that of the EE and LI zone (Fig. 7B,C). By contrast, restoration of the D1-SPN !SNr synaptic connection made the network activities different from as in both PD and healthy conditions (Fig. 7D,E). This is because in PD condition, in addition to the weakening of D1-SPN to SNr synapses, cortical inputs to D1-SPN were also weakened (Lindahl and Kotaleski, 2016), and therefore, cortical stimulation only evokes a weak response in D1-SPN. Thus, although restoration of the D2-SPN !GPe-TI or , and GPe-TA $GPe-TI connections after dopamine depletion. Cortical input was given to 50% of the striatal and STN neurons. A, Comparison of zone wise feature values for test network condition with respect to the healthy and PD states. Here, the test network refers to a PD network (tuned to generate triphasic response) in which synaptic weight of D2-SPN !GPe-TI was only restored to their healthy value (restored connection, D2-TI). Feature values for healthy condition is represented as unit circle (shown in green). L EE , D EE , A EE , H EE m ; H EE s , and H EE max denote the latency, duration, area, mean of bin heights, SD of bin heights, and maximum of bin heights for EE zone, respectively (for feature description, see Materials and Methods, Data analysis). L EI , D EI , A EI , H EI m ; H EI s , and H EI min denote the latency, duration, area, mean of bin heights, SD of bin heights, and minimum of bin heights for EI zone, respectively. L LE , D LE , A LE , H LE m ; H LE s , and H LE max denote the latency, duration, area, mean of bin heights, SD of bin heights, and maximum of bin heights for LE zone, respectively. L LI , D LI , A LI , H LI m ; H LI s , and H LI min denote the latency, duration, area, mean of bin heights, SD of bin heights, and minimum of bin heights for LI zone, respectively. B, Same as A but when GPe-TA $GPe-TI connections were restored (restored connection, TA-TI collaterals). C, Same as A but when STN $GPe-TI loop were restored (restored connection, STN-TI loop). D, Same as A but when D1-SPN !SNr was restored (restored connection, D1-SNr). E, Distance computed for different test network conditions based on A-D from healthy and PD network state. The green bars indicate the distance (calculated using Eq. 5) between healthy network activity state and test network activity state, whereas red bars indicate the distance between PD network activity state and test network activity state. Here, the test network refers to a PD network (tuned to generate triphasic response) in which individual synaptic weights (mentioned on the x-axis) were restored to their healthy value. STN $GPe-TI makes the transient response in PD condition more similar to the healthy condition, the early inhibitory phase is not restored.

Ongoing spontaneous activity of BG network
In the steady states, the b band oscillations and synchrony within and between different BG subnuclei during stimulus free ongoing spontaneous activity are two prominent hallmarks of PD condition in vivo (Brown et al., 2001;Mallet et al., 2006Mallet et al., , 2008. Therefore, next we tested whether the network parameters we used to generate the aberrant triphasic and biphasic responses could also induce b band oscillations. To this end, we tuned the BG network in PD condition when it showed either triphasic or biphasic transient response and measured the oscillations and synchrony in the ongoing (stimulus-free or spontaneous) activity.
We found that indeed the same set of parameters that generated aberrant transient responses were sufficient to elicit clear b band oscillations in both the biphasic and triphasic response modes (Fig. 8A-D). Next, we measured the phase relationship between different subnuclei of the BG. Mallet et al. (2008) reported that there exists an inphase relationship between activities of GPe-TA and STN neurons and anti-phase relationship between GPe-TA and GPe-TI neurons. In our model the phase relationships between GPe-TA and GPe-TI, GPe-TA and STN, GPe-TI and STN (Fig. 8E-G) were similar to that observed in experimental data. Thus, these results suggest that similar changes in the network connection could underlie the aberrant transient response and ongoing activity in PD conditions as well.
Effect of striato-pallidal and pallido-subthalamic pathways on the b oscillations While b band oscillations are a clear neural signature of PD, the mechanisms underlying the emergence of these oscillations are still debated. Both experimental data (Plenz and Kital, 1999;Hammond et al., 2007;Tachibana et al., 2011;de la Crompe et al., 2020) and computational models (Kumar et al., 2011;Holgado et al., 2010;Tachibana et al., 2011;Pavlides et al., 2015;Corbit et al., 2016;Bahuguna et al., 2020) have implicated essentially all the various network interactions in generating oscillations. Here, we have developed the BG model primarily to understand the transient response and found that the same model can also generate b band oscillations. Thus, we have a more constrained model of the BG than used previously, and this could help us narrow down on the key determinants of oscillations.
Based on our simulations and available experimental data (de la Crompe et al., 2020) GPe has emerged as a key network necessary to induce b band oscillations. Figure 8. b band oscillations in the ongoing activity of the BG. A, Spectrum of the GPe-TA activity in healthy (green), PD-triphasic (red), and PD-biphasic (brown) response conditions. B, Same as in panel A but for GPe-TI. C, Same as in panel A but for STN. D, Same as in panel A but for SNr. E, Phase relation between GPe-TA and GPe-TI shown in the range of 0 and p (in radians). F, Same as E but for the phase relation between GPe-TA and STN. G, Same as E but for the phase relation between GPe-TI and STN. The phase difference between GPe-TA and GPe-TI attains a peak around 2p 3 in the phase histogram shown in E. The in-phase relation between GPe-TA-STN and approximate anti-phase relation between GPe-TI-STN can be seen in F, G, respectively. However, it remains unclear which of its input and output connections are more crucial to generate oscillatory activity. Therefore, to quantify the relative contribution of GPe connectivity we computationally followed the path of lesion experiments usually done in vivo studies. To do this, we either removed striatal input to GPe-TI neurons, or GPe feedback to the striatal FSIs, or GPe-STN interactions. All these perturbations were performed in two different BG networks which showed biphasic or triphasic response in PD condition.
In both PD conditions (biphasic response and triphasic response), removal of D2-SPN input to GPe-TI neurons reduced the oscillations and synchrony GPe-TA, GPe-TI, STN, and SNr neurons nearly to a level observed in healthy state (Fig. 9, pale green bars). The convergence of a relatively large numbers of D2-SPNs onto a smaller set of the GPe-TI neurons greatly influenced the patterning of GPe-TI activity (Kovaleski et al., 2020). This supports the hypothesis that increase in D2-SPN activity in dopamine depleted state is responsible for unleashing oscillations in the BG (Mallet et al., 2006;Kumar et al., 2011;Sharott et al., 2017). Moreover, recent experiments also suggest that D2-SPN inputs control oscillations in the GPe-TI population (de la Crompe et al., 2020).
By contrast, removal of GPe feedback to striatal FSIs (Fig. 9, orange bars) did not affect oscillations or synchrony considerably. The effect of interactions within and between GPe-TA and GPe-TI neurons was dependent on the state of the network: removal of these connections (Fig. 9, yellow bars) reduced oscillations and synchrony of GPe-TA, GPe-TI, STN, and SNr neurons by a larger amount when the BG was tuned to exhibit triphasic response in PD condition (Fig. 9C,D, yellow bars). However, even after removal of collateral within the GPe neurons, both oscillations and synchrony were much higher than that observed in a healthy state.
Surprisingly, removal of STN $GPe-TI connections did not affect the oscillation, regardless of the type of transient responses the network showed in PD condition ( Fig.  9, blue bars). These findings imply that following dopamine loss, the abnormal patterning of GPe-TI by stronger striato-pallidal connection was opposed by STN-GPe transmission and thus, removal or inhibition of the STN will also not have any effect on the oscillations (Kovaleski et al., 2020). This is consistent with the recent findings by de la Crompe et al. (2020), who showed that optogenetic inhibition of STN does not quench oscillations (see also Gradinaru et al., 2009). Although experimental evidences suggest that the oscillatory dynamics in SNr in dopamine depleted state could be induced through either stronger D2-SPN !GPe-TI (Gradinaru et al., 2009;de la Crompe et al., 2020) or STN $GPe-TI loop (Holgado et al., 2010;Pavlides et al., 2015); however, in our model, the removal of the D2-SPN !GPe-TI connection is sufficient to suppress the oscillation in the SNr. The results emphasize the role the indirect pathway as opposed to the hyperdirect pathway in the manifestation of aberrant activity in PD.

Diversity of transient responses
As noted earlier, our BG network model is homogeneous and therefore, we could either generate biphasic or triphasic shaped transient response in the network. This approach however, allowed us to identify the key network interactions that are involved in changing the response shape from biphasic to triphasic (i.e., D2-SPN !GPe-TI, D1-SPN !SNr, and GPe-TI !STN). An inhomogeneous change in these connections could be one of the reasons for the observed diversity of transient responses in in vivo. However, oscillations in the ongoing activity could also contribute to the diversity of transient responses because the shape of transient response may depend on the oscillation phase at which cortical stimulation was delivered. In fact, recent experimental data suggests that when b band oscillations are weak or absent in PD, transient responses variability is reduced .
To test this hypothesis, we tuned the network in a PD state in which it responded with a biphasic shape (Table  10) and delivered the stimulus at different phases of oscillations. There were 48 trials of each cortical stimulus given at a specific phase of the SNr oscillation. We pooled the data for each of such 48 trials and observed various responses namely, "EE-EI-LE," "EE-EI," "EI-LE," "EE-LE," "EE," and "LE." Here, "EE-EI-LE" denotes a triphasic response which is observed mainly in healthy state.
This variation was primarily because of the differences in phase of the oscillation at which cortical stimulation was delivered. We found that when the input arrived at the trough (1.20p ; Fig. 10A) or during the falling edge of the SNr b oscillations (0.23p ; Fig. 10D), SNr responded with a biphasic transient response. The magnitude of the LE zone was stronger when the input arrived at the trough instead of the falling phase of the b oscillations (Fig. 10A,D).
By contrast, when the stimulation arrived during the rising phase, it resulted in triphasic responses in SNr (Fig.  10B,C), although the strength (area/time) and duration of the EE, EI, and LE varied with the actual phase of the stimulation (Fig. 10E,F). It is important to note that the triphasic response observed for some stimulation phases was still quantitatively different from that triphasic response seen in a healthy state.
In order to characterize the contribution of different BG nuclei to the transient response in healthy and PD conditions, we varied the strength of several connections in BG (e.g., see Figs. 5, 6). We pooled all those simulations together, where strengths of the synaptic connection were varied according to Figure 6, and estimated the variability of the transient responses. The rationale to do this was that each network simulation with different connection strength may represent a different SNr/GPi region or animal where the transient response was recorded. Indeed, such pooling of the data resulted in a high heterogeneity in the transient responses in both healthy and PD conditions (see Table  12), which closely matched with the experimental data.
These results while they do not explain the full diversity of the responses observed in in vivo, they show that the oscillation phase as well as diversity of synaptic connectivity are important variables in determining the shape of the response.

Discussion
Here, we have studied how the changes induced by low-dopamine affect both transient response (induced by cortical stimulation) as well as the ongoing spontaneous activity state of the BG network. Typically, a transient stimulation of the cortex results in a triphasic response in the SNr/GPi (the output of the BG). The shape of the response is impaired in chronic low-dopamine conditions such as PD. The different zones of the transient response can be associated with different aspects of initiation of voluntary movements. For instance, it has been hypothesized that EE zone resets the cortical activity, EI zone allows for the execution of movements and LE zone stops the movement (Nambu et al., 2002;Chiken et al., 2021). A weaker or completely absent EI zone in PD is thought to be related to akinesia. Indeed, L-dopa treatment or local inhibition of the STN both of which restore the EI zone also ameliorate motor deficits in PD . The triphasic response in the SNr/GPi is usually explained by difference in the relative timing of the hyperdirect, direct and indirect pathways of the BG which converge in the SNr/GPi.
Here, we show that changes in the shape of the transient response in PD state involve not only changes in the feed-forward connections between different subnuclei of the BG (D1-SPN !SNr) but also by interactions between STN and GPe (GPe-TI $STN; Fig. 6K) and to some extent by GPe-TA $GPe-TI (Fig. 7). Moreover, we show that same changes in the BG network (both synaptic and neuronal excitability) may underlie the impairment of transient response and emergence of induced population level oscillations and synchrony in the BG.
In the PD condition, neurons either show biphasic or triphasic transient responses (Sano and Nambu, 2019;Chiken et al., 2021;Wahyu et al., 2021); the latter is, however, quantitatively different from the triphasic response observed in healthy state. In our model, the aberrant biphasic response in PD condition appeared as we changed the parameters to a low-dopamine state (according to the model by Lindahl and Kotaleski, 2016). However, to obtain a triphasic response, we needed to reduce D2-SPN !GPe-TI, increase D1-SPN !SNr, and reduce GPe-TI !STN connections (see Table 10). This suggests that dopamine effects are not homogeneous within and between different subnuclei of BG. To restore a healthy state, it is important to experimentally characterize the heterogeneity of dopamine action. The diversity of dopamine action and phase of oscillations at which stimulation was delivered, together could explain the observed diversity of transient responses in in vivo.
Previously, Blenkinsop et al. (2017) suggested that in a healthy state, biphasic and triphasic responses in the SNr arise because of interactions among functionally segregated channels of competing inputs with different strengths. In a BG model with functionally segregated channels, local inhibition within the GPe and excitation from a small number of highly active STN neurons (presumably because of stronger cortical inputs) are responsible for emergence of LE zone, rendering a response biphasic or triphasic (Blenkinsop et al., 2017). Here, we have used a BG model without functionally segregated channels. Our results suggest that diversity of synaptic strengths within and between BG nuclei could give rise to some neurons responding in a triphasic and others in a biphasic manner. Consistent with the model by Blenkinsop et al. (2017), in our model, the magnitude of LE zone can be controlled by the strength of cortical stimulation (Fig. 5). Our work points out a strong influence of the indirect pathway (D2-SPN to GPe-TI) in controlling the shape of the transient response both in normal and PD condition. This observation is also consistent with the proposal of Blenkinsop et al. (2017). Figure 10. Diversity of transient response may depend on the phase of cortical stimulation. The network was tuned to operate in a PD state in which transient inputs result in a biphasic response in the SNr. Cortical input was given to 50% of the striatal and STN neurons. A, blue trace, Average transient response (average across 48 trials) when the stimulation arrived at trough of the b oscillations (U 1 = 1.20p ). It can be seen that during the stimulation, SNr is close to the trough of oscillation. The shaded areas above/ below the transient response (blue line) denote the 95% confidence interval. Green trace, The grand average (grand average) of the transient response averaged across all trials and stimulation phases. The corresponding transient response (shown in blue) with only EE and LE zones. B, Same as panel A; however, the stimulation arrived at a phase U 2 = 1.61p . The corresponding transient response (shown in dark red) is triphasic (EE-EI-LE). C, Same as A; however, the stimulation arrived at a phase U 3 = 1.81p , i.e., very close to the peak of SNr oscillation. The corresponding transient response (shown in magenta) is triphasic (EE-EI-LE); however, the EI zone is very weak. D, Same as A; however, the stimulation arrived at a phase U 4 = 0.23p , i.e., on the falling edge close to the trough of SNr oscillation. The corresponding transient response (shown in violet) is biphasic (EE-LE). E, Variation in the zone wise duration of the transient response for the stimulation corresponding to phase U 1 ; U 2 ; U 3 , and U 4 . F, Variation in the zone wise area per unit time (area/time) feature of the transient response for the stimulation corresponding to phase U 1 ; U 2 ; U 3 , and U 4 . ## The duration and area/time corresponding to phase U 1 ; U 4 and for the grand average are given for the complete excitatory response comprising of both EE and LE zones. In these cases, the EI was not detectable using statistical test, hence the two excitations (EE and LE) were merged during computation of the parameters. ND denotes that the zone was not detected using significance test (see Materials and Methods, Data analysis).
Although we managed to generate a triphasic response in PD condition, it was quantitatively different from the one observed in a healthy state. The differences were most clearly seen in the late excitation which was lasted longer in PD condition as compared with a healthy state. Moreover, these differences in the triphasic response could not be compensated by increasing the magnitude of the cortical stimulation, suggesting that impaired transient response also entails impaired recurrent interactions within and between BG subnuclei.
Here, we assumed that the GPe to SNr and STN to SNr synapses are static. However, experimental data suggests that synapses between GPe to SNr show shortterm depression (Connelly et al., 2010). Lindahl et al. (2013) has argued that when GPe to SNr synapses show short-term depression, the STN to SNr synapse should also show short-term depression to keep the SNr response small. Lindahl et al. (2013) further showed that short-term depression can have a big effect on the response of the BG when inputs last hundreds of milliseconds. Here, in our model we have only considered very short-lasting stimuli and therefore, short-term depression of synapses might not affect our results. In this work, so far, we have also ignored the effect of NMDA synapses. Such synapses can result in highly nonlinear synaptic integration (Du et al., 2017) and may affect the shape of the transient responses. The role of short-term dynamics and NMDA currents should be investigated in a more detailed model.
Dopamine has multiple effects on neuron's excitability, synaptic strength, and synaptic plasticity (see Table 10). To better understand which one of these are most detrimental for the shape of the transient response, we individually perturbed six of the most crucial parameters (Fig. 6). This analysis revealed that the connection D2-SPN !GPe-TI is the most crucial for the shape of the transient response as it controls both LE and LI zones (Fig. 6). In addition, D1-SPN !SNr connection is expected to be crucial for determining the EI zone. We further corroborated these results by restoring the strength of D2-SPN !GPe-TI connection to their normal level while keeping all other parameters to their low-dopamine levels. This single change was effective in bringing the triphasic response in PD state closer to the one observed in healthy state.
Here, we have essentially characterized the impulse response of the BG network. The impulse response of the SNr does not resembled with the response observed in behavioral tasks (for example, see SNr response in Basso and Wurtz, 2002;Gulley et al. 2002;Wichmann and Kliem, 2004;and GPi response in Schwab et al., 2020). The temporal structure of the BG response means that in behavioral tasks inputs to the BG have a complex temporal structure. Our model can be used to predict the responses in different BG nuclei given a certain output pattern in the SNr or GPi. To illustrate this, we assumed a polyphasic response in activity of SNr. This particular response shape is inspired from the shape of the GPi response Schwab et al. (2020). In our model, such a polyphasic change in the GPi/SNr activity would require that STN and D2-SPNs change their activity before D1-SPNs (see Fig. 11). Moreover, the model gives information about the relative time scale of the transient activity in the D1-SPNs, D2-SPNs, and the STN neurons (see Fig. 11B) in response to the cortical input as shown in Figure 11A. Simultaneous measurement of activity in the striatum, STN and GPi/SNr could verify this prediction of the model. While deriving such input patterns, we only considered the three pathways. Other excitatory pathways like cortical innervation of the GPe could also affect the response. However, more experimental data are needed to infer the potential impact of such inputs.
In network models of b band oscillations when both STN and GPe are included, invariably STN $GPe connections emerge as a key parameter in shaping the oscillations (Holgado et al., 2010;Pavlides et al., 2015). In the full model of BG with both striatum and cortico-BG loop, STN $GPe may not be as important. Indeed Leblois et al. (2006) showed that altered interactions among direct and hyperdirect pathways are sufficient to induce oscillations. However, in the model by Leblois et al. (2006) GPe plays no role in generating oscillationsthis is inconsistent with the experimental data (de la Crompe et al., 2020). In our model, consistent with the recent experimental data (de la Crompe et al., 2020) STN $GPe is not important for generating oscillations. In fact, in our model removal of STN $GPe-TI connections did not affect the oscillations (Fig.  9, blue bars). These observations combined with the experimental data (de la Crompe et al., 2020;Cruz et al., 2011) raise the question which network interactions generate oscillations, if not the STN-GPe loop. We have not explored this question in this work as the question requires a more systematic study. However, we speculate that besides the STN-GPe, the back-projections from GPe to striatum together with recurrent connections within the GPe can form an effective excitatory-inhibitory network necessary for generating oscillations. It is worth noting that previous experimental data (Mallet et al., 2006;de la Crompe et al., 2020;Sharott et al., 2017) and computational models (Kumar et al., 2011;Mirzaei et al., 2017;Bahuguna et al., 2020) provide a strong evidence that strengthening of D2-SPN !GPe-TI connection is also sufficient to induce b band oscillations/synchrony in the ongoing activity state of the BG. Thus, here, we provide a unified explanation of impaired transient response and ongoing activity in PD state. Our results highlight the importance of the GPe in controlling the dynamics and function of the BG.
In our model, emergence of b oscillations does not require synchrony among FSIs unlike some previous studies Damodaran et al. (2015). In our model, FSI synchrony could affect oscillations by modulating the D2-GPe-TI pathway and its effect could both strengthen or weaken the oscillations depending on the parameters. If synchronous firing of FSIs can reduce the firing rate of D2-SPNs, it will weaken the oscillations. If synchronized inhibition from FSIs can synchronize D2-SPNs then it can induce oscillations, regardless of the firing rate (Manferlotti et al., 2021).
Despite its simplicity, our model not only provides network interaction that shapes the properties of transient responses in the BG, but also clearly suggests that recurrent interactions within and between subnuclei of BG are crucial in shaping the transient response. We found that the duration of EE, EI, LI zones (and not LE zone) of the transient responses is largely robust to changes in the BG network interactions while the area/time of the different zones is not. This suggests that in in vivo data, we should find a narrow distribution of the duration of different zones and a wider distribution of the area/time of different zones. Next, our model predicts that by strengthening of cortical inputs, the normal shape of transient response cannot be restored in PD state. This prediction can be tested by either increasing the stimulus strength or by increasing the number of stimulated neurons (e.g., using optogenetic stimulation methods). Finally, the model predicts that by restoring the normal strength of D2-SPN !GPe-TI (or also by reducing the activity of D2-SPN), a near to healthy shape of transient response could be restored even in PD condition.