We analyze five data sets for spike rates from principal neurons under behavioral activation:

In monkey IT, single unit activity was recorded over 200ms for passive viewing of 77 different natural stimuli for 100 neurons, each stimulus shown 10 times ¹⁶ . This yielded 770 spike rate response data points per neuron. For these data, we show mean spike rate, standard deviation, max-min values, coefficient of variation (CV) and Fano factor (FF) ( Figure 1), (cf. 1). What is remarkable is that the dispersion for each neuron (variance related to mean, FF) is fairly constant, and not related to the rank of a neuron as high or low-frequency. In other words, neurons have roughly similar behavioral responses relative to their average spike rate. For this reason,many behavioral experiments have reported percentage of increase/decrease of spiking as the relevant parameter.

Additionally, we show data from primary auditory cortex (A1) from awake monkeys, which were recorded for spike responses to a 50ms, 100ms, or 200ms pure tone ( 17, Figure 2A and B) and data for spike rates from the primary auditory cortex of rats for four different conditions, which were recorded as cell-attached in vivo recordings ( 2, Figure 2C).

For midbrain nuclei neurons (IC), we re-analyzed spike rates in response to tones (for 200ms after stimulus onset) under variations of binaural correlation ²² . The frequency ranking of neurons by mean spike rate, standard deviation, min-max values, CV and FF are shown in ( Figure 3A–C). CV and FF are similar to the cortical data. Data from GABAergic cerebellar Purkinje cells offer some difficulty for this analysis since they have regular single spikes at high frequencies, and in addition, calcium-based complex spikes ²⁰ . Complex spike rates however are low (< 1 Hz). This can therefore be regarded as a form of multiplexing, with two separate codes, where single spike rates can be separately assessed in their distribution. Here we report data for single spikes from in vivo recordings in anesthetized rats ²⁰ ( Figure 4A)and data from spontaneous spiking (in the absence of synaptic stimulation) under in vitro conditions ( 18, 19, Figure 4B and C).

In order to show values for standard deviation and variance, data for two Purkinje cells from a behavioral experiment, ²¹ , i.e. single spike rates during arm movements from monkeys, have been added to the ranking of spontaneously firing neurons by mean spike rate from 18 ( Figure 4D).

The logarithmic (heavy-tailed) distribution of spike rates is evident under all conditions.

The distribution of spike rates for neurons spiking in the absence of synaptic input shows that there are differences in the intrinsic excitability of neurons. To further explore this we looked at three additional datasets, which report the action potential firing of a cell due to injected current, such as by a constant pulse. This can be defined as the neuronal gain parameter (spike rate divided by current, [Hz/nA]).

In 1, we already presented the data from striatum, which show that the spike response to a constant current follows a heavy-tailed distribution ³ . Figure 5A shows the spike rate in response to current pulses of different magnitude in two different areas, nucleus accumbens (NAcb) shell and dorsal striatum. Figure 5B and C show the distribution of rheobase (current-to-threshold) for dorsal striatal and NAcb shell neurons. Distributions appear mostly lognormal, with the exception of the 200pA current pulse response and the data in Figure 5C, which appear normally distributed.

We extend this dataset by recordings from different types of cortical neurons in cat area 17 in vivo ( 23, Figure 6A and B) and from GP in awake rats ( 5, Figure 6C).

A lognormal distribution of intrinsic gain is clearly apparent, except for fast-spiking interneurons, which, however, may be a sampling error (n=33).

Synaptic weight distributions have been investigated starting with ¹⁰ in hippocampus by measuring EPSC magnitude ^{6, 24– 27} ( Figure 7). There is also a review paper available ²⁸ to summarize the findings. Recently, the expression of AMPA receptor subunit GluA1, which is correlated with spine size, has also been measured ²⁹ , Figure 8. We used five datasets from cortex, hippocampus and cerebellum:

In 6, EPSP magnitude was measured for L5 pyramidal neurons in slices from rat visual cortex, averaged over 45–60 responses, and peak amplitude recorded, ( Figure 7A). Similar data were used in 26 for slices from a single barrel column in rats ( Figure 7B). A 30-fold variation of coupling strength was noted. In 10, EPSPs between CA3 and CA1 in hippocampal slices were recorded, by detecting somatic membrane potential changes in response to presynaptic neuron stimulation ( Figure 7C). There are also synaptic weight data on granule cell to Purkinje cell connections ^{11, 27, 28} , which show a similar distribution, but have an order of magnitude weaker connections than cortical connections ( Figure 7D). Finally, a different type of evidence was obtained in 29, namely labeling for a subunit of AMPA receptors in layer 2/3 mouse barrel cortex in vivo both before and after whisker stimulation. The AMPA intensity is distributed lognormally over the spines, corresponding to the observations on the strengths of EPSPs. It is noticeable that stimulation leads to increase of on average 200% (two-fold) in about 30% of spines ²⁹ . Yet as we know, over time the overall distribution of synaptic strengths remains stable.

For synaptic weights, just as for intrinsic gains and spike rates, lognormal distributions have been found for both EPSPs and AMPA receptor distribution in a highly consistent manner.

In many cases, the data were only available in the form of histograms. The parameters of the lognormal distribution were then obtained by fitting the data using a Nelder-Mead optimization method. A number of parameters were derived from these fits and reproduced in Table 1– Table 3, cf. 3.2.

Given are two neuron populations I, J with each n = 1000 neurons and variable random connectivity C between I and J. C determines the density between I and J. The input population I always has excitatory output onto J. Inhibitory input to J is modeled by a population H with n = 200. The output neuron population J may also have recurrent excitatory connectivity. Figure 9 shows the architecture for the generic neural network used. The model (GNN) was programmed in Matlab, and is available in GitHub https://github.com/gscheler/GNN . The input population I is modeled according to ² for pyramidal cortical neurons with a spike rate distribution of µ* = 4.95 and σ* = 1.98 ( Table 1). The goal is to generate a spike rate distribution R ^J for J, given a gain distribution G for the target neurons, and the weight distribution W _IJ such that R ^J is similar to R ^I .

For each neuron j in J the spike rate r _j is calculated by applying its gain g _j to the weighted sum of its connected excitatory inputs and its inhibition. C _j is the set of neurons from I that have excitatory connections to neuron j.

r j = g j ( ∑ i ∈ C j w i j r i − r j H ) ( 1 )

where the rate r _i is taken from the distribution R ^I , r j H from R ^H , w _ij from W _IJ , and g _j from G. g _j is modeled as a factor for a linear gain function. It is possible to use a sigmoidal gain function instead, but this makes no difference for the conclusions from the model (Section 3.3). The output R ^J may used as input to I with a matrix W _I,J for tests of the adaptation rules.

For the adaptation of weights W _IJ and gains G, we use Hebbian or homeostatic rules, as described in Section 3.3. The system described in this way is sufficient for all the calculations on the shape of distributions used in this paper.

We have documented the distribution of spike rates, gains, and weights for different types of neurons ( Figure 1– Figure 8). The distribution in all cases follows a lognormal shape. In some cases, we had data on the variability of spike rates and analyzed them for dispersion (CV, FF) under behavioral stimulation. While the fold-change from low spiking neurons to high spiking neurons is high, 5- to 7-fold, the variability for each neuron is comparatively low. It also seems to be adequately described by a percentage change over the whole population. This means that a low spiking neuron never reaches the same rate as a high spiking neuron, even when fully activated.

The similarities across neural systems are striking. For instance, in a midbrain nucleus (inferior colliculus) which is essentially an ’output’ site for auditory and somatosensory cortex, spike rates are high overall ²² , nonetheless the distribution of mean spike rate, and the variability are comparable to cortical data ( Figure 3). The distribution of mean spike rate is essentially the same under spontaneous and under behavioral conditions.

Lognormal rate distributions appear to be an essential property of neural tissues that occur in areas with very different neuron types and connectivity, and different absolute spike frequencies. They are present during spontaneous activity, and under activation of a network, in vivo as well as in vitro. They have a counterpart in a lognormal distribution of intrinsic excitability, and lognormal synaptic connectivity. This type of distribution seems to be an essential component of the functional structure of a mature network, which is not affected by learning, plasticity, or processing of information.

A lognormal distribution is characterized by parameters µ*and σ*. µ* = e ^µ is the median, a scale parameter, which determines the height of the distribution. σ* = e ^σ is the multiplicative standard deviation, a shape parameter which determines the width of the distribution. For distributions with small σ*(approximately σ* < 1.2 or σ < 0.182) a lognormal distribution is essentially identical to a normal distribution. (The coefficient of variation CV ∼ σ* − 1, so that for CV < 0.18, a lognormal equals normal distribution.) We collected data on spike rate, gain and synaptic weight distribution for a number of tissues in different experimental conditions ( Table 1– Table 3). For the height of the spike rate distribution, there are known differences, e.g. with lower values for cortex ( µ* ≈ 4.5) and higher values for Purkinje cell ( µ* ≈ 30) and midbrain nuclei (cf. Table 1). In other words, spike rates differ between brain areas such as cortex and cerebellum by a factor of 10.

In contrast to that, the width of spike rate distributions is more similar, with an average at σ* ≈ 2.2, with one outlier. The gain has a smaller σ*, i.e. a more normal, less heavy-tailed distribution than the spike rate. Minus the outlier (3.46), the mean for σ* is only 1.86, considerably lower than the width of the spike rate distribution ( Table 2). For weight distributions ( Table 3), the width σ* is consistently larger, with an average of almost 3 (2.91). The synaptic strength ( µ*) varies over at least one order of magnitude between cortex and cerebellum.

It turns out that σ* values are significantly different for rates, gains and weights, lowest for gains ( σ* ≈ 1.8), higher for rates ( σ* ≈ 2.2) and highest,( σ* ≈ 3), for weights. The data that we have are not precise enough to draw quantitative conclusions, but no large distinctions are apparent between the tissues ( Table 1– Table 3). We use a generic neural network to recreate lognormal distributions by adaptation rules and we will also show that distribution widths are structural properties which follow from general network properties.

Since not only mean spike rates but also both components, intrinsic excitability and synaptic weights, have lognormal distribution, this raises the question of how the functional system that we observe is generated. It is obvious, if the data are accurate, that these are basic parameters of any simulation and need to be reproduced in any model to make it biologically realistic.

We set up a generic neural network model (cf. Section 2.2) to explore the mechanisms of generating and maintaining rate, weight and gain distributions. The model consists of a source neuron group I, a target group J, a population of inhibitory neurons H, which are connected with J, and potentially recurrent excitation in the target group J. The spike rate distribution R ^I acts through a weight distribution W onto a gain distribution G, where inhibition H is subtracted, and a spike rate output distribution R ^J is produced ( Figure 9).

In the simplest case, we look at two sets of neurons, the source and the target. The source sends excitatory connections to the target, and exhibits variable weights at outgoing synapses. The input that a target neuron receives is fed through a linear filter G to produce an output rate R ^J according to Equation (1). The distribution for R ^J depends on G and W as well as on R ^I . The system is sufficient for calculations on the shape of distributions, as well as the effects of Hebbian and homeostatic plasticity.

We have explored the dependencies between gain, weight and rate distributions in simulations. First, we found that the width of the output spike rate distribution R ^J depends heavily on the gain distribution, but only slightly on the input weight distribution ( Figure 10). It does depend on the overall connectivity C, where σ R J * is wider for lower connectivity, but not very much ( Figure 10). Secondly, the width of the output distribution R ^J does not depend on R ^I or R ^H either ( Figure 11). The most important factor for a spike rate distribution remains the gain σ G * .

