Mathematical properties of Gaussian process neural code

These notes develop the mathematical properties of Gaussian processes that matter when they are used as models of neural representation. The motivating example is a population code for space. A smooth random field supplies the subthreshold response of each neuron, and a cellular threshold turns its excursions into place fields. This point of view explains much of the field structure seen across species and environments (Mainali, da Silveira & Burak, 2024). It also gives a tractable way to ask what such a code can represent and how accurately it can be read.

The note began as a chalk talk at Janelia on 18 November 2024. We begin with Gaussian variables and processes, pass through the geometry of threshold excursions, and finish with a simple calculation of local and global decoding error.

§Background

Probability enters theoretical neuroscience in several distinct ways. It can represent noise in neural recordings, pose a learning problem with stochastic data or plasticity, model random synaptic weights, or provide an analytical route into large statistical systems. Here randomness serves a different purpose. It models a neural representation itself as a random code over a continuous external variable such as position.

§Gaussian variables

§Conditional Gaussian

The most useful closure property of a multivariate Gaussian is that conditioning one subset of variables on another still gives a Gaussian. Let $\boldsymbol\xi$ and $\boldsymbol\eta$ be jointly Gaussian, with means $\mathbf m_\xi$ and $\mathbf m_\eta$, and covariance

§Gaussian processes

Moving from a Gaussian variable to a Gaussian vector replaces one value with a finite indexed collection. A Gaussian process continues that move to a continuously indexed random function. It is completely determined by its mean and covariance functions,

For a stationary process with constant mean, the covariance depends only on the separation $\Delta x=y-x$, so we write $r(x,y)=r(\Delta x)$. In particular, $r(0)=r_0$ is the pointwise variance.

§Origins of a random code

Consider $N$ spatially selective neurons presynaptic to a CA1 place cell. Let $u_i(x)$ be the bounded response of input neuron $i$ over an environment $x\in[0,L]$. For clarity, suppose all inputs share one tuning shape up to translation,

with preferred positions $x_i$ that tile the environment densely and uniformly. The subthreshold input to the CA1 cell is a random weighted sum,

where the independent weights $W_i$ have zero mean and finite variance. For any finite set of locations $\{x_k\}$, the multidimensional central limit theorem gives a jointly Gaussian vector $(f(x_1),\ldots,f(x_k))$ as $N$ grows (van der Vaart, 1998; Lehmann & Casella, 1998). The limiting random function is therefore a Gaussian process. The central limit theorem supplies the Gaussian finite-dimensional distributions. The tuning curves and uniform tiling assumptions supply smoothness and stationarity.

Uniform tiling also makes the covariance translation invariant. With input density $\rho=N/L$,

Choosing $\Var(W)=L/N$ leaves a covariance $r(\Delta x)$ that does not grow with population size.

Compare the random weighted sum with the cellular threshold that selects its connected excursions in Plate 1.1.

In this model, position $x$ is supplied by the environment. The process $f(x)$ represents subthreshold CA1 activity produced by many synaptic inputs. Cellular thresholding, together with any incident inhibition, reveals firing fields wherever $f(x)$ rises above a level of order its standard deviation.

§Excursion structure

The neural response can be idealized by a rectified threshold,

Its nonzero regions are the excursion sets of $f$ above $\theta\sqrt{r_0}$. Their number, length, separation, height, and boundary slopes expose much of the geometry of the code before any further output nonlinearity is chosen. Use the fixed path in Plate 2.1 to see how the threshold changes the crossings and connected excursions.

§The Kac-Rice formula

Kac and Rice gave a way to count level crossings of a smooth function, and then of a stochastic process, without locating each root directly (Rice, 1944; Kac & Slepian, 1959).

Compare shallow and steep crossings in Plate 2.2. Each contributes one to the count.

§Expected crossing count

We now apply the formula to a stationary, mean-zero, differentiable Gaussian process. At a fixed location, $f(x)$ and $f'(x)$ are jointly Gaussian. Stationarity gives

and

Thus

The zero covariance makes $f(x)$ and $f'(x)$ independent. Taking the expectation of (2.3) at level $T=\theta\sqrt{r_0}$ gives

Every bounded excursion has one upcrossing and one downcrossing. Its expected density is therefore half the crossing density,

Let $\Phi$ be the standard normal cumulative distribution function. The expected fraction of the domain above the threshold is $1-\Phi(\theta)$. Dividing the total occupied and unoccupied lengths by the expected number of excursions gives the approximate mean excursion length $s$ and mean gap $\bar s$,

At high thresholds, the gaps are much longer than the excursions. Consecutive events are then separated by several correlation lengths, and the crossing process approaches a Poisson process with the mean above.

§Spectral density

Define the even spectral moments

Then $\omega_0=r_0$ and $\omega_2=-r''_0$. The excursion density becomes

At the mean level, $\mu(0)=(2\pi)^{-1}\sqrt{\omega_2/\omega_0}$. It has units of inverse length and supplies a geometric correlation scale.

Repeating the argument for the derivative process gives equal densities of local maxima and minima,

The ratio between turning points and mean-level crossings defines the scale-invariant irregularity parameter

