We analyzed the population responses of 2343 cells from the AC of three awake mice recorded using calcium imaging techniques (data from Deneux et al., 2016). The neurons were recorded while the mice passively listened to randomized presentations of different auditory stimuli. In this study we consider a total of 16 stimuli, consisting of two groups of intensity modulated UP- or DOWN-ramping sounds. In each group, there were stimuli with different frequency content (either 8 kHz pure tones or white noise [WN] sounds), different durations (1 or 2 s) and different intensity slopes (either 50–85 dB or 60–85 dB and reversed, see Table 1 and Materials and methods, Section 'The data set').

We first illustrate the responses of single cells to the presentation of different auditory stimuli, focusing on the periods following the onset and offset of the stimulus. The activity of individual neurons to different stimuli was highly heterogeneous. In response to a single stimulus, we found individual neurons that were strongly active only during the onset of the stimulus (ON responses), or only during the offset (OFF responses), while other neurons in the population responded to both stimulus onset and offset, consistent with previous analyses (Deneux et al., 2016). Importantly, across stimuli some neurons in the population exhibited ON and/or OFF responses only when specific stimuli were presented, showing stimulus-selectivity of transient responses, while others strongly responded at the onset and/or at the offset of multiple stimuli (Figure 1A).

Figure 1

Download asset Open asset

Because of the intrinsic heterogeneity of single-cell responses, we examined the structure of the transient ON and OFF responses to different stimuli using a population approach (Buonomano and Maass, 2009; Saxena and Cunningham, 2019). The temporal dynamics of the collective response of all the neurons in the population can be represented as a sequence of states in a high-dimensional state space, in which the i-th coordinate corresponds to the (baseline-subtracted) firing activity $r_i(t)$ of the i-th neuron in the population. At each time point, the population response is described by a population activity vector $\mathbf{𝐫}(t)$ which draws a neural trajectory in the state space (Figure 1B left panel).

To quantify the strength of the population transient ON and OFF responses, we computed the distance of the population activity vector from its baseline firing level (average firing rate before stimulus presentation), corresponding to the norm of the population activity vector $||\mathbf{𝐫}(t)||$ (Mazor and Laurent, 2005). This revealed that the distance from baseline computed during ON and OFF responses was larger than the distance computed for the state at the end of stimulus presentation (Figure 1B right panel). We refer to this feature of the population transient dynamics as the transient amplification of ON and OFF responses.

To examine what the transient amplification of ON and OFF responses implies in terms of stimulus decoding, we trained a simple decoder to classify pairs of stimuli that differed in their frequency content, but had the same intensity modulation and duration. We found that the classification accuracy was highest during the transient phases corresponding to ON and OFF responses, while it decreased at the end of stimulus presentation (Figure 1C). This result revealed a robust encoding of the stimulus features during ON and OFF responses, as previously found in the locust olfactory system (Mazor and Laurent, 2005; Saha et al., 2017).

To further explore the structure of the neural trajectories associated with the population OFF responses to different stimuli, we analyzed neural activity using dimensionality reduction techniques (Cunningham and Yu, 2014). We focused specifically on responses within the period starting 50 ms before stimulus offset to 300 ms after stimulus offset.

By performing principal component analysis (PCA) independently for the responses to individual stimuli, we found that the dynamics during the transient OFF responses to individual stimuli explored only about five dimensions, as 80% of the variance of the OFF responses to individual stimuli was explained on average by the first five principal components (Figure 2A; see Figure 2—figure supplement 1 for cross-validated controls; note that, given the temporal resolution, the maximal dimensionality explaining 80% of the variance of the responses to individual stimuli was 9). The projection of the low-dimensional OFF responses to each stimulus onto the first two principal components revealed circular activity patterns, where the population vector smoothly rotated between the two dominant dimensions (Figure 2B).

Figure 2 with 4 supplements see all

Download asset Open asset

A central observation revealed by the dimensionality reduction analysis was that the OFF response trajectories relative to stimuli with different frequency content spanned orthogonal low-dimensional subspaces. For instance, the response to the 8 kHz sound was poorly represented on the plane defined by the two principal components of the response to the WN sound (Figure 2B), and vice versa, showing that they evolved in distinct subspaces. To quantify the relationship between the subspaces spanned by the OFF responses to different stimuli, we proceeded as follows. We first computed the first five principal components of the (baseline-subtracted) OFF response ${\mathbf{𝐫}}^{(s)}(t)$ to each individual stimulus s. Therefore, for each stimulus s these dimensions define a five-dimensional subspace. Then, for each pair of stimuli, we computed the relative orientation of the corresponding pair of subspaces, measured by the subspace overlap (Figure 2D; see Materials and methods, Section 'Correlations between OFF response subspaces').

This approach revealed an interesting structure between the OFF response subspaces for different stimuli (Figure 2D and Figure 2—figure supplement 2). Stimuli with different frequency content evoked in general non-overlapping OFF responses reflected in low values of the subspace overlap. Two clusters of correlated OFF responses emerged, corresponding to the 8 kHz UP-ramps and WN UP-ramps of different durations and intensity. Surprisingly, DOWN-ramps evoked OFF responses that were less correlated than UP-ramps, even for sounds with the same frequency content.

