QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

HOME  |   SEARCH  |   ARCHIVE  |   SUBSCRIBE  |   CONTACT  |   HELP

The Journal of Neuroscience, March 29, 2006, 26(13):3567-3583; doi:10.1523/JNEUROSCI.5050-05.2006

This Article Abstract Full Text (PDF) Submit an eLetter Alert me when this article is cited Alert me when eLetters are posted Alert me if a correction is posted Citation Map Services Email this article to a friend Related articles in J. Neurosci. Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (13) Citing Articles via Google Scholar Google Scholar Articles by Leblois, A. Articles by Hansel, D. Search for Related Content PubMed PubMed Citation Articles by Leblois, A. Articles by Hansel, D.

 Previous Article  |  Next Article 

Behavioral/Systems/Cognitive
Competition between Feedback Loops Underlies Normal and Pathological Dynamics in the Basal Ganglia

Arthur Leblois,^1,2,4 Thomas Boraud,^2,4 Wassilios Meissner,² Hagai Bergman,^3,4 and David Hansel^1,4

¹Laboratoire de Neurophysique et Physiologie du Système Moteur, Centre National de la Recherche Scientifique (CNRS) Unité Mixte de Recherche (UMR) 8119, Université René Descartes, 75270 Paris, France, ²Basal Gang, Laboratoire de Neurophysiologie, CNRS UMR 5543, Université Victor Segalen, 33076 Bordeaux, France, ³Interdisciplinary Center for Neural Computation, The Hebrew University, Jerusalem 91904, Israel, and ⁴Laboratoire Franco-Israelien de Neurophysiologie et Neurophysique des Systèmes, Université René Descartes, 75270 Paris, France

	Abstract

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
Experiments performed in normal animals suggest that the basal ganglia (BG) are crucial in motor program selection. BG are also involved in movement disorders. In particular, BG neuronal activity in parkinsonian animals and patients is more oscillatory and more synchronous than in normal individuals.

We propose a new model for the function and dysfunction of the motor part of BG. We hypothesize that the striatum, the subthalamic nucleus, the internal pallidum (GPi), the thalamus, and the cortex are involved in closed feedback loops. The direct (cortex–striatum–GPi–thalamus–cortex) and the hyperdirect loops (cortex–subthalamic nucleus–GPi–thalamus–cortex), which have different polarities, play a key role in the model. We show that the competition between these two loops provides the BG–cortex system with the ability to perform motor program selection. Under the assumption that dopamine potentiates corticostriatal synaptic transmission, we demonstrate that, in our model, moderate dopamine depletion leads to a complete loss of action selection ability. High depletion can lead to synchronous oscillations. These modifications of the network dynamical state stem from an imbalance between the feedback in the direct and hyperdirect loops when dopamine is depleted.

Our model predicts that the loss of selection ability occurs before the appearance of oscillations, suggesting that Parkinson's disease motor impairments are not directly related to abnormal oscillatory activity. Another major prediction of our model is that synchronous oscillations driven by the hyperdirect loop appear in BG after inactivation of the striatum.

Key words: neural network; models; action selection; oscillations; synchrony; Parkinson's disease

	Introduction

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
The basal ganglia (BG) are a complex network of subcortical nuclei involved in motor control, sensorimotor integration, and in cognitive and motivational processes (Gerfen and Wilson, 1996; Bolam et al., 2000). Their functions remain a subject of debate. Several experiments and anatomical considerations suggest that BG may be involved in the selection of motor programs (Chevalier and Deniau, 1990; Mink and Thach, 1991; Mink, 1996; Turner and Anderson, 1997). Other experiments lead to the view that BG play a critical role in reinforcement learning (Apicella et al., 1991; Schultz et al., 1993; Bar-Gad and Bergman, 2001).

Parkinson's disease (PD) and Huntington's disease involve BG dysfunctions. A model of motor symptoms of these diseases was proposed by Albin et al. (1989) and DeLong (1990). It relies on a segregation between the direct and indirect pathways going from the striatum to BG output structures. It also assumes that dopamine (DA) has opposing effects in these two pathways mediated through D₁ and D₂ receptors, respectively. This model predicts successfully that inactivation of internal pallidum (GPi) or subthalamic nucleus (STN) palliates PD symptoms (Bergman et al., 1990; Benazzouz et al., 1993). However, it is not sufficient to account for the modifications in firing patterns observed over the BG–cortical network after DA depletion (Bergman et al., 1994; Nini et al., 1995; Hutchison et al., 1997; Raz et al., 2001). Moreover, striatal neurons are involved in both the direct and indirect pathways (Bolam et al., 2000; Wu et al., 2000; Levesque and Parent, 2005), D₁ and D₂ receptors are colocalized on striatal neurons (Aizman et al., 2000), and the effects of DA mediated through D₁ and D₂ receptors are not always opposing (Calabresi et al., 2000; Kerr and Wickens, 2001; Nicola et al., 2004). Last but not least, a recent study concludes that lesions in the external pallidum in the monkey do not induce parkinsonian-like motor symptoms and do not affect the activity pattern in GPi (Soares et al., 2004). This raises additional questions regarding the primary involvement of the indirect pathway in PD symptoms.

