Linear dynamical systems and latent state

One neuron, one eigenvalue

Everything so far in this series has been static. We took a snapshot of neural activity at one moment, wrote it as a vector, and asked questions about its structure. But neural activity unfolds over time. The firing rate at time step $t$ depends on the firing rate at time step $t-1$. To say anything about how neural populations evolve, we need a model of dynamics. Start with the simplest possible case: a single neuron.

Let $x_t$ be the firing rate of one neuron at time step $t$. Suppose it evolves by a fixed rule: the rate at the next time step is some fraction $\lambda$ of the current rate.

The scalar $\lambda$ controls everything. Consider a neuron in prefrontal cortex during a working memory delay. The monkey saw the target, the target disappeared, and now the neuron holds its firing rate steady across several seconds of empty waiting. In the language of Equation 1, that is $\lambda \approx 1$. Each time step, the rate reproduces itself. The state persists.

Now consider a V1 neuron responding to a brief flash. The firing rate spikes and then decays back to baseline. That is $\lambda \approx 0.7$. Each time step, the rate shrinks to 70% of what it was. After ten steps, it is down to $0.7^{10} \approx 0.028$ of its initial value. The state fades.

If $\lambda > 1$, the rate explodes. Real neurons do not sustain this for long, but it matters for understanding stability. A network with even one unstable mode will, absent some nonlinearity to rein it in, blow up.

The eigenvalue controls the timescale. Define the time constant as $\tau = -1 / \log |\lambda|$. When $\lambda = 0.9$, $\tau \approx 9.5$ time steps: the rate falls to $1/e$ of its initial value in about 9.5 steps. When $\lambda = 0.99$, $\tau \approx 99.5$ time steps. As $\lambda$ approaches 1, the timescale diverges. The state persists indefinitely.*Working memory models in prefrontal cortex often use recurrent networks tuned so that the dominant eigenvalue sits near 1 [1, 2]. This gives persistent activity without external input. The eigenvalue is the knob that controls how long information is held. Continuous attractor models achieve $\lambda = 1$ exactly through fine-tuned symmetry in the weight matrix, while more robust implementations allow $\lambda$ slightly below 1 and periodically refresh the state.

This one-dimensional model is almost trivially simple. But every idea that matters for high-dimensional dynamics is already present. The eigenvalue $\lambda$ determines stability, timescale, and qualitative behavior. In the next section, those scalar eigenvalues become a spectrum, and the system becomes far richer.

Populations and the state equation

Scale up from one neuron to a population. Instead of a single firing rate $x_t$, we now have a $d$-dimensional latent state vector $\mathbf{x}_t$ that evolves according to a matrix:

The matrix $A$ is the rule that advances the state by one time step. It encodes the coupling between dimensions: how the activity along one direction at time $t$ influences every direction at time $t+1$. This is the state equation, the core of every linear dynamical system.

From Part 5, we know that if $A$ has a full set of eigenvectors, it factorizes as:

where $\Lambda$ is the diagonal matrix of eigenvalues and $V$ is the matrix of eigenvectors. This decomposition turns the coupled system into a set of independent modes. Define $\mathbf{z}_t = V^{-1} \mathbf{x}_t$, the state expressed in the eigenvector basis. Then $\mathbf{z}_{t+1} = \Lambda \, \mathbf{z}_t$, which is just $d$ copies of Equation 1 running in parallel. Each eigenvector defines a mode. Its eigenvalue controls whether that mode grows, decays, or oscillates, exactly as in the scalar case.

Real eigenvalues give modes that scale, one dimension at a time. A positive eigenvalue with $|\lambda| < 1$ decays monotonically. A negative eigenvalue oscillates between positive and negative values while decaying (or growing, if $|\lambda| > 1$).

Complex eigenvalues give modes that rotate. A complex conjugate pair $a \pm bi$ corresponds to a 2-by-2 rotation-scaling block:

The magnitude $|\lambda| = \sqrt{a^2 + b^2}$ controls whether the rotation spirals inward (decaying), outward (growing), or holds a fixed radius (neutral). The angle $\omega = \arctan(b/a)$ controls the rotation frequency. A pair of complex eigenvalues on the unit circle gives a perpetual oscillation. Inside it, a decaying spiral. Outside, an exploding one.

Churchland et al. [3] fit a linear dynamics model to motor cortex population activity during reaching. They found that the dominant eigenvalues were complex: the neural trajectory rotated in state space. The rotation frequency matched the reach duration. Their method, jPCA, was essentially finding the subspace where the dynamics matrix $A$ has complex eigenvalues with the largest imaginary parts. The rotational structure was not imposed; it was extracted from the data by looking for the eigenvalues that produce rotation.*The distinction between discrete time ($\mathbf{x}_{t+1} = A\mathbf{x}_t$) and continuous time ($d\mathbf{x}/dt = A\mathbf{x}$) matters more than it might seem. In discrete time, the stability boundary is the unit circle: $|\lambda| = 1$. In continuous time, it is the imaginary axis: $\text{Re}(\lambda) = 0$. Most neural data is sampled in discrete time (binned spike counts, calcium imaging frames), so this series uses the discrete formulation. The conversion is $A_{\text{discrete}} = \exp(A_{\text{continuous}} \cdot \Delta t)$, where $\Delta t$ is the time bin width.

So the eigendecomposition from Part 5 is not just a factorization trick. Applied to a dynamics matrix $A$, it gives you the independent modes of the system, their timescales, and their frequencies. Every linear dynamical system is, in the eigenvector basis, a collection of independent scalars and 2D rotators. The complexity comes from how many modes there are, where their eigenvalues sit, and how the initial state projects onto each one.

What we observe is not what evolves

The state vector $\mathbf{x}_t$ is latent. We never record it. What we actually measure is the activity of $N$ neurons, each of which reflects the latent state through its own linear weighting. If the latent state is $d$-dimensional and we record $N$ neurons with $N \gg d$, the observation equation is:

Here $C$ is an $N \times d$ observation matrix and $\mathbf{w}_t$ is noise. Each row of $C$ belongs to one neuron. Neuron $i$ fires at a rate proportional to $C_{i,:} \cdot \mathbf{x}_t$, plus noise. The matrix $C$ tells you how each neuron weighs the latent dimensions: how much of mode 1, how much of mode 2, and so on.

Consider what this means concretely. Suppose the latent state is a clean 2D spiral, two dimensions rotating with a slow decay, exactly the kind of trajectory Churchland et al. [3] found in motor cortex. Project that spiral through a $100 \times 2$ matrix $C$, and you get 100 correlated time series. Each neuron sees a different linear combination of the same two latent variables. Add noise, and the spiral is invisible in individual traces. Neuron 37 shows a bump then a dip. Neuron 84 shows two bumps. No single neuron reveals the structure. The structure lives in the relationships between neurons, encoded in $C$.

This is the same idea we encountered in Part 2 (a basis change hides structure) and Part 4 (projection loses dimensions), but now it plays out over time. The question is no longer "which subspace does the data live in" but "which subspace does the data evolve in." The latent state has dynamics; the observations are a high-dimensional shadow of those dynamics.

Together, Equations 2 and 4 form a linear dynamical system (LDS), also called a state-space model. The state equation (Equation 2) governs how the latent state evolves. The observation equation (Equation 4) governs what we see. Fitting an LDS to neural data means estimating $A$, $C$, and the noise statistics from the recorded spike trains or fluorescence traces.*GPFA [4] and LFADS [5] both use this state-space structure. GPFA assumes linear dynamics with Gaussian process smoothing of the latent trajectories. LFADS replaces the linear dynamics with a recurrent neural network, gaining the ability to model nonlinear evolution. But the observation equation is the same in both: $\mathbf{y} = C\mathbf{x} + \text{noise}$. The split between latent dynamics and observed projection is the shared scaffold.

The eigenvalue map

Every linear dynamical system is characterized by the eigenvalues of its dynamics matrix $A$. Plot them in the complex plane.

