# Random Compressed Coding with Neurons

Simone Blanco Malerba, Mirko Pieropan, Yoram Burak, Rava Azeredo da Silveira

► To cite this version:

Simone Blanco Malerba, Mirko Pieropan, Yoram Burak, Rava Azeredo da Silveira. Random Compressed Coding with Neurons. Cell Reports, 2025, 44 (3), pp.115412. 10.1101/2022.01.06.475186. hal-04305609v2

HAL Id: hal-04305609

https://hal.science/hal-04305609v2

Submitted on 10 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
**Cell Reports**

# Random compressed coding with neurons

## Graphical abstract

<figure><img src="image_2.png" /></figure>

## Authors

Simone Blanco Malerba, Mirko Pieropan, Yoram Burak, Rava Azeredo da Silveira

## Correspondence

rava@iob.ch

## In brief

Blanco Malerba et al. demonstrate that extended, irregular tuning curves, observed in the brain, efficiently compress information in a low-dimensional representation. Optimal smoothness of tuning curves that balances local and global errors yields an exponentially precise neural code in population size. This coding scheme is illustrated in monkey motor cortex.

## Highlights

* Extended irregular tuning curves compress stimulus information in a small neural population

* Optimal compression yields a tuning-curve smoothness that balances local and global errors

* This efficient code is achieved in a simple feedforward neural network with random weights

* Tuning curves in monkey motor cortex exhibit properties of highly efficient neural coding

<figure><img src="image_3.png" /></figure>
Article

# Random compressed coding with neurons

Simone Blanco Malerba,<fnref n="1" />,<fnref n="2" /> Mirko Pieropan,<fnref n="1" /> Yoram Burak,<fnref n="3" />,<fnref n="4" /> and Rava Azeredo da Silveira<fnref n="1" />,<fnref n="5" />,<fnref n="6" />,<fnref n="7" />,<fnref n="8" />,*

## SUMMARY

Classical models of efficient coding in neurons assume simple mean responses—"tuning curves"— such as bell-shaped or monotonic functions of a stimulus feature. Real neurons, however, can be more complex: grid cells, for example, exhibit periodic responses that impart the neural population code with high accuracy. But do highly accurate codes require fine-tuning of the response properties? We address this question with the use of a simple model: a population of neurons with random, spatially extended, and irregular tuning curves. Irregularity enhances the local resolution of the code but gives rise to catastrophic, global errors. For optimal smoothness of the tuning curves, when local and global errors balance out, the neural population compresses information about a continuous stimulus into a low-dimensional representation, and the resulting distributed code achieves exponential accuracy. An analysis of recordings from monkey motor cortex points to such "compressed efficient coding." Efficient codes do not require a finely tuned design—they emerge robustly from irregularity or randomness.

## INTRODUCTION

Neurons convey information about the physical world by modulating their responses as a function of parameters of sensory stimuli. Classically, the mean neural response to a stimulus—referred to as the neuron's tuning curve—is often described as a smooth function of a stimulus parameter with a simple monotonic or unimodal form.<fnref n="1" />⁻<fnref n="6" /> The deviation from the mean response—the “neural noise”—may lead to ambiguity in the identity or value of the encoded stimulus, and the coding performance of a population of neurons as a whole is dictated by the forms of the tuning curves and the joint neural noise. In the study of population codes, the efficient-coding hypothesis has served as a theoretical organizing principle. It posits that tuning curves are arranged in such a way as to achieve the most accurate coding possible given a constraint on the neural resources engaged.<fnref n="7" />⁻⁹ The latter is often interpreted as a metabolic constraint on the maximum firing rate of a single neuron or on the mean firing rate of the whole population.¹⁰⁻¹²

In order to tackle this constrained optimization problem in practice, tuning curves are parametrized and the corresponding parameters are optimized. A large body of literature addresses this constrained-optimization problem, in particular in the perceptual domain. For example, many studies model tuning

curves as Gaussian or other bell-shaped functions and obtain the values of their centers and widths that minimize the “perceptual” error committed when information is decoded from the activity of a population of model neurons.¹⁰,¹³⁻¹⁶ In the resulting optimal populations, and if noise among neurons is independent, the coding error typically scales like $1/\sqrt{N}$, where N is the number of model neurons.¹⁷ This behavior can be intuited based on the observation that the “signal” in the neural population grows in proportion to N while the noise grows as $\sqrt{N}$, yielding a signal-to-noise ratio that increases in proportion to the square root of the population size. (In some models of population neural coding of a one-dimensional parameter, the width of tuning curves can be further optimized to yield an additional factor of $1/\sqrt{N}$; the error then scales like $1/N^{18,19}$.)

Real neurons, however, can come with much more complex tuning curves than simple Gaussian or bell-shaped ones. Grid cells recorded in the entorhinal cortex offer a salient example.²⁰⁻²³ Their tuning curves in two-dimensional, open-field environments are multimodal and periodic as a function of spatial coordinates. When animal motion occurs in three dimensions, grid cells also exhibit irregular, locally periodic but globally aperiodic firing fields.²⁴,²⁵ An appreciable fraction of cells in the superficial layers of the medial entorhinal cortex are spatially selective, but, unlike grid cells, their firing fields are highly irregular

<figure><img src="image_4.png" /></figure>

<footnote n="1">¹Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, 75005 Paris, France</footnote>

<footnote n="2">²Institute for Neural Information Processing, Center for Molecular Neurobiology, University Medical Center Hamburg-Eppendorf, 20251 Hamburg, Germany</footnote>

<footnote n="3">³Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 9190401, Israel</footnote>

<footnote n="4">⁴Edmond and Lily Safra Center for Brain Sciences, Hebrew University of Jerusalem, Jerusalem 9190401, Israel</footnote>

<footnote n="5">⁵Institute of Molecular and Clinical Ophthalmology Basel, 4031 Basel, Switzerland</footnote>

<footnote n="6">⁶Faculty of Science, University of Basel, 4056 Basel, Switzerland</footnote>

<footnote n="7">⁷Department of Economics, University of Zurich, 8001 Zurich, Switzerland</footnote>

<footnote n="8">⁸Lead contact</footnote>

<footnote>*Correspondence: rava@iob.ch</footnote>

<footnote>https://doi.org/10.1016/j.celrep.2025.115412</footnote>
even during motion in two dimensions.²⁶,²⁷ Recent measure-
ments indicate that place cells in large environments are also
characterized by complex tuning curves, with multiple place
fields that are irregular both in location and in size.²⁸,²⁹
(By contrast, “classical” place cells exhibit only one or two place
fields.) Multimodal and irregular frequency selectivity was
observed in auditory neurons in mice and ferrets,³⁰ and a number
of other examples of irregular tuning curves were observed in
other cortical regions and in different species.³¹⁻³³ The growing
body of observations of complex tuning curves calls for a sys-
tematic theoretical understanding of their implications for cod-
ing. In the context of grid cells, the periodicity of their tuning
curves, as well as their functional organization in modules, was
shown to yield a spatial representation with an accuracy that
scaled exponentially (rather than algebraically, as above) in the
number of neurons.³⁴⁻³⁸

Here, we ask whether such strong boost in the efficiency of a
population code must rely on finely tuned properties, such as
the tuning curves' periodicity and the organization in modules,
or, alternatively, can arise more generally in populations of neurons
with complex, but irregular, tuning curves. We approach this ques-
tion by studying the case of a random neural code: a population
code that relies on tuning curves sampled from a stochastic pro-
cess. Tuning curves span the entirety of the stimulus space (i.e.,
they are spatially extended), and are characterized by their degree
of irregularity (or smoothness), parametrized by the correlation
length of the stochastic process. Both irregularity and spatial
extension of tuning curves lead to the emergence of two qualita-
tively distinct types of error—local errors, in which the encoding
of nearby stimuli is ambiguous, and global (or catastrophic) errors,
in which the identity of the stimulus is lost altogether. The efficient
coding problem then translates into a trade-off between these two
types of error: at the optimum, tuning curves have to be sufficiently
irregular to suppress local errors but smooth enough to prevent a
prohibitive rate of global errors. We show that, in the resulting
optimal population code, the coding error is suppressed exponen-
tially with the number of neurons in the population, allowing for an
accurate low-dimensional neural representation of a continuous
stimulus space, even in the presence of high-variance noise.

We first develop the theory for a one-dimensional input (e.g., a
spatial location along a line or an angle), then we generalize it to
higher-dimensional inputs. This generalization allows us to apply
our model to data from monkey motor cortex, where cells display
complex tuning curves.³³ We fit our model to the data and
discuss the merit of a “compressed code” that relies on irregular
tuning curves. Our results can be understood in terms of the ge-
ometry of the manifold of population activity: we relate the prop-
erties of individual neurons to the fine structure of the manifold,
which governs information processing. Overall, our study ex-
tends earlier theoretical work on grid cells and other finely de-
signed codes by proposing that efficient compression of infor-
mation can occur robustly even in the case of complex, but
irregular, tuning curves.

RESULTS

We organize the description of our results as follows. First, we
present, in geometric terms, the qualitative difference between

a code that uses simple, bell-shaped tuning curves and one
that uses more complex (periodic or irregular) ones. Second,
we introduce a simple model of neural population that can
interpolate between simple and complex tuning curves de-
pending on the values of its parameters. Third, we charac-
terize the accuracy of the neural code in the limiting case of
maximally irregular tuning curves. Fourth, we extend the dis-
cussion to the more general case in which an optimal code
is obtained from a trade-off between local and global errors.
All the above is done for the case of a one-dimensional input
space. Fifth, we generalize our approach to the case of a
multi-dimensional stimulus. This allows us, sixth, to apply
our model to recordings of motor neurons in monkey, and to
analyze the nature of population coding in that system. Sev-
enth, we give a quantitative description of the geometry of
the population response induced by our network as a function
of its parameters through a measure of dimensionality. In the
final sections, we consider various extensions of our approach
and how they impact the coding performance. These include
noise correlations, non-linearity in the response function of
second-layer neurons, and biologically motivated constraints
on the synaptic weight matrix.

The geometry of neural coding with simple vs. complex
tuning curves

A neural code is a mapping that associates given stimuli to a
probability distribution on neural population activity; in particular,
the code maps any given stimulus to a mean population activity.
In the case of a continuous, one-dimensional stimulus space, the
latter is mapped into a curve in the *N*-dimensional space of
the population activity, whose shape is dictated by the form of
the tuning curves of individual neurons. As an illustration, we
compare the cases of three neurons with Gaussian tuning curves
and three neurons with complex but regular (i.e., periodic, grid-
cell-like) tuning curves with three different periods (Figure 1A).
Simple tuning curves generate a smooth population-response
curve, implying that similar stimuli are mapped to nearby re-
sponses; by contrast, more complex tuning curves give rise to
a serpentine curve. The latter makes better use of the space of
possible population responses than the former and hence can
be expected to yield higher-resolution coding. Indeed, when
the population response is corrupted by noise of a given magni-
tude, it will elicit a smaller local error in the case of complex tun-
ing than in the case of simple tuning: by “stretching” the mean
response curve over a longer trajectory within the space of
possible population activities, complex tuning affords the code
higher resolution relative to the range of the encoded variable.
However, this argument does not capture in full the influence
of noise on the nature of coding errors. In the case of a winding
and twisting mean-response curve, two distant stimuli are some-
times mapped to nearby activity patterns. In the presence of
noise, this geometry gives rise to global (or catastrophic) errors.
The enhanced resolution of the neural code associated with the
occurrence of global errors was noted in the context of grid-cell
coding.³⁵,³⁹ Because of this trade-off, whether a simple or com-
plex coding scheme is preferable becomes a quantitative ques-
tion, which depends upon the details of the structure of the
encoding.
<figure><img src="image_5.png" /></figure>

(C) Examples of population activity (across the stimulus line, color indicates stimulus value) for three sample representation neurons (bottom), for increasing values of $\sigma$, and corresponding examples of typical tuning curves (top). When $\sigma \to 0$ (left, $\sigma = 0.001$), neurons produce uncorrelated random responses to different stimuli, generating a spiky curve made up by broken segments. As $\sigma$ grows ($\sigma = 0.015$, $\sigma = 0.03$) irregularities are smoothed out, and nearby stimuli evoke increasingly correlated responses. Ultimately, for large values of $\sigma$ (right, $\sigma = 0.15$), we recover a scenario similar to that with unimodal tuning curves.

<caption>Figure 1. Geometrical approach to coding, and the random feedforward neural network architecture</caption>

(A) Top: mean responses of three-neuron populations encoding a one-dimensional stimulus. Left: population of neurons with Gaussian tuning curves. Right: population of neurons with periodic tuning curves. Bottom: mean activity in the neural populations, parametrized by the stimulus value colored according to the legend, as one-dimensional curves in a three-dimensional space. Unimodal tuning curves (left) evoke a single-loop curve, which preserves the distances between stimuli in the evoked responses. Periodic tuning curves (right) evoke a more complex curve in which two distant stimuli may be mapped to nearby points in the joint-activity space; the curve is longer and fills up a larger portion of the activity space.

(B) Feedforward neural network as generative model of complex tuning curves. An array of $L$ sensory neurons with Gaussian tuning curves (one highlighted in purple) encodes a one-dimensional stimulus into an $L$-dimensional representation. These tuning curves determine the mean response of the population for a given stimulus, $x_0$ (dots). This layer projects onto a smaller layer of $N$ representation neurons with an all-to-all random-connectivity matrix, $W$, generating irregular responses. We plot the tuning curves of three sample neurons, highlighting their response to the stimulus $x_0$.

## A simple generative model of complex, irregular tuning curves

We consider irregular tuning curves that are instantiations of a stochastic, Gaussian process. This choice is unbiased in the sense that it is the maximum-entropy process given a two-point correlation function of the mean neural response across stimuli. We can thus think of this choice as a simple, unbiased model for the analysis of population coding in the presence of irregular tuning curves. The correlation function controls the degree of smoothness of tuning curves. It is a symmetric decaying function; for the sake of concreteness, we chose it to be Gaussian so that it is parametrized by only two numbers (variance and amplitude).

This kind of irregular tuning curve can be generated by a simple mechanism, namely, in a two-layer feedforward neural network with random synaptic weights. As we will see, the advantage of viewing the irregular tuning curves as being thus generated in the output layer of the network is that they can then be related to the properties of the classical, bell-shaped tuning curves in the input layer. Importantly, the irregular tuning curves do not rely on any finely tuned architecture or parameter tuning; they emerge from the variability in synaptic weights, and their statistics are controlled by the width of tuning curves in the input layer. Furthermore, recent work suggests that the detailed statistics of place fields in area CA1 are consistent with such an

architecture.40 We discuss the generality of our modeling choice in STAR Methods section “Generality and use of the network model” and further discuss possible biologically motivated extensions in STAR Methods section “Model Extensions.”

The architecture of the model network and the symbols associated with its various parts are illustrated in Figure 1B. In the first layer, a large population of $L$ sensory neurons encodes a one-dimensional stimulus, $x$, into a high-dimensional representation. Throughout, we assume that $x$ takes values between zero and one, without loss of generality. (If the input covered an arbitrary range, say $r$, then the coding error would be expressed in proportion to $r$. In other words, one cannot talk independently of the range of the input and of the resolution of the code. We set the range to unity in order to avoid any ambiguity.) Sensory neurons come with classical tuning curves; i.e., the mean activity of neuron $j$ in response to stimulus $x$ is given by a Gaussian with center $c_j$ (the preferred stimulus of that neurons) and width $\sigma$:

$$u_j(x) = A \exp \left( - \frac{(x - c_j)^2}{2\sigma^2} \right). \quad (\text{Equation 1})$$

Following a long line of models, we assume that the preferred stimuli in the population are evenly spaced, so that $c_j = j/L$. As a result, the response vector for a stimulus $x_0$, $u(x_0)$, can be represented as a Gaussian “bump” of activity centered at $x_0$.
Complex, irregular tuning curves appear in the second layer containing N representation neurons; we shall be interested in instances with $N \ll L$, in which efficient coding results in compression of the stimulus information from a high-dimensional to a low-dimensional representation. Each representation neuron receives random synapses from each of the sensory neurons; specifically, the elements of the all-to-all synaptic matrix, W, are i.i.d. Gaussian random weights with vanishing mean and variance equal to $1/L$ ($W_{ij} \sim N(0, 1/L)$). In the simple, linear case that we consider, the mean activity of neuron i in the second layer is thus given by

$$v_i(x) = \sum_{j=1}^{L} W_{ij} u_j(x). \qquad (\text{Equation 2})$$

In STAR Methods, we show that the tuning curves generated in this way, $v_i(x)$, are indeed samples from a Gaussian process with vanishing mean and Gaussian correlation function. The width of the correlation function is equal to $\sqrt{2\sigma}$, proportional to the width of the sensory tuning curves, $u_j$. Thus, $\sigma$ controls the smoothness of the tuning curves in the second layer. By the same token, for a fixed value of the prefactor A (Equation 1), the amplitude of the variations of $v_i$ with $x$ depends upon the value of $\sigma$. For a legitimate comparison of population codes in different networks, we set this amplitude to a constant on average,

$$\langle \int_{0}^{1} dx \left[ v_{i}(x) - \int_{0}^{1} dx' v_{i}(x') \right]^{2} \rangle_{W} = R, \qquad (\text{Equation 3})$$

by calibrating the value of the prefactor in Equation 1, A. Because of the averaging over the synaptic weights, indicated by the brackets $\langle \cdot \rangle_W$, A does not depend upon a specific realization of the synaptic weights. Equation 3 corresponds to the usual constraint of “resource limitation” in efficient coding models; it amounts to setting a maximum to the variance of the output over the stimulus space, as is commonly assumed in analyses of efficient coding in sensory systems.8,41–43

Returning to our geometric picture, we observe that, by changing the value of $\sigma$, we can interpolate between smooth and irregular tuning curves in the second layer (Figure 1C). In the limiting case of large $\sigma$, representation neurons come with smooth, unimodal tuning curves akin to classical ones; in the other limiting case of small $\sigma$, the mean population response curve becomes infinitely tangled. Thus, as the value of $\sigma$ is decreased, the mean-response curve “stretches out” and necessarily twists and turns in such a way as to fit within the allowed space of population responses defined by Equation 3. A longer population response curve fills the space of population responses more efficiently and represents the stimulus at a higher resolution, but its twists and turns may result in greater susceptibility to noise due to global errors.

To complete the definition of the model, we specify the nature of the noise in the neural response. We assume that the activity of neuron i in the second layer is affected by noise, which we denote by $z_i$, such that its response at each trial (in which stimulus x is presented) is given by $r_i = v_i(x) + z_i$. For the sake of

simplicity, we use Gaussian noise with vanishing mean and variance equal to $\eta^2$. In most of our analyses, we suppose that responses in the first layer are noiseless and that the noise in the second layer is uncorrelated among neurons. As a result, in the first layer, information is perfectly encoded, for any value of $\sigma$. Our focus is thus on the implications of noise in the second layer of neurons and how it affects the representation of a continuous stimulus perfectly encoded in the first layer of neurons. (In the STAR Methods section “Model extensions,” we discuss a generalization of the model in which first-layer neurons are noisy, inducing correlated noise among representation-layer neurons.) We emphasize that the coding properties of the representation neurons depend only on the statistics of their tuning curves, not on the specific mechanism with which these are generated. The advantage of a simple circuit mechanism is that it allows us to relate abstract parameters of the statistics of tuning curves to biophysical quantities; the nature of the circuit model, however, has no bearing on the results on information coding.

We quantify the performance of the code in the second layer through the mean squared error (MSE) in the stimulus estimate as obtained from an ideal decoder, “ideal” in the sense that it minimizes the MSE. (Throughout, in heuristic arguments and analytical calculations, we focus on the MSE. In a number figures, however, we plot its square root, the root-mean squared error (RMSE), so as to allow for a direct comparison with the stimulus range. The figure captions specify which of the two quantities is illustrated.) The use of an ideal decoder is an abstract device that allows us to focus on the uncertainty inherent to encoding (rather than to imperfections in decoding); it is nevertheless possible to obtain a close approximation to an ideal decoder in a simple neural network with biologically plausible operations (see STAR Methods).

## Compressed coding in the limiting case of tuning curves with extreme irregularity

It is instructive to study the properties of coding in our model in the limiting case of neurons with narrow tuning curves in the sensory layer ($\sigma \to 0$), because this limit yields the most irregular tuning curves in the representation layer of our network (Figure 1C). As we shall see, this limiting case also corresponds to that of a completely uncorrelated, random code, for which the mathematical analysis simplifies. When the value of $\sigma$ is much smaller than $1/L$, neurons in the sensory layers respond only if the stimulus coincides with the preferred stimulus of one of the neurons, and only that neuron is activated by the stimulus presentation; stimulus values that lie in between the preferred stimuli of successive sensory neurons in the first layer do not elicit any activity in the system. We can thus consider that any stimulus of interest is effectively chosen within a discrete set of L stimuli with values $x_j = j/L$, with $j = 1, \dots, L$.

Each of these stimuli elicits a mean response

$$v_i(x_j) = \tilde{A}W_{ij} \sim N(0, R) \qquad (\text{Equation 4})$$

in neuron i of the second layer. Here, the value of $\tilde{A}$ is chosen so as to set the amplitude of the variations of $v_i$ to be equal to the constant R (analogously to Equation 3 but for the case of discrete
<figure><img src="image_6.png" /><figcaption>Figure 2. Probability of error for narrow tuning curves in the sensory layer</figcaption></figure>

Simulation results are represented as mean ± SD.

(A) Joint mean responses of two neurons to $L = 50$ stimuli, colored according to the legend in Figure 1C. Noise is represented as a cloud of possible responses (in gray) around the mean. An error occurs when the noisy response, $r$, falls closer to a mean response corresponding to a stimulus, $\hat{x}$, different from the true one, $x_0$. Since mean responses are uncorrelated, $\hat{x}$ may be distant from $x_0$.

(B) Theoretical (solid curves, Equation 5) and numerical (dots) results for the probability of error as a function of the population size, for different values of $L$ ($\eta^2 = 0.5$). The probability of error scales exponentially with the number of neurons, $N$, with a multiplicative constant involving the number of stimuli, $L$.

(C) Theoretical (solid curves) and numerical (dots) results for the probability of error as a function of the population size for different values of $\eta^2$ ($L = 500$).

In stimuli). Geometrically, Equation 4 represents a mapping from $L$ stimulus values to a set of uncorrelated, random locations in the space of the population activity (as illustrated in Figure 2A for a two-neuron population). In any given trial, however, the responses in the representation layer are corrupted by noise (Figure 2A). The ideal decoder interprets a single-trial response as being elicited by the stimulus associated to the nearest possible mean response (Figure 2A). The outcome of this procedure can be twofold: either the correct or an incorrect stimulus is decoded; in the latter case, because the possible mean responses are arranged randomly in the space of population activity (Figure 2A; Equation 4), errors of all magnitudes are equiprobable. In other words, a model with narrow sensory tuning curves results in a second-layer code that does not preserve distances among inputs, and, consequently, the decoding error is either vanishing or, typically, on the order of the input range (set to unity here). The mean error is then simply proportional to the probability with which the ideal decoder makes a mistake, with a constant of proportionality of the order of the stimulus range.

In STAR Methods, we provide a derivation of this quantity. In the case of low-error coding, which interests us, we obtain the dependence of the probability of a decoding error as a function of the various model parameters, as

$$P_{\text{error}} \approx \frac{L}{\sqrt{2\pi N}} \exp\left(-\log\left(1 + \frac{R}{2\eta^2}\right)\frac{N}{2}\right). \quad \text{(Equation 5)}$$

The main dependence to note, here, is the exponentially strong suppression as a function of the number of neurons in the second layer (Figure 2B). By contrast, the probability of error scales merely linearly with the size of the stimulus space, $L$, as is expected in the low-error limit. This result implies that it is possible to compress information highly efficiently in a comparatively small representation layer ($N \ll L$) even though the code is completely random. The price to pay for the use of randomness is that any error is likely catastrophic (on the order of the stimulus range), but these large errors happen prohibitively rarely. It is also worth noting that