The fact that most of the stimuli evoked non-overlapping OFF responses is reflected in the number of principal components that explain 80% of the variance for all OFF responses considered together, which is around 60 (Figure 2C and Figure 2—figure supplement 1). This number is in fact close to the number of dominant components of the joint response to all 16 stimuli (see Table 1) that we would expect if the responses to individual stimuli evolved on uncorrelated subspaces (given by #PC per stimulus × #stimuli ≈80). Notably this implies that while the OFF responses to individual stimuli span low-dimensional subspaces, the joint response across stimuli shows high-dimensional structure.

We finally examined to what extent the structure observed between the OFF response trajectories ${\mathbf{𝐫}}^{(s)}(t)$ to different stimuli (Figure 2D) could be predicted from the structure of the population activity states reached at the end of stimulus presentation, corresponding to the initial states ${\mathbf{𝐫}}^{(s)}(0)$. Remarkably, we found that the initial states exhibited structure across stimuli that matched well the structure of OFF response dynamics (Figure 2E left panel and Figure 2—figure supplement 4), even though the component along the initial states was substracted from the corresponding OFF response trajectories (Figure 2E right panel and Figure 2—figure supplement 4). This suggests that initial states contribute to determining the subsequent dynamics of the OFF responses.

Our analyses of auditory cortical data identified three prominent features of population dynamics: (i) OFF responses correspond to transiently amplified trajectories; (ii) responses to individual stimuli explore low-dimensional subspaces; (iii) responses to different stimuli lie in largely orthogonal subspaces. We next examined to what extent these three features could be accounted for by a single-cell mechanisms for OFF response generation.

The AC is not the first stage where OFF responses arise. Robust OFF responses are found throughout the auditory pathway, and in subcortical areas the generation of OFF responses most likely relies on mechanisms that depend on the interaction between the inhibitory and excitatory synaptic inputs to single cells (e.g. post-inhibitory rebound or facilitation, see Kopp-Scheinpflug et al., 2018). To examine the possibility that OFF responses in AC are consistent with a similar single-cell mechanism, we considered a simplified, linear model that could be directly fit to calcium recordings in the AC. Adapting previously used models (Anderson and Linden, 2016; Meyer et al., 2016), we assumed that the cells received no external input after stimulus offset, and that the response of neuron i after stimulus offset is specified by a characteristic linear filter, which describes the cell’s intrinsic response generated by intracellular or synaptic mechanisms (Figure 3A). We moreover assumed that the shape of this temporal response is set by intrinsic properties, and is therefore identical across different stimuli. In the model, each stimulus modulates the response of a single neuron linearly depending on its activity at stimulus offset, so that the (baseline-subtracted) OFF response of neuron i to stimulus s is written as:

where ${\displaystyle r_{0,i}^{(s)}}$ is the modulation factor for neuron i and stimulus s. We estimated the single-cell responses $L_i(t)$ from the data by fitting basis functions to the responses of neurons subject to prior normalization by the modulation factors ${\displaystyle r_{0,i}^{(s)}}$ using linear regression (see Materials and methods, Section 'Single-cell model for OFF response generation'). We first fitted the single-cell model to responses to a single stimulus, and then increased the number of stimuli.

Figure 3 with 3 supplements see all

Download asset Open asset

The single-cell model accounted well for the first two features observed in the data: (i) because the single-cell responses $L_i(t)$ are in general non-monotonic, the distance from baseline of the population activity vector displayed amplified dynamics (Figure 3B); (ii) at the population level, the trajectories generated by the single-cell model spanned low-dimensional subspaces of at least two dimensions (Figure 3C) and provided excellent fits to the trajectories in the data when the response to a single stimulus was considered (Figure 3C). However, we found that the single-cell model could not account for the structure of responses across multiple stimuli. Indeed, although fitting the model to OFF responses to a single stimulus led to an excellent match with auditory cortical data (coefficient of determination $R^2=0.75$), increasing the number of simultaneously fitted stimuli led to increasing deviations from the data (Figure 3B,C), and strongly decreased the goodness of fit (Figure 3D). When fitting all stimuli at once, the goodness of fit was extremely poor (coefficient of determination $R^2=0.1$), and therefore the single-cell model could not provide useful information about the third feature of the data, the structure of subspaces spanned in response to different stimuli. A simple explanation for this inability to capture structure across stimuli is that in the single-cell model the temporal shape of the response of each neuron is the same across all stimuli, while this was not the case in the data (Figure 3E).

We contrasted the single-cell model with a network model that generated OFF responses through collective interactions between neurons. Specifically, we studied a linear network of $N$ recurrently coupled linear rate units with time evolution given by: 

The quantity $r_i(t)$ represents the deviation of the activity of the unit i from its baseline firing rate, while $J_{ij}$ denotes the effective strength of the connection from neuron j to neuron i (Figure 3F). As in the single-cell model, we assumed that the network received no external input after stimulus offset. The only effect of the preceding stimulus was to set the initial pattern of activity of the network at $t=0$. The internal recurrent dynamics then fully determined the subsequent evolution of the activity. Each stimulus s was thus represented by an initial state ${\mathbf{𝐫}}_0^{(s)}$ that was mapped onto a trajectory of activity ${\mathbf{𝐫}}^{(s)}(t)$. We focused on the dynamics of the network following stimulus offset, which we represent as $t=0$.

To quantify to what extent recurrent interactions could account for auditory cortical dynamics, we fitted the network model to OFF responses using linear regression (see Materials and methods, Section 'Fitting the network model'). For each stimulus, the pattern of initial activity ${\mathbf{𝐫}}_0^{(s)}$ was taken from the data, and the connectivity matrix $\mathbf{𝐉}$ was fitted to the responses to a progressively increasing number of stimuli.

Qualitatively, we found that the fitted network model reproduced well the first two features of the data, transient amplification and low-dimensional population trajectories (Figure 3G,H). Quantitatively, we evaluated the goodness of fit by computing the coefficient of determination $R^2$. While the goodness of fit of the network model was lower than the single-cell model when fitting the response to a single stimulus ($R^2=0.52$), for the network model the goodness of fit remained consistently high as the number of simultaneously fitted stimuli was increased (Figure 3D, $R^2=0.52$ when fitted on the responses to all stimuli). Computing the goodness of fit by taking into account the number of parameters of the network and single-cell models led to the same conclusion (Figure 3—figure supplement 1). When fitted to all stimuli at once, the fitted network model captured the structure of the subspace overlaps (Figure 3I left panel) and its correlation with the structure of initial conditions (Figure 3I right panel). In contrast to the single-cell model, the network model therefore accounted well for the third feature of the data, the structure of responses across stimuli. This can be explained by the fact that, in the network model, the temporal OFF-responses of single cells can in general differ across stimuli (see Figure 6A,D), similar to the activity in the AC (Figure 3E).

Having found that the fitted network model reproduced the three main features of the data, we next analyzed this model to identify the mechanisms responsible for each feature. The identified mechanisms provided new predictions that we tested on the effective connectivity matrix $\mathbf{𝐉}$ obtained from the fit to the data.

The first feature of the data was that the dynamics were transiently amplified, i.e. the OFF responses first deviated from baseline before eventually converging to it. A preliminary requirement to reproduce this feature is that dynamics are stable, that is, eventually decay to baseline following any initial state. This requirement leads to the usual condition that the eigenvalues of the connectivity matrix $\mathbf{𝐉}$ have real parts less than unity. Provided this basic requirement is met, during the transient dynamics from the initial state at stimulus offset, the distance from baseline can either monotonically decrease, or transiently increase before eventually decreasing. To generate the transient increase, the connectivity matrix needs to belong to the class of non-normal matrices (Trefethen and Embree, 2005; Murphy and Miller, 2009; Goldman, 2009; Hennequin et al., 2012; see Materials and methods, Section 'Normal and non-normal connectivity matrices'), but this is not a sufficient condition. More specifically, the largest eigenvalue of the symmetric part defined by ${\mathbf{𝐉}}_S=(\mathbf{𝐉}+{\mathbf{𝐉}}^T)/2$ needs to be larger than unity, while the initial state ${\mathbf{𝐫}}_0$ needs to belong to a subset of amplified patterns (Bondanelli and Ostojic, 2020; see Materials and methods, Section 'Sufficient condition for amplified OFF responses').

To test the predictions derived from the condition for amplified transient responses, that is, that $\mathbf{𝐉}$ belongs to a specific subset of non-normal matrices, we examined the spectra of the full connectivity $\mathbf{𝐉}$ and of its symmetric part ${\mathbf{𝐉}}_S$ (Bondanelli and Ostojic, 2020). The eigenvalues of the fitted connectivity matrix had real parts smaller than one, indicating stable dynamics, and large imaginary parts (Figure 4A). Our theoretical criterion predicted that for OFF responses to be amplified, the spectrum of the symmetric part ${\mathbf{𝐉}}_S$ must have at least one eigenvalue larger than unity. Consistent with this predictions, we found that the symmetric part of the connectivity had indeed a large number of eigenvalues larger than one (Figure 4B) and could therefore produce amplified responses (Figure 3G).

Figure 4

Download asset Open asset

The second feature of auditory cortical data was that each stimulus generated a population response embedded in an approximately five-dimensional subspace of the full state space. We hypothesized that low-dimensional responses in the fitted network model originated from a low-rank structure in the connectivity (Mastrogiuseppe and Ostojic, 2018), implying that the connectivity matrix could be approximated in terms of $R\ll N$ modes, that is, as

where each mode was specified by two $N$-dimensional vectors ${\mathbf{𝐮}}^{(r)}$ and ${\mathbf{𝐯}}^{(r)}$, which we term the right and left connectivity patterns respectively (Mastrogiuseppe and Ostojic, 2018). This set of connectivity patterns is uniquely defined from the singular value decomposition (SVD) of the connectivity matrix $\mathbf{𝐉}$ (see Materials and methods, Section 'Low-dimensional dynamics').

To test the low-rank hypothesis, we re-fitted the network model while constraining the rank of connectivity matrix (Figure 5A, Figure 5—figure supplement 1, Figure 5—figure supplement 2; see Materials and methods, Section 'Fitting the network model'). Progressively increasing the rank $R$, we found that $R=6$ was sufficient to capture more than 80% of the variance explained by the network model when fitting the responses to individual stimuli (Figure 5A).

Figure 5 with 3 supplements see all

Download asset Open asset

The dynamics in the obtained low-rank network model are therefore fully specified by a set of patterns over the neural population: patterns corresponding to initial states ${\mathbf{𝐫}}_0$ determined by the stimuli, and patterns corresponding to connectivity vectors ${\mathbf{𝐮}}^{(r)}$ and ${\mathbf{𝐯}}^{(r)}$ for $r=1,\mathrm{\dots },R$ that determine the network dynamics. Recurrent neural networks based on such connectivity directly generate low-dimensional dynamics: for a given stimulus s the trajectory of the dynamics lies in the subspace spanned by the initial state ${\mathbf{𝐫}}_0^{(s)}$ and the set of right connectivity patterns ${\{{\mathbf{𝐮}}^{(r)}\}}_{r=1,\mathrm{\dots },R}$ (Mastrogiuseppe and Ostojic, 2018; Bondanelli and Ostojic, 2020; Schuessler et al., 2020; Beiran et al., 2020). More specifically, the transient dynamics can be written as (see Materials and methods, Section 'Low-dimensional dynamics')

where $\mathbf{𝐔}$ and $\mathbf{𝐕}$ are $N\times R$ matrices that contain respectively the $R$ right and left connectivity vectors as columns, and ${\mathbf{𝐉}}^{ov}$ is the $R\times R$ matrix of overlaps ${\mathbf{𝐉}}_{kl}^{ov}={\mathbf{𝐯}}^{(k)T}{\mathbf{𝐮}}^{(l)}$ between left and right connectivity vectors. This overlap matrix therefore fully determines the dynamics in the network (see Materials and methods, Section 'Low-dimensional dynamics').

Inspecting the overlap matrix ${\mathbf{𝐉}}^{ov}$ obtained from the fitted connectivity matrix revealed a clear structure, where left and right vectors were strongly correlated within pairs, but uncorrelated between different pairs (Figure 5B left panel). This structure effectively defined a series of $R/2$ rank-2 channels that were orthogonal to each other. Within individual rank-2 channels, strong correlations were observed only across the two modes (e.g. between ${\mathbf{𝐯}}^{(1)}$ and ${\mathbf{𝐮}}^{(2)}$, ${\mathbf{𝐯}}^{(2)}$ and ${\mathbf{𝐮}}^{(1)}$, etc.; Figure 5B), so that the connectivity matrix corresponding to each rank-2 channel can be written as

where ${\mathrm{\Delta }}_1$ and ${\mathrm{\Delta }}_2$ are two scalar values. Rank-2 matrices of this type have purely imaginary eigenvalues given by $\pm i\sqrt{{\mathrm{\Delta }}_1{\mathrm{\Delta }}_2}$, reflecting the strong imaginary components in the eigenspectrum of the full matrix (Figure 4A). These imaginary eigenvalues lead to strongly rotational dynamics in the plane defined by ${\mathbf{𝐯}}^{(1)}$ and ${\mathbf{𝐯}}^{(2)}$ (Figure 6C; see Materials and methods, Section 'Dynamics of a low-rank rotational channel'), qualitatively similar to the low-dimensional dynamics in the data (Figure 2B). Such rotational dynamics however do not necessarily imply transient amplification. In fact, a rank-2 connectivity matrix as in Equation (5) generates amplified dynamics only if two conditions are satisfied (Figure 6B; see Materials and methods, Section 'Dynamics of a low-rank rotational channel'): (i) the difference $|{\mathrm{\Delta }}_2-{\mathrm{\Delta }}_1|/2$ is greater than unity; (ii) the initial state ${\mathbf{𝐫}}_0$ overlaps strongly with the left connectivity patterns ${\mathbf{𝐯}}^{(r)}$. A direct consequence of these constraints is that when dynamics are transiently amplified, the population activity vector at the peak of the transient dynamics lies along a direction that is largely orthogonal to the initial state at stimulus offset (Figure 6C,E, see Materials and methods, Section 'Correlation between initial and peak state').

Figure 6

Download asset Open asset

The theoretical analyses of the model therefore provide a new set of predictions that we directly tested on the connectivity fitted to the data. We first computed the differences $|{\mathrm{\Delta }}_2-{\mathrm{\Delta }}_1|/2$ for each channel using the SVD of the fitted matrix ${\mathbf{𝐉}}^{(s)}$, and found that they were sufficiently large to amplify the dynamics within each rank-$R$ channel (Figure 5C left panel). We next examined for each stimulus the component of the initial state ${\mathbf{𝐫}}_0$ on the ${\mathbf{𝐯}}^{(r)}$’s and found that it was significantly larger than the component on the ${\mathbf{𝐯}}^{(r)}$’s of a random vector, across all stimuli (Figure 5C right panel). Finally, computing the correlation between the peak state and the initial state at the end of stimulus presentation, we found that this correlation took values lower than values predicted by chance for almost all stimuli, consistent with the prediction of the network model (Figure 5D, Figure 5—figure supplement 3). The single-cell model instead predicts stronger correlations between initial and peak states, in clear conflict with the data (Figure 5E).

The third feature observed in the data was that population responses evoked by different stimuli were orthogonal for many of the considered stimuli. Orthogonal trajectories for different stimuli can be reproduced in the model under the conditions that (i) initial patterns of activity corresponding to different stimuli are orthogonal, (ii) the connectivity matrix is partitioned into a set of orthogonal low-rank terms, each individually leading to transiently amplified dynamics, and (iii) each stimulus activates one of the low-rank terms. Altogether, for $P$ stimuli leading to orthogonal responses, the connectivity matrix can be partitioned in $P$ different groups of modes:

with

where the vectors ${\mathbf{𝐯}}^{(r;s)}$ have unit norm. The set of modes indexed by s correspond to the s-th stimulus, so that the pattern of initial activity ${\mathbf{𝐫}}_0^{(s)}$ evoked by stimulus s overlaps only with those modes. Moreover, modes corresponding to different groups are orthogonal, and generate dynamics that span orthogonal subspaces (Figure 2D). Each term in Equation (6) can therefore be interpreted as an independent transient coding channel that can be determined from the OFF-response to stimulus s alone.

We therefore examined whether the fitted connectivity $\mathbf{𝐉}$ consisted of independent low-rank coding channels (Equation (6)), a constraint that would generate orthogonal responses to different stimuli as observed in the data. If $\mathbf{𝐉}$ consisted of fully independent channels as postulated in Equation (6), then it could be equivalently estimated by fitting the recurrent model to the OFF responses to each stimulus s independently, leading to one matrix ${\mathbf{𝐉}}^{(s)}$ for each stimulus. The hypothesis of fully independent channels then implies that (i) the connectivity vectors for the different connectivity matrices are mutually orthogonal, (ii) the full connectivity matrix can be reconstructed by summing the matrices ${\mathbf{𝐉}}^{(s)}$ estimated for individual stimuli. To test these two predictions, we estimated the connectivity matrices ${\mathbf{𝐉}}^{(s)}$ (the rank parameter for each stimulus was set to $R=5$). We then computed the overlaps between the connectivity vectors corresponding to different pairs of stimuli (see Materials and methods, Section 'Analysis of the transient channels') and compared them to the overlaps between the subspaces spanned in response to the same pairs of stimuli (see Figure 2D). We found a close match between the two quantities (Figure 7A,B): pairs of stimuli with strong subspace overlaps corresponded to high connectivity overlaps, and pairs of stimuli with low subspace overlaps corresponded to low connectivity overlaps. This indicated that some of the stimuli, but not all, corresponded to orthogonal coding channels. We further reconstructed the OFF responses using the matrix ${\mathbf{𝐉}}_{\mathrm{Sum}}={\sum }_s{\mathbf{𝐉}}^{(s)}$. We then compared the goodness of fit computed using the matrices ${\mathbf{𝐉}}_{\mathrm{Sum}}$, ${\mathbf{𝐉}}_{\mathrm{Full}}$ and multiple controls where the elements of the matrix ${\mathbf{𝐉}}_{\mathrm{Sum}}$ were randomly shuffled. While the fraction of variance explained when using the matrix ${\mathbf{𝐉}}_{\mathrm{Sum}}$ was necessarily lower than the one computed using ${\mathbf{𝐉}}_{\mathrm{Full}}$, the model with connectivity ${\mathbf{𝐉}}_{\mathrm{Sum}}$ yielded values of the goodness of fit significantly higher than the ones obtained from the shuffles (Figure 7C). These results together indicated that the full matrix ${\mathbf{𝐉}}_{\mathrm{Full}}$ can indeed be approximated by the sum of the low-rank channels represented by the ${\mathbf{𝐉}}^{(s)}$.

Figure 7

Download asset Open asset

So far we examined the activity obtained by averaging for each neuron the response to the same stimulus over multiple trials. We next turned to trial-to-trial variability of simultaneously recorded neurons and compared the structure of the variability in the data with the network and single-cell models. Specifically, we examined the variance of the activity along a direction ${\mathbf{𝐳}}_0$ of state-space (Figure 8A; Hennequin et al., 2012), defined as:

where $⟨\cdot ⟩$ denotes the averaging across trials. We compared the variance along the direction ${\mathbf{𝐫}}_{\mathrm{ampl}}$ corresponding to the maximum amplification (distance from baseline) of the trial-averaged dynamics with variance along a random direction ${\mathbf{𝐫}}_{\mathrm{rand}}$.

Figure 8

Download asset Open asset

Inspecting the variability during the OFF-responses in the AC data revealed two prominent features: (i) fluctuations across trials are amplified along the same direction ${\mathbf{𝐫}}_{\mathrm{ampl}}$ of state-space as trial-averaged dynamics, but not along random directions (Figure 8B,C,D); (ii) cancelling cross-correlations by shuffling trial labels independently for each cell strongly reduces the amplification of the variance (Figure 8C,D). This second feature demonstrates that the amplification of variance is due to noise-correlations across neurons that have been previously reported in the AC (Rothschild et al., 2010). Indeed the variance along a direction ${\mathbf{𝐳}}_0$ at time $t$ can be expressed as

where $\mathbf{𝐂}(t)$ represents the covariance matrix of the population activity at time $t$. Shuffling trial labels independently for each cell effectively cancels the off-diagonal elements of the covariance matrix and leaves only the diagonal elements that correspond to single-neuron variances. The comparison between simultaneous and shuffled data (Figure 8C,D) demonstrates that an increase in single-neuron variance is not sufficient to explain the total increase in variance along the amplified direction ${\mathbf{𝐫}}_{\mathrm{ampl}}$.

Simulations of the fitted models and a mathematical analysis show that the network model reproduces both features of trial-to-trial variability observed in the data (Figure 8E,F), based on a minimal assumption that variability originates from independent noise in the initial condition at stimulus offset (see Materials and methods, Section 'Structure of single-trial responses in the network model'). Surprisingly, the single-cell model also reproduces the first feature, the amplification of fluctuations along the same direction of state-spate as trial-averaged dynamics, although in that model the different neurons are not coupled and noise-correlations are absent (Figure 8G,H). Instead, the single-cell model fails to capture the second feature, as it produces amplified variability in both simultaneous and shuffled activity. This demonstrates that the amplification of variability in the single-cell model is due to an amplification of single-cell variances, that is, the diagonal elements of covariance matrix (see Materials and methods, Section 'Structure of single-trial responses in the single-cell model'). Somewhat counter-intuitively, the disagreement between the single-cell model and the data is not directly due to the lack of noise-correlations in that model, but due to the fact that the single-cell model does not predict accurately the variance of single-neuron activity during the OFF-responses, and therefore fails to capture variance in shuffled activity.

In summary, the network model accounts better than the single-cell model for the structure of single-trial fluctuations present in the AC data.

To perform PCA on the population responses to $C$ stimuli $s_{i_1},\mathrm{\dots },s_{i_C}$ we considered the matrix $\mathbf{𝐗}\in {ℝ}^{N\times TC}$, where $N$ is the number of neurons and $T$ is the number of time steps considered. The matrix $\mathbf{𝐗}$ contained the population OFF responses to the stimuli $s_{i_1},\mathrm{\dots },s_{i_C}$, centered around the mean over times and stimuli. If we denote by ${\lambda }_i$ the i-th eigenvalue of the covariance matrix ${\displaystyle \mathbf{X}{\mathbf{X}}^T/(TC)}$, the percentage of variance explained by the i-th principal component is given by:

while the cumulative percentage of variance explained by the first $M$ principal components (Figure 2A,C) is given by:

In Figure 2A, for each stimulus $s$ we consider the matrix ${\mathbf{𝐗}}^{(s)}\in {ℝ}^{N\times T}$ containing the population OFF responses to stimulus $s$. We computed the cumulative variance explained as a function of the number of principal components $M$ and then averaged over all stimuli.

We used cross-validation to estimate the variance of OFF responses to individual stimuli and across stimuli attributable to the stimulus-related component, discarding the component due to trial-to-trial variability (Figure 2—figure supplement 1). Specifically, we applied to the calcium activity data the method of cross-validated PCA (cvPCA) developed in Stringer et al., 2019. This method provides an unbiased estimate of the stimulus-related variance by computing the covariance between a training and a test set of responses (e.g. two different trials or two different sets of trials) to an identical collection of stimuli. Let ${\mathbf{𝐗}}^{(train)}$ and ${\mathbf{𝐗}}^{(test)}$ be the $N\times TC$ matrices containing the mean-centered responses to the same $C$ stimuli. We consider the training and test responses to be two distinct trials. We first perform ordinary PCA on the training responses ${\mathbf{𝐗}}^{(train)}$ and find the principal component ${\mathbf{𝐮}}_i^{(train)}$ (${\displaystyle i=1,...,N}$). We then evaluate the cross-validated PCA spectrum $\{{\lambda }_i\}$ as:

We repeat the procedure for all pairs of trials $(i,j)$ with $i\ne j$ and average the result over pairs. The cross-validated cumulative variance is finally computed as in Equation (13).

To quantify the degree of correlation between pairs of OFF responses corresponding to two different stimuli, termed subspace overlap, we computed the cosine of the principal angle between the corresponding low-dimensional subspaces (Figure 2D, Figure 3I). In general, the principal angle ${\theta }_P$ between two subspaces $U$ and $V$ corresponds to the largest angle between any two pairs of vectors in $U$ and $V$ respectively, and it is defined by Bjorck and Golub, 1973; Knyazev and Argentati, 2002:

To compute the correlations between the OFF responses to stimuli s₁ and s₂ we first identified the first $K=5$ principal components of the response to stimulus s₁ and organized them in a $N\times K$ matrix $\mathbf{𝐐}(s_1)$. We repeated this for stimulus s₂, which yielded a matrix $\mathbf{𝐐}(s_2)$. Therefore the columns of $\mathbf{𝐐}(s_1)$ and $\mathbf{𝐐}(s_2)$ define the two subspaces on which the responses to stimuli s₁ and s₂ live. The cosine of the principal angle between these two subspaces is given by Bjorck and Golub, 1973; Knyazev and Argentati, 2002:

where ${\sigma }_{\max }(\mathbf{𝐀})$ denotes the largest singular value of a matrix $\mathbf{𝐀}$. We note that this procedure directly relates to canonical correlation analysis (CCA; see Hotelling, 1936; Uurtio et al., 2018). In particular the first principal angle corresponds to the first canonical weight between the subspaces spanned by the columns of $\mathbf{𝐐}(s_1)$ and $\mathbf{𝐐}(s_2)$ (Golub and Zha, 1992; Bjorck and Golub, 1973).

In this section we describe the controls that we used to evaluate the statistical significance of the measures of orthogonality between subspaces spanned by neural activity during OFF responses to different stimuli (Figure 2D, Figure 3I), and between initial and peak activity vectors for a single stimulus (Figure 5D). We assessed the significance of the orthogonality between OFF response subspaces across stimuli using two separate controls, which aim at testing for two distinct null hypotheses (Figure 2D, Figure 3I, Figure 2—figure supplement 2). The first hypothesis is that small subspace overlaps (i.e. low correlations) between OFF responses to different stimuli may be due to the high number of dimensions of the state-space in which they are embedded. To test for this hypothesis we compared the subspace overlap computed on the trial-averaged activity with the subspace overlaps computed on the trial-averaged activity where the stimulus labels had been shuffled across trials for each pair of stimuli. We shuffled stimulus labels multiple times, resulting in one value of the subspace overlap for each shuffle. For each pair of stimuli, significance levels were then computed as the fraction of shuffles for which the subspace overlap was lower than the subspace overlap for the real data (lower tail test; Figure 2—figure supplement 2A).

Alternatively, small subspace overlaps could be an artifact of the trial-to-trial variability present in the calcium activity data. In fact, for a single stimulus, maximum values of the correlation between two different trials were of the order of 0.2 (Deneux et al., 2016). To test for the possibility that small subspace overlaps may be due to trial-to-trial variability, for each pair of stimuli we computed the values of the subspace overlaps by computing the trial-averaged activity on only half of the trials (10 trials), subsampling the set of 10 trials multiple times for each stimulus. This yielded a set of values $\cos {\theta }_{\mathrm{real}}(s_1,s_2,n)$, where s₁ and s₂ are the specific stimuli considered and $n=1,\mathrm{\dots },N_{\mathrm{shuffle}}$. We then computed the subspace overlaps between the trial-averaged responses to the same stimulus, but averaged over two different sets of 10 trials each, over multiple permutations of the trials, resulting in a set of values $\cos {\theta }_{\mathrm{shuffle}}(s,n)$, where $s\in \{s_1,s_2\}$ is the specific stimulus considered and $n=1,\mathrm{\dots },N_{\mathrm{shuffle}}$. For each pair of stimuli s₁ and s₂, significance levels were computed using two-tailed t-test and taking the maximum between the p-values given by $p(\cos {\theta }_{\mathrm{real}}(s_1,s_2,:),\cos {\theta }_{\mathrm{shuffle}}(s_1,:))$ and $p(\cos {\theta }_{\mathrm{real}}(s_1,s_2,:),\cos {\theta }_{\mathrm{shuffle}}(s_2,:))$ for those stimulus pairs for which ${⟨\cos {\theta }_{\mathrm{real}}(s_1,s_2,:)⟩}_n<{⟨\cos {\theta }_{\mathrm{shuffled}}(s_1,:)⟩}_n$ and ${⟨\cos {\theta }_{\mathrm{real}}(s_1,s_2,:)⟩}_n<{⟨\cos {\theta }_{\mathrm{shuffled}}(s_2,:)⟩}_n$, where the bar symbol indicated the mean over shuffles (Figure 2—figure supplement 2B).

The same null hypotheses have been used to test the significance of the orthogonality between initial state and peak state for individual stimuli (Figure 5D, Figure 5—figure supplement 3). A procedure analogous to the one outlined above was employed. Here, instead of the subspace overlaps, the quantity of interest is the correlation (scalar product) between the initial and peak state. For the first control shuffled data are obtained by shuffling the labels ‘initial state’ and ‘peak state’ across trials. Significance levels were evaluated as outlined above for the first control (Figure 5—figure supplement 3A). To test for the second hypothesis, we computed correlations between activity vectors defined at the same time point, but averaged over two different sets of 10 trials each. For each trial permutations we computed these correlations for all time points belonging to the OFF response (35 time points) and average over time points. Significance level was then assessed as outlined above for the second control (Figure 5—figure supplement 3B).

In the next section we describe the procedure used to fit the single-cell model to the auditory cortical OFF responses. We consider the single-cell model given by Equation (1), where ${\displaystyle r_{0,i}^{(s)}{L}_i(0)}$ represents the initial state of the response of unit i to stimulus s. Without loss of generality we assume that the dynamic range of the temporal filters, defined as $|{\max }_t{L}_i(t)-{\min }_t{L}_i(t)|$, is equal to unity, so that ${\displaystyle r_{0,i}^{(s)}}$ represents the firing rate range of neuron i for stimulus s. If that was not true, i.e. if the single-neuron responses were given by ${\displaystyle r_i^{(s)}(t)=r_{0,i}^{(s)}{K}_i(t)}$ with ${\alpha }_i=|{\max }_t{K}_i(t)-{\min }_t{K}_i(t)|\ne 1$ we could always write the responses as ${\displaystyle r_i^{(s)}(t)={\tilde{r}}_{0,i}^{(s)}{\tilde{L}}_i(t)}$, where ${\displaystyle {\tilde{r}}_{0,i}^{(s)}=r_{0,i}^{(s)}{\alpha }_i}$ and ${\tilde{L}}_i(t)=K_i(t)/{\alpha }_i$, formally equivalent to Equation (1).

To fit the single-cell OFF responses ${\displaystyle r_i^{(s)}(t)}$, we expressed the single-cell responses on a set of basis functions (Pillow et al., 2008).

where the shape and the number of basis function $N_{\mathrm{basis}}$ are predetermined. We choose the functions $f_i(t)$ to be Gaussian functions centered around a value ${\overline{t}}_i$ and with a given width ${\displaystyle w_i}$, that is, $f_i(t)=\exp (-{(t-{\overline{t}}_i)}^2/2w_i^2)$. The problem then consists in finding the coefficients $a_{ij}^{(s)}$ that best approximate Equation (17). By dividing the left- and right-hand side of Equation (17) by the firing rate range ${\displaystyle r_{0,i}^{(s)}}$ we obtain:

In general the coefficients $b_{ij}^{(s)}$ could be fitted independently for each stimulus. However, the single-cell model assumes that the coefficients $b_{ij}$ do not change across stimuli (see Equation (1)). Therefore, to find the stimulus-independent coefficients $b_{ij}$ that best approximate Equation (18) across a given set of stimuli, we minimize the mean squared error given by:

The minimization of the mean squared error can be solved using linear regression techniques. Suppose we want to fit the population responses to $C$ different stimuli simultaneously. Let ${\mathbf{𝐑}}^{(C)}$ be the matrix of size $N\times TC$ obtained by concatenating the $N\times T$ matrices ${\displaystyle ({\mathbf{r}}_{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}^{(s)}{)}_{it}=r_i^{(s)}(t)/r_{0,i}^{(s)}}$ ($i=1,\mathrm{\dots },N$, $t=1,\mathrm{\dots },T$, $s=1,\mathrm{\dots },C$) corresponding to the normalized responses to the $C$ stimuli. Let ${\mathbf{𝐅}}^{(C)}$ be the $N_{\mathrm{basis}}\times TC$ matrix obtained by concatenating $C$ times the $N_{\mathrm{basis}}\times T$ matrix ${(\mathbf{𝐟})}_{it}=f_i(t)$. Let $\mathbf{𝐁}$ be the $N\times N_{\mathrm{basis}}$ matrix given by ${\mathbf{𝐁}}_{ij}=b_{ij}$. Then Equation (18) can be written as:

which can be solved using linear regression.

In order to avoid normalizing by very small values, we fit only the most responding neurons in the population, as quantified by their firing rate range. We set $w_i=w=35ms$ for all $i=1,\mathrm{\dots },N_{\mathrm{basis}}$, and we set the number of basis functions to $N_{\mathrm{basis}}=10$, corresponding to the minimal number of basis functions sufficient to reach the highest cross-validated goodness of fit.

We consider a recurrent network model of $N$ randomly coupled linear rate units. Each unit $i$ is described by the time-dependent variable $r_i(t)$, which represents the difference between the firing rate of neuron $i$ at time $t$ and its baseline firing level. The equation governing the temporal dynamics of the network reads:

where $\tau$ represents the time constant (fixed to unity), and $J_{ij}$ is the effective synaptic strength from neuron $j$ to neuron $i$. The system has only one fixed point corresponding to $r_i=0$ for all $i$. To have stable dynamics, we require that the real part of the eigenvalues of the connectivity matrix $\mathbf{𝐉}$ is smaller than unity, that is, ${\left(\Re 𝔢\lambda (\mathbf{𝐉})\right)}_{\max }<1$. We represent each stimulus as the state of the system reached at the end of stimulus presentation, which we denote by the vector ${\mathbf{𝐫}}_0$. We moreover assume that during the OFF response the network receives no external input. The OFF response to the stimulus associated with the state ${\mathbf{𝐫}}_0$ therefore corresponds to the network dynamics starting from the initial condition ${\mathbf{𝐫}}_0$.

In this section we describe the procedure we used for fitting a linear recurrent network model (Equation (21)) to the auditory cortical data using ridge and reduced-rank ridge regression. To fit the network model to the responses to $C$ stimuli, we first concatenate the $C$ matrices ${\mathbf{𝐗}}^{(s)}\in {ℝ}^{T\times D}$ (where ${\mathbf{𝐗}}_{t,i}^{(s)}=r_i^{(s)}(t)$, with $s=1,\mathrm{\dots },C$ and $i=1,\mathrm{\dots },D$) containing the neural responses to stimulus s across ${\displaystyle T}$ timesteps and ${\displaystyle D}$ neurons (or PC dimensions) along the first dimension, obtaining a matrix $\mathbf{𝐗}\in {ℝ}^{CT\times D}$ contains the neural activity across $C$ stimuli, $T$ time steps, and $D$ neurons or PC dimensions. We then fit the network model $\dot{\mathbf{𝐗}}=\mathbf{𝐗}(\mathbf{𝐉}-\mathbf{𝐈})$ by first computing from the data the velocity of the trajectory as

Both ridge and reduced-rank ridge regression aim at minimizing the mean squared error ${||\dot{\mathbf{𝐗}}-\mathbf{𝐗}(\mathbf{𝐉}-\mathbf{𝐈})||}^2$ subject to different constraints, which are described below. Since ridge regression involves computationally expensive matrix inversion, we first reduced the dimensionality of the original data set by using PCA. Unless otherwise stated, we kept a number of PC components equal to K = 100 for the fit of individual stimuli (Figure 3G,H, Figure 5A–C, Figure 8, Figure 7) and for the fit of all stimuli together (Figure 3I, Figure 4, Figure 7C, Figure 3—figure supplement 2, Figure 3—figure supplement 3, Figure 5—figure supplement 1; K = 100 principal components sufficed to explain more than 90% of the variance of the responses to all stimuli together). As a result, the data matrix $\mathbf{𝐗}\in {ℝ}^{TC\times K}$ contains the activity across $T$ time bins and across all $K$ dimensions for a number $C$ of stimuli.

Ridge regression aims at determining the connectivity matrix $\mathbf{𝐉}$ that minimizes the mean squared error ${||\dot{\mathbf{𝐗}}-\mathbf{𝐗}(\mathbf{𝐉}-\mathbf{𝐈})||}^2$ with a penalty for large entries of the matrix $\mathbf{𝐉}$, so that the cost function to minimize is given by ${||\dot{\mathbf{𝐗}}-\mathbf{𝐗}(\mathbf{𝐉}-\mathbf{𝐈})||}^2+\lambda {||(\mathbf{𝐉}-\mathbf{𝐈})||}^2$, where $\lambda$ is a hyperparameter of the model. We can write the ridge regression optimization problem as:

where we defined ${\dot{\mathbf{𝐗}}}_{\lambda }=(\dot{\mathbf{𝐗}},\mathrm{})$ and ${\mathbf{𝐗}}_{\lambda }=(\mathbf{𝐗},\sqrt{\lambda }\mathbf{𝐈})$. The solution to Equation (23) is given by:

This procedure was used for fits shown in Figure 3, Figure 3—figure supplement 1, Figure 3—figure supplement 3, and Figure 4.

Reduced-rank regression aims at minimizing the mean squared error ${||\dot{\mathbf{𝐗}}-\mathbf{𝐗}(\mathbf{𝐉}-\mathbf{𝐈})||}^2$ under a rank constraint on the matrix $\mathbf{𝐉}$, that is, $\mathrm{rank}\mathbf{𝐉}\le R$, where $R$ is a hyperparameter of the model (Izenman, 1975; Davies and Tso, 1982). Here we combined reduced-rank and ridge regression in the same framework. The reduced-rank ridge regression optimization problem with hyperparameters $R$ and $\lambda$ can be therefore written as (Mukherjee et al., 2015):

To solve Equation (25) we consider the solution to the ridge regression problem given by Equation (24). If the matrix ${\mathbf{𝐗}}_{\lambda }{\mathbf{𝐉}}_{\lambda }^{*}$ has SVD given by ${\mathbf{𝐗}}_{\lambda }{\mathbf{𝐉}}_{\lambda }^{*}=\mathbf{𝐔}\mathrm{\Sigma }{\mathbf{𝐕}}^T$, then it can be shown that the solution to the reduced-rank ridge regression problem given by Equation (25) can be written as:

We note that each term in the sum of Equation (26) has unit rank, so that the resulting matrix ${\mathbf{𝐉}}_{R,\lambda }^{*}$ has rank equal to or less than $R$.

This procedure was used for fits shown in Figure 5, Figure 7, Figure 8, Figure 3—figure supplement 2, Figure 5—figure supplement 1, and Figure 5—figure supplement 2.

To select the value of the hyperparameter $\lambda$, we fitted the network model to the data and computed the goodness of fit as the coefficient of determination $R^2(\lambda )$ using cross-validation. We then selected the value of $\lambda$ which maximized the goodness of fit.

To select the values of the hyperparameters $\lambda$ and $R$, we fitted the network model to the data with hyperparameters $(\lambda ,R)$ and computed the goodness of fit as the coefficient of determination $R^2(\lambda ,R)$ using cross-validation. We repeated the process for a range of values of $(\lambda ,R)$. We observed that, independent of the value of $\lambda$, the goodness of fit as a function of the rank $R$ saturates at a particular value of the rank $R^{*}$, but does not exhibit a clear maximum. We took the value $R^{*}$ as the minimal rank hyperparameter, while we defined the best ridge parameter as ${\lambda }^{*}={\mathrm{argmax}}_{\lambda }{R}^2(\lambda ,R^{*})$ (Figure 5—figure supplement 1).

For both ridge and reduced-rank ridge regression, we used $K$-fold cross-validation, with $K=10$. When fitting multiple stimuli at once, for each stimulus we partitioned the temporal activity into $K$ chunks, resulting in a total of $KC$ chunks. At the $i$-th iteration of the cross-validation procedure, we leave out the $i$-th partition for each stimulus to construct the training set (consisting of $(K-1)C$ chunks) and test the trained model on the remaining $C$ folds.

In Figure 4A,B, the temporal window considered for the fit was extended from 350 ms to 600 ms to include the decay to baseline of the OFF responses, thus obtaining stable eigenvalues. The extension of the temporal window was possible only for the 1 s long stimuli ($n=8$), since the length of the temporal window following stimulus offset for the 2 s stimuli was limited to approximately 380 ms by the length of the neural recordings.

To evaluate whether the fitted network model captured nontrivial collective dynamics of auditory cortical OFF responses, we followed the approach introduced in Elsayed and Cunningham, 2017. We first computed the goodness of fit (as quantified by the coefficient of determination $R^2$), then repeated the model fitting on control data sets which kept no additional structure than the one defined by correlations across time, neurons and stimuli of the original neural responses. We found that the goodness of fit obtained from fitting the model to the original data set was significantly higher than the one obtained from fitting the control data sets (Figure 3—figure supplement 2, Figure 3—figure supplement 3), confirming that the recurrent network model captured nontrivial collective dynamics in the data beyond the correlations across time, neurons and stimuli.

We generated the control data sets using a recent method based on a maximum entropy model (Savin and Tkačik, 2017) and described in Elsayed and Cunningham, 2017. This method, termed Tensor Maximum Entropy, allowed us to build surrogate data sets that are maximally random (entropy maximization), but constrained in a way that their marginal means and covariances are equal to the means and covariances of the original data set.

Let the temporal activity along all $K$ dimensions for all $C$ stimuli be organized in a tensor $\mathbf{𝐙}\in {ℝ}^{T\times K\times C}$. The mean tensor $\mathbf{𝐌}$ is defined as the tensor that makes all the marginal means of $\mathbf{𝐙}$ vanish. Specifically, if $\overline{\mathbf{𝐙}}=\mathbf{𝐙}-\mathbf{𝐌}$, the tensor $\mathbf{𝐌}$ is such that:

The marginal covariances of the tensor $\overline{\mathbf{𝐙}}$ across times, neural dimensions, and stimuli are therefore defined as:

The method generates the desired number of surrogate data sets ${\mathbf{𝐒}}^{(i)}\in {ℝ}^{T\times K\times C}$. Each of these surrogates is randomly drawn from a probability distribution that assumes a priori no structure apart from that expected from the marginal means and covariances of the original data. Let µ, ${\mathrm{\Lambda }}^T$, ${\mathrm{\Lambda }}^K$, and ${\mathrm{\Lambda }}^C$ be the marginal means and covariances of surrogate $\mathbf{𝐒}$. The method computes the probability $P(\mathbf{𝐒})$ over the surrogates that maximizes the entropy function

subject to the constraints

where ${𝔼}_P[\cdot ]$ denotes the expectation over the probability density $P$. We used three types of surrogate data sets, denoted as T, TK, and TKC. All the three types of surrogates obey to the first constraint in Equation (30) on the marginal means. In addition, surrogates of type T obey the constraint on the time covariances, surrogates of type TK on time and dimension covariance, while surrogates TKC obey all the three covariance constraints.

In this section, we provide details on the analysis of the structure of the transient channels obtained from fitting the network model to OFF responses to individual stimuli using reduce-rank regression (Figure 7, see Materials and methods, Section 'Fitting the network model'). Let ${\mathbf{𝐉}}_{\mathrm{Full}}$ be the connectivity matrix resulting from fitting the network model to the responses to all stimuli at once, and ${\mathbf{𝐉}}^{(s)}$ the connectivity obtained from fitting the model to the response to stimulus s. Before fitting both matrices ${\mathbf{𝐉}}_{\mathrm{Full}}$ and ${\mathbf{𝐉}}^{(s)}$, we reduced the dimensionality of the responses to the first 100 principal components. We set the value of the rank parameter to $R=70$ and ${\displaystyle R=5}$ respectively for ${\mathbf{𝐉}}_{\mathrm{Full}}$ and ${\mathbf{𝐉}}^{(s)}$ (Figure 5). Using the SVD, we can write the matrices ${\mathbf{𝐉}}^{(s)}$ as:

where ${\widehat{\mathbf{𝐔}}}^{(s)}$ and ${\mathbf{𝐕}}^{(s)}$ contain respectively the $R_s$ right and left connectivity vectors as columns (see Equation (6)), while ${\mathrm{\Sigma }}^{(s)}$ contains the norms of ${\mathbf{𝐔}}^{(s)}$ on the diagonal (see Equation (6)). The transient dynamics elicited by stimulus s has a strong component along the state-space dimensions specified by the right connectivity vectors ${\mathbf{𝐮}}^{(r;s)}$ (columns of ${\mathbf{𝐔}}^{(s)}$; see also Equation (43)). We therefore define the overlap between the transient channels (connectivity overlaps) corresponding to pairs of stimuli s₁ and s₂ (Figure 7A) as the principal angle between the subspaces defined by ${\widehat{\mathbf{𝐔}}}^{(s_1)}$ and ${\widehat{\mathbf{𝐔}}}^{(s_2)}$ (see Equations (15 and 16)). We show that the structure of the connectivity overlaps (Figure 7A) matched well with the structure of the subspace overlaps (Figure 2D) across pairs of stimuli in Figure 7B.

To test whether the connectivity fitted on all stimuli at once ${\mathbf{𝐉}}_{\mathrm{Full}}$ consisted of the sum of low-rank transient coding channels, we defined the matrix ${\mathbf{𝐉}}_{\mathrm{Sum}}$ as the sum of the individual transient channels for all stimuli (see Equation 6):

We then compared the cross-validated goodness of fit of the population OFF responses using the full matrix ${\mathbf{𝐉}}_{\mathrm{Full}}$ and the matrix ${\mathbf{𝐉}}_{\mathrm{Sum}}$ (Figure 7C).

In Figure 7, when fitting individual stimuli, the ridge and rank hyperparameters have been set respectively to $\lambda =5$, ${\displaystyle R=5}$, which correspond to the values that maximize the agreement (coefficient of determination) between the connectivity overlaps (Figure 7A) and the subspace overlaps (Figure 2D). For these values the coefficient of determination between the transient channel overlaps and the subspace overlaps is >0.7.

In this section we provide a detailed mathematical analysis of the network model given by Equation (21)

Specifically, we derive the conditions on the network components (i.e. connectivity spectra, connectivity vectors, and initial conditions) to produce low-dimensional, transiently amplified OFF responses in networks with low-rank connectivity. We then focus on the specific case of rotational channels in the connectivity.

In this section, we summarize the relationship between amplified OFF responses and non-normal connectivity matrices in linear recurrent networks previously reported in Trefethen and Embree, 2005; Hennequin et al., 2012; Bondanelli and Ostojic, 2020. To characterize the amplification of OFF responses we focus on the temporal dynamics of the distance from baseline, defined as the norm of the population activity vector $||\mathbf{𝐫}(t)||$ (Hennequin et al., 2014). The network generates an amplified OFF response to the stimulus associated with the initial condition ${\mathbf{𝐫}}_0$ when the value of $||\mathbf{𝐫}(t)||$ transiently increases before decaying to its asymptotic value $||\mathbf{𝐫}(t\to \mathrm{\infty })||=0$. Note that having a transiently increasing value of the distance from baseline implies that the OFF response $r_i(t)$ of at least one unit displays non-monotonic temporal dynamics. Importantly, the transient behavior of $||\mathbf{𝐫}(t)||$ depends on the stimulus through ${\mathbf{𝐫}}_0$, and on the properties of the connectivity matrix $\mathbf{𝐉}$, in particular on the relationship between its eigenvectors (Trefethen and Embree, 2005).

Connectivity matrices for which the eigenvectors are orthogonal to each other are called normal matrices and they are formally defined as matrices $\mathbf{𝐉}$ that satisfy ${\mathbf{𝐉𝐉}}^T={\mathbf{𝐉}}^T\mathbf{𝐉}$. In networks with normal connectivity, any stimulus ${\mathbf{𝐫}}_0$ evokes an OFF response for which the distance from baseline decays monotonically to zero. Such networks thus cannot produce amplified OFF responses, as defined by a transiently increasing $||\mathbf{𝐫}(t)||$. Note that any symmetric matrix is normal.

On the other hand, connectivity matrices for which some eigenvectors are not mutually orthogonal are called non-normal (Trefethen and Embree, 2005), and they consist of all connectivity $\mathbf{𝐉}$ for which ${\mathbf{𝐉𝐉}}^T\ne {\mathbf{𝐉}}^T\mathbf{𝐉}$. It is well known that non-normal networks can lead to transiently increasing values of $||\mathbf{𝐫}(t)||$, therefore producing amplified OFF responses. However, the non-normality of the network connectivity $\mathbf{𝐉}$ constitutes only a necessary, but not a sufficient condition for the generation of amplified responses.

In this section we identify the necessary and sufficient condition for the generation of amplified OFF responses in linear recurrent networks. To find the necessary and sufficient condition for amplified responses, we start by writing the differential equation for the dynamics of the distance from baseline as (Neubert and Caswell, 1997; Bondanelli and Ostojic, 2020):

where ${\mathbf{𝐉}}_S$ denotes the symmetric part of the connectivity $\mathbf{𝐉}$. A linear recurrent network exhibits amplified responses when the rate of change of the distance from baseline, $\mathrm{d}||\mathbf{𝐫}||/\mathrm{d}t$, takes positive values at time $t=0$. The right-hand side of Equation (33) takes its largest value when the initial condition ${\mathbf{𝐫}}_0$ is aligned with the eigenvector of ${\mathbf{𝐉}}_S$ associated with the largest eigenvalue ${\lambda }_{\max }({\mathbf{𝐉}}_S)$. In this case, the rate of change of the distance from baseline at time $t=0$ takes the value ${\lambda }_{\max }({\mathbf{𝐉}}_S)-1$. From Equation (33) it can be shown that the necessary and sufficient condition for the generation of amplified responses in a recurrent networks with connectivity $\mathbf{𝐉}$ is given by

This criterion defines two classes of networks based on the properties of the connectivity matrix: networks in which amplified responses cannot be evoked by any stimulus, and networks able to generate amplified responses to at least one stimulus.

In the following section we examine OFF dynamics in networks with low-rank connectivity of the form given by Equation (3):

We first show that such connectivity directly constraints the network dynamics to a low-dimensional subspace. We then derive the conditions for the stability and amplification of OFF responses. Finally we apply these results to the specific case of low-rank rotational channels as in Equation (5).

Here we study the dynamics of the population activity vector for a unit-rank network ($R=1$) and for a general network with rank-$R$ connectivity structure. We consider low-rank networks in which the initial state is set to $\mathbf{𝐫}(0)={\mathbf{𝐫}}_0$ and no external input acts on the system at later times $t>0$. By construction, the autonomous dynamics generated by low-rank networks are constrained by the rank $R$ of the connectivity matrix and are therefore low-dimensional when $R\ll N$.

We first illustrate the linear dynamics in the case of a unit-rank connectivity ($R=1$), given by

with ${\mathbf{𝐯}}^{(1)T}$ of unit norm. The vectors ${\mathbf{𝐮}}^{(1)}$ and ${\mathbf{𝐯}}^{(1)}$ are respectively the right and left eigenvectors of $\mathbf{𝐉}$, and in the following analysis we refer to them respectively as the right and left connectivity vectors. In this case, the dynamics following the initial state ${\mathbf{𝐫}}_0$ can be explicitly computed as the product between the time-dependent propagator of the dynamics, given by ${\mathbf{𝐏}}_t=\exp (t(\mathbf{𝐉}-\mathbf{𝐈}))$, and the initial condition ${\mathbf{𝐫}}_0$ (Arnold, 1973; Bondanelli and Ostojic, 2020):

By expanding the exponential matrix $\exp (t(\mathbf{𝐉}-\mathbf{𝐈}))$ in power series, we can write the propagator for the unit-rank dynamics as:

where $\lambda$ is the only nonzero eigenvalue of $\mathbf{𝐉}$, given by

As a result, the full dynamics in the case of rank-1 connectivity structure can be written as:

where

Since ${\mathbf{𝐮}}^{(1)}$ is the right eigenvector of $\mathbf{𝐉}$ corresponding to $\lambda$, from Equation (40) we note that, when ${\mathbf{𝐫}}_0$ is fully aligned with ${\mathbf{𝐮}}^{(1)}$ the dynamics are one-dimensional and exhibit a monotonic decay along the same dimension. Instead, when the initial state is not fully aligned with ${\mathbf{𝐮}}^{(1)}$, the dynamics are confined to the plane defined by ${\mathbf{𝐫}}_0$ and ${\mathbf{𝐮}}^{(1)}$. In this case, while the component of the dynamics along ${\mathbf{𝐫}}_0$ decays exponentially as a function of time, the component along the direction of ${\mathbf{𝐮}}^{(1)}$ increases initially in proportion to the value of its norm, $||{\mathbf{𝐮}}^{(1)}||$ (since at time ${\displaystyle t=0}$, $\mathrm{d}{e}^{-t}a(t)/\mathrm{d}t=||{\mathbf{𝐮}}^{(1)}||$). Eventually, the activity decays to its asymptotic value given by $\mathbf{𝐫}(t\to \mathrm{\infty })=\mathrm{}$. Therefore in a unit-rank network the dynamics draw a neural trajectory that explores at most two-dimensions in the state space.

These observations can be extended to the general case of rank-$R$ connectivity matrices. For simplicity we rewrite Equation (3) as

where the matrices $\mathbf{𝐔}$ and $\mathbf{𝐕}$ contain respectively the $R$ right and left connectivity vectors as columns. Writing the connectivity matrix in this form is always possible by applying the SVD on $\mathbf{𝐉}$. The SVD allows us to write the connectivity matrix as $\mathbf{𝐉}=\widehat{\mathbf{𝐔}}{\mathbf{𝐒𝐕}}^T$, where $\mathbf{𝐔}=\widehat{\mathbf{𝐔}}\mathbf{𝐒}$, and ${\widehat{\mathbf{𝐔}}}^T\widehat{\mathbf{𝐔}}={\mathbf{𝐕}}^T\mathbf{𝐕}=\mathbf{𝐈}$. In particular, this implies the norm of each left connectivity vector ${\mathbf{𝐯}}^{(r)}$ is unity, while the right connectivity vectors ${\mathbf{𝐮}}^{(r)}$ are not normalized.

Following steps similar to Equation (38), the linear dynamics evoked by the initial state ${\mathbf{𝐫}}_0$ can be written as

where we defined the $R$-dimensional column vector

in analogy with Equation (40). Therefore, in the case of rank-$R$ connectivity matrix, the dynamics evolve in a $(R+1)$-dimensional space determined by the initial state ${\mathbf{𝐫}}_0$ and the $R$ right connectivity vectors ${\mathbf{𝐮}}^{(r)}$ (columns of $\mathbf{𝐔}$). Equation (44) shows that the dynamics of a rank-$R$ system are determined by the matrix of scalar products between left and right connectivity vectors, which we refer to as the overlap matrix (Schuessler et al., 2020; Beiran et al., 2020)

We conclude that low-rank connectivity matrices of the form given by Equation (3) with $R\ll N$ generate low-dimensional dynamics that explore at most $R+1$ dimensions.

In this paragraph we examine the conditions required to generate stable amplified dynamics in networks with low-rank connectivity $\mathbf{𝐉}$ and initial state ${\mathbf{𝐫}}_0$. Specifically, two sets of conditions need to be satisfied. The first set of conditions is directly derived by applying the general criterion for stable amplified dynamics given by Equation (34) to low-rank networks, which effectively constrains the connectivity $\mathbf{𝐉}$ through the relative arrangement of the connectivity vectors ${\mathbf{𝐮}}^{(r)}$ and ${\mathbf{𝐯}}^{(r)}$ ($r=1,\mathrm{\dots },R$). When this criterion is satisfied, amplified dynamics can be generated only if the initial condition ${\mathbf{𝐫}}_0$ is aligned with specific directions in the state-space. We thus examine a second set of conditions on the initial state ${\mathbf{𝐫}}_0$ for which amplified trajectories can be evoked, and express these conditions in terms of relationship between ${\mathbf{𝐫}}_0$ and the modes ${\mathbf{𝐮}}^{(r)}$-${\mathbf{𝐯}}^{(r)}$. As before, without loss of generality, we assume that the norm of the vectors ${\mathbf{𝐯}}^{(r)}$ is equal to one, while the norm of the vectors ${\mathbf{𝐮}}^{(r)}$ can vary.

We first consider a network with unit-rank connectivity $\mathbf{𝐉}={\mathbf{𝐮}}^{(1)}{\mathbf{𝐯}}^{(1)T}$, with ${\mathbf{𝐯}}^{(1)}$ of unit norm. In unit-rank networks, the dynamics is stable only if the nonzero eigenvalue $\lambda$ is smaller than one. From Equation (39) this yields the stability condition 

The condition for the generation of amplified responses given by Equation (34) can be derived in terms of ${\mathbf{𝐮}}^{(1)}$ and ${\mathbf{𝐯}}^{(1)}$ by computing the eigenvalues of the symmetric part of the connectivity, ${\mathbf{𝐉}}_S=({\mathbf{𝐮}}^{(1)}{\mathbf{𝐯}}^{(1)T}+{\mathbf{𝐯}}^{(1)}{\mathbf{𝐮}}^{(1)T})/2$. The matrix ${\mathbf{𝐉}}_S$ is of rank two, and has in general two nonzero eigenvalues, given by:

Therefore the condition for amplified OFF responses Equation (34) can be written in terms of the length of the vector ${\mathbf{𝐮}}^{(1)}$ as

Previous work has shown that amplification in a unit-rank network can take arbitrarily large values only if $0\le \cos \theta <1/||{\mathbf{𝐮}}^{(1)}||$ (Bondanelli and Ostojic, 2020). As a result, stable and unbounded amplified dynamics can be generated in a unit-rank network if and only if the norm of ${\mathbf{𝐮}}^{(1)}$ is sufficiently large (Equation (48)) and the correlation between the connectivity vectors is positive and sufficiently small (Equation (46)).

The conditions for stability and amplification derived for a rank-1 network can be implicitly generalized to a rank-$R$ network model, using the fact that the nonzero eigenvalues $\lambda$ of the low-rank connectivity matrix $\mathbf{𝐉}={\mathbf{𝐔𝐕}}^T$ (where $\mathbf{𝐉}$ is a $N\times N$ matrix) are equal to the nonzero eigenvalues of its overlap matrix ${\mathbf{𝐉}}^{ov}={\mathbf{𝐕}}^T\mathbf{𝐔}$ (where ${\mathbf{𝐉}}^{ov}$ is a $R\times R$ matrix) (Nakatsukasa, 2019):

Therefore, the condition for stability can be implicitly written in terms of the nonzero eigenvalues $\lambda$ of ${\mathbf{𝐉}}^{ov}$ as:

To derive the conditions for transient amplification, we need that the eigenvalues of the symmetric part of the connectivity be larger than unity, where the symmetric part ${\mathbf{𝐉}}_S$ is a matrix of rank $2R$ given by:

Using Equation (49) (Nakatsukasa, 2019), we can write the nonzero eigenvalues ${\lambda }_S$ of ${\mathbf{𝐉}}_S$ as:

where we used the fact that ${\mathbf{𝐕}}^T\mathbf{𝐕}=\mathbf{𝐈}$, and defined ${\mathbf{𝐒}}^2=\mathrm{diag}(||{\mathbf{𝐮}}^{(1)}||,||{\mathbf{𝐮}}^{(2)}||,\mathrm{\dots })$. By definition, the eigenvalues of the symmetric part multiplied by 2, that is, $2{\lambda }_S$, should satisfy:

where we used the expression of the determinant of a $2\times 2$ block matrix (Horn and Johnson, 2012). Since for two square matrices $\mathbf{𝐀}$ and $\mathbf{𝐁}$ we have that $det(\mathbf{𝐀𝐁})=det(\mathbf{𝐀})det(\mathbf{𝐁})$, $det({\mathbf{𝐀}}^T)=det(\mathbf{𝐀})$, and $det({\mathbf{𝐀}}^{-1})=1/det(\mathbf{𝐀})$, we can write: