Visual Abstract
Abstract
The striatum is the main input nucleus of the basal ganglia. Characterizing striatal activity dynamics is crucial to understanding mechanisms underlying action selection, initiation, and execution. Here, we studied the effects of spatial network connectivity on the spatiotemporal structure of striatal activity. We show that a striatal network with nonmonotonically changing distance-dependent connectivity (according to a gamma distribution) can exhibit a wide repertoire of spatiotemporal dynamics, ranging from spatially homogeneous, asynchronous-irregular (AI) activity to a state with stable, spatially localized activity bumps, as in “winner-take-all” (WTA) dynamics. Among these regimes, the unstable activity bumps [transition activity (TA)] regime closely resembles the experimentally observed spatiotemporal activity dynamics and neuronal assemblies in the striatum. In contrast, striatal networks with monotonically decreasing distance-dependent connectivity (in a Gaussian fashion) can exhibit only an AI state. Thus, given the observation of spatially compact neuronal clusters in the striatum, our model suggests that recurrent connectivity among striatal projection neurons should vary nonmonotonically. In brain disorders such as Parkinson’s disease, increased cortical inputs and high striatal firing rates are associated with a reduction in stimulus sensitivity. Consistent with this, our model suggests that strong cortical inputs drive the striatum to a WTA state, leading to low stimulus sensitivity and high variability. In contrast, the AI and TA states show high stimulus sensitivity and reliability. Thus, based on these results, we propose that in a healthy state the striatum operates in a AI/TA state and that lack of dopamine pushes it into a WTA state.
Significance Statement
Recent findings suggest that striatal activity is organized in spatially compact neuron clusters. Here, we show that striatal projection neurons should have a nonmonotonically changing distance-dependent connectivity to obtain spatially localized activity patterns in striatum. Among the different states a striatal network can show, asynchronous-irregular and transition activity states closely resemble striatal activity in the healthy state. In contrast, strong cortical inputs as observed in Parkinson’s disease drive the network into a winner-take-all state, in which the striatum loses its stimulus sensitivity. Thus, our model makes specific predictions about the spatial network connectivity in the striatum and provides new insights about how the striatum might make a transition from a healthy state to a Parkinson’s disease state.
Introduction
The basal ganglia (BG) are a collection of nuclei located at the base of the forebrain. BG are crucial for critical functions such as voluntary movement control, decision-making, procedural learning, and a variety of cognitive functions. The striatum as the main input stage of the BG plays a major role in BG-related brain function and dysfunction. Like many other subcortical structures, the striatum is a network of inhibitory neurons driven by excitatory inputs from multiple regions of the neocortex (Wall et al., 2013), thalamus (Smith et al., 2009a), and hippocampus. The output of the striatum directly projects to the globus pallidus (GP) and forms two major functional pathways in the BG, the “go and “no-go” pathways. The interaction of these pathways forms the basis of functions such as action selection, initiation, and execution (Albin et al., 1989).
Experimental data suggest that striatal firing rates in ongoing activity are quite low (≤1 Hz; Mahon et al., 2004; Gage et al., 2010), and ≈30% of striatal projection neurons respond with a markedly increased firing rate (≈10–20 Hz; Kimchi and Laubach, 2009; Gage et al., 2010; Adler et al., 2012). However, it remains unclear how the sparse and low-firing-rate task-related activity of the striatum is organized in space and time to shape the activity of the GP to initiate action selection.
Recent observations of temporal (Adler et al., 2013) and spatial (Barbera et al., 2016) clustering of spiking activity in the striatum suggest that striatal activity may be organized in transient coactivations of groups of neurons. Such coactivated neuron groups can be termed neuronal assemblies (NAs). Such NAs can effectively increase the inhibitory influence of striatal neurons on downstream networks. However, conditions under which NAs appear in the striatum are not well understood.
Within the framework of dynamical systems, the NA-type activity state could be an outcome of a “winnerless competition” (WLC; Rabinovich et al., 2001) or a noise driven transition state between the asynchronous-irregular (AI) activity state and the “winner-take-all” (WTA) state. Although the mechanisms underlying the emergence of such transient NAs are not well understood, heterogeneity in both the external input and the recurrent connectivity among striatal projection neurons is a common feature of the existing models of transient NAs (Humphries et al., 2009; Ponzi and Wickens, 2010; Angulo-Garcia et al., 2015).
Here, we investigate the existence of NAs in a large-scale network model of the striatum in which neurons are connected in a distance-dependent manner. We address two key questions: (1) How does the structure of the network connectivity define the repertoire of the ongoing activity dynamics? and (2) Under which conditions can task-related inputs modify these ongoing activity dynamics?
Specifically, we considered two qualitatively different spatial connectivity profiles in which the connection probability between any two neurons varied as a function of the distance between the neurons: decreased monotonically (according to a Gaussian distribution) or changed nonmonotonically (according to a gamma distribution; Rinzel, 1998). Using network simulations and mathematical analysis, we show that spatially structured activity in the network can emerge only for a nonmonotonically shaped connectivity kernel. In inhibitory networks with a monotonically changing connectivity kernel, network activity remains AI. Only networks with a nonmonotonically changing connectivity kernel can exhibit different activity regimes. Weak background inputs or high input variance induced unstructured AI activity, whereas stronger background inputs or low input variance induced stable dynamics (e.g., WTA). For moderately strong inputs, the network activity organized into unstable spatial bumps [transition activity (TA)] resembling the experimentally observed NAs (Adler et al., 2013; Barbera et al., 2016). Finally, we describe how stimulus response characteristics are affected by the ongoing activity state of the network. We show that only AI and TA states support a reliable stimulus response and that stimulus-induced changes in pairwise correlations are maximal in the TA state, especially for weak stimuli.
These results help us to better understand information processing in the striatum and other primarily inhibitory subcortical networks as in the amygdala and the cerebellum. Moreover, our model provides new insights into neuronal mechanisms underlying brain diseases such as Parkinson’s disease (PD). For instance, the strong resemblance between the properties of the WTA state and the activity observed in PD (Costa et al., 2006) leads to the hypothesis that PD may be a transition from AI or TA to a WTA state.
Materials and Methods
Neuron model
Neurons in the network were modeled as “leaky-integrate-and-fire” (LIF) neurons. We intentionally chose this simple neuron model to ensure that the observed network dynamics would be the result of the connectivity profiles studied and not the intrinsic dynamics of a more complex neuron model. The subthreshold membrane potential dynamics of LIF neurons are given by
(1)
Synapses were modeled as conductance transients. The temporal profile of the transient conductance change in response to a single presynaptic spike was modeled as an alpha function. The recurrent connectivity within an inhibitory network such as the striatum is weak and sparse (Tepper et al., 2004; Taverna et al., 2008; Planert et al., 2010). We adjusted the synaptic conductances to obtain weak synapses such that a unitary inhibitory postsynaptic potential (IPSP) had an amplitude of 0.8 mV at a holding potential of –44.0 mV and an excitatory postsynaptic potential (EPSP) had an amplitude of 1.6 mV at a holding potential of –70.0 mV. The various neuron and synapse parameters are summarized in Table 1. Whenever possible, we used parameters corresponding to the striatal network (Yim et al., 2011).
Parameter values for the neuron and synapse model
Network architecture
We considered a population of 10,000 inhibitory neurons. The neurons were placed on a 100 × 100 grid and folded as a torus to avoid boundary effects. The distance between neighboring nodes in the grid network amounted to 10 µm. Each neuron sent 1000 inhibitory recurrent inputs to other neurons within the network (i.e., 10% total connection probability), mimicking the sparse connectivity in most biological inhibitory networks. We implemented no self-connections within the network, and neurons were allowed to form multiple connections between them.
To implement a distance-dependent connectivity, we chose two qualitatively different spatial profiles (Rinzel, 1998). For the first type of connectivity, we assumed that the distance-dependent connection probability decreased monotonically as a function of distance. To implement such connectivity, we used the Gaussian distribution to model the distance-dependent decrease in connectivity (Fig. 1a, top). We will refer to this network type as a Gaussian network.
Spiking activity of the network with gamma-distributed (top) and Gaussian-distributed (bottom) connectivity. (a) Connection probability between neurons as a function of the distance between neurons normalized to the full size of the network. (b) Examples of the spiking activity of 200 of 10,000 neurons for different average firing rates (
) of the network. With increasing background activity, the spiking activity of the Gaussian networks remained irregular and the bursting behavior increased. In contrast, the gamma networks showed transient or persistent bursting behavior and local synchronization of spiking activities.
For the second type of connectivity, we assumed that the distance-dependent connection probability varied nonmonotonically as a function of distance. Recent experimental data suggest that the connectivity between pairs of medium spiny neurons (MSNs) in the striatum may indeed change nonmonotonically as a function of distance (Fujiyama et al., 2011; López-Huerta et al., 2013). To implement such connectivity, we used a gamma distribution to model the distance dependent variation in connectivity (Fig. 1a, bottom). We will refer to this network type as a gamma network.
The out-degree of individual neurons was fixed to 1000 neurons. In our network, target neurons were located on a uniformly spaced grid (100 × 100), and each neuron index can be calculated by its coordinate on the grid map. Moreover, each target neuron could also be identified by an angle and radius in polar coordinates:
where x and y are the coordinates of the grid locations, idx is the index of the target neurons, and ncol = 100 is the size of the network in one of two dimensions. To identify 1000 targets, we drew ϕ from a uniform distribution 
and r from either Gaussian or gamma distribution. This approach is much faster than deciding connectivity based on pairwise distances, but can be used only when neurons are arranged on a grid.
To appropriately compare these different network types, we ensured that the connection weights and the degree distributions were identical in the two network types.
Ongoing/background excitatory drive
All neurons received background excitatory inputs mimicking the global ongoing activity in striatum. The external drive was modeled by uncorrelated homogeneous Poisson-type spike trains with fixed firing rate. Each neuron received an independent realization of such spike train, obtained from the same underlying Poisson process. The rate of the Poisson-type spike trains was varied systematically to study the different dynamical states of the inhibitory networks.
Stimulus
To measure the impact of the an external stimulus on network activity in each dynamic state, we stimulated the network with an external stimulus, on top of the background inputs described above. The stimulus was provided to a subset of ≈45 neurons. To distribute the stimulated neurons in a spatial manner, we defined a region of interest (ROI) of size 30 × 30 neurons. In the ROI, we defined a location as center for stimulated neurons chosen according to a two-dimensional (2D) Gaussian profile (σ: 2 grid points). The external stimulus was modeled as an injection of direct current to the selected neurons. The stimulus presentation lasted for 1 s, and its amplitude was systematically varied between 0 and 150 pA. To collect sufficient data for further statistical analysis, each stimulus was presented 20 times to its assigned subset of neurons.
To estimate how the response magnitude and trial-by-trial variability changes when one stimulus is immediately followed by another, we extended the simulation of the network with two different external stimuli. Each of these two external stimuli was provided to two nonoverlapping subsets of ≈45 neurons. In the ROI, as described above, we defined two locations that were 10 grid points apart. Using these two locations as centers, stimulated neurons were chosen according to a 2D Gaussian profile (σ: 2 grid points). The two stimuli were presented in alternating fashion, and each stimulus presentation lasted for 1 s, separated by a pause of 2 s. Thus, the two stimuli differed only in terms of their target neurons; other properties (strength, duration) were identical.
Characterization of the network dynamics
The network activity states were characterized on the basis of the spiking activity using standard descriptors: firing rates and the coefficient of variation of interspike intervals (CVISI), the latter used to measure the degree of irregularity of the spiking of individual neurons Kumar et al. (2008b).
Count and duration of activity bumps
We constructed a series of binary matrices representing the 2D map (100 × 100 neurons) of the network activity over disjunct time windows of 100 ms each from the spike train recordings of all neurons. To enhance the spatial contrast, these matrices were convolved with the 2D “Mexican hat”–shaped kernel (Ricker wavelet). The size of the hat was taken from the estimation of the bump size (Figs. 4b and 5). The integral of the Mexican hat kernel was set to zero to ensure that the filtered activity was comparable for different inputs. Additionally, the integrals of the positive and the negative part of the convolution kernel were set to 1 and –1, respectively. After convolving the spatial activity of the network (measured in a 100 ms window) with the Mexican hat kernel, we thresholded the resulting image to identify patches of high activity. The center of each patch was taken as the position of individual bumps. The positions of individual bumps were then used to compute the number of bumps and determine the affiliation of neurons to individual bumps. Finally, trackers of time-merged bumps were used to analyze the lifespan of the bumps, determined by finding similar positions of bump centroids over subsequent time frames.
Characterization of the stimulus response
To estimate the influence of the external stimuli on striatal ongoing activity, we analyzed the activity of stimulated neurons in the ROI (of size 30 × 30 grid). To characterize the impact of the external stimulus on the network activity dynamics, we measured the average firing rate, the trial-by-trial variability of the network response, and the spectrum of pairwise correlations of the network spiking activity. To estimate whether the evoked response of the stimulated neurons was discernible from the background, we measured
for each stimulus. The trial-by-trial variability was quantified using the Fano factor (FF), defined as the ratio of the variance and the mean of the response of the stimulated neurons.
The correlation spectrum was estimated for the activity of individual neurons. Neuronal spiking activity was binned using a rectangular bin of 100 ms. We measured the correlation spectrum for the trial-by-trial average of both ongoing and stimulus-evoked activity.
Simulation tools
All simulations of the networks of spiking neurons were performed using NEST simulation software (Gewaltig and Diesmann, 2007). The dynamical equations were integrated at a fixed temporal resolution of 0.1 ms. Simulation data were analyzed with Python using the scientific libraries SciPy (http://www.scipy.org) and NumPy (http://www.numpy.org) and visualized using the plotting library Matplotlib (http://matplotlib.org).
Results
Neuronal activity in the striatum is maintained by thalamic and cortical excitatory inputs. In a randomly connected inhibitory network model of the striatum, the recurrent inhibition and the level of external excitatory inputs define the dynamical states and stimulus response properties of the network (Brunel and Hakim, 1999; Ponzi and Wickens, 2010; Yim et al., 2011; Angulo-Garcia et al., 2015). Here, we extend this earlier work by investigating the effects of spatial connectivity on the dynamical states and the stimulus-response properties of the striatal network activity. To this end, we used both analytical calculations using neural field equations and numerical simulations of biologically realistic, large-scale inhibitory networks models with 10,000 spiking neurons.
In our network models, the profile of the spatial connection probability between any two neurons could vary either monotonically or nonmonotonically as a function of the distance between the neurons. For this, we used two different kernels for the spatial connectivity: In Gamma networks, the nonmonotonic distance-dependent change in connection probability was modeled as a Gamma distribution [off-center inhibition (Rinzel, 1998); Fig. 1a, top], whereas in Gaussian networks, the connection probability decreased monotonically in a Gaussian manner as a function of distance between neurons [on-center inhibition (Rinzel, 1998); Fig. 1a, bottom].
Firing rates and spike train irregularity
First, we determined the firing rate and the degree of irregularity of the spiking activities in the networks connected according to the Gaussian-shaped (Figs. 1a, b, top) and gamma-shaped (Figs. 1a, b, bottom) connectivity profiles. As expected, the average firing rates in both network types increased as a function of the background excitatory input (νext; Fig. 1b). However, comparing the network activity at similar average output rate (
), the activity pattern of network differed in connectivity profiles (Fig. 1b). Gamma networks were clearly more sensitive to a change in the background input rate (νext) than Gaussian networks (Fig. 2a, black trace). Thus, the same range of average firing rates as in Gaussian networks was achieved in gamma networks already for a much lower range of external inputs (νext; Fig. 2a, black trace). This suggests that a gamma-shaped spatial connectivity kernel increased the sensitivity of the network activity to changes in the background input.
Analysis of the spike patterns in networks with different connectivity profiles. a, Distribution of average firing rates (
; top) and coefficient of variation (CVISI; bottom) as a function of background input strength (νext) for Gaussian (left) and gamma (right) networks. In Gaussian networks, increasing νext resulted in a steadily widening distribution of 
and CVISI, which for a large fraction of neurons tended to the value 3. In contrast, in gamma networks, the distribution of 
and CVISI was rapidly widely distributed from a relatively low νext (1.5 kHz) upward. Gamma networks were clearly more excitable than Gaussian networks. Green trace indicates the skewness of the firing rate in both networks. Black and blue traces refer to the average firing rate (
) and standard deviation (σλ) of the firing rate distribution, respectively. b, Relationship between the irregularity of the spiking pattern (CVISI) and the average firing rate (
). The color of the traces represents the background input rate (νext). In Gaussian networks (top), neurons with higher 
tended to exhibit a higher CVISI, whereas in gamma networks (bottom), they tended to exhibit a lower CVISI.
Moreover, these two network types differed considerably in the distribution of firing rates over neurons (Fig. 2a, top row) and spike timing irregularity (CVISI; Fig. 2a, bottom row). We compared the output firing rate distributions for a range of external inputs that resulted in the same range of output firing rates. The recurrent inhibition in the network ensured that even for very strong background input, only few neurons could spike at a high rate, and a major fraction of neurons spiked at relatively low rates, resulting in a positively skewed distribution of firing rates. For Gaussian networks, the average firing rate (
) and the standard deviation (σλ) and skewness (κλ) of the neuronal firing rate distribution monotonically increased as background input was increased (Fig. 2a, top right). In contrast, in Gamma networks, both 
and σλ increased nonmonotonically, but the skewness of the firing rate distribution changed in a nonmonotonic fashion as a function of background input rate (Fig. 2a, top right). This suggests that a nonmonotonically changing connectivity profile results in a larger diversity of firing rates in the network and that the diversity of firing rates is maximized for a moderate amount of background input.
Interestingly, in Gaussian networks, the neurons’ CVISI monotonically increased with an increase in their firing rate. An average CVISI > 1 indicates that most high-firing-rate neurons spiked in bursts (Fig. 2b, top). Because we used LIF neurons, the bursting in the network activity was caused by the recurrent inhibition and not intrinsic neuron properties. In the Gaussian networks, the connection probability peaked over only a small range in the vicinity of a given neuron, thereby reducing the effective recurrent inhibition and allowing the neuron to maintain its high firing rate for some time. However, fluctuations in the external input (modeled as Poisson spike trains) could rapidly switch the high-rate activity from any one neuron to another, thereby creating a spike pattern consisting of short-lived bursts and pauses.
In contrast, in the Gamma networks, only neurons with moderate output firing rates showed a CVISI > 1 (Fig. 2b, top). Neurons with very high firing rates (
≥100 Hz) spiked in a regular manner (CVISI ≤0.5). Such a small CVISI implies that the network was operating in the so-called mean-driven regime (Brunel, 2000), in which the background input is strong enough to keep the free-membrane potential of the neurons above spiking threshold. Within the physiologic range of output firing rates (
≤50 Hz), neurons in both network types elicited spike bursts, but the CVISI in Gaussian networks was nearly twice as high as that in Gamma networks.
Spatial structure of the network activity
Next, we included spatial information about neurons in our analysis and characterized the spatial activity patterns in both network types. Visual inspection of the spike rasters (Fig. 1b) suggested that in Gaussian networks, the structure of the network activity remained spatially homogeneous, even for very high background input rates (νext). In contrast, when increasing the background input rate in gamma networks, the activity got confined into local clusters, resulting in spatially periodic activity bumps. The spatial structure became more apparent when we rendered the neuronal activity as a 2D surface (Figs. 3 and 4a). Depending on the strength of the background input, three qualitatively different network activity states were observed. In the AI state, neurons spiked at a low rate and the activity was more or less homogeneously distributed across the network (Fig. 3, top row). For very strong background inputs, the network activity organized into a spatially periodic and temporally stable bump structure (Fig. 3, bottom row). This state resembles the k-WTA state (Hutt and Atay, 2005). In between, for moderate inputs, the bump structure was aperiodic and unstable (Fig. 3, middle row). We refer to this state as transition activity (TA).
Characterization of the bump activity states for the gamma network. Time series snapshots of the 2D pattern (100 × 100 neurons) of bump activity, contrast-enhanced by Mexican hat filtering. Each frame was measured by summing the neuronal activity over 100 ms. Three snapshots (columns) of activity were taken 1 s apart. Three representative gamma network states for three different amounts of external inputs are shown in each row. AI: external input = 1 kHz; no bump activity is observed and the network activity remains noisy. TA: external input = 1.5 kHz; the network is in an unstable state, with several bumps appearing transiently, in the company of noisy activity. WTA: external input = 3 kHz; the network forms mostly persistent bumps throughout the entire network.
To better characterize the emergence of activity bumps and these three dynamic states in gamma networks, we measured the numbers, sizes, and lifespans of activity bumps. To reliably identify bumps, we designed a bump detection algorithm (see Methods). First, we measured the average firing rate as a function of distance from the center of a bump. For weak inputs, the average network activity decayed and reached a baseline level with increasing distance from the center of the bump. However, for stronger inputs, not only did the firing rate in the bump increase, but the spatially periodic nature of the activity also became more apparent (Fig. 4a).
Quantification of bump activities in gamma networks. a, Spatial autocorrelation of the network spiking activity, showing the mean firing rate (
) of each bump as a function of the distance from the cluster centroid normalized to the full size of the network. Different colors represent the strength of the background excitation. b, Distance between bumps in various gamma distributions and its comparison between numerical simulations and mean field equations. The background color is measured bump distance from the network simulation data. Solid traces show results from analytical estimation of bump distance; dashed traces represent the estimation of bump distance from network simulations. The parameter for the gamma-distributed connection is used for the spiking network model (red circle). Spatial representations of activity bumps are also observed in Fig. 5 (c–e). These subfigures display the bump count and their relative lifespan during the entire simulation as a function of the background external excitation (νext), which modulates the nature of the bump activity. c, With increasing νext, the number of bumps increases in a sigmoidal fashion. For higher νext, the number of bumps saturates, owing to the limited capacity of the finite spatial map. The error bars for the bump counts indicate the standard deviation of bump counts over the simulation time. d, The lifespan of bumps reflects the dynamic state of the bumps: a shorter lifespan reflects TA dynamics, whereas a larger lifespan indicates stable bump activity reflecting WTA dynamics. By increasing νext, the distribution of lifespans shifts from short to long term. The ordinate indicates lifespan, normalized to the duration of the entire simulation (10 s), of individual bumps. Because the average bump counts are different in each dynamical state (TA, WTA), we normalized the color bar of bump count to the average bump count in individual states. e, The lifespan distribution is split into three groups (dashed lines in d). The long (red trace) appearance of bumps reflects the WTA state of bump activity, whereas the short (blue trace) appearance of bumps reflects the TA state of bump activity. Between these two states, the network is in a highly unstable state, characterized by a wider distribution of lifespans (light blue trace).
The number of bumps increased in a sigmoidal fashion with the background input rate (Fig. 4c). For strong inputs (νext ≥1.7 Kspikes/s), the bumps were stable: neurons that started to spike at the beginning of the simulation remained active for the entire simulation (Fig. 4d, e, red trace). For weak inputs (νext ≤1.2 Kspikes/s), the few bumps that appeared lasted only for short durations (Fig. 4d, e, blue trace). For the medium range of input rates (1.2 ≤ νext ≤ 1.7 Kspikes/s), the bumps showed a wide range of lifespans (Fig. 4d, e, light blue trace).
This analysis of the network activity showed that gamma networks exhibited random unstructured activity for weak background inputs and stable periodic bump activity for strong background inputs. The stable bump activity was similar to the periodic bump activity observed in networks of excitatory and inhibitory neurons interconnected according to a spatial Gaussian connectivity profile (Mehring et al., 2003; Roxin et al., 2005; Kumar et al., 2008a). This raises the question of why inhibitory networks with a gamma-type connectivity profile exhibit spatial bump activity patterns and why networks with Gaussian-distributed connectivity profile do not.
Necessary conditions for the emergence of bump states
To understand why a gamma-type connectivity profile exhibits spatial bump activity and a Gaussian connectivity profile does not, we investigated the dynamical states of spatially interconnected inhibitory neurons using neural field equations (Ermentrout, 1992; Hutt and Atay, 2005). For simplicity, we started by formulating the neural field equation for a one-dimensional network with circular boundary conditions. The results hold also for a 2D network with torus-folded boundary conditions because we are considering a homogeneous and isotropic scenario. The mean membrane potential v(x,t) in the continuum limit is given by
(2)
where W(x) denotes the spatial connectivity profile and f(v) the neuron transfer function, mapping the mean membrane potential to the output firing rate (
). Without loss of generality, we used normalized connectivity profiles W(x), and the scaling parameter was chosen such that the absolute connection strength was absorbed into f(v). The background input to the neurons is denoted by I(x,t). For the network, the stationary and spatially homogeneous solution v(x,t) = v0 = constant is a solution to Eq. (2) for a constant background input I(x,t) = I0. Here, v0 is given by the implicit equation
(3)
To investigate the stability of this homogeneous solution, we considered small perturbations 
around v0 and linearized the transfer function in Eq. (2). After subtracting Eq. (3), this yields
(4)
Here and in the remainder 
is used as shorthand notation for 
. In the Fourier domain, this expression simplifies to
(5)with 
denoting the Fourier transform of ν with respect to space. We can now obtain the eigenvalues:
(6)
When the eigenvalues are positive, small perturbations do not die out, indicating unstable dynamics. Assuming that the slope of the transfer function 
is always nonnegative, negative values in 
are a necessary condition for positive eigenvalues ζ and, hence, for spatially periodic activity bumps. In purely inhibitory networks, this condition can be fulfilled by off-center inhibition kernels, such as the gamma-kernel under investigation, or, e.g., a mixture of two Gaussian distributions arranged symmetrically around zero. Although for biologically plausible connection kernels, off-center inhibition (i.e., a nonmonotonic kernel) fulfills this condition, it should be noted that this is not a necessary condition. For instance, the Fourier transform of a box-shaped kernel around zero takes negative values at nonzero frequencies and therefore could, in principle, generate bump states.
For the gamma distribution–shaped connectivity kernel,
(7)with shape parameter n and scale parameter Θ, this condition of positive eigenvalue ζ > 0 can be fulfilled for n > 1, when the Fourier transform 
takes negative values (Fig. 6, right). In contrast, the Gaussian kernel Fourier-transforms into another Gaussian, which never takes negative values (Fig. 6, right), and therefore, purely inhibitory networks with a Gaussian-type connectivity profile do not show any spatially periodic bump activity. When the eigenvalues ζ are negative, the network activity remains spatially homogeneous, similar to the AI state observed in both the Gaussian networks and the gamma networks (Fig. 1c).
Eq. (6) also revealed that the slope of the transfer function at stationary rate 
is an important factor controlling the eigenvalues. For positive eigenvalues, it needs to be sufficiently large. That is, the following condition needs to be fulfilled:
(8)Here, 
is controlled by a number of factors. For instance, an increase in the synaptic weights (e.g., cortico-striatal, MSNs to MSNs, and fast spiking interneurons to MSNs) increases 
and can make the network cross the bifurcation. Similarly, increasing the background input rate will also increase the slope 
, because the transfer function of LIF neurons in the simulated network is convex for low firing rates (Burkitt, 2006). That is, both an increase in the strength of recurrent synapses and an increase in the background input rate can cause the transition from spatially homogeneous firing to periodically organized bump solutions.
This analysis shows that the activity of an inhibitory network driven by constant input I(x) has two stable solutions: for weak inputs, the stationary state is spatially homogeneous (AI state; Fig. 3, top row), whereas for strong inputs, the stationary state is spatially periodic (WTA state; Fig. 3, bottom row). When the network is driven with noisy inputs or the neurons have unequal numbers of synapses (owing to random connectivity), the transition between the two stable solutions is smoothened (TA state; Fig. 3, middle row). This result of the analytical calculations is consistent with the numerical network simulations (Figs. 1 and 4) where noise was introduced into the network by the connectivity and the Poisson type spike trains as external background input.
In addition to stating the condition for spatially periodic bump solutions to arise, Eq. (6) enables us to estimate the distance between bumps. The wave number of the emerging spatially periodic solution is approximately given by the wavenumber kc minimizing the Fourier transform of the gamma kernel (7). This critical wavenumber kc is given by
(9)
This means that the spacing between bumps will increase as the shape (n) and scale parameter (Θ) of the gamma function are increased (Fig. 4b, dashed traces). We confirmed this in numerical simulations by systematically varying the two parameters of the connectivity kernel (Fig. 5). We found that, indeed, the analytical approximation closely matched the estimates of the interbump distance measured in network simulations (Fig. 4b).
Spatial map of bump activity patterns for different gamma distributions. A snapshot, contrast-enhanced by Mexican hat filtering, of the 2D pattern (100 × 100 neurons) of bump activity for different parameters (shape, scale) of gamma connection profile defines the size of bumps and the distance between bumps. The rate of the external Poisson excitation (νext) was set to 5 kHz to obtain WTA states in networks with different gamma distributions.
Impact of the ongoing activity state on the stimulus response
The three different networks states we identified were obtained by changing the global input to the network, and therefore, these states can be thought as possible ongoing states of the network activity, i.e., when the animal is not engaged in a specific task. To understand how a task-related activity may be represented in the striatum when it operates in one of these three ongoing activity states, we selectively stimulated two different subsets of neurons located in spatially nonoverlapping regions in the network alternately (stimulus A and B; Fig. 7a; see Methods) and measured the stimulus response.
In the AI state, the network instantaneously responded to the stimulus and switched to a different spatial pattern as the input was changed (Fig. 7a, top row). In this state, there were no bumps in the ongoing activity; therefore the stimulus created new bumps. In the TA as the ongoing state (Fig. 7a, middle row), external stimuli modified the already existing bump pattern. The resulting patterns were different for the two stimuli. By contrast, in the WTA state, the ongoing activity already showed strong stable bumps and the external stimuli proved insufficient to alter the ongoing bump pattern (Fig. 7a, bottom row), unless it overlapped with an existing bump (e.g., the stimulus B in Fig. 7a, bottom row).
To quantify the striatal response in different ongoing activity states, we measured the change in the average firing rate of the stimulated neurons (Δresponse, Fig. 7b) and the trial-by-trial variability (measured over 20 trials) for the two stimuli (Fig. 7c). In AI and TA states, the external input elicited a strong response, well discernible from the background activity. As expected, the response magnitude increased with an increase in the stimulus strength (Fig. 7b). The variance of the evoked response (Δresponse measured across trials for the whole duration of the stimulus) in the TA state was higher than in the AI state, but for stronger stimuli both average response and trial-by-trial variability were similar (Fig. 7b). The time-resolved trial-by-trial variability of the evoked response (measured as FF) of the network response was maximal at stimulus onset (Fig. 7c). However, the trial-by-trial variability of the response in the steady state (in the presence of the stimulus) was smaller than that observed in the ongoing state. By contrast, in the WTA state, not only the Δresponse was smallest among the three states, it was also more variable across trials (depending on the location of the stimulus). For strong input it was possible to elicit a strong reliable response even in the WTA state, but that response depended strongly on the location of the stimulus, e.g., among the two stimuli we tested, only stimulus B resulted in a high Δresponse (Fig. 7b, bottom row). Besides, in the WTA state, trial-by-trial variability also depended on the stimulus location. For the stimulus A, FF in the WTA state was highest of the three states (Fig. 7c). Although average FF of stimulus B was much lower (≈7) than that of the stimulus A (≈25) in the WTA state, it was still higher than that measured in AI and TA states. These observations suggest that evoked responses are far more variable across trials in the WTA state than in both AI and TA states.
Taken together, the stimulus response properties of the three ongoing network activity states in gamma networks showed that only AI and TA states provide a suitable substrate to reliably encode different external stimuli. In contrast, in the WTA state, the response depends not only on the stimulus amplitude but also on the spatial location of the input. When the input coincides with an existing bump, the response is reliable and discernible from the background, but otherwise the response is weak and unreliable (e.g., stimulus B).
Modulation of pairwise correlations is maximal in the TA state
AI and TA states are similar in terms of stimulus sensitivity and response reliability (Fig. 7). Further analysis of the evoked activity, however, revealed a crucial difference between these two network states that renders the TA state as a more suitable ongoing activity state for stimulus encoding. In the AI state, a weak external stimulus affects only the rate of stimulated neurons and, thereby, induces only a small effect on neighboring neurons. Therefore, the correlation spectrum during the evoked activity is slightly positively skewed (Fig. 8, left, light blue trace). That is, although weak stimuli can evoke activity responses in the AI state, the spectrum of correlations remains largely unaffected (Fig. 8).
In contrast, in the TA state, even a weak external stimulus can create new activity bumps and, thereby, introduce both positive and negative correlations (Fig. 8, middle, light blue trace). The resulting distribution of correlations in the evoked activity is much wider than observed in both AI and WTA states. That is, in the TA state even weak stimuli are able to induce large changes in the structure of pairwise correlations (Fig. 8, middle, light blue trace). These correlations may be effective in carving out the selected action not only in the striatum, but also downstream in external and internal GP, both of which have high baseline activity and require coordinated striatal inhibition to be suppressed. Moreover, the change in the pairwise correlation distribution can also be useful in distinguishing stimulus-evoked bumps from spontaneously generated bumps, which otherwise cannot be distinguished based on firing rates alone.
For strong external stimuli, the correlation spectra are similar in both AI and TA states. In the WTA state, however, a weak stimulus did not induce any visible modulation in the correlation spectrum (Fig. 8, right, blue and light blue traces), and only very strong inputs could trigger a small change in the correlation structure (Fig. 8, right, red and orange traces). These results suggest that although both AI and TA states allow for reliable and discernible stimulus responses, the TA state may be more suitable to process weaker stimuli than the AI state, because in the TA state the correlation spectrum can already be modulated by very weak stimuli.
Discussion
Here, we investigated the activity dynamics and stimulus response properties of the striatum as a purely inhibitory network with different spatial connectivity profiles. We showed that a nonmonotonically changing spatial connectivity profile can lead to the emergence of spatially structured activity in purely inhibitory networks. In contrast, when the connectivity changes monotonically as a function of distance between neurons, the network activity is uniformly distributed over the network, respective of the background input. Specifically, we have shown that with a nonmonotonically shaped connectivity profile, the striatal network can exhibit three qualitatively different activity states: AI, TA, and WTA dynamics. Importantly, among these three different dynamical network states, both AI and TA states have the necessary properties for reliably encoding external stimuli. Between the AI and TA states, the average stimulus response is similar; however, the TA state has several interesting properties that make it a dynamically richer and more responsive state: (1) the overall firing rate distribution in the TA state is more skewed than in the other two states (Fig. 2a); (2) the lifespan of NAs in the TA state is more widely distributed that in the other two states (Fig. 4d); and (3) in the TA state, even weak stimuli are able to alter the spectrum of pairwise correlations (Fig. 8).
Transient NAs in striatum-like purely inhibitory networks with spatial connectivity structure (Humphries et al., 2009) or without any spatial connectivity structure (Ponzi and Wickens, 2010) have been defined as groups of neurons showing a conspicuous correlation in their temporal firing rate profiles. Such NAs and their member neurons were identified by offline analysis of the spiking activities of neurons in sparsely connected, random recurrent inhibitory network models with weak synapses. In such networks, NAs were found to be randomly distributed over the entire network and appeared to involve mutually unconnected neurons (Angulo-Garcia et al., 2015). However, such assemblies may not influence the downstream network, unless they are specifically wired to share their downstream targets. Moreover, experimental data showed that transient NAs exist as spatially compact clusters (Barbera et al., 2016). Here, we extend the previous work and show that the existence of spatially compact NAs in a striatum-like network requires that the connection probability between striatal neurons changes in a nonmonotonic fashion as a function of their distance.
Previously, a few studies have investigated the dynamics of the striatal network with distance-dependent connectivity (Wickens et al., 1995; Humphries et al., 2010), albeit without focusing on the role of the connectivity structure and the external input. Wickens et al. (1995) showed that when the connectivity of striatal neurons is confined to a small neighborhood, the network can exhibit multiple spatially localized, persistent bumps (akin to the WTA state). Such multiple bump–state activity required that the connections were symmetric. Any heterogeneity and asymmetry in connections led to traveling waves or AI-type activity. However, that network was very small, and the connection probability was fixed over a finite distance. Humphries et al. (2010) developed a more realistic model of the striatum using distance-dependent connectivity estimated from the 3D morphology of MSNs. Structurally, that network is similar to the gamma network we investigated here. However, Humphries et al. (2010) investigated neither the dynamics of the network as a function of the input nor the relationship between network structure and dynamics. Our work complements and extends these previous studies and, importantly, explains how external input and the structure of spatial connectivity profiles (Gaussian and gamma) interact to shape the dynamics of spatially compact activity clusters in the striatum.
Relevance for striatal network activity dynamics
In the ongoing network activity state, striatal neurons are relatively silent, and task-related activity can increase up to 20 Hz (Gage et al., 2010). Recently, advances in recording methods have enabled recording 10–100 neurons simultaneously, using calcium imaging in behaving animals. Analyses of such high-density sampling of striatal neurons showed that striatal activity is organized as spatially compact clusters of coactivated neurons (or NAs; Barbera et al., 2016; R. Costa, personal communication). Similar observations were made earlier in in vitro (Carrillo-Reid et al., 2011) and in task-related activity in behaving monkeys (Adler et al., 2013). In an unhealthy low-dopamine state (as in Parkinson’s disease), D2 type dopamine receptor expressing striatal projection neurons increase their firing rates (Mallet et al., 2006), and in general, striatal neuronal activity loses its diversity and ability to switch neuronal activity in a task-dependent fashion (Costa et al., 2006; Costa, 2011).
It is conceivable that the spatially compact clusters (or NAs) are a simple consequence of spatially localized cortico-thalamic inputs to the striatum. However, this is a trivial solution and indicates that the striatum acts as a simple transmitter of cortico-thalamic activity to downstream targets. Instead here we argue that intrinsic dynamics of striatum-like inhibitory networks are able to generate NAs even when cortico-thalamic inputs are not spatially compact.
The three network states that we identified in our study capture different aspects of the ongoing and evoked activity of the striatum in normal and in low-dopamine states. Both AI and TA states in our network models match some properties of the ongoing activity measured in vivo. In data from experiments in which animals wait for a cue to initiate a task (e.g., Gage et al., 2010; Adler et al., 2012), the ongoing activity in the striatum appears be similar to the AI state. When the animal engages in the task, after the cue presentation, the AI activity is transformed into a TA-like state (Gage et al., 2010; Adler et al., 2012). On the other hand, in data from freely moving animals, although the animals are not performing any goal-directed behavior, striatal activity shows spatially compact clusters of coactivated neurons, similar to the TA state (Barbera et al., 2016). Whether the striatal activity recorded in freely moving animals can be treated as an ongoing activity, a goal-directed activity, or a combination of both is an open question. We hypothesize that because animals were not engaged in specific goal-directed behavior, the data reported in Barbera et al. (2016) represent an ongoing activity state in the sense that any task-related cue or reward was absent. It will be interesting to explore whether these spatial activity clusters also emerge in a goal/cue-directed task, thereby possibly hinting at differences in the contextual input received by the striatum during ongoing and goal-directed states.
Thus, in our opinion, depending on the context (e.g., freely moving or cued goal-directed behavior), the ongoing activity of the striatum can be in an AI or TA state. As stated earlier, both states allow for a reliable and discernible response in the striatum. A subtle, but important, difference between the AI and TA states is that in TA-type ongoing activity, even weak stimuli can induce a large change in stimulus-induced correlations. Assuming that neuronal correlations form the basis of modification of synaptic strengths, we speculate that when an animal operates in a TA state, even weak stimuli can drive learning.
In the WTA state, neurons have higher firing rates and the spatial bumps are stable (Figs. 3 and 4). In this state, the ongoing activity shows only a very small diversity, and only very strong inputs, arriving in specific locations, can induce any discernible and reliable response. With these properties, the WTA state resembles the striatum dynamics in Parkinson’s disease (Costa et al., 2006; Costa, 2011). Moreover, the WTA state is observed when neurons are spiking at a high rate (Fig. 4), which could be achieved either by increasing the ongoing external excitation or by increasing the excitability of the neurons. This is also consistent with the fact that in Parkinson’s disease the cortical drive to the striatum is increased due to the potentiation of cortico-striatal synapses (Smith et al., 2009b; Fieblinger et al., 2014), and the lateral inhibition among MSNs is decreased or even disrupted (Taverna et al., 2008).
Here, we estimated the stimulus response for one particular size of the external stimulation (i.e., the number of stimulated neurons and their spatial distribution) and varied only the magnitude of the stimulation current. The size of the external stimulus may also influence the stimulus-response magnitude and reliability. We expect that increasing the number of neurons in a fixed region will increase the response reliability and magnitude. However, increasing the stimulated region with a fixed number of neurons may have a nonmonotonic effect on the response magnitude and reliability. As long as the stimulated region is smaller than the size of an individual bump, the response magnitude and reliability will increase. Increasing the size of the stimulated region beyond the single bump will recruit surrounding neurons, which should be inhibited by the bump itself, and hence, the stimulus response may decrease. In another scenario, distributing the stimulated neurons in small islands may have nonlinear effects on the response dynamics. A detailed analysis of the relationship between stimulus size and response dynamics is a complex topic and should be addressed in a separate study.
Despite the similarities and interesting insights we noted regarding the striatal activity in the healthy and in Parkinson’s disease states, we note that our model is highly simplified and ignores several key features of the striatal network. Specifically, it would be necessary to study how other components of the striatal network [e.g., interneurons, separation of the striatal network into D1- and D2-type medium spiny neuron populations (Taverna et al., 2008; Bahuguna et al., 2015), neuromodulators] affect the stability and stimulus responses in the AI and TA states. In addition, we have used only simple LIF model neurons in our networks. Voltage-dependent ion channels may introduce nonlinear effects and change the network dynamics qualitatively (Wickens et al., 1995; Rinzel, 1998). Hence, in future work, it will be important to understand how the three dynamical states identified in our study are possibly altered when neurons are endowed with voltage-dependent ionic currents or larger heterogeneity in neuron/synapse properties.
Model validation
Recent experimental data show that the ongoing activity of the striatum can be either in an AI state or a TA state, depending on the experimental conditions (Gage et al., 2010; Adler et al., 2012; Barbera et al., 2016). Furthermore, these data also suggest that the striatum exhibits a TA-like state as animals get engaged in a task. In addition to direct measurement of NAs using modern imaging tools (Barbera et al., 2016), our model suggests that even the relationship between the firing rate and spike time irregularity (Fig. 2) can provide further indirect evidence for the existence of NAs in the striatum.
The key to the emergence of spatially compact transient NAs in the TA state is the nonmonotonically shaped connection probability. Indirect estimates of the anatomic (from neuron morphologies) and functional connectivity both indicate that MSNs do not inhibit their nearest neighbors and that the connection probability peaks at a distance of ≈80 µm and then decays to zero beyond ≈200 µm (Fujiyama et al., 2011; López-Huerta et al., 2013). Similar estimates regarding the distance-dependent connectivity between MSNs have been drawn from computational analysis of the 3D morphologies of MSN axons and dendrites (Humphries et al., 2010). However, more experimental work is required to measure the spatial profile of not just the functional, but also the structural, connectivity within the striatum, in particular to measure the spatial connection profile and spatially compact neural clusters that would support our network model. In addition, our model predicts that neurons participating in the NAs should have a low connection probability and share their inputs. In addition, our model predicts that the CVISI values of MSNs should increase as their firing rate increases (Fig. 2b). Finally, we predict that the pairwise correlation spectrum should be more susceptible to weak inputs in the TA state than in the AI state (Fig. 6). High-density sampling of striatal neurons, when available, would be sufficient to check these predictions.
Analytical description of connectivity profiles. The graphs show the spatial connectivity profile (left) and its Fourier transform (right) as a function of the distance between neurons and wave numbers, respectively, for both gamma (red) and Gaussian (blue) connectivity kernels. Note that the spatial connectivity profile remains positive for both connectivity kernels. However, their Fourier transforms behave differently: the Gaussian kernel remains positive, whereas the gamma kernel takes negative values for larger (absolute) wave numbers.
Impact of the network dynamics on the stimulus response. a, Spatial distribution of the spiking activity displayed in time series for different dynamic states (rows). Each frame shows a spatial map of 30 × 30 neurons from the ROI (black squares) in a time window of 100 ms, with 500-ms intervals between successive frames. The gradient background is the probability area for stimulated neurons, and its color refers to a stimulus phase. b, The change of response activity of the stimulated neurons to their corresponding stimuli A (red) or B (orange) in different ongoing bump states (AI, TA, and WTA). Each row represents the strength of the stimuli (50–150 pA). A lower Δresponse indicates a weaker impact of the external stimuli on the network activity, and lower variance of activity reflects a higher reliability of the response. For each subpanel, the white lines are the median value of the data. The colored boxes extend from the 25% to 75% of the data, i.e., the box contains ≈50% of the data. Whiskers extend from minimum to maximum values of the data. c, The temporal variability (FF) of the response of the stimulated neurons as a function of time. A lower FF indicates a higher reliability of the stimulus response. A higher FF is observed at each stimulus onset, in both the AI and TA states. In contrast, in the WTA state, the network is not able to reliably respond to external stimuli. Both stimulus phases are displayed at the bottom of each subpanel.
Impact of the network dynamics on the modulation of the spectrum of pairwise correlations. Different subfigures show the pairwise correlation spectrum of the network activity in three dynamic network states (AI, TA, and WTA). Different colored traces represent different stimulus strengths on the selected neurons. Compared with the correlation spectrum in ongoing activity (blue), a higher excitation is required to modulate correlations in the AI state than in the TA state. With a stronger external excitation, the correlations are more widely distributed in the network activity in both AI and TA states.
Relationship with models of cortical networks
Recurrent networks with both excitatory and inhibitory (E-I) neurons, interconnected according to a monotonically decaying (e.g., Gaussian) connectivity profile for excitation and inhibition, can be tuned to exhibit spatially clustered or stationary bump type activity (Ben-Yishai et al., 1995; Roxin et al., 2005). A key feature of E-I networks that show spatial clusters is that the excitatory connectivity decays more rapidly with distance than the inhibitory connectivity. In such networks, the summation of excitatory and inhibitory connectivity kernels or excitatory and inhibitory synaptic strengths yield the well-known Mexican hat profile as the effective connectivity kernel with its characteristic, nonmonotonic shape. With such a connectivity profile, local recurrent excitation activates neighboring neurons which, in turn, inhibit the surrounding region because of the stronger distal inhibition. That is, in such E-I networks, a coactivated local group of neurons is brought together by their mutual, predominantly excitatory, connections and by their common field of surround inhibition. Examples of such behavior have been reported in the experimental literature, e.g., in the monkey prefrontal cortex (Vaadia et al., 1995). In contrast, we found that for purely inhibitory networks, using both network simulations and neural field equations, a nonmonotonic spatial connectivity kernel (such as the gamma distribution) generates spatially clustered bump type activity for high background input. As shown by our mean field analysis, purely inhibitory networks with a Gaussian spatial connectivity profile cannot possibly support any spatially periodic bump activity. Intuitively, that is because coactivation of neighboring neurons requires that these neurons do not inhibit each other while creating an inhibitory surround. Thus, in purely inhibitory networks, a coactivated local group of neurons is defined by the lack of mutual connectivity and the presence of a common range of surround inhibition. Mathematically, as we have shown, the key condition to have spatial clusters of activity is that the effective connectivity kernel has a nonmonotonic shape as a function of distance in the network (Fig. 6 and “Necessary conditions for the emergence of bump states”). Both the Mexican hat–shaped effective connectivity kernel of E-I networks and the gamma distribution shaped connectivity kernel of purely inhibitory networks fulfill that condition. We note that the stable bump state (WTA) observed in our network models closely resembles the grid patterns observed in the medial entorhinal cortex of rodents and that computational models of grid cells also use a nonmonotonic kernel to form connections between inhibitory neurons (Couey et al., 2013; Roudi and Moser, 2014).
In summary, we have shown how the shape of the distance-dependent recurrent connectivity profile and the strength of ongoing external background excitation together determine the state of the ongoing network activity as well as the stimulus-response properties in a purely inhibitory network, such as the striatum. These results, when properly adapted to the specific inhibitory network of interest, could provide important new insights into the functional characterization of the activity dynamics in inhibitory brain networks such as the striatum, globus pallidus, and central amygdala.
Acknowledgments
Acknowledgments: We acknowledge helpful discussions with Drs. Robert Schmidt, Ioannis Vlachos, and Alejandro Fernandez-Bujan.
Footnotes
Authors report no conflict of interest.
Partial funding of this work by the German-Isreali Foundation for Scientific Research and Development (GIF), the German Federal Ministry of Education and Research (BMBF) grant 01GQ0830 to the Bernstein Focus Neurotechnology (BFNT) Freiburg/Tuebingen, the NeuroTime EU-Erasmus Mundus Joint Doctorate, Parkinsonfonden Sweden, and the Carl Zeiss Foundation is gratefully acknowledged.
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
References
In this issue
