text:
# A theoretical model of neuronal population coding of stimuli with both continuous and discrete dimensions

## Introduction

The strategy used by populations of neurons to code for external stimuli or behavioural correlates is a major issue which has been recently investigated both through data analysis and theoretical modelling. The mutual information between external correlates and the spiking activity of the population is one way to assess such coding quantitatively . Several analyses have focused on the coding of a discrete set of stimuli (, see  for a review), which is the paradigm used in many experiments. In this situation the mutual information is bounded by the entropy of the stimulus set. Some theoretical studies have also considered the coding of stimuli varying in a continuous domain , which is interesting with respect to basic properties like orientation in visual stimuli, frequency in auditory stimuli, velocity and position in motor actions. In particular in  the authors have studied the asymptotic (large population) behaviour of the mutual information, with respect to a stimulus with a continuously varying dimension.

Yet no study has been proposed so far considering a mixture of both continuous and discrete features, which is obviously closer to real world stimuli or behavioural correlates. Moreover the initial rise of the mutual information for small but increasing population size is more relevant for a comparison with estimates from real data, at least as far as the possibility of having simultaneous recordings from very large populations of neurons is restricted to very few cases.

We have recently analyzed data recorded in the motor areas of behaving monkeys, in the laboratory of Eilon Vaadia. The monkeys moved a manipulandum in several possible directions (approximating a continuous correlate) and with different combinations of arms (4 types of movement, i.e. a discrete correlate). In trying to characterize the neural coding of these movements, we were particularly interested in whether, as it is reasonable to expect, different motor areas differ, at least quantitatively, in their coding properties. The results, obtained from records of activity in areas M1 and SMA, will be reported elsewhere ; but they suggest the importance of developing theoretical models of how populations of neurons might code simultaneously discrete and continuous correlates. For example, one clear conclusion has been that type and direction are not independent dimensions of the movement, in the sense for example that the information about direction, extracted from the activity recorded with all movement types, is much lower (roughly half) of the average information about direction, obtained with a single movement type.

Can we embody similar properties in a model of the scheme used by neurons to code movements? What would then be the dependence of the mutual information on the number of the possible types of movement? How would it depend on the resolution with which the continuous correlate is sampled? How on the level of noise affecting the firing patterns?

In a recent work  some of these questions have been investigated for a set of $`p`$ discrete stimuli, under the assumption that neurons fire independently of one another to each stimulus and that the distribution of the firing rates is gaussian. The linear and the quadratic approximations to the mutual information as a function of population size were studied analytically, in the limit of large noise, as well as the approach to the ceiling in the case of small noise. We generalize this study considering both a discrete and a continuous dimension in the stimulus, referring specifically to motor actions characterized by a direction and a “type”. Nonetheless, our model is equally applicable to other complex correlates. We introduce a more realistic conditional firing rate distribution, than the simple gaussian one, and find a simple resulting correction to the gaussian model: the analytical expression of the mutual information remains the same, except for a renormalization of the expansion parameter.

We then evaluate the information loss when the original activity distribution is averaged across the discrete correlate, as is sometimes the case in the analysis of real data, and the mutual information is evaluated solely with respect the continuously varying feature. Averaging out dimensions in the stimulus corresponds to losing accuracy in its description, and hence the information loss.

Our theoretical analysis allows a direct comparison with real curves; we present one comparison and discuss possible causes for the discrepancy between data and model. In particular, correlations between neurons, that are not included in the model, might play a relevant role, enhancing or decreasing redundancy in population coding . This issue will be the object of future work.

## The gaussian model

First, we consider a coding scheme where the distribution of the firing rates conditional to the movements is gaussian, similarly to the case examined in . This assumption implies that negative rates have a non zero probability to occur, but it allows an easier analytical treatment. We will examine a more realistic scheme later on.

Consider the following distribution:

$$P(\{\eta _i\}|\vartheta ,s)=\underset{i=1}{\overset{N}{}}\frac{1}{\sqrt{2\pi \sigma ^2}}exp\left[\left(\eta _i\stackrel{~}{\eta }_i(\vartheta ,s)\right)^2/2\sigma ^2\right];$$
(1)
$`\eta _i`$ is the firing rate of neuron $`i`$; $`\vartheta `$ and $`s`$ parameterize respectively the direction and the type of movement; $`\stackrel{~}{\eta }_i(\vartheta ,s)`$ is the average firing rate of the neuron with the movement parameterized by $`\vartheta `$,$`s`$.

In general, the directional tuning of real cells in motor cortices is modulated by the type of movement performed. We show an example of this modulation, with the typical shape of tuning curves, in fig.1 (data kindly provided by Eilon Vaadia). The modulation of the preferred direction looks weaker than the overall amplitude modulation.

For our model we have considered the following function:

$$\stackrel{~}{\eta }_i(\vartheta ,s)=\epsilon _s^i\overline{\eta }_i(\vartheta )(1\epsilon _s^i)\eta _i^f;$$
(2)
$$\overline{\eta }_i(\vartheta \vartheta _i^0)=\eta _i^0\mathrm{cos}^{2m}\left(\frac{\vartheta \vartheta _i^0}{2}\right)=\eta _i^0\left[\frac{1}{2^{2m}}\left(\begin{array}{c}2m\\ m\end{array}\right)+\frac{1}{2^{2m1}}\underset{\nu =0}{\overset{m1}{}}\left(\begin{array}{c}2m\\ \nu \end{array}\right)\mathrm{cos}\left\{\left(m\nu \right)\left(\vartheta \vartheta _i^0\right)\right\}\right];$$
(3)
$`\epsilon _s^i`$ is a quenched random variable distributed between $`0`$ and $`1`$. The meaning of eq.(1),(2),(3) is that the firing rate of neuron $`i`$ given the movement parameterized by ($`\vartheta `$,$`s`$) follows a gaussian distribution centered around a tuning curve $`\overline{\eta }_i(\vartheta )`$ whose flatness is modulated through the parameters $`\epsilon _s^i`$. If $`\epsilon _s^i`$ is zero for some particular $`s`$ the firing of the cell does not depend, for that movement type, on the direction of the movement. On the other hand if $`\epsilon _s^i`$ assumes a fixed value for all $`s`$ the directional tuning does not modulate with the type of movement. Tuning curves with a cosinusoidal shape have already been considered to model the directional selectivity of sensory neurons . Fig.2 shows the amplitude of the tuning curve.

We neglect the modulation of the preferred direction with the type, as it would burden the analytical calculations; it will be the object of future analyses.

## Evaluation of the mutual information for the gaussian model

We are interested in the mutual information  between the neuronal firing rates and the movements:

$$I(\{\eta _i\},\vartheta s)=\underset{s=1}{\overset{p}{}}𝑑\vartheta \underset{i}{}d\eta _iP(\vartheta ,s)P(\{\eta _i\}|\vartheta ,s)\mathrm{log}_2\frac{P(\{\eta _i\}|\vartheta ,s)}{P(\{\eta _i\})}_{\epsilon ,\vartheta ^0};$$
(4)
where the distribution $`P(\eta _i|\vartheta ,s)`$ is given in eq.(1) and $`.._{\epsilon ,\vartheta ^0}`$ is a short notation for the average across the quenched variables $`\{\epsilon _s^i\}`$,$`\{\vartheta _i^0\}`$. In fact we are not interested in a particular realization of the tuning, but in the average across all its possible realizations.

