Effective connectivity of the subthalamic nucleus–globus pallidus network during Parkinsonian oscillations

In Parkinsonism, subthalamic nucleus (STN) neurons and two types of external globus pallidus (GP) neuron inappropriately synchronise their firing in time with slow (∼1 Hz) or beta (13–30 Hz) oscillations in cortex. We recorded the activities of STN, Type-I GP (GP-TI) and Type-A GP (GP-TA) neurons in anaesthetised Parkinsonian rats during such oscillations to constrain a series of computational models that systematically explored the effective connections and physiological parameters underlying neuronal rhythmic firing and phase preferences in vivo. The best candidate model, identified with a genetic algorithm optimising accuracy/complexity measures, faithfully reproduced experimental data and predicted that the effective connections of GP-TI and GP-TA neurons are quantitatively different. Estimated inhibitory connections from striatum were much stronger to GP-TI neurons than to GP-TA neurons, whereas excitatory connections from thalamus were much stronger to GP-TA and STN neurons than to GP-TI neurons. Reciprocal connections between GP-TI and STN neurons were matched in weight, but those between GP-TA and STN neurons were not; only GP-TI neurons sent substantial connections back to STN. Different connection weights between and within the two types of GP neuron were also evident. Adding to connection differences, GP-TA and GP-TI neurons were predicted to have disparate intrinsic physiological properties, reflected in distinct autonomous firing rates. Our results elucidate potential substrates of GP functional dichotomy, and emphasise that rhythmic inputs from striatum, thalamus and cortex are important for setting activity in the STN–GP network during Parkinsonian beta oscillations, suggesting they arise from interactions between most nodes of basal ganglia–thalamocortical circuits.


Introduction
In idiopathic Parkinson's disease (PD) and its animal models, chronic loss of dopamine from basal ganglia circuits profoundly alters the firing rates and patterns of the neurons therein. Excessive synchronisation of oscillatory neuronal activity is an especially prominent functional change accompanying Parkinsonism (Bergman et al. 1998;Hammond et al. 2007). In unmedicated bradykinetic/rigid PD patients, synchronisation of neuronal activity within and between the basal ganglia and cerebral cortex preferentially occurs at beta (13-30 Hz) frequencies (Brown et al. 2001;Hammond et al. 2007). These clinically observed beta oscillations might be truly pathological (Brown, 2006) and have been recapitulated in studies of 6-hydroxydopamine (6-OHDA)-lesioned rats in vivo; excessive beta oscillations emerge throughout the basal ganglia including the output nuclei, subthalamic nucleus (STN), external globus pallidus (GP) and striatum, as well as in the cortex (Sharott et al. 2005;Mallet et al. 2008a,b;Degos et al. 2009;Avila et al. 2010;Moran et al. 2011). The appearance of such rhythms after chronic dopamine loss is associated with hyperactivity and hypoactivity in the STN and GP, respectively, and a shift from mostly uncorrelated ensemble activities to highly synchronised firing across the network (Mallet et al. 2008a,b).
The widespread occurrence of such abnormal neuronal rhythms across the Parkinsonian basal ganglia, together with the myriad interactions between these nuclei, challenges ready definition of the underlying cellular substrates. Experimental and computational modelling studies have highlighted the possibility that the STN-GP network is a 'pacemaker' for synchronised oscillations (Bevan et al. 2002;Gillies et al. 2002;Terman et al. 2002;Holgado et al. 2010). However, past work has taken the mainstream view that GP is populated by a single cell type with homogeneous functional properties, despite a host of contrary evidence (Kita, 2007). Importantly, functional diversity in GP comes to the fore in Parkinsonism, such that two main classes of GP neuron (termed here Type-I and Type-A neurons; see below) can be readily defined in vivo according to the distinct rates and preferred phases of their firing during ongoing slow (ß1 Hz) and beta oscillations (Magill et al. 2001;Zold et al. 2007;Mallet et al. 2008a). This temporal dichotomy persists across extreme brain states, suggesting it is largely subserved by immutable connections (Mallet et al. 2008a). Functional dichotomy within GP aside, neither STN nor GP alone readily generates synchronous pathological oscillations and, although the precise origin of abnormal beta rhythms is undefined, it is likely that they emerge in PD as a consequence of wider network interactions (Bevan et al. 2002;Holgado et al. 2010;Moran et al. 2011). The STN receives massive GABAergic input from GP, and major glutamatergic inputs from the cerebral cortex and intralaminar thalamus (Smith et al. 1998;Bevan et al. 2007). The GP receives most of its GABAergic inputs from striatal and GP neurons, and major glutamatergic input from STN (Smith et al. 1998;Kita, 2007). Which of these key glutamatergic or GABAergic connections of the STN-GP network are mechanistically imperative for generating synchronised oscillations with the hallmark phase differences in GP is unclear, as is the impact of distinct cell types within GP. A resolution would be facilitated by synergising experimental approaches and innovative computational models that embrace GP dichotomy (Jaeger, 2013).
To address these issues and better define the effective connectivity and other parameters underpinning both slow and beta oscillations in the Parkinsonian STN-GP network, we generated a series of firing rate-based computational models that differed in terms of structural and/or 'physiological' architecture. We used extracellular recordings of unit activity in the STN-GP network of 6-OHDA-lesioned rats in vivo to constrain the models and provide an empirical benchmark to quantitatively assess their utility and accuracy. Our models offer several new experimentally tractable predictions, including that the intrinsic properties and connections of GP-TI and GP-TA neurons are quantitatively different.

Ethical approval
All experiments on animals were conducted in accordance with the Animals (Scientific Procedures) Act, 1986 (UK).

Overview of experimental and computational strategy
We explored the effective connectivity underpinning abnormal oscillations in the Parkinsonian STN-GP network in a series of computational models that were successively more complex in terms of either structural and/or physiological architecture (see below). Each computational model was tightly constrained by data from our in vivo electrophysiological recordings of cortical oscillations and neurons in the STN-GP network in a clinically relevant rat model of PD, as well as by biological data in existing publications. As detailed below, our general strategy was to mathematically define several firing rate-based computational models with distinct combinations of fixed and unknown parameters, and then use machine learning techniques to estimate the unknown parameters of most interest, namely those describing the effective connectivity of the neural system. We quantitatively assessed the utility of each model with measures of an accuracy/complexity trade-off, which effectively gauge how well the experimental data are being fitted while penalising model complexity. We use the term 'effective connectivity' to distinguish the approach used here (a model-based characterisation of predicted causal influences) from studies elucidating 'functional connectivity' , which entail descriptive characterisations of the statistical dependence between two time series (see Moran et al. 2011). From the outset here, we follow the widely held view that neither the STN nor GP alone can generate synchronous pathological oscillations and thus, that they emerge in PD as a consequence of wider network interactions (Bevan et al. 2002;Holgado et al. 2010;Moran et al. 2011).

Electrophysiological data collection
Experimental procedures were carried out on adult, male Sprague-Dawley rats (Charles River). Some of the electrophysiological data that were gathered in 6-hydroxydopamine (6-OHDA)-lesioned Parkinsonian rats to constrain our computational models have been published (Mallet et al. 2008a), although we have executed new analyses of these data. The lesioning and recording methods have been documented in detail elsewhere (Mallet et al. 2008a). Briefly, unilateral 6-OHDA lesions were carried out on rats (190-305 g) anaesthetised by inhalation of 2-4% v/v isoflurane (Isoflo, Schering-Plough, Welwyn Garden City, UK) in O 2 . The neurotoxin 6-OHDA (hydrochloride salt; Sigma, Gillingham, UK) was injected (3 μl of a 4 mg ml −1 solution) adjacent to the medial substantia nigra. The extent of the dopamine lesion was assessed 14 or 15 days after unilateral 6-OHDA injection by challenge with apomorphine (0.05 mg kg −1 , S.C.; Sigma). The lesion was deemed successful in animals that made ࣙ80 net contraversive rotations in 20 min. Electrophysiological recordings were performed ipsilateral to successful 6-OHDA lesions in 27 anaesthetised rats (285-470 g at the time of recording) 21-45 days after surgery.
Anaesthesia for electrophysiology was induced with isoflurane (as above), and maintained with urethane (1.3 g kg −1 , I.P.; Sigma), and supplemental doses of ketamine (30 mg kg −1 , I.P.; Pfizer, Sandwich, UK) and xylazine (3 mg kg −1 , I.P.; Bayer, Newbury, UK). Frontal electrocorticograms (ECoGs), electrocardiographic activity and respiration rate were constantly monitored to ensure the animals' wellbeing. Extracellular recordings of the spikes fired by ensembles of single neurons in the STN and/or GP were typically made using multi-electrode arrays ('silicon probes'; NeuroNexus Technologies, Ann Arbor, USA). Because STN neurons are highly active and densely packed, silicon probe contacts within STN often registered multiunit activity (Mallet et al. 2008b). To increase our sample of well-isolated single units, some recordings of STN neurons were made using higher-impedance glass electrodes (Mallet et al. 2008a,b). Unit activity in the STN-GP network was recorded during two cortically defined brain states (see below). After the recording sessions, the animals were killed with ketamine (150 mg kg −1 , I.P.) and perfused via the ascending aorta with 100 ml of 0.01 M PBS (pH 7.4), followed by 300 ml of 0.1% w/v glutaraldehyde and 4% w/v paraformaldehyde in 0.1 M phosphate buffer (pH 7.4), and then by 100 ml of PBS. Brains were sectioned 24-72 h later for verification of recording locations using standard histological procedures (see Mallet et al. 2008a).

