Global and local excitation and inhibition shape the dynamics of the cortico-striatal-thalamo-cortical pathway

The cortico-striatal-thalamo-cortical (CSTC) pathway is a brain circuit that controls movement execution, habit formation and reward. Hyperactivity in the CSTC pathway is involved in obsessive compulsive disorder (OCD), a neuropsychiatric disorder characterized by the execution of repetitive involuntary movements. The striatum shapes the activity of the CSTC pathway through the coordinated activation of two classes of medium spiny neurons (MSNs) expressing D1 or D2 dopamine receptors. The exact mechanisms by which balanced excitation/inhibition (E/I) of these cells controls the network dynamics of the CSTC pathway remain unclear. Here we use non-linear modeling of neuronal activity and bifurcation theory to investigate how global and local changes in E/I of MSNs regulate the activity of the CSTC pathway. Our findings indicate that a global and proportionate increase in E/I pushes the system to states of generalized hyper-activity throughout the entire CSTC pathway. Certain disproportionate changes in global E/I trigger network oscillations. Local changes in the E/I of MSNs generate specific oscillatory behaviors in MSNs and in the CSTC pathway. These findings indicate that subtle changes in the relative strength of E/I of MSNs can powerfully control the network dynamics of the CSTC pathway in ways that are not easily predicted by its synaptic connections.

in specific nodes and cell types shapes the network dynamics of the CSTC pathway. Here we address this question using a theoretical approach, based on a mathematical model of the excitatory (glutamatergic) and inhibitory (GABAergic) connections in the CSTC pathway. By using coupled Wilson-Cowan models and bifurcation theory 16 , we determine how global and local changes in excitation and inhibition contribute to orchestrate the activity of D1-and D2-MSNs and of the entire CSTC pathway.