Eq.(4) can be written as:

$$I(\{\eta _i\},\vartheta s)=H(\{\eta _i\})H(\{\eta _i\}|\vartheta ,s)_{\vartheta ,s};$$
(5)
$$H(\{\eta _i\}|\vartheta ,s)_{\vartheta ,s}=\underset{s=1}{\overset{p}{}}𝑑\vartheta \underset{i}{}d\eta _iP(\vartheta ,s)P(\{\eta _i\}|\vartheta ,s)\mathrm{log}_2P(\{\eta _i\}|\vartheta ,s)_{\epsilon ,\vartheta ^0};$$
(6)
$$H(\{\eta _i\})=\underset{s=1}{\overset{p}{}}𝑑\vartheta \underset{i}{}d\eta _iP(\vartheta ,s)P(\{\eta _i\}|\vartheta ,s)\mathrm{log}_2\left[\underset{s^{}=1}{\overset{p}{}}𝑑\vartheta ^{}P(s^{},\vartheta ^{})P(\{\eta _i\}|\vartheta ^{},s^{})\right]_{\epsilon ,\vartheta ^0};$$
(7)
The calculation of the equivocation $`H(\{\eta _i\}|\vartheta ,s)_{\vartheta ,s}`$ is straightforward and the result is additive in the population size:

$$H(\{\eta _i\}|\vartheta ,s)_{\vartheta ,s}=\frac{N}{2\mathrm{ln}2}[1+\mathrm{ln}(2\pi \sigma ^2)];$$
(8)
The linearity in $`N`$ is standard whenever the conditional distribution of firing rates can be factorized across neurons, as in eq.(1).

The evaluation of the rate entropy $`H(\{\eta _i\})`$ can be carried out introducing $`n`$ replicas  for both the discrete and the continuous dimensions $`\vartheta `$ and $`s`$, which allows to get rid of the logarithm in eq.(7):

$`H(\{\eta _i\})={\displaystyle \frac{1}{\mathrm{ln}2}}\underset{n0}{lim}{\displaystyle \frac{1}{n}}({\displaystyle \underset{s_1..s_{n+1}=1}{\overset{p}{}}}{\displaystyle }d\vartheta _1..d\vartheta _{n+1}{\displaystyle }{\displaystyle \underset{i}{}}d\eta _i{\displaystyle \frac{1}{(2\pi p)^{n+1}}}`$
$`\text{x}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{k=1}{\overset{n+1}{}}}{\displaystyle \frac{1}{\sqrt{2\pi \sigma ^2}}}exp[(\eta _i\epsilon _{s_k}^i\overline{\eta }_i(\vartheta _k)(1\epsilon _{s_k}^i)\eta _i^f)^2/2\sigma ^2]_{\epsilon ,\vartheta ^0}1);`$ (9)
Integrating over $`\{\eta _i\}`$ and rearranging terms one obtains:

$$H(\{\eta _i\})=\frac{1}{\mathrm{ln}2}\underset{n0}{lim}\frac{1}{n}(\underset{s_1..s_{n+1}=1}{\overset{p}{}}d\vartheta _1..d\vartheta _{n+1}\frac{(n+1)^{\frac{N}{2}}}{(2\pi p)^{n+1}}\frac{1}{(\sqrt{2\pi \sigma ^2})^{nN}}\underset{i=1}{\overset{N}{}}exp(R_i)_{\epsilon ,\vartheta ^0}1);$$
(10)
$$R_i=\underset{k,l}{}\left(\epsilon _{s_k}^i(\overline{\eta }_i(\vartheta _k)\eta _i^f)\epsilon _{s_l}^i(\overline{\eta }_i(\vartheta _l)\eta _i^f)\right)^2/(4\sigma ^2(n+1));$$
(11)
### Large $`𝝈`$ limit

An exact analytical evaluation of eq.(10) is not possible without resorting to some approximation. In line with the analysis performed in  we assume now that the quenched randomness is uncorrelated and identically distributed across neurons:

$$\varrho (\{\epsilon _s^i\})=\left(\underset{s}{}\varrho (\epsilon _s)\right)^N;$$

$$\varrho (\{\vartheta _i^0\})=\left[\varrho (\vartheta ^0)\right]^N=\frac{1}{(2\pi )^N};$$

We assume also that $`\eta _i^f=\eta ^f`$ $`i`$. Then one can write:

$$\underset{i}{}exp(R_i)_{\epsilon ,\vartheta ^0}=exp(R)_{\epsilon ,\vartheta ^0}^N;$$

We consider now the limit of large noise $`\sigma `$; in this case, since $`R1/\sigma ^2`$ we can expand $`exp(R)`$; keeping only terms of order $`(N/\sigma ^k)^l`$, with $`k2`$ and $`l=1,2`$ we obtain:

$`exp(R)_{\epsilon ,\vartheta ^0}^N1NR_{\epsilon ,\vartheta ^0}+{\displaystyle \frac{N^2}{2}}R_{\epsilon ,\vartheta ^0}^2=`$
$`1{\displaystyle \frac{N}{4\sigma ^2(n+1)}}{\displaystyle \underset{k,lk}{}}\left(\epsilon _{s_k}(\overline{\eta }(\vartheta _k)\eta ^f)\epsilon _{s_l}(\overline{\eta }(\vartheta _l)\eta ^f)\right)^2_{\epsilon ,\vartheta ^0}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{N^2}{(4\sigma ^2(n+1))^2}}`$
$`{\displaystyle \underset{k,lk}{}}{\displaystyle \underset{\varrho ,\mu \varrho }{}}\left(\epsilon _{s_k}(\overline{\eta }(\vartheta _k)\eta ^f)\epsilon _{s_l}(\overline{\eta }(\vartheta _l)\eta ^f)\right)^2_{\epsilon ,\vartheta ^0}\left(\epsilon _{s_\varrho }(\overline{\eta }(\vartheta _\varrho )\eta ^f)\epsilon _{s_\mu }(\overline{\eta }(\vartheta _\mu )\eta ^f)\right)^2_{\epsilon ,\vartheta ^0}`$ (12)
To first order in $`N/\sigma ^2`$ we obtain:

$`H(\{\eta _i\})={\displaystyle \frac{1}{\mathrm{ln}2}}\underset{n0}{lim}{\displaystyle \frac{1}{n}}({\displaystyle \underset{s_1..s_{n+1}=1}{\overset{p}{}}}{\displaystyle }d\vartheta _1..d\vartheta _{n+1}{\displaystyle \frac{(n+1)^{\frac{N}{2}}}{(2\pi p)^{n+1}}}{\displaystyle \frac{1}{(\sqrt{2\pi \sigma ^2})^{nN}}}`$
$`(1{\displaystyle \frac{N}{4\sigma ^2(n+1)}}{\displaystyle \underset{k,lk}{}}(\epsilon _{s_k}(\overline{\eta }(\vartheta _k)\eta ^f)\epsilon _{s_l}(\overline{\eta }(\vartheta _l)\eta ^f))^2_{\epsilon ,\vartheta ^0})1);`$ (13)
This result is valid for a generic directional tuning curve $`\overline{\eta }(\vartheta )`$. We consider now our specific choice, eq.(3), examining the simplest case $`m=1`$ first. After averaging across direction selectivities $`\{\vartheta ^0\}`$ and integrating over continuous replicas $`\vartheta _1..\vartheta _{n+1}`$ we obtain

