An emerging paradigm proposes that neural computations can be understood at the level of dynamic systems that govern low-dimensional trajectories of collective neural activity. How the connectivity structure of a network determines the emergent dynamical system, however, remains to be clarified. Here we consider a novel class of models, gaussian-mixture, low-rank recurrent networks in which the rank of the connectivity matrix and the number of statistically defined populations are independent hyperparameters. We show that the resulting collective dynamics form a dynamical system, where the rank sets the dimensionality and the population structure shapes the dynamics. In particular, the collective dynamics can be described in terms of a simplified effective circuit of interacting latent variables. While having a single global population strongly restricts the possible dynamics, we demonstrate that if the number of populations is large enough, a rank network can approximate any -dimensional dynamical system.
A newly emerging paradigm posits that neural computations rely on collective dynamics in the state-space corresponding to the joint activity of all neurons in a network (Churchland & Shenoy, 2007; Rabinovich, Huerta, & Laurent, 2008; Buonomano & Maass, 2009; Saxena & Cunningham, 2019; Vyas, Golub, Sussillo, & Shenoy, 2020). Experiments in behaving animals have found that trajectories of neural activity are typically restricted to low-dimensional manifolds in that space (Machens, Romo, & Brody, 2010; Mante, Sussillo, Shenoy, & Newsome, 2013; Rigotti et al., 2013; Gao, Ganguli, Battaglia, & Schnitzer, 2015; Gallego et al., 2018; Chaisangmongkon, Swaminathan, Freedman, & Wang, 2017; Wang, Narain, Hosseini, & Jazayeri, 2018; Sohn, Narain, Meirhaeghe, & Jazayeri, 2019) and can therefore be described by a small number of collective, latent variables. It has been proposed that these collective variables form dynamical systems that implement computations through their responses to inputs (Eliasmith & Anderson, 2003; Hennequin, Vogels, & Gerstner, 2014; Rajan, Harvey, & Tank, 2016; Remington, Egger, Narain, Wang, & Jazayeri, 2018; Remington, Narain, Hosseini, & Jazayeri, 2018). How synaptic connectivity shapes the effective dynamics of collective variables, and therefore computations, however remains to be clarified.
Recurrent neural networks (RNNs) trained to perform neuroscience tasks are an ideal model system to address this question and further develop the theory of computations through dynamics (Sussillo, Churchland, Kaufman, & Shenoy, 2015; Rajan et al., 2016; Barak, 2017; Wang et al., 2018; Yang, Joglekar, Song, Newsome, & Wang, 2019). A recently introduced class of models, low-rank RNNs, directly embodies the idea of low-dimensional collective dynamics, opens the door to relating connectivity and dynamics, and provides a framework that unifies a number of specific RNN classes (Mastrogiuseppe & Ostojic, 2018). Low-rank RNNs rely on connectivity matrices that are restricted to be low rank, which directly generate low-dimensional dynamics. The rank of the network determines the number of collective variables needed to provide a full description of the collective dynamics. While previous works have shown that other specific classes of RNNs can approximate arbitrary dynamical systems (Doya, 1993; Maass, Joshi, & Sontag, 2007), the range of collective dynamics that can be implemented by low-rank RNNs remains to be clarified.
In this work, we focus on low-rank RNNs in which neurons are organized in distinct populations that correspond to clusters in the space of low-rank connectivity patterns. Each population is defined by its statistics of connectivity, described by a multivariate gaussian distribution, so that the full network is specified by a mixture of gaussians. The total number of populations in the network is a hyperparameter distinct from the rank of connectivity. Previous works have considered low-rank networks consisting of a single, global gaussian population (Mastrogiuseppe & Ostojic, 2018, 2019; Schuessler, Dubreuil, Mastrogiuseppe, Ostojic, & Barak, 2020). In the opposite limit, by increasing the number of populations, a gaussian mixture model can approximate any arbitrary low-rank connectivity distribution. Here we examine how the number of populations and their structure determine and limit the resulting collective dynamics in the network.
We first derive three general properties of gaussian-mixture low-rank networks: (1) in an autonomous network of rank , dynamics are characterized by collective variables that form a dynamical system; (2) the dynamics are determined by an effective circuit description, where collective variables interact through gain-modulated effective couplings; and (3) the resulting low-dimensional dynamics can approximate any arbitrary -dimensional dynamical system if the number of populations is large enough. We then illustrate how increasing the number of populations in a network extends its dynamical range. For that, we specifically focus on fixed points of the dynamics. While a network consisting of a single population can generate at most a pair of stable fixed points, independently of its rank, we show that adding populations allows the network to implement arbitrary numbers of stable fixed points embedded in a subspace determined by the rank of the connectivity matrix. Finally, we propose a general algorithm to approximate a given -dimensional dynamical system with a multipopulation network of rank and show one example network that is designed to implement complex temporal dynamics.
2 Model Class: Gaussian Mixture Low-Rank Networks
In this section, we introduce the class of models we study and define the key underlying quantities.
We assume that for each neuron, the set of pattern loadings is generated independently from a multivariate probability distribution . We moreover restrict ourselves to a specific class of loading distributions: mixtures of multivariate gaussians. This choice is motivated by the fact that gaussian mixtures can approximate any arbitrary multivariate distribution, afford a natural interpretation in terms of populations, and allow for a mathematically tractable and transparent analysis of the dynamics as shown below.
3 Dynamics in Gaussian Mixture Low-Rank Networks
In this section, we present three key properties of dynamics in mixture of gaussian low-rank networks: (1) in a network of rank , dynamics can be characterized by collective variables that form a dynamical system; (2) for loadings drawn from gaussian mixture distributions, the dynamics can be further described as an effective circuit in which collective variables interact through gain-modulated effective couplings; and (3) with a sufficient number of populations, the resulting low-dimensional dynamics can approximate an arbitrary -dimensional dynamical system.
Details of the derivations are provided in appendixes A and B.
3.1 Low-Dimensional Dynamics
We focus in the following on networks receiving a constant input, so that there is only one collective variable along the input dimension, the value of which is constant. The recurrent connectivity contributes to the dynamics of through the term .
The dynamics of collective variables in equation 3.2 are valid for any finite-size, low-rank network, without any assumption on the values of pattern loadings. We next turn to networks where the pattern loadings are generated from specific distributions.
3.2 Dynamics in Multipopulation Networks
3.3 Universal Approximation of Low-Dimensional Dynamical Systems
Note that the universal approximation theorem does not state how many populations are required to implement a given dynamical system and does not provide an algorithm for finding the statistics of the different populations.
4 Dynamics in Networks with a Single Population
Having shown that a rank network with an arbitrary number of populations can approximate any -dimensional dynamical system, we now illustrate how having a small number of populations in contrast limits the possible dynamics.
We focus first on the case of networks consisting of a single gaussian population. For simplicity, we focus on autonomous networks with zero-mean connectivity patterns. Specifically, we show that independent of their rank, the range of dynamics such networks can implement is restricted. This case was previously studied for connectivities that combined a rank-one or rank-two structure and a full-rank random component (Mastrogiuseppe & Ostojic, 2018; Schuessler, Dubreuil et al., 2020). Here we provide an overview of those results and extend them to single-population networks of arbitrary rank. The fact that we focus on networks whose connectivity is low-rank allows us to provide a deeper analysis of the dynamics.
4.1 Fixed Points
When the eigenvectors of are not orthogonal to each other, the eigenvectors of are not eigenvectors of the linear stability matrix . However, the eigenvalues of are still given by equation 4.5 (see appendix D), so that the same stability properties hold (see appendix D): every fixed point is stable in the direction toward the origin, and the fixed point in the direction given by the largest eigenvalue is stable while the other ones become unstable.
In summary, if all eigenvalues of the covariance matrix are real and nondegenerate, only the pair of nontrivial fixed points corresponding to the largest eigenvalue is stable. All the other nontrivial fixed points of the dynamics are saddle points. This implies that low-rank networks consisting of a single gaussian population can have at most two stable fixed points independent of their rank.
4.2 Limit Cycles
Complex eigenvalues of the covariance matrix , if they exist, always appear in conjugate pairs. They lead to spiral dynamics around the origin, in the plane spanned by the real and imaginary part of the corresponding eigenvectors. If the real part of the complex eigenvalues is smaller than unity, , the spiral dynamics decay back to the origin. Otherwise, if , there is a limit cycle on the plane, around the origin. Similar to the case with only real eigenvalues of the covariance matrix, if the real part of the complex eigenvalue is larger than the real part of any other eigenvalue of , any trajectory will converge to the plane defined by the real and imaginary parts of the corresponding eigenvectors. On this plane, we then find that the limit cycle is stable.
In this analysis, equation 4.10 is derived for the particular covariance matrix in equation 4.6, which is the sum of an isotropic matrix (proportional to the identity) and an antisymmetric matrix. However, this equation is valid more generally for any connectivity matrix with a pair of complex eigenvalues (see appendix D). When the covariance matrix is not antisymmetric but still has complex eigenvalues, the limit cycle is no longer a circle but resembles an ellipse, while the frequency of oscillation appears to still be given by equation 4.10.
Complex eigenvalues can be combined with real eigenvalues in networks with rank larger than two. The same stability properties are kept: only the fixed point or the limit cycle generated by the eigenvector with the largest real part is stable. The unstable fixed points or limit cycles remain attractive within the dimensions spanned by the corresponding eigenvector (see Figure 9 in appendix D) for an example of a rank-three network combining an unstable limit cycle and two stable fixed points).
4.3 Slow Manifolds
In finite-size simulations, the sampling of random loadings introduces spurious correlations in the matrix , breaking the degeneracy of the eigenvalues. As a consequence, only a small number of points on the continuous attractor predicted by the mean-field theory give rise to actual fixed points. While the rest of the points on the predicted continuous attractor are not fixed points of the finite-size network, the dynamics around them are typically slow. More specifically, any trajectory of activity quickly converges toward the predicted continuous attractor and then slowly evolves along it until it reaches a fixed point (see Figure 3L) (Mastrogiuseppe & Ostojic, 2018). In finite-size networks, the continuous attractor predicted by the mean-field analysis therefore gives rise to a low-dimensional manifold in state space, along which the dynamics are slow.
When degenerate and nondegenerate real and complex eigenvalues are combined, the global stability appears to be given by the criterion in equation 4.5: each eigenvalue generates its corresponding nontrivial dynamics (fixed points, continuous attractors or limit cycle) independently. The stability of these dynamical phenomena depends on the global eigenspectrum: the eigenvalues with the largest real part generate stable attractors, while the other eigenvalues lead to repellers.
In a low-rank network consisting of a single gaussian population, the possible nontrivial steady states are a pair of fixed points, a limit cycle, or a continuous attractor that gives rise to a small number of fixed points in finite networks. On top of this limited range of stable solutions, increasing the rank leads to additional unstable fixed points and limit cycles that can potentially be used to control the dynamics, a point we do not further explore here. We instead proceed to show that increasing the number of gaussian populations allows networks to implement a larger range of stable dynamics.
5 Dynamics in Networks with Multiple Populations
As described in the previous section, a major limitation of rank- networks consisting of a single gaussian population is that they cannot give rise to more than two stable fixed points, symmetrically arranged around the origin. We next show that networks consisting of several gaussian populations can exhibit a larger number of stable fixed points. We specifically describe two mechanisms by which multiple fixed points can be generated and controlled and show that classical Hopfield networks (Hopfield, 1982) correspond to a particular limit of gaussian-mixture low-rank networks.
5.1 Nonlinear Gain Control
The first mechanism for generating multiple fixed points is based on having several populations that reach saturation in different regions of the collective space. The local dynamics in these different regions are then controlled solely by the statistics of the nonsaturated populations.
For concreteness, we consider a rank-one network consisting of two populations, defined by different statistics of pattern loadings. Within population , for , the joint distribution of and values over neurons is specified by a covariance matrix . For simplicity, we take the mean of the distribution to be zero. In the two-dimensional loading space defined by and , the two populations correspond to different gaussian clusters, both centered at zero but with different shapes and orientations (see the green and purple dots in Figure 4A).
In particular, the network can have three stable fixed points: one at the origin and a pair of symmetrical nontrivial fixed points. First, the origin is always a fixed point of dynamics in equation 5.1. The origin is, moreover, a stable fixed point if the effective feedback at zero, which is given by , is smaller than 1. Second, at large values of , the effective feedback should be positive to cancel the contribution of the leaky term and generate a nontrivial fixed point. Therefore, one of the populations, which we define to be the first one (), must have a strong, negative overlap, . Given equation 5.2, this implies that the gain of the positively correlated population two should be large, whereas the gain of the negatively correlated population one should be close to zero. A small gain is achieved in the first population by having a large value , so that the second condition reads . Figures 4C and 4D show the dynamics of such a network given by the mean-field equation and in finite-size networks.
More generally, with more than two populations, this mechanism can be extended to produce a larger number of stable fixed points in rank-one networks. The two key components of this mechanism are (1) an independent control of the gain of the different populations, so that the contribution of each population to the effective feedback takes place in different ranges of the collective variable , and (2) covariances of different signs, so that the effective feedback can flexibly take both positive and negative values in different ranges of . These mechanisms can also be applied to networks with rank higher than one. In that case, the overlap between loadings is given by a matrix instead of a scalar, while the gain of each population is a scalar value. Populations with different covariance matrices and gains that vary at different ranges of the collective variables are able to generate multiple fixed points in different regions of the collective space, or combinations between stable limit cycles and stable fixed points (Dubreuil, Valente, Beiran, Mastrogiuseppe, & Ostojic, 2020).
5.2 Symmetries in Loading Space
A second mechanism for generating multiple fixed points is to exploit symmetries in the distribution of loadings . Such a symmetry in the connectivity induces a symmetry in the dynamics of the collective variables. In consequence, if a network generates one nontrivial stable fixed point, additional stable fixed points appear at symmetric points in the collective space.
We focus here on networks where the overlap between the connectivity patterns is given by the nonzero means of the loadings, which is complementary to the previous section where the connectivity patterns had zero mean and the recurrent dynamics is determined by the covariances between the loadings. We introduce a symmetry in the distribution of loadings by arranging the means of the loadings accordingly. Note, however, that symmetrical distributions of loadings can also be generated in the zero-mean case.
As an illustration, we consider first a rank-two network, with units evenly split into populations. In each population, the loadings have a different set of means , , , , and the covariances are zero. The variance, of the loadings, and , are identical in all populations. As a consequence, different populations correspond to clusters of identical spherical shape, but centered at different points in the four-dimensional loading space.
Because of the symmetry, if , there are at least stable fixed points arranged symmetrically on a circle (see Figures 5A to 5C for a network with populations and Figures 6A to 6C for a network with populations). If the number of population pairs is odd, there are stable fixed points symmetrically arranged on a circle, because there is also a symmetry with respect to the origin, imposed by the symmetry in the transfer function. Otherwise, if is even, stable fixed points are generated by the network.
In general, any -dimensional discrete symmetry in the loadings will generate a dynamical system with stable fixed points on a -dimensional sphere, arranged with the symmetry of the dual polytope. A regular polytope is defined as the generalization of a regular polyhedron generalized to more than three dimensions.
5.3 Relation to Hopfield Networks
A rank-two network with four populations , characterized by equation 5.3 (see Figures 6A and 6C), is therefore equivalent to a two-pattern Hopfield network in the limit of no dispersion around the mean of each cluster, . In this limit, saddle points are located at the midpoints between neighboring stable fixed points. In the more general rank-two networks in equation 5.3, where , the saddle points between stable fixed points move farther away from the origin (such as in Figure 6B, where ), but the four stable fixed points remain on the vertices of a square along the axes and . In the limit of very large , the saddle points between stable fixed points approach the circle that circumscribes the stable fixed points.
The rank-three network presented in equations 5.10 and 5.11 also becomes a classical Hopfield network in the limit of . Allowing for values , as illustrated in Figures 6D and 6G, does not change the number of fixed points generated by the Hopfield network or their direction in collective space. These networks generate pairs of stable fixed points along the directions , , and . When is large, additional fixed points become stable along directions . These additional fixed points correspond to well-known spurious mixture states in Hopfield networks (Amit, Gutfreund, & Sompolinsky, 1987).
6 Approximating Dynamical Systems with Gaussian-Mixture Low-Rank Networks
In the previous section, we focused on generating multiple fixed points in an autonomous network by means of a few gaussian populations in the connectivity. More generally, as shown in section 2.3, multipopulation rank- networks can approximate any -dimensional dynamical system. In this section, we propose an algorithm for that purpose.
Previous works have developed algorithms for training recurrent networks to implement given dynamics that effectively used low-rank connectivity (Eliasmith & Anderson, 2003; Rivkind & Barak, 2017; Pollock & Jazayeri, 2020). These methods rely on tuning the loadings of individual neurons, given fixed external inputs and connectivity loadings . Here we focus instead on networks based on mixtures of gaussian populations, in which the couplings between individual neurons are not precisely set but instead sampled from a distribution. We extend previous methods to find the first- and second-order moments of multiple gaussian populations that approximate a given dynamical system.
Our algorithm proceeds as follows. We first fix the number of gaussian populations in the network and the fraction of neurons included in each population, . Depending on the complexity of the target dynamics and the required accuracy, a smaller or larger number of populations is required. Second, we set the mean and variance of the vectors in each population, and , together with the mean and variance of the external input, and . We assign these parameters randomly according to a certain probability distribution. Finally, we determine the statistics of the vectors, the only unknown in the network, using linear regression.
The number of populations, the level of regularization, together with the distributions chosen to fix the mean and covariance values , , , and are hyperparameters of the algorithm. These hyperparameters can be tuned progressively by running several iterations of the algorithm. For example, a possible goal is to search for the minimal number of populations required for approximating a given dynamical system within some accuracy limits. In general, we observe empirically that the distribution underlying the fixed parameters does not play a crucial role in the accuracy of the algorithm as long as they span a wide range of values.
We first analyze the performance of the algorithm as a function of the number of populations, for two different levels of regularization. We find that regularization is necessary to avoid diverging values of the parameters in (see Figure 7B, black versus red curve). This is due to the fact that the matrix is close to singular, so that its inverse reaches very high values. Interestingly, the error in the approximation increases only slightly with a strong level of regularization (see Figure 7C). As the number of populations is increased, the error in the approximation for all levels of regularization monotonically decreases.
The algorithm is based on the mean-field description in equation 6.2, which holds in the limit of a very large number of neurons per population. We next study whether the obtained mean-field networks describe well the dynamics in networks with a finite number of neurons per population. We therefore sampled the connectivity and input loadings from the multivariate gaussian distribution that characterizes each population. For that purpose, it is necessary to set the variances of the right connectivity patterns , which do not influence the mean-field dynamics. As a general approach, we set those variance' parameters to be as low as possible but high enough so that the correlation matrix of the multivariate gaussian distribution is positive-definite. Sampling a finite number of values from the distribution introduces deviations in the sampled mean and sampled covariance matrix that introduce additional finite-size errors in the approximated dynamics. We find that when the parameters obtained in reach very large values, the approximation error in networks with per population is very large (see Figure 7D, black curve). However, when the algorithm is constrained by a strong level of regularization, the finite-size error remains approximately constant as the number of populations increases (see Figure 7D, red curve).
We then show the approximated dynamics for two different numbers of populations, (see Figures 7E and 7F) and (see Figures 7G and 7H), in the regularized case. Fifteen populations are enough to obtain a limit cycle with similar features to the Van der Pol oscillator, although there are deviations in the shape and frequency of oscillation (see Figure 7E). The finite-size effects remain small when units are considered in each population, so that the dynamics of a simulated network resemble the mean-field approximation (see Figure 7F). For a larger number of populations, the mean-field approximation increases the accuracy considerably (see Figure 7G). However, finite-size effects increase with the number of populations, when the number of units per population is kept constant, so that the approximation error in finite-size simulations is not necessarily reduced (see Figure 7H).
The finite-size deviations that contribute to the approximation error depend on the number of units per population that we chose to fix in Figures 7D, 7F, and 7H. We expect a decrease in finite-size errors as increases, so that the network dynamics resemble more closely the mean-field description. Accordingly, with a high enough number of units per population , a finite-size network with populations (as in Figure 7F) could approximate better a Van der Pol oscillator than a finite-size network with populations (as in Figure 7H). A full characterization of finite-size effects is beyond the scope of this study.
This algorithm can be applied to generate any given dynamics in collective space within a finite domain. Beyond this finite domain sampled through the chosen set points, if the target vector field does not follow the required asymptotic behavior (see equation 3.14), as is the case for the Van der Pol oscillator, the network will not extrapolate to the target dynamics (region outside square in Figures 7E to 7H, left). However, in practice, it may produce qualitatively similar dynamics: in the example of the Van der Pol oscillator, if the network is initialized at a point outside the limit cycle, the resulting trajectories still converge to the limit cycle.
In this letter, we have examined the dynamics in gaussian-mixture low-rank recurrent neural networks, a class of models in which the connectivity is defined by a low-rank matrix, with connectivity patterns consisting of several populations with distinct gaussian statistics. In these networks, the collective dynamics can be described by collective variables, where is the rank of the connectivity matrix and the dimensionality of the input patterns. These collective variables form a dynamical system, the evolution of which is determined by the connectivity statistics of the populations forming the network. The rank of the network and the population structure therefore play different roles: the rank of the network sets the internal dimensionality of the dynamics and defines the corresponding collective variables, while individual populations shape the dynamics of these collective variables but do not contribute new ones. We specifically showed that in the limit of a large number of populations, this class of network displays a universal approximation property and can therefore implement a large range of dynamical systems. Having a small number of populations instead imposes constraints and limits the achievable range of dynamics.
We focused here on a specific family of distributions for the connectivity patterns, mixtures of multivariate gaussians. This choice was motivated by several considerations. First, this family of distributions can be used to approximate any multivariate distribution for the pattern loadings. Second, it leads to a particularly simple form of dynamics for the collective variables, where the time evolution is formulated in terms of a simple effective circuit (see equation 3.11). Remarkably, in this description of the dynamics, which is exact and nonlinear, the collective variables appear to interact linearly through effective couplings and effective inputs. This allows for a particularly transparent interpretation of dynamics in terms of gain modulation. Several of our results are, however, independent of the specific assumption for the type of distribution; this is in particular the case for the influence of symmetry in the connectivity on the dynamics. When a large number of populations is needed to approximate the connectivity structure, other parametric distributions may be more suitable, and the interpretation in terms of discrete populations may not be appropriate.
Previous work (Mastrogiuseppe & Ostojic, 2018, 2019; Schuessler, Dubreuil et al., 2020) has studied the generation of low-dimensional dynamics by recurrent neural networks with a low-rank component in the connectivity, focusing on fixed points and limit cycles. In these studies, the connectivity matrix contained a random full-rank term and a low-rank term, where the low-rank loadings are randomly drawn from a single gaussian population. The full-rank term in the connectivity does not allow for the expression of the emerging dynamics directly as a low-dimensional dynamical system, as we have presented in this work (see section 3), and the study of the potential constraints on the low-dimensional dynamics.
Low-rank networks with arbitrary pattern distributions form a rich and versatile framework that encompasses a number of previously studied types of recurrent neural networks. As shown in the last part of the results, Hopfield networks storing patterns can be seen as a particular limit of gaussian-mixture low-rank networks, in which pattern loadings are binary and exhibit a specific type of symmetry. The neural engineering framework (Eliasmith & Anderson, 2003) and the manifold embedding approach (Pollock & Jazayeri, 2020) provide algorithms that implement specific low-dimensional dynamics by controlling the structure of fixed points and Jacobians using linear-regression methods. These algorithms generate recurrent networks with low-rank connectivity, in which the individual couplings between individual neurons are set to specific values. In contrast, we focus on low-rank networks in which individual couplings are sampled from an underlying distribution, and our algorithm determines the statistics of this distribution rather than individual couplings. This approach automatically endows our networks with strong robustness with respect to the deletion of individual neurons.
The dynamics of recurrent networks of rate units that include several populations have been studied in recent years (Aljadeff, Stern, & Sharpee, 2015; Aljadeff, Renfrew, Vegué, & Sharpee, 2016; Kadmon & Sompolinsky, 2015). In these works, the connectivity consisted of full-rank, random matrices in which populations define the pair-wise connectivity statistics. For example, in a two-population model of excitatory and inhibitory neurons, the inhibitory population is defined by the connectivity statistics (mean and variance of synaptic strengths) between inhibitory neurons and the connectivity statistics toward excitatory neurons. In contrast, here we focus on low-rank structure in the connectivity matrices and specify populations in terms of the statistics for the low-rank structure vectors within each population. This allows us to go beyond random connectivity and study the effects of the rank and the number of populations independently. Directly relating the two descriptions of connectivity is a topic of ongoing work.
Our framework is also closely related to echo-state (Jaeger, 2001) and FORCE networks (Sussillo & Abbott, 2009), which rely on randomly connected recurrent networks controlled by feedback loops. Each feedback loop is mathematically equivalent to adding a unit-rank component to the connectivity matrix. Echo-state and FORCE networks therefore correspond to low-rank networks with an additionnal full-rank, random term in the connectivity (Mastrogiuseppe & Ostojic, 2018, 2019). Because the feedback loops are trained to produce specific outputs, the low-rank part of the connectivity is typically correlated to the random connectivity term (but see Mastrogiuseppe & Ostojic, 2019). Such correlations increase the dimensionality and the range of the dynamics (Schuessler, Dubreuil et al., 2020; Logiaco, Abbott, & Escola, 2019), although the low-rank connectivity structure and the number of populations still generate strong constraints. For instance, for rank-one networks with a random term in the connectivity but consisting of a single population, the fixed points are restricted to lie on a one-dimensional but nonlinear manifold, and typically at most two nontrivial stable fixed points can be generated (Schuessler, Dubreuil et al., 2020). More generally, random components in the connectivity can strongly influence learning dynamics during training (Schuessler, Mastrogiuseppe, Dubreuil, Ostojic, & Barak, 2020).
Gaussian-mixture low-rank networks, the neural engineering framework, and echo-state networks all exhibit universal approximation properties (Eliasmith, 2005; Maass, Natschläger, & Markram, 2002). It is, however, important to distinguish between several variants of this property. In our case, in analogy with the NEF, we started from an -dimensional dynamical system fully specified by its flow function and showed that gaussian-mixture low-rank networks can approximate this flow function, provided a large number of populations is available and the flow function satisfied specific constraints. Echo-state and FORCE networks instead start by specifying a target readout, and universal approximation means that any such readout can be generated by training the feedback (Maass et al., 2007). This readout corresponds to a low-dimensional projection of a large dynamical system, and echo-state networks are free to implement any dynamical system consistent with the specified output projection. This is a major distinction with our approach and that of the NEF, where the overall dynamical system is more tightly constrained.
In this work, we have examined only networks with fixed inputs. Varying the inputs instead modifies the low-dimensional dynamics, an effect that can be understood through modulations of effective couplings that govern the interactions between collective variables. In a companion paper (Dubreuil et al., 2020), we have used gaussian-mixture, low-rank RNNs to reverse-engineer networks trained on a range of neuroscience tasks and found that gain modulation through input control underlies complex computations, such as flexible input-output mappings (Fusi, Miller, & Rigotti, 2016). Varying inputs while keeping connectivity fixed therefore has the potential of implementing a large range of dynamical systems and computations (Pollock & Jazayeri, 2020), but the full capacity of this mechanism still remains to be fully elucidated.
Appendix A: Dynamics in Multipopulation Networks
In this appendix, we derive the equation for the dynamics of a multipopulation, low-rank network, equation 3.4. We consider a low-rank network that consists of populations, where each population is defined by different statistics of the probability distribution . We assume that the external input is constant in time and uncorrelated with the left connectivity patterns at the level of each population. Each neuron in the network is assigned to a population according to the probability . In the following, we set the statistics of each population to be drawn from a multivariate gaussian with mean vector , as defined in equation 2.5 and covariance matrix (see equation 2.6).
Appendix B: Universal Approximation of Low-Dimensional Dynamics
There is a direct mapping between the second term of equation B.1 and the recurrent dynamics of a low-rank RNN. The recurrent dynamics in equation 3.3 can be directly mapped to equation B.1: the variables correspond to , to , and to . This implies that the recurrent dynamics can approximate any flow function within a finite domain. The parameters , , and can be adjusted independently. In particular, is independent of all the other collective variables , although the opposite is not true: the collective variables obviously depend on the external tonic input . We are assuming that the temporal profile of the background input is constant, so that is also constant after discarding a possible initial transient (see equation 3.2).
so that maps to , maps to , and is mapped to the bias term . The transfer function is, however, different. In equation B.1, the nonlinear function used is , while in equation B.2, the nonlinear function used is . Both functions are nonlinear and nonpolynomial, so that the theorem applies in each case. The contribution given by the disorder in the population loadings, and , are not required for the universal approximation. However, quadratic terms like the one introduced by the variance of loadings improve the approximation in terms of expressibility and efficiency (Fan, Xiong, & Wang, 2020). Overall, this means that a low-rank network with a finite number of populations can approximate any dynamical system within a bounded domain.
Appendix C: Linear Stability Matrix at Fixed Points in Networks with a Single Population
Furthermore, it can be shown analytically that is never zero for any finite value , by studying the minima of its primitive function. The primitive function of is proportional to . The primitive function has no local minima because it is bounded between zero and one, that is, a monotonically decreasing function of . We can show that the function is monotonically decreasing because for sigmoidal transfer functions, for any fixed value , if and only if . Thus, this property is still conserved when calculating the gaussian average: if and only if .
Putting these analyses together, we conclude that the function is for , is always smaller than zero and tends asymptotically to this upper bound as approaches infinity. This result is used to study whether the linearized dynamics around fixed points in low-rank networks with a single gaussian population are stable (see equation 4.5).
Appendix D: Stability Analysis of Rank-Two Networks with Nonnormal Covariance
In section 4, we analyzed the dynamics generated by a rank-two network with one single gaussian population when the covariance matrix is normal, that is, the eigenvectors are mutually orthogonal to each other. We extend here the analysis to nonnormal matrices and show that the main features of the dynamics are conserved when the correlation matrix is nonnormal.
We first studied the case of normal matrices with real eigenvalues (see Figures 3A and 3D and equations 4.2 to 4.5). The analysis showed that each real eigenvector of the covariance matrix leads to a pair of fixed points when the associated eigenvalue is larger than one.
The linear dynamics around fixed points is given by the Jacobian in equation 4.4, where . The eigenvectors of coincide with the eigenvectors of if the covariance matrix has real eigenvectors mutually orthogonal to each other. Using this property, we can then determine the eigenvalues of the linearized dynamics around each fixed point in equation 4.5. As a conclusion, this analysis showed that all fixed points are saddle points, except for the two fixed points in the direction of the eigenvector with the largest associated eigenvalue, where the fixed points are stable.
Figures 8A and 8C show how the dynamics vary when the correlation matrix is nonnormal but the eigenvalues remain constant. The parameter determines the degree of correlation between eigenvectors: as increases, the two eigenvectors become more strongly correlated (see Figure 8A). This correlation moves the direction of the associated fixed point, while keeping the radial distance constant (see Figure 8B). The eigenvalues of the linear dynamics around the fixed points do not vary. Figure 8C shows an example of the dynamical landscape in collective space. It is interesting to note that the trajectories that go from the saddle points toward the stable fixed points resemble ellipses, similar to the case of normal (see Figures 3A and 3D). However, the stable fixed points are no longer located at the farthest point of this curve from the origin.
In summary, we showed that the number of fixed points and the stability in a low-rank network with one single population depends only on the eigenvalues of , and not on its eigenvectors. When the eigenvectors are all mutually orthogonal, we can calculate analytically the eigenvalues and eigenvectors of the linearized dynamics around each fixed point. When the eigenvectors are not all mutually orthogonal, the eigenvectors of the linearized dynamics change, but not the eigenvalues. This suggests that increasing the correlation between eigenvectors introduces a continuous deformation of the collective space while the timescales of the linearized dynamics around fixed points are preserved.
We observed numerically that the previous result holds as well for covariance matrices with complex eigenvalues. Each pair of complex conjugate eigenvalues with eigenvector generates a limit cycle on the plane spanned by –h if . If is a normal matrix, such that vectors and have the same norm and are orthogonal to each other, the limit cycle is a circle. In that case, we showed in section 4 that the frequency of the limit cycle is given by .
Adding correlations between the real and imaginary part of the eigenvectors and , or changing their relative norm, will change the shape of the limit cycle to a curve resembling an ellipse. Figures 8D an 8F show an example of a covariance matrix where the nonnormality is controlled by parameter without affecting the eigenvalues. As increases, the relative norm of the imaginary and real part of the eigenvector changes (see Figure 8D). This correlation changes the shape of the limit cycle, from a circle when is one (normal case; see the gray line in Figure 8E) to a closed curve resembling an ellipse as increases (see Figures 8E and 8F, the gray line). The frequency of oscillation along the limit cycle, however, stays invariant as is varied (see Figure 8E bottom).
To sum up, also in the case of complex eigenvalues, the correlations between eigenvectors deform the collective space, introducing a continuous mapping between the circular limit cycle and the new limit cycle, which resembles an ellipse. However, the temporal features of the dynamical landscape, such as the frequency of the limit cycle, remain constant.
Real and complex eigenvalues are combined in networks with rank three or larger. The same stability rules hold: the dynamic structure (limit cycle or fixed point) generated by the eigenvalue with the largest real part is stable, while the other ones are not stable in all directions. Figures 9A and 9B show an example of a rank-three network, whose connectivity matrix has a real eigenvalue and a pair of complex conjugate eigenvalues and . The real part of all eigenvalues is larger than one, so the real eigenvalue leads to a pair of fixed points, and the complex eigenvalues generate a limit cycle. Given that the real eigenvalue is larger than the real part of the other eigenvalues, the fixed points are stable. The limit cycle is marginally stable in the plane spanned by the real and imaginary parts of the complex eigenvector of but unstable in any other direction. Therefore, trajectories starting in the plane converge to the limit cycle in the mean-field equation (see the gray trajectory in Figure 9C). Small perturbations, such as those introduced by finite-size effects, make these trajectories deviate from the limit cycle and converge to one of the two stable fixed points (see the gray trajectories, Figure 9D).
The project was supported by the Ecole de Neurosciences de Paris, the ANR project MORSE (ANR-16-CE37-0016), the CRCNS project PIND, and the program Ecoles Universitaires de Recherche launched by the French government and implemented by the ANR, with the reference ANR-17-EURE-0017. There are no competing interests. We thank Mehrdad Jazayeri and Eli Pollock for discussions.
Code and trained models are available at https://github.com/emebeiran/low-rank2020.