the rate of exponential suppression depends on the variance of the noise, $\eta^2$, or, more precisely, on the single-neuron signal-to-noise ratio, $R/\eta^2$ (where $R$ is the variance of the signal, Equation 3). In numerical simulations, we set $R = 1$ and we vary $\eta^2$ to explore different noise regimes. Interestingly, even when this signal-to-noise ratio becomes small, i.e., when the noise in the activity of individual neurons is comparable to modulations of their mean responses, the exponential suppression as a function of $N$ of the probability of error remains valid, with a rate approximately equal to $R/4\eta^2$.

We note that these results are in line with the central result of compressed sensing, stating that the number of random projections, $N$, needed to decode an $L$-dimensional signal that is $K$-sparse (meaning that it can be represented with at most $K$ non-zero components in some domain, with $K = 1$ in our case) with negligible error, scales logarithmically with $L$, $N > O(K \log(L/K))$.44,45

## Compressed coding with general irregular tuning curves: Trade-off between local and global errors

As we saw in the previous section, in the case of infinitely narrow tuning curves the coding of a stimulus in a given trial is either perfect or indeterminate; that is, any error is typically a global error, on the order of the entire stimulus range. In the more general case of irregular tuning curves generated from sensory neurons with arbitrary tuning width, the picture is more complicated: in addition to global errors, which result from the twisting and turning of the mean-response curve, the population code is also susceptible to local errors (Figure 3A). This is because broad tuning curves in the sensory layer partly preserve distances: locally, nearby stimuli are associated with nearby points on the mean-response curve. As a result, the coding of any given stimulus is susceptible to local errors due to the response noise. As the tuning width in the sensory layer, $\sigma$, decreases, two changes occur in the mean-response curve: it becomes longer (it stretches out) and it becomes more windy (Figure 1C). Stretching increases the local
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 3. Trade-off between local and global errors</figcaption></figure>

Simulation results are represented as mean ± SD. (A) Different types of error in an irregular curve of mean population activity (joint response of two neurons, colored according to the legend in Figure 1C). Here, $r'$ and $r''$ are two possible noisy responses to the same stimulus, extracted from the Gaussian cloud surrounding the mean response, $v(x_0)$. An ideal decoder outputs the stimulus corresponding to the closest point on the curve. In one case, $r'$ results in a local error by selecting a point on the curve that represents a nearby stimulus, $\hat{x}'$. In the other case, $r''$ is closer to a point on the curve that represents a stimulus distant from the true one, $\hat{x}''$, causing a global error.

(B) Normalized histogram of absolute error magnitudes, $\Delta x = |\hat{x} - x|$, made by an ideal decoder, for different values of $\sigma$ ($N = 25$). For better visualization, we consider a stimulus with periodic boundary conditions. The contribution of the two types of error varies with $\sigma$. For small $\sigma$, coding is precise locally (fast drop of the purple curve for small errors), but many global errors occur (tail of the distribution is high). For large $\sigma$ (green curves) local accuracy is poorer but global errors are suppressed.

(C) Theoretical prediction for the two contributions to the MSE as a function of $\sigma$ ($N = 30$). The magnitude of local errors increases with larger $\sigma$ (solid curve), while the number of global errors decreases (dashed curve).

(D) RMSE as a function of $\sigma$: comparison between numerical simulations (dots) and theoretical prediction of Equation 6 (solid curve).

(E) RMSE as a function of $\sigma$ for different population

sizes $N$ (increasing from violet to yellow). The smallest RMSE occurs at an optimal value of $\sigma$, $\sigma^*(N)$, which decreases with increasing $N$.

(F) Same data, but the error is displayed as a function of $N$, for a fixed value of $\sigma$. The MSE decreases exponentially rapidly until global errors are suppressed, then the local errors are linearly reduced. A smaller value of $\sigma$ implies a larger value of $N$ at which the crossover occurs, as well as a smaller MSE at this crossover value.

resolution of the code (because it allows for two nearby stimuli to be mapped to two more distant points in the space of population activity), while windiness increases the probability of global errors. This trade-off is apparent when we plot the histogram of error magnitudes as a function of $\sigma$: for larger values of $\sigma$, global errors are less frequent, but local errors are boosted (Figure 3B). Also noticeable, here, is that the large-error tails of the histograms are flat, consistent with the observation that global errors of all sizes are equiprobable. (Strictly speaking, this happens if the stimulus has periodic boundary conditions, such that, picking two random points, the probability that they are at a given distance does not depend on the location of one or the other point.)

For a more quantitative understanding, we carried out an approximate analytical calculation. First, we approximated the mean-response curve by a linear function locally. Second, we considered that the distance between two segments of the curve representing the mean responses to two stimuli distant by more than $\sigma$ is random and independent of the stimulus values. Using these two assumptions, we obtained the MSE as a sum of two terms (see STAR Methods for mathematical details) corresponding to local and global errors, as

$$ \epsilon^2 = \langle E^2 \rangle_W \approx \epsilon_l^2 + \epsilon_g^2 \approx \frac{2\sigma^2\eta^2}{RN} + \frac{\epsilon_g^2}{\sigma\sqrt{2\pi N}} \exp\left(-\log\left(1 + \frac{R}{2\eta^2}\right)\frac{N}{2}\right), \quad (\text{Equation 6}) $$

where $\epsilon_g^2$ is a term of $\mathcal{O}(1)$ that depends upon the choice of stimulus boundary conditions (see STAR Methods). This expression quantifies the MSE for a “typical” set of tuning curves, obtained by averaging over possible choices of synaptic weights, as indicated by the brackets $\langle \cdot \rangle_W$. The first term on the right-hand-side of Equation 6 represents the contribution of local errors, while the second term corresponds to global errors (Figure 3C). Their form can be intuited as follows. The magnitude of local errors is proportional to $\eta^2$ and inversely proportional to $N$, as in classical models of population coding with neurons with bell-shaped tuning curves (see, e.g., Zhang and Sejnowski$^{10}$). Furthermore, decreasing $\sigma$ stretches out the mean-response curve, which increases the local resolution of the code and explains the factor $\sigma^2$ in Equation 6. (The form of this first term can also be understood as the inverse of the Fisher information.$^{17,46}$) The second
<figure><img src="image_13.png" /><figcaption>Figure 4. Scaling of the optimal width and the optimal MSE as a function of population size and signal-to-noise ratio</figcaption></figure>

Simulation results are represented as mean ± SD. (A) The optimal $\sigma^*$ decreases exponentially rapidly with the number of representation neurons, saturating the lower bound imposed by the finite number of neurons of the first layer (corresponding to the spacing of the preferred positions, 1/L). Simulations (dots) show good agreement with analytical results (solid curve).

(B) The optimal RMSE is suppressed exponentially rapidly with N. Simulations (dots) agree with analytical results (solid curve).

(C and D) Optimal width (C) and RMSE (D) as a function of the parameters N and $\eta^2$. The color coding is in log scale in order to highlight the exponential scaling.

term on the right-hand-side of Equation 6 is obtained as an extension of Equation 5: instead of considering the probability that two mean-response points are placed nearby, we consider the probability that two segments of the mean-response curve with size $\sigma$ each fall nearby. There are $1/\sigma$ such segments (since we have set the stimulus range to unity), and this explains why the factor L in Equation 5 is replaced by a factor $1/\sigma$ in Equation 6. Importantly, the two terms in Equation 6 are modulated differently by the two parameters N and $\sigma$. Depending upon their values, either local or global errors dominate (Figure 3C).

We tested the validity of Equation 6: it agrees closely with results from numerical simulations (Figure 3D, see STAR Methods for details). The non-trivial dependence is illustrated by the observation that the MSE may decrease or increase as a function of $\sigma$, around a given value of $\sigma$, depending upon the value of N (Figure 3E). Furthermore, the strong (exponential) reduction in MSE with increasing N occurs only up to a crossover value that depends on $\sigma$ (Figure 3F); beyond this value, global errors disappear, and the error suppression is shallower (hyperbolic in N, due to improved local resolution). For small values of $\sigma$, the crossover values of N are larger and occur at lower values of the MSE.

As is apparent from Figures 3D and 3E, for any value of N there exists a specific value of $\sigma = \sigma^*(N)$ that balances the two contributions to the MSE such as to minimize it. This optimal width can be thought of as the one that stretches out the mean-response curve as much as possible to increase local accuracy but that stops short of inducing too many catastrophic errors. The MSE is asymmetric about the optimal width, $\sigma^*$: smaller values of $\sigma$ cause a rapid increase of the error due to an increased probability of global errors, while larger values of $\sigma$ mainly harm the code's local accuracy, resulting in a milder effect. From Equation 6, we obtain the dependence of the optimal width upon the population size, as

$$ \sigma^* \approx \left( \frac{\bar{\epsilon}_g^2}{4\eta^2} \sqrt{\frac{N}{2\pi}} \right)^{1/3} \exp \left( -\log \left( 1 + \frac{R}{2\eta^2} \right) \frac{N}{6} \right), \quad (\text{Equation 7}) $$

and the optimal MSE as a function of N, as

$$ \epsilon^{2*} = \langle E^2(\sigma^*) \rangle W \approx \left( \frac{\eta \bar{\epsilon}_g^2}{\sqrt{2\pi N}} \right)^{2/3} \exp \left( -\log \left( 1 + \frac{R}{2\eta^2} \right) \frac{N}{3} \right). \quad (\text{Equation 8}) $$

Both these analytical results agree closely with numerical simulations (Figures 4A and 4B). Equation 8 and Figure 4B show that the optimal MSE is suppressed exponentially with the number of representation neurons in the second layer. Thus, highly efficient compression of information and exponentially strong coding also occurs when tuning curves of representation neurons are smooth: furthermore, a degree of smoothness in the tuning curves is advantageous. With the optimal choice of the sensory tuning width, the rate of scaling with N of the argument within the exponential in Equation 8 depends upon the noise variance, $\eta^2$; in Figures 4C and 4D, we illustrate the dependence of $\sigma^*$ and $\epsilon^*$ upon N and $\eta^2$.

## Compressed coding of multi-dimensional stimuli

We generalized our framework to include higher, K-dimensional stimuli. Complex, irregular multi-dimensional tuning curves can be defined in multiple ways, depending on how neural responses are jointly affected by different stimulus dimensions. In one limiting case, neurons respond independently to different stimulus dimensions, with multi-dimensional tuning curves obtained as simple superpositions of one-dimensional ones. In the other
A

<figure><img src="image_14.png" /></figure>

B

<figure><img src="image_15.png" /></figure>

D

<figure><img src="image_16.png" /></figure>

C

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

E

<figure><img src="image_19.png" /></figure>

lower in the pure case. As soon as *N* is sufficiently large and *σ* is sufficiently large to allow for coverage of the stimulus space, the conjunctive case outperforms the pure case. This effect is stronger in the small-*σ* region due to the slower scaling of the global errors in the pure case. When *σ* is large, the ratio saturates at the value given by the ratio of the local errors.

(D and E) Optimal tuning width (D) and relative RMSE (E), for pure (blue, red) and conjunctive (green, violet) cases. The global error decreases more slowly in the pure case. For *N* ≥ 75 the optimal width in the conjunctive case saturates, due to loss of stimulus coverage, while the pure population does not suffer from this limitation. Thus, the RMSE in the conjunctive case stops decreasing exponentially and starts decreasing only linearly with *N*.

limiting case, tuning curves are not additively separable func-
tions of stimulus dimensions. We model this case by describing
tuning curves as samples from a stochastic processes with a
correlation function that jointly depends on all stimulus dimen-
sions. In our model network, one limiting case is achieved
when each sensory neuron is tuned to a single stimulus dimen-
sion and insensitive to others. As a result, tuning curves of repre-
sentation neurons are obtained as sums of independent realiza-
tions of Gaussian processes for each stimulus dimension. In the
other limiting case, when sensory neurons are each sensitive to
all stimulus dimensions, tuning curves in the second layer are
samples from multi-dimensional Gaussian processes. We refer
to these two coding schemes in the sensory layer as “pure”
and “conjunctive,” following Finkelstein et al⁴⁷, who examined
them in the context of head-direction neurons in bats. We

emphasize that neurons in the representation layer receive
random inputs from all sensory neurons; as such, representation
neurons are sensitive to all stimulus dimensions no matter what
coding scheme is implemented in sensory neurons. For the sake
of simplicity, we consider here only the two limiting cases of pure
and conjunctive coding in the sensory layer of our model.

By extending the geometric picture (illustrated in Figure 1
for the case of a one-dimensional stimulus), we can analyze dif-
ferences in coding properties between pure and conjunctive
coding schemes; in Figure 5A, we illustrate the case of a two-
dimensional stimulus. In this case, the mean response of repre-
sentation neurons corresponds to a mapping from a two-dimen-
sional stimulus space to a random “sheet” (a two-dimensional
surface) in the *N*-dimensional space of the population activity.
In the pure case, the activity of a given sensory neuron is

<caption>Figure 5. Compressed coding with multi-dimensional stimuli</caption>

Simulation results are represented as mean ± SD.
We illustrate the case of a three-dimensional
stimulus, with L = 3375, η² = 1, and R = 1.

(A) Mapping of a multi-dimensional stimulus space into neural population activity as obtained from a two-layer coding scheme. Top: two-dimensional stimulus space; colors serve as stimulus legend for subsequent plots. Middle: mean activity (z-axis) of a sample sensory neuron, for two cases, as a function of the two stimulus coordinates (x- and y-axis). In the pure case (left), a single sensory neuron folds the two-dimensional sheet across a direction, specified by its preferred position and dimension (here, x₂). In the conjunctive case (right), a sensory neuron creates a bump in the sheet. Bottom: joint activity of three representation neurons as a function of the stimulus. Each of these neurons randomly sums the output of sensory neurons, producing a randomly folded sheet in the pure case (left) and a crumpled sheet in the conjunctive case (right).

(B) RMSE as a function of σ for different population
sizes *N* (increasing from violet to yellow), when the
first layer consists of pure (left) or conjunctive
(right) cells. The optimal σ, which decreases with
*N*, optimizes the balance between local and global
errors, similarly to the one-dimensional case. In
the conjunctive case, the rapid increase of the
RMSE below σ = 0.05 is due to the sensory neu-
rons not tiling the stimulus space, and it is inde-
pendent of *N*.

(C) Ratio of the RMSE in the two cases, ε<sub>c</sub>/ε<sub>p</sub>, as a function of σ and N. The yellow (violet) region indicates an outperformance of the pure (conjunctive) population. To aid visualization, the yellow region indicates all the values greater than 2. This regime of small σ is characterized by a better coverage of the pure population, independently of N. Values greater than one occur also when N is small, due to the prefactor of the global error being
maximally modulated when the stimulus varies along a particular dimension, the one to which the neuron is sensitive. Variations of the stimulus along orthogonal directions have no effect on the mean neural activity. Neurons in the representation layer compute a linear sum of these responses, and therefore their activity can be decomposed as a sum of one-dimensional functions. This implies that the response sheet is maximally curved along each of the stimulus dimensions; geometrically, this results in a folded structure, with creases along the directions of mild sensitivity. By contrast, in the conjunctive case, the activity of a sensory neuron is modulated by variations of the stimulus along any direction. As a result, the response sheet that represents the joint mean activity of neurons in the second layer comes with (random) curvature equally along all stimulus dimensions: rather than folded, it behaves like a “crumpled” sheet (Figure 5A).

This geometric picture offers an intuitive explanation of the behavior of the MSE in the two coding schemes. (For the corresponding mathematical treatment, see STAR Methods.) The local error is determined by how much the response sheet is stretched out; in turn, the more the response sheet is stretched out, the more it has to fold (or crumple) to fit in the allowed range of neural activity. Folding allows for a more modest stretching of the sheet than crumpling, and, as a result, the pure scheme incurs a larger local error than the conjunctive scheme (see Equations 71 and 76). The behavior of the global error is also different in the two coding schemes; there are two mechanisms at play here. First, in the pure scheme, as tuning curves are superpositions of one-dimensional functions, global errors occur primarily in a single stimulus dimension. This is also apparent in Figure 5A: the folded structure of the response sheet induces global errors in a single stimulus dimension. By contrast, in the conjunctive scheme global errors occur in an arbitrary number of stimulus dimensions. Second, the total variance of the tuning curve across the stimulus space is fixed (and, in particular, set to the same value for the pure and conjunctive schemes), but the signal-to-noise ratio, which governs the rate of error suppression with N, scales differently as a function of K. Both mechanisms, in a regime in which N is large enough to suppress the contribution of global errors, enhance the probability of global error in the pure scheme as compared to the conjunctive scheme (compare Equations 81 and 83 in STAR Methods). Intuitively, this is because a folded sheet has a larger surface area of contact with itself than a crumpled sheet. Thus, for sufficiently large values of N, the conjunctive scheme is more favorable than the pure one. The corresponding crossover value of N, however, depends on K, and large values of K impose a stringent constraint in the conjunctive case.

We illustrate these conclusions with numerical results in the case of a three-dimensional stimulus (K = 3), relevant to the data analysis we present in the next section. In Figure 5B, we illustrate the behavior of the RMSE as a function of N and σ for the pure and conjunctive coding schemes. In order to quantify the relative advantage of one scheme with respect to the other, we plot the ratio of the RMSE in the two schemes as a function of N and σ (Figure 5C). The resulting, relatively intricate pattern can be understood by considering different regimes. If

the population size is small, the pure scheme slightly outperforms the conjunctive one (because of a difference in the prefactors that affect the probability of error in the two cases); in this regime, global errors dominate and coding is poor overall. At larger values of N, the contribution of local errors becomes non-negligible. If local errors dominate relative to global errors (which occurs for large N and sufficiently large σ), then the conjunctive scheme outperforms the pure one, and the ratio of the RMSEs approaches the ratio between local errors only (Equation 77 for K = 3, implies $ε_{1,c}/ε_{1,p} \approx 1/\sqrt{3}$). In the non-trivial regime in which local and global errors are balanced (for large N and intermediate values of σ), the advantage of the conjunctive scheme is further boosted. As explained above, this is due to a stronger suppression of global errors as a function of N in the conjunctive case. Finally, if σ becomes smaller than a crossover value that depends on the number of sensory neurons, the latter no longer cover the stimulus space sufficiently densely, and the conjunctive scheme breaks down; in this regime, thus, the pure scheme is favored.

As illustrated in Figure 5B, similar to the one-dimensional case, there exists in each of the two coding schemes an optimal value of the tuning-curve width, σ, which achieves a balance between local and global errors, and it decreases with N. This dependence is somewhat different in the two coding schemes (Figure 5D) and contributes to the form of the suppression of the RMSE in the two schemes (Figure 5E). Both quantities, the optimal tuning-curve width and the RMSE, decrease more rapidly as a function of N in the conjunctive scheme. This results from the fact that global errors are suppressed more strongly with N in the conjunctive case (as explained above), and therefore a smaller σ, yielding a lower local error, is preferable. At the same time, the requirement that sensory neurons cover the stimulus space yields a more stringent constraint on σ in the conjunctive scheme, yielding a bound on the extent of the regime of exponential error suppression.

## Compressed coding in monkey motor cortex

The activity of neurons in the primary motor cortex (M1) of monkey is correlated with the location and movement of the limbs. Here, we consider spatial tuning in the context of a “static task.”⁴⁸ In this task, the monkey is trained to keep its hand motionless during a given delay after having placed it at one of a set of preselected positions on a three-dimensional grid labeled by the vector $\mathbf{x} = (x_1,x_2,x_3)$. Tuning curves of hand-position selectivity can be extracted from recordings in M1,⁴⁸,⁴⁹ and it has been customary to model these as a linear projection of the hand position onto a so-called preferred vector or positional gradient, $\mathbf{p}$, which thus points in the direction of maximal sensitivity.⁴⁹ The tuning curve of neuron i is then written as

$$v_i(\mathbf{x}) = a_i + p_{1,i}x_1 + p_{2,i}x_2 + p_{3,i}x_3 = a_i + \mathbf{p}_i \cdot \mathbf{x}. \quad (\text{Equation 9})$$

A recent study³³ noted, however, that a model of tuning curves that includes a form of irregularity yields an appreciably superior fit to the simple linear behavior of Equation 9. This more elaborate model bears similarity with our model of irregular tuning curves, and this naturally led us to ask about potential coding advantages that a complex coding scheme may afford M1.
<figure><img src="image_20.png" /><figcaption>Figure 6. Irregular vs. linear tuning</figcaption></figure>

Simulation results are represented as mean ± SD. (A) Mean fractional improvement for irregular tuning, as compared to linear tuning, as a function of population size and tuning-curve width, for a fixed noise variance, $\eta^2 = 1$. The black line indicates the critical values of *N* and $\sigma$ at which the two coding schemes perform equally well. In the region below (violet), global errors penalize the irregular case, making a smoother code more efficient. With increasing *N*, global errors become rarer while irregularities improve the local accuracy of the code (yellow region). This advantage increases at smaller values of $\sigma$, but so does the value of *N* required for the irregular case to be advantageous.

(B) Mean fractional improvement in the irregular case, generated with the data-fitted model, compared to the linear one, as a function of *N* and $\eta^2$. At small population sizes, irregular tuning curves produce global errors, and smoother tuning curves perform better (violet region, $\Delta\epsilon < 0$). By increasing *N*, global errors are suppressed and irregularities improve the local accuracy (yellow region, $\Delta\epsilon > 0$). The black line marks the transition values.

(C and D) Mean fractional improvement (C) and RMSE (D) in the irregular case as a function of population size for the noise model extracted from data. A noise variance is assigned to each neuron according to the distribution extracted from the data, shown in the inset of (C). For small *N*, linear tuning yields a better coding performance. At *N* ~ 40, the higher local accuracy compensates for global errors, and the irregular code starts to perform better, although the error is still substantial. The improvement saturates to a finite value of ~ 0.4 at a value of *N* ~ 100, when global errors are fully suppressed; the scaling of the error as a function of the population size is no longer exponential but only hyperbolic.

To be more specific, one can interpret the first layer in our network, made up of neurons with three-dimensional Gaussian tuning curves, as representing neurons in the parietal reach area (or premotor area), which are known to display spatially localized tuning properties.⁵⁰ This population of neurons projects onto a smaller population of M1 neurons that display spatially extended and irregular tuning profiles. As we have shown in the previous sections, whether a coding scheme with complex, irregular tuning curves is preferable to one with simple tuning curves depends on the degree of irregularity, on the population size, and on the magnitude of the noise. We quantify the merit of a compressed code making use of irregular tuning curves by computing the MSE, $\epsilon_{irr}^2$, and comparing the latter with the corresponding quantity in a coding scheme making use of smooth tuning curves defined in Equation 9, $\epsilon_{lin}^2$. We plot our results in terms of the mean fractional improvement, $\Delta\epsilon \equiv (\epsilon_{lin} - \epsilon_{irr})/\epsilon_{lin}$. $\Delta\epsilon$ is positive when irregularity favors coding, and is at most equal to unity (in the extreme case in which irregularity allows for error-free coding). First, we explore the performance of the two coding schemes for different values of the parameters *N* and $\sigma$, in the case in which all neurons have the same noise variance (Figure 6A). We note the existence of a crossover value of *N*, $N^*$, such that $\Delta\epsilon(N^*) = 0$. When $N < N^*$, small values of $\sigma$ induce prohibitively frequent global errors in the compressed (irregular) coding scheme, and linear (smooth) tuning curves are more efficient. For $N > N^*$, however, irregularities are always advantageous, and the more so the smaller the value of $\sigma$.

Next, we set the degree of irregularity in tuning curves by fitting our model to recordings from M1 neurons by Lalazar at al.³³ (see STAR Methods for details). Following these authors, we did not

aim at fitting individual tuning curves; instead, we allowed for randomly distributed synaptic weights (as in our original model) and we fitted a single parameter, the width of the tuning curves in the first layer, $\sigma$. The fit was aimed at reproducing specific summary statistic of the data referred to as “complexity measure” (a discrete version of the Lipschitz derivative that quantifies the degree of smoothness of a curve; see STAR Methods and Lalazar at al.³³). The complexity measure varies from neuron to neuron, and we chose $\sigma$ so as to minimize the Kolmogorov-Smirnov distance (see Equation 120 in STAR Methods) between the distribution implied by our model and the one extracted from the data. While our model is somewhat simpler than a model of irregular M1 tuning curves employed previously,³³ it yields a comparable fit. Given the best-fit value of $\sigma$ ($\approx 23$), the crossover population size, $N^*$, takes a moderate value that depends on the magnitude of the noise; the larger the noise, the larger $N^*$ (Figure 6B).

Finally, for a more realistic modeling of M1 neurons, we analyzed the performance of a model in which each neuron's noise variance is extracted from data (Figures 6C and 6D). For each recorded neuron, we computed the variance of the signal as the variance, across different stimuli, of the mean firing rate (left-hand side of Equation 3). Then, we estimated the variance of the noise by averaging the trial-to-trial variability of responses to the same stimulus. These two quantities allowed us to define a signal-to-noise ratio for each neuron of the population (see Equation 121 in STAR Methods). As in simulations we set the variance of the signal for each neuron to a constant value, we modeled the heterogeneity in the signal-to-noise ratio as a heterogeneous noise variance; the resulting distribution is skewed,
with an appreciable fraction of neurons exhibiting low signal-to-noise ratios (Figure 6C, inset). For each value of *N*, we sampled eight different pools of *N* neurons from the population, and we averaged the corresponding mean fractional improvement, $\Delta\epsilon$. We found, again, that the relative merit of compressed coding (with irregular tuning curves) grows with the population size; interestingly, when compressed coding becomes advantageous ($\Delta\epsilon > 0$ in Figure 6C), the error magnitude is still appreciable (Figure 6D). This means that, even though local and global errors are balanced, both contributions are substantial. $\Delta\epsilon$ continues to grow with *N* until global errors are suppressed; beyond this second crossover value, $N_{local}$, $\Delta\epsilon$ saturates because, in both coding schemes (with irregular and linear tuning curves), local errors dominate. Correspondingly, the MSE scales differently for *N* above or below $N_{local}$. When $N < N_{local}$, the MSE decreases exponentially with *N* due to the suppression of global errors, while, when $N > N_{local}$, the suppression of the MSE is hyperbolic in *N*, reflecting the behavior of local errors only (Figure 6D). This second crossover occurs at $N_{local} \approx 100$, a figure comparable to the number of neurons that control individual muscles in this specific task, as estimated from decoding electromyography (EMG) signals corresponding to individual muscles from subsets of M1 neurons.³³

## Dimensionality of a compressed neural code

The set of mean responses of a population of *N* neurons to a continuous stimulus can be viewed as a geometric object in an *N*-dimensional space, often referred to as neural manifold.⁵¹,⁵² In the STAR Methods section “Dimensionality of a compressed neural code” and Figure S2, we relate the picture put forth by our model to recent work that quantifies the dimensionality of neural manifolds inferred from neural recordings.⁵³,⁵⁴ In our model, optimal coding of a one-dimensional stimulus obtains when the neural manifold has a dimensionality of approximately *N*/2, as defined by a dimensionality reduction method⁵⁵ (see STAR Methods). This observation, that a one-dimensional stimulus can give rise to an effectively high-dimensional neural manifold, highlights the subtlety of the relation and the challenge of characterizing a manifold’s dimensionality from data when the neural code is non-linear.

## Model extensions

We proposed a simple, two-layer, feedforward neural network as a circuit mechanism that can generate downstream irregular tuning curves from upstream simple, bell-shaped ones. Beyond this simple scheme, we have investigated how biologically realistic features—noise correlations, synaptic sparsity, Dale’s law, non-linear responses, and recurrent connectivity—affect the coding performance (see STAR Methods section “Model extensions” and Figures S3–S6). We find that noise correlations inherited from sensory noise and produced through divergent synaptic connections modestly degrade coding accuracy, while random noise correlation structures generically enhance coding accuracy (see also section “Limitations of the study”). Imposing synaptic sparsity or Dale’s law yields a quantitative suppression of the coding accuracy in the small tuning-width regime but does not alter the qualitative behavior of the model. A rectifying non-linearity applied to the response function of

representation neurons mildly reduces the effective size of the neural population and thereby slightly increases the optimal tuning width. In particular, a variant of our simple, linear model with a thresholding non-linear transfer function (that sets the neural response to zero when the summed input is lower than a given value) has been used to fit recordings from the CA1 area of hippocampus; the quality of the fit indicates that our model, with minor refinements, captures the behavior of actual neurons⁴⁰ (see also section “Limitations of the study”). Finally, we show that recurrent connectivity in the representation layer affects noise correlations but preserves the coding performance in steady-state conditions.

## DISCUSSION

### Summary

We analyzed the properties of a neural population encoding a one-dimensional, continuous stimulus by means of random, spatially extended and irregular tuning curves, which might emerge from simple ones in a neural network with random synaptic weights. This model can interpolate between an extremely irregular coding scheme, highly efficient but prone to catastrophic errors, and a smoother one, more robust in the face of noise. Optimality is achieved at an intermediate level of irregularity, which depends on the population size and on the variance of the noise. At optimality, the mean error is suppressed exponentially with population size; as a result, irregular neural codes allow to compress the representation of a low-dimensional, continuous stimulus from a large first layer of neurons to a small second layer. We extended these results to the case of multi-dimensional, continuous stimuli, more intricate because each sensory neuron can be sensitive to one, some, or all stimulus dimensions. We considered two extreme cases: a pure coding scheme, in which multi-dimensional tuning curves are superpositions of one-dimensional ones, and a conjunctive coding scheme, in which neurons are sensitive to the combination of all stimulus dimensions. These two coding schemes correspond, in our model network, to different arrangements of tuning curves in the first layer. We examined the relative advantage of one scheme with respect to the other⁴⁷,⁵⁶ and elucidated its dependence on the number of representation neurons and on the tuning parameters. These analyses enabled us to revisit data from M1 neurons in monkey³³ and to discuss the benefits of an irregular code in the context of the representation of hand position. Finally, we broadened the picture of compressed coding by considering various model extensions and by relating our picture to analysis of the dimensionality of population activity.

### “Exponentially strong” neural population codes

Our results on the exponential scaling of the mean error with population size are similar to results obtained in the context of the representation of position by grid cells.³⁴⁻³⁶,³⁸,⁵⁷ According to the terminology adopted in this literature, the random compressed coding presented here is an exponentially strong population code. Grid-cell-tuning is a particular instance of exponentially strong codes making use of periodicity; the model presented here offers another example, in which tuning curves are random.
The notion of an exponentially strong code predates work in computational neuroscience: Shannon introduced it in the context of communication systems and analog signals.⁵⁸ In his framework, a sender maps a “message” (a continuously varying quantity analogous to our stimulus) into a “signal” (a higher-dimensional continuous quantity analogous to the output of our representation layer), which is transmitted over a noisy channel and then decoded by a receiver. The specific illustration he provides is that of a one-dimensional message mapped into a higher-dimensional signal (Figure 4 in Shannon⁵⁸), analogous to the mapping illustrated in Figure 1C; this mapping corresponds to a curve that wraps around in a higher-dimensional space. Shannon argues that an efficient code is obtained by stretching this curve to make it as long as possible up to the point at which the winding and twisting causes the curve to pass too close to itself, thereby generating catastrophic errors.

However, Shannon went further and showed that such a code need not to be carefully designed. His calculation corresponds, in our framework, to the case of infinitely narrow tuning curves in the sensory layer (Figure 2): he demonstrated that it is possible to send a discrete set of messages with an error suppressed exponentially in the dimensionality of the signal. Our work proposes an extension of this “fully random” scenario for the representation of a continuous variable based on a smooth, but irregular, mapping in a higher-dimensional encoding space. By varying the width of tuning curves in sensory neurons, $σ$, one can modulate the smoothness of the mapping and trade off global errors with local errors. In this more general, “correlated random” scenario, it is optimal to choose a non-vanishing value of $σ$ that depends on the population size and other model parameters.

## Coding with complex, irregular tuning curves

A large body of literature has addressed the problem of coding low-dimensional stimuli in populations of neurons with simple tuning curves. The most common assumption is that of bell-shaped tuning curves; these are often chosen to model sensory coding in peripheral neurons. Various studies set in this context discussed the shape of optimal tuning curves as a function of population size and stimulus dimensionality,¹⁰ stimulus geometry,⁵⁹ and the timescale on which coding operates.¹¹,¹⁴ More recent work analyzed the influence of a (non-uniform) prior distribution of stimuli on the optimal arrangement and shapes of tuning curves across a population of neurons; a particular prediction is that the tuning-curve width is narrower for neurons with a preferred stimulus over-represented in the prior.¹⁵,⁶⁰,⁶¹ Beyond simple, bell-shaped (or monotonic) tuning curves, and, especially, motivated by the discovery of grid cells, a number of studies have discussed the implications of periodic tuning curves for coding and their optimal allocation in modules with different periodicities.³⁴⁻³⁶,⁵⁷,⁶² In particular, the benefit of periodic tuning curves, as compared to spatially localized ones, on the coding accuracy was analyzed.

Our study falls in this line of work: we consider tuning curves that are more complex and span all the stimulus space, as opposed to bell-shaped, localized ones, and irregular curves (in particular, aperiodic, as opposed to two-dimensional grid cells). In recent years, there has been an increasing number of

experimentally observed tuning curves that exhibit these characteristics.²⁴,²⁸,³⁰,³³ Irregular tuning curves have also been used to fit data from neural recordings. Recent studies⁶³⁻⁶⁵ have aimed at fitting neural responses by defining a set of tuning curves as functions of a low-dimensional, latent representation of the stimulus. Crucially, in contrast to tuning curves commonly defined as functions of observable stimulus features, and fitted with simple parametric forms, here tuning curves (which are functions of a latent quantity to be identified) are multimodal and irregular.

Given the ubiquity of these observations, we chose a simple model to investigate the implications of complex tuning curves on coding properties in general terms: namely, a family of random tuning curves sampled from a Gaussian process in the stimulus space. A multivariate Gaussian distribution of neural responses across stimuli is indeed the maximum-entropy distribution given the mean neural responses and their covariance matrix, and, up to a non-linear transformation, it has been used as a prior in fitting neural data recordings.⁶³ Thus, our approach, while simple enough to enable us to derive explicit results, in fact rests on a random process that is used to fit actual neural activity.

Our work also points to the benefits of heterogeneity for coding. This general observation has been made before, but the emphasis was on the heterogeneity among neurons.¹⁶,⁶⁶⁻⁶⁸ In our case, it is each neuron's tuning curve that displays heterogeneity. The irregularity of tuning curves makes the response of any single neuron highly ambiguous; the resulting code is thus distributed, and the neural population as a whole is viewed as the relevant unit of computation.⁶⁹

Distributed codes have been argued to come with high capacity. An early example was developed in the context of face coding in the superior temporal sulcus of monkey.⁷⁰ Data analysis indicated that single-neuron sensitivity was heterogeneous and uninformative, but the number of distinguishable face stimuli grew exponentially with the population size. Our work provides an example of a random distributed code for continuous stimuli that exhibits similar scaling properties. The main difference is that, in the case of continuous stimuli, the precise identity of the stimulus cannot be recovered in presence of noise, and what matters is the magnitude of the distance between the decoded stimulus and the true one, quantified by an appropriate metric. In other words, both the probability of occurrence of an error and its magnitude matter. The requirement of minimizing the MSE then yields a particular coding scheme that balances small (local) and large (global) errors.

## Geometry and dimensionality of population responses

In the past decade, the progress in experimental methods has allowed for the recording of neural populations on a large scale.⁵⁵,⁶⁹ In an effort to interpret the way in which information is represented in population activity, various approaches have been focusing on the geometric properties of population responses to a battery of stimuli.⁵²⁻⁵⁴,⁷¹ Points in a high-dimensional space, each corresponding to the neural population response to a stimulus, are often interpreted as being located on a manifold that describes the space of possible population activity. Quantifying the geometry and, more specifically, the dimensionality of this manifold offers a characterization of neural
population activity. This geometric element is eminently relevant in our work. In particular, our theoretical framework illustrates the connection between single-neuron tuning properties and a non-trivial structure of the population response, bringing together two perspectives that are often considered separately.⁷² This, in turn, allows us to study how the fine structure of the neural response manifold determines the coding properties of the population, a feature that can be easily missed when the geometry of population activity is characterized through dimensional reduction methods applied to a dataset.⁷³ Indeed, even with non-linear methods, the observation of a limited number of neurons and the noise in neural responses might hinder the identification of a low-dimensional structure of the neural manifold.

A specific geometrical question is that of the dimensionality of the population response in the representation layer. We showed that the spectrum of the covariance of the population activity in the representation layer, across the stimulus space, comes with a band-pass structure; by decreasing the width of tuning curves in the sensory layer, the band-pass profile acquires additional modes. Stringer et al.⁵⁴ discussed a similar picture in analyzing recordings from a large population of visual neurons responding to a large, but discrete, set of images. In their case, the spectrum of the covariance matrix of population responses exhibits an algebraic (power-law) tail, and the authors argue that this property allows for a high-dimensional population activity while retaining smoothness of the code. Our work presents a different, and more elementary, mechanism by which a large number of modes can be accommodated by the population activity (while retaining smoothness). The non-trivial point, in our case, is that it is not beneficial for coding to be poised in the limiting case in which the number of modes is maximal but the code becomes singular (non-smooth), as, in this limit, global errors proliferate. The optimal effective dimensionality of the response manifold, as defined by the participation ratio, lies at an intermediate value at which intersections of the manifold with itself are rare and local and global errors are balanced (Figures S2A–S2C).

## Compressed sensing

We studied a network in which the information encoded in a high-dimensional activity pattern is compressed into the activity of a comparatively small number of neurons, a setting that exhibits analogies with the one of compressed sensing.⁴⁴ Compressed sensing is a signal-processing approach for reconstructing *L*-dimensional signals, which are *K*-sparse in some basis (i.e., they can be expressed as vectors with only *K* non-vanishing elements), from *N* linear, noisy measurements, with *K* ≪ *L* and *N* ≪ *L*.⁷⁴ In our study, the low dimensionality of the stimulus, *x*, implies sparsity of the *L*-dimensional activity of the sensory layer as long as the tuning curves in the sensory layer are not too wide.