Values near one indicate a regular process with roughly one maximum and one minimum between mean-level upcrossings. Small values indicate additional high-frequency structure.

Compare the normalized densities for $\omega_0$, $\omega_2$, and $\omega_4$ in Plate 2.3. Dividing each weighted spectrum by its own moment gives unit area, so the three shapes share one vertical scale even though their unnormalized moments have different units.

§Excursion shape

Counting fields is only the first step. To describe a typical field, we condition on the event that the process crosses the threshold from below (Kac & Slepian, 1959). Palm conditioning selects the observation location from the crossings of the process. A location sampled uniformly in space follows a different distribution.

§Slope at an upcrossing

Although $f'$ is Gaussian at a uniformly chosen point, it is not Gaussian at an upcrossing. Steep slopes create more crossings per unit length and are sampled more often. The derivative factor in the Kac-Rice formula expresses exactly this selection bias. From this point onward, $\theta$ denotes the threshold itself. Replace it with $\theta\sqrt{r_0}$ to restore the variance scale.

§The Slepian process

For a Gaussian process, the conditional law inside this integral is Gaussian. Its mean and covariance therefore determine it completely. Take

Bochner’s representation gives the remaining cross-covariance,

The joint covariance is therefore

With $\omega_0=r_0$ and $\omega_2=-r''(0)$, Gaussian conditioning gives

and residual covariance

The last two terms are the reduction in uncertainty caused by fixing the height and slope at the crossing. They vanish far from the origin when the covariance and its derivative decay.

Compare the conditional mean and residual process with the high-threshold parabolic limit in Plate 2.4.

§High excursions

At high threshold, the length and excess height of an excursion are both of order $\theta^{-1}$. Expand the covariance around the crossing:

The conditioned residual satisfies $\kappa(x/\theta)=o(x/\theta)$. Substituting these terms into the Slepian process yields

The excursion is locally parabolic. Its excess height is approximately $\zeta^2\omega_0/(2\theta\omega_2)$ and its length is approximately $2\omega_0\zeta/(\theta\omega_2)$. The length inherits the Rayleigh law of the crossing slope, while the height inherits its square.

In an isotropic $D$-dimensional process, make the simplifying approximation that the $D$ principal lengths share this Rayleigh scaling. If a normalized linear size $S$ has density

and volume is $V=S^D$, then the change of variables $s=v^{1/D}$ gives

For $D=2$, this normalized volume, and likewise the squared Rayleigh height variable, is exponentially distributed.

§Coding properties

So far the Gaussian process has described the structure of a neural code. We now ask how well a population of these fields can encode position. The calculation follows a simplified random-code model without thresholding (Blanco Malerba et al., 2022).

Let $N$ neurons independently sample a stationary Gaussian process over a one-dimensional environment of length $L$,

Suppose the response at the true position $x=0$ is corrupted by independent Gaussian noise $z_i\sim\mathcal N(0,\eta^2)$. Under that noise model, maximum-likelihood decoding minimizes squared distance:

We measure performance by mean squared position error,

Two mechanisms contribute. A local error moves the single posterior peak away from the true position. A global error makes a distant, weakly correlated codeword appear closer than the true one. Their dependence on field roughness points in opposite directions.

§Local error

Suppose the error $\Delta x=\widehat x$ is small. Linearizing each field around zero gives

Conditioned on the field derivatives, averaging over the readout noise gives

The normalized sum is $\chi_N^2$. Since $\EE[(\chi_N^2)^{-1}]=1/(N-2)$ for $N>2$,

Rougher fields have larger $\omega_2$ and steeper local gradients, so they reduce this error.

§Global error

Now assume global errors are rare and occur at positions separated by more than one correlation length. Their count is then approximately Poisson. If $\lambda$ is the rate per independent segment,

At a distant position $x$, define the vector difference

A global error occurs when the distant noisy code is closer to the observation than the true code. The noise is symmetric, so changing the sign of the cross term leaves the error probability unchanged.

or equivalently when

Conditioned on $\widetilde{\mathbf f}$,

Dividing by the positive norm gives a Gaussian variable with mean $\norm{\widetilde{\mathbf f}}$ and variance $4\eta^2$. The conditional error probability can be written as

For points beyond the correlation length, each component of $\widetilde{\mathbf f}$ has variance approximately $2r_0$. Averaging the exponential over the resulting norm distribution, using the noncentral $\chi^2$ moment-generating function, gives the approximation used in the talk,

Completing the remaining half-Gaussian integral yields

To turn this single-comparison probability into an error rate, the note uses the simplifying estimate that an environment of length $L$ contains $L\omega_2$ effectively independent comparisons. This is a rough correlation-length approximation rather than an exact identity. It gives

Finally, approximate the signed displacement by a uniform variable on an interval of length $L$ centered at zero. Its second moment is $L^2/12$, so

Adding the local and global contributions gives

Compare the two error terms and their minimum in Plate 3.1 as field roughness changes.

The two terms assign opposite roles to $\omega_2$. More roughness improves local resolution but creates more independent opportunities for a distant confusion. Increasing $N$, by contrast, suppresses the global term exponentially (Blanco Malerba et al., 2022). The balance gives an optimal roughness within this simplified model.

§References