The conductance-based neural model of a striatal medium spiny neuron is described in detail in 6. The membrane voltage V _m is modeled using the equation

V ˙ m = − 1 C [ μ 1 ( I 1 ) + μ 2 ( I 2 ) + ... + μ n ( I n ) − I s y n ] , ( 1 )

where the I _i are the currents, induced by the individual ion channels. Variability of the neuron is modeled by modifications to µ _i . This model includes ion channels for Na (INa), K (IK), slow A-type K channels (IAs), fast A-type K channels (IAf), inward rectifying K channels (IKir), L-type calcium channels (ICaL), and the leak current (I _leak). The definition of all parameters and the dynamics of the ion channels can be found in 6.

For the experiments in this paper, we only focus on variability induced by changes in the strength of the slow A-type K channels. The total current contribution for this channel is µ _IAs where µ was selected between 1.0 and 1.5.

In order to illustrate the variability in neuron behavior, we excited the neuron model by input signals, resembling two kinds of synaptic input: uncorrelated and correlated. These signals were generated by superposition of excitatory and inhibitory spikes from individual Poisson-distributed spike trains (50 excitatory and 10 inhibitory), and biased Gaussian background noise.

The amount of pairwise correlation in these spike trains governs the type of input signal. A high correlation factor was used in order to generate sequences which have short periods (10–15ms) of high activity.

In order to do large-scale simulation we needed to employ a simple, computationally tractable neuron model. We used a two-dimensional model of a neural oscillator (cf. 7), and employed an instantiation of the model with parameters fitted to the general properties of cortical pyramidal neurons ¹³ as a generic model (g). The model consists of an equation for the membrane model v ( Equation 2), fitted to experimental values for cortical pyramidal neurons, and an equation for a gating parameter u ( Equation 3).

                                                                                         v ˙ = 0.04 v ² + 5 v + 140 − u − I _syn                                         (2)

                                                                                                            u ˙ = a( bv − u)

                                                                                                            b = 0.2                                                              (3)

                                                                                                            a = 0.02

When the neuron fires a spike (defined as v( t) = 30 mV), v is set back to a low membrane potential v := c; c = −65.8 mV and the gating variable u is increased by a fixed amount d ( u := u + d; d = 8) (cf. 13). This formulation allows for a very simple neuron model, which avoids the explicit modeling of the downslope of the action potential, and rather resets the voltage. Time-dependence is modeled by the gating variable u.

Neuronal heterogeneity is achieved by systematic variation of inactivation parameters. By varying d, we can vary the inactivation dynamics of the model after a spike, by varying a we vary the inactivation dynamics throughout the computation. In this way, we can attempt to model neuronal variability in activation/inactivation dynamics, which is sufficient to model frequency-selectivity as a stored intrinsic property. The parameters used in this paper for different neuron types are listed in Table 1.

We created graphs of N (= 210) excitatory neurons, and K ( ∼ 1900) excitatory connections. For the excitatory neurons, we use randomly connected graphs (N,K) with different dispersion σ ^*. This corresponds to normal (Gaussian) to lognormal graphs with different widths and length of the heavy tail. We model inhibition by Poisson-distributed inhibitory synaptic input directly onto excitatory neurons.

We use specific instantiations of these graphs (RG, LG1) for the simulations. Table 2 shows global graph characteristics for the Gaussian graph (RG), the lognormal graphs (LG1), and intermediate graphs LG2, LG3, and LG4.

We define synchronization s in a network by pairwise correlations: for each neuron n _i , we count, for each other neuron n _j , the number of spikes which occur within a window W ( W = 10 ms) of n _i ’s spike events, divided by the total number of spikes for n _i . More precisely, for each neuron, we bin all firing events into 5ms bins. We then count the number of other neurons, which fire in a 10ms window around the (start) of the bin. The synchronization s is then the average over all neuron pairs in the network:

s = 1 N ( N − 1 ) ∑ i 1 S j ∑ i B i j ,

where B _ij is the number of spikes that neurons i and j have in common within a moving window of W = 10 ms during the entire measuring time. S _j is the number of spikes of neuron j during the entire measuring time.

The rewiring algorithm used to change the properties of a graph G is a greedy algorithm, which iteratively selects the node s with the highest degree. One of its edges is then rewired to random nodes with lower degrees, decreasing σ ^*. The algorithm terminates, when the value of σ ^* falls below a given threshold.

All simulations were performed with the software tool CNeuroSim, which is implemented in Matlab (R2016b) and C, and available at 14.

We show how we can model gain as a stored intrinsic property, defined as a spike rate in response to constant or fluctuating input of fixed strength (A/Hz). We use a full ion channel based model (the MSN model ⁶ ), with variation in the slow A-type potassium channel (IAs).