Here, we propose a new model for the motor part of BG, which does not require a segregation between the direct and indirect pathways. Instead, the model assigns a primary role to the hyperdirect pathway, which goes from the cortex to the STN and the GPi (Nambu et al., 2000). It also assumes that the direct and hyperdirect pathways operate as competing feedback loops (Deniau et al., 1996; Hoover and Strick, 1999; Kelly and Strick, 2004). Hence, it is a dynamical model of the BG–cortical network. We show that the competition between these two loops may be a basis for BG functions and dysfunctions. Specifically, it provides a natural mechanism for action selection and its loss after DA depletion. It also suggests that pathological synchronous oscillations can stem from a dynamical imbalance between the direct and hyperdirect pathways. Part of this work has been presented previously in abstract form (Leblois et al., 2002, 2003).

	Materials and Methods

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
The model
In the present study, we focus on the implications of the direct and hyperdirect pathways in BG physiology and pathophysiology. Therefore, in most of our study, the indirect pathway and the external segment of the globus pallidus (GPe) are not included. Hence, our model of the BG network contains five populations of neurons. Three of these populations are excitatory. These are the motor cortex (Ctx), the ventral anterior and ventral lateral nucleus of the thalamus (Th), and the STN. Two populations are inhibitory. They correspond to the sensory–motor territories of the striatum (Str) and the GPi. The impact of adding a GPe population to the network is briefly addressed in Discussion.

Anatomical and electrophysiological studies have shown that the connectivity along the direct pathway is topographically organized and that the flow of information through this pathway displays functional and somatotopic segregation (Alexander et al., 1986; Deniau et al., 1996; Nakano, 2000). In line with these features, our BG network model consists of two parallel circuits that control two motor programs (Fig. 1). It can be thought of as one "somatotopic channel" involved in the execution of these two motor programs. We also assume that the two circuits interact at the level of the STN to GPi connection. This hypothesis is supported by anatomy and electrophysiological evidence, which shows that the connectivity between STN and GPi is more divergent than the connectivity along the direct pathway, and that the functional segregation observed along the direct pathway is only partly maintained along the hyperdirect pathway (Parent and Hazrati, 1995b). The architecture of our model is summarized in Figure 1. In particular, each STN population interacts with the GPi populations in the two circuits. A motor program is executed only if the average activity of the cortical population crosses some defined threshold in the corresponding circuit.

View larger version (19K): [in this window] [in a new window]   Figure 1. Architecture of the model. The network consists of two circuits, each comprising a cortical, a striatal, a thalamic, a subthalamic, and a pallidal population. The two circuits interact via diffused subthalamic–pallidal connections. Arrows, Excitatory connections. Dots, Inhibitory connections. The substancia nigra pars compacta is not explicitly represented in the model.

 
It is also clear from this figure that each circuit can be thought of as being composed of two feedback loops as follows: (1) A global positive feedback loop: CtxStrGPiThCtx. In this loop, the cortex acts on the thalamus through a sequence of one excitatory and two inhibitory sets of connections. Therefore, the cortex disinhibits the thalamus. Because the CtxStr GPi path is usually called the direct pathway, this loop will be referred to in the following as the "direct loop." (2) A global negative feedback loop: CtxSTNGPiTh Ctx. In this loop, the cortex acts on the thalamus through a sequence of three sets of connections: two excitatory (CtxSTNGPi) and one inhibitory (GPiTh). Therefore, it inhibits the thalamus. Following Nambu et al. (2000), we call this loop the "hyperdirect loop."

Neuronal dynamics
The number, N, and the cell properties of the neurons in a given population ( = Ctx, Str, GPi, Th, STN) are assumed to be identical in the two circuits. The single neuron dynamics are described by a rate model (Wilson and Cowan, 1972; Hopfield, 1984; Shriki et al., 2003). A neuron in population in circuit k is characterized by its instantaneous activity A_ik(t), where i = 1, ..., N, and k = 1, 2. If I_ik is the total input to the neuron and S_ik(x) is its nonlinear input–output transfer function, its instantaneous activity is given by the following: A_ik(t) = S_ik(I_ik). For simplicity, we will consider threshold linear transfer functions, S_ik(x) = [x – T_i]₊, where [x]₊ = x for x > 0 and 0 otherwise; T_i is the threshold of neuron i in population , and > 0 is the gain of the input–output transfer function of the neurons in that population. For simplicity, we take T_i = T independently of i for = Ctx, GPi, Th, STN. In contrast, we take a broad Gaussian distribution for the thresholds of the striatal neurons (average, T_Str; SD, T_Str/2) to reproduce the distribution of spontaneous activity of the striatal projection neurons observed in experiments (Sandstrom and Rebec, 2003; Slaght et al., 2004), as shown in Results.

The synaptic output of neuron (i, , k) to a neuron in population is characterized by a smooth variable m_ik(t), which is a low-pass filtered version of its instantaneous level of activity (Ermentrout, 1996; Shriki et al., 2003). This variable follows the dynamical equation as follows:

where , which only depends on and , is the time constant of the synapses between neurons in population and in population . We will assume that = for (, ) (Ctx, STN). To account for the presence of NMDA receptors on STN neurons, _CtxSTN is assumed to be larger than . The ratio _CtxSTN/ will be denoted by µ.

Connectivity
The connectivity from population to population in circuit k = 1, 2 is random. The connectivity matrix, J_ijk, is such that J_ijk = 1 if neuron j, , k is connected presynaptically to neuron i, , k, and J_ijk = 0 otherwise. The average number of inputs from population to neuron i, , k depends solely on and and will be denoted by K. The strength of the interaction between two neurons in the same circuit depends solely on the populations to which they belong. It will be denoted G. As stated above, populations belonging to different circuits do not interact, except for the STN and the GPi. The strength of the cross-connection between the STN in one circuit and the GPi in the other circuit is G_GPiSTN, where is a constant and G_GPiSTN is the strength of interaction between the STN and the GPi in the same circuit. The connectivity along the cross-connection from circuit k to circuit k' is denoted by J_ij^crossk'k.

For the sake of simplicity, we do not introduce any lateral interactions in the various populations of the model.

Total inputs to the neurons
The total synaptic inputs I_ik received by neuron i in population (for = Ctx, Str, GPi, Th, STN) in circuit k, are given by the following:

where denotes the synaptic delay from population to population . The cortical and striatal populations receive external inputs (see below) denoted by H_ik^Ctx and H_ik^Str, respectively. The last term in these equations, _ik(t), represents additional Gaussian white noisy inputs with a zero mean and a SD, .

External inputs to the network in relation to the planning and execution of movements. We assume that, in relation to the planning and execution of a movement, the neurons in the cortical and the striatal populations receive additional external inputs. These inputs represent synaptic entries to cortical and striatal neurons coming from other cortical areas and brain regions that are not explicitly incorporated in the model. For simplicity, the inputs to the cortical and striatal populations in circuit k (k = 1, 2) are taken to be homogeneous over the populations. In the text and figures, external inputs are displayed in arbitrary units.

The inputs to cortical neurons. The two cortical populations in our model represent neurons in the motor cortex encoding for two motor programs. Because other brain regions sending inputs to the motor cortex are not explicitly represented in the model, we take these inputs into account by including an external input, H_ik^Ctx, in the dynamics. We assume that the temporal profile of this external input has the following form:

for t_m – D_mvt/2 < t < t_m + D_mvt/2 and zero otherwise. Hence, D_mvt is the total duration of the input related to the movement and t_m is the time at which this input takes its maximal value, H_k^Ctx.

The degree of selectivity of the external inputs to the cortex (i.e., the extent to which the external inputs to the two cortical populations in the two circuits are different) will be measured by the following:

A weakly selective external input to the cortex corresponds to a small . In reality, the input to the cortex related to action execution is presumably selective to some extent. However, our idea is that this selectivity may not be sufficient to induce action selection in normal situations but that the BG network is required for it to occur. To highlight this role of the BG, we assume in most of this study that the cortical input is nonselective [i.e., identical in the two cortical populations ( = 0)]. We will show that even in this case the BG can make the response of the cortex selective.

The inputs to striatal neurons. During motor planning, information coming from sensorimotor cortical area, and related to the motor program to be executed, is sent to the sensory–motor striatum (Gardiner and Nelson, 1992; Boussaoud and Kermadi, 1997; Lee and Assad, 2003). This sensory information is represented in the model by adding a transient and weakly selective external input to the striatum. This input starts at the same time as the cortical input, but its duration, D^Str, is much shorter than the duration of the cortical input (D_mvt). We take the following for its profile:

for t_m – D_mvt/2 < t < t_m – D_mvt/2 + d_Istr, where d_Istr = 200 ms is the duration of this input, and zero otherwise.

Modeling the effect of DA
Dopaminergic SNc neurons project to several BG nuclei (for review, see Graybiel, 2000), but these projections and the dopaminergic receptors are more numerous in the striatum (Haber and Fudge, 1997). Thus, we have assumed that the DA primarily affects the BG dynamics via its effects in the striatum. The effects of DA on synaptic and cellular properties in the striatum is still a topic of debate (for review, see Calabresi et al., 2000; O’Donnell, 2003; Nicola et al., 2004). Most of the available data regarding the presynaptic and postsynaptic effects of DA were obtained in vitro. These data are often contradictory, presumably because of the different experimental conditions under which they were obtained, and whether the changes in DA concentration are phasic and/or attributable to local fluctuations or whether they are pathological. Here, we focus on the pathological changes in the dynamics of the cortex–BG network induced by a depletion in the average extracellular DA concentration. Previous in vitro studies suggest that the threshold of striatal neurons increases with DA (Calabresi et al., 2000; O’Donnell, 2003). Consistent with this idea, in vivo studies have shown that the activity of striatal projection neurons increases after DA depletion both in anesthetized and awake animals (Kish et al., 1999; Tseng et al., 2001). We modeled this increase in activity after DA depletion phenomenologically by the following dependence of the average threshold of striatal neurons, T_Str, on the level of DA as follows:

where D is the relative level of striatal DA compared with its normal level (D = 100%). This function is plotted in Figure 2A. In addition, corticostriatal transmission is altered by DA depletion (Calabresi et al., 2000) and the signal-to-noise ratio in the striatum increases with the level of DA (O’Donnell, 2003; Nicola et al., 2004). Although little is known about the change in corticostriatal synaptic strength with DA, it seems reasonable to assume that the DA potentiates this synaptic connection. In fact, the combined increase in corticostriatal synaptic strength and average striatal threshold would result in a higher signal-to-noise ratio in the striatum. Moreover, the phasic activation induced by glutamate in striatal neurons is amplified by DA (Kiyatkin and Rebec, 1996), suggesting a postsynaptic potentiation of the corticostriatal synapses. We thus consider variations of the corticostriatal synaptic strength with DA of the following form:

as shown on Figure 2B.

View larger version (14K): [in this window] [in a new window]   Figure 2. The effects of DA. A, The effect on the threshold of striatal neurons (Eq. 10). B, The effect on the effective strength of the corticostriatal synapses (Eq. 11). D = 100% corresponds to the "normal" physiological level of DA.

 
Note that the two effects of DA depletion captured by Equations 10 and 11 are antagonistic because the first increases the response of the striatal neurons to an elevation of cortical activity, whereas the second effect decreases it.

Numerical simulations of the detailed model
The complexity of the model described above makes it very hard to study analytically. For this reason, we will investigate it with numerical simulations.

The BG–cortical circuitry is characterized by high probabilities of connections between interacting neurons. Recent estimates of connectivity from cortex to striatum and STN indicate that the number of cortical neurons connected to given neurons in these two structures are in the order of K_StrCtx = 5000 (Kincaid et al., 1998) and K_STNCtx = 100 (Mink, 1996), respectively. The connectivity from the striatum to the GPi, the STN to GPi, and the GPi to the thalamus are less well known. In our simulations, we took K_GPiSTN = 800 and K_ThGPi = 500. For the connectivity between the striatum and GPi, we chose K_GPiStr = 50. We also checked that all the results were unchanged if K_GPiStr = 500. To account for the large connectivity from the thalamus to cortex, we took K_CtxTh = 1000.

Obviously, it is unfeasible to simulate a model of the cortex–BG network that would incorporate a number of neurons in each of the populations that would be similar to reality. The issue is thus how best to scale the network to a size suitable for simulations with reasonable time and memory consumption. In the present study, we used the approach developed by Golomb and Hansel (2000). The idea is that a network in which all the populations have the same size, N_sim, will display properties very close to those of the original network, provided that the connectivity between population and population is scaled as follows:

where N and K are the number of neurons in population and the average connectivity between populations and in the original system, and K^sim is the average connectivity between the two populations in the simulated scaled network.

In our simulations, we took N_sim = 1000. The scaled connectivity parameters are subsequently K_CtxStr^sim = 909, K_CtxSTN^sim = 92, K_STNGPi^sim = 446, K_GPiTh^sim = 333, K_ThCtx^sim = 500, and K_StrGPi^sim = 48.

Integration method. The simulations of the model were performed by integrating the dynamics using a standard first-order Euler algorithm with fixed time step t = 0.5 ms (Press et al., 1993). We also ran several simulations with a smaller time step (0.05 ms) to check that the dynamics were unchanged.

Single-unit and population activity. Spike trains for a given neuron (i, , k) were generated according to an inhomogeneous Poisson process with an average rate given by its instantaneous firing rate, A_ik. Population activities were calculated by averaging the instantaneous firing rates over one set of 20 units taken randomly in the populations.

PETH analysis of movement related activity. To characterize the movement-related activity in our model, we generated the perievent time histogram (PETH) of a single unit around the movement initiation. To that end, we repeated N_trial = 30 times the simulation of the system when the movement related external inputs, H^ctx(t), H^Str(t), were applied. The initial conditions as well as the noise were different from trial to trial. The PETH was computed by averaging the activities in these trials centered on the onset of movement related input, t_m – D_mvt/2 (see above), as is usually done to analyze experimental data. We proceeded as follows. For each trial i, the neuron firing-rate time course F_i(t) was first determined with a time in 10 ms by a kernel estimator in which the spike times T_i^j were convolved with a kernel function K(t): F_i(t) = _j=1ⁿ K(t – T_i^j). We used a Gaussian kernel K(t) = , where s determined the kernel width, controlling the degree of smoothing. We took s = 0.25/F, where F is the mean firing rate of the neuron over the recording period (Baker and Gerstein, 2001). The mean firing rate of the neuron across the trials, aligned on the onset of cortical additional input, was then calculated, and yielded a smoothed version of the standard PETH. The mean and SD of the mean rate were determined over a baseline region (during the 500 ms preceding onset of additional input to cortex). The onset time was defined as the first bin where the estimate rate was modified by 10% from the mean in the same direction (i.e., elevation or suppression). Neurons were classified according to the polarity of the response: inhibition or activation.