Materials and Methods
We use a system of seven nonlinear equations to describe the major excitatory (glutamatergic) and inhibitory (GABAergic) synaptic connections of the CSTC pathway, thought to be equivalent in rodents and humans of either sex or age ( Fig. 1) 2 . The general structure of this circuit is conserved across rodents and humans of either sex and age ( Fig. 1) 2 , allowing our findings to be applicable to a variety of contexts and animal species. In this system, each equation provides a quantitative description of the strength of excitation (c e ) and inhibition (c i ) onto each of the seven nodes of the CSTC pathway (Eqs 3-10). We then examine the network dynamics of the CSTC pathway using a modeling approach based on the use of analytic equations originally developed by Wilson and Cowan 16 and extensively used in a number of modeling contexts [17][18][19] . According to this model, the level of activity in each node (representing the mean activity of a population of neurons) can be described as: Figure 1. Synaptic connections of the CSTC pathway in humans and rodents. (A) Anatomical distribution of the synaptic connections of the CSTC pathway in the human brain. The CSTC pathway is a multi-synaptic neuronal circuit that connects the cortex, striatum and thalamus. The striatum receives glutamatergic inputs (green) from the cortex and the thalamus and sends out GABAergic (red) inputs to the sub-thalamic nucleus. (B) Anatomical distribution of the synaptic connections of the CSTC pathway in the rodent brain. (C) Schematic representations of the synaptic connections in the CSTC pathway. The two most abundant cell types in the striatum, D1-and D2-MSN, receive glutamatergic inputs from the cortex and the thalamus and inhibit each other via GABAergic synaptic connections. D1-and D2-MSNs send GABAergic signals to other basal ganglia nuclei via direct and indirect connections to the substantia nigra (SNr), respectively. Abbreviations: GPe (globus pallidus external part), GPi (globus pallidus internal part), STN (sub-thalamic nucleus), SNs (substantia nigra pars compacta). where X represents the mean level of activity in the cortex (C), striatum (S), thalamus (T), external and internal capsules of the globus pallidus (E and I, respectively), D1 and D2-MSNs (D1 and D2, respectively). We define = ∑ + ∑ Z cI cI e e i i as the weighted sum of all network inputs received by a node X. Here I e and I i represent the activity in all nodes from which X receives excitatory/inhibitory inputs, whereas c e and c i represent average synaptic strengths. The original Wilson-Cowan model obtains theoretical conditions for the synaptic parameters that produce transitions in behavior and then interprets them in a physiological context. Then specific ranges are used to illustrate these conditions and the corresponding bifurcations. Since, for the coupled dynamics, one would expect the relevant dynamic behavior to be encountered in similar parameter ranges as in the original model, we use an extended range that includes all values used in the Wilson-Cowan paper and in subsequent work 16,20 to describe the synaptic strength values c e and c i (0 ≤ c e , c i ≤ 40). In this context, values of c e = 20, c i = 20 are considered to be indicative of a normal, physiological state. The function S θ,b expresses the integrated response of a node to the input at time t and has sigmoidal shape, which represents the typical physiological neural response to stimuli of increasing strength. Accordingly, a gradual increase in input strength elicits first a low firing response (initial portion of the sigmoidal curve), followed by a window of stimulus strengths that induce profound changes in the firing response (steep portion of the sigmoidal curve), and by a region where the firing rate does not change any more (asymptotic portion of the sigmoidal curve). Equation 1 also incorporates the time constant (identical for all nodes) and the activity history of the node/population. As in the Wilson-Cowan model, in each node, the state X = 0 corresponds to a state of low-level background neuronal activity, with small negative values of X representing depression of resting activity. Hence the state in which X = 0 for all nodes must be a stable steady-state solution for our equations in the absence of external inputs, in order to have physiological significance. In the Wilson-Cowan model, this requirement is fulfilled by setting S θ,b (0) = 0, therefore considering it of the form: The sigmoidal parameters θ and b give the value and position of the maximum slope of the sigmoid function in each node. Varying these parameters controls the position and impact of the steep portion of the curve, but the results of the model remain qualitatively similar regardless of their specific value (in terms of number and stability of steady states, hysteresis effects, presence of limit cycles, etc.) as long as the general shape of the sigmoidal curve remains the same. Therefore, in this study, we fix the sigmoidal parameters describing excitation (θ e = 4, b e = 1.2) and inhibition (θ i = 2, b i = 1) to values that are within the range of those originally considered by Wilson 16 ), and consider them as indicative of the physiological control state. Therefore, the state describing the physiological ongoing activity in different nodes of the CSTC pathway is defined as the state where: c e = 20, The model exhibits dynamic phenomena typical of a coupled system of non-linear oscillators, which are described in more detail in the Results section. The time constants of all nodes are identical (i.e. all nodes are synchronized, with no aperiodic behavior/chaos). The strength of this modeling framework is that it approaches simultaneously the temporal behavior of seven coupled nodes of the CSTC pathway. However, because the model represents a seven dimensional phase space combined with a multi-dimensional parameter space, it is difficult to perform a global analysis of the system's dynamics in both phase and parameter spaces. To overcome this potential limitation, one can either perform computer-assisted searches for all types of behaviors that the system can achieve 21 , or focus on specific behaviors of interest, which is the approach that we use here. By constructing a sequence of progressively refined models, we can tease apart the effects of gradually finer differences in E/I onto D1-and D2-MSNs in regulating the network dynamics of the CSTC pathway without having to vary more than two key parameters at the same time (i.e. we work in a 2D parameter space). We then shift our attention to the phase space, to study the effects of specific perturbations on the activity of selected nodes.
The mathematical (numerical) analysis of the network dynamics of the system (i.e. the CSTC pathway) is performed as follows. First, for fixed parameter values, we analyze the phase space to identify the geometry and position of asymptotic attractors, such as stable equilibria and cycles. The basin of attraction of the stable equilibria is determined by the range of values of the seven-dimensional initial conditions that evolve asymptotically towards that particular attractor. The position of the asymptotic attractors describes the activity of the nodes when they have reached their long-term stable states. Some initial conditions may converge to stable equilibria with higher activity for D1-MSNs than for D2-NSNs. Other initial conditions may converge to stable equilibria with higher activity for D2-MSNs than for D1-MSNs. Some initial conditions may converge to states of perpetual oscillations with characteristic amplitude and duty cycle. Second, we analyze how the phase space, with the position and geometry of its attractors, evolves in response to changes in E/I, which we refer to as perturbations. Third, we perform a numerical search for bifurcations, which are critical states where the system exhibits sharp transitions between qualitatively different phase space dynamics (e.g. from a non-oscillatory to oscillatory activity). We find co-dimension one bifurcations by studying the sensitivity of the system to changes of one parameter at a time. We find co-dimension two bifurcations by simultaneously changing two parameters at a time. In this study, we specifically analyze how the system responds to independent perturbations of either E or I (and search for co-dimension one bifurcations) or to concurrent perturbations of both E and I (and search for co-dimension two bifurcations). The perturbations are performed on excitatory and inhibitory connections of the entire CSTC pathway or only on those targeting D1-or D2-MSNs.
We use the Matcont package 22,23 to search specifically for four types of co-dimension one bifurcations: branch points (BPs, or transcritical bifurcations), limit points (LPs, or saddle node bifurcations), Hopf bifurcations (H) and limit point cycles (LPCs, or fold bifurcations). We use this nomenclature in the remainder of the paper for consistency with the Matcont package terminology. BPs occur when two equilibrium curves intersect and simply swap stability. LPs occur when two equilibria collide and vanish. When one of these equilibria is attracting, a drastic rearrangement of the phase space is necessary so that a different attractor may take over the basin of the lost stable equilibrium (hysteresis). In our system, LPs may occur when a perturbation causes a change from high activity state for D1-MSNs to a high activity state for D2-MSNs, even when the system is initiated at the same original state for both variables. H bifurcations occur when an equilibrium suddenly changes stability, giving rise to a limit cycle. We are particularly interested in supercritical Hopf bifurcations, where a stable equilibrium may give birth to a stable cycle that acts as a new asymptotic attractor. At a supercritical H bifurcation, the system may be redirected from a long-term stable state (e.g. high activity for D1-MSNs) to a regime of perpetual oscillations (e.g. between high and low D1-MSN activity). LPCs occur when two cycles collide and disappear (similarly to when two equilibria collide at an LP). If either cycle is an attractor for a set of initial conditions, its loss requires redirecting these solutions towards a different attractor (e.g. a stable equilibrium). An LPC may be, for that reason, a plausible mechanism for cessation of oscillations by small parameter perturbations. In addition to asymptotic attractors, we also search for saddle equilibria and cycles, which only attract solutions originating in lower-dimensional basins ("shallow" basins). For example, at a supercritical Hopf bifurcation, a saddle equilibrium may lose two attracting dimensions to give birth to a saddle cycle with an attracting basin contained in a two-dimensional sub-space of the original seven-variable space. The presence of shallow basins may strongly influence the geometry of other attractors and their own basins, affecting a wide array of trajectories. Shallow basins are hard to identify by randomly choosing initial conditions that eventually converge to a specific saddle point/cycle. In our numerical investigations, we identify saddle equilibria and cycles and their attracting basins whenever possible and we specify when computations become intractable. We consider the system to be in a state of multistability if multiple local attractors co-exist in the phase space, for fixed parameter values, situation which leads to initial conditions that evolve towards different long-term dynamics, depending on the basin of which attractor they belong to. Initial conditions close to the boundary between two basins can easily be pushed towards either attractor by very small perturbations in the initial state of the system (without even changing the system parameters). The ability of the system to swiftly converge to a different steady state from slightly different original states supports the idea that the long-term behavior of the network relies on the E/I balance acting in conjunction with the current state of the system. As we identify all of these different behaviors, we point out their relevance for our study. Our models are not independent, but rather developed incrementally, with increasingly subtler phenomena brought under focus with each section.

Results
Effect of network-wide changes in E/I. To provide a mathematical representation of the network activity of the CSTC pathway ( Fig. 1), we use the following system of seven nonlinear equations (see Methods): (1 ) (1 ) ( ) We first assume that all inputs mediating either excitation or inhibition behave identically throughout the entire CSTC pathway and therefore have the same sigmoidal parameters (excitation: θ e = 4 and b e = 1.2; inhibition: θ i = 2 and b i = 1) and strength (excitation: c e ; inhibition: c i ). We ask how changing the strength of E/I throughout the entire system affects its network dynamics. To address this question, we first set c e = 8 and c i = 10 ( Fig. 2, bottom left). We then vary c e and c i progressively until c e = 30 and c i = 40 (Fig. 2, top right). When c e and c i are small, the system initiated at zero stabilizes to an asymptotic attractor of low cortical, striatal and thalamic activity. In contrast, when c e and c i are large, the system converges to a state of high activity for the cortex, striatum and thalamus (Fig. 2, top right). The transition between these two distinct states with small or large (c e , c i ) values can occur through different paths, along which c e is either larger (Fig. 2, bottom right) or smaller than c i (Fig. 2, top left). These paths can cross regions of the parameter space where the system generates oscillations, with amplitude and duty cycles that depend on the specific values of c e and c i (Fig. 3). That is to say that the equilibrium curves describing the activity of the cortex, D1-and D2-MSNs and the thalamus show different profiles depending on whether c e is varied while keeping c i constant ( Fig. 3A-D) or whether c i is varied while keeping c e constant ( Fig. 3E-H). When the system reaches a supercritical Hopf bifurcation point, it gains access to a regime of periodic oscillations. In our simulations, network oscillations that arise at the Hopf bifurcation point disappear via a fold bifurcation. Varying only c e for a specified fixed value of c i , like varying only c i for a specified value of c e , does not significantly affect the position of the Hopf bifurcation point (Fig. 3). In contrast, increasing both c e and c i alters more efficiently the position of the Hopf bifurcation point and therefore changes the dynamics of the system more profoundly (Fig. 3I). Values of (c e , c i ) that bring the system inside the Hopf loop (blue area in the (c e , c i ) parameter plane illustrated in Fig. 3I) allow the system to show oscillations. Values of (c e , c i ) that bring the system out of the Hopf loop, prevent it from showing stable oscillations (white and green areas in Fig. 3I). The results of this set of simulations show that global changes in either excitation or inhibition control the steady state level of activity of the CSTC pathway. In contrast, concurrent changes in global E/I control the stability of the system, allowing it to show network oscillations with specific amplitudes and duty cycles.

Effect of local changes in E/I onto D1-and D2 -MSNs.
We next ask how local changes in E/I specifically targeted onto D1-and D2-MSNs affect the network dynamics of the CSTC pathway. To address this question, we add two new parameters to the system of seven nonlinear equations, which we name ′ c e and ′ c i . We begin our analysis by altering the strength of excitation onto D1-and D2-MSNs. To do so, we change ′ c e , while keeping ′ c i fixed (Section and Fig. 4A-D). Then, we alter the strength of inhibition onto D1-and D2-MSNs. This is accomplished by changing ′ c i , while keeping ′ c e fixed (Section and Fig. 4E-H). As in the previous model, we analyze the evolution of the system's asymptotic behavior in the cortex, D1-and D2-MSNs and the thalamus.     The most notable feature in this set of simulations is that the system shows a bi-stable behavior along the entire range of tested values for ′ c e (Fig. 4). When ′ c e is small ( < ′ < c 5 10 e ), the system has two stable equilibrium curves (blue and green curves) and one saddle equilibrium curve (orange curve). If all variables are initiated at rest (zero initial conditions), the system converges to the orange equilibrium; if the rest initial conditions are slightly perturbed in the D1-MSN component, then the system converges to the blue equilibrium (Fig. 4B); if they are perturbed in the D2-MSN component, then the system converges to the the green equilibrium (Fig. 4C). The position of all these equilibrium curves remains approximately constant, allowing the system to remain either in a stable state of high activity for the cortex, D1-MSNs and thalamus and low activity for D2-MSNs, or in a stable state of low activity for the cortex, D1-MSNs and thalamus and high activity for D2-MSNs. As ′ c e increases ( < ′ < c 10 1 5 e ), the stable (green) and saddle (orange) branches do not change appreciably, while the stable (blue) branch increases steadily in its D2-MSNs component. At the BP, when ′ ∼ c 15 e , the upper part of the pink branch takes over the stability from the blue branch, and maintains, for values of ′ c e higher than BP, the asymptotic state of high activity in all four regions of the system, attracting all initial conditions with higher activity of D1-versus D2-MSNs. A decrease in the initial activity of D1-versus D2-MSNs can lead the system to converge towards the other stable state (green), of low activity for the cortex, D1-MSNs and thalamus and high activity for D2-MSNs . Hence different types of unbalanced activity of D1-versus D2-MSNs can bring the system into one or the other attraction basins. This new bistable asymptotic regime survives without state transitions for all values of ′ > c 15 e . A perfect balance in the activity of D1-and D2-MSNs allows the system to reach a saddle equilibrium curve (orange curve) for all values of ′ c e . This saddle equilibrium curve corresponds to a low activity state in the cortex, striatum and thalamus, which persists until the Hopf bifurcation point, where the system gains a saddle cycle.

Effect of local changes in inhibition onto D1-and D2-MSNs.
We use the following system of equations to determine the effect of changing the strength of local inhibition onto D1-and D2-MSNs (i.e. ′ c i ).
(1 ) (1 ) ( ) A first notable feature in these simulations is that multiple equilibrium curves remain similar along the entire range of values for ′ c i (Fig. 4E-H). When < ′ < c 0 7 i , the system converges to a steady state characterized by high activity in the cortex, D1-and D2-MSNs and thalamus (pink curve). A second notable feature is the presence of stable limit cycles. When ′ ∼ c 7 i , the system undergoes a fold bifurcation and generates a stable cycle (not shown in Fig. 4), which later collides with the blue equilbrium curve. More precisely, this cycle coexists with the blue equilibrium branch within a narrow window of ′ c i values ( < ′ < c 7 9 i ). Then it disappears through a Hopf bifurcation when ′ ∼ c 9 i , changing stability of the blue equilibrium curve from a saddle to an attracting equilibrium. When ′ c i increases beyond the Hopf bifurcation value ( ′ > c 9 i ), the system reaches a regime with two stable equilibria (blue and pink curves). One stable branch that emerges from the Hopf point (blue curve) corresponds to a state of high activity in the cortex, thalamus and D2-MSNs, coupled with low activity in D1-MSNs. The second ScieNtific REpoRTS | 7: 7608 | DOI:10.1038/s41598-017-07527-8 branch (pink curve) corresponds to state of high activity in all brain regions, which preserves high activity levels for both D1-and D2-MSNs until the following bifurcation point ( ′ = c 26 i ). Here the blue branch loses stability and leads to a dramatic drop in in the activity of D2-MSNs. When ′ > c 26 i , the D2-MSN component of the blue equilibrium curve decays rapidly, while all other system components retain high activity levels. These findings indicate that changing the strength of local inhibition onto D1-and D2-MSNs allows these cells to show not only concurrent or non-concurrent high activity levels, but also stable oscillations between these states.