$`H(\{\eta _i\})={\displaystyle \frac{1}{\mathrm{ln}2}}\underset{n0}{lim}{\displaystyle \frac{1}{n}}({\displaystyle \underset{s_1..s_{n+1}=1}{\overset{p}{}}}{\displaystyle \frac{(n+1)^{\frac{N}{2}}}{(p)^{n+1}}}{\displaystyle \frac{1}{(\sqrt{2\pi \sigma ^2})^{nN}}}`$
$`\text{x}(1{\displaystyle \frac{N(\eta ^0)^2}{4\sigma ^2(n+1)}}{\displaystyle \underset{k,lk}{}}[(\alpha (\alpha 1)+{\displaystyle \frac{1}{4}})(\epsilon _{s_k}\epsilon _{s_l})^2_\epsilon +{\displaystyle \frac{1}{8}}\epsilon _{s_k}^2+\epsilon _{s_l}^2_\epsilon ])1);`$
where we have defined $`\alpha =\eta ^f/\eta ^0`$. To perform the average across the quenched variables $`\{\epsilon _s\}`$ we assume that they are equally distributed across $`p`$ movement types, namely:

$$\varrho (\{\epsilon _s\})=\underset{s}{}\varrho (\epsilon _s)=[\varrho (\epsilon )]^p;$$

In this case the sum over indexes $`k`$ and $`l`$ generates $`n(n+1)`$ identical terms. The summation on discrete replicas yields a factor $`p^{n+1}`$ multiplying the term $`(\epsilon _{s_k})^2_\epsilon `$ and a factor $`p^n(p1)`$ multiplying the term $`(\epsilon _{s_k}\epsilon _{s_l})^2_\epsilon `$, since this last term is non zero only when $`s_ks_l`$.

Taking the limit $`n0`$ yields

$`H(\{\eta _i\}){\displaystyle \frac{N}{2\mathrm{ln}2}}[1+\mathrm{ln}(2\pi \sigma ^2)]`$
$`+{\displaystyle \frac{1}{\mathrm{ln}2}}{\displaystyle \frac{N(\eta ^0)^2}{4\sigma ^2}}\left[{\displaystyle \frac{p1}{p}}2\left(\alpha (\alpha 1)+{\displaystyle \frac{1}{4}}\right)\lambda _1+{\displaystyle \frac{1}{4}}\lambda _2\right];`$ (14)
$$\lambda _1=𝑑\epsilon \varrho (\epsilon )\epsilon ^2\left[𝑑\epsilon \varrho (\epsilon )\epsilon \right]^2;$$
(15)
$$\lambda _2=𝑑\epsilon \varrho (\epsilon )\epsilon ^2;$$
(16)
From eqs.(14) and (8) the final expression for the mutual information can be written

$$I(\{\eta _i\},\vartheta s)\frac{1}{\mathrm{ln}2}\frac{N(\eta ^0)^2}{4\sigma ^2}\left[\frac{p1}{p}2\left(\alpha (\alpha 1)+\frac{1}{4}\right)\lambda _1+\frac{1}{4}\lambda _2\right];$$
(17)
In the more general case of a power $`2m`$ of the cosine, in eq.(3), it is easy to show that the final result can be expressed as

$$I(\{\eta _i\},\vartheta s)\frac{1}{\mathrm{ln}2}\frac{N(\eta ^0)^2}{4\sigma ^2}\left[\frac{p1}{p}2\left(\alpha A_1\right)^2\lambda _1+2\left(A_2(A_1)^2\right)\lambda _2\right];$$
(18)
where:

$$\eta ^0A_1=\overline{\eta }(\vartheta _k)_{\vartheta ^0}=\eta ^0\frac{1}{2^{2m}}\left(\begin{array}{c}2m\\ m\end{array}\right);$$
(19)
$$\left(\eta ^0\right)^2A_2=\overline{\eta }^2(\vartheta _k)_{\vartheta ^0}=\left(\eta ^0\right)^2\frac{1}{2^{4m}}\left(\begin{array}{c}4m\\ 2m\end{array}\right);$$
(20)
The calculation of the coefficient of the second order term, which multiplies $`N^2/\sigma ^4`$ in eq.(12), is slightly more complex and implies integration of terms with 4-replica interaction. The detailed analytical evaluation is given in the appendix. The final result up to the quadratic approximation reads:

$`I(\{\eta _i\},\vartheta s){\displaystyle \frac{1}{\mathrm{ln}2}}\{{\displaystyle \frac{N(\eta ^0)^2}{4\sigma ^2}}[{\displaystyle \frac{p1}{p}}2(\alpha A_1)^2\lambda _1+2(A_2(A_1)^2)\lambda _2]{\displaystyle \frac{N^2(\eta ^0)^4}{2(4\sigma ^2)^2}}`$ (21)
$`\{{\displaystyle \frac{p1}{p^2}}2\left(2(\alpha A_1)^2\lambda _1\right)^2+[{\displaystyle \frac{p1}{p}}(\lambda _1\lambda _2)^2+{\displaystyle \frac{(\lambda _2)^2}{p}}]\left[\left({\displaystyle \frac{1}{2^{2m1}}}\right)^4{\displaystyle \underset{\nu =0}{\overset{m1}{}}}\left[\left(\begin{array}{c}2m\\ \nu \end{array}\right)\right]^4\right]\}\};`$ (24)
where the expressions of $`\lambda _1`$,$`\lambda _2`$,$`A_1`$,$`A_2`$ are given respectively in eq.(15),(16),(19),(20) In the limit of large $`p`$ we have

$`I(\{\eta _i\},\vartheta s){\displaystyle \frac{1}{\mathrm{ln}2}}\{{\displaystyle \frac{N(\eta ^0)^2}{4\sigma ^2}}[2(\alpha A_1)^2\lambda _1+2(A_2(A_1)^2)\lambda _2]`$
$`{\displaystyle \frac{N^2(\eta ^0)^4}{2(4\sigma ^2)^2}}\left\{(\lambda _1\lambda _2)^2\left[\left({\displaystyle \frac{1}{2^{2m1}}}\right)^4{\displaystyle \underset{\nu =0}{\overset{m1}{}}}\left[\left(\begin{array}{c}2m\\ \nu \end{array}\right)\right]^4\right]\right\}\};`$ (27)
Fig.3 shows the linear and the quadratic approximations of eq.(24) for different values of the expansion parameter $`g^2=\left(\eta ^0/2\sigma \right)^2`$. It is easy to see that, for very small values of $`g^2`$, linear and quadratic approximations roughly coincide, while when $`g^2>0.8`$ the quadratic approximation begins to fail and one should add higher orders in perturbation theory.

Fig.4(a) shows the dependence of the mutual information on the number of types $`p`$. The dependence on $`p`$ is weak for the linear approximation and somewhat stronger for the quadratic one. In both cases an increase in the number of discrete correlates $`p`$ produces an increase in the mutual information. The distance between the linear and the quadratic approximation remains asymptotically finite (for $`p\mathrm{}`$), contrary to what happens in the case of discrete stimuli alone .

