Randomly connected networks generate emergent selectivity and predict decoding properties of large populations of neurons

Advances in neural recording methods enable sampling from populations of thousands of neurons during the performance of behavioral tasks, raising the question of how recorded activity relates to the theoretical models of computations underlying performance. In the context of decision making in rodents, patterns of functional connectivity between choice-selective cortical neurons, as well as broadly distributed choice information in both excitatory and inhibitory populations, were recently reported [1]. The straightforward interpretation of these data suggests a mechanism relying on specific patterns of anatomical connectivity to achieve selective pools of inhibitory as well as excitatory neurons. We investigate an alternative mechanism for the emergence of these experimental observations using a computational approach. We find that a randomly connected network of excitatory and inhibitory neurons generates single-cell selectivity, patterns of pairwise correlations, and indistinguishable excitatory and inhibitory readout weight distributions, as observed in recorded neural populations. Further, we make the readily verifiable experimental predictions that, for this type of evidence accumulation task, there are no anatomically defined sub-populations of neurons representing choice, and that choice preference of a particular neuron changes with the details of the task. This work suggests that distributed stimulus selectivity and patterns of functional organization in population codes could be emergent properties of randomly connected networks.


Introduction
With the deluge of neural recordings made possible by modern recording methods [2], theoretical neuroscience must address the challenge of relating complex activity patterns to the algorithms thought to underlie brain function [3,4]. For example, evidence accumulation tasks explore how the brain makes decisions based on the temporal integration of incoming sensory information. One class of models for performing this discrimination is that of attractor networks. In an attractor network model of decision making, pools of neurons fire selectively for a particular choice. A transient activation is prolonged through slow recurrent excitation within the pool while inhibiting other pools of neurons that are selective in their firing for other choices through non-specific it is unclear how the specific pattern of connectivity would be generated. In this paper, we explore an alternative mechanism for these experimental observations.
More specifically, in this task [1], a rodent is presented with an irregular train of either visual or auditory impulses and must determine whether the average frequency with which those pulses arrived is above or below some internally remembered threshold.
Recordings of population activity in the posterior parietal cortex (PPC) during the task revealed weak choice selectivity in single cells, with a fraction of individual cells showing significant selectivity for one of the two choices. A linear classifier operating on the activity across the population decoded the choice with high accuracy. Both excitatory and inhibitory neurons were selective for choice, and noise correlations between pairs of neurons reflected whether stimulus selectivity was shared or opposing. A straightforward mechanism for these observations is that some specific pattern of connectivity exists in the cortical network that separates inhibitory cells into "pools" that are selective for specific choices.
We explore an alternative mechanism for the emergence of these observations. We hypothesize that a randomly connected network can produce patterns of activity that are sufficiently distinct to differentiate input frequency. The eigenvalue spectrum of a randomly drawn connectivity matrix often has a tail of eigenvalues with large real parts [9], and this remains true even when networks follow Dale's principle of separate excitatory and inhibitory neurons [10,11]. We reason that input modes overlapping with fast-growing modes are amplified by network dynamics, and nonlinear input-output function at the level of single cells could result in an the activation pattern that depends on the temporal frequency of the input. Thus, the network would produce different patterns of activation as a function of the input frequency, resulting in emergent selectivity for frequency at the level of the population. We are primarily interested in whether such a generic random network of firing rate units could generate the observed patters of selectivity, functional connectivity, and cell-type-specific readout weight September 24, 2019 3/25 distributions simply through randomly arising heterogeneity in synaptic inputs at the level of single neurons. A thorough theoretical derivation of selectivity and readout weight distributions would be of interest but rather complex, so instead we explore this question computationally.
Our simulation results support this hypothesis. Our main findings are (i) that heterogeneity in connectivity generates differences in inputs to single cells that are dependent on stimulus frequency and (ii) that these differences are sufficient to distinguish between low-and high-frequency inputs. Our model reproduces experimental findings, including the distribution of single-neuron selectivity, patterns of pairwise noise correlations, the performance of a classifier, and the distribution of readout weights. Our theory makes the verifiable experimental predictions that, if the mechanism is through emergence rather than specific connectivity, then (1) there is no anatomical basis for sub-populations tuned to a particular choice and (2) when task parameters, such as input frequency are changed, neural selectivity also changes. These results suggest a mechanism for how cortical networks could exhibit functional organization without specific patterns of cortical connectivity.