Neuronal activity
Unit activity was recorded during cortical slow-wave activity (SWA), which is similar to activity observed during natural sleep, and during episodes of spontaneous 'cortical activation' , which contain patterns of activity that are more analogous to those observed during the awake, behaving state (Steriade, 2000). Both brain states were defined according to the frontal ECoGs recorded simultaneously with unit activity in the STN-GP network (Mallet et al. 2008a). Beta oscillations are only exaggerated in cortico-basal ganglia circuits after a profound and chronic (>4 days) disruption of dopaminergic transmission (Mallet et al. 2008b). Moreover, in anaesthetised J Physiol 592.7 6-OHDA-lesioned animals, exaggerated beta oscillations emerge in these circuits during spontaneously activated brain states (Mallet et al. 2008a,b), thus accurately mimicking the oscillatory activity recorded in awake, unmedicated PD patients (Brown, 2006).
When recorded in epidural ECoGs, the slow (ß1 Hz) oscillation that dominates SWA contains surface-positive 'active' components, which are correlated with the synchronous discharges of cortical projection neurons, and surface-negative 'inactive' components, which are correlated with a widespread and profound reduction in cortical neuron activity (Fig. 1). We first quantified the temporal relationship between the cortical slow oscillation and single-unit activity in the STN and GP as detailed previously (Mallet et al. 2008a). In doing so, we were able to define one type of STN neuron and two main types of GP neuron in 6-OHDA-lesioned rats (Fig. 1). Type-A GP neurons (GP-TA) were statistically defined as those preferentially discharging during the active component of the slow oscillation in ECoGs, whereas Type-I GP neurons (GP-TI) preferentially discharged during the inactive component (Fig. 1). Note that STN neurons preferentially fire during the active component of the slow oscillation (Magill et al. 2001;Walters et al. 2007;Mallet et al. 2008a). Neurons of the same type tend to fire together, with small phase differences, whereas different types of neuron tend not to do so (Mallet et al. 2008a). This diversity in temporal coupling persisted across SWA and activated brain states, suggesting it is strongly governed by 'hard wiring' , that is, the connections we endeavoured to model here (Mallet et al. 2008a;see Fig. 1). From the recorded single-unit activities, we calculated the average maximum, mean and minimum firing rates (during an oscillatory cycle) for each neuronal population for direct comparison with data from our computational models (see below). In our recordings, the peak frequency of the cortical slow oscillation during SWA was 0.92 ± 0.03 Hz (Mean ± SEM), whereas the peak oscillation frequency in the beta band during cortical activation was 20.3 ± 0.3 Hz. These oscillation frequencies were used to guide extrinsic rhythmic inputs to the STN-GP network in our models (see below).

Population-level computational models
A key innovation of our computational modelling approach here is that we consider that GP neurons belong to one of two major populations, that is, either GP-TA or GP-TI, as supported by electrophysiological studies in vivo (Magill et al. 2001;Zold et al. 2007;Mallet et al. 2008a). While the axonal projections of electrophysiologically defined GP-TA and GP-TI neurons have recently been elucidated in part (Mallet et al. 2012), their synaptic inputs are unknown. However, we started the current study knowing that published work strongly suggested that all GABAergic GP neurons in the rat emit local axon collaterals (Sadek et al. 2007) and also innervate STN (Baufreton et al. 2009). Moreover, all GABAergic GP neurons are thought to receive inputs from both striatum and STN (Smith et al. 1998;Kita, 2007). It is also now widely acknowledged that STN receives substantial monosynaptic inputs from ipsilateral frontal neocortex (Fujimoto & Kita, 1993;Nambu et al. 2002;Magill et al. 2004). We thus incorporated these fundamental structural elements as likely effective connections of GP-TA, GP-TI and STN neurons in all our computational models. Another innovation in our approach is the inclusion (in only some of the studied models) of a direct thalamic projection to the STN-GP network. We explore this connection because glutamatergic neurons of the intralaminar thalamus, and particularly those of the parafascicular nucleus (Pfn), innervate both STN and GP (Bevan et al. 1995;Yasukawa et al. 2004), and because intralaminar thalamic projections to the basal ganglia may play a key role in the pathophysiology of PD (Smith et al. 2009). Our full model architecture, which has all connections that are considered in this study, is shown in Fig. 2. This architecture captures all the major intrinsic/extrinsic glutamatergic and GABAergic connections of the STN-GP network, although we note that the extrinsic afferents to these nuclei also include monoaminergic and other inputs with functions that are poorly understood. Although dopaminergic inputs to basal ganglia-thalamocortical circuits are important for controlling the activity patterns therein (Hammond et al. 2007), we did not include the connections of midbrain dopaminergic neurons in our computational models because we sought to investigate only the pathophysiological sequelae of profound and chronic dopamine depletion, a key feature of Parkinsonism. In short, we explored only the effective connectivity of major glutamatergic and GABAergic inputs that directly impinge on STN and/or GP neurons. Corticostriatal and thalamostriatal projections were thus not explicitly modelled with weights, but are included in Figs 2-5 for completeness. For our purposes here, it was not imperative to simulate the STN and GP connections to other basal ganglia structures (i.e. output nuclei and striatum).
Another novel aspect of our modelling strategy is that we sought to reproduce the firing rates (oscillations) of STN and GP neurons that prevail across two brain states, namely SWA and cortical activation. To characterise the firing rates of neuronal populations in STN and GP we use the well-described 'mean firing rate model' (Dayan & Abbott, 2001;Vogels et al. 2005). We have already used population-level mean firing rate models to successfully reproduce some key electrophysiological features of the STN-GP network, including beta oscillations, albeit with a single cell type in GP (Holgado et al. 2010). Mean firing rate models have the additional advantage of providing simple mathematical descriptions, i.e. a first-order differential equation to simulate firing in each neural population. This allowed us to both reproduce our extensive electrophysiological data sets, and to analytically study the model afterwards. All firing rates in this study are expressed as the number of action potentials ('spikes') fired per second (spk s −1 ). The mean firing rate model describes the changes in the firing rate of a neural population by the following equation:  In eqn (1), v is the firing rate of the neural population being modelled,v denotes the rate of change in the firing (i.e. derivative of v), u is the vector of firing rates of the presynaptic neural populations, w is the vector of weights describing strengths of connections from the presynaptic populations, and τ is the firing-rate time constant describing how rapidly the postsynaptic population reacts to its inputs. F(·) is the input-output relationship (or 'activation function') of the neurons in the steady state; this is often termed the frequency-current (f-I) relationship in experimental studies. Using eqn (1) to describe the activity in the STN population and activity in the two distinct GP populations, as mapped onto the full model architecture of Fig. 2, we obtain the following set of differential equations: Note that only some models included connections from the thalamic parafascicular nucleus (Pfn) to the STN-GP network. Italic text pertains to variables used in model equations. Adjacent to each modelled connection are the associated weight (w) and, where appropriate, transmission delay ( t). Corticostriatal and thalamostriatal projections (dashed lines) were not explicitly modelled with weights but are included here for completeness. Each population of GP-TI, GP-TA and STN neurons was described by a differential equation (see eqn (2)), which modelled a variable mean firing rate that tends to a stable constant value (F I (in), F A (in) and F S (in), respectively). The input-output function of each population is diagrammatically shown as a sigmoidal f-I curve, where M n and B n are the maximum and basal firing rates, respectively, of neural population n. B n effectively captures autonomous activity in the absence of inputs. τ n is the membrane time constant. Other parameters: w nm , weight of connection from population n to population m; t nm , delay in transmission from population n to m. Indices G , I , A , S , C , X and P pertain to the GP, GP-TI, GP-TA, STN, cortex, striatum and Pfn neuronal populations, respectively.
In eqn (2), STN is the firing rate of the STN neural population, GP A is the firing rate of the GP-TA neural population, and GP I is the firing rate of the GP-TI population. SṪN,ĠP A andĠP I denote rates of change in the firing of these populations. Ctx, Str and Pfn are the inputs from cortex, striatum and thalamic parafascicular nucleus, respectively. τ S , τ A and τ I are the firing-rate time constants for STN, GP-TA and GP-TI populations, respectively. F S (·), F A (·) and F I (·) are the input-output relationships for STN, GP-TA and GP-TI populations. w nm denotes the weights of the connections from neural population n to neural population m. t nm denotes the transmission delays of the connections from population n to population m. The indices n and m indicate different neural populations and they can be: S for STN, A for GP-TA, I for GP-TI, C for cortex, X for striatum and P for Pfn.
The firing rates of STN, GP-TA and GP-TI neurons are simulated with differential equations (eqn (2)), while cortical, striatal and thalamic firing rates are treated as external inputs to these differential equations. The firing rates of most neurons in cortex, striatum and parafascicular nucleus are oscillatory (at ß1 Hz) during SWA in 6-OHDA-lesioned rats (Mallet et al. 2006(Mallet et al. , 2008aBallion et al. 2008;Parr-Brownlie et al. 2009). During the activated brain state, synchronised beta oscillations (ß20 Hz) are instead observed at the population level (as evinced in local field potentials) in the cortex and striatum of these Parkinsonian rats (Moran et al. 2011). We assume that the firing rates of the neurons in these brain regions oscillate accordingly. Thus, the firing rates of cortical, striatal and thalamic projection neuron populations are described in the computational model by sinusoidal waves, with means, amplitudes and frequencies inferred from the appropriate experimental studies (see below). The equation of this sinusoidal wave for each neural population is: in spikes per second. In eqn (3), f s represents the frequency of oscillations during the specific brain state s (i.e. 1 Hz for SWA and 20 Hz for the activate state), θ n represents the phase of activity in nucleus n with respect to ongoing cortical oscillations, and R sn the mean firing rate of nucleus n during the specific brain state s. The amplitude of oscillations in eqn (3) is equal to α s R sn , thus α s is a scaling constant that describes how large the amplitude of the oscillation is with respect to the mean firing rate during the specific brain state s. The indices s and n indicate, respectively, the brain state and the nucleus related to the parameters α, θ, f and R. The index n can take the values C (cortex), X (striatum) or P (parafascicular nucleus), while s can take the values Act (cortical activation) or Swa (slow-wave activity). Using this notation, the activity of striatum during cortical activation, for instance, becomes: Str(t) = α Act R ActX sin(2πf Act t + θ X )/2 + R ActX (4) in spikes per second. How the firing rates of STN neurons vary as a function of their inputs, as mimicked by somatic current injections (that is, their f-I relationships), has been studied in detail in vitro (Nakanishi et al. 1987;Bevan & Wilson, 1999;Hallworth et al. 2003;Wilson et al. 2004;Barraza et al. 2009). Most STN neurons show sigmoidal f-I relationships ('curves'). Some GP neurons have f-I curves that approximate a sigmoidal function (Deister et al. 2009), while other GP neurons probably do not (Kita & Kitai, 1991;Nambu & Llinas, 1994;Bugaysen et al. 2010). We explore the impact of heterogeneity in the intrinsic properties of GP neurons in some of our computational models (see Results). However, with respect to f-I relationships, theoretical studies have shown that a population of neurons with varied f-I curves can be approximated in a population-level mathematical model by a unit with a sigmoidal input-output relationship (Wilson & Cowan, 1972). Therefore, in all our models, populations of both STN and GP neurons have sigmoidal f-I curves, which we mathematically describe with the function F n (·). To unambiguously determine the F n (·) function, we specify the values of the maximum firing rate (maximum point of F n (·), called M n ), the 'basal' firing rate (F n (0), called B n ), which is equivalent to the mean firing rate when deprived of all synaptic inputs, and the slope of the sigmoid function (determined by S n ). Note that because most STN and GP neurons are able to fire autonomously, i.e. fire at a steady rate in the absence of applied current or synaptic inputs (Surmeier et al. 2005), B n is not zero. Using these three parameters to define the activation function, its equation becomes: In eqn (5), the index n represents the neural population being modelled, and 'in' the input that this population receives.

Model parameters
There are a large number of parameters in the differential equations of the mathematical models. Following the same procedure used by Holgado et al. (2010), we fixed the values of those parameters that could be accurately determined on the basis of published experimental studies, while we estimated as best as possible the values of the parameters that could not be fixed. Table 1 lists the parameter values that were fixed in most models. Note that, to fix model parameters, all values used were derived from experimental studies carried out in rodents (rats), partly because our own electrophysiological data were obtained in rats but also because the basal ganglia of rats is the most thoroughly studied of all species (endowing a large resource of data). In the case of some parameters, it can be debated whether it is better to fix a value to the available experimental data or to estimate it.
To study which parameters should be fixed and which ones should not be fixed, we designed different computational models based on our full architecture (see Fig. 2), and these models differed slightly from each other in terms of which parameters are fixed and estimated. These different models will be described in detail below. Table 2 shows the various parameters that were fixed and estimated for each model. During the description of the different model parameters below, the index notation introduced in the previous subsection will be followed. These indices may be presented with the wild-card letters n or m to indicate neuronal population names, and they can take the values S for STN, A for GP-TA, I for GP-TI, C for cortex, X for striatum and P for parafascicular thalamus. We use an additional index to represent brain state, which is presented with the wild-card letter s and may take the values Act for cortical activation or Swa for slow-wave activity.
The following is a description of the parameters used in our computational models (also see Tables 1 and 2).
Transmission delays ( t nm ). This parameter describes the delay in information transmission between two connected neuronal populations n and m, and is a sum of axonal conduction and synaptic transmission times. More specifically, once a change in firing rate occurs in the presynaptic population (n), this parameter describes the time required for this change to affect the firing rate of the postsynaptic population (m). These parameters are fixed at the same values in all but one of our models, and the values for connections between STN and GP neurons can be readily extracted from electrophysiological studies in vivo (Kita et al. 1983;Kita & Kitai, 1991). However, because the transmission delays between GP neurons themselves are unknown, we approximated the parameter as half of the transmission delay between STN and GP neurons; the spread of the local axon collaterals of GP neurons is less than half the separation between STN and GP in the rat (Sadek et al. 2007).

Firing-rate time constants (τ n ).
Once the input to a group of postsynaptic neurons changes, this parameter measures the time required for the mean firing rate of these postsynaptic neurons to adapt to such a change, i.e. to acquire the new stable firing rate indicated by its sigmoidal input-output relationship. The time constant of a firing rate model is approximately equal to the experimental membrane time constant of the neurons being modelled (Dayan & Abbott, 2001) and therefore we set this parameter to the membrane time constant in all models, as in Holgado et al. (2010). In our models there are three such time constants: for STN (τ S ), GP-TA (τ A ) and GP-TI (τ I ). Their values were estimated from published intracellular recordings of STN and GP neuron membrane potentials (Kita et al. 1983;Nakanishi et al. 1987;Kita & Kitai, 1991), that is, the time needed for these neurons to reach approximately two-thirds of a final 'stable' firing rate in response to somatic injection of a step current. In the absence of definitive empirical evidence, we assume that τ A and τ I are equal.
Synaptic weights (w nm ). These parameters describe the synaptic influence that each presynaptic neuronal population (n) has on each postsynaptic population (m). The set of values of all the synaptic weights represents the effective connectivity of the neuronal network. Because deriving estimates of effective connectivity to and within the STN-GP network was the principal goal of this study, all w nm parameters were estimated in all models.
Oscillation frequencies (f s ). These parameters determine the frequencies of the firing rate oscillations taking place in cortex, striatum and thalamus, which are treated as external inputs to our computational models of the STN-GP network (variables Ctx, Str and Pfn in eqn (2)). We used two such frequency parameters, one when simulating cortical activation (f Act ) and the other when simulating SWA (f Swa ), with values of 20.0 Hz and 1.0 Hz, respectively, as approximated from our electrophysiological recordings (see above).
Relative oscillation amplitudes (α s ). These parameters determine the amplitude of the oscillations taking place in cortex, striatum and thalamus, as a proportion of the mean firing rate of the neuronal population in each of these brain regions. There are two amplitude parameters, one used when simulating cortical activation (α Act ) and other one used when simulating SWA (α Swa ). In all but one of the models, α s parameters were estimated; they were fixed to a value of 1.0 in the exceptional case (see Table 2).
Mean firing rates (R sn ). These parameters represent the mean firing rates of oscillating neurons in cortex, striatum and thalamus during each brain state (SWA or activation). For consistency with our own experimental data, we fixed most of these parameters using values reported in published studies of unit activity recorded during defined brain states in anaesthetised Parkinsonian rats (see Table 1).

Parameter
Value Source of biological data Transmission t SA and t SI 2.8 ms (Kita & Kitai, 1991) delays ( t) t AS and t IS 1.3 ms (Kita et al. 1983) t AA , t IA , t AI and t II 1 ms Based on the above and the assumed closer proximity between GP neurons Time constants (τ) τ S 10 ms (Kita et al. 1983;Nakanishi et al. 1987) τ A and τ I 15 ms (Kita & Kitai, 1991 The value of each fixed model parameter is given together with the literature source from which it was derived. In the rare case that direct experimental evidence was not available, an assumption was made as described. Not all these parameter values were used in all models. Moreover, in some models, a number of the parameters were instead estimated with the fitting algorithm. Table 2 shows the parameters that were estimated as well as those fixed to the values of this table. Swa, slow-wave activity; Act, cortical activation; M, maximum firing rate; B, basal firing rate; S, slope of the f-I curve; B Gdiff , basal firing rate difference between GP-TA and GP-TI neurons. Indices G , I , A , S , C , X and P pertain to the GP (as a whole), GP-TI, GP-TA, STN, cortex, striatum and Pfn neuronal populations, respectively. Note that B A and B I take into account the relative proportions of GP-TA neurons to GP-TI neurons (ratio of 1:3.5), and were set to be symmetrically distributed around the mean basal firing rate of all GP neurons (B G ), i.e. if one population had a basal firing rate higher than B G then the other population had a basal firing rate lower than B G . spk s −1 , spikes fired per second (a measure of firing rate).
Phase difference with cortex (θ n ). These parameters represent the phase difference of oscillations in cortex, striatum and thalamus with respect to cortex. There are three phase parameters (θ C , θ X and θ P ), one for each of these neural populations. By definition then, θ C is equal to 0. In all but one of the models, θ X was estimated; it was fixed to a value of 0 in the exceptional case (see Table 2). In the absence of experimental data to argue otherwise, θ P was also conservatively fixed to 0. Although the time taken for action potential propagation and synaptic transmission could introduce appreciable phase differences, both experimental data and dynamical system theory suggest that widely distributed neuronal populations can synchronise their activity with no phase difference (Roelfsema et al. 1997;Izhikevich, 2007).
Intrinsic properties or f-I curve parameters (M n , B n and S n ). These unambiguously define the shape of the sigmoidal activation functions F n (·) for STN, GP-TA and GP-TI neurons. To fix parameters of maximum and basal firing rates (M n and B n ) we chose, whenever possible, values detailed in published studies of rat neurons that were recorded using 'near physiological' techniques in vitro. This was particularly important for B n to ensure that our models best captured the autonomous firing rates of STN and GP neurons. Thus, we preferentially used data from cell-attached or perforated-patch recordings of neurons that were performed at ࣙ35°C (Abbott et al. 1997;Hallworth et al. 2003;Baufreton et al. 2009;Deister et al. 2009) to minimise the effects of cell dialysis that often cause gradual changes in activity profiles over time (see Barraza et al. 2009;Bugaysen et al. 2010). For B n , we focused on data recorded in the presence of antagonists of glutamate and GABA receptors. We chose the values of S n on the basis of experimental data (Nakanishi et al. 1987;Kita & Kitai, 1991;Hallworth et al. 2003;Deister et al. 2009), such that the F n (·) function reached the firing rates B n and M n evident in published f-I curves J Physiol 592.7 We studied 8 models, each of which incorporated a unique set of parameters that were either estimated with the fitting algorithm or fixed to the values arising from the experimental findings listed in Table 1. The estimated and fixed parameters for each model are marked as 'EST' and 'Fix', respectively. Parameters not used by a model are marked with dashes. For clarity, the parameters with EST/Fix classifications that differed as compared to those of Model 3-Optimal (the model that best fitted the experimental data) are in bold and italic. For definitions of symbols and other abbreviations, see Table 1.
derived from recordings of neuron responses to somatic current injections (in pA); with appropriate input scaling, the activation functions shown in Figs 3-5 mimic these published f-I curves. Because the intrinsic membrane properties and f-I relationships of physiologically defined GP-TA and GP-TI neurons are unknown (they have not been intracellularly recorded), we were conservative in initially estimating M n and S n to be equal for these different cell types. Accordingly, in some of our models, we assumed that both types of GP neuron shared the same basal firing rate (B G ). In other models, we assumed that GP-TA and GP-TI neurons have different basal firing rates, that is, B A and B I , and the difference between them is described by the parameter B Gdiff . This is supported by experimental evidence that not all GP neurons can autonomously fire in a robust manner (Nakanishi et al. 1987;Nambu & Llinas, 1994;Cooper & Stanford, 2000;Bugaysen et al. 2010;Chuhma et al. 2011;Miguelez et al. 2012), including those recorded in vitro after 6-OHDA lesions (Chan et al. 2011;Miguelez et al. 2012).
Modification of B A and B I to take B Gdiff into account will produce horizontal shifts in the two f-I curves. We did not bias the estimation of B Gdiff with respect to the polarity of any difference; B A could thus be higher than B I , or vice versa. Nevertheless, B A and B I were set such that not only was their average equal to the experimentally observed mean basal firing rate of all GP neurons (B G ), but also that they were symmetrically distributed, i.e. if one population had basal firing rate above B G then the other population had basal firing rate below B G . In ensuring that the average of B A and B I was equal to B G , we also accounted for the relative proportions of GP-TA neurons to GP-TI neurons, which is approximately 1 to 3.5 (Mallet et al. 2008a), by weighting B Gdiff according to cell type (see Table 1).

Candidate models of effective connectivity
To investigate the model assumptions about physiological properties and connectivity patterns that lead to the best explanation of (i.e. closest fit to) our in vivo electrophysiological data, we studied eight different computational models based on the mathematical descriptions presented above. All studied models were described by eqn (2), and each is based on different assumptions that define which of the parameters listed in Table 2 were fixed and which were estimated. All models followed a similar division between fixed and estimated parameters, and each model classified only one or a few parameters differently from those commonly employed by all other models (see Table 2). For clarity, the specifics of each model are described with respect to the model that showed the best accuracy/complexity trade-off when fitting the experimental data; this model is referred to as Model 3-Optimal (see Results). Therefore, the first model we present in detail in the Results is called Model 1-NoTB (see Fig. 3) because it differs from Model 3-Optimal in that it has neither thalamus nor basal firing rate differences in GP (see Table 2).

Estimating effective connectivity
Our principal goal was to derive estimates of the effective connections between neurons of the STN-GP network and their extrinsic inputs. To achieve this, we used a fitting technique to determine the values of all estimated parameters in all models (Table 2), while placing special emphasis on those describing the synaptic weights of the connections (w nm ). This fitting technique was designed to find, for a given model, the combination of parameter values that simultaneously best reproduced two experimental data sets, namely the unit activity recorded in the STN-GP network during SWA and cortical activation. Each data set contained measures that could be used to unambiguously describe the oscillations in firing rate taking place in each modelled population of STN, GP-TA and GP-TI neurons. These activity measures are listed in Table 3, and were derived from the experimental data as follows. First, activity histograms and linear phase histograms of single-unit activity with respect to cortical slow and beta oscillations, respectively, were calculated for each qualifying STN, GP-TA and GP-TI neuron as previously described in detail (Mallet et al. 2008a). Next, average phase histograms were computed for each of the three neuronal populations in each of the two brain states. For each population, we repeatedly concatenated the average phase histogram to form a synthetic sequence of periodic firing rates, and then smoothed these firing rates by low-pass filtering. Thus, all components of the fast Fourier transform with frequencies higher than a threshold value were eliminated; these thresholds were set at 1.75 Hz and 22 Hz when filtering SWA data and cortical activation data, respectively, i.e. just above the dominant frequencies of oscillatory firing. Subsequently, we divided the filtered firing rates into individual oscillation cycles, from which we calculated the average 'firing-rate profile' for each neuronal population in each brain state. These profiles express the average firing rate for each phase value and are convenient for the qualitative visual assessment of a given model (see panels A and B in Figs 3-5). For each neuronal population (STN, GP-TA and GP-TI) and brain state (SWA and cortical activation), we also determined the following statistical activity measures to be fitted with the models: the mean, minimum and maximum of the average firing rate during an oscillation cycle (3 per population per brain state); the average phase difference of peak firing with respect to cortex (1 per population per brain state); and the dominant oscillation frequency (1 per brain state). This gives a total of 26 experimental measurements fitted with each model (see Table 3).
To determine the combination of estimated parameter values that best reproduced or best fitted the experimental data sets, we used a computational search technique that employed an advanced, adaptive 'genetic algorithm' , so-called because such algorithms are inspired by the optimisation of genes during the course of evolution (Goldberg, 1989;Mitchell, 1998;Eiben & Smith, 2007). The details of this fitting method are described in Section A of the Supplemental material (available online). Briefly, for each estimated parameter, this technique explores possible values that lie within the boundaries listed in Table 4, and seeks the combination of these values that best reproduces all experimental data sets simultaneously, that is, STN and GP unit activity recorded during both the slow (ß1 Hz) and beta oscillations (ß20 Hz) prevalent in SWA and cortical activation, respectively. The boundaries of most of the parameter values in Table 4 were iteratively selected to give a sufficiently wide exploratory range around values listed in published experimental studies (also see Table 1), while also taking into consideration that the fitting algorithm will fail when boundaries are set too wide (see Section A of Supplemental material). With respect to effective connections, the lower and upper boundaries for the synaptic weights between STN and GP neurons were set at 0 and 5, respectively, because the firing rates and numbers of these neurons are of the same order of magnitude (Oorschot, 1996;Mallet et al. 2008a). The boundaries of synaptic weights of striatal efferents were set wider (0 and 20) to take into account that, first, the firing rates of striatal projection neurons are an order of magnitude lower than those of GP neurons (Mallet et al. 2006(Mallet et al. , 2008a, and secondly, that the number of striatal neurons (and, presumably, their efferents) is several orders of magnitude larger than the number of GP neurons (Oorschot, 1996). The genetic algorithm we used here is essentially an extension of the one used for the same purpose by Holgado et al. (2010), but with the refinement of additional Multiobjective Optimisation, using Pareto sets, and a new measure of model fit, using fuzzy logic operands J Physiol 592.7 In the genetic algorithm, each partial objective function E i assesses how closely the experimental measurements are being fitted by the model. Here we show all 26 activity measurements (note that a single oscillation frequency is listed three times for a given brain state for completeness), indicating in superscript the index of the partial objective function that assesses each quantity. Experimental data sets used here were obtained after additional analyses of the cortical oscillations and STN-GP unit activity reported in Mallet et al. (2008a,b). Note that phase differences here are given with respect to the modelled input rather than the ECoG waveform. deg, degrees.  Lower and upper boundaries of the parameter values searched by the fitting algorithm are indicated. In searching for the combination of free parameter values that allowed each model to best reproduce the experimental data, the genetic algorithm tested only values within these boundaries. a.u., arbitrary units. (Mitchell, 1997;Eiben & Smith, 2007). This combination of techniques attempts to maximise several measures of fit accuracy simultaneously, where each measure is calculated through a partial objective function (E i ). These fuzzy-Pareto functions compare the previously described experimental data with model simulations produced for a given set of parameters. There are three partial objective functions in total, numbered from 1 to 3, and each one of them monitors a series of experimental measurements as indicated by the superscripts in Table 3. Each one of these functions can take any value from 0.0 to 1.0, where 1.0 indicates that the corresponding combination of parameters perfectly reproduces all experimental data (zero error in all experimental measures). To find the combination of parameter values that best reproduces the experimental data, the algorithm seeks the parameter values with the highest scores in these three partial objective functions.
To confirm that we were using sufficient amounts of experimental data to estimate the value of the unknown parameters, we executed the genetic algorithm 500 times over the same experimental data, and obtained 10,000 potential parameter values, which we refer to as 'fittings' . We report on our studies of the most accurate fittings in the Results section. The algorithm code was implemented in object oriented C++ (Deitel & Deitel, 2007) and executed over a Linux operating system on an IBM BlueCrystal supercomputer with 2500 processors in the Advanced Computing Research Centre, University of Bristol. Once the algorithm was executed 500 times, the final Pareto front obtained by each execution was stored in a file. These Pareto fronts were subsequently included into a global collection of possible parameter sets. From this global collection, we selected the 66 parameter sets that produced the best fits (i.e. the ones giving highest fuzzy value E 1 × E 2 × E 3 , where in the case of Model 3-Optimal all E i values were higher than 0.8). These were designated as the 'best sets' . The mean and standard deviation were then calculated for each estimated parameter across the best sets. The resulting mean of each parameter value was considered to be a correct estimation, as the corresponding standard deviations were small (i.e. generally <5% of the exploration boundary size of the algorithm on that parameter). These means and standard deviations are shown in panel C of Figs 3-5. The simulations with the parameters that gave the best fit were used to produce, after graphing in MatLab (MathWorks), the data in panels A, B and E of Figs 3-5; similar results were observed when using some of the other 66 best sets, suggesting that no local maxima with similar fitting measures existed in the parameter space. The absence of other maxima further indicates that the model is not ambiguous, i.e. it could not equally well explain the experimental data with a substantially different combination of parameters.

Comparing the models
The relative adequacy of the models to describe the experimental data was determined by first calculating the sum of squared errors for each model. We then used the sum of squared errors to derive the finite-size-corrected Akaike information criterion (AIC C ), which additionally penalises model complexity, before finally transforming the AIC C into so-called 'Akaike weights' that estimate the probability of each model being the best with respect to an accuracy/complexity trade-off and all other considered models. The AIC C (Hurvich & Tsai, 1989) is a modification of the traditional Akaike information criterion (Akaike, 1974); it is designed to correct for biases that arise when the number of simulated parameters is in the same order of magnitude as the number of fitted data points. The AIC C thus contains a bias-correction term and measures the sum of squared errors of each computational model with respect to the experimentally derived data, while penalising each model depending on the number of parameters that it contains: In eqn (6), j is an index indicating the model for which AIC C is being calculated, k is the number of parameters of the mathematical model, n is the number of experimental values that the model tries to reproduce (as listed in Table 3, i.e. n = 26), e i are the experimental values that the model tries to reproduce (Table 3), and s ij are the corresponding simulation values for model j. Although AIC C values give a relative measure of which model carries the best accuracy/complexity trade-off, the related Akaike's weights are more intuitive quantities (Wagenmakers & Farrell, 2004). Akaike's weights are derived through a simple transformation of AIC C values, and attempt to approximate the probability of each model being the best among all the considered alternatives (8 models were compared here, with higher Akaike's weights indicating better accuracy/complexity trade-offs). These quantities, denoted AW(j), were calculated from the AIC C as follows: In eqn (7), j again indicates the model for which the Akaike's weight is being calculated, while i and l are indexes ranging through all considered models.

Experimental data sets
We recorded frontal electrocorticograms as well as single-unit activity in the STN-GP network of 6-OHDA-lesioned rats during two brain states, slow-wave activity and cortical activation (Fig. 1). Initially focusing on SWA, we quantified the temporal relationships between the dominant cortical slow (ß1 Hz) oscillation and the spike firing of STN and GP neurons. We could define one type of STN neuron and two main types of GP neuron by their temporal signatures; STN neurons (n = 32) and GP-TA neurons (n = 83) preferentially fired during the active component of the cortical slow oscillation, whereas GP-TI neurons (n = 350) preferentially discharged during the inactive component ( Fig. 1A and C). Importantly, this cell-type-dependent diversity in temporal coupling with cortex persisted across brain states, such that the distinct phase-locked firing of STN and GP neurons with slow oscillations ( Fig. 1A and C) was maintained with the excessive beta oscillations (ß20 Hz) that emerged during cortical activation ( Fig. 1B and D). Thus, during activation, STN neurons (n = 49) and GP-TA neurons (n = 62) tended to fire with relatively small phase differences ('in-phase'), whereas GP-TI neurons (n = 280) fired in 'anti-phase' to both. We derived several measures of neuronal activity from these experimental data sets to test the validity of the model data (see below). In short, the suitability of each computational model was quantitatively assessed by its ability to match these experimental in vivo data.

Model architecture
We studied eight firing-rate-based computational models that each incorporated various features of the full architecture shown in Fig. 2. In this architecture, the STN-GP network is embodied as a reciprocally coupled system; glutamatergic STN neurons project to and exert excitatory postsynaptic effects on GP neurons, whereas GABAergic GP neurons project to and exert inhibitory postsynaptic effects on STN neurons (and other GP

. Data fits and effective connectivity obtained with Model 1-NoTB
This model excluded thalamus and also assumed that all GP neurons have identical input-output functions. A, left, average firing-rate profiles of STN, GP-TA and GP-TI neurons within a single cycle of the slow oscillation (peak at 1 Hz) during slow-wave activity. Grey bins indicate experimental data and black lines indicate model data. All data sets were individually normalised according to the minimum average firing rates (0%) and the maximum average firing rates (100%). Right, comparison of average firing rates, network oscillation frequencies and phase relationships in experimental data (crosses) with the best data fits from this model (circles). Note the models attempted to fit these statistical activity measures rather than the firing-rate profiles per se. B, firing-rate profiles of neurons (left) and their activity measures (right) during the beta oscillations (peak at 20 Hz) that prevail in cortical activation. This model reproduced the activity of STN neurons, but not that of either GP-TI or GP-TA neurons, with reasonable accuracy for both SWA (A) and activation (B). C, the estimated and fixed parameter values in this model. Data are means (±SEM) of the estimated parameter values in the 66 best fits derived from multiple executions of the genetic algorithm. Parameters fixed to a value of zero are represented as stars. D, the input-output functions (f-I curves) that were fixed for each neuron population. Input is scaled to mimic experimental f-I curves of neuron responses to somatic current injections (in pA). Note that curves for GP-TI and GP-TA neurons are the same in this model. E, graphical representation of the effective connectivity estimated with neurons). The STN (only) receives additional excitatory inputs from cortex, and the GP (only) receives additional inhibitory inputs from striatum. Both STN and GP can also receive excitatory inputs from the thalamus, embodied by connections originating in Pfn. In all our models, GP neurons are explicitly divided into two major populations (GP-TA and GP-TI). We explored only the effective connectivity of those glutamatergic and GABAergic inputs that directly impinge on STN and/or GP neurons (Fig. 2). All our computational models contained a group of parameters whose values were fixed on the basis of published experimental data (see Methods and Table 1), and another group of parameters whose values were estimated. The models differed from one another in terms of the parameters that were fixed and estimated (Table 2).

Model data fits and effective connectivity
We used our recordings of ECoGs and unit activities made during both SWA and activation to calculate a set of statistical measurements for STN, GP-TA and GP-TI neurons, according to brain state. These statistical measures summarise the suprathreshold activity patterns of each of the modelled neuronal populations, and the computational models were therefore fitted to them. The activity measures we used were: the mean, minimum and maximum of the average firing rate during an oscillation cycle; the average phase difference of peak firing with respect to cortex; and the dominant oscillation frequency (see Table 3). We also generated 'firing-rate profiles' that embodied the average firing rate during each phase of the dominant ongoing network oscillation (peak frequencies of either ß1 Hz or ß20 Hz; see Methods). For each model, we used a genetic algorithm (see Section A of Supplemental material) to find the combination of parameter values that allowed the model to best fit or reproduce the activity measures arising from these experimental data. Among the estimated parameters, the 'synaptic weights' (w nm ) are of special interest because they define the effective connectivity of the STN-GP network as predicted by each model.
The first, and simplest, of the three models presented in detail here is Model 1-NoTB (Fig. 3). It includes one STN population and two GP populations, which receive extrinsic inputs from cortex and striatum, respectively, and is based on two key assumptions: (1) that the parafascicular nucleus of the thalamus does not significantly contribute to the patterned firing rate fluctuations of neurons in the STN-GP network during slow and beta oscillations, and (2) that GP-TA and GP-TI neurons have the same intrinsic properties, i.e. the same input-output (f-I) curves (Fig. 3D). To implement the first assumption, the synaptic weights coming from Pfn to the STN-GP network (w PS , w PA and w PI ) were ignored by setting them to 0. Likewise, the parameters describing the firing rates of Pfn during SWA (R SwaP ; see Table 1) and cortical activation (R ActP ) are irrelevant because they have no effect on the behaviour of this model. The second assumption constrains the basal firing rates of both GP populations to be the same, and therefore B G from Table 1 is used in the model instead of B A and B I (B Gdiff was thus set at 0). In summary, this first model had 29 fixed parameters and 15 estimated parameters, as indicated in Table 2. Model 1-NoTB provided reasonably good fits of the firing-rate profiles and activity measures of STN neurons during SWA (Fig. 3A) and cortical activation (Fig. 3B). However, this first model struggled to accurately reproduce the activity measures of both GP-TA and GP-TI neurons, evident as substantial deviations of the model data from the experimental data ( Fig. 3A and B). Despite the low accuracy of the best fits, the genetic algorithm gave stable values for the estimated parameters (Fig. 3C). When synaptic weight parameters were graphically represented, the resulting map of the effective connectivity of the STN-GP network (Fig. 3E) suggested a number of salient features, including: connections from striatum are much stronger to GP-TI neurons than to GP-TA neurons (ratio of striatal connection weights of ß20:1); connections within and between the two populations of GP neuron are not equivalent; and connections from both GP populations to the STN are relatively weak. Model 1-NoTB also suggested that STN connections are stronger to GP-TA neurons than to GP-TI neurons.
Because the ability of Model 1-NoTB to accurately reproduce the experimental data was not entirely satisfactory, we designed and tested a second model called Model 2-NoB, which explored in detail whether the inclusion of thalamic connections to STN and GP improved the ability of the model to reproduce the experimental data (Fig. 4). This second model, in common with the previous one, is based on the assumption that the intrinsic properties of GP-TA and GP-TI neurons are the same (i.e. no basal firing rate differences in GP; Fig. 4D). However, in contrast to the previous model,     Model 2-NoB assumed that the thalamus contributes to the experimentally observed firing rates of neurons in the STN-GP network. Because Pfn was included in this model, the synaptic weights coming from the thalamus were not equal to 0, and they were included in the list of estimated parameters that the fitting algorithm had to calculate (see Table 2). Likewise, the parameter describing the firing rates of Pfn neurons during cortical activation (R ActP ) became part of the list of estimated parameters, as its value could not be derived from existing literature but should nevertheless affect the behaviour of the model. Compared to the first model, Model 2-NoB provided relatively poor fits of the activity measures of STN neurons during SWA (Fig. 4A) and cortical activation (Fig. 4B). Notably though, Model 2-NoB was able to more accurately reproduce the activity measures of GP-TA neurons and GP-TI neurons ( Fig. 4A and B). The first and second models were similar in other respects, with the latter recapitulating many of the previously estimated parameter values (Fig. 4C) and salient features of effective connectivity (Fig. 4E). Moreover, the inclusion of thalamic connections served to reinforce the finding of Model 1-NoTB that STN connections to GP-TA neurons are of a higher weight than those to GP-TI neurons (ratio of STN connection weights of ß4:1; Figs 3E and 4E). Perhaps most striking, Model 2-NoB suggested that Pfn connections to STN neurons and GP-TA neurons are much stronger than Pfn connections to GP-TI neurons (ratio of thalamic connection weights of ß8:1; Fig. 4E).
Although the accuracy of Model 2-NoB in reproducing the experimental data was better than that of Model 1-NoTB, neither could adequately reproduce all the experimental measures of activity in all three neuronal populations in the STN-GP network. Because of this shortcoming, we designed and tested a third model, called Model 3-Optimal, which explicitly accounted for the possibility that the two types of GP neuron have intrinsic physiological differences (Fig. 5). This seemed a reasonable way forward given that: (1) not all GP neurons can fire autonomously in a robust manner, i.e. at 20 spk s −1 , which was the basal firing rate assumed for all GP neurons in the absence of inputs (fixed parameter B G ) in Models 1 and 2 (Nakanishi et al. 1987;Nambu & Llinas, 1994;Cooper & Stanford, 2000;Bugaysen et al. 2010;Chuhma et al. 2011;Miguelez et al. 2012); and (2) GP-TA and GP-TI neurons have significantly different firing rates in vivo (Mallet et al. 2008a), which might be underpinned by differences in their intrinsic membrane properties. Thus, Model 3-Optimal differs from Model 2-NoB in that the basal firing rates of GP-TA and GP-TI neurons were assumed to be different from each other. This difference was introduced into the mathematical definition of the model by using the parameters B A and B I which, through the estimated parameter B Gdiff , give a different basal firing rate to each GP population (Table 1). Like in Model 2, however, the Pfn connections were included, which meant that the parameters w PS , w PA , w PI and R ActP were retained in the list of estimated parameters (Table 2). Compared to the previous two models, Model 3-Optimal provided much better fits of the firing-rate profiles and activity measures of STN, GP-TA and GP-TI neurons; it almost perfectly reproduced the experimental data recorded during either SWA (Fig. 5A) or cortical activation (Fig. 5B). These better fits were associated with clear horizontal shifts in the f-I curves of GP-TA and GP-TI neurons (with respect to each other), such that B A was 6 spk s −1 and B I was 24 spk s −1 (Fig. 5D). Thus, the genetic algorithm estimated that GP-TA neurons have substantially lower autonomous firing rates than GP-TI neurons. Model 3-Optimal further reiterated many of the salient features of effective connectivity that had emerged in the previous two models. For example, the connection from striatum to GP-TI neurons was ß20 times stronger than that to GP-TA neurons, whereas the connection from Pfn to GP-TA neurons (and to STN neurons) was more than 10 times stronger than that to GP-TI neurons ( Fig. 5C  and E). A further contrast between GP-TA and GP-TI neurons was evident at the level of their inputs from STN, with the subthalamic connection to the GP-TA population being ß4 times stronger. Different connection weights between and within the two populations of GP neuron were also a consistent feature: the connection from GP-TI neurons to GP-TA neurons was approximately twice the weight of the reciprocal connection; the local connections between GP-TA neurons were more modest; and the weight of the local connections between GP-TI neurons was almost negligible (Fig. 5C and E). Model 3-Optimal also revealed two more notable features of effective connectivity. First, although the connection from GP-TI neurons to STN was relatively modest, the connection from GP-TA neurons to STN was tiny. Thus, while the reciprocal connections between GP-TI and STN and GP-TI curves. E, scaled graphical representation of the effective connectivity estimated with this model, with thicker lines and larger terminations indicating relatively higher connection weights. In recapitulating and refining some salient features of effective connectivity of Models 1 and 2, this model suggested that the connection from striatum was much stronger to GP-TI neurons than that to GP-TA neurons, whereas the connection from thalamic Pfn to GP-TA neurons (and to STN neurons) was much stronger than that to GP-TI neurons. Moreover, the reciprocal connections between GP-TI and STN neurons were well matched in weight, but those between GP-TA and STN neurons were not. Different connection weights between and within the two populations of GP neuron were also evident. This figure follows the conventions of Figs 2-4. neurons were well matched in weight, those between GP-TA and STN neurons were far from being equivalent. Second, the connection weight from Pfn to STN became larger than the connection weight from cortex to STN (Fig. 5E).

Further model comparison
The findings above indicate that, among the three models described in detail, Model 3-Optimal most accurately reproduced the electrophysiological data sets. We next compared the models systematically, and assessed which showed the best accuracy/complexity trade-off for our experimental data set. The fitting algorithm we used here employs a specialised measure of fitting accuracy, the so-called 'fuzzy-Pareto estimation functions' (see Methods). This measure evaluates how accurately each combination of estimated parameter values reproduces the experimental data, and it is designed on the basis of heuristic considerations, which maximise the capacity of the algorithm to find good fits. For this reason, fuzzy-Pareto functions are not entirely suitable for formal comparisons of model accuracy. We therefore used more standard error measures for model comparison: the sum of squared errors, which directly measures the accuracy of each fit; and the so-called Akaike's weights, which are derived from the sum of squared errors and additionally take into account the mathematical complexity of each model (see Methods). Akaike's weights are thus a useful shown together with the corresponding Akaike's weight (B). A smaller squared error indicates that a given model can better describe or fit the experimental data, while the Akaike's weight approximates the probability of a given model being the one with best accuracy/complexity trade-off as compared to the other models tested. Note that Model 3-Optimal provides the best fit, even though it is not the most complex model tested. metric for selecting models with the best trade-off between accuracy and complexity, with higher values indicating better matches between a model and the constraining (experimental) data as compared to the other models being considered. As expected, Model 3-Optimal had a lower summed squared error (Fig. 6A) and a higher Akaike's weight (Fig. 6B) than Model 2-NoB and Model 1-NoTB. Furthermore, Model 2 had a substantially lower squared error and a higher Akaike's weight than Model 1 (Fig. 6), indicating that the Pfn connections in the former were particularly important for more accurate fits. To explore this further, we constructed a fourth model (Model 4-NoT) that differed from Model 1 only in that it assumed that the two types of GP neuron have intrinsic differences (see Table 2). Model 4-NoT did not perform better than Model 1-NoTB (Fig. 6). Together with the comparison of Model 2 with Model 3, this indicates that the freedom to estimate parameters B A , B I and B Gdiff only improved model fits when thalamic inputs were already included. However, Models 1-4 clearly suggest that intrinsic differences between the two GP populations, and Pfn connections, both contribute to the experimentally observed activity in the STN-GP network. Equally important, these comparisons also provide a first indication that model accuracy did not improve by simply adding complexity (i.e. extra estimated parameters), but rather, that fits were highly dependent on the connections present and the precise combination of estimated and fixed parameters.
To explore further whether the assumptions (i.e. combinations of estimated and fixed parameters) of Model 3-Optimal were necessary and sufficient to explain the experimental data, we built on Model 3 to characterise four additional models (Table 2) and compared their accuracies as above. The first of these additional models, Model 5-Noθ, was only different from Model 3 in that it assumed that the phase difference between the oscillations taking place in cortex and striatum (see Table 4) is not necessary to explain the experimental data, and therefore sets parameter θ X to 0. Model 6-Noα was similar to Model 3 but assumed that the firing rate oscillations taking place in cortex, striatum and Pfn have the same amplitude during SWA and cortical activation, which is equal to the mean firing rate in each structure. This constrained the parameters α Swa and α Act to be equal to 1 (Table 4). Model 7-NoR assumed that the firing rates of striatal and Pfn neuronal populations do not significantly change upon transition from SWA to the activated state. This constrained R ActX and R ActP to be equal to R SwaX and R SwaP , respectively (see Table 1). Finally, we constructed an eighth model that had many more estimated parameters than Model 3-Optimal in order to re-examine whether adding complexity brought further improvements in the fitting accuracy. This Model 8-Do had the freedom to explore different transmission delays ( t nm ) from those J Physiol 592.7 fixed in all other models on the basis of published electrophysiological data (see Table 1), i.e. it allowed the genetic algorithm to find the values of transmission delays that best fitted the experimental data. Again, Model 3-Optimal had a lower summed error and a higher Akaike's weight as compared to Models 5, 6 or 7 (Fig. 6). In summary, these results suggest that a small phase difference between cortex and striatum is of some importance (Fig. 5C), that the firing rate oscillations in cortex, striatum and Pfn have different amplitudes during SWA and cortical activation (Fig. 5C); and that the firing rates of striatal and Pfn neuronal populations change upon transition from SWA to the activated state (Fig. 5C). Indeed, our models predict that, as compared to activity during SWA (Table 1), their firing rates approximately double (Fig. 5C). Although Model 3-Optimal and Model 8-Do had similarly low summed errors, the Akaike's weight of Model 3 was much higher than that of Model 8 (Fig. 6). The transmission delays estimated by Model 8 were ß2 ms between STN and GP (both directions), and 1 ms between GP neurons, which are in good agreement with the values fixed from the literature (Table 1) and used in all other models. This also implies that the additional complexity in Model 8 is not statistically justified, and again highlights that more complex models are not necessarily more accurate.

Elements required for experimentally observed phase relationships
During Parkinsonian network oscillations, STN and GP-TA neurons tend to oscillate in-phase with each other, but anti-phase to GP-TI neurons (Fig. 1). Our best candidate computational model accurately reproduced the experimental data, including these persistent phase relationships. However, although this model is optimal with respect to explaining the totality of our experimental data, it does not reveal which of its elements are responsible for the in-phase and anti-phase firing observed across the STN-GP network. In a final investigation, we thus sought to analytically identify the key elements underlying these oscillatory phase relationships. In essence, we mathematically derived the set of minimal conditions that the model (of six neuronal populations; see Fig. 5E) needed to satisfy in order to reproduce the inverted phase of GP-TI neuron oscillations. We provide a detailed explanation of our analysis in Section B of the Supplemental material. Briefly, as an initial step, we were able to eliminate from the main differential equation of our model (see eqn (2)) the elements that are unlikely to play major roles in the phase inversion of GP-TI neurons. These elements were the sigmoidal activation functions and time constants of the modelled STN, GP-TA and GP-TI neurons, as well as the transmission delays between them; this elimination step reduced eqn (2) to a linear equation. Then, using complex numbers to represent the oscillations present in each neuronal population (i.e. 'phasors' algebra), we were able to derive the minimal conditions in this linear equation that inverted the phase of GP-TI neuron oscillations. We observed, first, that these conditions were critically dependent on the connection weights of eqn (2) and, secondly, that they clearly divided the neuronal populations into two groups. Each group was characterised by the predominant phase at which the member population(s) oscillated; the first group was GP-TI neurons, and the second group consisted of all other neuronal populations. Figure 7 graphically represents the derived conditions, and shows how excitatory connections were prominent between member populations of the second group, all of which oscillated in-phase with each other. In contrast, inhibitory connections were prominent (and excitatory connections were low) between this second group and the GP-TI neurons that oscillated in anti-phase to them (Fig. 7). In summary, the excitatory connections tend to maintain in-phase firing, while the inhibitory connections tend to maintain anti-phase firing. This seems intuitive given that neurons receiving strong excitatory inputs will tend to fire most when the firing of these afferents is maximal, whereas neurons receiving relatively strong inhibitory inputs will tend to fire least during the period of the oscillation at which the firing rate of the afferent neurons is maximal. It is important to note, however, that this derivation with phasors algebra shows it is these connection features, rather than more exotic dynamics, that are responsible for the phase relationships observed in the STN-GP network in vivo. Thus, strong excitatory connections along an axis of cortex, thalamus, STN and GP-TA would promote firing therein at ß0 deg, whereas strong inhibitory inputs to and from GP-TI neurons would promote their firing at ß180 deg to that of the other neuronal populations (Fig. 7). This is especially pertinent for cases in which the transmission delays between neuronal populations are much shorter than the oscillation period, as would be the case for the slow oscillations that dominate SWA. However, as the oscillation period becomes shorter and approaches the transmission delays, additional phase shifts between nuclei may appear. This helps explain why, during beta oscillations, the phases of STN, GP-TA and GP-TI neurons are not as close to the values 0 and 180 deg.

Discussion
We generated a series of computational models to explore the effective connections and physiological parameters underlying the distinct firing rates and phases of STN neurons and two populations of GP neurons that were experimentally observed during excessively synchronised oscillations in the Parkinsonian brain in vivo. Our best candidate model (Model 3-Optimal) was able to reproduce the experimental data across two brain states with great accuracy, and predicted how the afferent and efferent connections of each type of neuron in the STN-GP network are quantitatively different.

Limitations of the computational models
We studied several firing rate-based computational models of the STN-GP network, each of which contained dozens of fixed and estimated parameter values. Because we used a genetic algorithm to find the combination of estimated values that endowed the best model fits to the experimental data, it could be argued that our approach is relatively slow and computationally extensive. However, it gave us the opportunity to incorporate a large amount of biological data to fix or guide model parameter values. Dynamic causal modelling (Friston et al. 2003) of effective connections in Parkinsonian brain networks has proven utility (Moran et al. 2011), and would probably be quicker to implement, but this approach cannot yet be constrained by key experimental measurements used here (i.e. unit activity data). Models simulating large networks of spiking neurons would provide an extra step towards biological reality, but would be a challenge to implement for our purposes here, not least because a substantial number of extra parameters of unknown empirical value would arise. Our firing rate-based models thus provided a reasonable and tractable 'middle ground' . In them, we assumed that the populations of glutamatergic neurons in STN, cortex and thalamus exert only excitatory effects on their direct targets, and that the populations of GABAergic GP-TA, GP-TI and striatal neurons exert inhibitory effects. Studies in vitro show that synaptic dynamics in the STN-GP network extend beyond classical excitation and inhibition, and inputs may cause phase locking, 'rebound bursting' , and/or the 'shunting' or enhancement of other inputs (Bevan et al. 2007;Atherton et al. 2010). Our models were not designed to explore such dynamics, but rather, because our strategy is innovative, we chose to constrain them using a physiologically well-characterised data set comprising hundreds of single units recorded extracellularly in vivo (Mallet et al. 2008a,b). Although imperfect, Model 3-Optimal still generated excellent fits to these experimental data. That said, modelled connection weights do not reveal the precise biological substrates, e.g. a strong connection weight could be subserved by a few powerful synapses and/or by many weak synapses. We can nevertheless interpret our estimates of effective connectivity in light of the known structural and physiological properties of STN and GP neurons and their inputs.

Predicted effective connectivity of the Parkinsonian STN-GP network
In order for Model 3-Optimal to accurately reproduce the patterns of unit activity recorded in the STN-GP network, it incorporated certain assumptions about the effective connectivity and electrophysiological properties of the constituent neurons. By using the term effective connectivity, we distinguish our approach (a model-based characterisation of predicted causal influences) from studies of functional connectivity. In accepting that the computational model is a useful tool for reproducing STN-GP activity, we can conclude that the model assumptions are not only necessary to reproduce STN-GP network dynamics in silico, but are also important to explain the experimental data. Some of these assumptions and their corresponding predictions are unintuitive, and their necessity would be difficult to infer without model simulations. We generated four sets of novel predictions about the effective connectivity of the STN-GP network during the excessively synchronised oscillations that arise in idiopathic PD and its animal models.
First, our modelling data provide a compelling theoretical basis to assume that thalamic neurons innervating STN and GP, and nominally those in the parafascicular nucleus in the caudal intralaminar group, play important roles in setting activity levels in the STN-GP network during Parkinsonian oscillations. Indeed, adding Pfn connections substantially increased the ability of the models to accurately reproduce our experimental data (compare Models 1-4). Importantly, our best candidate model also suggested that Pfn J Physiol 592.7 innervation of the STN-GP network is highly selective in strength. Thus, while the weights of Pfn connections to STN and GP-TA neurons were similarly large, the Pfn connection to GP-TI neurons was tiny. Anatomical data suggest that intralaminar thalamic projections to STN are modest, at least when compared to thalamostriatal projections (Smith et al. 2004). However, our models indicate that the weight of this connection to STN is large (even greater than that of the corticosubthalamic projection), suggesting a paradoxically strong impact. Electrophysiological analyses combined with selective activation of intralaminar thalamic efferents, e.g. through optogenetic approaches (Ellender et al. 2012), will be required to test this prediction. For the GP, a selective or specific innervation of GP-TA neurons, the smaller of the two neuron populations, could be one substrate by which a relatively small thalamopallidal input (Smith et al. 2004) could still powerfully impact on local activity and, thence, the wider STN-GP network. With respect to STN-GP activity during slow oscillations, our predictions are supported by the demonstration in 6-OHDA-lesioned rats that Pfn neurons fire in time with STN neurons (Parr-Brownlie et al. 2009) and thus, GP-TA neurons. Whether and how Pfn neurons, or any other type of thalamic neuron, fire individually or collectively during Parkinsonian beta oscillations is unknown (see below). However, our model predictions help explain why surgical interventions targeted to the intralaminar thalamus have some therapeutic benefits in Parkinsonism (Smith et al. 2009;Jouve et al. 2010), as is the case for those targeted to STN or GP (Limousin et al. 1998;Vitek et al. 2004).
Secondly, all of our models identified a mismatch in the weight of striatal connections to GP-TI and GP-TA neurons. Ultimately, the connection to GP-TI neurons was an order of magnitude higher than the connection to GP-TA neurons. This might partly explain why GP-TI neurons (but not GP-TA neurons) fire anti-phase to striatopallidal neurons during SWA (Mallet et al. 2006(Mallet et al. , 2008a. Our model data concurred that striatopallidal neurons fire infrequently (Mallet et al. 2006) but their huge numbers, convergent projections and synchronised firing could still result in powerful postsynaptic effects. The key prediction then is that rhythmic striatal outputs are better suited to directly influence the activities of GP-TI neurons than those of GP-TA neurons during network oscillations. This prediction has recently been validated for slow (ß1 Hz) oscillations by pharmacological manipulations of striatum in 6-OHDA-lesioned rats (Zold et al. 2012). Moreover, an optogenetics study in vitro identified two classes of GP neuron that can be readily distinguished by the strengths of their striatal inputs, amongst other electrophysiological differences (Chuhma et al. 2011). It is thus tempting to speculate that these cell classes receiving strong and weak striatal inputs in vitro are GP-TI and GP-TA neurons, respectively. Thirdly, our models predicted that the reciprocal connections between GP-TI and STN neurons are modest and well matched in weight, but those between GP-TA and STN neurons are not balanced. The connection from GP-TI neurons to STN is stronger than that from GP-TA neurons, whereas the connection from STN is stronger to GP-TA neurons. Thus, while connections from STN neurons might directly and powerfully impact on GP-TA neuron firing, the reverse might not be true. An entirely different modelling approach has also suggested that GP-TI neurons are more important than GP-TA neurons for sculpting STN neuron firing during beta oscillations (Cruz et al. 2011). This surprising prediction appears valid because it has very recently been shown that only GP-TI neurons substantially innervate STN (Mallet et al. 2012). Importantly then, our current models corroborate these previous theoretical and experimental findings, thus providing further 'face validity' . There is insufficient experimental evidence to validate the prediction that connections from STN are relatively stronger to GP-TA neurons, although this idea is supported by the fact that GP-TA neurons (but not GP-TI neurons) tend to fire in-phase with STN neurons. Yet, this could also result from disparate striatal innervation and/or a common input such as that from Pfn (see above). Recording the responses of GP-TA and GP-TI neurons to selective stimulation of STN efferents would help resolve this.
Fourthly, our models pointed to substantial differences in the strengths of the local connections within and between the two populations of GP neurons. Of these local connections, that from GP-TI to GP-TA was relatively strong, while the reciprocal was modest. The connection between GP-TA neurons was similarly modest, while that between GP-TI neurons was negligible. Recent anatomical data establish the precedent that GP-TI and GP-TA neurons innervate other GP-TI neurons, though quantification is lacking (Mallet et al. 2012). Whether GP-TI and GP-TA neurons innervate other GP-TA neurons is unknown. Nevertheless, our models make the novel prediction that there is considerable cell-type selectivity in the innervation patterns and postsynaptic impact of local axon collaterals of GP neurons. Moreover, the current models and past theoretical studies (Cruz et al. 2011) concur that inhibitory interactions between GP-TI and GP-TA neurons could promote anti-phase firing across the two populations (also see Miguelez et al. 2012).

Additional predictions
Further predictions arise from estimated parameter values other than those describing connection weights. As discussed above, our models suggested that oscillatory outputs from striatum (striatopallidal neurons) and thalamus (Pfn neurons) are important for the slow and beta oscillations recorded in the STN-GP network in Parkinsonism. How striatopallidal and Pfn neurons fire during excessive beta oscillations is unknown, but our models predict not only that they will rhythmically synchronise at beta frequencies but also that their firing rates will approximately double with respect to those during SWA. Because excessive beta oscillations emerge in striatal field potentials after 6-OHDA lesions (Moran et al. 2011), it is likely that the firing of some striatal projection neurons will be phase-locked accordingly. Moreover, the fact that neurons in basal ganglia output nuclei, some of which will innervate thalamus, are engaged by excessive beta oscillations (Brown et al. 2001;Avila et al. 2010), together with anecdotal thalamic recordings in unmedicated PD patients (Kempf et al. 2009), lend support to the prediction of recruitment of Pfn activity to abnormal beta oscillations. It remains to be determined whether striatal and thalamic neuron firing rates increase as predicted (but see Zold et al. 2012).
Model 3-Optimal also predicted that, in the absence of synaptic inputs, the rate of basal activity of GP-TA neurons is markedly lower than that of GP-TI neurons. Whether the autonomous firing of GP-TA and GP-TI neurons differs is unknown, partly because their definition currently requires electrophysiological recordings in vivo rather than in vitro. However, this prediction is corroborated by reports that not all GP neurons can produce robust (>10 spk s −1 ) autonomous activity (Nakanishi et al. 1987;Nambu & Llinas, 1994;Cooper & Stanford, 2000;Bugaysen et al. 2010;Chuhma et al. 2011;Miguelez et al. 2012). The idea that GP-TA neurons are intrinsically less active would also fit well with their predicted connectivity; these cells might not need to generate robust activity themselves because, as compared to GP-TI neurons, they receive stronger excitatory connections from STN and thalamus, and weaker inhibitory connections from striatum. In line with this, in vitro work has revealed a class of unusually slow-firing GP neurons that receive only weak striatal inputs (Chuhma et al. 2011). Because the autonomous firing of some (but not all) GP neurons might be compromised after 6-OHDA lesions (Chan et al. 2011;Miguelez et al. 2012), the basal firing differences predicted here could reflect pathological changes in the intrinsic membrane properties of GP neurons in Parkinsonism. Such pathological alterations might be cell-type selective, with the excitability/activity of GP-TA and GP-TI neurons being differentially regulated by dopaminergic transmission or its disruption (also see Billings & Marshall, 2004).

Further functional considerations
Here, we used computational modelling to gain new insights into the factors underlying the complex activity dynamics of the Parkinsonian STN-GP network in vivo. Our models are the first to explicitly divide the GP into two neuronal populations, and to reproduce the firing rates/phases of STN-GP neurons that prevail in two brain states and the distinct oscillations that define them. We formally demonstrate using mathematics tools that the connections of GP-TA and GP-TI neurons are quantitatively different. GP-TA neurons receive stronger excitatory connections from STN and thalamus, whereas GP-TI neurons receive stronger inhibitory connections from striatum. The inhibitory connections from each type of GP neuron to STN are also of different weights, as are their local connections. Together with intrinsic physiological differences in GP-TA and GP-TI neurons, these distinct connectivity patterns are predicted to be of causal significance for the observed STN-GP network activity, particularly for the in-phase vs. anti-phase oscillatory firing relationships. Our data reiterate that the GP should not be considered a functionally homogeneous structure in circuit-level descriptions (Mallet et al. 2012), and argue for a cell-type-specific diversity of connections being incorporated into future schemes of the functional organisation of the basal ganglia.
We have previously established some conditions that are theoretically necessary for the STN-GP network to independently generate beta oscillations (Holgado et al. 2010;Pavlides et al. 2012). The current models do not address the issue of whether the connections between STN neurons and the two types of GP neurons could subserve a central pattern generator function in vivo (this would not apply for slow oscillations anyway because they originate in cortex; see Chauvette et al. 2010). Instead, our new results emphasise the importance of rhythmic inputs to the STN-GP network, and particularly those arising in cortex, striatum and thalamus, for setting the activity of STN, GP-TA and GP-TI neurons during Parkinsonian oscillations. Excessive beta oscillations are thus likely to arise from the interactions of most nodes of basal ganglia-thalamocortical loop circuits.
Although we focused our investigations on the excessive beta oscillations that arise after severe and chronic dopamine depletion (Mallet et al. 2008b), it is important to note that lower levels of beta oscillations can be observed in basal ganglia-thalamocortical circuits when dopamine is intact (Mallet et al. 2008b;Leventhal et al. 2012). The neural basis of these ostensibly 'normal' beta oscillations, which are presumably beneficial (Leventhal et al. 2012) and well controlled in time and space by midbrain dopaminergic outputs, is unknown. However, the differences in connections and intrinsic properties of neurons that we have elucidated here could still apply, at least qualitatively, to the STN-GP network in the dopamine-intact brain. Assuming sufficient experimental data are available for model comparison and constraint, our computational strategy and models can be used in the future to investigate the effective connectivity of the STN-GP network in health and, as part of this, extended J Physiol 592.7 to incorporate the activity patterns and connections of midbrain dopaminergic neurons.

Translational perspective
Dopamine loss from the basal ganglia in idiopathic Parkinson's disease (PD) leads to abnormally synchronised, oscillatory neuronal activity in these brain circuits. The disturbed activities of neurons in the subthalamic nucleus (STN) and external globus pallidus (GP) are causally important for the cardinal motor symptoms of PD, but how the connection patterns and intrinsic physiological properties of these neurons might support the emergence of widespread pathological oscillations is still largely unknown. Here, we used electrophysiological recordings of STN and GP neurons in the 6-hydroxydopamine-lesioned rat model of PD to constrain a series of computational models designed to address this issue. Importantly, we took a novel approach that was founded on the premise that GP contains two major types of projection neuron that have dissimilar input-output connections and firing properties. Our optimal computational model, which accurately reproduced our in vivo electrophysiological data, predicted that the major input-output connections of the two types of GP neuron are indeed quantitatively distinct. We further demonstrate that, together with the intrinsic physiological differences of these GP neurons, such mismatched connections explain some of the complex patterns of STN-GP network activity observed in experimental Parkinsonism. These data collectively reiterate the need to incorporate GP functional dichotomy into contemporary schemes of basal ganglia organisation. They additionally suggest that inappropriate interactions between most nodes of basal ganglia-thalamocortical circuits underlie the excessive neuronal oscillations seen in Parkinsonism. This study also provides further rationale for exploring the external globus pallidus as a target for improved surgical interventions in PD.

Supporting Information
Supplementary material