Spectral analysis. Spectral analysis was performed on spike trains formed over 100 s of unitary activity in the full-DA (D = 100%) and DA-depleted (D = 20%) states. Single neuron oscillatory behavior was analyzed using power spectral analysis of the spike trains, whereas synchronized oscillatory activity between neurons was detected using coherence spectra between pairs of spike trains. Both autospectra and coherence functions were calculated with Hanning windows and mean detrending of data segments. Coherence is a function of frequency and is calculated from the cross-spectral density between the two waveforms normalized by the power spectral density of each waveform. Coherence values can range from 0 if the spike trains are not linearly related to a value of 1 if the spike trains have a perfectly linear relationship. A statistically significant coherence value between the discharge of two neurons was used to indicate the presence of oscillatory synchronization. A 95% confidence level was determined by calculating a coherence value given by 1–(1– , where a = 0.95, L is the number of windows used, and the factor 0.375 results from Hanning window weighting (Halliday and Rosenberg, 2000). This value or greater was considered to indicate a significant probability (p < 0.05) of synchronized oscillatory activity between two cells. For statistical analysis of the autospectra, the elements of the interspike intervals of each spike train were randomly shuffled, and the power spectrum of the resulting (random) spike train was calculated (Soares et al., 2004). The mean and SD of the resulting power spectra of 20 iterations of shuffled data were computed. Peaks in the power spectra of the original data stream were considered significant if they were at least 5 SD greater than the mean of the shuffled data. An autospectrum (resp. coherence function) displaying one or more significant peak(s) is referred to an oscillatory spectrum [respectively (resp.) coherence].

The reduced model
We will also consider a simplified version of the model described above, which has the advantage that several aspects of its dynamics can be investigated analytically. The knowledge of the properties of this reduced model will guide us in the investigation of the detailed model.

In the reduced model, the connectivity between two interacting populations and ( presynaptic to ) is all-to-all (e.g., J_ijk = J = 1 for any neuron i and j in population and , respectively). To make the analysis more tractable, we also assumed that there is no noise in the system. The thresholds are assumed to be the same for all the neurons in a given population. Therefore, all the neurons in a given population receive the same total synaptic input. Hence, the dynamics of the system can be fully described in terms of 12 activity variables, m_k(t).

The steady states of the network for time-independent external inputs as well as their stability can be determined analytically as a function of the synaptic weights, the external inputs, synaptic delays, and synaptic time constants.

The steady states of the network
Steady states of the network are obtained by solving the fixed point equations for the dynamics. Here, we focus on the steady states when the external inputs are nonselective (i.e., the external inputs are the same for the two circuits). Obviously, by symmetry, solutions of the fixed point equations in which the activities in the two circuits are the same always exist. We will call such solutions "symmetric" steady states. In particular, states in which the total inputs are above threshold for all the populations are symmetric. In these states, referred below as "linear steady states," the activities are determined by a set of linear equations that can be solved straightforwardly (see Appendix, Eqs. 26272829–30). To be consistent, such a solution must fulfill the condition that all the inputs are above threshold (five conditions).

Other symmetric steady states, in which the total input to at least one of the neuronal populations in each circuit is under threshold, may exist. In such states, at least one of the feedback loops is open. For a given set of synaptic parameters, several such steady states may coexist. They can also coexist with the linear steady state.

Under certain conditions that are analyzed below, symmetry breaking occurs and the response of the network is asymmetric (populations from the two circuits display different levels of activity) despite the fact that the external inputs are the same on both circuits. This second type of solution exists because of the threshold nonlinearities of the transfer functions. As we will see below, these solutions play a central role in the ability of the BG network to perform action selection.

Stability of steady states
A fixed point solution of the dynamical equations is stable if any small perturbation around it eventually decays at large time. If certain perturbations increase with time, the fixed point is unstable.

To investigate the stability of a steady state, we need to study the equations of the dynamics linearized around that state and look for solutions of the following form (Strogatz, 1994):

This is a solution for the linearized dynamics, provided that satisfies a transcendental equation, P() = 0, which depends on the coupling strengths, delays, and synaptic time constants. The solutions to this equation are in general complex numbers. The steady state is stable provided that Re() < 0 for all the solutions. It is unstable if at least one solution with Re() > 0 exists. If, for this solution, Im() = 0, the system undergoes a nonoscillatory instability. If Im() 0, the instability is an Hopf bifurcation (Strogatz, 1994) at a frequency = Im(). In general, at the onset of an instability, the spatial structure of the unstable mode is obtained by determining the prefactors _k.

In the case of linear steady states, one finds that (for details, see Appendix):

with:

and

Note that G₊ and G_– are the products of the synaptic strength along the direct (positive feedback) and hyperdirect (negative feedback) loops respectively and that ₊ and _– are the total delays along these two loops.

In other steady states that play an important role in our analysis, the cortical population is silent in one of the circuits, whereas all of the populations in the other circuit are active. These states are nonsymmetric. For these states, the function P() is as follows (for details, see Appendix):

Finally, note that steady states in which all the feedback loops are open are always stable.

	Results

Top Abstract Introduction Materials and Methods Results Discussion Conclusions Appendix References

 
The properties of the reduced model
The phase diagram of the reduced model for nonselective time-independent external inputs
In this section, we consider the case of nonselective external inputs [i.e., we assume that the inputs to the cortex are the same in the two circuits ( = 0) and that there is no external input to the striatum (H_k^Str = 0 for k = 1,2)].

When the total synaptic input to the cortical populations is smaller than the threshold T_Ctx, they are inactive. This happens in particular if the external input to the cortex is too small, for H^Ctx < H_rest^Ctx, where H_rest^Ctx is given by Equation 32 in the Appendix. Note that, if H^Ctx satisfies this condition, the state of the network does not depend on H^Ctx. Note also that the rest state is always stable, because all of the feedback loops are open in that state.

For sufficiently large external input to the cortex, the feedback is suppressed in the direct and hyperdirect loops, because either the GPi or the thalamus becomes inactive. If G₊ > G_–, the GPi is silent in the two circuits for H^Ctx > H_max1^Ctx, where H_max1^Ctx is given by Equation 33 in the Appendix. Similarly, if G_– > G₊ and H^Ctx > H_max2^Ctx, where H_max2^Ctx is given by Equation 34 in the Appendix, the thalamus is strongly inhibited by the GPi and it becomes inactive in the two circuits.

For intermediate values of the external input to the cortex, H^Ctx, (H_rest^Ctx < H^Ctx < H_max1^Ctx if G_– < G₊ or H_rest^Ctx < H^Ctx < H_max2^Ctx if G₊ < G_–), the linear steady state in which all the populations are active in the network exists. As explained in Materials and Methods, in this state, the activities of the various populations are determined by a set of linear equations (see Appendix, Eqs. 26272829–30). This state is symmetric (i.e., the activities in the populations are the same in the two circuits). In particular, the activities of the cortex in the two circuits are as follows:

where I₀ is given by Equation 31 in the Appendix.

We now consider the stability of the linear steady state (see Materials and Methods).

Nonoscillatory instabilities. Such an instability occurs when a real eigenvalue, , changes its sign from negative to positive as a control parameter is varied. This happens if P(0) = 0 [for the definition of the function P(), see Eqs. 14 and 20]. Two nonoscillatory instabilities may occur, one for the following:

and the other for the following:

Calculation of the coefficients _k (Eq. 13) shows that, for the first of these instabilities, the unstable mode is inhomogeneous [i.e., this instability breaks the symmetry between the two circuits (Fig. 3A)]. This symmetry-breaking instability requires strong feedback in the direct loop but also sufficiently strong negative feedback in the hyperdirect loops. Its occurrence can be intuitively understood as follows. If activity is increased in the cortical population of one of the two circuits, the strong positive feedback in this circuit amplifies this increase. In contrast, the cortical populations in the two circuits tend to inhibit each other via the negative polarity "cross-path" Ctx–STN–GPi–Th–Ctx, which goes from the cortical population in one circuit to the other one in the other circuit. As a result, the increase in the activity of one cortical population reduces the total input received by the other one, and hence its activity is reduced. If this effect is sufficiently strong, it breaks the symmetry between the two circuits.

View larger version (13K): [in this window] [in a new window]   Figure 3. Instabilities of the symmetric fixed point solution of Equations 12345–6. In A and B, the network settles at the unstable symmetric fixed point at t = –500 ms. It remains there until t = 0 when the cortical population in circuit 1 is perturbed by a brief external current. After a short while, an instability develops. Red, Cortex; black, thalamus; blue, GPi. Solid (resp. dashed) lines, The activities in circuit 1 (resp. circuit 2). A, Symmetry-breaking instability for G+ = 2.47, G– = 2.85. The network ends in an asymmetric state in which the cortical population is active in one circuit alone. B, Oscillatory instability for G+ = 0.18 and G– = 2.85. The network settles in a homogeneous oscillatory state. C, Solid line, The frequency of the unstable mode at instability onset as a function of the overall delay derived from Equation 36. Circles, The frequency of the oscillations as a function of at the instability onset for G– = 2.85.

 
Oscillatory instabilities. An oscillatory Hopf instability occurs if the real part of a complex eigenvalue changes sign when a control parameter is varied. Therefore, at the onset of an oscillatory instability, = i is pure imaginary. The complex equation P(i) = 0 determines the values of and the conditions that G₊, G_–, ₊, _–, µ satisfy at the instability onset.

The instability condition simplifies if one assumes that ₊ = _– and µ = 1. Then one finds that the linear steady state undergoes a Hopf instability for the following:

where G₀ and are determined by Equations 36 and 37 (see Appendix). In particular, the frequencies of the unstable modes depend only on the synaptic delays. For a given value of , Equation 36 has an infinite discrete set of solutions, _n, n = 1, 2, ...; ₁ < ₂ < .... The solution _n corresponds to a mode that is unstable when G_– is larger than the critical value G_–^c = G₊ + G₀(, _n). Hence, the stationary fixed point first becomes unstable via the mode n = 1, for which G_–^c is the smallest. The frequency of this mode at instability onset is plotted in Figure 3C as a function of . It decreases monotonically with . For = 0, = 31.8 Hz. For = 20 ms, = 12.8 Hz.

Using simulations, we determined the frequency of the actual oscillations that developed at the instability onset for G_– = 2.85 for several values of . As shown in Figure 3C (circles), in the range of delays we explored, the frequency is extremely well approximated by the imaginary part of the instability eigenvalue.

The oscillatory instability (Fig. 3B) is driven by the negative feedback in the hyperdirect loops, which induces oscillations if it is sufficiently strong. This is because an increase in the activity of the cortical populations in the two circuits leads, through this negative feedback, to a decrease in the input that feeds back from the network to both cortical populations. This decrease opposes the initial increase in cortical activities. If it is strong enough, oscillations occur. The positive feedback present in the direct loop compensates for this tendency. This is the meaning of Equation 24.

Clearly, in the unstable mode, oscillations are present in all of the populations in each circuit. The oscillations of the two circuits are in-phase, but different populations in the same circuit oscillate with some phase shifts that are functions of the various synaptic delays.

The phase diagram of the reduced model. The results described above concerning the stability of the linear steady state can be summarized in a phase diagram in the plane (G₊, G_–) like the one shown in Figure 4. This phase diagram corresponds to the case ₊ = 26 ms, _– = 20 ms, = 5 ms, and µ = 4. The linear steady state is stable in the region between the dashed and the solid line.

View larger version (19K): [in this window] [in a new window]   Figure 4. The phase diagram for the various dynamical regimens of the reduced model as a function of G+ and G– for = 0.4, + = 26 ms and – = 20 ms. The synaptic time constant is = 5 ms for all of the synapses except synapses from the cortex to STN for which STNCtx = 20 ms.

 
On the solid line, the linear steady state undergoes a symmetry-breaking instability. Below this line, G₊ > G_–, and the symmetric linear steady state exists if the external input to the cortex is such that H_rest^Ctx < H^Ctx < H_max1^Ctx, where H_rest^Ctx and H_max1^Ctx are given by Equations 32 and 33 in the Appendix. However, because of the symmetry-breaking instability in this region, any small asymmetric perturbation leads the system to an asymmetric state in which the cortical population is active in one of the circuits but is inactive in the other. Because of the symmetry between the two circuits, there are two such states. They differ in terms of the circuit which has an active cortical population. The activity of the cortex in the circuit where it is active is as follows:

where T_Ctx is the threshold of the cortical neurons and I₀ is given by Equation 31 in the Appendix. These two states also differ in term of the activity in the other populations. In particular, activity in the GPi is smaller in the circuit in which the cortex is active than in the circuit in which it is inactive. As a consequence, activity in the thalamus is greater in the former circuit than in the latter. The stability boundary of these latter states, given by Equation 23, is represented by the dotted line.

On the dashed line, the linear steady state undergoes an oscillatory Hopf bifurcation driven by the negative feedback on the hyperdirect loop. This instability to oscillations depends on ₊, _–, and µ. In the special case in which ₊ = _– and µ = 1, the corresponding line in the phase diagram is straight. This is not the case in general. For instance, if ₊ > _–, or if µ > 1, this line is shifted toward larger values of G_– and it is convex as shown in Figure 4.

When G₊ = 1 + G_–, A_k^Ctx in Equation 25 diverges. This indicates that the asymmetric state becomes unstable on this line (Fig. 4, dotted line). In fact, in the region G₊ > 1 + G_– (below the dotted line), a detailed analysis shows that depending on the external input, H^Ctx, one, two, or three stable steady states exist. In all of these states, the two feedback loops are open. If H^Ctx < H_max1^Ctx, only one steady state exists in which the cortical populations are inactive in the two circuits. If H_max1^Ctx < H^Ctx < H_rest^Ctx, this steady state still exists, but it coexists with two other states in which the cortical population is active in one of the two circuits and the GPi population is silent, whereas in the other circuit, the cortical population is silent and the GPi population is active. Thus, for these parameters, the network displays tristability. If H_rest^Ctx < H^Ctx, only two asymmetric states coexist.

The phase diagram of Figure 4 was obtained with the values of the synaptic delays given in Table 1. These parameters are sufficient to reproduce recent experimental results on response latencies in BGs after cortical stimulation (Nambu et al., 2000) (see below). In particular, the latency from the STN to the GPi is smaller than from the Str to GPi, and the overall latency is shorter in the hyperdirect pathway than in the direct pathway (₊ = 26 ms; _– = 20 ms). We also studied the phase diagram as a function of the parameters. We found that the location of the symmetry-breaking bifurcation does not depend on the delays, but the location of the bifurcation to oscillations does. However, the phase diagram remains qualitatively the same in a range of delays compatible with the data reported by Nambu et al. (2000) (data not shown).

View this table: [in this window] [in a new window]   Table 1. Values of the various parameters of the model

 
Response of the network to an external transiently selective input in the various regimens
With the network parameters given in Table 1 and in the absence of external input, the average firing rates of the GPi, thalamus, STN, striatum, and cortex are 80, 20, 25, 0, and 0 spikes/s, respectively. Moreover, the negative feedback of the hyperdirect loop generates oscillations if it is not compensated for by a sufficiently strong positive feedback in the direct loop. Hence, as the strength of the corticostriatal interactions decreases, the parameter G₊ decreases and the network moves in the phase diagram from the symmetry-breaking regimens, to the linear regimen and eventually to the oscillatory regimen.

In this section, we describe how the network responds in these various regimens to a time-dependent input that is nonselective on the cortex ( = 0) and transiently and slightly selective on the striatum. This input is shown in Figure 5A. The nonselective input to the cortex and the slightly selective input to the striatum start at the same time. The input to the striatum has a duration of 200 ms.

View larger version (24K): [in this window] [in a new window]   Figure 5. The responses of the reduced network model to a weakly and transiently selective external input in the various regimens of the phase diagram. Network parameters are given in Table 1. Activities in and inputs to circuit 1 (resp. circuit 2) are plotted with solid lines (resp. dashed lines). A, The external input. Top, The input to the cortical populations in the two circuits. Bottom, The input to the striatum. The responses of the cortex (red), the thalamus (black), and the GPi (blue) are plotted in B, the symmetry-breaking regimen (GStrCtx = 0.7); C, the linear regimen (GStrCtx = 0.4); and D, the oscillatory regimen (GStrCtx = 0.05). Note that only one cortical population is activated in the symmetry-breaking regimen (B), whereas both populations are activated in the other regimens. E, The response of the network to a strong cortical input in the oscillatory regimen (GStrCtx = 0.05). The striatal input and the profile of the cortical input are as in A. The amplitude of the latter has been increased by a factor of 3. The oscillation is suppressed because the thalamus is silent, and subsequently the feedback loops are open. F, The response of the network to the input displayed in A in the multistability regimen (GStrCtx = 0.9). As in the symmetry-breaking regimen, one cortical population is activated, whereas the other is silent.

 
Selectivity in the symmetry-breaking regimen. The response of the network in this regimen is shown in Figure 5B. At the input onset, the cortical populations in the two circuits are activated in the same way. Hence, initially, the striatal neurons in both circuits receive the same excitation from the cortex. However, the transient input to the striatum, which is slightly larger in circuit 1, induces a slight difference between the activities of the striatum in the two circuits, with the striatum in circuit 1 being more active. Just like the brief external input in Figure 3A, this transient striatal input induces a symmetry breaking in the system. As a result, the cortex becomes more active in circuit 1 than in circuit 2. After the transient component of the external input is over, the difference between the activities in the two circuits persists and even amplifies. This is the outcome of the symmetry-breaking instability that drives the network toward an asymmetric state in which the cortical population in one of the circuits is quiescent. Therefore, in this regimen, a transient weakly biased input to the striatum is sufficient to induce a strongly asymmetric response of the cortex even if the two cortical populations are activated exactly the same way. In this sense, the cortex–BG network is able to perform action selection.

Linear regimen. The network is in this regimen if the cortex–striatum synaptic efficacies are not strong enough to achieve symmetry breaking but are sufficiently strong to compensate for the oscillations induced by the negative feedback of the hyperdirect loop. In this regimen, the nonselective input to the cortex leads to coactivation of both cortical populations (Fig. 5C). The GPi neurons also increase their activity in relation to this additional input. This is because the excitatory input on the GPi coming via the hyperdirect pathway increases more than the inhibitory input coming via the direct pathway. The striatal selective input induces a transient, small asymmetry between the activities in the two circuits. This asymmetry is proportional to H^Str (data not shown). However, once the striatal input is over, the symmetry between the two circuits is rapidly restored, and the activities in the cortex are the same in the two circuits. Hence, in this regimen, once the transient input to the striatum is over, the response of the cortex is qualitatively similar to what would happen in the absence of the cortex–BG–thalamus network. The effect of the latter is only to change the value of the gain of the effective input–output relation of the cortical neurons.

Oscillatory regimen. If the weight of the corticostriatal synapses is too small, the direct loop cannot prevent the oscillations from appearing. Hence, when the cortex is activated, the network tends to develop oscillatory activity. As in the linear regimen, the input to the striatum induces a transient asymmetry in the response of the two circuits, which rapidly decays once the input to the striatum is over. The pattern of the response is shown in Figure 5D. Interestingly, if H^Ctx is sufficiently strong, the oscillations are suppressed because the activity of the GPi becomes so strong that it completely inhibits the thalamus. The feedback is then suppressed in the hyperdirect loop and the system cannot display oscillations. This is shown in Figure 5E.

The multistability regimen. The response of the network in the multistability regimen is shown in Figure 5F. The external input induces an asymmetric response of the network because the symmetric steady state is unstable. One cortical population is activated, whereas the other is inhibited, as in the symmetry-breaking regimen. However, in contrast to what happens in this latter regimen, the strong positive feedback sustains the activity in the active cortical population and the network remains in an asymmetric state even after removal of the external input. Hence, the network displays multistability between the rest state and the two asymmetric states. This regimen can be relevant for modeling the physiology of the substantia nigra pars reticulata, where selective delayed activity is found (Hikosaka et al., 2000). However, a detailed study of the behavior of the network in this region of the phase diagram is beyond the scope of the present study.

The physiological regimen. Symmetry breaking occurs in the whole region of the phase diagram below the solid line. This region can be divided into two parts. Between the solid and the dotted line, the network is monostable. Below the dotted line, it displays multistability. In this regimen, the positive feedback sustains the activity and prevents the network to return to rest state after the removal of the external input, inducing persistent activity. We interpret the first of these regions as the one in which the system (which represents the motor part of the BG–cortical network) performs normally. The fact that the size of this region is reduced means that normal function can only be achieved for a right balance between the hyperdirect and direct feedbacks. We will show in the framework of our more detailed model that dopamine depletion disturbs this balance, because it reduces G₊, leaving G_– constant. As a consequence, selectivity is impaired, because the DA level is reduced. If G_– is large, oscillations can emerge for sufficiently strong DA depletion. Interestingly, the range of G₊ for which the system behaves "normally" is broader for larger G_–. Hence, the propensity of the network to display pathological oscillations increases with the robustness of its functional state.

The detailed model
A major asset of the reduced model is that its phase diagram can be derived analytically. This allows us to show how the BG network can perform action selection and how this ability is lost if the feedback on the direct loop is too weak. We also show that oscillatory activity can emerge in the BG network if the feedback in hyperdirect pathway is strong and the feedback in the direct pathway does not compensate for it.

However, our reduced model is too simplified to account for the resting properties of the BG network, because in this model, in the absence of external input the cortex is not active. This prompted us to investigate the more realistic model of the cortex–basal–ganglia network described above (see Materials and Methods, Numerical simulations of the detailed model). In this model, the connectivity is random, the thresholds of striatal neurons are distributed, and noise is present in all of the populations. The parameters (synaptic strength, delays, thresholds) of the model are given in Table 1. Note that the synaptic strength, G, and the thresholds of the neurons have been modified compared with their values in the reduced model. With these synaptic coupling values, the feedback in the hyperdirect loop is strong enough to induce oscillatory activities driven by this loop (see below).

The resting state in the normal situation (D = 100%)
Neuronal activities at rest for D = 100%. We define the rest state of the network as the state in the absence of external inputs. In contrast to the reduced model, in the detailed model, the neurons in the cortical and striatal populations display some low but nonzero level of activity at rest. This is because noise is present in their input. For the chosen parameters, the average activity of the cortical neurons is 5 spikes/s. In the GPi, STN, thalamus, and striatum, the activities are 80, 20, 25, and 0.6 spikes/s, respectively. In the GPi and striatum, the average activities are broadly distributed as shown in Figure 6A. In particular, the distribution of the activity in the striatum is very skewed, with one-half of the neurons completely quiescent.

View larger version (27K): [in this window] [in a new window]   Figure 6. Neuronal activities in the network at rest. A, The distributions of the average spontaneous firing rates of striatal and GPi neurons. B, Population activities of groups of 20 neurons randomly chosen in the GPi (blue lines), the thalamus (black lines), the STN (yellow lines), the cortex (red lines), and the striatum (cyan lines). C, Correlation matrix of three GPi neurons. Diagonal, Autocorrelograms. Off-diagonal, Cross-correlations. D, The effective input–output transfer function of a neuron in the striatum in the absence (gray line) and in the presence of noise (black line, SD of the noise given in Table 1). In the presence of noise and for low input, the effective gain of striatal neurons is reduced.

 
In this rest state, there is no symmetry breaking and the average and distribution of activity in populations of the same type are the same in the two circuits. This is because the activity in the striatal populations is very low. In fact, as shown in Figure 6D, the effective transfer function of the striatal neurons in the presence of noise is nonlinear, and the effective gain of the striatal neurons is decreased for low activities. Hence, the effective feedback of the direct loop is too small for symmetry breaking to occur in the rest state. Nevertheless, in the normal rest state, the feedback in the direct loop is sufficiently strong to compensate for the oscillatory instability driven by the negative feedback in the hyperdirect loop. This is clear from Figure 6, B and C, which shows that the activities in cortex, thalamus, and GPi do not oscillate and are not synchronous across neurons.

Response of the neurons to brief cortical stimulations. The emergence of oscillations and their frequency depend on the delays in the various pathways (for values, see Table 1). Therefore, it is important to verify that the values we have selected for these parameters are compatible with experimental data.

One way to characterize these delays is by looking at the dynamics of the activity changes in the different populations following a brief and sufficiently strong stimulation in the cortex. This approach has been used in several experimental studies (Deniau et al., 1978; Maurice et al., 1998; Nambu et al., 2000). In particular, Nambu et al. (2000) showed that a cortical stimulation induces triphasic response in the GPi with an early excitation followed by a suppression and by a late excitation. They argued that the early excitation was mediated by the hyperdirect pathway and that the suppression was mediated by the direct pathway. They also argued that the hyperdirect pathway contributes to the late excitation, because inactivation of the GPe did not completely eliminate this phase of the response.

To compare the behavior of our model with the results of Nambu et al. (2000), we studied the response of the network to a brief (1 ms duration) stimulation in the cortical population in one of the circuits. The population response of a group of 20 neurons in the GPi is shown in Figure 7. It displays a triphasic pattern (Fig. 7, top graph) with an early excitation followed, 10 ms