Effect of local changes in the relative E/I of D1-versus D2 -MSNs.
We next ask how the activity of the system changes when altering the relative excitation or the relative inhibition of D1-and D2-MSNs. To do this, we keep constant the excitation and inhibition in the rest of the CSTC pathway (c e = c i ) and then introduce new variables to vary either: (1) the relative strength of excitation onto D1-MSNs (c e1 ) and D2-MSNs (c e2 ) or (2) the relative strength of inhibition onto D1-MSNs (c i1 ) and D2-MSNs (c i2 ). In both cases, we work with a two-dimensional parameter space, as described in the Methods section.
Effect of local changes in the relative excitation of D1-versus D2-MSNs. We use the following system of equations to determine the effect of changing the relative strength of excitation onto D1-and D2-MSNs. and the orange curve that remains constant for the whole c e 1 range (Fig. 5A,B). All other branches that we identified numerically are unstable saddles, and the supercritical Hopf bifurcation marked along the pink curve creates saddle cycles with shallow basins (not captured in the subsequent phase planes). Phase plane plots confirm that, when = c 10 (co-existence of this existing node with a new stable node of high activity for D1-MSNs; Fig. 5D). As mentioned above, this transition occurs through an LP bifurcation (Fig. 5A,B). Similarly, when setting = c 20 e 2 , the system has access to a stable steady state with high D2-MSNs activity levels throughout the range of c e 1 (Fig. 5E,F). The system exhibits a similar phase transition between mono-and bi-stability, but in this case the transition appears for larger values of c e 1 ( ∼ c 15 ). When < c 15 e 2 , the only possible activity level of D1-MSNs is low (Fig. 5G). When > c 15 e 2 , the system gains access to an alternative stable state where both D1-and D2-MSNs have high activity levels (Fig. 5H). Taken together, these findings indicate that increasing excitation of D1-MSNs generates bistability irrespective of the levels of excitation of D2-MSNs. When excitation of D2-MSNs is low, the activity level of D1-MSNs is opposite to the activity level of D2-MSNs: high activity in D1-MSNs causes low activity in D2-MSNs and vice versa (Fig. 5C,D). When excitation of D2-MSNs is high, the activity level of D2-MSNs remains high regardless of the activity level of D1-MSNs. Therefore, high levels of excitation of D2-MSNs can be coupled with either high or low activity levels of D1-MSNs (Fig. 5G,H).

Effect of local changes in the relative inhibition of D1-versus D2 -MSNs.
We use the following system to study the effect of altering the relative strength of inhibition onto D1-and D2-MSNs:  . We let c i 1 vary between = c 0 The bifurcation diagrams show hysteresis regardless of the specific value of c i 2 (Fig. 6A,B,F,G). We first considered the case of = c 7 i 2 ( Fig. 6A,B). For ). In fact, this portion of the equilibrium curve is almost identical to its counterpart for the case of = c 7 However, in our current case, the second stable regime emerges sooner, when > .  ), in the sense that the saddle equilibria and cycle that coexist with these attractors have shallow attraction basins, and are thus unlikely to be reached by generic initial conditions. The phase diagrams in Fig. 6C-E,H-J illustrate some of these behaviors in more detail for = c 7 i 2 , in a (D1, D2) phase space slice. When < . c 7 58 i 1 , the system has a unique attracting equilibrium with high D1-and D2-MSN activity (Fig. 6C,D). When = .
c 7 58 , it crosses the LPC bifurcation, leading to the birth of a stable cycle and the onset of network oscillations which co-exists with the original attracting equilibrium (Fig. 6E). At the Hopf bifurcation point ( ∼ . H 10 15), the cycle collapses into a low D1-and D2-MSN activity equilibrium, which continues to coexist with the original high D1-and D2-MSN equilibrium (Fig. 6H). Through a short hysteretic window, flanked by LP = 19.97 and LP = 20.77, a new low equilibrium of low D1-and D2-MSN activity is born (Fig. 6I) and replaces the pre-existing low D1-and D2-MSN activity equilibrium (Fig. 6J).
In order to illustrate more comprehensively the simultaneous dependence of the system's behavior on both inhibition levels c i 1 and c i 2 , we plot bifurcation curves in the c c ( , ) i i 1 2 parameter plane (Fig. 7A). We focus in a region around a Zero-Hopf co-dimension two bifurcation point (ZH), where three bifurcation curves intersect, separating three parameter regions that identify three distinct dynamic behaviors, illustrated in Fig. 7B-D. The green region, analyzed in Fig. 7B, is characterized by the presence of two stable equilibria, one with high D1-and D2-MSN activity, the other with low D1-and D2-MSN activity. The blue region, analyzed in Fig. 7C, is characterized by the formation of a stable cycle when entering from the green region by crossing the H (green) curve. The cycle vanishes via a fold bifurcation when crossing the LPC (blue) curve, so that in the pink region only the stable equilibrium of high D1-and D2-MSN activity survives (Fig. 7D). Finally, one can re-enter the green region by crossing the LP curve (pink), with formation of two new equilibria (a node and a saddle), bringing the system back to the initial bistability (Fig. 7B).
Taken together, these findings indicate that increasing inhibition onto D1-MSNs when inhibition onto D2-MSNs is low introduces a bistable behavior where both cell types show either high or low levels of activity. In contrast, increasing inhibition onto D1-MSNs when inhibition onto D2-MSNs is high allows to change the activity level only in D1-MSNs, while the activity of D2-MSNs remains high. In the previous section we saw that altering the relative excitation of D1-and D2-MSNs allows the system to transition through dynamics regimes with saddle limit cycles, which attract only a shallow subspace of initial conditions. Here we show that altering the relative inhibition of D1-and D2-MSNs pushes the system through dynamics windows with stable cycles with attraction basins of significant size. Therefore, altering the relative inhibition of D1-and D2-MSNs provides a more powerful mechanism than altering the relative excitation of these cells, which can push the system into a regime of stable oscillations. Effect of local and coupled changes in the E/I of D1-and D2 -MSNs. In the previous sections, we considered simplified scenarios where we changed either excitation or inhibition onto D1-MSNs, while keeping excitation or inhibition onto D2-MSNs constant at arbitrary values. Here we consider the more complex scenario in which the two effects are coupled. First, we change the magnitude of either local excitation or local inhibition proportionately onto D1-and D2-MSNs (using the parameter h). Second, we scale these effects differently in D1and D2-MSNs (using the parameters a and b).
Effect of local and coupled changes in the excitation of D1-and D2 -MSNs. We use the following system to study the effect of a coupled modulation of excitation onto D1-and D2-MSNs: (1 ) (1 ) ( ) (1 ) ( ) . We analyze these effects at two different values of a. By increasing a, we scale the effects induced by varying h differently in D1-and D2-MSNs.
Here we analyze two scenarios in which we vary h: when a = 1 (Fig. 8A-D) and when a > 1 (Fig. 8E-H). We plot a series of bifurcation curves obtaine by varying h over an extended (and partly unbiological) range (−1 < h < 2), to visualize the origin of the bifurcation phenomena. Consistent with our previous model (Fig. 5), increasing h when a = 1 moves the system between two bi-stability windows, where a low D1-and high D2-MSN activity state (green curve) coexists with a high D1-and low D2-MSN activity state (orange curve; Fig. 8A-D). This is to be expected, since this model introduces only a re-parametrization of the excitation strengths from the c c ( , ) e e 1 2 parameter plane to the new (h, a) parameter plane. Similar findings are observed for a larger scaling factor a = 1.2 (Fig. 8E-H). These findings suggest that local changes in excitation onto D1-and D2-MSNs can alter the level of activity in these cells regarless of the fact that they may be more or less pronounced in one of the two cell populations.

Effect of local and coupled changes in the inhibition of D1-and D2 -MSNs.
We use a similar approach to determine the network effects of local changes in inhibition. To do this, we use the following system of equations: . When b = 0.8, progressively increasing h allows the system to transition from a bi-stable regime of high D1-and D2-MSN activity coexisting with low D1-and high D2-MSN activity, to a regime of high D1-and D2-MSN stable activity, via an LP bifurcation (Fig. 9A-D). When b = 1, changes in h allow the system to transition between the same two regimes described for b = 0.8 (i.e. bistability when h is small, mono-stability when h is large; Fig. 9E-J). In addition, when b = 1, increasing h can trigger oscillations via two supercritical Hopf bifurcations (Fig. 9F,G). The first Hopf point (at ∼ . h 0 04) gives birth to a saddle cycle (Fig. 9E) which survives until ∼ .
h 0 17 and disappears via a LPC. This saddle cycle remains symmetric for the whole duration of its life (i.e. D1 = D2), and has an attraction basin that also lies in the sub-space D1 = D2 (Fig. 9E). The second Hopf point (at h = 0.82) moves the system from bi-stability (Fig. 9H) into stable oscillations (Fig. 9I), which disappear at ∼ .
h 0 92 via a LPC. This cycle attractor coexists throughout its life with a high D1and D2-MSN stable equilibrium (Fig. 9I), before disappearing for larger values of h (Fig. 9J). Setting b = 1.2 generates qualitatively similar results, only in a more efficient way. Therefore, under these conditions, a progressively larger reduction of inhibition onto D1-than D2-MSNs evoked by increasing h allows the system to transition from bi-stability to small stable oscillations to a single high-activity equilibrium (Fig. 9K-O).
Effect of local and coupled changes in the E/I of D1-and D2 -MSNs. As a last step in our analysis, we ask how simultaneous, coupled changes of both excitation and inhibition onto D1-and D2-MSNs can alter the activity level of these cells, using the following system of equations:   Figure 10 shows how the system can be pushed through a sequence of bifurcations and a succession of regimes by increasing the parameter h, for two instances of (a, b), shown side by side. In Fig. 10A-F, a = 1 and b = 1. In Fig. 10G-L, a = 1 and b = 1.2. Both sets of models show qualitatively similar results. Accordingly, for small values of h, the system resides in a bi-stable regime of differentially high/low activity for D1-and D2-MSNs. As h increases, the system shows stable out-of-phase oscillations until both D1-and D2-MSNs show high levels of activity ( Fig. 10A-F). Increasing b from b = 1 to b = 1.2 does not change the occurrence of these transitions, but causes subtler, yet important effects ( Fig. 10G-L). First, for the higher b, the transition between the bi-stability and mono-stability regimes for increasing values of h is interrupted by two LPs, which cause an abrupt transition in the equilibrium curves. During this transition, for 0.25 < h < 0.35, the system goes through a small window with a unique equilibrium, characterized by low activity in both D1 and D2-MSNs. This suggests that when inhibition decreases more in D1than D2-MSNs, the system passes through a transient state of basal D1-and D2-MSN activity before reaching a mono-stability regime. Second, the stable oscillation window between the LPC and the H bifurcation points is prolonged. This is important, since the presence of oscillations prevents the system from converging to a state of steady hyperactivity for D1-and D2-MSNs.

Discussion
The goal of this study is to detemine how global and local changes in excitation and inhibition in distinct domains of the CSTC pathway affect the activity of this brain circuit and of specific classes of striatal neurons. The two major pieces of evidence that motivate our study are: (1) the CSTC pathway controls important physiological functions, like movement execution; (2) the CSTC pathway is a major site of synaptic dysfunction in neurospychiatric disorders like OCD and Tourette's syndrome, which are characterized by the repeated execution of involuntary movements 2,24 . In this pathway, cortical and thalamic neurons projecting to the striatum form glutamatergic synapses onto D1-and D2-MSNs. These cells form local inhibitory GABAergic synapses onto each other and send long-projection GABAergic afferents to the substantia nigra pars reticulata (SNr) and the internal capsule of the globus pallidus (GPi) via two anatomically distinct pathways (the direct and indirect pathways, respectively). Information from D1-and D2-MSNs is ultimately relayed to the thalamus and the motor cortex 25 . Coordinated activation of D1-and D2-MSNs is crucial for action selection and inhibition of unwanted behaviors as well as to control the overall level of activity in the CSTC pathway 26,27 . Accordingly, reduced activity of D1-MSNs and increased activity of D2-MSNs generate a powerful inhibition of the motor cortex. In contrast, increased activity of D1-MSNs and reduced activity of D2-MSNs promote motor cortex activity and movement execution, which is thought to contribute to the motor symptoms of OCD [28][29][30][31] . Consistent with these findings, clinical sudies in humans indicate that patients with OCD show functional hyperactivity in the CSTC pathway. These findings are largely based on the use of positron emission tomography (PET) and single photon emission CT (SPECT) [32][33][34][35][36][37] and on resting-state functional magnetic resonance imaging (rs-fMRI) [38][39][40][41] , all of which lack single-cell resolution. A recent connectomics study based on the combined use of magnetic resonance scans and deterministic fiber tracking based on diffusion tensor imaging data indicates that these functional alterations in the CSTC pathway are associated with structural alterations in the connectivity strength between the orbito-frontal cortex and the striatum, two important nodes of the CSTC pathway 42 . What is not known is whether these functional and structural abnormalities in the CSTC pathway can be attributed to changes in excitatory and/or inhibitory transmission onto MSNs.
In this work, we model the activity of the CSTC pathway using a system of seven non-linear equations, to provide a mathematical description of the currently known wiring diagram of the CSTC pathway. Our modeling approach, originally proposed by Wilson and Cowan 16 , allows us to describe the mean field activity in neural populations, not the firing rate of individual cells (consistently with the level of resolution of neuroimaging works). According to this definition, a given population firing rate can be due to stereotyped firing rate in all neurons or to the presence of a small number of high firing neurons. The spontaneous activity of medium spiny neurons in the striatum involves recurring shifts in membrane potential from a hyperpolarized (down) state to a depolarized (up) state accompanied by irregular spikes at low frequency and burst firing [43][44][45][46] . This is different from the activity . of other types of neurons (e.g. dopaminergic neurons in the ventral tegmental area of the midbrain 47 , which exhibit ongoing (tonic) firing activity interleaved by brief periods of bursting (phasic) activity 48 . In this context, the Wilson-Cowan type model used here allows us to accurately capture the temporal evolution of the mean level of activity of distinct neuronal populations without adding further information on the specific firing activity of individual cells.
Although some works suggest that glutamatergic synaptic dysfunction contributes to the etiology of OCD 12 , others indicate that changes in GABAergic inhibition disrupt the activity of the CSTC pathway and lead to the onset of OCD-like behaviors 11,[13][14][15] . Our network analysis sheds light on this conundrum. It indicates that changes in excitation or inhibition can both cause an increased activity in D1-and reduced activity in D2-MSNs (see red cells in Fig. 11). In particular, this pattern of activity can be generated by: (1) increasing excitation on both D1and D2-MSNs (Fig. 4); (2) varying excitation onto D1-MSNs at low levels of excitation of D2-MSNs (Figs 5 and 8); (3) decreasing inhibition of D1-MSNs when excitation of both MSNs is low (Fig. 10). Varying GABAergic inhibition onto MSNs exerts a more powerful control on the oscillatory activity (than the overal level of activity) of D1 and D2-MSNs. Power spectra analysis from striatal local field potentials and electrocorticograms highlight the presence of two major types of spontaneous oscillatory activities in the low (2-9 Hz; δ, θ) and high frequency range (30-80 Hz; γ) 49 . A dynamic progression in the dominant mode of oscillatory activity is thought to occur in the striatum when learning (ventral striatum) and when performing a specific motor repertoire (dorsal striatum) [50][51][52] . According to these findings, the oscillatory activity in the ventral striatum evolves from transient bouts in the γ range to transient bouts in the β range (15-28 Hz) as a specific movement is consolidated into a habit through learning 52 . The ventral striatum is thought to direct the activity of the dorsal striatum through corresponding changes in its oscillatory activity 50 . For this reason, changes in the relative inhibition of D1-and D2-MSNs provide a valuable candidate mechanism that contributes to disrupted habit formation and movement execution in OCD 5, 53-60 .
In our analysis, we point out that there are wide parameter ranges that confer bi-stability to the system. In particular, the asymptotic trajectory of the system depends crucially on the positions of its initial state. In most cases, the hyperspace D1 = D2 delimits a boundary between two different attracting regimes: one with high D1-and low D2-MSN activity, a second with low D1-and high D2-MSN activity (Figs 5,8,9 and 10). A small initial bias in the activity of either D1-or D2-MSNs would lead the system to one of these attractors. In many instances, the hyperspace D1 = D2 contains a saddle equilibrium or a saddle cycle. Only initial conditions in which the activity of D1-MSNs is identical to that of D2-MSNs would maintain the system in this hyperspace. However, regimes with perfectly identical D1-and D2-MSN activity are extremely unlikely to be reached, although other factors that have not included in this model (e.g. spike timing) could make this regime more likely to be maintained.
Taken together, these findings allow us to analyze the effects of small perturbations on the activity of D1-and D2-MSNs and of the entire CSTC pathway. They highlight that local inhibition, among all other factors, can powerfully regulate the activity of the entire system, possibly contributing to the hyperactivity patterns observed in patients with OCD.