In Figure 1A, we show the response of MSN model neurons with a scaling of µ _IAs = 1.0, 1.3, 1.5 to a noisy signal, derived from uncorrelated Poisson-distributed synaptic input. The top panel shows the development of the membrane potential, V _m , over time for all neurons. The middle panel shows the spike-train for each neuron with the mean ISI and its standard deviation; the total number of spikes is shown on the right. The total number of spikes includes bursts, which were excluded from ISI calculation. The bottom panel shows the synaptic input. The dots correspond to the spiking events for a single neuron 3 ( µ _IAs = 1.5). The resulting mean ISIs are 25, 37, and 45 ms. With a standard deviation of 6, 11, and 8, they are clearly distinguishable. This is also shown by the Gaussian distribution for the mean ISIs for each neuron type ( Figure 1B).

This model shows frequency-specificity as read-out of the relative contribution of the slow A-type potassium channel, indicated by the scaling factor µ _IAs. The relative contribution of an ion channel corresponds to its density or distribution on the somato-dendritic membrane, or in some cases its specific localization at dendritic branch points. Experimental evidence has shown that this is a plastic feature for neurons.

We then employ highly correlated synaptic input, defined as in 6 (see Methods). We stimulate the same neurons with the correlated input and observe the spike pattern ( Figure 2). We can show that the frequency-specificity of the neuron disappears. Instead we see a time-locked spike pattern which is expressed by a similar spike frequency ( Figure 2A) and an overlap of the mean ISIs ( Figure 2B). What this experiment showed is that a stored intrinsic property, the gain, is available to the processing network in a conditional manner. The property is continually expressed, the differences in ion channel density persist. Depending on the mode of stimulation, however, this property is manifested as intrinsic gain, or it is obscured when a neuron is driven by strongly correlated input.

To continue with exploring this property of model neurons in a networked context, we switched to a simplified model neuron ⁷ and created a set of variations for this model (see Methods). We show the response of two-dimensional model neurons to asynchronous input in Figure 3, and to regular, synchronous input in Figure 4. In the first case, we have clearly separated frequencies, and in the second case, the ISIs are nearly identical with a narrow distribution. When we stimulate the neurons with irregular, but synchronous input, the ISIs become identical, but with a wider distribution to reflect the different duration of pauses between the synchronous stimulation ( Figure 5).

We may also consider the question of whether a neuron can simultaneously respond to an input and read out its stored spike frequency. If there are single synchronous events, which interrupt ongoing spiking, can we recover the intrinsic properties for each neuron? In Figure 6 it is shown that this is possible. Figure 6A shows the input and the synchronous responses, and in Figure 6B we still see a clear separation of frequencies.

We conclude that we can multiplex asynchronous and synchronous input. It is also apparent that there needs to be a lower limit on the intervals between synchronous events that can be processed without disrupting intrinsic properties. This interval needs to be defined as functionally dependent on the intrinsic frequencies. In this case, it is 3/s for the synchronous events, with 10Hz for the slowest neuron.

The simplified model neurons allow the creation of large networks of heterogeneous neurons and exploration of different topologies. We hypothesized that a lognormal graph, because of its hierarchical topology and the existence of hub neurons would lead to synchronization of action potentials – even with heterogeneous neurons – while a Gaussian topology would support asynchronous spiking behavior ^{15, 16} . We define synchronization s in a network by pairwise correlation ( Methods). The spike frequency for each neuron type is assessed by the mean and standard deviation for ISIs, as before.

We first use a Gaussian connected graph (RG) with N = 210 and K = 1800 and use 7 different neuronal types (1–6, plus the generic neuron g) with 30 Neurons each ( Methods). Figure 7A shows an excerpt of the graph structure. We can see that the graph is connected such that all neurons have a comparable number of connections. This is also apparent in Figure 8, where we can see a (narrow) normal distribution for connectivity for the Gaussian graph RG. Table 2 contains the graph characteristics for both graphs.

N = 210 is about the size of a minicolumn or ensemble unit within a larger network with presumably dense interconnections ¹⁷ . The maximal density d = K/( N×N−1) in a cortical microcircuit is estimated at 0.1 for 10 ⁴ neurons, and 10 ⁷ synaptic connections, ¹⁸ . With ˜50% of synapses internal to the network, d = 0.04 − 0.07 is a realistic value for internal connectivity ¹⁷ . There is also a small background inhibition to all neurons present, implemented by 10% inhibitory neurons with Poisson-distributed firing and complete connectivity to excitatory neurons.

We now stimulate the graph by an initial stimulation to 10 excitatory neurons (for about 1s). In Figure 9A, we see highly asynchronous neuronal activity after 1s of stimulation. The pairwise correlation values s is low ( s = 0.11). Figure 9B shows that each neuronal type retains its own frequency, i.e., has its own typical ISI, separated from other neuronal types. We also notice that some neurons fire with low frequencies (5Hz) and others with higher frequencies (20Hz). Very low firing neurons (2Hz) which are typical for cortex are not represented in this model.