We may now ask, where do lognormal spike rate distributions come from? How is the system set up, i.e. what rules of adaptation generate lognormal distributions in weights and gains?

In the case of cortical networks, there are excitatory recurrent interactions that constitute a significant part of total input. In the case of cerebellar or striatal neurons, there are no recurrent excitatory interactions, only inhibitory interneurons and excitatory input. The generation of lognormal distributions must therefore be independent of recurrent excitation. It requires a system where continuous input shapes the weights and gains of a target network J. We start with the system that we described before, with random assignment of weights and gains. We employ adaptivity for weights, and also for gains, by positive Hebbian learning, or by negative homeostatic learning. The output of I is fed into J, and W and G are adaptive. Additionally, J may have excitatory recurrent connectivity, and learning takes place within the network J.

From any given initial spike rate distribution (Gaussian, uniform, lognormal) for I, we calculate W assuming a positive learning (Hebbian) adjustment rule, which is dependent on input and output frequencies. Each individual weight w _ij is updated by

w i j ′ = w i j + λ w i j ( r i I r j J − μ )

We use parameters λ and µ, such that the generated spike rate output R ^J is compatible in strength with the input rate R ^I .

Using Hebbian learning, we generate a weight distribution W that is lognormally distributed, independent of the initial configuration or the distribution of the gains in the system ( Figure 12). The lognormal distribution also develops independently of the rate distribution of the inputs, it only develops faster with lognormal rather than normally distributed spike rate input (not shown). It makes no difference whether we use a recurrent system J, or a non-recurrent population J with input from a population I with a spike rate distribution, as long as we use a Hebbian weight adaptation rule. For the shape of the distribution, it also does not matter whether we route the output of J back to I, or whether we use local or no recurrence. To show the effect of the adaptation rule, we also used homeostatic synaptic plasticity to adjust the weights. This means that the weight is adjusted inversely to the spike rate of input and output neurons.

w i j ′ = w i j − λ w i j ( r i I r j J − μ )

In this case, it is very clear that with any input or initial configuration and any gain distribution, only a normal distribution of weights results ( Figure 13). Again, a lognormal input spike rate slows the process of adaptation, but the end result is the same, a normal distribution.

Since gain distributions are also lognormal, we may ask in the same way how they develop and are maintained by plasticity rules. We adapt the linear gain G by either Hebbian or homeostatic learning. Each gain can be adjusted by a Hebbian rule

g j ′ = g j + λ g j ( r j J − μ )

or a homeostatic rule

g j ′ = g j − λ g j ( r j J − μ )

with parameters λ and µ.

We start with uniform or normally distributed G in an environment where W is lognormal, normal or uniform, and R ^I is normal or lognormal. If we adapt only G for any initial configuration, using any distribution for R ^I , including the lognormal distribution, and a lognormal or normal weight distribution, we arrive at a normal distribution for G with homeostatic learning and a lognormal distribution with Hebbian learning ( Figure 14).

Lognormal distributions develop from Hebbian plasticity, and homeostatic plasticity generates only normal distributions. The explanation lies in the nature of random statistical events, which generate normal distributions when the underlying mechanisms are sums of many small events, but lognormal distributions when the underlying mechanisms are multiplicative ³⁰ . We also wanted to understand the observed widths of the distributions. We hypothesized that the differences for σ* between W, R and G result from the network structure. Accordingly we started a simulation with initial uniform values for G and W and Hebbian update rules using the same learning rate λ for both ( Figure 15). We find that gain, rate and weight distributions match the experimental values, and that this is true for any tested constellation. We also found that Hebbian learning alone quickly escalates values, which develop exponentially, and that additional rounds of homeostatic adaptation are required to stabilize the system. Homeostatic learning pushes the system back towards a normal distribution.