Simulation of a temporal evidence accumulation task
We first review key experimental results for a particular version of an evidence accumulation task that can be performed by rodents [1,6,7]. In these tasks, rodents are trained to discriminate low-and high-frequency inputs (pulses, either visual or auditory) and to report their choice. The exact timing of the inputs is random, so the rodent is not reporting the inter-pulse interval but deciding based on the stimulus history over a short period, typically one second. They do this accurately. Recordings of neurons in PPC during the task show a distribution of selectivity for choice, meaning that the reported outcome of the evidence accumulation is predictable from the activity of some fraction of the neurons. A classifier trained to determine choice from the recorded population activity performed accurately. Interestingly, it performed equally well for excitatory or inhibitory sub-populations, and an analysis of the classifier weights showed that excitatory and inhibitory cell activities were weighted similarly [1].
Here we study a randomly connected rate-based network of N E excitatory cells and N I inhibitory cells that performed the evidence accumulation task described above (Fig. 1A). While similar in some aspects, this network is not set up as required for reservoir computing [12,13], in that our network operates in a regime in which spontaneous activity is low in the absence of inputs. The firing rate of individual units is the sum of any external input, representing the stimulus, and synaptic inputs from the rest of the network, passed through a saturating non-linearity so that activity is between 0 and 1. Connectivity in the network is sparse, with a 20% probability of connection between two cells. Excitatory and inhibitory synaptic weights are drawn from truncated normal distributions (∼ N (g E,I , σ 2 )) ( Fig. 1B, Methods), and excitatory and inhibitory synaptic weights to neurons are balanced on average, such that N E g E − N I g I = 0.
These are set to generate a broad spectrum of eigenvalues in the connectivity matrix [10]. We fix g E , g I and σ, and we focus on a single combination of parameters, but the essential results are not dependent on making these particular choices (see S2 This network was conceived as a model of the posterior parietal cortex (PPC) of rodents, which does not receive direct sensory inputs but rather receives inputs that have passed through upstream networks [14]. Moreover, the average population firing rates in PPC during such a task are not directly related to stimulus frequency [6]. To account for this effect in our simulation, we scaled the amplitude of inputs such that the network firing rate is, on average, equal for each frequency ( Fig. 1D and S1 Fig). From a computational standpoint, this choice makes the task of decoding from the network activity more difficult, as the average firing rate does not encode stimulus identity.
We present the simulation results by following the experimental observations presented by Najafi and colleagues [1]. First, we analyze simulated choice selectivity at the level of single cells and pairs of cells. Given that we know the network connectivity, we analyze both noise correlations and the underlying connectivity pattern. Next, using  Overview: Simulation of an input frequency discrimination task in a recurrent network of excitatory and inhibitory neurons with random connectivity. A: Recurrent E-I network of randomly connected neurons. Synaptic weights are heterogeneous (line thickness), and connectivity in the network is random. Inputs to the network have the same pattern of spatial activation but different temporal features. Readouts of network activity (from all neurons, or restricted to excitatory or inhibitory neurons) discriminate between stimulus categories. We then analyze the distributions of readout weights. B: Distributions of non-zero inhibitory (shaded red) and excitatory (shaded blue) synaptic weights. 20% of the weights are non-zero. C: Inputs to the network are pulses arriving at random times with low (top row, yellow) or high (bottom row, red) average frequency. D: Histogram of population average network activity over trials for low-and high-frequency stimuli. Average network activity is not informative of input categories.

Emergent selectivity for stimulus category in single cells
We first examine choice selectivity in single cells from the network simulations. For this analysis, choice was defined to be the correct stimulus label (i.e., low vs. high frequency of input pulses) on each trial, which is equivalent to analyzing the correct trials only in a behavioral experiment. Single-neuron activity was averaged over time, yielding a single number per trial for each cell. An ideal observer analysis was used to discriminate between low-and high-frequency inputs (see Methods). Choice selectivity [1] was defined as the area under the receiver-operator curve (AUC, [15]), which will be less (greater) than 0.5 when the cell is selective for the low-frequency (high-frequency) stimulus. For each network realization, this generated a distribution of AUC values across all cells ( Fig. 2A). To assess the significance of a single AUC calculation, we computed the 95% confidence interval of AUC values obtained with shuffled trial labels.
Values exceeding these bounds are significant. Approximately 30% of excitatory and 30% of inhibitory cells in the example network in (

Noise correlations reflect relative selectivity of pairs of cells
We next asked whether this simple, unstructured network also explained pairwise relationships observed in the experimental recordings. Specifically, pairs of cells in PPC that shared the same selectivity had higher noise correlations [1] than pairs that had opposite selectivity. To compare this to our simulation results, we computed noise correlation as the cross-correlation between neurons of stimulus-specific responses, averaged across stimuli. Because the only source of noise in the simulation, the variable timing of inputs, is a shared input to all neurons, we expect noise correlation in the model to be higher in the model than in the data, but cells with the same selectivity are expected to have higher noise correlation. Restricting analysis to the cells that exceeded the significance criterion for AUC values, we categorized cells by selectivity and compared average noise correlation between pairs of cells with same and opposite preference. As observed experimentally, noise correlation in the simulation was higher between pairs of cells with the same selectivity than between cells with opposite selectivity (Fig. 2D). Thus, organization of functional connectivity in the network emerged without setting up distinct clusters of connections in the network.
We further examined the patterns of connectivity between these sets of cells (Fig. 2E). The overall probability of a synapse was set to 20% for all simulations. There was a nearly identical probability of connection from an excitatory cell to an inhibitory cell with the same selectivity, 20% ± 1% (SD), as with opposite selectivity, 19% ± 3% (SD). The inhibitory to excitatory connection between cells with the same selectivity was similar as well, 18% ± 2% (SD). Among cells with opposite selectivity, the probability of connection from an inhibitory to an excitatory cell was 24% ± 4% (SD).
Even this small amount of excess connectivity between inhibitory and excitatory cells of opposite selectivity was sufficient to reproduce the experimentally observed trends in functional connectivity and stimulus selectivity patterns. We emphasize that this bias in connectivity was not put in by hand, but rather was uncovered by the dynamics that it shaped.

Decoding from population activity
To determine how accurately the simulated network represented choice, we trained linear classifiers to discriminate between low and high stimulus frequency. We fit a classifier using activity near the end of the stimulus period (circles, Fig. 3A) and tested the classifier over the full stimulus period on a set of reserved trials (Methods). In this September 24, 2019 9/25 simulation, classifier accuracy reached maximal performance within 5τ (time constant of the network) and decoded with 83% ± 2% accuracy over the last half of the stimulus period. A classifier fit using only the activity of the inhibitory cells performed with 76% ± 2% accuracy, and a randomly drawn subset of excitatory cells equal in number to the number of inhibitory cells decoded with 78% ± 2% accuracy. Across all instances of the network, the accuracy of decoders was comparable for inhibitory cells and excitatory ones (Fig. 3B). The range of performance we observed across randomly drawn networks (69% ± 2% to 86% ± 2%) was highly consistent with the population decoding accuracy observed in experiments (about 70% to 85%, [1]).  significantly different (Fig. 3C). In all simulated networks, the distributions of weights for excitatory and inhibitory cells were overlapping (Fig. 3D-E)  either gained or lost selectivity as the task parameters changed (Fig. 4C). The fraction of all cells that changed selectivity (in either direction) varied across networks but was always greater than zero, averaging 22% of cells (+/-9%, range 8% to 42%). To summarize, we predict that the set of selective cells depends on the temporal features of the input and will change if the task changes.

Discussion
We presented a set of simulation results that account for several key features of population recordings in PPC during an evidence accumulation task. This simulated patterns of activity, we measure single-cell selectivity for input rate, patterns of noise correlations between same and opposite selectivity cells, and readout performance and readout weight distributions. We find that these measures are highly consistent with those observed experimentally. Importantly, our simulations do not include specific patterns of connectivity between excitatory and inhibitory neurons; any biases that emerge are the result of dynamics shaped by random heterogeneity in connectivity patterns. We suspect this is a fairly generic effect of a network with a broad spectrum of eigenvalues, such as the random networks that we studied, but the theoretical connection between the spectra of a connectivity matrix in a nonlinear network and the emergence of selectivity and distribution of readout weights among excitatory and inhibitory populations remains to be explored.
Experimental work has shown that PPC is integral to performing a sensory evidence integration task [6,7]. showed that an unstructured network produces all of these effects as well.
Additionally, previous results have argued that, while the posterior parietal cortex is involved in performing a visual sensory evidence integration task [6,7] in rodents, PPC does not itself integrate evidence [4,7]. One possibility, consistent with our simulations, is that the population response is subtly distinct across input frequencies, and these distinctions are learned through a reinforcement mechanism in some other area, which then feeds back to PPC to enhance the distinction between sensory inputs. This feedback mechanism may interact with other biases, such as motivation and trial history [16,17], which modulate choice. Such a feedback loop could further enhance apparent selectivity in connectivity in PPC because it emphasizes dynamics that were shaped by heterogeneity in connectivity.
Finally, we found that stimulus information could be decoded very early in the stimulus window. The technical reason for this is that due to the input scaling, the first pulse of the input series carried stimulus information. Scaling was used to make the decoding task more challenging and to reproduce the experimental observation that PPC responses to low-frequency inputs are not lower than the PPC responses to high-frequency inputs. We did not implement more realistic adaptation dynamics, which could be the mechanism underlying such scaling. One would expect that such adaptation mechanisms (e.g., short-term depression or facilitation) could be transformed into a population code in a spiking network, readily decoded downstream [18], and it is an interesting question for future study whether spiking networks such as these replicate the experimentally observable quantities that we focused on here. Adding this feature to the model would slow the ramp-up in decoding accuracy, but would also require more choices about adaptation rates and tuning. Moreover, such elaboration of the model is superfluous to answering the question of whether a heterogeneous network of excitatory and inhibitory units can distinguish between inputs with different temporal frequencies and produce statistical features comparable to experimental observations. instance, spatially distinct neural representations could trigger distinct neural trajectories [24], matching the spatio-temporal multineuronal dynamics on single trials more closely. We did not pursue this here, as our goal was to show that heterogeneity in network connectivity could explain many features of population recordings during a simple discrimination task.
One of our key results is that selectivity can emerge for task parameters from an unstructured, random network. Several theoretical studies have previously examined the emergence of strong selectivity in random networks. For example, it has long been recognized that orientation selectivity in primary visual cortex could emerge from random projections from geniculate inputs [25,26] and, in balanced state networks, this selectivity would be robust due to the dynamic cancellation of non-selective inputs [23,[27][28][29]. Moreover, selectivity from random projections could be enhanced through learning mechanisms either in the feed-forward projections [25] or in the recurrent cortical network [23]. The mechanism for selectivity presented here adds to these by demonstrating selectivity for parameters that reflect temporal, rather than spatial, patterns. In this model, there is no spatial heterogeneity in the inputs, and the source of selectivity is not from heterogeneity in feed-forward projections, but from the dynamics of a recurrently connected non-linear network.
Finally, our theory makes the following prediction, which should be verifiable experimentally. Suppose, the decision boundary for reporting low-versus high-frequency stimuli changed. If the network is structured as separate pools of excitatory and inhibitory neurons representing choice for one or the other stimulus category, then the representation of choice in PPC will not change with the task. If instead functional organization is generated by emergent network properties, when the task changes, the September 24, 2019 14/25 selectivity of individual cells will shift, as different pools of neurons represent the lowversus high-frequency stimuli. These pools of neurons would be overlapping, as the frequencies being discriminated change.
As ever larger populations of neurons are simultaneously recorded, and experiments frequently focus on a variants of well-controlled sensory discrimination tasks, we face a tremendous challenge in inferring mechanism from observations. Very generally, there are two mechanisms that could account for diverse and mixed selectivity along with patterned functional connectivity across a population of neurons engaged in some experimental task. The first is that the network is wired specifically to achieve this, and if that is the case, then one must also explain the developmental or learning process that produced such intricate topology in the network. The alternative, which we explored here, is that the observed patterns of selectivity can be explained as an emergent phenomena from simple patterns of statistical connectivity. In the particular case examined here, we were able to reproduce distributions of selectivity, functional connectivity, and population readout weight patterns that were experimentally observed.
Even though we analyzed the emergence of selectivity in a specific experiment, we believe that similar conclusions could hold in other applications, in which a broad distribution of selectivity is observed.

Random recurrent firing rate network
We simulated a recurrent neural network of N = 500 neurons (N E = 400 excitatory neurons and N I = 100 inhibitory neurons) using standard firing rate equations: where x i is the "membrane potential" of neuron i and r i is its firing rate, obtained from

The frequency discrimination task
In the simulated task, the network was driven by stimuli consisting of irregular impulses with average frequency f (Fig. 1A). The length of the sampling period in the simulated September 24, 2019 16/25 task was 50τ , corresponding to an effective τ of 20 ms. In the first set of simulations, either f = 8 Hz or f = 16 Hz. In simulations in Fig. 4, we analyzed f = 10 Hz to f = 20 Hz.
For each trial with average frequency f , the input times (t k ) were selected randomly with uniform probability from the stimulus interval (50τ , or 1000ms) by drawing a fixed number of time points (e.g., 8). Temporal precision of impulse timing was 0.01τ , so impulse times were drawn from the integers 1 to 5000 without replacement. For simulated trials with the same impulse frequency, the trial-by-trial variation in input times was the only source of variability across trials.
We assumed that upstream sensory networks filtered the pulse inputs, so the overall input current to the network was described by where a is a filtering timescale of pre-processing networks and the summation was taken over all t k < t. We set a = 0.5 for all simulations (full width at half max of a single pulse is 1.7τ ).
Input currents were scaled by a frequency-dependent factor α to match the average firing rate across the network between conditions (Fig. 1, S1 Fig). To implement this scaling, 25 sets of input parameters (frequency and amplitude) were simulated with input frequencies from 0.1 to 0.5 (per τ ), and amplitudes from 0.3 to 15, a range spanning parameters that generated a range of firing rates in the network. From each of these 25 simulations, we computed the average network firing rate, and we interpolated this surface to find contours of equal firing rate. For the simulated experiment, we used combinations of frequency and amplitude that fall on a fixed contour. Across networks, the typical amplitude ratio between low-and high-frequency inputs was 2.1. We verified post-simulation that the average firing rates match across frequency conditions (see, e.g., Fig. 1D).

Noise correlations
We computed noise correlations from the activity at the end of the stimulus period by subtracting the stimulus-averaged response and then computing neuron-neuron correlation coefficients.

Classifier analysis
The goal of the classification was to discriminate between two input frequencies using the simulated activity patterns. Simulated network activity was temporally down-sampled by averaging over time windows of size 1 (τ ). Neurons that received direct inputs were excluded from the decoding analysis, leaving 320 excitatory neurons and 100 inhibitory neurons. We fit classifiers separately on the full population ("all," 420 neurons), a subset of excitatory neurons ("exc-sub", 100) and a subset of inhibitory neurons ("inh", 100). We split the 800 simulated (400 in each condition) into "test" and "training" sets. We trained a classifier (linear kernel, SVM) to predict the stimulus label based on the activity (in "all," "exc-sub" and "inh") at each time point. We trained the classifier on the z-scored activity from each cell [7]: The classifier finds a rule ξ, η ξ * z > η.
September 24, 2019 18/25 We also calculated the weights (w) and bias (b) that operate on firing rates directly Classifier accuracy was calculated on the reserved test set. In Fig. 3A, a classifier was fit at a single time point at the end of the stimulus window and tested at all other time points. Uncertainty in classifier accuracy was estimated by fitting the classifier using different cuts of the data: train/cross-val/test sets were drawn randomly, a classifier fit, and weights and accuracy on the test set recorded. This was repeated 50 times, and the error bar on classifier accuracy is the standard deviation across test set accuracy generated in this way. Lines show contours at fixed average firing rate. For each firing rate (indicated by line color), we extract amplitudes and frequencies on the corresponding contour. B: Amplitudes for the high-frequency (freq2=2*freq1) input, plotted against the lower frequency (freq1). C: Amplitudes for the low-frequency input, plotted against the lower frequency (freq1). The highest firing rate contour that is defined over freq1 = 0.14 to 0.2 is simulated. Units of frequency are per τ ; multiplication by 50 converts to Hz.  Figure. Additional parameter combinations. Fixing network topology (i.e., which elements of J are non-zero), we simulated three networks: with original weights (i), with all synaptic weights scaled by a factor of 1.5, (ii) and with homogeneous excitatory and homogeneous inhibitory synaptic weights (iii). A: Image of network connectivity for 25 (of 500 total) neurons showing that the topology was kept the same for each simulation. B: Distribution of non-zero excitatory and inhibitory weights in the network. Note that there are approximately four times as many excitatory weights, but they are on average a quarter of the strength of inhibitory weights. C: Histogram of stimulus 1 (yellow) and stimulus 2 (purple) population firing rates for each parameter scaling. Firing rates are matched for each simulation individually; these are operating over a similar population firing rate range. D: Decoding accuracy of the full population in each network. E: Decoding accuracy of excitatory and inhibitory cells in each network. D and E show that decoding accuracy persists after a drastic parameter change, for this network topology. F: Decoding performance for all simulated networks.