Next we changed the topology of the network to a lognormal graph (LG1), with N = 210 and K = 1924 and used the same neurons as before.

Figure 7B shows an excerpt of a lognormal graph structure. The connectivity structure seems much denser, because of ’hub’ neurons in the center of the graph. In Figure 8, we can see a much wider distribution of degrees for the lognormal graph (blue), with a number of nodes with high connectivity. Presumably, those nodes are capable of synchronizing the network, because they can reach many neurons simultaneously. What is the effect on the presence of neural heterogeneity?

Figure 10 shows that a high amount of synchronization can be achieved in spite of heterogeneity of intrinsic frequency of model neurons. The rasterplot ( Figure 10A) shows the activity in LG1 with the same neurons and the same stimulation as before. The overall correlation, defined by pairwise correlation of neurons, is much higher (s=0.32). The distribution of ISIs in this case is strongly overlapping ( Figure 10B), similar to Figure 4B, where neurons were explicitly driven by highly synchronous input. This means that synchronization is dependent on the network topology, and a lognormal graph exhibits a higher tendency for pairwise synchronization. Also, that neuronal heterogeneity is apparent in an asynchronous network mode but is repressed in a synchronous firing mode.

We could show that differences in intrinsic properties appear or become more prominent when there is less synchronicity in a network. In our model, the pairwise synchronicity s is dominated by the network topology, more precisely by the width of the degree distribution (dispersion) ranging from Gaussian to lognormal with a heavy tail.

To confirm this observation we used a number of intermediate graphs and mapped the pairwise synchronization dependent on the dispersion σ ^* ( Figure 11). The graphs RG and LG1 that we used have values of σ ^* = 1.44 and σ ^* = 2.89 ( Methods). They have the same density, i.e., the same number of connections and neurons ( 0 ˜ .05). Additionally, we analysed the dependence of synchronicity on the density of the graph between 0.01 and 0.1 ( Figure 11).

There is higher synchronization in the lognormal region, especially with higher σ ^* > 2.5 but no synchronization for Gaussian graphs. For heavy-tail graphs, synchronization depends linearly on the density between d = 0.03 − 0.08 ( s = 0.2 − 0.5).

How are the different graphs related? We hypothesized that fast synaptic switching ¹⁹ by neuromodulation could change the network topology sufficiently to switch from a synchronous to an asynchronous regime. In Figure 12, we plot the number of edges that were changed to achieve different dispersions σ ^* of a graph. The algorithm used was a simple greedy algorithm ( Methods), which is suboptimal, i.e. overestimates the number of edges required. It appears that 30–50% of edges changed would be sufficient.

We employ a parameterizable two-dimensional neural oscillator model to encode different intrinsic excitability manifested by different frequency responses to constant input. What the experiments show is that a stored intrinsic property, the gain, is available to the processing network in a conditional manner: The gain is continually present, the differences in ion channel density persist. Depending on the mode of stimulation, however, this property is manifested as spike rate, or it is obscured when a neuron is driven by strongly correlated input. This is interesting because it shows a property of memory that synaptic plasticity lacks: the memory is not always ’read-out’ in any processing step. It is conditional, it can be accessed or ignored depending on the state of the network. This seems to be an essential property of memory in any intelligent system.

Different statistical properties of synaptic input can be modeled by a variability in the correlation properties of input neurons. In a network model, this means that the overall correlation in the network determines what input a neuron receives. With a Gaussian degree distribution topology, correlation is low and neurons fire irregularly with their own preferred frequency. With a heavy-tailed, lognormal topology, correlation is higher, and neurons fire when they receive correlated input, irrespective of intrinsic properties. I.e. driving neurons by correlated vs. uncorrelated input leads to uniform spiking behavior vs. read-out of stored differences in ion channel conductances.

A restriction of the present model with respect to a biological simulation model is the simplified treatment of inhibition. However, experimental work shows that cortical parvalbumin-expressing (PV+), fast-spiking interneurons have no connection specificity to pyramidal neurons, rather they present as an ’unspecific, densely homogeneous matrix covering all nearby pyramidal cells’ ( 20, p. 13260), which corresponds to our model.

Conditions for neuronal read-out may include the activity of inhibitory neurons. Inhibition and excitation are tightly linked by feedback interaction. 21 suggested that the close coupling of inhibition and excitation in cortical tissue cancels out purely input-dependent, i.e. not network generated synchrony. 13 suggested that with highly correlated input, both inhibitory and excitatory, the neuron may receive less input which allows it to be driven only by strong synaptic input, with distributed input it may receive a barrage of excitatory and inhibitory inputs where the membrane voltage remains close to firing threshold and the neuron fires continuously. In our sense, it is reading out stored properties. Both inhibitory and excitatory synaptic input will conform to be asynchronous or synchronous, to drive neurons or it causes a neuron to emit spikes according to its own stored intrinsic properties.