The unit circle is the stability boundary. An eigenvalue inside the circle corresponds to a decaying mode. Outside, a growing mode. On the circle, a mode that neither grows nor decays, sustained indefinitely.

Position on the real axis matters too. An eigenvalue sitting on the real axis produces no oscillation. The mode simply scales up or down each time step. An eigenvalue off the real axis oscillates. The angle from the positive real axis gives the frequency: larger angles mean faster rotation. The distance from the origin gives the growth or decay rate per step.

Translate this to neural populations. Slow modes (eigenvalues near 1 on the real axis) correspond to persistent activity. A prefrontal population holding a stimulus identity across a delay period lives in a mode like this. Fast modes (eigenvalues near 0) correspond to transient responses that die out within a few time steps: a brief sensory-evoked volley that fades before the next stimulus arrives. Complex eigenvalues correspond to rhythmic or rotational structure. The motor cortex reaching dynamics from Section 2, where the neural trajectory traces a spiral in state space, arise from a conjugate pair sitting just inside the unit circle. The angle sets the rotation speed; the radius sets how quickly the spiral decays.

The eigenvalue map is the complete signature of a linear dynamical system. Two systems with the same eigenvalues, even when expressed in entirely different coordinate systems with different observation matrices and different neuron counts, have the same qualitative dynamics. This is the dynamical analogue of the basis-change story from Part 2: the eigenvalues are invariant, the coordinates are not.*For neural data, you estimate $A$ from recordings and then inspect its eigenvalues to learn the population's dynamical repertoire. Macke et al. [7] used this approach to characterize the timescales and oscillation frequencies of cortical population dynamics, reading off the decay rates and rotational periods directly from the eigenvalue spectrum.

What comes next

The state-space model has two unknowns: the dynamics matrix $A$ and the observation matrix $C$. Given recorded neural data $\mathbf{y}_1, \ldots, \mathbf{y}_T$, can we recover them? In practice we record only the observations (the spiking of hundreds of neurons across time), and the latent state is never directly accessible. The question of recovering the model from data is not just a technical convenience; it is the only way to test whether a linear dynamical system is actually a good description of how a population evolves.

This is the system identification problem. The key tool turns out to be the Hankel matrix, built from time-lagged observations, where row $i$and column $j$ holds the observation at lag $i + j - 2$. Its SVD reveals the latent dimensionality and the subspace structure of the system: the singular values tell you how many dimensions the dynamics actually occupy, and the singular vectors encode the observation and state subspaces. The next post derives subspace identification[8] from this starting point, connecting the Hankel construction back to the SVD and low-rank approximation tools developed earlier in this series.

References

1. Romo, R., Brody, C. D., Hernández, A., and Lemus, L. "Neuronal correlates of parametric working memory in the prefrontal cortex," Nature, vol. 399, pp. 470-473, 1999.
2. Goldman, M. S. "Memory without feedback in a neural network," Neuron, vol. 61, no. 4, pp. 621-634, 2009.
3. Churchland, M. M., Cunningham, J. P., Kaufman, M. T., et al. "Neural population dynamics during reaching," Nature, vol. 487, pp. 51-56, 2012.
4. Yu, B. M., Cunningham, J. P., Santhanam, G., et al. "Gaussian-process factor analysis for low-dimensional single-trial analysis of neural population activity," Journal of Neurophysiology, vol. 102, no. 1, pp. 614-635, 2009.
5. Pandarinath, C., O'Shea, D. J., Collins, J., et al. "Inferring single-trial neural population dynamics using sequential auto-encoders," Nature Methods, vol. 15, no. 10, pp. 805-815, 2018.
6. Shenoy, K. V., Sahani, M., and Churchland, M. M. "Cortical control of arm movements: a dynamical systems perspective," Annual Review of Neuroscience, vol. 36, pp. 337-359, 2013.
7. Macke, J. H., Buesing, L., Cunningham, J. P., et al. "Empirical models of spiking in neural populations," Advances in Neural Information Processing Systems, vol. 24, 2011.
8. Van Overschee, P. and De Moor, B. Subspace Identification for Linear Systems. Kluwer Academic Publishers, 1996.