A central result in the field of compressed sensing is that random measurements can yield near-optimal reconstructions. Furthermore, for near optimality to be achieved, the required number of measurements scales approximately linearly in *K* and only logarithmically in the dimensionality of the signal: *N* > O(*K* log(*L*/*K*)).⁴⁴,⁴⁵ In effect, in our network, the representation layer extracts a limited number of random measurements from

the sensory layer. We obtain an analog scaling form by inverting Equation 5: the number of random projections, *N*, necessary to decode *L* different stimuli with negligible error scales logarithmically with the number of stimuli. We note, however, that our framework differs from that of compressed sensing as the objective is to decode the identity of the stimulus rather than a high-dimensional signal vector (in our case, the activity pattern of the sensory layer). Thus, our optimal decoding mechanism differs from the one typically employed in these settings.

## Encoding vs. decoding

We focused in this study exclusively on the properties of encoding in a neural population. For this aim, we assumed an ideal decoder; in principle, this is not a limitation: we show in STAR Methods that an ideal decoder can be implemented by a simple, two-layer neural network. The first layer computes a discretized approximation of the posterior distribution over stimuli, and the second layer computes the mean of this distribution in such a way as to minimize the MSE (similar decoders have been employed by others¹⁵,⁷⁵). Furthermore, all the operations carried out by this two-layer network—linear filtering, non-linear response, and normalization—are plausible biological operations.¹³,⁷⁶,⁷⁷ The parameters involved, however, have to be chosen with the knowledge of the tuning curves and noise model.

One can ask whether biologically plausible learning rules can result in a decoder that approximates the ideal one. A closely related question has been examined by Bordelon et al.,⁷⁸ who analyzed how the generalization error in a deep neural network trained with gradient descent depends on the number of training samples and on the structure of the decomposition of a target function into a set of modes (e.g., Fourier modes). Bordelon and Pehlevan⁷⁹ found that learning the high-frequency Fourier components of a target function requires a larger number of training samples, as compared to learning its low-frequency components. Similarly, in the context of our network, one expects that learning a decoder in the case of narrow tuning curves in the sensory layer is more laborious than in the case of broad tuning curves. Noise in the training samples may also hamper learning severely in the presence of global errors. Furthermore, one can ask how our results might be modified if decoding is carried out by a decoder different from the ideal one, for example by a decoder obtained through adequately chosen learning rules. We leave these questions for future work.

## Limitations of the study

This study assumes tuning curves sampled from a Gaussian process and uncorrelated noise. The simplicity of the framework makes it tractable analytically, but it ignores possibly important aspects of neural systems, such as correlated noise and recurrent connectivity, which can potentially affect information coding.⁸⁰,⁸¹ Furthermore, while irregular tuning curves offer theoretical advantages, their prevalence in the brain and functional role require further exploration. Recent evidence from hippocampal CA1 neurons demonstrates that a version of the proposed model with rectifying response functions captures neural data with quantitative accuracy.²⁸,⁴⁰ The implication is that hippocampal area CA1 compresses inputs from area CA3 into an efficient, sparse representation. Finally, while this work focuses on
encoding, it raises the issue of decoding: complex encoding requires a flexible decoding machinery. Whether optimal decoding is realistic for the kind of encoding we propose here, and whether, alternatively, sub-optimal but effective decoders can be learned from experience are questions left for future research. A complete theory would consider encoding and decoding jointly and explain encoding schemes from functional demands and restrictions imposed by possibly limited decoding resources.⁸² (See STAR Methods section “Limitations of the study (extended)” for an extended discussion of the limitations of the study.)

## RESOURCE AVAILABILITY

### Lead contact

Requests for further information and resources should be directed to and will be fulfilled by the lead contact, Rava Azeredo da Silveira (rava@iob.ch).

### Materials availability

This study did not generate new materials.

### Data and code availability

* The monkey motor cortex data were reported and discussed by Lalazar et al.³³ Data are publicly available at OSF: https://osf.io/u57df/.

* All original code has been deposited at Zenodo: https://doi.org/10.5281/zenodo.14210005 and is publicly available as of the date of publication at GitHub: https://github.com/simonebmalerba/RandomCompressedCoding.git.

* Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

## ACKNOWLEDGMENTS

We thank Anatoly Khina and Trang-Anh Nghiem for helpful discussions. Support for this research was provided by the CNRS through Unité Mixte de Recherche (UMR) 8023, the Swiss National Science Foundation Sinergia Project (CRSII5_173728), the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program grant agreement no. 866376, the European Union's Horizon 2023 research and innovation program under the Marie Skłodowska-Curie Actions grant 101152984 ASTRONET, the Simons Foundation Autism Research Initiative (SFARI) grant No. 982347, the Israel Science Foundation (grant nos. 1319/13, 1978/13, and 1745/18), the Gatsby Charitable Foundation, and the William N. Skirball Chair in Neurophysics.

## AUTHOR CONTRIBUTIONS

Conceptualization, S.B.M., Y.B., and R.A.d.S.; methodology, S.B.M., Y.B., and R.A.d.S.; software, S.B.M.; investigation, S.B.M. and M.P.; formal analysis, S.B.M., Y.B., and R.A.d.S.; writing – original draft, S.B.M.; writing – review & editing, S.B.M., M.P., Y.B., and R.A.d.S. funding acquisition, Y.B. and R.A.d.S.; supervision, R.A.d.S.

## DECLARATION OF INTERESTS

The authors declare no competing interests.

## STAR★METHODS

Detailed methods are provided in the online version of this paper and include the following:

* KEY RESOURCES TABLE

* METHOD DETAILS

* Network model

* Population coding and optimal decoder

* Analytical derivations

* Influence of correlated output noise on population coding

* Influence of recurrent connectivity on population coding

* Data analysis and model fitting

* Dimensionality of a compressed neural code

* Model extensions

* Generality and use of the network model

* Limitations of the study (extended)

* QUANTIFICATION AND STATISTICAL ANALYSIS

* Details of numerical simulations

* EXPERIMENTAL MODEL AND STUDY PARTICIPANT DETAILS

## SUPPLEMENTAL INFORMATION

Supplemental information can be found online at https://doi.org/10.1016/j.celrep.2025.115412.

Received: December 31, 2021

Revised: November 20, 2023

Accepted: February 18, 2025

Published: March 19, 2025

## REFERENCES

1. Georgopoulos, A.P., Kalaska, J.F., Caminiti, R., and Massey, J.T. (1982). On the relations between the direction of two-dimensional arm movements and cell discharge in primate motor cortex. J. Neurosci. 2, 1527-1537. https://doi.org/10.1523/jneurosci.02-11-01527.1982.

2. Taube, J.S., Muller, R.U., and Ranck, J.B. (1990). Head-Direction Cells Recorded from the Postsubiculum in Freely Moving Rats. I. Description and Quantitative Analysis. J. Neurosci. 10, 420-435. https://doi.org/10.1523/jneurosci.10-02-00420.1990.

3. Miller, J.P., Jacobs, G.A., and Theunissen, F.E. (1991). Representation of Sensory Information in the Cricket Cercal Sensory System. I. Response Properties of the Primary Interneurons. J. Neurophysiol. 66, 1680-1689. https://doi.org/10.1152/jn.1991.66.5.1680.

4. Bremmer, F., Ilg, U.J., Thiele, A., Distler, C., and Hoffmann, K.P. (1997). Eye Position Effects in Monkey Cortex. I. Visual and Pursuit-Related Activity in Extrastriate Areas MT and MST. J. Neurophysiol. 77, 944-961. https://doi.org/10.1152/jn.1997.77.2.944.

5. Dayan, P., and Abbott, L.F. (2001). Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems (MIT Press).

6. Kayaert, G., Biederman, I., Op De Beeck, H.P., and Vogels, R. (2005). Tuning for shape dimensions in macaque inferior temporal cortex. Eur. J. Neurosci. 22, 212-224. https://doi.org/10.1111/j.1460-9568.2005.04202.x.

7. Barlow, H.B. (2012). Possible Principles Underlying the Transformations of Sensory Messages. Sensory Communication 1, 216-234. https://doi.org/10.7551/mitpress/9780262518420.003.0013.

8. Atick, J.J., and Redlich, A.N. (1990). Towards a Theory of Early Visual Processing. Neural Comput. 2, 308-320. https://doi.org/10.1162/neco.1990.2.3.308.

9. Lewicki, M.S. (2002). Efficient coding of natural sounds. Nat. Neurosci. 5, 356-363. https://doi.org/10.1038/nn831.

10. Zhang, K., and Sejnowski, T.J. (1999). Neuronal Tuning: To Sharpen or Broaden? Neural Comput. 11, 75-84. https://doi.org/10.1162/089976699300016809.

11. Bethge, M., Rotermund, D., and Pawelzik, K. (2002). Optimal Short-Term Population Coding: When Fisher Information Fails. Neural Comput. 14, 2317-2351. https://doi.org/10.1162/08997660260293247.

12. Wang, Z., Stocker, A.A., and Lee, D.D. (2016). Efficient Neural Codes That Minimize Lp Reconstruction Error. Neural Comput. 28, 2656-2686. https://doi.org/10.1162/NECO_a_00900.
13. Deneve, S., Latham, P.E., and Pouget, A. (1999). Reading population codes: A neural implementation of ideal observers. Nat. Neurosci. 2, 740-745. https://doi.org/10.1038/11205.

14. Yaeli, S., and Meir, R. (2010). Error-based analysis of optimal tuning functions explains phenomena observed in sensory neurons. Front. Comput. Neurosci. 4, 130-146. https://doi.org/10.3389/fncom.2010.00130.

15. Ganguli, D., and Simoncelli, E.P. (2014). Efficient Sensory Encoding and Bayesian Inference with Heterogeneous Neural Populations. Neural Comput. 26, 2103-2134. https://doi.org/10.1162/NECO_a_00638.

16. Fiscella, M., Franke, F., Farrow, K., Müller, J., Roska, B., da Silveira, R.A., and Hierlemann, A. (2015). Visual Coding with a Population of Direction-Selective Neurons. J. Neurophysiol. 114, 2485-2499. https://doi.org/10.1152/jn.00919.2014.

17. Seung, H.S., and Sompolinsky, H. (1993). Simple models for reading neuronal population codes. Proc. Natl. Acad. Sci. USA 90, 10749-10753. https://doi.org/10.1073/pnas.90.22.10749.

18. Berens, P., Ecker, A.S., Gerwinn, S., Tolias, A.S., and Bethge, M. (2011). Reassessing optimal neural population codes with neurometric functions. Proc. Natl. Acad. Sci. USA 108, 4423-4428. https://doi.org/10.1073/pnas.1015904108.

19. Kim, J.H.J., Fiete, I., and Schwab, D.J. (2020). Superlinear Precision and Memory in Simple Population Codes. Preprint at arXiv. http://arxiv.org/abs/2008.00629.

20. Hafting, T., Fyhn, M., Molden, S., Moser, M.B., and Moser, E.I. (2005). Microstructure of a spatial map in the entorhinal cortex. Nature 436, 801-806. https://doi.org/10.1038/nature03721.

21. Doeller, C.F., Barry, C., and Burgess, N. (2010). Evidence for grid cells in a human memory network. Nature 463, 657-661. https://doi.org/10.1038/nature08704.

22. Yartsev, M.M., Witter, M.P., and Ulanovsky, N. (2011). Grid cells without theta oscillations in the entorhinal cortex of bats. Nature 479, 103-107. https://doi.org/10.1038/nature10583.

23. Killian, N.J., Jutras, M.J., and Buffalo, E.A. (2012). A map of visual space in the primate entorhinal cortex. Nature 491, 761-764. https://doi.org/10.1038/nature11587.

24. Ginosar, G., Aljadeff, J., Burak, Y., Sompolinsky, H., Las, L., and Ulanovsky, N. (2021). Locally ordered representation of 3D space in the entorhinal cortex. Nature 596, 404-409. https://doi.org/10.1038/s41586-021-03783-x.

25. Grieves, R.M., Jedidi-Ayoub, S., Mishchanchuk, K., Liu, A., Renaudineau, S., Duvelle, É., and Jeffery, K.J. (2021). Irregular distribution of grid cell firing fields in rats exploring a 3D volumetric space. Nat. Neurosci. 24, 1567-1573. https://doi.org/10.1038/s41593-021-00907-4.

26. Diehl, G.W., Hon, O.J., Leutgeb, S., and Leutgeb, J.K. (2017). Grid and Nongrid Cells in Medial Entorhinal Cortex Represent Spatial Location and Environmental Features with Complementary Coding Schemes. Neuron 94, 83-92.e6. https://doi.org/10.1016/j.neuron.2017.03.004.

27. Hardcastle, K., Maheswaranathan, N., Ganguli, S., and Giocomo, L.M. (2017). A Multiplexed, Heterogeneous, and Adaptive Code for Navigation in Medial Entorhinal Cortex. Neuron 94, 375-387.e7. https://doi.org/10.1016/j.neuron.2017.03.025.

28. Eliav, T., Maimon, S.R., Aljadeff, J., Tsodyks, M., Ginosar, G., Las, L., and Ulanovsky, N. (2021). Multiscale representation of very large environments in the hippocampus of flying bats. Science 372, eabg4020. https://doi.org/10.1126/science.abg4020.

29. Harland, B., Contreras, M., Souder, M., and Fellous, J.M. (2021). Dorsal CA1 hippocampal place cells form a multi-scale representation of megaspace. Curr. Biol. 31, 2178-2190.e6. https://doi.org/10.1016/j.cub.2021.03.003.

30. Gaucher, Q., Panniello, M., Ivanov, A.Z., Dahmen, J.C., King, A.J., and Walker, K.M. (2020). Complexity of frequency receptive fields predicts tonotopic variability across species. Elife 9, e53462. https://doi.org/10.7554/eLife.53462.

31. Kadia, S.C., and Wang, X. (2003). Spectral Integration in A1 of Awake Primates: Neurons with Single- and Multipeaked Tuning Characteristics. J. Neurophysiol. 89, 1603-1622. https://doi.org/10.1152/jn.00271.2001.

32. Sofroniew, N.J., Vlasov, Y.A., Hires, S.A., Freeman, J., and Svoboda, K. (2015). Neural coding in barrel cortex during whisker-guided locomotion. Elife 4, e12559. https://doi.org/10.7554/eLife.12559.

33. Lalazar, H., Abbott, L.F., and Vaadia, E. (2016). Tuning Curves for Arm Posture Control in Motor Cortex Are Consistent with Random Connectivity. PLoS Comput. Biol. 12, e1004910-e1004927. https://doi.org/10.1371/journal.pcbi.1004910.

34. Fiete, I.R., Burak, Y., and Brookings, T. (2008). What Grid Cells Convey about Rat Location. J. Neurosci. 28, 6858-6871. https://doi.org/10.1523/JNEUROSCI.5684-07.2008.

35. Sreenivasan, S., and Fiete, I. (2011). Grid cells generate an analog error-correcting code for singularly precise neural computation. Nat. Neurosci. 14, 1330-1337. https://doi.org/10.1038/nn.2901.

36. Mathis, A., Herz, A.V.M., and Stemmler, M.B. (2012). Resolution of Nested Neuronal Representations Can Be Exponential in the Number of Neurons. Phys. Rev. Lett. 109, 018103. https://doi.org/10.1103/PhysRevLett.109.018103.

37. Burak, Y. (2014). Spatial coding and attractor dynamics of grid cells in the entorhinal cortex. Curr. Opin. Neurobiol. 25, 169-175. https://doi.org/10.1016/j.conb.2014.01.013.

38. Wei, X.-X., Prentice, J., and Balasubramanian, V. (2015). A principle of economy predicts the functional architecture of grid cells. Elife 4, e08362. https://doi.org/10.7554/eLife.08362.

39. Welinder, P.E., Burak, Y., and Fiete, I.R. (2008). Grid cells: The position code, neural network models of activity, and the problem of learning. Hippocampus 18, 1283-1300. https://doi.org/10.1002/hipo.20519.

40. Mainali, N., da Silveira, R.A., and Burak, Y. (2025). Universal statistics of hippocampal place fields across species and dimensionalities. Neuron, in press. https://doi.org/10.1101/2024.06.11.597569.

41. Van Hateren, J.H., and Ruderman, D.L. (1998). Independent component analysis of natural image sequences yields spatio-temporal filters similar to simple cells in primary visual cortex. Proc. Biol. Sci. 265, 2315-2320. https://doi.org/10.1098/rspb.1998.0577.

42. Doi, E., Gauthier, J.L., Field, G.D., Shlens, J., Sher, A., Greschner, M., Machado, T.A., Jepson, L.H., Mathieson, K., Gunning, D.E., et al. (2012). Efficient Coding of Spatial Information in the Primate Retina. J. Neurosci. 32, 16256-16264. https://doi.org/10.1523/JNEUROSCI.4036-12.2012.

43. Zhaoping, L. (2014). Understanding Vision: Theory, Models, and Data (Oxford University Press).

44. Candes, E.J., and Tao, T. (2006). Near-optimal signal recovery from random projections: Universal encoding strategies? IEEE Trans. Inf. Theor. 52, 5406-5425. https://doi.org/10.1109/TIT.2006.885507.

45. Baraniuk, R., Davenport, M., DeVore, R., and Wakin, M. (2008). A simple proof of the restricted isometry property for random matrices. Constr. Approx. 28, 253-263. https://doi.org/10.1007/s00365-007-9003-x.

46. Brunel, N., and Nadal, J.P. (1998). Mutual Information, Fisher Information, and Population Coding. Neural Comput. 10, 1731-1757. https://doi.org/10.1162/089976698300017115.

47. Finkelstein, A., Ulanovsky, N., Tsodyks, M., and Aljadeff, J. (2018). Optimal dynamic coding by mixed-dimensionality neurons in the head-direction system of bats. Nat. Commun. 9, 3590. https://doi.org/10.1038/s41467-018-05562-1.

48. Kettner, R.E., Schwartz, A.B., and Georgopoulos, A.P. (1988). Primate motor cortex and free arm movements to visual targets in three-dimensional space. III. Positional gradients and population coding of movement direction from various movement origins. J. Neurosci. 8, 2938-2947. https://doi.org/10.1523/jneurosci.08-08-02938.1988.
49. Wang, W., Chan, S.S., Heldman, D.A., and Moran, D.W. (2007). Motor Cortical Representation of Position and Velocity During Reaching. J. Neurophysiol. 97, 4258-4270. https://doi.org/10.1152/jn.01180.2006.

50. Andersen, R.A., Essick, G.K., and Siegel, R.M. (1985). Encoding of spatial location by posterior parietal neurons. Science 230, 456-458. https://doi.org/10.1126/science.4048942.

51. Seung, H.S., and Lee, D.D. (2000). The manifold ways of perception. Science 290, 2268-2269. https://doi.org/10.1126/science.290.5500.2268.

52. Gallego, J.A., Perich, M.G., Miller, L.E., and Solla, S.A. (2017). Neural Manifolds for the Control of Movement. Neuron 94, 978-984. https://doi.org/10.1016/j.neuron.2017.05.025.

53. Fusi, S., Miller, E.K., and Rigotti, M. (2016). Why neurons mix: high dimensionality for higher cognition. Curr. Opin. Neurobiol. 37, 66-74. https://doi.org/10.1016/j.conb.2016.01.010.

54. Stringer, C., Pachitariu, M., Steinmetz, N., Carandini, M., and Harris, K.D. (2019). High-dimensional geometry of population responses in visual cortex. Nature 571, 361-365. https://doi.org/10.1038/s41586-019-1346-5.

55. Cunningham, J.P., and Yu, B.M. (2014). Dimensionality reduction for large-scale neural recordings. Nat. Neurosci. 17, 1500-1509. https://doi.org/10.1038/nn.3776.

56. Harel, Y., and Meir, R. (2020). Optimal Multivariate Tuning with Neuron-Level and Population-Level Energy Constraints. Neural Comput. 32, 794-828. https://doi.org/10.1162/neco_a_01267.

57. Mosheiff, N., Agmon, H., Moriel, A., and Burak, Y. (2017). An efficient coding theory for a dynamic trajectory predicts non-uniform allocation of entorhinal grid cells to modules. PLoS Comput. Biol. 13, e1005597. https://doi.org/10.1371/journal.pcbi.1005597.

58. Shannon, C.E. (1949). Communication in the Presence of Noise. Proc. IRE 37, 10-21. https://doi.org/10.1162/NECO_a_00991.

59. Montemurro, M.A., and Panzeri, S. (2006). Optimal Tuning Widths in Population Coding of Periodic Variables. Neural Comput. 18, 1555-1576. https://doi.org/10.1162/neco.2006.18.7.1555.

60. Wei, X.-X., and Stocker, A.A. (2012). Efficient coding provides a direct link between prior and likelihood in perceptual Bayesian inference. Adv. Neural Inf. Process. Syst. 25, 1304-1312.

61. Yerxa, T.E., Kee, E., DeWeese, M.R., and Cooper, E.A. (2020). Efficient sensory coding of multidimensional stimuli. PLoS Comput. Biol. 16, e1008146. https://doi.org/10.1371/journal.pcbi.1008146.

62. Wei, X.-X., Prentice, J., and Balasubramanian, V. (2013). The Sense of Place: Grid Cells in the Brain and the Transcendental Number e. Preprint at arXiv. http://arxiv.org/abs/1304.0031%0A.

63. Wu, A., Roy, N.A., Keeley, S., and Pillow, J.W. (2017). Gaussian process based nonlinear latent structure discovery in multivariate spike train data. Adv. Neural Inf. Process. Syst. 30, 3496-3505.

64. Wu, A., Pashkovski, S.L., Datta, S.R., and Pillow, J.W. (2018). Learning a latent manifold of odor representations from neural responses in piriform cortex. Adv. Neural Inf. Process. Syst. 31, 5378-5388.

65. Jensen, K.T., Kao, T.C., Tripodi, M., and Hennequin, G. (2020). Manifold GPLVMs for discovering non-Euclidean latent structure in neural data. Adv. Neural Inf. Process. Syst. 34, 22580-22592.

66. Wilke, S.D., and Eurich, C.W. (2002). Representational Accuracy of Stochastic Neural Populations. Neural Comput. 14, 155-189. https://doi.org/10.1162/089976602753284482.

67. Shamir, M., and Sompolinsky, H. (2006). Implications of Neuronal Diversity on Population Coding. Neural Comput. 18, 1951-1986. https://doi.org/10.1162/neco.2006.18.8.1951.

68. Berry, Ii, M.J., Lebois, F., Ziskind, A., and da Silveira, R.A. (2019). Functional Diversity in the Retina Improves the Population Code. Neural Comput. 31, 270-311. https://doi.org/10.1162/neco_a_01158.

69. Saxena, S., and Cunningham, J.P. (2019). Towards the neural population doctrine. Curr. Opin. Neurobiol. 55, 103-111. https://doi.org/10.1016/j.conb.2019.02.002.

70. Abbott, L.F., Rolls, E.T., and Tovee, M.J. (1996). Representational Capacity of Face Coding in Monkeys. Cerebr. Cortex 6, 498-505. https://doi.org/10.1093/cercor/6.3.498.

71. Kobak, D., Pardo-Vazquez, J.L., Valente, M., Machens, C.K., and Renart, A. (2019). State-dependent geometry of population activity in rat auditory cortex. Elife 8, e44526-e44527. https://doi.org/10.7554/eLife.44526.

72. Kriegeskorte, N., and Wei, X.-X. (2021). Neural tuning and representational geometry. Nat. Rev. Neurosci. 22, 703-718. https://doi.org/10.1038/s41583-021-00502-3.

73. De, A., and Chaudhuri, R. (2023). Common population codes produce extremely nonlinear neural manifolds. Proc. Natl. Acad. Sci. USA 120, e2305853120. https://doi.org/10.1073/pnas.2305853120.

74. Donoho, D.L. (2006). Compressed sensing. IEEE Trans. Inf. Theor. 52, 1289-1306. https://doi.org/10.1109/TIT.2006.871582.

75. Sharpee, T.O., and Berkowitz, J.A. (2019). Linking neural responses to behavior with information-preserving population vectors. Curr. Opin. Behav. Sci. 29, 37-44. https://doi.org/10.1016/j.cobeha.2019.03.004.

76. Kouh, M., and Poggio, T. (2008). A Canonical Neural Circuit for Cortical Nonlinear Operations. Neural Comput. 20, 1427-1451. https://doi.org/10.1162/neco.2008.02.07-466.

77. Carandini, M., and Heeger, D.J. (2011). Normalization as a canonical neural computation. Nat. Rev. Neurosci. 13, 51-62. https://doi.org/10.1038/nrn3136.

78. Bordelon, B., Canatar, A., and Pehlevan, C. (2020). Spectrum Dependent Learning Curves in Kernel Regression and Wide Neural Networks. In International Conference of Machine Learning, 37, pp. 1024-1034. https://doi.org/10.48550/arXiv.2002.02561.

79. Bordelon, B., and Pehlevan, C. (2022). Population codes enable learning from few examples by shaping inductive bias. Elife 11, e78606. https://doi.org/10.7554/eLife.78606.

80. Moreno-Bote, R., Beck, J., Kanitscheider, I., Pitkow, X., Latham, P., and Pouget, A. (2014). Information-limiting correlations. Nat. Neurosci. 17, 1410-1417. https://doi.org/10.1038/nn.3807.

81. Rumyantsev, O.I., Lecoq, J.A., Hernandez, O., Zhang, Y., Savall, J., Chrapkiewicz, R., Li, J., Zeng, H., Ganguli, S., and Schnitzer, M.J. (2020). Fundamental bounds on the fidelity of sensory cortical coding. Nature 580, 100-105. https://doi.org/10.1038/s41586-020-2130-2.

82. Blanco Malerba, S., Micheli, A., Woodford, M., and Azeredo da Silveira, R. (2024). Jointly efficient encoding and decoding in neural populations. PLoS Comput. Biol. 20, e1012240. https://doi.org/10.1371/journal.pcbi.1012240.

83. Rasmussen, C.E., and Williams, C.K.I. (2005). Gaussian Processes for Machine Learning (MIT press). https://doi.org/10.7551/mitpress/3206.001.0001.

84. Tarmast, G. (2001). Multivariate log-normal distribution. In International Statistical Institute: Seoul 53rd Session 210, 1-2. https://doi.org/isiproceedings.org/pdf/329.PDF.

85. Livan, G., Novaes, M., and Vivo, P. (2018). Introduction to Random Matrices - Theory and Practice (Cham: Springer). https://doi.org/10.1007/978-3-319-70885-0.

86. Arakaki, T., Barello, G., and Ahmadian, Y. (2019). Inferring neural circuit structure from datasets of heterogeneous tuning curves. PLoS Comput. Biol. 15, e1006816. https://doi.org/10.1371/journal.pcbi.1006816.

87. Everitt, B.S., and Skrondal, A. (2018). The Cambridge Dictionary of Statistics, 4 ed. (Cambridge University Press).

88. Gao, P., Trautmann, E., Yu, B., Santhanam, G., Ryu, S., and Ganguli, S. (2017). A theory of multineuronal dimensionality, dynamics and measurement. Preprint at bioRxiv. https://doi.org/10.1101/214262.
<figure><img src="image_21.png" /></figure>

89. Azeredo da Silveira, R., and Rieke, F. (2021). The Geometry of Information Coding in Correlated Neural Populations. Annu. Rev. Neurosci. 44, 403-424. https://doi.org/10.1146/annurev-neuro-120320-082744.

90. Thorpe, S., Fize, D., and Marlot, C. (1996). Speed of processing in the human visual system. Nature 381, 520-522. https://doi.org/10.1038/381520a0.

91. Masquelier, T., and Thorpe, S.J. (2007). Unsupervised learning of visual features through spike timing dependent plasticity. PLoS Comput. Biol. 3, e31. https://doi.org/10.1371/journal.pcbi.0030031.

92. Barak, O., Rigotti, M., and Fusi, S. (2013). The sparseness of mixed selectivity neurons controls the generalization-discrimination trade-off. J. Neurosci. 33, 3844-3856. https://doi.org/10.1523/JNEUROSCI.2753-12.2013.

93. Babadi, B., and Sompolinsky, H. (2014). Sparseness and Expansion in Sensory Representations. Neuron 83, 1213-1226. https://doi.org/10.1016/j.neuron.2014.07.035.

94. Tootoonian, S., and Lengyel, M. (2014). A Dual Algorithm for Olfactory Computation in the Locust Brain. Adv. Neural Inf. Process. Syst. 27, 2276-2284.

95. Zhang, Y., and Sharpee, T.O. (2016). A robust feedforward model of the olfactory system. PLoS Comput. Biol. 12, e1004850.

96. Litwin-Kumar, A., Harris, K.D., Axel, R., Sompolinsky, H., and Abbott, L.F. (2017). Optimal Degrees of Synaptic Connectivity. Neuron 93, 1153-1164.e7. https://doi.org/10.1016/j.neuron.2017.01.030.

97. Stettler, D.D., and Axel, R. (2009). Representations of Odor in the Piriform Cortex. Neuron 63, 854-864. https://doi.org/10.1016/j.neuron.2009.09.005.

98. Van Vreeswijk, C., and Sompolinsky, H. (1996). Chaos in neuronal networks with balanced excitatory and inhibitory activity. Science 274, 1724-1726. https://doi.org/10.1126/science.274.5293.1724.

99. Renart, A., De La Rocha, J., Bartho, P., Hollender, L., Parga, N., Reyes, A., and Harris, K.D. (2010). The Asynchronous State in Cortical Circuits. Science 327, 587-590. https://doi.org/10.1126/science.1179850.

100. Takács, V.T., Klausberger, T., Somogyi, P., Freund, T.F., and Gulyás, A.I. (2012). Extrinsic and local glutamatergic inputs of the rat hippocampal CA1 area differentially innervate pyramidal cells and interneurons. Hippocampus 22, 1379-1391. https://doi.org/10.1002/hipo.20974.

101. Bezanson, J., Edelman, A., Karpinski, S., and Shah, V.B. (2017). Julia: A fresh approach to numerical computing. SIAM Rev. 59, 65-98. https://doi.org/10.1137/141000671.

102. Datseris, G., Isensee, J., Pech, S., and Gál, T. (2020). DrWatson: the perfect sidekick for your scientific inquiries. J. Open Source Softw. 5, 2673. https://doi.org/10.21105/joss.02673.
## STAR★METHODS

### KEY RESOURCES TABLE

<table><thead><tr><th>REAGENT or RESOURCE</th><th>SOURCE</th><th>IDENTIFIER</th></tr></thead><tbody><tr><td colspan="3">Deposited data</td></tr><tr><td>Multi-electrode recordings from primate M1 during arm posture task</td><td>Lalazar et al.<sup>33</sup></td><td>https://osf.io/u57df</td></tr><tr><td colspan="3">Software and algorithms</td></tr><tr><td>Julia 1.0</td><td>Bezanson et al.<sup>101</sup></td><td>https://julialang.org/</td></tr><tr><td>DrWatson</td><td>Datseris et al. 2020<sup>102</sup></td><td>https://juliadynamics.github.io/DrWatson.jl/dev/</td></tr><tr><td>Custom code</td><td>https://github.com/simonebmalerba/RandomCompressedCoding.git</td><td>https://doi.org/10.5281/zenodo.14210005</td></tr></tbody></table>

### METHOD DETAILS

Throughout, we denote vectors by bold letters, e.g., $\mathbf{r} = (r_1, r_2, ..., r_N)$, and the $L_2$ norm as $|| \mathbf{r} ||_2^2 = \sum_i r_i^2$. Capital bold letters, e.g., $\mathbf{W}$, refer to matrices. We denote the derivative of a function as $f'(x) = \partial f / \partial x$.

#### Network model

#### Network model for one-dimensional stimuli and constraints on its parameters

We consider a two-layer feedforward network. The first, sensory layer is made up of $L$ neurons, each responding to a continuous scalar stimulus, $x \in [0, 1]$, according to a Gaussian tuning curve. The mean activity of neuron $j$ in response to a stimulus, $x$, is given by

$$u_j(x) = A \exp \left( - \frac{(x - c_j)^2}{2\sigma^2} \right), \qquad (\text{Equation 10})$$

where $c_j$ is the preferred stimulus of neuron $j$, $\sigma$ is the tuning-curve width, and $A$ is a fixed response amplitude. The preferred stimuli are evenly spaced, $c_j = j/L$. Each neuron in the first layer projects onto all $N$ neurons in the second, representation layer. The response function is assumed to be linear, and the random synaptic weights are independent realizations of a Gaussian random variable, $W_{ij} \sim N(0, 1/L)$; hence, the mean activity of representation neuron $i$ can be written as

$$v_i(x) = \sum_{j=1}^{L} W_{ij} u_j(x). \qquad (\text{Equation 11})$$

The value of the amplitude, $A$, is chosen so as to set the 'dynamic range' of representation neurons to a fixed value; more precisely, we choose the value of $A$ so that the variance of each neuron's response across the stimulus range is invariant under variations in the other parameter of the model, $\sigma$, on average over network realizations. This quantity is calculated as

$$ 
\begin{aligned}
R &= \left\langle \int_{0}^{1} dx \left[ v_{i}(x) - \int_{0}^{1} dx' v_{i}(x') \right]^{2} \right\rangle_{W} \\
&= \left\langle \int_{0}^{1} dx v_{i}(x)^{2} - \left( \int_{0}^{1} dx v_{i}(x) \right)^{2} \right\rangle_{W} \\
&= \left\langle \sum_{j=1}^{L} W_{ij} W_{ij'} \left[ \left( \int_{0}^{1} dx u_{j}(x) u_{j'}(x) \right) - \left( \int_{0}^{1} dx u_{j}(x) \right) \left( \int_{0}^{1} dx u_{j'}(x) \right) \right] \right\rangle_{W} \\
&= \int_{0}^{1} dx u_{j}(x)^{2} - \left( \int_{0}^{1} dx u_{j}(x) \right)^{2},
\end{aligned} 
\qquad (\text{Equation 12})
$$
where $\langle \cdot \rangle_W$ indicates an average over the distribution of synaptic weights. Here (and below), we approximate Gaussian integrals on a bounded domain as

$$ \int_{0}^{1} dx u_{j}(x) \approx \int_{-\infty}^{\infty} dx u_{j}(x) = A \sqrt{2\pi\sigma^{2}}; \qquad (\text{Equation 13}) $$

this approximation is valid when $\sigma$ is small with respect to the stimulus range and $c_j$ is separated from the boundaries (0 and 1) by a distance that exceeds $\sigma$. As we will consider a large number of neurons in the sensory layer and relatively small values of $\sigma$ (up to 1/10 th of the stimulus range), errors introduced by this approximation will be negligible. By inserting Equation 13 and a similar approximation for $\int_0^1 dx u_j(x)^2$ into Equation 3, we obtain $A$ as a function of $\sigma$, as

$$ A^2 = \frac{R}{\sqrt{\pi\sigma^2} - 2\pi\sigma^2}. \qquad (\text{Equation 14}) $$

## Tuning curves as samples from a Gaussian process

The response of each neuron in the second layer to a stimulus, $x$, is a sum of realizations of Gaussian random variables; as a result, it is also a realization of a Gaussian random variable, with mean

$$ \langle v_i(x) \rangle_W = \sum_{j=1}^{L} \langle W_{ij} \rangle_W u_j(x) = 0, \qquad (\text{Equation 15}) $$

and its covariance is calculated as

$$ \begin{aligned} \langle v_i(x) v_j(x') \rangle_W &= \sum_{j' = 1}^{L} \langle W_{ij} W_{j'k'} \rangle_W u_j(x) u_{j'}(x') \\ &= \sum_{j=1}^{L} \frac{1}{L} u_j(x) u_j(x') \\ &\approx A^2 \int_0^1 d c_j \exp \times \left( - \frac{(x - c_j)^2 + (x' - c_j')^2}{2\sigma^2} \right) \\ &\approx A^2 \sqrt{\pi\sigma^2} \exp \left( - \frac{\Delta x^2}{4\sigma^2} \right), \end{aligned} \qquad (\text{Equation 16}) $$

where $\Delta x = x - x'$. The first approximation is obtained by replacing a sum by an integral $\sum_{j=1}^L \frac{1}{L} f(c_j) \approx \int_0^1 f(c_j) dc_j$ and the second approximation consists in extending the integration domain to the entire real line. The first approximation is valid if the spacing between the centers is small relatively to the width of the Gaussian, that is $L\sigma \gg 1$, while the second is valid if the arithmetic mean of $x$ and $x'$ is far from the stimulus boundaries. According to Equations 15 and 16 each neuron's tuning curve can be viewed as a sample from a one-dimensional Gaussian process with vanishing mean and Gaussian correlation function with standard deviation equal to $\sqrt{2\sigma}$. As the expected number of zero-crossings of a Gaussian process in the unit interval decreases as the inverse of the width of the correlation function (specifically, as $1/2\sqrt{2\sigma}^{83}$), for large values of $\sigma$ the tuning curves become unimodal.

## Tuning curves in the presence of an exponential response non-linearity

We consider the case in which an exponential non-linear function is applied to the linear combination of activity of the sensory neurons. The tuning curve of neuron $i$ is then obtained as

$$ v_i(x) = A_{nl} \exp \left( \sum_{j=1}^{L} W_{ij} u_j(x) \right). \qquad (\text{Equation 17}) $$

Since the arguments of the non-linearity, the linear tuning curves $v'_i(x) = \sum W_{ij} u_j(x)$, are samples from a Gaussian process (and, thus, neural responses are distributed according to a multivariate Gaussian), the joint distribution of non-linear neural responses is a multivariate log-normal distribution.⁸⁴ Its means is obtained as

$$ \langle v_i(x) \rangle_W = \exp \left( \mu + \frac{\langle v_i'(x)^2 \rangle_W}{2} \right), \qquad (\text{Equation 18}) $$

where $\mu = \log(A_{nl})$. The covariance is calculated as

$$ \begin{aligned} \langle (v_i(x) - \langle v_i(x) \rangle_W) (v_j(x') - \langle v_j(x') \rangle_W) \rangle_W &= \langle A_{nl}^2 \exp (v_i'(x) + v_j'(x')) \rangle_W - \langle A_{nl} \exp (v_i'(x)) \rangle_W^2 \\ &= \exp \left( 2\mu + \langle v_i'(x)^2 \rangle_W \right) \left( \exp \left( \langle v_i'(x) v_j'(x') \rangle_W \right) - 1 \right), \end{aligned} \qquad (\text{Equation 19}) $$

where have used the fact that $v'_i(x) + v'_j(x') \sim N(0, 2\langle v_i'(x)^2 \rangle_W + \langle v_j'(x) v_i'(x') \rangle_W)$, and the analog of Equation 18 to calculate the mean of the log-normal distribution. The covariance of the linear tuning curves, $\langle v'_i(x) v'_j(x') \rangle_W$, is obtained as in Equation 16; by setting $\Delta x = 0$
in Equation 16, we obtain the variance, $\langle v_i^\prime(x)^2 \rangle_W \approx 1$ (in the limit of small $\sigma$). The computation of $A_{nl}$ in the general case is difficult; however, when $\sigma \to 0$ and responses are close to independent, we obtain

$$A_{nl}^2 = \frac{R}{e(e - 1)}. \qquad (\text{Equation 20})$$

The correlation function, Equation 19, is plotted in Figure S6 for different values of $\Delta x = x - x'$.

**Network model for multi-dimensional stimuli**

We denote by $K$ the stimulus dimensionality, such that the stimulus is a $K$-dimensional vector, $\mathbf{x} = \{x_1, x_2, \dots, x_K\}$. Analogously to the one-dimensional case, each stimulus dimension can assume values in a bounded interval, $x_k \in [0, 1]$. We consider the two cases of pure and conjunctive tuning for sensory neurons. In both cases, the sensory layer is made up of $L$ neurons, which project onto all $N$ representation neurons. Similarly to the one-dimensional case, synaptic weights are independent realizations of a Gaussian random variable, $W_{ij} \sim N(0, 1/L)$.

**Sensory neurons with pure tuning.** The $L$ neurons are divided in $K$ sub-populations of $Q = L/K$ neurons. Neurons in the sub-population $k$ are sensitive to the single stimulus dimension $x_k$. The mean activity of neuron $j$ assigned to sub-population $k$ is given by the one-dimensional Gaussian

$$u_{j,k}^{p}(\mathbf{x}) = u_{j,k}^{p}(x_k) = A_p \exp \left( - \frac{(x_k - c_j^k)^2}{2\sigma^2} \right), \qquad (\text{Equation 21})$$

with preferred stimuli evenly spaced, $c_j^k = j/Q$ for $j = 1, \dots, Q$. The mean activity of representation neuron $i$ can be written as a superposition of one-dimensional tuning curves,

$$\begin{aligned} v_i^p(\mathbf{x}) &= \sum_{k=1}^{K} \sum_{j=1}^{Q} W_{ij,k} u_{j,k}(\mathbf{x}) \\ &= \sum_{k=1}^{K} v_{i,k}^p(x_k), \end{aligned} \qquad (\text{Equation 22})$$

each of which is an independent sample from a one-dimensional Gaussian process with statistics defined in Equations 15 and 16 and covariance scaled by a factor of $1/K$. Imposing the resource constraint, Equation 3, we obtain $A_p^2 = R/((\pi\sigma^2)^{1/2} - 2\pi\sigma^2)$.

**Sensory neurons with conjunctive tuning.** Neurons are sensitive to all stimulus dimensions. The mean activity of sensory neuron $j$ is given by the multi-dimensional Gaussian function

$$u_j^c(\mathbf{x}) = A_c \exp \left( - \frac{\| \mathbf{x} - \mathbf{c}_j \|_2^2}{2\sigma^2} \right), \qquad (\text{Equation 23})$$

with preferred stimuli, $\mathbf{c}_j$, arranged on a $K$-dimensional square grid with mesh size $L^{-1/K}$. The mean activity of representation neuron $i$ is obtained as

$$v_i^c(\mathbf{x}) = \sum_{j=1}^{L} W_{ij} u_j^c(\mathbf{x}). \qquad (\text{Equation 24})$$

By adapting Equations 15 and 16 to this case, we obtain that tuning curves are samples from a $K$-dimensional Gaussian process with vanishing mean and covariance

$$\langle v_i(\mathbf{x}) v_j(\mathbf{x}') \rangle_W = A_c^2 (\pi \sigma^2)^{K/2} \exp \left( - \frac{\Delta \mathbf{x}^2}{4\sigma^2} \right), \qquad (\text{Equation 25})$$

where $\Delta \mathbf{x} = \| \mathbf{x} - \mathbf{x}' \|_2^2$. Imposing the resource constraint, Equation 3, we obtain

$$A_c^2 = R / \left( (\pi\sigma^2)^{K/2} - (2\pi\sigma^2)^K \right)$$

## Population coding and optimal decoder

### Noise model

We assume that the response of representation neurons is corrupted by noise. The vector of responses to a given stimulus, $\mathbf{x}$, is

$$\mathbf{r} = \mathbf{v}(\mathbf{x}) + \mathbf{z}, \qquad (\text{Equation 26})$$
where **z** is a noise vector of independent Gaussian entries with vanishing mean and fixed variance, $z_i \sim N(0, \eta^2)$. Here, $\mathbf{v}(\mathbf{x}) = \{v_1(\mathbf{x}), v_2(\mathbf{x}), \dots, v_N(\mathbf{x})\}$ is the vector of mean responses of second-layer neurons to a stimulus, $\mathbf{x}$ (see Equation 2). The probability density of a response vector, $\mathbf{r}$, given a stimulus, $\mathbf{x}$, is written as

$$p(\mathbf{r}|\mathbf{x}) = \frac{1}{(2\pi\eta^2)^{N/2}} \exp\left(-\frac{\|\mathbf{r} - \mathbf{v}(\mathbf{x})\|_2^2}{2\eta^2}\right). \qquad (\text{Equation 27})$$

Below, we will furthermore consider an extension that takes into account a generic noise covariance matrix, $\Sigma$, resulting in the more general form

$$p(\mathbf{r}|\mathbf{x}) = \frac{1}{(2\pi)^{N/2}[\det(\Sigma)]^{1/2}} \exp\left(-\frac{1}{2}(\mathbf{r} - \mathbf{v}(\mathbf{x}))^T\mathbf{\Sigma}^{-1}(\mathbf{r} - \mathbf{v}(\mathbf{x}))\right). \qquad (\text{Equation 28})$$

### Loss function and ideal decoder

We quantify the coding performance of the neural population by the mean squared error (MSE) in the stimulus estimate,⁵ as obtained from the ideal decoder, or estimator, $\hat{\mathbf{x}} = f_{\text{dec}}(\mathbf{r})$, expressed as

$$E^2 = \int_0^1 dx \int dr p(\mathbf{r}|x)(\hat{\mathbf{x}} - x)^2, \qquad (\text{Equation 29})$$

where we have assumed a uniform prior over stimuli, $p(\mathbf{x}) \sim U(0,1)$. We consider the average of this quantity over network realizations, $\epsilon^2 = \langle E^2 \rangle_W$; in some figures, we plot the square root of this quantity, the RMSE, $\epsilon = \sqrt{\langle E^2 \rangle_W}$.

For multi-dimensional stimuli, the ideal decoder outputs a vector estimate of the stimulus, $\hat{\mathbf{x}} = f_{\text{dec}}(\mathbf{r})$. In this case, we define the MSE as the average squared norm of the difference between the stimulus and the decoder output,

$$E^2 = \int d\mathbf{x} \int dr p(\mathbf{r}|\mathbf{x}) \| \hat{\mathbf{x}} - \mathbf{x} \|_2^2, \qquad (\text{Equation 30})$$

where $\| \hat{\mathbf{x}} - \mathbf{x} \|_2^2 = \sum_{k=1}^K (\hat{x}_k - x_k)^2$.

The estimator that minimizes the MSE (Minimum-MSE or MMSE) is given by the mean of the posterior density. We can write the optimal estimator as

$$\hat{x}_{\text{MMSE}} = \int_{0}^{1} dx p(x|\mathbf{r}) x = \frac{\int_{0}^{1} dx p(\mathbf{r}|x)x}{\int_{0}^{1} dx p(\mathbf{r}|x)}. \qquad (\text{Equation 31})$$

We note that a simple neural network can output the MMSE estimate. Indeed, if we approximate the integrals in Equation 31 with a discrete sum over $M$ values and we use Equation 27, we obtain

$$\begin{aligned}
\hat{x}_{\text{MMSE}} & \approx \frac{\sum_{m=1}^{M} x_m p(\mathbf{r}|x_m) \Delta x_m}{\sum_{m=1}^{M} p(\mathbf{r}|x_m) \Delta x_m} \\
& = \frac{\sum_{m=1}^{M} x_m \exp\left(-\frac{1}{2\eta^2}\left(\sum_{i=1}^{N} r_i^2 + \sum_{i=1}^{N} v_i^2(x_m) - 2 \sum_{i=1}^{N} v_i(x_m) r_i\right)\right)}{\sum_{m=1}^{M} \exp\left(-\frac{1}{2\eta^2}\left(\sum_{i=1}^{N} r_i^2 + \sum_{i=1}^{N} v_i^2(x_m) - 2 \sum_{i=1}^{N} v_i(x_m) r_i\right)\right)} \\
& = \frac{\sum_{m=1}^{M} x_m \exp\left(\frac{1}{2\eta^2}\left(\sum_{i=1}^{N} 2v_i(x_m)r_i - \sum_{i=1}^{N} v_i^2(x_m)\right)\right)}{\sum_{m=1}^{M} \exp\left(\frac{1}{2\eta^2}\left(\sum_{i=1}^{N} 2v_i(x_m)r_i - \sum_{i=1}^{N} v_i^2(x_m)\right)\right)} \\
& = \sum_{m=1}^{M} x_m \tilde{h}_m,
\end{aligned} \qquad (\text{Equation 32})$$

where the terms $\sum_{i=1}^N r_i^2$ in both numerator and denominator cancel and we assumed a constant spacing,

$\Delta x_m = \Delta x_0$. The approximate estimate specified by Equation 32 can be produced by a two-layer neural network: a first layer of $M$ neurons, whose activities are given by

$$\tilde{h}_m = \frac{\exp\left(\sum_{i=1}^{N} \lambda_{mi} r_i + b_m\right)}{\sum_{m'=1}^{M} \exp\left(\sum_{i=1}^{N} \lambda_{m'mi} r_i + b_{m'}\right)}, \qquad (\text{Equation 33})$$

computes a normalized, discrete approximation of the posterior, $\tilde{h}_m \approx p(x_m|\mathbf{r})$, such that $\sum_{m=1}^M \tilde{h}_m = 1$. The unnormalized activity of neuron $m$, $h_m = \exp(\sum_{i=1}^N \lambda_{mi} r_i + b_m)$, is obtained as a linear combination of the activities of the representation neurons plus a bias
term, transformed through an exponential non-linearity. The ‘synaptic weight’ from the ith representation neuron to the mth decoder neuron is a function of the true mean response of neuron i to stimulus $x_m$ and of the variance of the noise, $\lambda_{mi} = v_i(x_m)/\eta^2$. Similarly, the bias term is obtained as $b_m = -\sum_i v_i(x_m)^2/2\eta^2$. Finally, to obtain the MMSE estimate, a single output neuron weights the activity of these M neurons according to their ‘preferred stimulus’, $x_m$.

In what follows, we will also use the maximum a posteriori (MAP) estimator, defined as

$$ \hat{x}_{MAP} = \arg \min_{x_m} \| \mathbf{r} - \mathbf{v}(x_m) \|_2^2. \quad (\text{Equation 34}) $$

It is equal to the maximum likelihood (ML) estimator given the uniformity of the stimulus prior, and it has a simple geometric interpretation: it identifies the stimulus which corresponds to the vector of mean responses closest to the noisy population output. In numerical simulations, the MSEs calculated with the MMSE and the MAP estimators are very similar.

The MMSE estimator can be extended to the case of non-diagonal noise covariance matrix, $\Sigma$, by combining Equations 28 and 31. The decoder weights and biases are then correlated, $\lambda_m = \mathbf{v}^T(x_m)\Sigma^{-1}$ and $b_m = \mathbf{v}^T(x_m)\Sigma^{-1}\mathbf{v}(x_m)$, where $\lambda_m$ denotes the vector with elements corresponding to the mth row of $\lambda$.

The MMSE estimator can also be extended to the case of multi-dimensional stimuli. In this case, the integrals of Equation 31 are K-dimensional and the output layer is made up by K neurons, which compute a vector estimate of the stimulus, $\hat{\mathbf{x}}$.

## Analytical derivations

In the calculations that follow, we consider the limit of $L \to \infty$ and we assume $N \ll L$.

### Narrow tuning curves

In the limiting case with $\sigma \to 0$, sensory neurons respond only to their preferred stimulus. Therefore, we consider the case of L discrete stimuli corresponding to the neurons' preferred stimuli, $x_j = j/L$. The mean activity of representation neuron i is written as $v_i(x_j) = \hat{A}W_{ij}$, with $\hat{A}^2 = LR$, such that $v_i(x_j) \sim N(0, R)$. The constant of proportionality is computed with the analog of Equation 3 for discrete stimuli, in the limit of large L.

The MSE in the case of narrow tuning curves, $\epsilon_n^2$, is obtained as

$$ \epsilon_n^2 = \langle E^2 \rangle_W = \left\langle \frac{1}{L} \sum_{j=1}^{L} \int d\mathbf{r} p_e(\mathbf{r}, \mathbf{W}, x_j) (\hat{x}_j - x_j)^2 \right\rangle_W, \quad (\text{Equation 35}) $$

where $p_e(\mathbf{r}, \mathbf{W}, x_j)$ denotes the probability, given a synaptic matrix $\mathbf{W}$ and noise, of having an incorrect estimate of $x_j$, i.e., $\hat{x}_j \neq x_j$. For every choice of $\mathbf{W}$ and $\mathbf{r}$, there are $L-1$ equiprobable realizations of the synaptic matrix which correspond to permutations of the identity of the decoded stimulus, such that $\hat{x}_j = x_j$ with $j' \neq j$. Therefore, the MSE can be written as

$$ \epsilon_n^2 = \langle P(E) \rangle_W \frac{1}{L(L-1)} \sum_{j=1}^{L} \sum_{j' \neq j} (x_j - x_{j'})^2, \quad (\text{Equation 36}) $$

where $\langle P(E) \rangle_W = \langle \int d\mathbf{r} p_e(\mathbf{r}, \mathbf{W}, x_j) \rangle_W$ is the probability of error averaged over the noise and the synaptic weights realizations. The MSE is the product of two terms: the mean probability of error and the average squared magnitude of the error. We now compute these two terms.

**Error probability.** An error occurs if there exists a $j'$ such that $\mathbf{r}$ is closer to $\mathbf{v}(x_{j'})$ than to $\mathbf{v}(x_j)$, where $x_j$ is the presented stimulus. We calculate the probability of error as a function of the probability of the complementary event, as

$$ P(E) = 1 - P(\| \mathbf{r} - \mathbf{v}(x_{j'}) \|_2^2 > \| \mathbf{r} - \mathbf{v}(x_j) \|_2^2 \quad \forall j \neq j'), \quad (\text{Equation 37}) $$

By averaging over different realizations of $\mathbf{W}$, $\langle P(E) \rangle_W$, the probabilities that an error is not committed on the possible values of $j'$ are independent; thus, we can express the probability of error, as a function of the mean responses, $v_i$, and the noise, $z_i$, as

$$ \begin{aligned} \langle P(E) \rangle_W &= 1 - \left(1 - \langle P(\| \mathbf{r} - \mathbf{v}(x_{j'}) \|_2^2 < \| \mathbf{r} - \mathbf{v}(x_j) \|_2^2) \rangle_W \right)^{L-1} \\ &\approx L \langle P(\| \mathbf{r} - \mathbf{v}(x_{j'}) \|_2^2 < \| \mathbf{r} - \mathbf{v}(x_j) \|_2^2) \rangle_W \\ &= L \langle P \left( \sum_{i=1}^N (v_i(x_j) - v_i(x_{j'}))^2 - \sum_{i=1}^N 2(v_i(x_j) - v_i(x_{j'}))z_i < 0 \right) \rangle_W. \end{aligned} \quad (\text{Equation 38}) $$

The approximation comes from the assumption that the probability of error is small, and $L-1 \approx L$ is large, while the last equality is obtained from the definition of noisy responses, Equation 26. The difference between the mean activity of the same neuron to two different stimuli is sampled according to a Gaussian distribution, $\hat{v}_i \equiv v_i(x_j) - v_i(x_{j'}) = \hat{A}(W_{ij} - W_{ij'}) \sim N(0, 2R)$. The mean probability of error is calculated as

$$ \langle P(E) \rangle_W \approx L \int \prod_{i=1}^{N} d\hat{v}_i \prod_{i=1}^{N} dz_i p(\hat{v}_i) p(z_i) \Theta \left( - \sum_{i=1}^{N} \hat{v}_i^2 + 2 \sum_{i=1}^{N} \hat{v}_i z_i \right). \quad (\text{Equation 39}) $$
This quantity is the probability that the random variable $\rho = \sum_{i=1}^{N} \hat{v}_i^2 - \sum_{i=1}^{N} 2\hat{v}_i z_i$, where $\hat{v}_i \sim N(0, 2R)$ and $z_i \sim N(0, \eta^2)$, is negative. With $\zeta = \sum_{i=1}^{N} \hat{v}_i^2$, the distribution of $\rho$ conditional on $\zeta$ is Gaussian with mean $\zeta$ and variance $4\zeta\eta^2$. Thus,

$$ \langle P(E) \rangle_W \approx L \int_0^\infty d\zeta p(\zeta) \int_{-\infty}^0 dp p(\rho|\zeta) \\ = \frac{L}{2} \int_0^\infty d\zeta p(\zeta) \text{erfc}\left(\sqrt{\frac{\zeta}{8\eta^2}}\right), \qquad (\text{Equation 40}) $$

where erfc is the complementary error function and

$$ p(\zeta) = \frac{\left(\frac{\zeta}{2R}\right)^{N/2-1} \exp\left(-\frac{\zeta}{4R}\right)}{2^{N/2+1} \Gamma(N/2)} \qquad (\text{Equation 41}) $$

is the probability density function of a chi-squared distribution. Computing the integral, we obtain

$$ \langle P(E) \rangle_W \approx L \frac{\left(\frac{\eta^2}{2R}\right)^{\frac{N}{2}} \Gamma(N)}{\Gamma\left(\frac{N}{2}\right)} {}_2F_1\left(\frac{N}{2}, \frac{1+N}{2}, \frac{2+N}{2}, -2\frac{\eta^2}{R}\right) \\ = L \frac{\left(\frac{\eta^2}{2R}\right)^{\frac{N}{2}} \Gamma(N)}{\Gamma\left(\frac{N}{2}\right) \Gamma\left(\frac{2+N}{2}\right)} \sum_{n=0}^{\infty} \frac{\left(\frac{N}{2}\right)_n \left(\frac{N+1}{2}\right)_n}{\left(\frac{N+2}{2}\right)_n n!} \left(-2\frac{\eta^2}{R}\right)^n, \qquad (\text{Equation 42}) $$

where ${}_2F_1(a, b, c, x)$ is the regularized 2F1 Hypergeometric function; we provide its definition in the last equality. The Pochammer symbol can be defined through Gamma functions, $(x)_n = \frac{\Gamma(x+n)}{\Gamma(x)}$. In order to obtain a closed-form expression, we make the approximation $(N+1)/2 \approx (N+2)/2$. By using the identity $\sum_{n=0}^{\infty} \frac{(x)_n}{n!} a^n = (1-a)^{-x}$ and the Stirling approximation for Gamma functions, we obtain the expression of the error probability that appears in the main text, Equation 5:

$$ \langle P(E) \rangle_W \approx L \left(\frac{\eta^2}{2R}\right)^{\frac{N}{2}} \frac{\Gamma(N)}{\Gamma^2\left(\frac{N}{2}\right) \frac{N}{2} \left(1+2\frac{\eta^2}{R}\right)^{\frac{N+1}{2}}} \qquad (\text{Equation 43}) \\ \approx \frac{L}{\sqrt{2\pi N}} \exp\left(-\log\left(1+\frac{R}{2\eta^2}\right)\frac{N}{2}\right). $$

### Average squared magnitude of error

We denote by $\bar{\epsilon}_{n,g}^2$ the second factor in Equation 36, which can be written as

$$ \bar{\epsilon}_{n,g}^2 = \frac{1}{L(L-1)} \sum_{j=1}^{L} \sum_{j' \neq j}^{L} \left(\frac{j}{L} - \frac{j'}{L}\right)^2 \\ = \frac{1}{L(L-1)} \left( \sum_{j=1}^{L} \sum_{j' \neq j}^{L} j^2 + \sum_{j=1}^{L} \sum_{j' \neq j}^{L} \frac{j'^2}{L^2} - 2 \sum_{j=1}^{L} \sum_{j'=1}^{L} \frac{j j'}{L^2} \right). \qquad (\text{Equation 44}) $$

These sums can be computed through the identities for the sum of the first $n$ squared numbers, $\sum_{i=1}^{n} i^2 = n(n+1)(2n+1)/6$, and for the sum of the first $n$ numbers, $\sum_{i=1}^{n} i = n(n+1)/2$. The first two sums in Equation 44 are identical, and yield

$$ \sum_{j=1}^{L} \sum_{j' \neq j}^{L} \frac{j^2}{L^2} = (L-1) \frac{L(L+1)(2L+1)}{6L^2}, \qquad (\text{Equation 45}) $$
while the last term is calculated as

$$ 
\begin{aligned}
& \frac{1}{L^2} \sum_{j=1}^{L} j \sum_{j' \neq j}^{L} j' \\
& = \frac{1}{L^2} \sum_{j=1}^{L} j \left( \frac{L(L-1)}{2} - j \right) \\
& = \frac{L^2(L-1)^2}{4L^2} - \frac{L(L+1)(2L+1)}{6L^2} \\
& = \frac{3L^3 - 10L^2 - 3L - 2}{12L}.
\end{aligned}
\quad (\text{Equation 46})
$$

Finally, combining Equations 45 and 46 into Equation 44, we obtain

$$ 
\begin{aligned}
\bar{\epsilon}_{n,g}^2 &= \frac{(L+1)(2L+1)}{3L^2} - \frac{3L^3 - 10L^2 - 3L - 2}{6L^2(L-1)} \\
&= \frac{1}{6} \left( 1 + \frac{12}{L-1} + \frac{1}{L} \right).
\end{aligned}
\quad (\text{Equation 47})
$$

This is a term of order 1, the size of the stimulus range, plus corrections of order 1/L.

## Broad tuning curves

In the case of broad tuning curves, we consider the regime of smooth response curves on the scale of the noise amplitude, such that the mean population activity can be approximated locally by a linear function of the stimulus. This regime obtains when the second-order term in the Taylor expansion is negligible with respect to the first-order one:

$$ \frac{1}{4} \| \mathbf{v}''(x) \|_2^2 \Delta x^4 \ll \| \mathbf{v}'(x) \|_2^2 \Delta x^2, \quad (\text{Equation 48}) $$

where $\Delta x^2 \approx \eta^2$. In order to express this condition in terms of model parameters, we impose it on average over network realizations; we note that this leads to the same result as imposing the condition on average over stimuli, but it requires a simpler calculation. From the identity $\langle W_{ij}W_{ij'} \rangle_W = \frac{1}{L}\delta_{ij'}$, the average of the left-hand-side of Equation 48 is obtained as

$$ \langle \| \mathbf{v}''(x) \|_2^2 \rangle_W = \langle \sum_{i=1}^N \sum_{j=1}^L W_{ij} W_{ij'} u_j''(x) u_j''(x) \rangle_W = \sum_{i=1}^N \sum_{j=1}^L \frac{1}{L} u_j''(x)^2 \approx \frac{3\sqrt{\pi N A^2}}{4\sigma^3}, \quad (\text{Equation 49}) $$

where the approximations consists in replacing the sum with an integral and in extending the integration domain to the real line. A similar calculation can be performed for the right-hand-side of Equation 48:

$$ \langle \| \mathbf{v}'(x) \|_2^2 \rangle_W = \sum_{i=1}^N \sum_{j=1}^L \frac{1}{L} u_j'(x)^2 \approx \frac{\sqrt{\pi N A^2}}{2\sigma}. \quad (\text{Equation 50}) $$

By combining Equations 50 and 49, and substituting $\Delta x^2$ by the variance of the noise, $\eta^2$, in Equation 48, we obtain the smoothness condition as

$$ \frac{3\eta^2}{8\sigma^2} \ll 1. \quad (\text{Equation 51}) $$

In the case of broad tuning curves the error can be of two qualitatively different types: local or global (Figure 3A). The width of the Gaussian correlation function in Equation 16 gives a measure of the distance in the stimulus space at which population responses are correlated; we refer to a global error when the distance between the stimulus and its estimate is greater than this ‘correlation length’, $\sigma$. We write the MSE as

$$ \epsilon^2 = \epsilon_f^2 + \epsilon_g^2, \quad (\text{Equation 52}) $$

and we compute these two terms.

## Local error

According to the ML decoder, Equation 34, the stimulus estimate corresponds to the value $x'$ that minimizes the distance between $\mathbf{v}(x')$ and $\mathbf{r}$; if the error is local, this is obtained by projecting the noise vector onto the curve of mean population activity. By expanding the response curve around $\mathbf{v}(x)$, we obtain, to linear order,

$$ \| \mathbf{z} \cdot \hat{\mathbf{v}}'(x) \|_2^2 \approx \| \mathbf{v}(x + \Delta x) - \mathbf{v}(x) \|_2^2 \approx \| \mathbf{v}'(x) \|_2^2 \Delta x^2, \quad (\text{Equation 53}) $$

where $\hat{\mathbf{v}}'(x) = \mathbf{v}'(x) / \|\mathbf{v}'(x)\|_2$. The local error can then be calculated as $\Delta x^2 = (\hat{x} - x)^2 \approx \|\mathbf{z} \cdot \hat{\mathbf{v}}'(x)\|_2^2 / \|\mathbf{v}'(x)\|_2^2$. By averaging over the noise and the synaptic weights, we obtain the mean local error as

$$ \epsilon_f^2 = \langle \int_0^1 dx \int dz \rho(\mathbf{z}) \frac{\| \mathbf{z} \cdot \hat{\mathbf{v}}'(x) \|_2^2}{\| \mathbf{v}'(x) \|_2^2} \rangle_W = \langle \int_0^1 dx \frac{\eta^2}{\| \mathbf{v}'(x) \|_2^2} \rangle_W. \quad (\text{Equation 54}) $$
The squared norm of the derivative of the tuning curves is the realization of the random variable

$$ \| \mathbf{v}'(x) \|_2^2 = \sum_{i=1}^{N} \left( \sum_{j=1}^{L} W_{ij} u_j'(x) \right)^2. \qquad (\text{Equation 55}) $$

The terms of the inner sum, $W_{ij}u_j'(x)$, are realizations of independent Gaussian random variables with variable variance; as a result, the outer sum is also the realization of a Gaussian random variable with mean equal to the sum of the means of its terms and variance equal to the sum of the variances. The sum of the variances can be calculated as

$$ \sum_{j=1}^{L} \frac{u_j'(x)^2}{L} \approx \int_{-\infty}^{\infty} dc_j u_j'(x)^2 = \frac{\sqrt{\pi} A^2}{2\sigma}, \qquad (\text{Equation 56}) $$

where the approximation consists in replacing the sum with an integral and in extending the integration domain to the real line. Therefore, the inner sum is distributed according to

$$ \bar{W}_i \equiv \sum_{j=1}^{L} W_{ij} u_j'(x) \sim \mathcal{N} \left( 0, \frac{\sqrt{\pi} A^2}{2\sigma} \right). \qquad (\text{Equation 57}) $$

As a result, the quantity $1/\|\mathbf{v}'(x)\|_2^2 = 1/\sum_{i=1}^N \bar{W}_i^2$ is sampled according to a scaled inverse chi-squared distribution, with mean given by

$$ \left( \frac{1}{\sum_{i=1}^{N} \bar{W}_i^2} \right)_W = \frac{2\sigma}{\sqrt{\pi(N-2)A^2}} \approx \frac{2\sigma}{\sqrt{\pi NA^2}}. \qquad (\text{Equation 58}) $$

The local error is then obtained, from Equations 54 and 58, as

$$ \varepsilon_i^2 = \frac{2\sigma\eta^2}{\sqrt{\pi NA^2}}. \qquad (\text{Equation 59}) $$

Finally, if we approximate the response amplitude for small $\sigma$ by $A^2 \approx R/\sqrt{\pi\sigma^2}$, we obtain the expression that appears in the main text (first term of Equation 6). We note that this expression is equal to the inverse of the Fisher information averaged over network realizations; the Fisher information in case of neural responses corrupted by independent Gaussian noise is given by

$$ \langle J(x) \rangle_W = \left\langle \frac{\| \mathbf{v}'(x) \|_2^2}{\eta^2} \right\rangle_W = \frac{\sqrt{\pi} N A^2}{2\sigma \eta^2}. \qquad (\text{Equation 60}) $$

## Global error

Here, we extend the calculation performed in the case of discrete stimuli. The analog of Equation 36 in the case of broad tuning curves is

$$ \varepsilon_g^2 = \langle P_b(E) \rangle_W \bar{\varepsilon}_g^2, \qquad (\text{Equation 61}) $$

where the two factors are the probability of a global error and the average squared magnitude of a global error, respectively.

We approximate the probability of a global error by considering a division of the curve of mean population activity into $\sigma$ 'segments'. These segments are roughly uncorrelated and appear in random locations in the space of population activity; as a result, we can replace $L$ by the number of segments in Equation 5 to obtain the probability of a global error, as

$$ \langle P_b(E) \rangle_W \approx \frac{1}{\sigma \sqrt{2\pi N}} \exp \left( -\log \left( 1 + \frac{R}{2\eta^2} \right) \frac{N}{2} \right). \qquad (\text{Equation 62}) $$

Similarly to the discrete case, when a global error occurs, the decoded stimulus is uniformly sampled from all the other stimuli belonging to incorrect segments. We illustrate the calculation for the case $x - \sigma > 0$ and $x + \sigma < 1$; similar results can be obtained for stimuli close to the boundaries of the stimulus range. In this case, the output of the decoder is distributed uniformly in the interval $\hat{x} \in [0, x - \sigma] \cup [x + \sigma, 1]$. The average magnitude of global errors is therefore

$$ \bar{\varepsilon}_g^2 = \left\langle \int_0^1 dx (\hat{x} - x)^2 \right\rangle_W \approx \int_0^1 dx \frac{1}{(1 - 2\sigma)} \left[ \int_0^{x-\sigma} d\hat{x} (\hat{x} - x)^2 + \int_{x+\sigma}^1 d\hat{x} (\hat{x} - x)^2 \right] = \frac{1(1 - 4\sigma^3)}{6(1 - 2\sigma)}, \qquad (\text{Equation 63}) $$

which is a term of the order 1, the size of the stimulus range, plus corrections of order $\sigma$. We obtain the expression for the global error which appears in the main text (second term of Equation 6), by combining Equations 61, 62, and 63.
## Local and global errors in the case of multi-dimensional stimuli

The MSE for multi-dimensional stimuli, Equation 30, averaged over synaptic weights realizations, is defined as the sum of the MSEs along each stimulus dimension,

$$ \epsilon^2 = \sum_{k=1}^{K} \epsilon_k^2 = \sum_{k=1}^{K} \left\langle \int_{0}^{1} dx_k \int d\mathbf{r} \rho(\mathbf{r}|\mathbf{x}) (\hat{\mathbf{x}}_k - x_k)^2 \right\rangle_W \quad (\text{Equation 64}) $$

The local error along stimulus dimension k can be calculated, similarly to Equation 54, as

$$ \epsilon_{i,k}^2 = \left\langle \int d\mathbf{z} \rho(\mathbf{z}) (\hat{\mathbf{x}}_k - x_k)^2 \right\rangle_W \approx \left\langle \int d\mathbf{z} \rho(\mathbf{z}) \frac{\|\mathbf{z} \cdot \hat{\mathbf{v}}'_k(\mathbf{x})\|_2^2}{\|\hat{\mathbf{v}}'_k(\mathbf{x})\|_2^2} \right\rangle_W \approx \left\langle \frac{\eta^2}{\|\hat{\mathbf{v}}'_k(\mathbf{x})\|_2^2} \right\rangle_W \quad (\text{Equation 65}) $$

where the noise is projected onto the direction parallel to the partial derivative of the mean activity with respect to stimulus dimension $k$, $\mathbf{v}'_k(\mathbf{x}) = \partial\mathbf{v}(\mathbf{x})/\partial x_k$.

Local error—sensory neurons with pure tuning. The derivative of the tuning function with respect to stimulus dimension k is given by

$$ \mathbf{v}'_{i,k}(\mathbf{x}) = \frac{\partial v_i(\mathbf{x})}{\partial x_k} = \sum_{j=1}^{Q} W_{ijk} \frac{\partial u_{j,k}^p(x_k)}{\partial x_k} \quad (\text{Equation 66}) $$

Similarly to the one dimensional case, this is a sum of realizations of independent Gaussian random variables. Dropping the superscript p for the sake of clarity, the sum of the variances of these terms is calculated as

$$ \sum_{j=1}^{Q} \frac{1}{L} \left( \frac{\partial u_{j,k}(x_k)}{\partial x_k} \right)^2 \approx \frac{1}{K} \int_{-\infty}^{\infty} dc_j^k \left( \frac{\partial u_{j,k}(x_k)}{\partial x_k} \right)^2 = \frac{\sqrt{\pi A^2}}{2K\sigma}, \quad (\text{Equation 67}) $$

where the approximation consists in replacing the sum $\sum_{j=1}^{Q} \frac{1}{L} f(c_j^k)$ with an integral $\int dc_j^k f(c_j^k)$, and in extending the integration domain to the real line. The sum is distributed as

$$ \bar{W}_{i,k}^p \equiv \sum_{j=1}^{Q} W_{ijk} \frac{\partial u_{j,k}(x_k)}{\partial x_k} \sim \mathcal{N} \left( 0, \frac{\sqrt{\pi A_p^2}}{2K\sigma} \right) \quad (\text{Equation 68}) $$

Finally, by calculating the mean of the scaled inverse chi-squared distribution,

$$ \left\langle \frac{1}{\sum_{i}^{N} (\bar{W}_{i,k}^p)^2} \right\rangle_W \approx \frac{2K\sigma}{\sqrt{\pi N A_p^2}}, \quad (\text{Equation 69}) $$

in Equation 65, we obtain the local error along a single stimulus dimension, as

$$ \epsilon_{p,j,k}^2 = \frac{2K\sigma\eta^2}{\sqrt{\pi N A_p^2}}, \quad (\text{Equation 70}) $$

the total local error is then obtained by summing over dimensions,

$$ \epsilon_{p,j,k}^2 = \sum_{k=1}^{K} \epsilon_{p,j,k}^2 = \frac{2K^2 \sigma \eta^2}{\sqrt{\pi N A_p^2}} \quad (\text{Equation 71}) $$

Local error—sensory neurons with conjunctive tuning. The derivative of the tuning function with respect to stimulus dimension k is given by

$$ \begin{aligned} \mathbf{v}'_{i,k}(\mathbf{x}) &= \frac{\partial v_i(\mathbf{x})}{\partial x_k} = \sum_{j=1}^{L} W_{ij} \frac{\partial u_j^c(\mathbf{x})}{\partial x_k} \\ &= -\sum_{j=1}^{L} W_{ij} \frac{(x_k - c_{jk})}{\sigma^2} u_j^c(\mathbf{x}), \end{aligned} \quad (\text{Equation 72}) $$

where $c_{jk}$ is the kth component of the preferred stimulus of neuron j, $c_j$. Similarly to the previous calculations, this is a sum of realizations of independent Gaussian random variables of different variances. Dropping the superscript c for the sake of clarity, the sum of the variances of these terms is calculated as
$$ \sum_{j=1}^{L} \frac{(x_k - c_{j,k})^2}{L \sigma^4} u_j(\mathbf{x})^2 \approx \int d\mathbf{c}_j \frac{(x_k - c_{j,k})^2}{\sigma^4} u_j(\mathbf{x})^2 \quad \text{(Equation 73)} \\ \approx \frac{\pi^{K/2} A_c^2}{2\sigma^{(2-K)}} $$

where the approximation consists in replacing the sum $\sum_{j=1}^{L} \frac{1}{L} f(\mathbf{c}_j)$ with a K-dimensional integral $\int d\mathbf{c}_j f(\mathbf{c}_j)$, and in extending the integration domain. The sum is therefore distributed as

$$ \overline{W}_{j,k}^{c} = \sum_{j=1}^{L} w_{ij} \frac{\partial u_j(\mathbf{x})}{\partial x_k} \sim \mathcal{N} \left( 0, \frac{\pi^{(K/2)} A_c^2}{2\sigma^{(2-K)}} \right). \quad \text{(Equation 74)} $$

Finally, by calculating the mean of the scaled inverse chi-squared distribution,

$$ \left\langle \frac{1}{\sum_{i}^{N} (\overline{W}_{i,k}^{c})^2} \right\rangle_W \approx \frac{2\sigma^{2-K}}{\pi^{(K-2)}NA_c^2}, \quad \text{(Equation 75)} $$

in Equation 65, and by summing over dimensions, we obtain the total local error as

$$ \epsilon_{c,j}^2 = \sum_{k=1}^{K} \epsilon_{c,j,k}^2 = \frac{2\sigma^{(2-K)}\eta^2}{\pi^{K/2}NA_c^2}. \quad \text{(Equation 76)} $$

If we approximates $A_c^2$ and $A_p^2$ for small values of $\sigma$, we obtain that the ratio of the local errors in case of sensory neurons with pure and conjunctive tuning is

$$ \frac{\epsilon_{c,j}^2}{\epsilon_{p,j}^2} \approx \frac{1}{K}. \quad \text{(Equation 77)} $$

*Global error—sensory neurons with pure tuning.* In the case of sensory neurons with pure tuning, the tuning function of a representation neuron is obtained as the superposition of one-dimensional tuning curves (Equation 22). According to the ML decoder, Equation 34, the decoder output can be written as

$$ \begin{aligned} \hat{\mathbf{x}} &= \arg \min_{\mathbf{x}'} \| \mathbf{v}(\mathbf{x}) + \mathbf{z} - \mathbf{v}(\mathbf{x}') \|_2^2 \\ &= \arg \min_{\mathbf{x}'} \left\| \sum_{k=1}^{K} (\mathbf{v}_k(\mathbf{x}_k) - \mathbf{v}_k(\mathbf{x}'_k) + \mathbf{z}_k) \right\|_2^2, \end{aligned} \quad \text{(Equation 78)} $$

where $\mathbf{z}_k$ is the projection of the noise vector onto the direction parallel to the partial derivative of the mean activity with respect to stimulus dimension $k$. For most realizations of the random tuning curves, if $K \ll N$, the K vectors summed in Equation 78 are likely orthogonal. Thus, minimizing the squared norm of the sum is equivalent to minimizing the sum of the squared norms of each of the vectors. This, in turn, the stimulus estimate can be obtained independently for each stimulus dimension, as

$$ \hat{\mathbf{x}}_k = \arg \min_{\mathbf{x}'_k} \| \mathbf{v}_k(\mathbf{x}_k) - \mathbf{v}_k(\mathbf{x}'_k) + \mathbf{z}_k \|_2^2. \quad \text{(Equation 79)} $$

Therefore, a global error can occur in one or several stimulus dimensions; it requires that $|\hat{\mathbf{x}}_k - \mathbf{x}_k| > \sigma$ for some $k$. If the probability of a global error on more than one stimulus dimension is negligible, the total probability of a global error can be approximated as the sum of probabilities over dimensions, $\langle P(E_{p,g}) \rangle_W \approx \sum_{k=1}^{K} \langle P(E_{k,g}) \rangle_W$. We calculated the probability of a global error in the one-dimensional case in the previous section. In order to extend the formula to this case, we have to take into account that the variance of the tuning function along one stimulus dimension is

$$ \langle \int_0^1 dx_k [v_i(\mathbf{x}) - (\int_0^1 dx_k v_i(\mathbf{x}))] \rangle_W = \sum_{j=1}^Q \frac{1}{L} \left( \int_0^1 dx_k u_{j,k}(\mathbf{x}_k)^2 - \left( \int_0^1 dx_k u_{j,k}(\mathbf{x}_k) \right)^2 \right) \approx \frac{R}{K}. \quad \text{(Equation 80)} $$

This quantity is the signal variance which governs the rate of exponential suppression of the probability of global error; replacing $R$ by $R/K$ in Equation 62, multiplying by the average squared magnitude of global errors and summing over dimensions, we obtain the global error as

$$ \epsilon_{p,g}^2 \approx \frac{K \bar{\epsilon}_g^2}{\sigma \sqrt{2\pi N}} \exp \left( - \log \left( 1 + \frac{R}{2K\eta^2} \right) \frac{N}{2} \right). \quad \text{(Equation 81)} $$
Global error—sensory neurons with conjunctive tuning. The correlation of the responses of neuron i to two stimuli, **x** and **x'**, reads

$$ \langle v_i(\mathbf{x})v_i(\mathbf{x}') \rangle_W \approx A_c^2 (\pi\sigma^2)^{K/2} \exp\left(-\frac{\Delta\mathbf{x}^2}{4\sigma^2}\right), \qquad (\text{Equation 82}) $$

where $\Delta\mathbf{x}^2 = \|\mathbf{x} - \mathbf{x}'\|_2^2$; it is exponentially suppressed if $\|\mathbf{x} - \mathbf{x}'\|_2 > \sigma$. By analogy to the one-dimensional case, we divide the surface described by the population activity as a function of the stimulus, $\mathbf{v}(\mathbf{x}) = \{v_1(\mathbf{x}), \dots, v_N(\mathbf{x})\}$, into $1/\sigma^K$ uncorrelated regions. We calculate the global error by replacing $L$ with the number of uncorrelated regions in Equation 5, obtaining

$$ \epsilon_{c,g}^2 \approx \frac{\epsilon_{c,g}^2}{\sigma^K \sqrt{2\pi N}} \exp\left(-\log\left(1 + \frac{R}{2\eta^2}\right) \frac{N}{2}\right), \qquad (\text{Equation 83}) $$

where $\epsilon_{c,g}^2$ is the average squared magnitude of a global error, a term of the order of the stimulus range.

## Influence of correlated output noise on population coding

### Correlated output noise due to independent noise in sensory neurons

We consider the case in which the activity of sensory neurons is affected by independent Gaussian noise: $\tilde{u}(x) = u(x) + z^u$, with $z_i^u \sim N(0, \xi^2)$. This results in a multivariate Gaussian noise in the responses of representation neurons, with covariance matrix $\Sigma = \eta^2 I + \xi^2 WW^T$. The matrix $WW^T$ is sampled according to a Wishart distribution, with mean $I$ and variance of the matrix elements of order $1/L^{85}$. We write the covariance matrix as the identity plus a perturbation, as

$$ \begin{aligned} \Sigma &= \tilde{\eta}^2 I + \xi^2 (WW^T - I) \\ &= \tilde{\eta}^2 \left(I + \frac{\xi^2}{\tilde{\eta}^2} (WW^T - I)\right), \end{aligned} \qquad (\text{Equation 84}) $$

where $\tilde{\eta}^2 = \eta^2 + \xi^2$. In order to quantify the effect of input noise on the coding performance, we calculate the inverse of the Fisher information (FI) as a lower bound to the MSE. The FI is written as

$$ \begin{aligned} J(x) &= v'(x)^T \Sigma^{-1} v'(x) \\ &= u'(x)^T W^{-1} u(x). \end{aligned} \qquad (\text{Equation 85}) $$

We expand the inverse of the noise covariance matrix to second order in $\xi^2/\tilde{\eta}^2$, as

$$ \Sigma^{-1} \approx \frac{1}{\tilde{\eta}^2} \left( I - \frac{\xi^2}{\tilde{\eta}^2} (WW^T - I) + \frac{\xi^2}{\tilde{\eta}^4} (WW^T - I)^2 \right). \qquad (\text{Equation 86}) $$

In this approximation, the FI becomes $J(x) \approx J_{\text{ind}}(x) + \delta J(x)$, with

$$ J_{\text{ind}}(x) = \frac{1}{\tilde{\eta}^2} u'(x)^T B u'(x), \qquad (\text{Equation 87}) $$

and

$$ \delta J(x) = \frac{1}{\tilde{\eta}^2} \left( - \frac{\xi^2}{\tilde{\eta}^2} u'(x)^T (\mathbf{B}^2 - \mathbf{B}) u'(x) + \frac{\xi^4}{\tilde{\eta}^4} u'(x)^T (\mathbf{B}^3 - 2\mathbf{B}^2 + \mathbf{B}) u'(x) \right), \qquad (\text{Equation 88}) $$

where $\mathbf{B} = WW^T$. The first term, $J_{\text{ind}}(x)$, is the FI for independent Gaussian output noise with variance $\tilde{\eta}^2$; by averaging over synaptic weights realizations, we obtain the expression in Equation 60,

$$ \langle J_{\text{ind}}(x) \rangle_W = \frac{\sqrt{\pi N A^2}}{2\sigma\tilde{\eta}^2}. \qquad (\text{Equation 89}) $$

The average of the second term, $\delta J(x)$, over network realizations depends on the moments of the matrix $\mathbf{B}$, which can be computed using Wick's theorem: from the identity $\langle W_{ij}W_{mn} \rangle_W = \frac{1}{L}\delta_{im}\delta_{jn}$, we obtain

$$ \langle B_{mn} \rangle_W = \left\langle \sum_{j=1}^{N} W_{jm} W_{jn} \right\rangle_W = \frac{N}{L} \delta_{mn}, \qquad (\text{Equation 90}) $$

$$ \langle B_{mn}^2 \rangle_W = \left\langle \sum_{i=1}^{L} \sum_{j=1}^{N} W_{im} W_{j'j} W_{i'j} W_{j'm} \right\rangle_W = \left( \frac{N}{L} + \frac{N^2}{L^2} + \frac{N}{L^2} \right) \delta_{mn}, \qquad (\text{Equation 91}) $$
$$ \langle B_{mn}^3 \rangle_W = \left\langle \sum_{i=1}^{L} \sum_{j'=1}^{N} W_{ji} W_{j'i} W_{ji'} W_{j'i'} W_{j'n} W_{j'n'} \right\rangle_W = \left( \frac{N}{L} + 3 \frac{N^2}{L^2} + 3 \frac{N}{L^2} + \frac{N^3}{L^3} + 3 \frac{N^2}{L^3} + 4 \frac{N}{L^3} \right) \delta_{mn}. \quad (\text{Equation 92}) $$

From now on, we consider the terms up to $O(N^2/L^2)$; the mean of the perturbation term in the FI becomes

$$ \langle \delta J(x) \rangle_W \approx \frac{1}{\tilde{\eta}^2} \mathbf{u}'(x)^T \mathbf{u}'(x) \left( - \frac{N^2 \xi^2}{L^2 \tilde{\eta}^2} + \frac{N^2 \xi^4}{L^2 \tilde{\eta}^4} \right). \quad (\text{Equation 93}) $$

Finally, we compute the first factor by approximating the discrete sum with the integral, similarly to previous calculations, obtaining

$$ \begin{aligned} \frac{1}{\tilde{\eta}^2} \mathbf{u}'(x)^T \mathbf{u}'(x) &= \frac{1}{\tilde{\eta}^2} \sum_{j=1}^{L} \frac{(x - c_j)^2}{\sigma^4} u_j(x)^2 \\ &\approx \frac{L}{\sigma^4 \tilde{\eta}^2} \int_{-\infty}^{\infty} dc_j (x - c_j)^2 u_j(x)^2 \\ &= \frac{\sqrt{\pi L A^2}}{2 \sigma \tilde{\eta}^2}. \end{aligned} \quad (\text{Equation 94}) $$

This quantity is proportional to the mean of the FI in the case of independent noise, Equation 89, by a factor $N/L$. Combining Equations 89, 93, and 94, we obtain

$$ \langle J(x) \rangle_W \approx \frac{\sqrt{\pi N A^2}}{2 \sigma \tilde{\eta}^2} \left( 1 - \frac{N \xi^2}{L \tilde{\eta}^2} + \frac{N \xi^4}{L \tilde{\eta}^4} \right). \quad (\text{Equation 95}) $$

We approximate the local error as the inverse of the FI; including only corrections up to $O(N^4/L^4)$, we obtain

$$ \epsilon_f^2 \approx \frac{1}{\langle J(x) \rangle_W} \approx \epsilon_{f, \text{ind}}^2 \left( 1 + \frac{N \xi^2}{L \tilde{\eta}^2} - \frac{N \xi^4}{L \tilde{\eta}^4} \right). \quad (\text{Equation 96}) $$

## Correlated output noise with random covariance structure

Similar calculations can be carried out for a noise covariance matrix that obeys the same statistics as those of $WW^T$, but that does not derive from the structure of synaptic weights. We consider

$$ \Sigma_{\text{rand}} = \eta^2 I + \xi^2 \mathbf{X}\mathbf{X}^T, \quad (\text{Equation 97}) $$

with $X_{ij} \sim N(0, I)$, such that $\langle X_{ij}W_{mn} \rangle_{W,X} = 0$ and $\langle X_{ji}X_{mn} \rangle_{X} = \frac{1}{L}\delta_{im}\delta_{jn}$. In this case, by expanding the inverse of the covariance matrix to second order in $\xi^2/\tilde{\eta}^2$ in Equation 85, we obtain the perturbation term in the FI as

$$ \delta J(x) = \frac{1}{\tilde{\eta}^2} \left( - \frac{\xi^2}{\tilde{\eta}^2} \mathbf{u}'(x)^T (\mathbf{W}^T \mathbf{X}\mathbf{X}^T \mathbf{W} - \mathbf{B}) \mathbf{u}'(x) + \frac{\xi^4}{\tilde{\eta}^4} \mathbf{u}'(x)^T (\mathbf{W}^T (\mathbf{X}\mathbf{X}^T)^2 \mathbf{W} - 2\mathbf{W}^T \mathbf{X}\mathbf{X}^T \mathbf{W} + \mathbf{B}) \mathbf{u}'(x) \right). \quad (\text{Equation 98}) $$

We compute the mean of these matrices over realizations of the noise covariance matrix and of the synaptic matrix using Wick's theorem. We obtain

$$ \langle (W^T XX^T W)_{mn} \rangle_{W,X} = \left\langle \sum_{i=1}^{L} \sum_{j'=1}^{N} X_{ji} W_{jm} X_{j'i} W_{j'n} \right\rangle_{W,X} = \frac{N}{L} \delta_{mn}, \quad (\text{Equation 99}) $$

$$ \langle (W^T (XX^T)^2 W)_{mn} \rangle_{W,X} = \left\langle \sum_{i=1}^{L} \sum_{j'=1}^{N} W_{ji} X_{ji} X_{j'i} X_{j'n} W_{j'n} \right\rangle_{W,X} = \left( \frac{N}{L} + \frac{N^2}{L^2} \right) \delta_{mn}. \quad (\text{Equation 100}) $$

Therefore, the first order correction vanishes, and the FI is increased,

$$ \langle J(x) \rangle_W x \approx \frac{\sqrt{\pi N A^2}}{2 \sigma \tilde{\eta}^2} \left( 1 + \frac{N \xi^4}{L \tilde{\eta}^4} \right), \quad (\text{Equation 101}) $$

yielding a negative correction to the MSE (Equation 125).
Influence of recurrent connectivity on population coding

We define the dynamic of the synaptic activity (see also Equation 127) as

$$ \tau \frac{ds_i}{dt} = -s_i + \sum_j M_{ij}s_j + \sum_j W_{ij}u_j(x) + \zeta_i(t), \qquad (\text{Equation 102}) $$

where we assume that the noise is uncorrelated in time and among different neurons, so that

$$ \langle \zeta_i(t) \zeta_j(t') \rangle = \Delta \delta_{ij} \delta(t-t'), \qquad (\text{Equation 103}) $$

where $\Delta$ is a constant. We consider a readout scheme based on a weighted averaged of the synaptic activity over a time scale $T$. For the sake of simplicity, we assume that this is done by convolving $s_i$ with an exponential kernel, to obtain an integrated variable $r_i$ according to

$$ T \frac{dr_i}{dt} + r_i = s_i, \qquad (\text{Equation 104}) $$

In the steady state the mean of $r$ obeys (in vector notation)

$$ \langle \mathbf{r} \rangle = \langle \mathbf{s} \rangle = (\mathbf{I} - \mathbf{M})^{-1} \mathbf{W} \mathbf{u}(\mathbf{x}). \qquad (\text{Equation 105}) $$

We next define the difference of $r_i$ from its mean through

$$ r_i = \langle r_i \rangle + \xi_i(t) \qquad (\text{Equation 106}) $$

and our goal is to evaluate the equal-time correlation matrix of $\xi_i$,

$$ \Sigma_{ij} = \langle \xi_i(t) \xi_j(t') \rangle. \qquad (\text{Equation 107}) $$

Using Equations 102 and 104, we can express $\xi_i(t)$ in terms of $\zeta_i(t)$: in Fourier space, this relation takes the form

$$ \tilde{\xi}_i(\omega) = \sum_k \Gamma_{ik}(\omega) \tilde{\zeta}_k(\omega), \qquad (\text{Equation 108}) $$

where $\tilde{f}(\omega)$ denotes the Fourier transform of $f(t)$, and the functions $\Gamma_{ik}(\omega)$ are computed below. By using the fact that

$$ \langle \tilde{\xi}_i(\omega) \tilde{\xi}_j(\omega') \rangle = 2\pi\Delta\delta_{ij}\delta(\omega + \omega'), \qquad (\text{Equation 109}) $$

we obtain the equal-time correlation as

$$ \Sigma_{ij} = \frac{\Delta}{2\pi} \int_{-\infty}^{\infty} \sum_k \Gamma_{ik}(\omega) \Gamma_{jk}(-\omega) d\omega. \qquad (\text{Equation 110}) $$

To evaluate $\Gamma_{ij}(\omega)$, we first reexpress Equation 102 in Fourier space and we isolate the term $\tilde{s}(\omega)$. In vector notation, we obtain for any $\omega \neq 0$

$$ \tilde{s}(\omega) = (i\omega\tau + \mathbf{I} - \mathbf{M})^{-1} \tilde{\zeta}(\omega). \qquad (\text{Equation 111}) $$

Next, by similarly reexpressing Equation 104 in Fourier space and substituting Equation 111 for $\tilde{s}$, we obtain

$$ \tilde{r}(\omega) = \frac{1}{i\omega T + 1} (i\omega\tau + \mathbf{I} - \mathbf{M})^{-1} \tilde{\zeta}(\omega). \qquad (\text{Equation 112}) $$

As a result,

$$ \Gamma(\omega) = \frac{1}{i\omega T + 1} (i\omega\tau + \mathbf{I} - \mathbf{M})^{-1} \qquad (\text{Equation 113}) $$

and, using Equation 110 we can write the covariance matrix of the noise as

$$ \Sigma = \frac{\Delta}{2\pi} \int_{-\infty}^{\infty} \frac{1}{\omega^2 T^2 + 1} (i\omega\tau + \mathbf{I} - \mathbf{M})^{-1} (-i\omega\tau + \mathbf{I} - \mathbf{M})^{-T} d\omega. \qquad (\text{Equation 114}) $$

The behavior for $T \gg \tau$ can be identified by rewriting this integral in terms of $\Omega = \omega T$,

$$ \Sigma = \frac{\Delta}{2\pi T} \int_{-\infty}^{\infty} \frac{1}{\Omega^2 + 1} \left( \frac{i\Omega\tau}{T} + \mathbf{I} - \mathbf{M} \right)^{-1} \left( -\frac{i\Omega\tau}{T} + \mathbf{I} - \mathbf{M} \right)^{-T} d\Omega, \qquad (\text{Equation 115}) $$
from which we obtain, to leading order in $\tau/T$,

$$ \Sigma \approx \frac{\Delta}{2\pi T} \mathbf{H}^{-1} \mathbf{H}^{-T}, \qquad (\text{Equation 116}) $$

where $\mathbf{H} = \mathbf{I} - \mathbf{M}$. We note, first, that the covariance matrix of the noise is proportional to $1/T$, as expected. Second, the structure of the covariance matrix is of the form obtained by distorting the uncorrelated noise among neurons through the transformation $\mathbf{H}^{-1}$.

## Data analysis and model fitting
### Description of the data and summary statistics

The data consist of the responses (firing rates) of $N \sim 500$ neurons, recorded during an arm posture ‘hold’ task including 27 different positions, with 2 hand orientations each, arranged in a virtual cube of size $40 \times 40 \times 40$ cm. The response of each neuron for each hand position is recorded in several trials (~10 trials per hand position). Tuning curves are computed by averaging over trials. In order to quantify the degree of irregularity of a tuning curve in a non-parametric form, the authors used a ‘complexity measure’: for neuron $i$, it is defined as the standard deviation of a discretized derivative of the mean response:

$$ c(D_{\min})_i = \text{std} \left( \frac{\| v_i(\mathbf{x}) - v_i(\mathbf{x} + \Delta\mathbf{x}) \|}{\sqrt{\| \Delta\mathbf{x} \|^2}} \quad \text{s.t.} \quad \| \Delta\mathbf{x} \|_2^2 < D_{\min} \right), \qquad (\text{Equation 117}) $$

where $v_i(\mathbf{x})$ is the mean response, $D_{\min}$ is the distance between two neighboring hand positions, and in the experiment is equal to 35.33 evaluated also another summary statistics, the distribution of $R^2$ values resulting from a fit of the tuning curves with a linear model (see Equation 9, originally proposed by$^{48}$):

$$ R_i^2 = 1 - \frac{\sum_x (v_i(\mathbf{x}) - v(\mathbf{x}))^2}{\sum_x v(\mathbf{x})^2}, \qquad (\text{Equation 118}) $$

where $v_i(\mathbf{x})$ is the response predicted by a linear regression of the data, and the sum is over hand positions used in the experiment. The distribution of these two quantities across neurons is a measure of the irregularity of the neural population response; if the population were perfectly described by a linear model, the $R^2$-distribution would be a constant for all neurons and equal to 1, while the complexity measure would exhibit low values.

### Model fitting and comparison between irregular and linear tuning curves

We consider neurons responding with at least 5 spikes/s at more than two target positions and we compute their tuning curves by averaging the firing rates over trials. Then, we shift and normalize the tuning curves to cancel their means and set their variances across hand positions to unity. We use a version of our shallow network model to produce three-dimensional mean-response profiles. The sensory layer is made up of $L = 100^3$ neurons; the preferred stimuli (here, hand positions) are arranged so as to cover a space of $100 \times 100 \times 100$ cm, in such a way that hand positions used in the experiment are placed far from the boundaries of the stimulus space. To limit computation load, we choose $\mathbf{W}$ as a sparse random matrix, with sparsity equal to 0.1, with Gaussian-distributed elements, similarly to the model of.33 The sparsity of the matrix does not affect our results, as long a proper normalization of the synaptic weights is taken into account and the representation neurons receive a sufficient number of inputs from the sensory layer, i.e., as long as the matrix is not too sparse and the tuning width is not too narrow. The tuning curves are normalized to have zero mean and unit variance across hand positions. With respect to the model of,33 there are two main differences: in their case the random weights were distributed according to a uniform distribution, and a rectifying non-linear function was applied to the random sum of the activity of first-layer neurons to enforce a positive activity of the representation neurons. Their model thus had two tunable parameters: the tuning width of first-layer neurons, $\sigma$, and the threshold of the non-linear response function in the second layer. The only tunable parameter in our model is $\sigma$.

In order to fit our model, we generate neural responses of a number of representation neurons equal to the number of recorded neurons, using the same set of hand positions to as used in the experiment. We then computed the distribution of the complexity measure for different values of $\sigma$; we denote by $\sigma_f$ the tuning-curve width which minimizes the Kolmogorov-Smirnov (KS) distance between the distribution produced by the model and that extracted from the data (Figure S1A). The KS distance is a measure of discrepancy between two probability distributions. We denote by $F_{\text{data/model}}(c)$ the empirical cumulative distribution function of the complexity measure across data/model, that is, the empirical probability of finding a neuron with complexity measure less than c,

$$ F_{\text{data/model}}(c) = \frac{\# \text{ neurons in the data/model with complexity measure} < c}{N}, \qquad (\text{Equation 119}) $$

where N is the total number of neurons. The KS distance is defined as the maximum absolute difference between the $F_{\text{data}}$ and $F_{\text{model}}$:

$$ KS \equiv \max_{c} |F_{\text{data}}(c) - F_{\text{model}}(c)|. \qquad (\text{Equation 120}) $$

Figure S1C compares the distribution of the complexity measure across neurons for our model with $\sigma = \sigma_f$ with the one found in data and the one calculated for a population with linear tuning curves. For the sake of completeness, we also computed the KS distance between the distributions of $R^2$ corresponding to model and data (Figure S1A, red line). We mention that the model of33 with
two tunable parameters did not reproduce the distributions of complexity measure and of $R^2$, and only the complexity measure was taken into account in the fitting procedure. A better fit can be obtained in a heterogeneous model, at the cost of tracking many more parameters (two per neuron): see$^{86}$ for a more detailed discussion of the fitting procedure in such a model.

We also extract a noise model from the data, as follows. We define the variance of the mean response of neuron i across hand positions as the variance of the average responses across hand positions, $\hat{R}_i = \langle(\tilde{v}_i(\mathbf{x}) - \langle\tilde{v}_i(\mathbf{x})\rangle_i)^2\rangle_\mathbf{x}$, where $\tilde{v}_i(\mathbf{x})$ is the un-normalized tuning curve. Similarly, we average the trial-to-trial variability across different stimuli to obtain the variance of the noise, $\hat{\eta}_i^2 = \langle(r_i^t - \tilde{v}_i(\mathbf{x})_t)\rangle_i$, where $r_i^t$ is the response at trial t. In the model, we set the variance of the signal to unity and we rescaled the noise variance correspondingly, as

$$ \eta_i^2 = \frac{\hat{\eta}_i^2}{\hat{R}_i}. \qquad (\text{Equation 121}) $$

In principle, the noise may depend on the stimulus. To control for this effect, we preprocess the data with a variance stabilizing transformation. We substitute $r_i(\mathbf{x})$ by $\sqrt{r_i(\mathbf{x})}^{87}$, and we recalculated the variance of the noise accordingly. In this way, if the noise were proportional to the mean, one would obtain a constant estimate of the variance of the responses for different hand positions. The distribution of noise variances across neurons calculated in this way does not differ substantially from the one obtained without this data transformation.

For numerical simulations (Figure 6), the tuning curves are computed at a finer scale than in the data (cubic grid of 21 × 21 × 21 points instead of 3 × 3 × 3). We illustrate three examples of tuning curves obtained with $\sigma = \sigma_f$, measured at these hand positions in Figures S1D–S1F, together with the prediction obtained from a linear regression (Equation 9). We note that there are some neurons which are well described by the linear model while others are not compatible with it. We generated the tuning curves for a number of neurons equal to the number of neurons analyzed in the fitting procedure ($N_{tot} = 400$). Results for a given population size, N, are obtained by averaging over 8 different pools of size N sampled with replacement from $N_{tot}$. In Figures 6A–6C, we compare the MSE as obtained in a population in which neurons respond according to the irregular tuning curves generated by our model and a population in which the tuning curves are linear, Equation 9. The latter are generated according to Equation 9, by sampling the preferred directions, $\mathbf{p}_i$, uniformly on the unit sphere; the tuning curves are shifted and normalized to have zero mean and unit variance across hand positions. The comparison is quantified through the mean fractional improvement, defined as

$$ \Delta\epsilon \equiv \frac{\epsilon_{\text{lin}} - \epsilon_{\text{irr}}}{\epsilon_{\text{lin}}}, \qquad (\text{Equation 122}) $$

where $\epsilon_{\text{lin}}$ is the RMSE as obtained in the population with linear/irregular tuning curves.

## Dimensionality of a compressed neural code

A great deal of recent work aimed to characterize the response of a population of N neurons to an ensemble of stimuli, viewed as a geometric object made up of all the N-dimensional vectors of the mean responses. In the case of continuous stimuli, this geometric object is often referred to as a ‘neural response manifold’.51,52 Our theoretical framework has allowed us to discuss how the response properties of individual neurons relate to the structure of the neural response manifold and, in turn, how the structure of the latter impacts information processing. In this section, we relate the picture put forth by our model to recent work that quantifies the dimensionality of neural activity from neural recordings as a way to characterize its nature and to infer strategies used by the brain to represent (sensory) information.53,54 In doing so, we point out that simple methods of dimensionality reduction might miss the fine structure in neural response manifolds, which, in the case of complex and irregular neural responses, is relevant to information coding.73

With smooth tuning curves, a continuous K-dimensional stimulus is represented as an M-dimensional hypersurface embedded in the N-dimensional space of neural activity, with $M \le K$ (and $M = K$ when there are no invariances in the neural response). However, the complicated geometry of such an M-dimensional manifold may make its identification difficult. In practice, one is faced with a data set, namely a noisy sample from a population of tuning curves, and from it one would like to make statements on the geometry of the population activity. Fitting a low-dimensional manifold to a neural population data set is not a trivial task, and is the focus of a large number of studies on dimensional reduction.55 A simple approach is to consider the eigenvalue spectrum of the covariance matrix of the neural responses across the stimulus range or, equivalently, the variance carried by the different modes in a principal component analysis (PCA). If we apply this approach to the population response in our model, for different values of $\sigma$, we find a spectrum that exhibits a band-pass structure, which plateaus up to a cut-off value before a sharp suppression; the cut-off value is larger for smaller values of $\sigma$ (Figure S2A). From this analysis one would conclude that the population activity occupies a low-dimensional subspace embedded in the N-dimensional space of neural activity, with dimensionality controlled by $\sigma$. As the value of $\sigma$ falls to zero, the population responses fill an increasingly large fraction of the N available dimensions, until they fill the space entirely for $\sigma \to 0$.

A quantification of the ‘intrinsic dimensionality’ of the population activity based on this PCA analysis is offered by the participation ratio, defined as $d = (\sum_{i=1}^N \lambda_i^2) / (\sum_{i=1}^N \lambda_i^2)$, where $\lambda_i$ denotes the i-th eigenvalue of the covariance matrix of the neural responses across the stimulus range.88 Loosely speaking, the participation ratio measures the number of eigenvalues (principal components) which are much larger than the others; for example, if M eigenvalues are of comparable size and much larger than any others, then $d \approx M$.
In our model, while *d* is close to unity for large values of *σ*, it becomes larger for smaller values of *σ* and approaches *N* when *σ* → 0 (Figure S2B). It is interesting to examine the manifold dimensionality of the optimal neural code. In Figure S2C, we display the fractional dimensionality (i.e., the participation ratio divided by the number of neurons, *d*/*N*) corresponding to the population activity at the optimal value of *σ*, *σ*<sup>*</sup>, as a function of the population size, for a fixed level of noise. As expected, *d* increases with *N*: larger populations allow for more irregular tuning curves which benefit the local accuracy without generating prohibitive global errors. Quantitatively, the value of *d* hovers around *N*/2. A possible interpretation of this result is that it corresponds to the largest value beyond which a random manifold embedded in *N* dimensions comes close to intersect itself; thus, this value of *d* ensures that global errors do not proliferate. While a naive interpretation of the value of the participation ratio would suggest that the neural population encodes an *N*/2-dimensional stimulus, in the context of our model it results from the efficient coding of a one-dimensional stimulus. Thus, simple methods of dimensionality reduction might miss the fine structure in neural response manifolds. This observation also points to the difficulty of using a simple criterion to define the dimensionality of a manifold when the latter is highly non-linear.

## Model extensions

We proposed a two-layer feedforward neural network as a circuit mechanism that can generate downstream irregular tuning curves from upstream simple, bell-shaped ones. Samples from a Gaussian process with a Gaussian covariance function were chosen as a simple model of complex, irregular tuning curves because of analytical tractability and for the purpose of generality; as mentioned, a Gaussian process is the maximum-entropy, i.e., least-unbiased, stochastic process under the constraint of a fixed covariance function. Here, we discuss how different extensions to the model, motivated by features of real neurons, affect the coding performance in the population of representation neurons. The extensions we consider also point to the robustness of our qualitative results. For the sake of comparison, in all model extensions we always normalize tuning curves to have unit variance, to impose the resource constraint defined in Equation 3.

## Compressed coding in the presence of noise correlations

Real neurons exhibit correlated variability in their responses. In general, the coding accuracy can be either suppressed or enhanced by such noise correlations<sup>89</sup>; here, we illustrate two examples of correlations statistics which have qualitatively different, in fact opposite, effects on the coding accuracy in the population of representations neurons. In the formulation of the model used in previous sections, we assumed that sensory neurons were noiseless. (The rationale for this assumption was that we were interested in studying the impact of noise in the activity of representation neurons.) It is more natural, however, to assume that sensory neurons are also noisy. Noise in the activity of sensory neurons not only limits the information about the stimulus represented in the activity of the first layer, it also induces correlated noise in the activity of second-layer neurons. If neurons in the first layer are affected by independent, Gaussian noise with variance equal to ξ<sup>2</sup>, the covariance of the noise in the second layer becomes Σ = η<sup>2</sup>I + ξ<sup>2</sup>WW<sup>T</sup>.

As we shall show, in the compression regime (N ≪ L) on which we focus, the kind of correlations generated by noise in the sensory layer has a negligible effect on the coding performance. The presence of sensory noise degrades coding, so a comparison of noisy and noiseless systems is not very telling. Instead, we compare population coding in the presence of the full noise covariance matrix, Σ, and in the presence of a diagonal covariance matrix (i.e., independent noise) with elements chosen as follows. Given the distribution of synaptic weights, the matrix WW<sup>T</sup> is sampled from a Wishart distribution with mean given by the identity matrix and fluctuations of order 1/L (see Method details); in the limit of L → ∞, the covariance matrix becomes

$$ \Sigma_{\text{ind}} \equiv (\eta^2 + \xi^2)I = \tilde{\eta}^2 I, \qquad (\text{Equation 123}) $$

i.e., the population noise becomes independent, with single-neuron variance $\tilde{\eta} = \eta^2 + \xi^2$. Hereafter, we compare the two cases of populations with covariance matrices Σ and Σ<sub>ind</sub>.

In numerical studies, we observe, first, that the MSE depends only weakly on the noise correlations, as a function of *σ*. This behavior obtains because noise correlations primarily affect local errors, not global errors. (As noise correlations reduce the noise entropy — they ‘shrink the cloud of possible noisy responses’ — with respect to the independent case, one expects that correlations reduce the probability of occurrence of global errors. Numerical simulations however indicate that this effect is quantitatively negligible.)

We can show analytically that in our model, if noise correlations are due to independent noise in the sensory layer, local errors are enhanced. By computing a correction to the diagonal behavior of the covariance matrix in the limit *L* → ∞ (see Method details), we obtain the local contribution to the MSE as

$$ \epsilon_l^2 = \epsilon_{l, \text{ind}}^2 \left( 1 + \frac{N \xi_l^2}{L \tilde{\eta}^2} - \frac{N \xi_l^4}{L \tilde{\eta}^4} + \dots \right), \qquad (\text{Equation 124}) $$

where $\epsilon_{l, \text{ind}}^2$ is the corresponding quantity calculated for the matrix Σ<sub>ind</sub> rather than the full covariance matrix Σ. From Equation 124, it appears that the effect of noise correlations on the MSE is deleterious, but scales proportionally to the ratio between the two population sizes, which we supposed to be small. We checked this behavior numerically (Figure S3A), and found a good match with the analytical result. We also compared the impact of different values of ξ<sup>2</sup>, while keeping the effective noise variance, $\tilde{\eta}^2$, fixed (i.e., varying the relative contribution of input noise and output noise). Both Equation 124 and Figure S3B indicate that there exists a regime in which increasing the relative contribution of input noise, ξ<sup>2</sup>, in fact mitigates the deleterious effect of the correlated noise (this is seen in Equation 124 as a partial cancellation of the second- and fourth-order terms).
The harmful impact of noise correlation discussed above results from the form with which sensory noise invests it. To illustrate this point, we examine a case in which noise correlations in representation neurons come with an arbitrary structure unrelated to the network architecture. Specifically, we consider a noise covariance matrix that has the same statistics as those of $\Sigma$, but in which the form of correlations is not inherited from the network structure through the synaptic matrix $W$: a random covariance matrix, $\Sigma_{rand} = \eta^2 I + \xi^2 XX^T$, where $X_{ij} \sim N(0, 1/L)$. In this case, noise correlations suppress the MSE as compared to the independent case (with $\Sigma_{ind}$), because the ‘cloud of possible noisy responses’ is reoriented randomly with respect to the curve of mean responses. Analytically, the analog of Equation 124 for the case of a covariance matrix $\Sigma_{rand}$ (instead of $\Sigma$) is similar, but skips the lowest-order, deleterious term:

$$ \epsilon_{f,rand}^2 \approx \epsilon_{f,ind}^2 \left( 1 - \frac{N\xi^4}{L\tilde{\eta}^4} \right). \qquad (\text{Equation 125}) $$

This result, as well as numerical simulations (Figure S3B), demonstrates that generically coding is improved by random noise correlations, and that this improvement increases with $N$ and also increases with the relative contribution of $\xi^2$ with respect to $\eta^2$. In sum, noise correlations in representation neurons are deleterious if they are inherited from independent noise in sensory neurons—yet, the effect is quantitatively modest.

In a feedforward network such as the one we considered, correlated noise in representation neurons results from joint inputs from sensory neurons. In a more general network, recurrent connections can generate noise correlations. We have studied a generalization of our model that includes recurrence, to examine the impact of correlations thus generated. In the context of the network we analyzed, noise correlation that was induced by recurrent connectivity did not affect the coding accuracy in representation neurons (see below for details).

### Compressed coding with sparse and constrained synaptic weight matrices

The choice of a synaptic weight matrix with Gaussian entries is biologically unrealistic because (i) biological synaptic weights obey Dale's law, and (ii) real neurons are sparsely connected. In an extension of the model, we impose the condition that sensory neurons can be either excitatory or inhibitory by setting to 0 the negative (respectively, positive) elements of each column $i$ of the matrix $W$ with probability $f$ (respectively, $1-f$), corresponding to the fraction of excitatory neurons. This procedure has the additional effect of sparsifying the connectivity matrix, by setting to 0 half of its entries (sparsity = 0.5). We find for this case that the qualitative behavior of the error as a function of $\sigma$ is unchanged, with an optimal value of $\sigma$, that balances local and global errors, that decreases exponentially with the population size (Figure S4B). Quantitatively, we observe an increase of the MSE in the low-$\sigma$ regime. By analyzing a different extended model, in which the constraint on the sign of the weights is omitted, but the sparsity of the matrix is maintained, we can see that the effect is due to Dale's law (Figure S4C). The reason is that, under Dale's law, the correlation among the responses of different neurons to the same stimulus increases: in the extreme case in which only one sensory neuron responds to each stimulus value, neural responses in the second layer all have the same sign. Similarly, a rapid increase of the decoding error with decreasing tuning width is observed for appreciable sparsity of the matrix (Figure S4D). In this case, representation (second-layer) neurons become unresponsive over large stimulus intervals, due to the lack of inputs from the first layer. Correspondingly, the exponential scaling of the optimal decoding error saturates at a comparatively higher value of the tuning width of sensory neurons, as compared with the simpler (non-sparse) case.

To summarize, sparsity and sign constraint on synaptic weights do not affect quantitatively the coding properties of the network in the large-$\sigma$ regime. In this regime, neurons in the second layer receive a large number of random inputs, and hence their activities behave as Gaussian random variables; furthermore, the correlation length of the Gaussian process governing the tuning curves of representation neurons remains unchanged (compare Figures S4A, S4C and S4D). We note that when $\sigma$ is sufficiently small, the number of co-active inputs to each neuron becomes small, and Gaussianity no longer holds. This leads in (Figure S4) to some quantitative differences between the sparse and the non-sparse cases, in the small-$\sigma$ regime. Specifically, the optimum value of $\sigma$ is slightly greater in the sparse case, but the qualitative features of the encoding scheme are unchanged.

### Compressed coding with non-linear neural responses

In the network model we proposed, tuning curves have zero mean, resulting in half of neural responses which are negative. These responses can be interpreted as deviations away from a baseline firing rate. Alternatively, we can interpret that quantity as the membrane potential, which is then passed through a non-linear response function, $g$, to obtain non-negative firing rates. The tuning curve of neuron $i$ is then obtained as $v_i(x) = A_{nl}g(v_i^f(x))$, with $v_i^f(x)$ the ‘linear’ tuning curve, Equation 2, and the normalization constant, $A_{nl}$, computed numerically to impose the resource constraint, Equation 3. We consider the effect of two common choices of the form of $g$ on the coding performance of the representation neurons.

In order to model neurons that respond only to a fraction of stimuli, as is the case for example of the place cells observed by,²⁸ one possibility is to impose a threshold-linear response function, as

$$ v_i(x) = A_{nl} \left[ \sum_{j=1}^{L} W_{ij} u_j(x) + c_i \right], \qquad (\text{Equation 126}) $$

where $[\cdot]_+$ = max(0, ·) and $c_i$ is a firing bias set so as to fix the ‘coding level,’ defined as the fraction of stimuli which evoke a non-vanishing neural activity. For a given stimulus, a lower coding level is accompanied by a larger fraction of silent neurons; this results in
lowering the effective number of neurons encoding the stimulus, and in turn enhances the MSE in the regime in which it scales exponentially with N (Figures S5B and S5C). The normalization constraint has the effect of increasing the amplitude of neural responses, to compensate for the fraction of the non-responsive neurons. In the regime of sufficiently large σ and large N, in which local errors dominate, this results in a lower MSE (reduced by a factor of ~ 0.85 and ~ 0.9 with respect to the linear case when the coding level is = 0.3 and 0.5, respectively).

An alternative common choice of a non-linear response function is an exponential, $g(\cdot) = \exp(\cdot)$; this choice yields the prior distribution of tuning curves used in$^{63,64}$ to fit data from neural recordings. Numerically, we observe that this choice has a similar effect as setting a low coding level above (Figure S5D). We obtain an intuitive understanding of this behavior by inspecting the resulting statistics of the tuning curves: neural responses are distributed according to a multivariate log-normal distribution (see Method details, Equations 18 and 19). The distribution of responses is 'squeezed' in a comparatively narrow interval of values, resulting in an effective lower number of neurons; moreover, the width of the correlation function of the tuning curves is reduced compared to the linear case (Figure S5). This yields a lower MSE in the regime in which local errors dominate, similarly to the case with threshold non-linearity and low coding level.

Overall, the introduction of non-linearity in the neural responses preserves the qualitative the coding properties of the bare model, characterized by an optimal σ that balances local and global errors, and yields an exponential scaling of the MSE as a function of the population size. Quantitatively, it has the effect of suppressing the effective number of neurons encoding the stimulus, and of enhancing the optimal tuning width of tuning curves in the first layer (Figure S5).

## Compressed coding in a population of neuron with recurrent connectivity

Our focus has been on feedforward processing, relevant to a number of brain computations, as in hippocampal area CA1 and in the context of fast processing in the visual cortex.$^{90,91}$ Yet, it is interesting to examine also how our framework can apply to a situation in which recurrent connections affect the neural activity in the representation layer. A model that incorporates recurrent connections needs to be formulated in a dynamical setting which allows for the reverberation of activity induced by recurrence. We adopt a standard description of the dynamics of the form (see Method details Equation 102)

$$ \tau \frac{ds_i}{dt} = -s_i + \sum_j M_{ij}s_j + \sum_j W_{ij}u_i(x) + \zeta_i(t), \qquad (\text{Equation 127}) $$

where $s_i$ is the synaptic activity of neuron i in the representation layer, $\tau$ the synaptic time constant of integration, M the matrix of recurrent synaptic weights, and $\zeta_i$ denotes neural noise that arises from input stochasticity, i.e., dynamic variability in the neural activity in the sensory layer. For the sake of simplicity, we assume that this noise is uncorrelated among neurons, but the argument we provide below can be extended to the case of correlated noise arising, e.g., from noise in the activity of the input neurons. We assume furthermore that the recurrent connectivity does not destabilize neural responses under the linearized dynamics; this implies that the matrix $H = I - M$ is invertible.$^5$

We consider a scheme in which coding and decoding occurs on a time scale much longer than the synaptic time scale, in order to examine the effect of recurrence. (On short time scales, recurrence has less relevance; we return to this point below.) Thus, we focus on the steady-state activity in the recurrent representation layer induced by a constant input to the sensory layer. The activity in the representation layer settles in the steady state on a mean response that depends on its inputs and the feedforward and recurrent weights: $\bar{s} = H^{-1}Wu(x)$, in vector notation, where $\bar{s}$ denotes the vector of mean responses (constant in time). The correlation structure of the noise is affected by the recurrent weights as well, due to the reverberation of the noise through the recurrent network dynamics. The coding properties depend on both the transformation of the inputs into mean steady-state responses and the correlation structure of the noise. Since the mean responses are constant (in the steady state) and the noise fluctuates over time, it is advantageous for a decoder to integrate the neural activity in the representation layer over a time scale $T \gg \tau$. We thus consider the quantity defined as

$$ r_i(t) = \int_{-\infty}^{t} \frac{dt'}{T} \exp \left( - \frac{(t - t')}{T} \right) s_i(t') \qquad (\text{Equation 128}) $$

as the relevant informative variable. The signal to noise ratio in $r_i$ scales proportionally to T if the time scale of the temporal correlations in the noise is smaller than T. As we show in Methods details, for large values of T this variable takes on a simple expression,

$$ r = H^{-1}Wu(x) + \xi, \qquad (\text{Equation 129}) $$

in vector notation. Here, $\xi$ is a Gaussian correlated term, which is obtained by passing uncorrelated noise through the transformation $H^{-1}$ (another way to state this is that $H\xi$ is independent among neurons). Hence, the noise is linearly transformed by the recurrent activity in the same way the signal is. Since we assume that $H$ is invertible, neural responses can be transformed deterministically into a signal governed by Equation 26 by applying the linear transformation $H$ to $r$, and all of our results on the decoding precision hold without modification.

If the time window of readout is short, the signal and the noise do not transform in the same way. Applying $H$ to $r$ still yields the correct mean readout, but $\xi$, the noise term in Equation 129, takes on a non-trivial correlation structure that depends on the specifics of the recurrent weight matrix. In this case, the impact of recurrence on coding properties depends on specific assumptions about the
structure of the recurrent architecture. For fast readout, however, the situation is closer to that of a feedforward network. We comment further on the issue of recurrence at different time scales, as well as on the inclusion of non-linear dynamics in the presence of recurrence, in the Discussion and in the next section, but leave a detailed analysis for future work.

## Generality and use of the network model

In sum, what matters in our derivations is that the tuning curves are sampled from a Gaussian process. The simple network we put forth serves as an example of a mechanism that can generate such tuning curves. We proposed a two-layer feedforward neural network as a simple model mechanism that generates complex tuning curves with the desired statistics from the response of sensory neurons with the usual bell-shaped tuning curves. In our simple model mechanism, the irregularity of tuning curves arises from random synaptic weights. Besides the fact that random synaptic weights have been used in many theoretical studies of information processing through layers of neural networks,⁹²⁻⁹⁶ random connectivity has also been proposed as a computational device in the context of observed neural populations (see, e.g.,⁹⁷). In sum, what matters in our derivations is that the tuning curves are sampled from a Gaussian process. The simple network we put forth serves as an example of they can be generated mechanistically. It is also useful, however, to think about possible extensions of the model (see the last section of Results and Supplementary Information for a discussion of some of these). Quantitatively, we find that some of the features with which our simple model can be enriched, such as sparsity in the connectivity matrix and non-linear response functions, lower the ‘effective’ number of coding neurons; they also rule out the unbiological prediction of extremely narrow tuning curves with $σ^* → 0$ for large population sizes (Figures S4 and S5).

In the context of the two-layer architecture, we considered an extension to a network with recurrent synaptic connectivity within the representation (output) layer. We showed that our results on coding apply in a ‘slow coding’ regime in which a downstream decoder integrates the output of the network over a time scale longer than the synaptic time scale. We have, however, made a number of simplifications. The main one is that we analyzed a linearized description of the recurrent dynamics. In general, however, linear networks may exhibit instabilities, depending on the strength and structure of the recurrent connections, and, in such cases, nonline- arities play an important role in stabilizing the dynamics. In this more involved case, both the effective structure of the noise and the interplay between mean response and noise, which govern the coding properties, can become intricate. For example, in balanced networks with random connectivity,⁹⁸,⁹⁹ noise may become uncorrelated in different neurons, but this is not the case in other classes of recurrent network models. At a conceptual level, one still expects the mean population activity to span a low-dimensional manifold as the low-dimensional stimulus varies, and efficient coding to arise when this manifold densely fills the population output space in a self-avoiding manner. Thus, even in recurrent and non-linear networks, random synaptic weights can support efficient encoding. It will be of interest to investigate specific instances of network dynamics that can instantiate such a scenario in future work.

## Limitations of the study (extended)

Our main results are derived under the assumptions that tuning curves are realizations of a Gaussian process and that neurons in the sensory layer exhibit uncorrelated noise. We examined a few extensions of the model but others remain to be considered, especially in applying the theory to specific neural systems. For example, in many brain areas, neural noise is correlated. While some patterns of noise correlations can have only a marginal effect on information coding, as discussed above, others can yield a more drastic ef- fect⁸⁰,⁸¹ which warrants study in the context of the coding scheme analyzed in the present paper. In this context, it would be partic- ularly relevant to study the patterns of correlation arising from correlated noise in the sensory layer or from recurrent, non-linear cou- plings in the representation layer. Indeed, while we consider here a purely feedforward network, sufficient for our purposes, a heterogeneous recurrent architecture is another, natural candidate for the emergence of extended, irregular tuning curves. In that more involved scheme, a thorough examination of the role of noise correlation would be relevant.

Given the theoretical advantages of extended, irregular tuning curves, it is natural to wonder about their prevalence in the brain. The historical focus on simple tuning curves might be partly due to a measurement bias: they are easier to identify and their corresponding sensitivity is more directly amenable to function. Recently, however, neurons with complex selectivity have been recorded in several different regions of the brain, as we mention in the Introduction, and we can ask whether they share common principles of information processing. We provide some elements of an analysis of the representation of information in monkey motor cortex in the light of our theoretical proposal, but there is a long way to go before we can answer the above, wide-ranging question. We note here, however, that while much remains to be discovered as to the functional role of extended, irregular tuning curve, strong evidence has now become available that the tuning curves of a class of real neurons are described in quantitative detail by the model proposed in the present paper. The statistics of sizes, arrangements, shapes, and topological properties of hippocampal CA1 neurons are captured with precision by our model (⁴⁰, see also²⁸). These findings support the feedforward architecture put forward here as a possible generator of realizations of a Gaussian process, since there is very little excitatory connectivity between CA1 pyramidal cells.¹⁰⁰ Furthermore, recent measurements in bats indicate that CA3 place cells, which provide the main excitatory drive to CA1, are sparsely tuned even in very large environments in which CA1 tuning curves are multimodal. These findings suggest that CA1 com- presses the inputs from CA3, in accordance with our theoretical framework, to generate a neural representation that can be reliably readout from a relatively small subset of the cells.

An important question we leave for future study pertains to decoding from extended, irregular tuning curves. In our work, we have assumed an optimal decoder—an unrealistic assumption as it requires the decoder to have access to the full knowledge of the tuning curve shapes. As we show in work to be published, this requirement can be circumvented because sub-optimal, but still faithful,
decoders can be learned from experience. It remains, though, that encoding with extended, irregular tuning curves requires heavier decoding machinery than encoding with simple, unimodal tuning curves, and learning it can be much more data-intensive. In a more complete functional theory, one would like to consider coding performance in a joint encoder-decoder system.⁸² In that richer picture, encoding optimality depends on decoding, and the encoding scheme might switch from one mode to another—simple vs. complex—according to the specific function of the neural assembly under scrutiny and the availability of decoding resources.

## QUANTIFICATION AND STATISTICAL ANALYSIS

Numerical simulations and data analyses were carried out with custom codes written in Julia,¹⁰¹ with the DrWatson project assistant.¹⁰²

### Details of numerical simulations

In numerical simulations, we compute the MSE with standard Monte Carlo methods. At each step, we sample a stimulus, we generate a noisy population response and we decode it using the ideal decoder; the squared difference between the stimulus and its estimate is used to update the MSE. This process is repeated and the MSE estimate is updated until convergence, defined as the point for which the variance of the MSE estimates in the last 500 steps, after a burn-in period of 5000 steps, is less than a tolerance threshold, set to $10^{-8}$. We set the number of decoder neurons equal to the number of sensory neurons, $M = L$, with uniformly spaced preferred stimuli, $x_m = m/M$. Unless otherwise stated, $L = 500$, $R = 1$ and $\eta^2 = 0.5$. The results are averaged over 8 network realizations and shaded regions correspond to one SD.

## EXPERIMENTAL MODEL AND STUDY PARTICIPANT DETAILS

This is a computational study.