Abstract
The external pallidum (globus pallidus pars externa [GPe]) plays a central role for basal ganglia functions and dynamics and, consequently, has been included in most computational studies of the basal ganglia. These studies considered the GPe as a homogeneous neural population. However, experimental studies have shown that the GPe contains at least two distinct cell types (prototypical and arkypallidal cells). In this work, we provide in silico insight into how pallidal heterogeneity modulates dynamic regimes inside the GPe and how they affect the GPe response to oscillatory input. We derive a mean-field model of the GPe system from a microscopic spiking neural network of recurrently coupled prototypical and arkypallidal neurons. Using bifurcation analysis, we examine the influence of dopamine-dependent changes of intrapallidal connectivity on the GPe dynamics. We find that increased self-inhibition of prototypical cells can induce oscillations, whereas increased inhibition of prototypical cells by arkypallidal cells leads to the emergence of a bistable regime. Furthermore, we show that oscillatory input to the GPe, arriving from striatum, leads to characteristic patterns of cross-frequency coupling observed at the GPe. Based on these findings, we propose two different hypotheses of how dopamine depletion at the GPe may lead to phase-amplitude coupling between the parkinsonian beta rhythm and a GPe-intrinsic γ rhythm. Finally, we show that these findings generalize to realistic spiking neural networks of sparsely coupled Type I excitable GPe neurons.
SIGNIFICANCE STATEMENT Our work provides (1) insight into the theoretical implications of a dichotomous globus pallidus pars externa (GPe) organization, and (2) an exact mean-field model that allows for future investigations of the relationship between GPe spiking activity and local field potential fluctuations. We identify the major phase transitions that the GPe can undergo when subject to static or periodic input and link these phase transitions to the emergence of synchronized oscillations and cross-frequency coupling in the basal ganglia. Because of the close links between our model and experimental findings on the structure and dynamics of prototypical and arkypallidal cells, our results can be used to guide both experimental and computational studies on the role of the GPe for basal ganglia dynamics in health and disease.
Introduction
The basal ganglia (BG) are a set of interconnected subcortical nuclei that form different feedback loops with cortex and thalamus (Alexander and Crutcher, 1990; Bolam et al., 2000). Because of its recurrent connections with nearly all other BG nuclei, the globus pallidus pars externa (GPe) plays a major role in information transmission through the BG (Kita, 2007). In Parkinson's disease (PD), synchronized oscillations have been reported throughout all major BG nuclei (Wichmann, 2019), including the GPe (Wichmann and Soares, 2006; Mallet et al., 2008). These oscillations are characterized by transient power increases in the β frequency band (12-30 Hz) and an increased phase-amplitude coupling (PAC) between the phase of a β signal and the amplitude of a high-frequency γ signal (50-250 Hz) (Jenkinson et al., 2013; Lofredi et al., 2019; Gong et al., 2021). Computational models of BG phase transitions in PD suggest that the GPe is involved in the oscillation generation, either via its recurrent coupling with the subthalamic nucleus (STN) or via its processing of inputs from striatum (STR) (Pavlides et al., 2015; Schroll and Hamker, 2016; Rubin, 2017). Most of these computational models regarded the GPe as a homogeneous population of neurons. However, two major cell types have been identified within the GPe, which differ in their electrophysiological properties, firing rates, and firing patterns: prototypical (GPe-p) and arkypallidal (GPe-a) cells (Cooper and Stanford, 2000; Abdi et al., 2015; Hegeman et al., 2016). Regarding their efferent synapses, it has been shown that GPe-p cells preferentially project to subthalamic nucleus (STN) and BG output nuclei, whereas GPe-a cells provide feedback to striatum (STR) (Mallet et al., 2012; Hernández et al., 2015; Fujiyama et al., 2016). Furthermore, recent studies found that STN and STR differentially affect GPe-p and GPe-a in mice (Pamukcu et al., 2020; Ketzef and Silberberg, 2021). Regarding cell type-specific differences in GPe-intrinsic axon collaterals, there is evidence from mice experiments that prototypical cells express more numerous axon collaterals than arkypallidal cells (Mallet et al., 2012; Higgs et al., 2021; Ketzef and Silberberg, 2021). Still, a substantial number of arkypallidal axon collaterals was identified that targeted prototypical GPe cells (Mallet et al., 2012).
Depending on the pattern of mutual inhibition between those two major GPe cell populations, different modes of GPe internal dynamics may exist. Asymmetric connections between the two cell types may give rise to a feedforward inhibition scenario, where an excitatory input to one population could silence the other population. Alternatively, winner-takes-all (WTA) dynamics can arise in scenarios of mutual inhibition between two populations (Schmidt et al., 2018). Either of these two scenarios could explain the recent findings from optogenetic stimulation in mice, where silencing the GPe-p led to increased GPe-a firing (Aristieta et al., 2021). The same work suggests that STN and STR express asymmetric projections to GPe-p versus GPe-a that allow for transient switching between the two different output pathways of the GPe. To be able to infer the neurodynamic mechanism underlying such data, we argue for an improved understanding of the relationship between synaptic coupling and neural dynamics in a GPe composed of arkypallidal and prototypical cells. Whereas some computational models have considered the role of arkypallidal feedback to STR (Nevado-Holgado et al., 2014; Corbit et al., 2016; Jȩdrzejewski-Szmek et al., 2018), a neurodynamic understanding of the structure-function relationships inside a dichotomous GPe is still missing. In this study, we examine the effects of GPe coupling patterns on GPe behavior. For this purpose, we derive and analyze a mean-field description of two fully coupled inhibitory populations, following the approach by Luke et al. (2013) and Montbrió et al. (2015). Importantly, this mean-field description captures the exact macroscopic dynamics of the underlying, heterogeneous spiking neural network and can thus capture population-intrinsic spike resonance phenomena that classic mean-field approaches would miss. This in itself makes our modeling approach interesting for the understanding of synchronization processes inside the GPe. Using bifurcation analysis of the two-population GPe model, we identify monostable, bistable, and oscillatory regimes, the existence of which depends on the GPe-intrinsic coupling pattern. We then show that the GPe expresses distinct responses to periodic input when initialized in either of these regimes. Finally, we analyze how the macroscopic phase transitions found in the GPe mean-field model translate to spiking neural networks with realistic numbers of neurons and axons.
Materials and Methods
Model definition
Mathematical formulation of population dynamics
We consider the GPe as a nucleus of two distinct populations of GABAergic projection neurons (Kita, 2007; Hegeman et al., 2016). Whereas prototypical neurons express high average spontaneous firing rates of 50-70 Hz (DeLong, 1971; Wichmann and Soares, 2006; Jaeger and Kita, 2011), arkypallidal neurons fire with considerably reduced average firing rates of 5-15 Hz (Abdi et al., 2015; Dodson et al., 2015; Hernández et al., 2015; Aristieta et al., 2021). To model synaptic influences on the spike timings of GPe neurons, it is important to know their type of excitability. This can be inferred from their phase-response curve (Gutkin et al., 2005). Neurons can either express Type I excitability, meaning that the direction in which the excitability of a neuron is changed by extrinsic input is not dependent on the intrinsic phase of the neuron, or express Type II excitability, meaning that the direction in which the excitability of a neuron is changed by extrinsic input does depend on the intrinsic phase of the neuron (Izhikevich, 2000). While computational studies demonstrated that both Type I and Type II excitability can be identified in single-cell models of GPe neurons (Schultheiss et al., 2010; Fujita et al., 2012), experimental investigations only revealed Type I excitability so far (Wilson, 2013). Furthermore, it has been shown that coupled networks of Type I excitable neurons can express Type II excitability on the network level (Dumont and Gutkin, 2019). Thus, as a base neuron model, we use the quadratic integrate-and-fire neuron (QIF), which is the canonical form of Type I excitable neurons and expresses a quadratic and thus nonlinear input-output relationship (Izhikevich, 2000). This choice also accounts for the nonlinear input-output relationship reported in prototypical and arkypallidal cells (Kita, 2007; Abdi et al., 2015). The evolution equation of the jth QIF neuron embedded within either the GPe-p or GPe-a is given by the following:
Here, the synaptic activation m takes the form of a simple convolution of the average firing rate r with the synaptic response kernel g, henceforth abbreviated by the convolution operator ⊗. The parameters η and Δ are the center and half-width at half-maximum of a Lorentzian distribution over the single-neuron parameters ηj. Thus, η and Δ allow to control the average and heterogeneity of the firing rates inside the QIF population, respectively. Spontaneous firing rates of GPe cells cannot be explained by glutamatergic input alone, since brain slice recordings still showed autonomous activity of up to 26 Hz after synaptic transmission was blocked pharmacologically (Günay et al., 2008). In other words, GPe cells are strong pacemaker cells that show regular firing at a cell-specific frequency under synaptic isolation (Mercer et al., 2007; Abrahao et al., 2017). Across GPe cells, a substantial amount of heterogeneity of the intrinsic firing frequencies has been reported (Wilson, 2013). By considering the background excitabilities ηj as distributed quantities, we account for these findings.
We are aware that the all-to-all coupling and infinite population sizes are in contrast to the actual GPe structure (Wilson, 2013; Hegeman et al., 2016). However, it has been recently shown that the mean-field model predictions can generalize to a fairly wide range of network sizes and coupling probabilities (Gast et al., 2020). Even for QIF networks with recurrent coupling probabilities of 1%, the authors found that population sizes of N = 8000 neurons were sufficient to reproduce the macroscopic dynamics predicted by the mean-field model accurately. Given that population sizes of primate GPe are on the order of 105 and recurrent coupling probabilities are ∼5% (Wilson, 2013), we expect that this mean-field model is sufficient to capture the macroscopic dynamics of QIF populations with realistic cell counts and coupling probabilities.
Mathematical formulation of axonal propagation and synaptic dynamics
In a next step, we define the coupling function g which, in our model, acts as a lumped representation of axonal propagation and synaptodendritic integration. In other words, g serves to link single spikes emitted by neuron j to changes in the membrane potential of any other neuron. GPe to GPe connections have been suggested to express axonal transmission delays of ∼1.0 ms (Jaeger and Kita, 2011) and make use of GABAergic synapses (Kita, 2007). Since axon collaterals can express a substantial variability in individual axon diameters and myelination properties (Schmidt and Knösche, 2019), we modeled the axonal transmission delays via γ distributions (Smith, 2011). The probability density function of the γ distribution can be written as follows:
In addition to the axonal delays, we also included a dynamic model of the electrochemical processes that lead to a change in the postsynaptic potential after a presynaptic action potential traveled down the axon. A popular choice to express these dynamics is via a convolution with a biexponential synaptic response kernel, for which the rise and decay time constants are specific to the type of presynapse and postsynapse (Deco et al., 2008). Such a biexponential synaptic response function is given by the following:
Specification of the two-population GPe model
Based on these dynamic equations for neural populations and synaptic transmission, we can now introduce the full set of equations of our GPe model. Since the number of equations of the ODE approximation (Eq. 7) to the γ kernel convolution (Eq. 5) depends on the parameter β of (Eq. 6), we chose to provide a set of integro-differential equations for generality and brevity. However, for our results, each γ kernel convolution was formulated as a set of coupled ODEs of the form (Eq. 7) and each convolution with a synaptic response kernel of the form (Eq. 8) was formulated as the ODE system given by Equations 9 and 10. The following set of coupled integro-differential equations describes the average firing rate and average membrane potential dynamics at GPe-p and GPe-a as follows:
Mathematical formulation of extrinsic model inputs
Extrinsic input can generally be applied via the extrinsic forcing parameters Ip(t) and Ia(t) to GPe-p and GPe-a, respectively. In our simulations, we applied step function inputs to each of the populations. These are defined as follows:
Here, α defines the input strength, whereas
Additionally, to account for the bursting characteristics of typical striatal inputs arriving at the GPe (Jaeger et al., 1995), we applied a sigmoidal transformation to the Stuart-Landau oscillator, giving us the final input as follows:
Model analysis
To analyze the behavior of the model given by Equations 11–14, we used the open-source Python toolbox PyRates (Gast et al., 2019). We chose PyRates' interface to the SciPy Runge-Kutta solver with adaptive integration step-size (Virtanen et al., 2020) for numerical integration of the model dynamics for a given initial condition. For bifurcation analysis, we used PyRates' interface to Auto-07p (Doedel et al., 2007) to perform numerical parameter continuation and automatic bifurcation detection (for an in-depth explanation of these techniques, see Kuznetsov, 2004; Meijer et al., 2009). To analyze the behavior of the spiking neural networks corresponding to our mean-field models, we used custom MATLAB code. Numerical integration of the spiking neural network dynamics was performed via an explicit Euler algorithm with an integration step-size of 0.001 ms, which we found to be sufficiently small to capture all model dynamics. The scripts and configuration files for all simulations and parameter continuations are available at the following public Github repository: https://github.com/Richert/GPe_Dynamics.
Spectral analysis
We also analyzed the GPe model behavior in the frequency domain. To this end, we used time series of 320 s of simulated GPe-p firing rate dynamics sampled at 1 ms and cut off the first 20 s to remove initial transients from the time series. Power spectral densities were calculated from the raw simulation data using Welch's method. We used FFT segments of length 2048 and an overlap between segments of 1024 time steps. For quantification of PAC and phase-phase coupling (PPC) between different frequency components of the GPe-p firing rate dynamics, we followed the procedure described by Gong et al. (2021). PAC measures the amount of modulation of the amplitude of a high-frequency signal by the phase of a low-frequency signal and was evaluated by means of the Kullback-Leibler-based modulation index (KL-MI) (Tort et al., 2010). Both the low- and high-frequency signals were acquired by bandpass filtering the GPe-p firing rate time series. Following the procedure described by Gong et al. (2021), we evaluated the KL-MI for multiple pairs of phases at frequencies
Model parameters
The dynamics at GPe-p and GPe-a are each governed by membrane time constants τ and two parameters η and Δ that determine the center and half-width at half-maximum of the distribution of single-cell firing rates inside the populations. Additionally, the four synaptic connections between GPe-p and GPe-a are each parameterized via a lumped synaptic strength J, two axonal delay parameters μ and σ, and the synaptic rise and decay time constants τr and τd. For all model analysis, we used re-scaled synaptic strength parameters
Model parameters
Configuration of the GPe model. A, Distributions of steady-state QIF firing rates under synaptic isolation, i.e., with
Results
In this section, we report the results of our analysis of the relationship between model parameters and neural dynamics for the GPe model given by Equations 11–14. We focus on parameters that contribute to a difference between GPe-p and GPe-a, which include the coupling strengths within the GPe as well as additional inputs to the two populations. All coupling parameters are reported as
Effects of GPe-intrinsic coupling
As the first part of our analysis, we performed a bifurcation analysis of the GPe mean-field model given by Equations 11–14 to investigate whether different coupling patterns between prototypical (GPe-p) and arkypallidal (GPe-a) cells promote different macroscopic states and phase transitions. We started out from the GPe coupling strengths listed in Table 1, which represent a coupling pattern where GPe-p inhibition of GPe-a is strongest and GPe-a axon collaterals are weaker than GPe-p axon collaterals, but still exist (for a comparison to experimental findings, see Mallet et al., 2012; Ketzef and Silberberg, 2021). In this default state, input to the GPe populations led to changes in their firing rates but did not induce any phase transitions (Fig. 2A). Still, inhibition of GPe-p via decreases in ηp caused a fast increase in GPe-a firing rates. This interaction pattern between GPe-p and GPe-a changes when the intrinsic connections of the GPe are altered. As can be seen in Figure 2C, D, increases in kpp lead to the emergence of a stable limit cycle via a supercritical Andronov-Hopf bifurcation. The emerging oscillations express a frequency in the γ range (≈50-60 Hz) and can be induced by changes in kpp as well as changes in ηp (Fig. 2D,E). The Hopf curve in the ηp-kpp plane shows that this emergence of synchronized oscillations critically depends on kpp > 0 as well as a sufficient excitatory drive as given by ηp > 0. We encountered a different phenomenon when increasing kpa, the connection strength from GPe-a to GPe-p. For sufficient increases in kpa, the system expresses twofold bifurcations that mark the outer boundaries of a bistable regime, in which transient inputs to GPe-p (or GPe-a) allow to switch between two stable states (Fig. 2F,H). One of those two stable states is a focus for which the GPe-p is in a high-activity regime and forces the GPe-a to a low-activity regime. The other stable state is also a focus where the GPe-a is in a high-activity regime and forces the GPe-p to a low-activity regime. These two stable equilibria are separated by a saddle focus. Thus, we found that strong bidirectional coupling between prototypical and arkypallidal GPe populations allows for the existence of a bistable activity regime, where the two populations compete over a high-activity state. In Figure 2G, we show the curve of the fold bifurcations in the ηp-kpa plane, which collapse in a cusp bifurcation when kpa becomes small. Furthermore, we find another Hopf curve that touches the fold curve at a zero-Hopf bifurcation. This Hopf curve covers a relatively small parameter range, however, and only appears on the lower branch within the bistable regime. The oscillations emerging from this limit cycle are small-amplitude γ oscillations.
Phase transitions in the GPe. A, One-parameter bifurcation diagram varying the background input of GPe-p ηp for the default connectivity. B, Depiction of the GPe circuit and the bifurcation parameters. C, 1D parameter continuation in ηp for kpp = 5. Green circles represent Andronov-Hopf bifurcations. Lines starting at the Hopf bifurcation indicate the minima and maxima of the emerging limit cycle. Solid (dotted) lines indicate stable (unstable) equilibria. D, 2D bifurcation diagram in ηp and kpp. Black curve represents the continuation of the 1D Hopf bifurcation from C in the 2D parameter space. Shaded regions represent the parameter space where stable oscillations exist. E, GPe-p and GPe-a firing rates for kpp = 5 and ηp = 30 show that additional input Ip(t) = 10.0, applied between 2200 ≤ t < 2600, forced the system over the Hopf bifurcation. F, 1D parameter continuation in ηp for kpa = 5. Gray triangles represent fold bifurcations. G, 2D bifurcation diagram in ηp and kpa. Gray rhombus represents a cusp bifurcation. Black star represents a zero-Hopf bifurcation. H, GPe-p and GPe-a firing rates for kpa = 5 and ηp = 21.0 show switching between the two stable branches via Ip(t) = 5 applied between 2200 ≤ t < 2400 and Ip(t) = – 5 applied between 2600 ≤ t < 2800.
As a next step, we investigated whether the frequency of the oscillations identified for increased self-inhibition of GPe-p critically depends on the axonal delays in our model, since estimates of these delays differ substantially in the literature (Jaeger and Kita, 2011; Ketzef and Silberberg, 2021). To this end, we simulated the behavior of our model for increasing values of
Effect of increasing axonal delays on GPe oscillations. Time series represent average firing rates of GPe-p, abbreviated as r. Reported values of μ and σ refer to changes the parameter μxy and
GPe response to periodic forcing
By now, we have established an understanding of the intrinsic, coupling-dependent GPe response to static, afferent inputs. We found that a dichotomous organization of the GPe with two distinct populations GPe-a and GPe-p results in coupling-dependent dynamic behavior that situates the GPe either near a bistable or near an oscillatory regime. These two different scenarios may have substantially different consequences for the transmission and amplification of periodic input arriving at the GPe. Hence, as a next step, we analyzed the response of the GPe to periodic inputs when initialized in the bistable regime, in the oscillatory regime, or in the healthy steady-state regime. To this end, we applied periodic striatal input with period ω and amplitude α to the prototypical population. The bursting properties of striatal input were generated by applying a sigmoidal transformation to a Stuart-Landau oscillator (see Eqs. 17–20). In each regime, we performed numerical simulations of the model behavior for different values of ω and α. We then evaluated the average PAC between the phase of low-frequency signal components (2-30 Hz) and the amplitude of high-frequency signal components (50-250 Hz) of the GPe-p firing rate dynamics. Furthermore, we evaluated the PPC, that is, the phase dependency of the high-frequency components on the phase of the dominating low-frequency component. A detailed description of these measures is provided in Materials and Methods. As can be seen in Figure 4, we find that the GPe responds differently to periodic input depending on its dynamic regime.
PAC and PPC in the periodically inhibited GPe. The input Ip(t) was applied with different frequencies
In Figure 4A, the GPe response to periodic input is depicted for the default coupling pattern. In this case, a stable focus is the only equilibrium, and the input perturbs the system around that equilibrium at the input frequency. After a perturbation, the system relaxes back to the focus via damped oscillations. The amplitude of these oscillations scales with α, as can be seen by comparing time series 1 and 2 of Figure 4A. Thus, stronger inputs generate stronger modulation of high-frequency amplitudes by low-frequency phases, resulting in increased PAC. Since the high-frequency focus dynamics are directly elicited by the low-frequency perturbations, increases in PAC always co-occur with increased PPC.
Figure 4B depicts the response of a bistable GPe to periodic input. In this regime, periodic inhibition forces the GPe-p toward the low-activity regime, if sufficiently strong (time series 1 and 2 of Fig. 4B for inputs that are too weak and sufficiently strong, respectively). If forced toward the low-activity regime, the GPe-p attempts to relax back to its natural high-activity regime (time series 2 and 4 of Fig. 4B). In this relaxation process, the system is affected by the strong focus dynamics of the saddle that separates the two stable states. For most combinations of α and ω, this behavior creates oscillations with interleaved large- and small-amplitude oscillations, where the large-amplitude oscillations act as an amplification of the low-frequency input and cause cross-frequency coupling with the high-frequency, small-amplitude focus dynamics. Stronger inputs generate stronger modulation of high-frequency amplitudes by low-frequency phases, as evaluated by PAC. Such increases in PAC occur together with increased phase locking between low- and high-frequency components. This can be observed by the generally high PAC-PPC correlations in Figure 4B. Interestingly, there also exists a relatively narrow window in ω, where the periodic inhibition of the GPe-p forces the system to stay within the domain of influence of the unstable saddle focus, thus causing periodic oscillations with strongly reduced PAC (see time series 3 of Fig. 4B).
More complex, resonant behavior can arise for periodic forcing of the GPe, if the GPe already expresses oscillations autonomously (Fig. 4C,D). When increasing α, the system undergoes a torus bifurcation that emerges from the interaction between the intrinsic limit cycle and the extrinsic, periodic input. As can be seen from the firing rate dynamics in Figure 4C, strong amplitude modulations of the intrinsic limit cycle exist in the vicinity of this torus bifurcation. A continuation of the torus bifurcation in the ω-α plane reveals that the system expresses various resonances of the intrinsic limit cycle with the extrinsic input. Close to regimes of 1:2 resonances, we were able to identify small loci of period doubling bifurcations, suggesting the existence of chaotic regimes. These bifurcations are also reflected in the PAC and PPC profiles of the system. PAC values are low before the system undergoes the torus bifurcation and increase near and after the torus bifurcation. In the vicinity of the torus bifurcation, we find regimes where strong PAC can coexist with low PPC values. These regions express negative correlations between PAC and PPC and are clearly separated from regions where increased PAC and PPC coexist (Fig. 4C).
Model generalization to GPe spiking neural networks
In this section, we report how the above-described findings generalize to spiking neural networks (SNNs) of coupled GPe-p and GPe-a cells with realistic cell counts and coupling probabilities. To this end, we attempted to replicate the mean-field model dynamics shown in time series 3 of Figure 4B and time series 4 of Figure 4C and in SNNs with (1) different network sizes and (2) different coupling probabilities. We created a total of four SNNs with (1) either all-to-all coupling or only 5% of all possible connections, and (2) either Np = 4000 (Na = 2000) GPe-p (GPe-a) cells or Np = 40 000 (Na = 20 000) GPe-p (GPe-a) cells. We then repeated our simulations of the GPe response to periodic stimulation for spiking neural networks initialized near the bistable and in the oscillatory regime (same parameterizations as reported in Fig. 4 for the mean-field model). The dynamics of all four SNNs can be seen compared with the mean-field predictions in Figure 5. As expected, we find that an all-to-all coupled SNN of large size behaves nearly identical to the mean-field prediction, where the remaining difference in the oscillation amplitude is an effect of the network size and would vanish if we increased the network size even further (see difference between SNNs with N1 and N2). Interestingly, we find that reducing the number of synaptic connections to p = 5% of all possible connections attenuates synchronized oscillations in the network for small network sizes. However, for a sufficiently large network, the SNN follows the macroscopic dynamics predicted by the mean-field model, even when p = 5%. This holds for both the bistable as well as the oscillatory regime.
Comparison between mean-field model (dashed lines) and spiking neural networks (solid lines) under periodic stimulation. The first and second columns represent the average firing rate and average membrane potential across the population, whereas the last column represents spike timings from 200 randomly chosen neurons of the population. SNNs are composed of
Discussion
In this work, we investigated the implications of a dichotomous GPe structure on its intrinsic dynamic regimes and on its response to oscillatory input. This was motivated by experimental results that strongly suggest that there exist two neuron types (GPe-p and GPe-a) inside GPe with different projection targets and electrophysiological features (Mallet et al., 2012; Abdi et al., 2015; Hernández et al., 2015; Hegeman et al., 2016). Our investigations were based on populations of coupled QIF neurons, for which we derived exact mean-field equations describing the low-dimensional dynamics of average membrane potentials and firing rates. This two-population model was set up such that it accounts for the I-f curves of GPe-p and GPe-a cells (Abdi et al., 2015), firing rate heterogeneity inside each population (Miguelez et al., 2012; Aristieta et al., 2021; Ketzef and Silberberg, 2021), as well as the response of GPe-p and GPe-a to extrinsic stimulation (Kita et al., 2004; Aristieta et al., 2021). Next, we investigated the dependence of the macroscopic neurodynamics of this model on the underlying connectivity between GPe-p and GPe-a. We found that strong GPe-p projections and weak GPe-a projections produced realistic steady-state firing rates (rp ≈ 60 Hz and ra ≈10 Hz) and responses to STR and STN stimulation (Kita et al., 2004; Dodson et al., 2015; Aristieta et al., 2021). In this regime, the model was relatively robust to changes in its extrinsic input. While alterations of the background input led to changes in the average firing rates, they were not sufficient to induce phase transitions in the model.
Under dopamine depletion, GABAergic postsynaptic currents of GPe-to-GPe synapses have been reported to increase in strength (Miguelez et al., 2012). However, it is unclear which GPe-to-GPe synapses are affected by this and how these changes may influence GPe dynamics. Thus, we investigated whether increases in the strength of specific GPe connections induce changes in the behavior of the system. To this end, we focused on the synapses inhibiting the GPe-p, since increased inhibition of GPe-a would merely reduce its already low steady-state firing rates and, hence, not induce any phase transitions to the GPe dynamics.
When the strength of GPe-p inhibition by GPe-a was increased, we found that a bistable regime emerged in which GPe-p and GPe-a competed over a high-activity state. This regime could provide a form of network memory and a reliable way to switch GPe output between the different projection targets of GPe-p (STN) and GPe-a (STR). Furthermore, if the GPe was situated in the bistable regime, periodic inputs from STR were amplified because of the existence of two different attracting states. Indeed, when we applied periodic input from STR to GPe-p in the β frequency range characteristic for PD (Brown, 2003; Hammond et al., 2007), the GPe-p was periodically forced from a high-activity state down to a low-activity state, thus resonating at the input frequency. Importantly, this led to strong phase-amplitude as well as PPC between β and γ components of the GPe-p firing rate dynamics. Interestingly, the synchronized neural activity that has been detected in recordings of STN and GPe activity from PD patients expressed not only increased power in the β frequency band, but also increased PAC between the phase of β components and the amplitude of γ components (López-Azcárate et al., 2010). Our model can explain these findings as follows: dopamine depletion at the GPe leads to an increased strength of GPe-a to GPe-p synapses in PD (Miguelez et al., 2012). This structural change moves the system closer to a bistable regime. By moving closer to the boundaries of the bistable regime, oscillatory inputs from STN or STR become more likely to elicit switching between the two stable states of the GPe. At both of these input sites, increased β oscillations have been reported in PD (Brown, 2003; Belluscio et al., 2014). Periodic forcing of the bistable GPe would then perturb the system in a phase space controlled by multiple stable and unstable foci with focus frequencies in the γ range, thus causing damped γ oscillations with input-triggered amplitude modulations. Hence, in this scenario, PD-related intrinsic changes can cause increased susceptibility of the GPe to periodic inputs, but not autonomous GPe oscillations.
When we instead increased the GPe-p to GPe-p self-inhibition, we found that stable oscillations in a γ frequency range (≈50-60 Hz) could emerge. These oscillations were driven by the dynamic interactions between the pace-making properties of the GPe-p and its delayed self-inhibition. Most likely, these oscillations reflect the same synchronization mechanism as reported for a single population with delayed self-inhibition (Luccioli et al., 2019). According to their results, oscillations are counteracted by neural heterogeneity. This way, our results can be linked to the considerations in (Wilson, 2013), which suggest that strong firing rate heterogeneity together with recurrent inhibition inside GPe may serve to desynchronize GPe activity under healthy conditions. In accordance with experimental data, we modeled the GPe-p with highly heterogeneous single-cell firing rates (Miguelez et al., 2012; Hernández et al., 2015; Aristieta et al., 2021; Ketzef and Silberberg, 2021). This way, inhibitory feedback from the GPe-p provides the means to suppress synchronized oscillations inside the GPe, which supports these considerations. Experimental evidence from animal models of PD suggest that GPe activity shows increased synchronization in PD (Wichmann and Soares, 2006; Mallet et al., 2012). Our model can explain these findings as follows: dopamine depletion causes the strength of GPe-p self-inhibition to increase in PD (Miguelez et al., 2012). This moves the GPe system closer to or even across the boundaries of an oscillatory regime. The emerging limit cycle leads to narrow-band γ oscillations and thus cannot explain the emergence of β oscillations in the parkinsonian BG (Brown, 2003; Hammond et al., 2007). However, assuming that burst-like afferent inputs drive the GPe at a β frequency, our findings predict that GPe-intrinsic γ oscillations can resonate with the input, leading to a waxing and waning of the γ oscillations. Such waxing-and-waning behavior also implicates increased PAC, which may occur together with decreased PPC, according to our simulations. The latter finding was unique for the oscillatory regime of the GPe. Similarly, complex patterns of cross-frequency coupling have been reported previously in an instantaneously coupled two-population QIF model with sinusoidal forcing in the α frequency range (10 Hz) (Ceni et al., 2020). Thus, our results show under which conditions the GPe system can express the characteristic dynamics that have been identified in more abstract models of two populations with mutual inhibition.
The mathematical model presented in this paper can serve as a basis for future BG models. We have shown that a mean-field model derived under the assumptions of all-to-all coupling and an infinite number of neurons accurately describes the dynamics of a SNN with realistic cell counts and coupling densities (Hegeman et al., 2016). Thus, our model lends itself to multiscale approaches. It can easily be extended by additional biological details, such as plasticity mechanisms (Gast et al., 2020) or gap junctions (Pietras et al., 2019). The latter may be of particular interest to investigate parkinsonian conditions inside the GPe (Schwab et al., 2014). Furthermore, our model can be used as a basis for experimental investigations of GPe structure-function relationships in healthy and diseased states. First, the I-f curves reported at the single-cell and population level allow a comparison of our model with experimental data on the input-output relationship of GPe cells. Second, our model predicts specific phase transitions of GPe dynamics to emerge for increased strengths of GPe-p to GPe-p or GPe-a to GPe-p synapses. Experimental investigation of the GPe response to GPe-p stimulation under dopamine-depleted conditions could serve to infer whether GPe dopamine depletion increases these synaptic strengths sufficiently to induce these phase transitions. Finally, our results indicate that the GPe expresses distinct responses to periodic input from STR, based on its dynamic regime. Specifically, PAC and PPC profiles of GPe firing rates can distinguish between the healthy steady-state regime, the bistable regime, and the oscillatory regime of the periodically forced GPe. Thus, extrinsic periodic stimulation of indirect pathway STR projection neurons could be used in combination with measurements of PAC and PPC of GPe firing rates to infer the underlying dynamic regime and synaptic connectivity of the GPe in healthy versus diseased states.
Footnotes
R. Gast was supported by Studienstiftung des deutschen Volkes. H.S. was supported by German Research Foundation DFG KN 588/7-1 awarded to T.R.K. via Priority Program 2041, “Computational Connectomics.” We thank Bastian Pietras and Ernest Montbrió for helpful discussions.
The authors declare no competing financial interests.
- Correspondence should be addressed to Richard Gast at rgast{at}cbs.mpg.de