Fig.4(b) shows the dependence of the mutual information on the width of the directional tuning (see fig.2,(b)). Since we are considering the case when the noise $`\sigma `$ is large, a very narrow tuning in $`\vartheta `$ corresponds to a larger overlap in the conditional probabilities $`p(\eta |\vartheta ,s)`$, for most angles $`\vartheta `$. A consequence is that, especially in the linear approximation, the mutual information is a (slowly) decreasing function of $`m`$.

### The limit of large $`N`$

We consider now the case when the number of neurons is large. Since we deal with an infinite number of stimuli the mutual information is unbounded. Thus we expect that when the number of neurons becomes large and the noise is finite the mutual information tends asymptotically to infinity.

In order to study this limit we discretize the $`\{\vartheta \}`$ space into a finite set of $`M=2\pi /\mathrm{\Delta }\vartheta `$ angles $`\vartheta _1..\vartheta _M`$, and then we take the limit $`\mathrm{\Delta }\vartheta 0`$.

The entropy of the responses $`H(\{\eta _i\})`$ can be written

$$H(\{\eta _i\})=\underset{s=1}{\overset{p}{}}\underset{k=1}{\overset{M}{}}\underset{i}{}d\eta _iP(\vartheta _k,s)P(\{\eta _i\}|\vartheta _k,s)\mathrm{log}_2\left[\underset{s^{}=1}{\overset{p}{}}\underset{k^{}=1}{\overset{M}{}}P(s^{},\vartheta _k^{})P(\{\eta _i\}|\vartheta _k^{},s^{})\right]_{\epsilon ,\vartheta ^0}$$
(28)
where in analogy with eq.(1) we define

$$P(\{\eta _i\}|\vartheta _k,s)=\underset{i=1}{\overset{N}{}}\frac{1}{\sqrt{2\pi \sigma ^2}}exp\left[\left(\eta _i\epsilon _s^i\overline{\eta }_i(\vartheta _k)(1\epsilon _s^i)\eta _i^f\right)^2/2\sigma ^2\right];$$
(29)
and we discretize the average across the directional selectivities $`\{\vartheta _i^0\}`$ as well

$$𝑑\vartheta ^0\varrho (\vartheta ^0)\underset{k=1}{\overset{M}{}}\varrho (\vartheta _k^0).$$
(30)
This situation corresponds to the case when each neuron can discriminate across different angles $`\vartheta _1..\vartheta _M`$ with a resolution $`\mathrm{\Delta }\vartheta `$.

The calculation can be carried out introducing replicas as in the previous case. One gets

$$H(\{\eta _i\})=\frac{1}{\mathrm{ln}2}\underset{n0}{lim}\frac{1}{n}\left(\underset{s_1..s_{n+1}=1}{\overset{p}{}}\underset{k_1\mathrm{}k_{n+1}=1}{\overset{M}{}}\frac{(n+1)^{\frac{N}{2}}}{(Mp)^{n+1}}\frac{1}{(\sqrt{2\pi \sigma ^2})^{nN}}exp(R)^N_{\epsilon ,\vartheta ^0}1\right)$$
(31)
where

$$R=\left[\underset{l,m}{}\left(\epsilon _{s_l}(\overline{\eta }(\vartheta _{k_l})\eta ^f)\epsilon _{s_m}(\overline{\eta }(\vartheta _{k_m})\eta ^f)\right)^2/(4\sigma ^2(n+1))\right];$$
(32)
and we have assumed symmetry across neurons in the quenched randomness and in the parameters characterizing the conditional distribution $`P(\eta _i|s,\vartheta _k)`$.

Now we take the limit $`N\mathrm{}`$. As it is evident from eq.(32), $`exp(R)1`$ and $`exp(R)=1`$ when $`s_l=s_m`$ and $`k_m=k_l`$ for each pair of indexes ($`m`$,$`l`$). Thus when $`N\mathrm{}`$ the only terms which survive in the sum on replicas are the ones with $`s_1=s_2..=s_{n+1}`$ and $`k_1=k_2..=k_{n+1}`$. Since we have $`p`$ stimuli $`s`$ and $`M`$ stimuli $`\vartheta _k`$ the total number of terms is $`Mp`$. Substituting this value in the sums over replicas in eq.(31) and putting $`exp(R)=1`$ one obtains an expression for the entropy of the responses $`H(\{\eta _i\})`$, which summed to the equivocation as in eq.(8) gives the final result for the mutual information:

$$I(\{\eta _i\},\vartheta s)=\mathrm{log}_2(p)+\mathrm{log}_2(M);$$
(33)
Now we remember that $`M=2\pi /\mathrm{\Delta }\vartheta `$. Taking the limit to continuous angles, $`\mathrm{\Delta }\vartheta 0`$, it is easy to see that asymptotically the mutual information tends logarithmically to infinity.

## Beyond the gaussian assumption: the tg model

So far we have considered the case where the rate distribution for each neuron is normal. This assumption implies that negative rates have a non zero probability to occur; the more the average rate is small and close to zero, the more this probability becomes large. The bias introduced by the inclusion of negative rates in the space of possible states might be even more serious since we have considered the limit of large noise, where the tail of the distribution in the domain $`\eta <0`$ acquires a significant weight.

Cutting the distribution at zero is not enough to assign the proper weight to under-threshold activity: each time the summation of the inputs coming from other units is lower than threshold the neuron remains silent, and this occurs with a well defined probability.

A natural choice for the rate distribution $`P(\eta _i|\vartheta ,s)`$ is a thresholded gaussian plus a $`\delta `$ peak in zero (tg model):