Our data, in the most general way, allowing for various conditions and architectures, show that Hebbian learning is required both for intrinsic gain and for weights in order to generate the attested lognormal distributions. This is an interesting result, because it shows that we need prominent Hebbian intrinsic learning to explain the gain distributions that we find experimentally. Intrinsic learning is not just homeostatic adaptation, it follows the same rules as synaptic weight learning.

Spike rates of neurons seem to be universally distributed according to a lognormal distribution, with many neurons at low spike rates, and a small number at successively higher spike rates (heavy-tail) ¹ . The same distributions are found for synaptic weights ⁶ , and intrinsic properties associated with excitability (gain). The neurons that we reported on are of very different types, and they are embedded in different kinds of connectivity. Medium spiny neurons and Purkinje cells are GABAergic (inhibitory), while cortical and IC neurons are glutamatergic (excitatory), but this is not reflected in a distinct spike rate distribution. They also fire with very different average spike rates. IC neurons operate at very high frequencies, and Purkinje neurons at much higher frequencies than cortical or striatal projection neurons. But they all have the same spike rate distribution. It has been suggested ² that lognormal spike distributions are only a feature of cortical tissue and arise from the strong excitatory recurrent connectivity, but this is neither experimentally nor theoretically correct. While cortical pyramidal neurons exist in a heavily excitatory recurrent environment, medium spiny neurons, cerebellar Purkinje cells and IC neurons act mostly in a feed-forward environment, i.e. have no significant recurrent connectivity.

Beyond spike rate distribution, we also gathered data on weight and gain distributions. Again the observation of lognormal distributions is ubiquitous. We find synaptic weight distributions for cortex ⁶ and cerebellum that are lognormal, with characteristic width of distributions. For intrinsic properties, striatal projection neurons and cortical neurons ²³ show responses to constant current and current-to-threshold (gain) distributions, which again appear lognormally distributed, with smaller widths than spike rate distributions.

Our models show that lognormal distributions arise even in a purely input-output environment, and that they are a result of Hebbian learning of weights and gains, quite independent of the overall magnitude of the spike rates.

Mean spike rates, as well as intrinsic excitability and synaptic weights, have lognormal distributions.

It has often been assumed that variability in intrinsic excitability is a source of noise in neural computation ³¹ , even though others have argued that intrinsic variability contributes to neural coding ^{32, 33} and that intrinsic plasticity follows certain rules ³⁴ . An excellent overview of the experimentally attested forms of intrinsic plasticity is contained in 35, cf. 36– 38. Many other detailed observations are contained in 39– 42.

Recently, Mahon and Charpier ⁴ have shown that intrinsic excitability is stable in individual neurons under control conditions, while stimulation protocols (e.g. in barrel cortex of anesthetized rats) change intrinsic excitability by at least 50–100%. However, the conclusions drawn from the experimental research are often contradictory. Intrinsic plasticity is sometimes assumed to act in a negative, homeostatic way, i.e. opposite to synaptic plasticity ⁴ , but sometimes in a ’Hebbian’, positive way, i.e. cooperative with synaptic plasticity ^{41, 43} . There is evidence for (short-term) negative or homeostatic plasticity, which has been previously investigated ⁴ .

Our work has now shown that any kind of neural system with linear gains requires positive, Hebbian intrinsic plasticity to produce and maintain a lognormal distribution of gains. We also could show that the observed widths of the distributions, i.e. the differences for σ* between W, R and G, naturally result from the network structure and are built into the system simply by Hebbian adaptation.

Lognormal distributions may arise as stable properties of the system during early development (the set-up of the system), i.e. before actual pattern storage or event memory develops, and they are maintained during processing by a Hebbian type of positive adaptation events. Homeostatic plasticity consists in downregulating gains or weights with increases in firing rates. Purely homeostatic learning results in normal distributions, and erases existing lognormal distributions. By combining homeostatic and Hebbian adaptation we can achieve and maintain stable lognormal distributions.