However, neuromodulation has effects on inhibitory neurons as well ^{22, 23} , which we have not modeled. Further simulations will show whether the I-E coupling is altered during enhanced neuromodulation, or whether the effects are synergistic to the present results.

Neuromodulation influences both intrinsic properties and synaptic connectivity ⁵ , e.g. acetylcholine, (via nucleus basalis stimulation), noradrenaline (via LC stimulation) or dopamine (via VTA stimulation) ^{5, 24} . Experimental estimates on the distribution of synaptic neuromodulatory receptors are at approximately 30%–50% of connections ¹⁹ . That is sufficient to transform the topological properties of a graph, such as the dispersion of its degree distribution from heavy-tailed graph to a more Gaussian, less clustered graph without requiring tight optimization for the positions of neuromodulatory receptors ( Figure 12).

Neuromodulation disables or enhances various ion channels, such as Sk-channels which guide reset times after a spike, or A-type potassium channels which influence latency to spike ⁶ . In this way, neuromodulation influences intrinsic properties. If neuromodulation reduces synchrony by acting at synaptic receptors, it uncovers intrinsic heterogeneity, and induces a mode of processing that allows read-out and storing of intrinsic properties. Depending on the neuromodulator used, and the amplitude and duration of the signal, different soma-dendritic ion channel profiles would emerge.

In the synchronous mode, intrinsic heterogeneities are reduced in the presence of tightly correlated input which drives neurons reliably. This invariance of neuronal intrinsic properties in synchronous mode allows synaptic transmission and information processing independent of neuronal heterogeneity.

The idea of introducing synchronous events by common input to an asynchronous background, and in this way use reliable synaptic transmission without affecting the state of the system (multiplexing) has also been documented in experimental results. For instance, ( 25, Figure 4A) shows a case of multiplexing in response to behavioral stimuli. In this case, intrinsic read-out can continue, and single events are transmitted reliably through driven activations.

Why should synchronization properties be switched by neuromodulation? Increased correlation in the network supports population-coded information to be propagated effectively ²⁶ . Turning on neuromodulation would decorrelate an area and increase the capacity for information coding in an ensemble or a cortical microcolumn ²⁷ . This area would become an information source to surrounding areas. When turned off, increased correlation would allow this area to transmit information and to disregard the stored neural memory.

Basal forebrain stimulation, which results in increased acetylcholine release and muscarinic/nicotinic receptor activation, decreases correlation between cortical neurons ( 28– 30 ( Figure 3.C)). Likewise, 31 ( Figure 3 and Figure 4) shows reduced noise (internal) correlations with cholinergic stimulation, while inactivation of the basal forebrain caused more synchronized activity. 32 shows reduction of correlation for task-relevant perception, where presumably task-relevance causes neuromodulatory activity. 33 provides evidence for the involvement of noradrenaline in desynchronization of cortical state and the enhancement of sensory coding.

There is considerable evidence ^{34– 37} showing that several neuromodulators, including at least noradrenaline and acetylcholine, modulate pairwise spike correlation, such that strongly synchronized states (anesthesia, slow wave sleep) have high correlation and low neuromodulation, while asynchronous states (normal waking), with higher neuromodulation, have lower pairwise correlation.

37 observed intrinsic fluctuations in synchronization of cortical networks during wakefulness which correlated with the amount of encoded perceptual information and perceptual performance. Their results showed a mean decrease in correlations from synchronized to desynchronized state corresponding to perceptual performance by approximately 20%, similar to values observed during attention ³⁸ , and after adaptation ³⁹ . We have shown ( Figure 11) that correlation changes are continuous with network topology and a 20% correlation change is well within the range of the current simulations. Importantly, the results in 37 point to fluctuations in synchronization that reflect local changes in network activity rather than just global cortical state dynamics which have traditionally been associated with central neuromodulatory release.

The role of presynaptic neuromodulation in suppressing cortical connections ^{11, 12} and changing attractor states ⁴⁰ , as well as allowing rapid synaptic weight changes ¹⁹ has previously been assessed. Theoretical work has also emphasized the connection between correlations and information content ^{27, 41– 43} .

Here we bring these observations together to suggest that neuromodulation of synapses may alter network topology and in this way bring about an increased decorrelation of spiking, and a more asynchronous state, with a higher informational capacity. It may provide a general explanation (a) on how fluctuations in synchrony can be engineered rapidly and in small cortical areas and (b) why intrinsic memory may be conditional, accessible only at certain times and in a localized fashion.