$$P(\eta _i|\vartheta ,s)=\frac{1}{\sqrt{2\pi \sigma ^2}}exp[(\eta _i\stackrel{~}{\eta }_i(\vartheta ,s))^2/2\sigma ^2]\mathrm{\Theta }(\eta _i)+2(1\mathrm{erf}(\stackrel{~}{\eta }_i(\vartheta ,s)/\sigma )\delta (\eta _i)\mathrm{\Theta }(\eta _i)$$
(34)
where $`\mathrm{\Theta }(x)`$ is the Heaviside step function, $`\stackrel{~}{\eta }_i(\vartheta ,s)`$ is the same as defined in eq.(2) and $`\mathrm{erf}(x)`$ is the error function:

$$\mathrm{erf}(x)=\frac{1}{\sqrt{2\pi }}_{\mathrm{}}^x𝑑te^{t^2/2}.$$
(35)
The factor multiplying the $`\delta `$ function ensures a correct normalization and it assigns the proper weight to the peak in zero, which is larger the more the average rate $`\stackrel{~}{\eta }_i(\vartheta ,s)`$ is close to zero. A similar distribution has already been considered in networks of threshold linear neurons .

The analytical evaluation of the mutual information is obviously more difficult than in the case of the simple gaussian, because of the presence of the error function, which cannot be integrated exactly. Nonetheless in the limit of large $`\sigma `$ it is possible to evaluate both the linear and the quadratic approximation in $`N`$ and thus quantify the impact of the correction with respect to the gaussian case, eq.(24).

### The limit of large $`\sigma `$ for the tg model: the equivocation

We remind the expression of the equivocation, eq.(6):

$$H(\{\eta _i\}|\vartheta ,s)_{\vartheta ,s}=\underset{s=1}{\overset{p}{}}𝑑\vartheta \underset{i}{}d\eta _iP(\vartheta ,s)P(\{\eta _i\}|\vartheta ,s)\mathrm{log}_2P(\{\eta _i\}|\vartheta ,s)_{\epsilon ,\vartheta ^0}.$$
(36)
Assuming independence among neurons in the conditional probability $`P(\{\eta _i\}|\vartheta ,s)`$, eq.(36) can be written

$$H(\{\eta _i\}|\vartheta ,s)_{\vartheta ,s}=\frac{N}{\mathrm{ln}2}\underset{s=1}{\overset{p}{}}𝑑\vartheta P(\vartheta ,s)𝑑\eta P(\eta |\vartheta ,s)\mathrm{ln}P(\eta |\vartheta ,s)_{\epsilon ,\vartheta ^0}.$$
(37)
In the specific case of the distribution (34) it is easy to show that

$`{\displaystyle 𝑑\eta P(\eta |\vartheta ,s)\mathrm{ln}P(\eta |\vartheta ,s)}_{\epsilon ,\vartheta ^0}=`$
$`{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \frac{\stackrel{~}{\eta }(\vartheta ,s)}{2\sigma }}e^{\left[\stackrel{~}{\eta }(\vartheta ,s)\right]^2/2\sigma ^2}_{\epsilon ,\vartheta ^0}{\displaystyle \frac{1}{2}}(1+\mathrm{ln}(2\pi \sigma ^2))\mathrm{erf}(\stackrel{~}{\eta }(\vartheta ,s)/\sigma )_{\epsilon ,\vartheta ^0}`$
$`+{\displaystyle _0^{\mathrm{}}}𝑑\eta \delta (\eta )\left(1\mathrm{erf}(\stackrel{~}{\eta }(\vartheta ,s)/\sigma )\right)\left[\mathrm{ln}\delta (\eta )+\mathrm{ln}2+\mathrm{ln}\left(1\mathrm{erf}(\stackrel{~}{\eta }(\vartheta ,s)/\sigma )\right)\right]_{\epsilon ,\vartheta ^0}.`$ (38)
To proceed with the calculation we have to be careful with the integration of the delta function. In fact it is easy to show that the integration of the product $`\delta (x)\mathrm{ln}\delta (x)`$ yields a logarithmic divergence. Since the mutual information must remain finite with a finite number of neurons $`N`$, we expect this divergence to cancel exactly with an analogous term in the rate entropy and, in fact, in the next section, we will show that this is the case. For the moment we use the equality

$$_{\mathrm{}}^+\mathrm{}𝑑x\delta (x)F(x)=\underset{ϵ0}{lim}_{ϵ/2}^{ϵ/2}𝑑x\frac{1}{ϵ}F(x).$$
(39)
Assuming as usual that the quenched disorder is identically distributed across neurons and stimuli and that $`\stackrel{~}{\eta }(\vartheta ,s)`$ is like in eq.(1) we can write:

$`H(\{\eta _i\}|\vartheta ,s)_{\vartheta ,s}={\displaystyle \frac{N}{2\mathrm{ln}2}}\{(1+\mathrm{ln}(2\pi \sigma ^2))\mathrm{erf}(\stackrel{~}{\eta }(\vartheta ,s)/\sigma )_{\epsilon ,\vartheta ^0}{\displaystyle \frac{\stackrel{~}{\eta }(\vartheta ,s)}{\sqrt{2\pi }\sigma }}e^{\left[\stackrel{~}{\eta }(\vartheta ,s)\right]^2/2\sigma ^2}_{\epsilon ,\vartheta ^0}`$
$`+2[1\mathrm{erf}(\stackrel{~}{\eta }(\vartheta ,s)/\sigma )]_{\epsilon ,\vartheta ^0}\mathrm{ln}{\displaystyle \frac{ϵ}{2}}2[1\mathrm{erf}(\stackrel{~}{\eta }(\vartheta ,s)/\sigma )]\mathrm{ln}[1\mathrm{erf}(\stackrel{~}{\eta }(\vartheta ,s)/\sigma )]_{\epsilon ,\vartheta ^0}\}.`$ (40)
The average across quenched disorder cannot be performed if we do not resort to some approximation. Since we have already focused on the limit of large $`\sigma `$ it is natural to consider an expansion of the error function in eq. (35) for a small value of its argument:

$$\mathrm{erf}(x)\frac{1}{2}+\frac{1}{\sqrt{2\pi }}x+o(x^2);$$
(41)
Approximating all the error functions in eq.(40) we obtain:

$`H(\{\eta _i\}|\vartheta ,s)_{\vartheta ,s}={\displaystyle \frac{N}{2\mathrm{ln}2}}\{{\displaystyle \frac{1}{2}}(1+\mathrm{ln}(2\pi \sigma ^2))+\mathrm{ln}(ϵ)`$
$`+{\displaystyle \frac{\stackrel{~}{\eta }(\vartheta ,s)}{\sqrt{2\pi }\sigma }}_{\epsilon ,\vartheta ^0}(2+\mathrm{ln}(2\pi \sigma ^2)2\mathrm{ln}ϵ){\displaystyle \frac{1}{\pi }}{\displaystyle \frac{\stackrel{~}{\eta }(\vartheta ,s)^2}{\sigma ^2}}_{\epsilon ,\vartheta ^0}\}+o(1/\sigma ^2)`$ (42)
where in line with the approximation used in the case of the simple gaussian distribution we have omitted terms of order $`N/\sigma ^k`$ with $`k>2`$.

### Evaluation of the mutual information

We reconsider eq.(7). Using replicas and assuming that the quenched randomness is identically distributed across neurons we obtain

$$H(\{\eta _i\})=\frac{1}{\mathrm{ln}2}\underset{n0}{lim}\frac{1}{n}(\underset{s_1..s_{n+1}=1}{\overset{p}{}}d\vartheta _1..d\vartheta _{n+1}\frac{1}{(2\pi p)^{n+1}}d\eta \underset{k=1}{\overset{n+1}{}}P(\eta |\vartheta _k,s_k)_{\epsilon ,\vartheta ^0}^N1),$$
(43)
where $`P(\eta |\vartheta ,s)`$ is given in eq.(34). Integrating over $`d\eta `$ yields

$`{\displaystyle }d\eta {\displaystyle \underset{k=1}{\overset{n+1}{}}}P(\eta |\vartheta _k,s_k)_{\epsilon ,\vartheta ^0}^N=(n+1)^{\frac{1}{2}}{\displaystyle \frac{1}{(\sqrt{2\pi \sigma ^2})^n}}exp(R)\mathrm{erf}({\displaystyle \frac{1}{\sqrt{n+1}}}{\displaystyle \underset{k}{}}\stackrel{~}{\eta }(\vartheta _k,s_k)/\sigma )`$
$`+\left({\displaystyle \frac{2}{ϵ}}\right)^n{\displaystyle \underset{k}{}}[1\mathrm{erf}\left({\displaystyle \frac{\stackrel{~}{\eta }(\vartheta _k,s_k)}{\sigma }}\right)]_{\epsilon ,\vartheta ^0}^N`$ (44)
where we have used eq.(39) to integrate the $`\delta `$ function in eq.(34), and the expression of $`R`$ is like in eq.(11). Using the approximation (41) for the error function and considering the expansion for small $`n`$

$$a^n1+n\mathrm{ln}a$$
(45)
we obtain

$$H(\{\eta _i\})=\frac{1}{\mathrm{ln}2}\underset{n0}{lim}\frac{1}{n}(\underset{s_1..s_{n+1}=1}{\overset{p}{}}d\vartheta _1..d\vartheta _{n+1}\frac{1}{(2\pi p)^{n+1}}(1nC\underset{k,lk}{}G_{kl})^N1);$$
(46)
where

$$C=\frac{1}{2}\left[\frac{1}{2}\left(1+\mathrm{ln}(2\pi \sigma ^2)\right)+\mathrm{ln}(ϵ)+\underset{k}{}\frac{\stackrel{~}{\eta }(\vartheta _k,s_k)}{\sqrt{2\pi }\sigma }_{\epsilon ,\vartheta ^0}\left(2+\mathrm{ln}(2\pi \sigma ^2)2\mathrm{ln}ϵ\right)\right];$$
(47)
$$G_{kl}=\frac{1}{8(n+1)\sigma ^2}\left[\left(\stackrel{~}{\eta }(\vartheta _k,s_k)\stackrel{~}{\eta }(\vartheta _l,s_l)\right)^2_{\epsilon ,\vartheta ^0}\frac{1}{2\pi }\stackrel{~}{\eta }(\vartheta _k,s_k)\stackrel{~}{\eta }(\vartheta _l,s_l)_{\epsilon ,\vartheta ^0}\right].$$
(48)
It is simple to verify that $`C`$ remains finite when $`n`$ goes to zero. Now we expand in powers of $`N`$ up to the second order:

$$\left(1nC\underset{k,lk}{}G_{kl}\right)^N1N\left(nC+\underset{k,lk}{}G_{kl}\right)+\frac{N^2}{2}\underset{k,lk}{}\underset{\varrho ,\mu \varrho }{}G_{kl}G_{\varrho \mu };$$
(49)
where we have omitted terms that are $`o(n)`$ when $`n0`$. This quantity has to be summed over continuous and discrete replicas, after having explicitly performed the average across quenched disorder in eqs.(47) and (48). It is easy to show that $`C`$ cancels exactly with analogous terms in the equivocation, eq.(42). The evaluation of the linear and quadratic term in the mutual information can be performed with a similar technique to the one used in the case of the gaussian distribution and it involves averaging terms with $`2`$ and $`4`$replica interactions.

The final expression for the mutual information in the linear and quadratic approximation reads

$`I(\{\eta _i\},\vartheta s)`$
$`{\displaystyle \frac{1}{\mathrm{ln}2}}\{{\displaystyle \frac{N(\eta ^0)^2}{4\sigma ^2}}({\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{\pi }})[{\displaystyle \frac{p1}{p}}2(\alpha A_1)^2\lambda _1+2(A_2(A_1)^2)\lambda _2]{\displaystyle \frac{N^2(\eta ^0)^4}{2(4\sigma ^2)^2}}({\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{\pi }})^2`$ (50)
$`\{{\displaystyle \frac{p1}{p^2}}2\left(2(\alpha A_1)^2\lambda _1\right)^2+[{\displaystyle \frac{p1}{p}}(\lambda _1\lambda _2)^2+{\displaystyle \frac{(\lambda _2)^2}{p}}]\left[\left({\displaystyle \frac{1}{2^{2m1}}}\right)^4{\displaystyle \underset{\nu =0}{\overset{m1}{}}}\left[\left(\begin{array}{c}2m\\ \nu \end{array}\right)\right]^4\right]\}\}.`$ (53)
Comparing eq.(24) and (53) it is evident that modifying the gaussian model into the more realistic tg model has no effect on the analytical expression of the mutual information, except for a renormalization of the expansion parameter $`g^2=(\eta ^0/2\sigma )^2`$:

$$g^{}=g\sqrt{\frac{1}{2}+\frac{1}{\pi }}.$$
(54)
Fig.5 shows the effect of the renormalization for different values of $`g^2`$. The mutual information is lower in the tg model than in the gaussian approximation, as expected.

We have explored whether the two models can fit the information rise estimated from real data. Since the analytical expression of the mutual information is the same in both cases the fit does not change between the two models.

Fig.6 shows the comparison between the information estimated from a sample of cells recorded in the right primary motor cortex  and the prediction given by either of the two theoretical distributions. In the limit of large $`\sigma `$ the gaussian model fails to provide a good fit, but we can conclude that this failure is not due to the inclusion of negative rates in the distribution.

## Information loss in averaging activity distributions

Fig.1 suggests that the directional tuning of real cells is modulated, albeit moderately, by the type of movement. In fact the analysis of real data has proved that the coding of the direction is not unique, but it is specific to the complex correlate which is being considered  (here, which arm moves).

More in general, distinct features characterizing a complex stimulus are not expected to be coded independently of one another. This raises the question of how central representations of external correlates are constructed and which are the basic featural components of these representations. Of course the categorization of natural stimuli is arbitrary and the more accurate a description is provided, the higher the dimensionality of the stimulus set. Since an infinite number of different descriptors could be chosen to characterize a stimulus, any (finite) categorization has the effect of emphasizing some relevant features and averaging out other irrelevant features. An obvious consequence is that we end up, even involuntarily, evaluating how some features are coded on average, with respect to the dimension we have chosen, explicitly or implicitly, to neglect.

Thus, with correlates which have one continuous and one discrete dimension, one might wonder which are the relationships among the information carried about the total number of continuous+discrete dimensions, the information carried about the continuous dimension, disregarding the discrete dimension and, finally, the information carried about the continuous dimension, if a single value of the discrete dimension had been fixed when recording neural activity. In other words, suppose that we investigate how the direction is coded on average across different types of movement. This corresponds to averaging the full distribution $`P(\{\eta _i\}|\vartheta ,s)`$ on $`s`$:

$$P(\{\eta _i\}|\vartheta ,s)\underset{s}{}P(s)P(\{\eta _i\}|\vartheta ,s)=P_s(\{\eta _i\}|\vartheta ).$$
(55)
The resulting expression of the mutual information is:

$$I(\{\eta _i\},\vartheta _s)=𝑑\vartheta \underset{i}{}d\eta _iP(\vartheta )P_s(\{\eta _i\}|\vartheta )\mathrm{log}_2\frac{P_s(\{\eta _i\}|\vartheta )}{P(\{\eta _i\})}_{\epsilon ,\vartheta ^0}.$$
(56)
The analytical evaluation is very similar to the cases already discussed.

As usual, the mutual information can be expressed as the difference between the entropy of the responses $`H(\{\eta _i\})`$ and the equivocation $`H(\{\eta _i\}|\vartheta )_\vartheta `$, where

$$H(\{\eta _i\}|\vartheta )_\vartheta =𝑑\vartheta \underset{i}{}d\eta _iP(\vartheta )P_s(\{\eta _i\}|\vartheta )\mathrm{log}_2P_s(\{\eta _i\}|\vartheta )_{\epsilon ,\vartheta ^0};$$
(57)
$$H(\{\eta _i\})=𝑑\vartheta \underset{i}{}d\eta _iP(\vartheta )P_s(\{\eta _i\}|\vartheta )\mathrm{log}_2\left[𝑑\vartheta ^{}P(\vartheta ^{})P_s(\{\eta _i\}|\vartheta ^{})\right]_{\epsilon ,\vartheta ^0}.$$
(58)
We focus on the more realistic tg model, eq.(34); the entropy of the responses is obviously independent on the chosen categorization of the stimuli; in the limit of large $`\sigma `$ the procedure is precisely the one followed in the previous section.

The difference with respect to the cases already discussed is in the evaluation of the equivocation: since $`P_s(\{\eta _i\}|\vartheta )`$ is obtained averaging the distribution (34) across $`s`$, we need to introduce discrete replicas $`s_1..s_n`$, as when evaluating the entropy of the responses. Then the calculation is straightforward and the basic steps are given in the previous section and in the appendix.

Since in one case (for the entropy of the responses) we sum both over discrete and over continuous replicas and in the other case (for the equivocation) we sum only on discrete replicas, it is clear that all terms that do not involve 2 or more replica interacting cancel out.

More in detail, if the evaluation of the entropy of the responses requires the analytical calculation of averages like $`\overline{\eta }(\vartheta _k)\overline{\eta }(\vartheta _l)_{\vartheta ^0}`$ (see the appendix), these averages disappear in the evaluation of the equivocation, since replica indexes are only for the discrete variable $`s`$.

The final result for the mutual information up to the quadratic approximation reads:

$`I(\{\eta _i\},\vartheta _s){\displaystyle \frac{1}{\mathrm{ln}2}}\{{\displaystyle \frac{N(\eta ^0)^2}{4\sigma ^2}}(1+{\displaystyle \frac{2}{\pi }})(A_2(A_1)^2)[\lambda _2{\displaystyle \frac{p1}{p}}\lambda _1]`$
$`{\displaystyle \frac{N^2(\eta ^0)^4}{2(4\sigma ^2)^2}}({\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{\pi }})^2\{{\displaystyle \frac{p1}{p^2}}8[(\alpha A_1)^4(A_22\alpha A_1+\alpha ^2)^2]\left(\lambda _1\right)^2`$
$`+[{\displaystyle \frac{p1}{p}}(\lambda _1\lambda _2)^2+{\displaystyle \frac{(\lambda _2)^2}{p}}]\left[\left({\displaystyle \frac{1}{2^{2m1}}}\right)^4{\displaystyle \underset{\nu =0}{\overset{m1}{}}}\left[\left(\begin{array}{c}2m\\ \nu \end{array}\right)\right]^4\right]\}\}.`$ (61)
Fig.7(a) compares the averaged information, $`I(\{\eta _i\},\vartheta _s)`$, with the full information $`I(\{\eta _i\},\vartheta s)`$, where in the case of $`I(\{\eta _i\},\vartheta s)`$ we have put $`p=1`$. The full information calculated with $`p=1`$ in fact gives the curve one would obtain, on average, by considering only one movement type (or value of the discrete correlate) at a time. As one could expect, averaging the distribution across the discrete correlates results in an information loss. Moreover, the the full information with $`p=4`$ movement types is obviously above the specific one, obtained by setting $`p=1`$ (compare fig.5 and 7(a)).

Fig.7(b) shows the dependence of the full and averaged information on the number $`p`$ of discrete correlates. Contrary to the full information, the averaged information decreases monotonically with $`p`$, both in the linear and in quadratic approximation.

As one would expect, averaging the distribution across a large number $`p`$ of correlates is equivalent to a regularization of the activity distribution, which results in a lower mutual information.

## Discussion

We have studied a model of the coding of discrete+continuous stimuli by a population of $`N`$ neurons, referring to the specific case of movements categorized according to a direction and a type. We have shown that, asymptotically in the limit of large populations of neurons, the mutual information tends to infinity logarithmically with the resolution in the continuous dimension. This result aside, we have focused on the initial rise of the mutual information with the population size, which may offer a more direct comparison with the analysis of real data.

In the limit of large noise we have therefore derived an analytical formula for such initial rise of the information, up to the quadratic approximation. We have examined the dependence on the number of discrete correlates and on the width of the directional tuning. A comparison with the information estimated from real data, in which the linear term is used as a fit coefficient, has shown that the quadratic approximation fails to capture the deviation from linearity.

We have then considered a more realistic model for the conditional firing distribution, a thresholded gaussian with a $`\delta `$-peak. We have shown that this more realistic distribution simply renormalizes the expansion parameter applicable to the gaussian distribution. Therefore, the discrepancy in the fit to real data does not originate in the firing rate distribution. There are several possible reasons for this discrepancy:

* the value of the expansion parameter $`g^2=(\eta ^0/\sigma )^2`$ corresponding to the best fit in fig.(6) is quite high ($`g^2=0.78`$). This value is in the range where the quadratic approximation is expected to fail on its own (see fig.(3),(5)). Adding higher orders in perturbation theory might improve the fit. Moreover, we have neglected terms of order $`N/\sigma ^k`$ with $`k>2`$, which might be non negligible when $`g^2`$ becomes large and $`N`$ is not too small.
* information estimates from real data are often biased because of poor sampling. Several procedures have been proposed to correct the bias , but the improvement given by the correction is not precisely quantifiable, and sampling biases cannot be discarded for good.
* both with the tg model and in the gaussian approximation, we have assumed that neurons fire independently of one another to each movement, but the analysis of extracellular recordings has shown that correlations may play a non-negligible role in the coding .

Finally, we have examined the effect of averaging the distribution across the discrete correlates, evaluating the mutual information with respect to the continuously varying dimension alone. As expected, averaging the distribution across $`s`$ results in an information loss, which is more serious the larger the number $`p`$ of discrete correlates.

Further developments of this work include: the introduction of weak correlations in the signal of different neurons and the analysis of the information transfer to other stages of processing, given this as the coding scheme in the input layer of the network. In the specific case of movement coding these research directions might help model how information is transmitted down the motor system, in the planning and execution of motor tasks.

## Appendix A Detailed evaluation of the second order coefficient

We show here how to evaluate the coefficient of the quadratic term in eq.(12), that we write again:

$$\underset{k,lk}{}\underset{\varrho ,\mu \varrho }{}\left(\epsilon _{s_k}(\overline{\eta }(\vartheta _k)\eta ^f)\epsilon _{s_l}(\overline{\eta }(\vartheta _l)\eta ^f)\right)^2_{\epsilon ,\vartheta ^0}\left(\epsilon _{s_\varrho }(\overline{\eta }(\vartheta _\varrho )\eta ^f)\epsilon _{s_\mu }(\overline{\eta }(\vartheta _\mu )\eta ^f)\right)^2_{\epsilon ,\vartheta ^0}.$$
(62)
This quantity has to be integrated over continuous and discrete replicas. First we perform the average across the quenched variable $`\vartheta _0`$. Expanding the products it is easy to see that the quantities to be averaged are all like $`\overline{\eta }(\vartheta _k)_{\vartheta ^0}`$ and $`\overline{\eta }^2(\vartheta _k)_{\vartheta ^0}`$, that we have already calculated in eq.(19),(20). Another average that we need is $`\overline{\eta }(\vartheta _k)\overline{\eta }(\vartheta _l)_{\vartheta ^0}`$, with $`kl`$:

$$\overline{\eta }(\vartheta _k)\overline{\eta }(\vartheta _l)_{\vartheta ^0}=\left(\eta ^0\right)^2\left\{\left(\frac{1}{2^{2m}}\right)^2\left[\left(\begin{array}{c}2m\\ m\end{array}\right)\right]^2+\frac{1}{2}\left(\frac{1}{2^{2m1}}\right)^2\underset{\nu =0}{\overset{m1}{}}\left[\left(\begin{array}{c}2m\\ \nu \end{array}\right)\right]^2\mathrm{cos}\left\{\left(m\nu \right)\left(\vartheta _k\vartheta _l\right)\right\}\right\}.$$
(63)
Thus only terms like $`\overline{\eta }(\vartheta _k)\overline{\eta }(\vartheta _l)_{\vartheta ^0}`$ still depend on continuous replicas.

Since terms like $`\mathrm{cos}(m\nu )(\vartheta _k\vartheta _l)`$, with $`kl`$ are zero when integrated on $`\vartheta _k`$,$`\vartheta _l`$ the only term which requires a careful evaluation in performing the integration on continuous replicas is the product $`\overline{\eta }(\vartheta _k)\overline{\eta }(\vartheta _l)_{\vartheta ^0}\overline{\eta }(\vartheta _\varrho )\overline{\eta }(\vartheta _\nu )_{\vartheta ^0}`$.

After integration on continuous replicas this term yields

$$\left(\eta ^0\right)^4\left\{\left(\frac{1}{2^{2m}}\right)^4\left[\left(\begin{array}{c}2m\\ m\end{array}\right)\right]^4+\frac{1}{4}\left(\frac{1}{2^{2m1}}\right)^4\underset{\nu =0}{\overset{m1}{}}\left[\left(\begin{array}{c}2m\\ \nu \end{array}\right)\right]^4\left(\delta _{k\varrho }\delta _{\mu l}+\delta _{k\mu }\delta _{\varrho l}\right)\right\}.$$

Rearranging all terms one has, finally

$`{\displaystyle \underset{k,lk}{}}{\displaystyle \underset{\varrho ,\mu \varrho }{}}\left[\left(\alpha A_1\right)^2(\epsilon _{s_k}\epsilon _{s_l})^2_\epsilon +\left(A_2(A_1)^2\right)\epsilon _{s_k}^2+\epsilon _{s_l}^2_\epsilon \right]x`$
$`x\left[\left(\alpha A_1\right)^2(\epsilon _{s_\varrho }\epsilon _{s_\mu })^2_\epsilon +\left(A_2(A_1)^2\right)\epsilon _{s_\varrho }^2+\epsilon _{s_\mu }^2_\epsilon \right]`$
$`+\epsilon _{s_k}\epsilon _{s_l}_\epsilon \epsilon _{s_\varrho }\epsilon _{s_\mu }_\epsilon \left({\displaystyle \frac{1}{2^{2m1}}}\right)^4{\displaystyle \underset{\nu =0}{\overset{m1}{}}}\left[\left(\begin{array}{c}2m\\ \nu \end{array}\right)\right]^4\left(\delta _{k\varrho }\delta _{\mu l}+\delta _{k\mu }\delta _{\varrho l}\right)`$ (66)
where $`A_1`$ and $`A_2`$ are defined in eq.(19) and (20). To correctly perform the summation across the discrete replicas and to keep only terms of order $`n`$ in the limit $`n0`$, we consider separately the different contributions to the sum:

$$\underset{k,lk}{}\underset{\varrho ,\mu \varrho }{}\underset{kl\varrho \mu }{}+\underset{\varrho ,kl\mu }{}\left(\delta _{\varrho k}+\delta _{\varrho l}\right)+\underset{\mu ,kl\varrho }{}\left(\delta _{\mu k}+\delta _{\mu l}\right)+\underset{\mu ,\varrho ,kl}{}\left(\delta _{\mu k}\delta _{\varrho l}+\delta _{\varrho k}\delta _{\mu l}\right).$$
(67)
Since in eq.(66) terms like $`(\epsilon _{s_k}\epsilon _{s_l})^2_\epsilon `$ are zero if $`s_k=s_l`$, we have to distinguish between different cases for each term in eq.(67) in performing the summation on discrete replicas:

$`{\displaystyle \underset{s_1}{}}..{\displaystyle \underset{s_{n+1}}{}}{\displaystyle \underset{s_1}{}}..{\displaystyle \underset{s_ks_l}{}}{\displaystyle \underset{s_\varrho s_m}{}}..{\displaystyle \underset{s_{n+1}}{}}+{\displaystyle \underset{s_1}{}}..{\displaystyle \underset{s_k=s_l}{}}{\displaystyle \underset{s_\varrho s_m}{}}..{\displaystyle \underset{s_{n+1}}{}}`$
$`+{\displaystyle \underset{s_1}{}}..{\displaystyle \underset{s_ks_l}{}}{\displaystyle \underset{s_\varrho =s_m}{}}..{\displaystyle \underset{s_{n+1}}{}}+{\displaystyle \underset{s_1}{}}..{\displaystyle \underset{s_k=s_l}{}}{\displaystyle \underset{s_\varrho =s_m}{}}..{\displaystyle \underset{s_{n+1}}{}};`$ (68)
With a bit of combinatorics and keeping only the terms order $`n`$, which give a finite contribution in the limit $`n0`$, we obtain the final result for the second order term

$$\frac{N^2(\eta ^0)^4}{2(4\sigma ^2)^2}\left\{\frac{p1}{p^2}2\left(2\left(\alpha A_1\right)^2\lambda _1\right)^2+\left[\frac{p1}{p}\left(\lambda _1\lambda _2\right)^2+\frac{\lambda _2^2}{p}\right]\left[\left(\frac{1}{2^{2m1}}\right)^4\underset{\nu =0}{\overset{m1}{}}\left[\left(\begin{array}{c}2m\\ \nu \end{array}\right)\right]^4\right]\right\};$$

where the expressions of $`\lambda _1`$,$`\lambda _2`$ are given in eqs.(15), (16). Considering the result obtained at the first order, eq.(18), it is easy to derive the expression for the mutual information up to the second order in $`N/\sigma ^2`$:

$`I(\{\eta _i\},\vartheta s){\displaystyle \frac{1}{\mathrm{ln}2}}\{{\displaystyle \frac{N(\eta ^0)^2}{4\sigma ^2}}[{\displaystyle \frac{p1}{p}}2(\alpha A_1)^2\lambda _1+2(A_2(A_1)^2)\lambda _2]`$
$`{\displaystyle \frac{N^2(\eta ^0)^4}{2(4\sigma ^2)^2}}\{{\displaystyle \frac{p1}{p^2}}2\left(2(\alpha A_1)^2\lambda _1\right)^2+[{\displaystyle \frac{p1}{p}}(\lambda _1\lambda _2)^2+{\displaystyle \frac{\lambda _2^2}{p}}]\left[\left({\displaystyle \frac{1}{2^{2m1}}}\right)^4{\displaystyle \underset{\nu =0}{\overset{m1}{}}}\left[\left(\begin{array}{c}2m\\ \nu \end{array}\right)\right]^4\right]\}\}.`$
### Acknowledgements

We have enjoyed extensive discussions with Inés Samengo, Eilon Vaadia, Edmund Rolls and Elka Korutcheva. Partial support from Human Frontier Science Programme grant RG 0110/1998-B.