A lognormal distribution means that values are normally distributed on a logarithmic scale. From an engineering perspective, basic Hebbian plasticity for synapses and intrinsic properties is sufficient to generate stable logarithmic distributions. If there is random variation of multiplicative events, as in Hebbian plasticity, a lognormal distribution will be the result ³⁰ .

This is related to principles of sensory coding, where logarithmic scale signal processing enhances perception of weak signals, while also being able to respond to large signals - effectively increasing the perceptual range compared to linear coding ¹⁴ . In an interconnected network logarithmic coding may turn into a property for the access of representations. Feature clusters, or event traces could be accessed by targeted connections to the top-level neurons, which then activate lower level neurons in their immediate vicinity. By accessing high frequency neurons preferentially, a whole feature area can be reached, and local diffusion will provide any additional computation. Similarly, the results of a local computation can be efficiently distributed by high frequency neurons to other areas. Fast point-to-point communication using only high frequency neurons may be sufficient for fast responses in many cases. Scale-free networks in general support synchronization, which is also a useful feature for rapid information transfer and access ⁴⁴ .

Recently, publications ^{45, 46} have shown that there is indeed a difference between high-frequency and low-frequency neurons in their connectivity: high-frequency neurons have short delays, strong connections, and directed targets, while low-frequency neurons have long delays, weak connections and diffuse targets.

The lognormal distribution of spike rates has significant implications for neural coding. Logarithmic spike rates are coupled with linear variance for responses to behavioral stimulation. In other words, the greatest part of the coding results already from the frequency rank of the neuron itself, such that high frequency neurons have the largest impact. A fixed mean rate for each neuron allows stable expectation values for network computations.

Logarithmic, hierarchical coding does not need to be sparse. The low frequency neurons may matter the most in terms of input response. With lognormal synaptic weight distributions, if strong synapses are kept stable, they may transmit an input neuron’s mean firing rate to targets and in this way provide stability to the system. All other synapses could be arbitrary. This would allow for continued pattern learning to be implemented by the bulk of low weight synapses, while the framework of neuronal interactions, e.g., the ensemble structure, could be unchanged. Such a division of labor between strong synapses and weaker ones could have many advantages in a complex, modular network.

Experimental data have often shown that sampling of neuronal responses from a large population (10 ⁵ or more neurons), which become activated at 30% or more, yields accuracy for a stimulus already for small samples (100–200 neurons or 1–2%) (e.g., 47). We may suggest that this happens when we sample from a highly modular structure, and we have been able to replicate the effect with lognormal networks ⁴⁸ .

In our earlier work ¹ , we found that intrinsic excitability manifested by spike response to current injection and rheobase in vitro for dorsal striatal and nucleus accumbens neurons seems to have the same distribution as the firing rate in cortex under in vivo conditions. Approximately at the same time ⁶ , had observed a heavy-tailed distribution of synaptic weights in cortical tissue.

In this article, we have done three things: (a) collected data to show that rate, weight and gain distributions in different brain areas all follow a heavy-tailed, specifically a lognormal distribution; (b) created a generic neural network model to show that these distributions arise from Hebbian learning, and specifically that intrinsic plasticity must be Hebbian as well; and (c) shown that the width of the distributions, as experimentally attested, arise naturally from the network structure and the role of its components, in a very robust way. We have also discussed what the lognormal distribution means for neural coding: a division of labor between fast transmission by high-frequency neurons and low-level computation by low-frequency neurons in a modular structure, and possibly a division of labor between stable components (strong synapses, high-frequency neurons) and more variable components (weak synapses, low-frequency neurons).

The GNN simulation software was programmed in Matlab, and is available in GitHub: https://github.com/ gscheler/GNN/tree/v0.1

Archived source code as at time of publication: https://doi.org/10.5281/zenodo.829949 ⁴⁹

OSS approved license: Apache 2.0.

All the data required for re-analysis of the study have been referenced throughout the manuscript.