Functional connectivity modules in recurrent neural networks: function,
  origin and dynamics
Understanding the ubiquitous phenomenon of neural synchronization across
species and organizational levels is crucial for decoding brain function.
Despite its prevalence, the specific functional role, origin, and dynamical
implication of modular structures in correlation-based networks remains
ambiguous. Using recurrent neural networks trained on systems neuroscience
tasks, this study investigates these important characteristics of modularity in
correlation networks. We demonstrate that modules are functionally coherent
units that contribute to specialized information processing. We show that
modules form spontaneously from asymmetries in the sign and weight of
projections from the input layer to the recurrent layer. Moreover, we show that
modules define connections with similar roles in governing system behavior and
dynamics. Collectively, our findings clarify the function, formation, and
operational significance of functional connectivity modules, offering insights
into cortical function and laying the groundwork for further studies on brain
function, development, and dynamics.
Jacob Tanner, Sina Mansour L., Ludovico Coletta, Alessandro Gozzi, Richard F. Betzel
# Functional connectivity modules in recurrent neural networks: function, origin and dynamics

###### Abstract

Understanding the ubiquitous phenomenon of neural synchronization across species and organizational levels is crucial for decoding brain function. Despite its prevalence, the specific functional role, origin, and dynamical implication of modular structures in correlation-based networks remains ambiguous. Using recurrent neural networks trained on systems neuroscience tasks, this study investigates these important characteristics of modularity in correlation networks. We demonstrate that modules are functionally coherent units that contribute to specialized information processing. We show that modules form spontaneously from asymmetries in the sign and weight of projections from the input layer to the recurrent layer. Moreover, we show that modules define connections with similar roles in governing system behavior and dynamics. Collectively, our findings clarify the function, formation, and operational significance of functional connectivity modules, offering insights into cortical function and laying the groundwork for further studies on brain function, development, and dynamics.

## I Introduction

The brain is a complex adaptive network with dynamically interacting parts that somehow come together to form the locus of control for our behavior [1; 2; 3]. In the quest to unravel the intricacies of brain function, one observation that has garnered considerable attention is the large-scale coordination of neurons [4; 5; 6]. The synchronization of neural activity across neuronal populations is a pervasive phenomenon, observed across multiple species and levels of neural organization. From the coordinated firing of neurons in cellular assemblies in fish [7; 8; 9], mice [10; 11; 12] and monkeys [13], to large-scale systems like the default mode network [14; 15; 16; 17], it is likely that synchronization plays a pivotal role in brain function. To explore this, correlation-based network approaches serve as an important tool for probing the statistical relationships between neural units at different spatiotemporal scales[18; 19].

Analyses of correlation-based networks often reveal modular structures, characterized by sets of neural units whose activity is highly correlated [20; 21; 22; 23; 24]. Previous studies have shown that the boundaries of modules in correlation networks circumscribe meta-analytic task co-activation patterns [25; 26] and early investigations suggest that correlation-based network modules emerge in the fetal brain and align with areas that will later support vision, movement, and language [27; 28]. In addition, research indirectly related to the origin and function of correlation-based modules at the micro-scale of neurons suggest that modules can form from activity-based plasticity mechanisms [11], and that stimulating these modules can elicit relevant task behavior [10]. Additionally, previous work has modeled this correlation structure as emerging from dynamic interaction with the underlying structural connectome [29; 30; 31].

These correlation-based networks are often referred to as _functional connectivity_ (FC) networks, a term that has sparked some controversy given that the name suggests that these correlations imply a functional relationship (e.g. [32; 33; 34]). This is a critique that can doubly be leveled at the analysis of modules in functional connectivity networks, given the temptation to relate them to the rich history of theory and research on the function of modules in structural networks [35; 36; 37; 38; 39; 40; 41; 42; 43]. In this way, the functional relationship implied by functional connectivity modules is contested. Additionally, although previous research has tracked modules in functional networks across development and, indeed, the entire human lifespan [44; 45; 46], the specific origin of functional connectivity modules as well as their specific relationship with the dynamical features of a system, remain unclear.

Recurrent neural networks offer trainable systems that operate based on a network of dynamically interacting parts and as such are ideal model organisms for investigating such questions [47; 48; 49; 50; 51; 52; 53; 54; 55]. Not only do these artificial systems provide us with complete access to all the information important for their function, but they also offer a
safe and ethically neutral platform for perturbation; an important tool for revealing the underlying causal relationships in a system. While the technology for recording and perturbing neural systems in living brains has been advancing at an increasing pace (e.g. [56; 57; 58; 59; 60]), complete access to _all_ of the information important to the function of a living brain is still an aspirational, and arguably distant goal, and the ethical implications of perturbing these systems remain complex.

Our study leverages recurrent neural networks trained on canonical systems neuroscience tasks [61] to investigate three critical questions about functional connectivity modules: 1) What specific functional role do functional connectivity modules serve? 2) What drives their formation? 3) What insights can they provide into system dynamics? Importantly, our research suggests these modules are not merely statistical artifacts but are functional units that encapsulate and specify dynamics that uniquely transform similar types of information. Additionally, our findings indicate that one origin of empirically observed FC modules may lie in the asymmetries of input projections. Finally, this result hints that input projections, particularly from subcortical areas like the thalamus, could serve as a genetically efficient means to encode module specialization in the cortex. Collectively, our findings offer key insights into functional connectivity modules in recurrent neural networks and lay the groundwork for future studies to take advantage of these insights in studying biological brains.

## Results

The following results are organized into three sections, speaking to the 1) function, 2) origin, and 3) dynamics of functional connectivity (FC) modules in recurrent neural networks. Each section contains subsections that support the general conclusions for that section. The conclusions for each section are as follows. The _function_ of modules in recurrent neural networks is to hold onto task relevant information. An important _origin_ of these functional modules involves asymmetries in the input projections. Finally, the _dynamics_ of these functional connectivity modules involves the transformation of input information into task relevant information.

### Function

In this section, we ask the following question: what is the _function_ of FC modules in recurrent neural networks (RNNs)? That is, how do functional modules, specifically, contribute to the RNNs ability to perform its prescribed task? In answering this question, we show that FC modules in RNNs accumulate task relevant information, which the output layer then reads out in order for the network to make its decision. We also show that the organization of these modules can be used to identify semantically similar information in feed-forward neural networks.

#### Functional connectivity modules take on task-relevant functional significance

In this section we show that FC modules in the recurrent layer of RNNs take on task-relevant functional significance after training on a perceptual decision-making task and a go _vs_ no-go task. These two tasks are common systems neuroscience tasks used to probe information transformation, decision-making, and memory [10; 55; 61]. Additionally, we extend this result by showing a similar effect in feed-forward neural networks trained on the recognition of hand-written digits, and a transformer deep neural network previously trained by openAI (weights from GPT-2 [62]).

The perceptual decision-making task presents the RNN with two stimuli drawn from normal distributions with different means (same variance). The task of the RNN is to compute which of the two stimuli come from the distribution with the greater mean. This task requires the RNN to track previous values of both stimuli and to compare them. After a fixation period, the RNN makes its decision based on the relative activity of two output nodes. The output node with the greatest activity corresponds to the RNNs decision about which stimuli came from the distribution with the greater mean (see Fig. 1a for a schematic of this task).

Here, we found that FC modules in the recurrent layer of these RNNs store information about the difference between the means of both these stimuli across the fixation period. We found that, before the dynamics in the recurrent layer, the modules defined by the input layer carry information about the current value of each stimulus (Fig. 1e,g; \(r=0.77,p<10^{-15}\); \(r=0.54,p<10^{-15}\)). After the recurrence however, this information is transformed into FC modules that carry information about the current difference between the cumulative means of each stimulus during the fixation period (Fig. 1f,h; \(r=0.78,p<10^{-15}\); \(r=-0.86,p<10^{-15}\)). Indeed, the activity of the output nodes that make the decision for this network are highly correlated with the mean activity of these FC modules (Fig. 1b-d; \(r=0.77,p<10^{-15}\); \(r=0.80,p<10^{-15}\)) suggesting that these modules are causally involved in the decision-making process.

To verify this, we performed a lesioning analysis wherein we selectively removed connections from the output layer that received information from neurons within a module ("output lesions"; Fig. 1i). This lesioning effectively removes the information within a functional connectivity module from consideration during the decision process. By performing an output lesion on module 1, we saw that task accuracy dropped to nearly 0% for trials where stimulus 1 was the correct decision, while trials where stimulus 2 was the correct decision were not effected (Fig. 1j; two-sample \(t\)-test \(p<10^{-15}\)). In con
trast, when we performed an output lesion on module 2 task accuracy was high for trials where stimulus 1 was the correct decision, while dropping to nearly zero for trials where stimulus 2 was the correct decision (Fig. 1k; two-sample \(t\)-test \(p<10^{-15}\)). The specificity of the effects of these output lesions on different trial types confirms that the information from these FC modules is being used by the network to make its decisions.

We find a similar result with an RNN trained to perform a go _vs_ no-go task. With this task, the RNN receives either a go signal, or a no-go signal, and after some delay must indicate which signal it received (see Fig. S1a for a schematic of this task). We find that the FC of the recurrent layer during this task maintains four modules (Fig. S1b) whose activity changes according the which signal was presented (Fig. S1c). Similar to the previous task, we found that the mean activity with the first two FC modules is highly correlated to the output activity in the RNN where the final decision is made by the network (Fig. S1d; \(r=0.72\); \(r=0.73\)). Finally, output lesions

Figure 1: **Functional connectivity modules take on task-relevant functional significance**, (_a_) Schematic describing the perceptual decision-making task and the architecture of the RNN trained to perform it. The input nodes are given stimulus information and a fixation input. The stimuli come from two distributions with different means. The RNN must determine which of the two stimuli come from the distribution with the greater mean while the fixation input has a value of 1. When the fixation input is zero, the decision is made based on which of the output nodes have the greater activity. (_b_) Functional connectivity of the recurrent layer of an RNN trained on the perceptual decision-making task reorganized according to a modular partition found using modularity maximization with Louvain. This algorithm found four modules, two of which are labeled ’module 1’ and ’module 2’. The other two modules are associated with the fixation input (one activates at the beginning of the fixation period, and the other activates at the end of the fixation period). (_c_) Mean activity of each module for different fixation periods as well as the difference in the cumulative mean between the stimuli. Notice how the mean activity in each module tracks with this value. (_d_) Two plots showing the correlation between activity in two of the modules and the output activity of the neurons that make the decision for the RNN.(_e_) Input projections create functional connectivity modules. (_f_) Functional connectivity modules are transformed by the recurrent layer. (_g_) Two plots showing that functional connectivity modules created by input projections are related to the current input stimuli values. (_h_) Two plots showing that the recurrent layer of the RNN transforms the information in modules so that they represent the cumulative stimulus difference. (_i_) Schematic showing the process of information lesioning to modules. Information lesions to modules were performed by removing weights in the output layer that sent information from nodes in that module. (_j_/_k_) Two sets of boxplots showing the accuracy and loss on different trials following an information lesion to module 1/2. Accuracy is considered separately for trials where stimulus 1 had the higher mean than stimulus 2, and vice versa.
to one of the modules significantly dropped the accuracy for trials of the task when the go signal was given, but not for trials where the no-go signal was given to the network. This suggests that this module holds information regarding the presence of the go signal(Fig. S1e-g; two-sample \(t\)-test \(p<10^{-15}\)).

Taken together, these results suggest that FC modules track task-relevant information that is then used by the network to make decisions. With both RNNs trained to perform a perceptual decision-making task, as well as RNNs trained to perform a go _vs_ no-go task, the mean activity of FC modules was found to correlate highly with the activity of output neurons (where the decision of the network is made). In addition, in both networks output lesions to these modules had specific effects on decision-making related to the information they were tracking.

In addition to our analysis of the function of FC modules in recurrent neural networks, we performed a supplemental analysis to investigate whether or not FC modules played a similar role in feed-forward neural networks. We found that the specific organization of neurons into FC modules in the feed-forward neural network would carry information about the categorical or semantic content of the input (see **Supplemental Section** 1, and Figs. S2 and S3 for more details on our analysis).

### Origin

In the previous section we investigated the function of modules in RNNs, showing that modules carry task-relevant information and relay this information to the output layer. Here, we investigate the origin and development of these modules. First, we demonstrate that FC modules emerge spontaneously in random networks and show that these modules emerge due to asymmetries in the sign and weight of projections from the input layer (input projections) that can be approximated using cosine similarity. Furthermore, we use this approximation to demonstrate a relationship between FC modules in the cortex of mice and humans and the cosine similarity of connections from the thalamus. Finally, we investigate the potential developmental role of initial asymmetries in input projections to the cortex from areas like the thalamus showing that initial weights can be used to guide the development of FC modules across learning.

Input projections influence the modular structure of functional connectivity in both synthetic and real brains

In previous sections, we showed that FC modules in artificial neural networks carry task-relevant information, but what are the origins of these FC modules? Surprisingly, we found that FC modules are present prior to and throughout training these networks. In fact, we found that training these networks - as well as feed-forward neural networks - can be interpreted as a search for the correct set of modules. Indeed, the rate of exploration is related to network performance across training (\(r=-0.87\); see **Supplemental Section** 2 and Fig. S6).

The presence of these FC modules before the networks were trained presented us with an insight about their origins. Before training, the weights of the neural network were randomly assigned. This meant that FC modules could be produced by randomly assigned weights. Upon further investigation, we found that these FC modules were being produced by the input layer of the RNN (Fig. S7b,c). That is, before the activity was transformed by the dynamics of the recurrent layer, this activity was already modular. Furthermore, we could modulate the number of FC modules created by the input layer by changing the number of input neurons (Fig. S7d, \(r=0.84\)).

In the parlance of linear algebra, the input layer is simply a projection from an \(N\)-dimensional space to an \(M\)-dimensional space, where \(N\) is the number of input neurons, and \(M\) is the number of neurons in the next layer. For this reason, these FC modules could be recovering the low-dimensional structure of the \(N\) input dimensions. Although appealing, this suggestion is incomplete given that these modules are still found when \(N>M\) as in the case of our supplemental analyses on feed-forward neural networks (Fig. S2).

Instead, we hypothesized that FC modules are created through weight-based "competition" for output activity. That is, some sets of outputs will receive more weight from input neuron 1, and another set will receive more weight from input neuron 2, and so on. Each of these sets will correspond to different modules (and each modules activity will be related to the "winning" input; see Fig. 1e). When the number of input neurons \(N\) get large enough, some of the input neurons do not "win" any of the output neurons. In a supplemental analysis, we demonstrate this, and show how the outcome of this competition, and therefore the structure of these modules, is also dependent upon the statistics of the input (see **Supplemental Section** 4 and Fig. S14).

In other words, this "competition" involves asymmetries in the weights of the connections from input neurons. If one neuron has more weight than the other neuron, it "wins". These asymmetries come in two basic types: asymmetries in sign (+/-), and asymmetries in weight. In a supplemental analysis, we develop generative models of these asymmetries in sign and weight to show how they modulate the modularity of the resulting FC (see **Supplemental Section** 3 and Fig S4). In this analysis we also show that in a simple model system with only two input neurons, the modules produced by the input layer can be nearly perfectly circumscribed by six connectivity-based rules (Fig. S4n-p). Importantly, we find that these connectivity-based rules can be approximated for larger systems (where \(N>2\)) by taking the cosine similarity of the connection weights (Fig. S8).

More specifically, for an \(N\times M\) weight matrix, we can compute the cosine similarity between the input weights
of all \(i\) by \(j\) outputs to produce a connection similarity matrix of size \(M\times M\) (see Fig. 2a for a schematic). By clustering this matrix of connection similarity we find modules that can approximate the connectivity-based rules (Fig. 2b and Fig. S8d). Additionally, this similarity matrix is nearly a perfect estimate of the FC of activity produced by sending Gaussian noise through the input projections (Fig. 2c-d; \(r=0.99\)).

That said, in RNNs this projection of activity from the input layer acts as a perturbation to the current state of the recurrent layer. This means that these functional connectivity modules will be altered by the _dynamics_ of the recurrent layer. As an example, we have seen how the input projection modules in an RNN trained on the perceptual decision-making task carry information about the current input stimulus value (Fig. 1e,g). However, as this input activity is transformed by the dynamics of the recurrent layer, the new FC modules begin to accumulate information about the current difference between the cumulative means of each stimulus (Fig. 1f,h). It is therefore an open question whether or not the FC modules defined by input projections (_before dynamics_) offer a reasonable partition of the modules found in the FC of the recurrent layer (_after dynamics_).

Here, we directly tested if the modular partition of the connection similarity matrix can be used to partition the FC of the recurrent layer of an RNN trained on a perceptual decision-making task. We can estimate the quality of a partition of FC into modules using the measure of modularity \(Q\). By imposing the modular partition of the connection similarity matrix on the FC of the recurrent layer during task trials, we induced a modularity of \(Q=0.33\) (Fig. 3c). We tested if this value was greater than chance by generating a null distribution of induced \(Q\) values. Briefly, we randomly permuted the order of this partition 1000 times, while preserving the number and size of each community and calculated the \(Q\) induced by imposing the permuted partition on the network for every permutation. We found that the real induced \(Q\) value was significantly greater than expected by our null model (Fig. 3d; \(p<10^{-15}\)). Indeed, we found a strong positive relationship between the input projection's connection similarity and the FC of the recurrent layer (Fig. 3e; \(r=0.51,p<10^{-15}\)) suggesting that the input projections greatly influence both the overall FC and the modules that it forms even in the presence of recurrent dynamics.

Taken together, these results suggest that input projections greatly influence the development of FC modules in recurrent neural networks, but it remains unclear whether or not the same can be said for the development of FC modules in _real_ brains. Here, we explore this using structural and functional connectivity data from the brains of both mice and humans. Given our interest in the development of FC modules in the cortex, we hypothesized that the thalamus could be modeled as the input projections to the cortex. In addition to its varied roles in information propagation and modulation [63, 64, 65], the thalamus is also the primary hub for relaying sensory information from the periphery to the cortex. As such, the thalamus might be a good analog of the input projections found in our model. Sensory information from the eyes, ears, and body often passes through the thalamus on its path towards the cortex [65, 66, 67].

First, we tested this in mice (Fig. 3f). We obtained weights for the structural connections from the thalamus to the cortex using publicly available tract tracing data from the Allen Brain Institute [68]. We took the cosine similarity of all in-weights from the thalamus to the left hemisphere of the cortex to create an \(M\times M\) matrix, where \(M\) is the number of cortical regions in the left hemisphere (\(M=2166\)). We then used modularity maximization to partition this similarity matrix into modules (Fig. 3g). We also collected functional magnetic resonance imaging (fMRI) data from lightly anaesthetized mice [69]. When we applied the modular parti

Figure 2: **Before dynamics, functional connectivity modules are created from the input projections.** (_a_) Schematic showing how to create the cosine similarity matrix representing input connection similarity. The input weights of each \(i\times j\) output neuron are related to one another using cosine similarity. (_b_) The process described in the previous panel results in a cosine similarity matrix representing the similarity of the input weights between all pairs of output neurons. We used modularity maximization with the Louvain algorithm to find modules in this matrix. (_c_) Using the same modular partition from the cosine similarity matrix, we reordered the functional connectivity that emerges from the input projections when inputting Gaussian noise. The modules are a near perfect match. (_d_) The cosine similarity and the functional connectivity that emerges from the input layer are nearly identical.
tion of the thalamocortical connections similarity matrix onto FC of the mouses cortex it induced a modularity of \(Q=0.17\) (Fig. 3h) We then generated a null distribution of induced \(Q\) values that we should expect given spatial autocorrelation in fMRI cortical data [70]. This involved randomly reordering nodes in a way that approximately preserved the variogram of the original data [70]. In this way, we produced 100 random partitions. We found that we could induce more modularity (\(Q\)) in mouse cortical FC using the connection similarity matrix partition than expected by chance (Fig. 3i; \(p<10^{-15}\)). We also found a strong positive relationship between thalamocortical similarity and FC values in the cortex (Fig. 3j; \(r=0.5,p<10^{-15}\)). We performed a similar analysis with other subcortical areas (subplate, pallidum, hypothalamus, pons, medulla, midbrain, hpf, striatum, olf, cerebellum) and we found that the relationship with the thalamus was the greatest (Fig. S15; two-sample \(t\)-test \(p=0.0057\)). This suggests that thalamocortical connections contribute to the FC modules found in the mouse cortex.

Next, we tested this in humans using dense structural and functional connectivity data from the human connectome project [71; 72; 73]. Using the same method that we used with the mice data, we again found that when we applied the modular partition of the thalamocortical pro

Figure 3: **After dynamics, input projections contribute to functional connectivity modules**, (_a_) Schematic of input projections onto the recurrent layer. (_b_) Input projection similarity matrix reordered using modularity maximization. (_c_) Functional connectivity of the recurrent layer reordered by the partition of the input projection similarity matrix. (_d_) Boxplot showing the null distribution of induced modularity (Q) values that we should expect by chance. This was produced by randomly permuting the partition labels from \(b\) and applying them to \(c\). The real induced modularity (Q) value is in blue. (_e_) Plot showing the relationship between the cosine similarity of input projections and the functional connectivity of the recurrent layer. Dot color and size indicates the number of points that fell in this bin. (_f_/_k_) Schematic showing the thalamus and the cortex in mice/humans. (_g_/_l_) Thalamocortical input projection similarity matrix reordered using modularity maximization (for mice/humans; human matrix down-sampled for plotting). (_h_/_n_) Functional connectivity of the cortex reordered by the partition of the thalamocortical input projection similarity matrix (for mice/humans; human matrix down-sampled for plotting). (_i_/_n_) Boxplot showing the null distribution of induced modularity (Q) values that we should expect by chance. This was produced by randomly permuting the partition labels from _g_/_l_ in a way that maintains the spatial autocorrelation in fMRI data and applying them to _h_/_m_. The real induced modularity (Q) value is in blue. (_j_/_o_) Plot showing the relationship between the cosine similarity of input projections and the functional connectivity of the recurrent layer. Dot color and size indicates the number of points that fell in this bin (for mice/humans).
jections onto FC of the human cortex (left hemisphere; number of nodes \(M=29696\); Fig. 3k-m) we could induce more modularity (\(Q\)) than expected by chance (Fig. 3i). We also found a strong positive relationship between thalamocortical similarity and FC values in the cortex in humans (Fig. 3o; \(r=0.26,p<10^{-15}\)).

In a supplemental analysis, we found that the thalamocortical connection similarity values were significantly concentrated in 7 brain systems defined based on the correlational structure of resting-state fMRI across 1000 subjects [15; 74] (Fig. S12a-b; two-sample \(t\)-test all \(p<10^{-15}\)), but the brain systems for which the concentrations were highest were primary sensory systems (visual, and somatomotor; Fig. S12c, two-sample \(t\)-test \(p<10^{-15}\)).

Taken together these results suggest that thalamocortical connections contribute to the structure of FC modules found in the cortex of both mice and humans. But what - if anything - does this say about the developmental origin of FC modules in the cortex?

We hypothesize that initial weights for the connections from the thalamus to the cortex can bias the development of FC modules. In this way, the thalamus might determine the broad placement of FC modules in the cortex. Importantly though, mice and humans learn, and this learning involves synaptic changes. Therefore, this relationship between initial weights and modules should be robust to weight changes that occur with learning. In a supplemental analysis we find evidence of this robustness of the development of FC modules to weight changes that occur with learning in recurrent neural networks (see **Supplemental Section**5 and Fig. S9).

### Dynamics

In this section we report our findings investigating the intersection between functional connectivity (FC) modules and the dynamics of the recurrent layer in RNNs. In the previous sections, we have shown that FC modules can emerge from asymmetries in the weights of input projections (Fig. S4) and that these asymmetries also produce a relationship between the current value of stimulus input information and the activity of the FC modules that emerge from the input projection (see Fig. 1e,g & Fig. S5). That said, RNNs have _dynamics_ and the tasks that they are well-suited to perform involve the memory and/or transformation of information across time. Therefore, the information regarding the current value of the input stimulus must either be maintained or transformed by the dynamics depending on the nature of the task. Here, we investigate this process directly by comparing the structure of FC modules with the effects of lesioning neurons (or their connections) in the recurrent layer.

In the first section on dynamics, we describe the result of lesioning connections within _input projection modules_ to show that the dynamics within these modules specialize in accumulating information within these modules. In the next section, we describe the result of lesioning neurons within _functional connectivity modules in the recurrent layer_ to show that lesions to neurons within the same FC module are likely to have similar effects on the phase portrait of the system.

Lesioning recurrent connections within input projections modules has circumscribed effects on behavior

In this section, we show that the input projection modules, while concentrating stimulus information into certain nodes in the recurrent layer also circumscribe the specific effects caused by lesions to the weights in the recurrent layer.

In order to test the specificity of lesions to tracking task relevant information from each stimulus, we set up a perturbation paradigm where we artificially sent a large amount of input into only one of the input neurons (corresponding to one of the stimuli (see Fig. 4a) for a schematic). Given such a perturbation, the RNN will infer that the artificial stimulus comes from a distribution with a much larger mean. We tested the hypothesis that lesions to weights in the recurrent layer that were associated with each input projection module would have specific effects on this inference. One input projection module is defined by having more weight from one stimulus, and the other module is defined as having more weight from the other stimulus. Our hypothesis was that lesions to each module would have specific effects on tracking the stimulus that it was associated with.

Previous work has shown that RNNs trained to perform this perceptual decision-making task set up a dynamical object referred to as a line attractor [55; 76]. A line attractor is a line formed by many fixed-point attractors such that the systems long term behavior will end up somewhere on this line in state space (Fig. 4b). When such a system is resting on this line attractor and is perturbed, the state of the system will return to another location along the line attractor. These line attractors have previously been shown to be used often by RNNs when they are trained to track continuous variables. In the case of this perceptual decision-making task, stimulus from each input neuron will perturb the state of the RNN in one of two general directions along the line attractor corresponding to accumulating information about the current difference between the stimulus means.

Along the center of this line attractor is a decision-boundary (see Fig. 4c for a schematic). When the state of the system is on one side of this boundary it will make one decision, and it will make a different decision on the other side of this boundary. In this way, movement away from the decision boundary corresponds to increasing evidence for one decision and against another decision. Indeed, when we perturbed the RNN at the input for stimulus 1, the state of the system traveled away from the decision boundary and fixed itself on a leftward portion
of the line attractor. In contrast, when we perturbed the RNN at the input for stimulus 2, the state of the system traveled away from the decision boundary in the other direction (Fig. 4c).

For our lesioning analysis, we used distance from the decision boundary as a proxy for the network's ability to accumulate relevant evidence for each stimulus. We operationalized the decision boundary based on the difference between the activity of the output neurons that made the decision for the RNN. We then trained 100 RNNs and for each RNN ran it through four lesioning conditions: _1)_ lesion module 1 and perturb stimulus 1, _2)_ lesion module 1 and perturb stimulus 2, _3)_ lesion module 2 and perturb stimulus 1, and _4)_ lesion module 2 and perturb stimulus 2. For each condition, we gradually lesioned a larger number of positively weighted connections (ordered from largest to smallest). We found that lesions to module 1 moved the system closer to the decision boundary when stimulus 1 was perturbed, but had a limited effect on perturbations to stimulus 2 (sometimes causing the system to move further from the decision boundary; Fig. 4f). In contrast, lesions to module 2 moved the system closer to the decision boundary when stimulus 2 was perturbed, but had a limited effect on perturbations to stimulus 1 (Fig. 4g). We found this general feature of the specificity of lesioning to modules across all 100 RNN models (Fig. 4h and Fig. S10c; two-sample _t_-test between stimulus 1 distances and stimulus 2 distances after 200 lesions to module 1: \(p<10^{-15}\), and to module 2: \(p=5.25\times 10^{-13}\)). We also tested this effect on modules defined using the

Figure 4: **Lesioning recurrent connections within module projections has circumscribed effects on behavior**, (_a_) Schematic showing how we perturb different input neurons during our lesioning trials. (_b_) Plot of the activity in the recurrent layer of a trained RNN projected into the first two principal components. In red, we plot the fixed/slow points approximated using a gradient-descent based method [75]. The activity trajectory is colored according to the correct decision about which stimulus came from the distribution with the greater mean. Note that these colors sometimes overlap given that the cumulative mean can be artificially higher for the incorrect stimulus early in the trial due to sampling variability. (_c_) Same as b, but instead of plotting the fixed points, we plot the decision boundary for the network (this is a visual estimate; see main text for how the boundary was calculated). (_d_) Same as b, but instead of plotting the fixed points, we plot the trajectories of two perturbation trials. In the green trial we perturbed stimulus 1. In the purple trial we perturbed stimulus 2. The state of the system starts at the black star and perturbations result in activity that stably rests at the colored stars. (_e_) Schematic showing how we lesioned the weights of the recurrent layer of the network based on the modules defined by the input projections. (_f_) Lesions to module 1 of the input projection cause the end points of stimulus 1 perturbations to move closer to the decision boundary, whereas the end points for the stimulus 2 perturbation in this example move further away from the decision boundary. (_g_) Showing a similar but opposite effect as j when lesioning module 2. (_h_) These plots show results from our four lesioning conditions (as described in main text) when applied based on modules defined by the input node giving the most weight to each output node in the input projection. When increasingly lesioning input projection module 1, stimulus 1 perturbations move closer to the decision boundary, but stimulus 2 perturbations do not move closer to the decision boundary. The opposite is shown for increasingly lesioning input projection module 2. These lines represent the average distance from the decision boundary across 100 trained RNNs.
sign of the input projections and we found the same effect (Fig. S10a,b; two-sample _t_-test between stimulus 1 distances and stimulus 2 distances after 200 lesions to module 1: \(p<10^{-15}\), and to module 2: \(p=2.01\times 10^{-14}\)).

Taken together, these results suggest that the FC modules defined by the input projections circumscribe functionally relevant effects of lesions to the weights of the recurrent layer. This further suggests that the symmetry breaking represented by these weight and sign asymmetries in the input layer are relevant for the development

Figure 5: **Similarity of lesioning effects on flow is related to functional connectivity modules**, (_a_) Quiver plots showing the flow of the activity in the recurrent layer of an RNN trained on the go vs. no-go task (projected onto the first two principal components). recurrent activity states start in a grid of points arrayed within the total state space that is explored during task trials. Arrows show where the system will end up after a single time step (from the current time step). The arrows eventually settle onto a limit cycle. Blue and red trajectories in the first panel indicate trajectories for different task conditions. (_b_) Plot of the effect of lesions on the flow shown in panel \(a\). A scatterplot of colored points is superimposed on a quiver plot describing the flow. The colors indicate the Euclidean distance between the non-lesioned flow after one time step and the flow after lesioning a neuron in the recurrent layer. Although plot is shown in two principal component dimensions, the distances are calculated in the original dimensions of recurrent layer activity. (_c_) Functional connectivity of the recurrent layer (in an RNN which has _not_ been lesioned). Modularity maximization was used to find modules in this matrix. This partition is also used to reorder the matrix in the next panel representing lesion similarity. (_d_) We individually lesioned the weights from and to every _i_-th neuron in the recurrent layer. This produced a grid of flow distance values for every neuron describing the effects of lesions on the dynamic flow. We then flattened this grid and compared the flow distances between all \(i\times j\) neurons producing a matrix of similarity values telling us how similar the effects of lesions were between all neurons in the recurrent layer. This matrix was reordered by the partition of functional connectivity into modules found in the previous panel. (_e_) Boxplot comparing the induced modularity (Q) of the lesion similarity matrix when using the modular partition of functional connectivity to a null model that randomly permuted this partition 1000 times. (_f_) Plot showing the relationship between lesion similarity and the functional connectivity of the recurrent layer. Dot color and size indicates the number of points that fell in this bin.
of task-based dynamics.

### Functional connectivity modules circumscribe sets of neurons with similar contributions to dynamics

In the previous section, we showed that the FC modules defined by the input projections delimit dynamics that are specific to the transformation of the input information concentrated in each input projection module. However, these modules are defined before the recurrent layer transforms this input activity via dynamics. An important open question is how the FC modules defined _after dynamics_ relate to the dynamics produced by the weights in the recurrent layer. Here we use a novel lesioning analysis to show that FC modules defined _after dynamics_ circumscribe sets of neurons with similar contributions to the dynamics.

We begin our analysis by visualizing the phase portrait of an RNN trained on the go _vs._ no-go task (Fig. 5a). A phase portrait can be used to visualize the dynamics of a system by initializing the state of that system in many different locations in state space(typically in a grid) and plotting the direction that each of these points moves in the phase space after some period of time. Here, we created a grid of initial points in the two-dimensional state space defined by the first two principal components of the activity in the recurrent layer. We then propagated this activity by stepping the RNN forward for a single time step and plotted the direction and magnitude of the resulting movement through state space using a quiver plot (Fig. 5a). After a sufficient period of time, all initial points fell onto a limit cycle (Fig. 5b,c).

In order to quantify the effects of lesions to the weights of neurons in the recurrent layer, we analyzed how the dynamics shown in these phase portraits changed following each lesion. Briefly, after lesioning a single neuron, we created a new phase portrait for the lesioned RNN and took the Euclidean distance between the lesioned and non-lesioned location in state space after one time step. Note that this distance was calculated in the original dimensions of the recurrent layer, where \(N=100\)(we refer to this as _flow distance_, see Fig. 5d-f). After lesioning the input and output connection weights for every neuron in the recurrent layer separately, we had an array of flow distance values for lesions to each neuron. We then measured the similarity of the dynamical effects of each neuron on the phase portrait by taking the Pearson correlation between every \(i\times j\) pair of flow distance arrays, resulting in an \(M\times M\) matrix of similarity values, where \(M\) is the number of neurons in the recurrent layer.

We found that this matrix of lesion similarity values describing the similarity of the dynamical contributions of different neurons to the phase portrait was highly related to the structure of the FC matrix _after dynamics_. Not only was there a positive linear relationship between lesion similarity and FC values (Fig. 5f, \(r=0.59,p<10^{-15}\)), but we also found that a modular partition of the FC matrix could also be used to identify modules in the lesion similarity matrix (Fig 5d,e; \(p<10^{-15}\)).

We replicated these results in an RNN trained on the perceptual decision-making task (Fig. S11) where again we found a positive relationship between FC and lesion similarity (Fig. S11c,d,f; \(r=0.26,p<10^{-15}\)), and a significant relationship between the modules in the FC matrix and the modules in the lesion similarity matrix (Fig. S11d,e; \(p<10^{-15}\)).

These results suggest that FC modules _after dynamics_ circumscribe sets of neurons with similar contributions to the dynamics of the system. Taken together with the results of lesioning input projection FC modules (_before dynamics_) these results suggest that FC modules can be used to identify meaningful task relevant dynamics that not only hold onto task-relevant information, but also circumscribe dynamics responsible for transforming input information into its task-relevant counterpart.

## Discussion

Through our research, we offer preliminary answers to three fundamental inquiries about functional connectivity modules, addressing their 1) function, 2) origin, and 3) dynamics. In the ensuing section, we explore the broader implications of these results for understanding functional connectivity modules in brains, suggesting prospective research paths, while acknowledging the limitations of our study.

### The origins of functional connectivity modules

In this paper we have shown that functional connectivity modules in recurrent neural networks are created by asymmetries in the weights of projections from the input layer. We found that these asymmetries concentrate input information into different areas of the recurrent layer in what we refer to as _"input projection modules"_. Although this information is then transformed via the dynamics of the recurrent layer, we found that we still see a correspondence between the modules defined by the weights of the input layer and the functional connectivity modules of the recurrent layer.

We then replicated these results in empirical structural and functional data from both mice and humans, finding agreement between the asymmetries of thalamocortical projections and cortical FC modules. Here, we argue for the plausibility of these results and their potential implications for the development of functional connectivity modules in biological brains.

If we model the cortex as the recurrent layer of this RNN, then a good analog for the input projections would be the thalamus. In tandem with its varied and complex role in information processing and modulation [64; 65; 77; 78], the thalamus is also a critical hub for
the communication of sensory information from the periphery to the cortex [65; 66]. Visual information from the retina passes through the lateral geniculate nucleus of the thalamus before being passed on to the visual cortex [66; 67]. Auditory information from the ears first passes through the medial geniculate body before being passed to the auditory cortex [79]. The ventral posterolateral nucleus receives information about pain, crude touch and temperature from the spine before passing it on to the somatosensory cortex [67; 80]. The thalamus is crucial for relaying information from the senses to the cortex for further processing.

Our results suggest that this important role of the thalamus as a relay also gives it a special place in the development of functional connectivity modules in the cortex. Given the clear functional roles these modules play in information processing and dynamics in these RNNs, and empirical results suggesting that functional connectivity modules in the cortex are highly related to the functional specialization of different regions (visual, somatosensory, etc.) [21; 25; 26], this might also have bearings on the development of functional specialization.

By asymmetrically delivering sensory information to different areas of the cortex, the thalamus could perform an important role in _symmetry-breaking_ during the development of cortical functional specialization. Symmetry-breaking is a concept from dynamical systems theory wherein wells or basins form in a previously flat attractor landscape [81]. Where before every state was equally likely, symmetry-breaking biases the evolution (or development [82]) of the system into a certain area of state space.

In the context of the development of functional specialization in the cortex, the symmetry-breaking performed by the thalamus could bias the development of areas like the visual, auditory, and somatosensory cortices by coding for asymmetries in the weights of thalamocortical connections. To developmental neurobiologists this is not an unfamiliar idea. Indeed, research from as early as the 1980s has suggested that redirecting thalamic projections can redirect functional specialization to new areas of the cortex [83; 84]. Indeed, more recent research involving transcriptomic identification of cell types in mice found that cortical modules identified by cell type were not only highly similar to those identified by connectivity-based methods, but thalamocortical connections refined cell type composition within these cortical modules [85] Our work enriches and supplies further evidence for this general hypothesis regarding the importance of thalamic connections for the development of functional specialization in the cortex.

Not only do we find that the asymmetries in thalamocortical projections predict modules in the functional connectivity of the cortex, but we also find evidence that could suggest that coding for the weights of these connections genetically would result in similar modules following development and learning. More specifically, we demonstrate that when the weights of input projections in RNNs are initialized with the same template prior to training, the resulting trained RNNs exhibit similar modular partitions. This suggests a possible developmental pathway - from genes to phenotype - for functional connectivity modules. That is, if some simple genetically coded mechanism, perhaps a chemical gradient, could bias asymmetries in the weights of thalamocortical projections, it might also have an outsized effect on the formation and organization of functional specialization in the cortex.

Indeed, recent empirical results involving the development of thalamocortical projections in premature human infants support this general idea. Thalamocortical connections develop from a transient brain area just below the cortex known as the cortical subplate [86; 87]. During a period known as the _"waiting period"_(between 20-38 weeks), thalamic projections embed themselves in the cortical subplate and _wait_ for their cortical postsynaptic neurons or target neurons [86; 87].

Because of this periods proximity to birth at 40 weeks, premature infants are particularly vulnerable to interruptions to this period [86; 88]. Many studies have found that premature birth can result in cognitive deficits (e.g.[89; 90]), but some work also suggests that these cognitive deficits are directly due to malformed thalamocortical projections [88; 91]. Perhaps one reason that malformation to these thalamocortical projections is detrimental to cognition is that it interrupts the development of functional specialization in the cortex. Indeed, more recent work from our lab has actually shown that the development of functional connectivity modules in premature human infants is a predictor of better cognitive outcomes one year later (unpublished).

Together with recent work on the community structure of thalamocortical projections [92; 93], the recapitulation of cortical gradients by thalamic anatomy [63], recent work showing a genetic and connectomic axis of variation in the thalamus relates to an anterior-posterior pattern in cortical activity [94], and general theory about the relationship between the thalamus and cortical functional connectivity [95], our work supports an emerging theory about the importance of the thalamus for the development of functional connectivity modules, as well as the development of functional specialization in the cortex.

Although we primarily focus our discussion here on the role of the thalamus in providing input projections to the cortex, it is important to note that our results on the origins of functional connectivity modules are not limited to the thalamus/cortex. Instead, our results could apply broadly to the connections impinging on any population of neurons and suggest that you could predict the functional connectivity modules that emerge in that population by measuring the asymmetries in the weights of those connections. This general result is consistent with recent work showing that low-rank recurrent neural networks produce low-dimensional dynamics [96]. Indeed, when our measure of similarity for grouping asymme
tries in input projections (cosine similarity) is high, it suggests that the columns of the weight matrix are not independent and therefore exhibit lower rank.

### The function and dynamics of functional connectivity modules

Functional connectivity modules are nearly universally present in neural activity. From calcium imaging data of single cell activity in the larval zebrafish [7], to fMRI imaging of regional activity in the cortex of human beings [15; 20; 21]. Although previous research has attempted to disambiguate the functional roles played by functional connectivity modules, their specific roles remain unclear.

For example, meta-analytic work on functional connectivity modules during the resting-state suggests that they recapitulate the co-activation of brain regions during different tasks such that categories of tasks map onto different modules [25; 26]. Indeed, common parcellations of cortical activity into large scale systems using functional connectivity analysis suggest a mapping exists between FC modules and cytoarchitectural features of primary sensory cortices as well as functional systems defined by neuropsychological and animal ablation studies [21; 97]. Additionally, recent work has used task and resting-state fMRI in tandem to infer that large-scale functional connectivity modules can be sub-divided into smaller functional connectivity modules with task-specific functional domains [98; 99]. That said, it is unknown whether or not a complete picture of the function of these modules would require cell-level imaging (beyond the current spatial resolution capability of fMRI) and perturbation (beyond that which is ethically permissible in humans).

At the microscale of neurons, some recent research involving imaging and perturbation in mice provides an even clearer picture of the possible function of functional connectivity modules. For example, after training mice on a go _vs._ no-go task researchers found cellular assemblies related to the go signal and no-go signal in layer 2/3 of the mouse primary visual cortex. This suggested that these assemblies might be keeping track of the presence of either signal. In fact, they found that optogenetic stimulation of the go signal assembly was sufficient to produce go behavior [10].

Taken together, this research is suggestive of various roles that functional connectivity modules _might_ play in functional specialization, but it leaves open the question of the _specific_ function of these modules. Recurrent neural networks provide a good platform for exploring such functional questions. Not only do recurrent neural networks maintain functional dynamics that emerge from network interactions in a manner similar to complex adaptive systems like the brain, but with recurrent neural networks we have complete access to all of the information responsible for behavior. Additionally, with recurrent neural networks there are no ethical issues with perturbing/lesioning the system to further probe the nature of its behavior.

Here, we directly explored the _function_ of functional connectivity modules in recurrent neural networks. We found that the function of these modules is to hold onto, or accumulate task-relevant information. For example, by restricting the network from using information from these modules for its decision, we were able to show that the network exhibits specific behavioral deficits. Furthermore, these deficits could be predicted by its lack of access to the information found in these modules.

Similarly, previous research has used the dynamical systems lens to probe the low-dimensional dynamics of neural data [4; 6; 100]. Such research has shown that axes in these low-dimensional spaces can be shown to represent task-relevant variables [6; 55; 100; 101] and that dynamics in these spaces are responsible for transforming this data and using it for behavior [55; 101]. While functional connectivity modules imply that a system is exhibiting low-dimensional behavior, it remains unclear how these dynamics are related to functional connectivity modules and their function.

Here, we directly explored the _dynamics_ of functional connectivity modules in recurrent neural networks. We found that these functional connectivity modules delimit neurons with similar dynamic effects. For example, we found that lesions to neurons within a module produced similar changes to the system's phase space. Additionally, we found that lesions within input projection modules result in specific disruptions to the dynamics responsible for transforming input information from different stimuli.

Taken together with the previous results on the origins of functional connectivity modules, our results suggest an intriguing story about their _function_ and _dynamics_. When the RNN begins it's task, input information is projected preferentially into different neurons in the recurrent layer. Where this information is preferentially sent can be determined by the input projections modules. Now, this information must be transformed or maintained by the dynamics of the system. In the case of the perceptual decision-making task, these initial modules hold information about the current values of each stimulus. The dynamics of the system then transform this information into an estimate of the difference between the cumulative mean's of the two input stimuli such that _after dynamics_ functional connectivity modules hold onto this task-relevant variant of the input information.

So, what might these results suggest about functional connectivity modules when we see them in brains? At the scale of the entire cortex, our results suggest that the thalamus concentrates input information onto different areas of the cortex, and then the cortex transforms this information via dynamics. The functional connectivity modules of the cortex maintain a similar profile to the input projection modules defined by the thalamocortical projections, but they begin to reflect the dynamical signature of each node, and the information that they carry
is likely some task-relevant transformation of sensory information.

At the scale of neurons in cellular assemblies, these results suggest something similar. Input projections to a given population of neurons reflects how input information is concentrated in that population. Functional connectivity modules in this population could also delimit neurons with similar dynamical contributions to the phase space of the system. These functional connectivity modules also likely hold onto task-relevant information.

One important open question is how these scales relate to one another. While the general dynamical and information-processing perspective described here might apply well to both micro- and macro-scales, the relationship between these scales is likely a hierarchical one. That is, functional connectivity modules at the cortical scale defined by regional BOLD activity are likely composed of the functional connectivity modules at the neural scale. Indeed, modules at the macro scale have been consistently shown to have a hierarchical structure [24; 102; 23]. Although a true bridge between these spatial scales has not been directly explored, this research suggests that large modules are composed of smaller modules, perhaps all the way down to the level of cell assemblies. Future research should explore how such a hierarchical structure might emerge in recurrent neural networks and how the dynamics and information-processing across these scales relate and interact.

In closing, it's important to note that, although it is unlikely, functional connectivity modules could _just_ be the signature of dynamics in the system and play no specific functional role. In this way, modules in functional connectivity networks would be epiphenomenal, present as a result of the systems dynamics, but ultimately not used by the system. In our work, we diminished this possibility by allowing the system dynamics to persist while not allowing the information from these modules to be used for decision-making. Importantly, this destroyed the systems ability to make the correct decisions.

### Low-dimensional manifolds and functional connectivity modules

Our results suggest that functional connectivity modules hold onto, or accumulate, task-relevant information. Importantly, this same thing can be said for low-dimensional manifolds in neural activity [100; 101; 6; 55]. Here, we take the opportunity to clarify that we see these as complementary perspectives on a similar underlying phenomenon.

Both functional connectivity modules and low-dimensional manifolds share the feature of being about the statistical dependence between sets of neural units. Indeed, low-dimensional manifolds are often based on the eigendecomposition of a covariance matrix, which itself is a non-standardized version of a functional connectivity matrix. Although the mathematical relationship between the two is more difficult to work out given that functional connectivity modules are often found through optimization, the relationship is nonetheless non-trivial.

Many of our results extend the literature on the dynamic behavior of recurrent neural networks by showing that the familiar properties of low-dimensional manifolds can also be seen at the level of functional connectivity modules. Additionally, we also show a link between such dynamics and network properties of the weights in these systems. As such, our work offers a bridge between various disciplines that model the brain as a network and those that model the brain as a dynamical system. It is our belief that these two important frameworks can be used together to tell a richer story about the function of recurrent neural networks, and hopefully also, brains.

### Limitations

Our study has a number of important limitations to consider. First, recurrent neural networks are abstract models of brains, and as such remove some biological details from consideration in the model. This biological detail could prove to be fundamental to understand the system in question (e.g.[103; 104; 105; 106]). Additionally, these systems might solve problems differently than brains. As such, it is essential that the insights provided by recurrent neural networks be tested on empirical data from brains. We were able to test a portion of our results on the origin of functional connectivity modules, but future research should attempt to replicate our results on the function and dynamics of functional connectivity modules in biological brains.

In addition, the tasks that these RNNs solved, although based in classic systems neuroscience tasks, are also highly simplified. Biological organisms evolve in a complex environment replete with uncertainty, many different ill-defined tasks, and serious stakes regarding survival and reproduction. Future work should explore the function, origin, and dynamics of functional connectivity modules in recurrent neural networks trained to solve more ethologically relevant tasks.

Additionally, one of the primary claims made here is that modules carry information relevant to a task. One unaddressed concern relates to whether these two properties - modular structure and task-relevance - tend to be correlated. That is, in the universe of all possible RNNs trained on the tasks used here or elsewhere, if a network exhibits modules, must the modules also carry task-relevant information? Are there networks that solve the task equally well whose correlation structure is modular (note that correlation structure, itself, may artificially inflate estimates of modularity)? If so, how prevalent are they and how likely is it that RNNs trained using backpropagation will discover such networks as solutions to systems neuroscience tasks? An important direction for future work is to investigate the interrelationship of modular correlation structure and the task-relevance of
the detected modules.

Finally, one might argue that we see functional connectivity modules in these recurrent neural networks because the number of neurons in the recurrent layer is too large. Perhaps if the network only needs to keep track of three task-relevant variables, then three neurons would be sufficient to track this information. While we acknowledge that this could be the case, it is important to clarify that there is a difference between neurons holding onto information, and neurons producing dynamics that can transform input information into task-relevant information. In fact, it may be the case that the number of neurons needed to easily find dynamics to solve a task, might be greater than the number of neurons whose state could theoretically hold onto task-relevant information. Indeed, although the perceptual decision-making task only requires keeping track of roughly three variables, in a supplemental analysis we found that recurrent neural networks with a recurrent size of three had a more difficult time finding the dynamics that led to high accuracy with this task than networks with larger recurrent sizes (Fig. S13).

Additionally, the same critique can be leveled at research on the low-dimensional manifolds found in recurrent neural networks [5; 101; 5], and yet we find both functional connectivity modules and low-dimensional manifolds in neural data from biological brains, suggesting again that the number of neurons required to easily find and produce task-related dynamics might be larger than the number of task-relevant variables being tracked.

That being said, after these dynamics transform this information, it can then be read off of the relevant functional connectivity modules by reader cells similar to the way in which our output layer uses single read-out neurons to make decisions based on the current state of the recurrent layer. Indeed, this could also be the reason why neuroscience will continue to find single cells encoding complex variables (e.g.[107]). These cells then would be reading out information from upstream functional connectivity modules [12], but nonetheless the original functional connectivity modules would hold onto the information being transformed because they are the statistical signature of the dynamics transforming this information across time.

## Materials and Methods

### Recurrent Neural Network

We used a specific type of recurrent neural network (RNN) referred to as continuous time recurrent neural networks (CTRNNs). CTRNNs are defined by the following equation [108; 109; 110; 47]:

\[\tau\frac{d\mathbf{r}}{dt}=-\mathbf{R}(t)+f(W_{rec}\mathbf{R}(t)+W_{in}\mathbf{ Input}(t)+\mathbf{b}).\]

Here, \(\tau\) is a time constant (set to 100), \(\mathbf{R}\) is the recurrent state, \(W_{rec}\) are the weights of the recurrent layer, \(W_{in}\) are the weights of the input layer (input projection), and \(\mathbf{b}\) is a bias. For all implementations of RNNs in this project we chose CTRNNs with a recurrent layer size of \(M=100\). Finally, a linear transformation of the recurrent state \(\mathbf{R}\) is used to create output activity:

\[\mathbf{O}(t)=g(W_{out}\mathbf{R}(t)+\mathbf{b}_{out}).\]

For more information on the design and implementation of this CTRNN see the github repository [https://github.com/gyyang/nn-brain/blob/master/RNN_tutorial.ipynb](https://github.com/gyyang/nn-brain/blob/master/RNN_tutorial.ipynb) [108; 47]. For more information on the dynamics and parameter space of CTRNNs generally, see [110; 50; 111].

### System neuroscience tasks for recurrent neural networks

We used a machine learning toolbox called neurogym: ([https://github.com/neurogym](https://github.com/neurogym)). The toolbox is built upon the popular machine learning framework _gym_ from OpenAI and provides a variety of systems neuroscience tasks that have been designed for easy presentation to recurrent neural networks (see [109]). We chose two tasks from this toolbox to focus on: 1) the perceptual decision-making task, and 2) the go vs no-go task.

The perceptual decision-making task implements a simplified version of a random-dot motion task wherein the subject is presented with randomly moving dots that have coherent motion in some direction [112; 55]. The task of the subject is to indicate the direction of coherent motion during some trial period. In the neurogym version of this task, the trial period is determined by a fixation input to the RNN. Two stimuli are presented to the RNN during this fixation period and the task is to indicate which of the two stimuli has the larger mean value. Each of these stimuli comes from a Gaussian distribution with different means (equal standard deviation). The means of this distribution represent coherent motion, and the random distribution of values around the mean represent the random motion.

The go _vs_ no-go task also has a trial structure that is determined by a fixation input. At the beginning of the fixation period, the RNN is presented with one of two signals: a go signal or a no-go signal. These two signals are represented by two different input neurons. Following presentation of this signal, there is a delay period wherein the RNN only receives the fixation signal. The delay period ends when the fixation signal is no longer present, and the RNN must now determine which of the two input signals it received before the delay period (the go signal, or the no-go signal).

RNNs were trained on both of these tasks in PyTorch, an open-source machine learning library for the implementation and training of machine learning models using
backpropagation [113]. Pytorch simplifies the training of such models by automatically tracking the gradient of all computations in the forward pass of a neural network, and storing them in a computational graph that can later be used to backpropagate (via the chain-rule).

### Modularity maximization

In this paper we analyzed network modules in the functional connectivity network of RNNs. In order to detect these modules, we used a method referred to as _modularity maximization_. The method of modularity maximization[114] is based on a straightforward principle: comparing the observed connectivity data with what is expected by chance. This entails comparing the observed connectivity matrix \(A\), with another matrix of the same dimensions \(P\). The elements \(P_{ij}\) of matrix \(P\) represent the expected weight of the connection between nodes \(i\) and \(j\) under a null model. From these two matrices we define the modularity matrix, \(B\), as:

\[B=A-P.\]

Each element \(B_{ij}\) in the modularity matrix represents whether the observed connection between nodes \(i\) and \(j\) is stronger (\(B_{ij}>0\)) or weaker (\(B_{ij}<0\)) than expected under the null model. Modularity maximization then uses the modularity matrix to assess the quality or "goodness" of a modular partition, which is a division of the network's nodes into non-overlapping modules or communities. The quality of a partition is quantified using the modularity function, \(Q\), calculated as:

\[Q=\sum_{ij}B_{ij}\delta(\sigma_{i},\sigma_{j}).\]

where \(\delta\) is the Kronecker delta function, and \(\sigma_{i}\) and \(\sigma_{j}\) are the community labels of nodes \(i\) and \(j\), respectively.

In addition to simply assessing the quality of a given partition, the variable \(Q\) can be optimized outright to identify high-quality partitions of a network's nodes into modules. We used the Louvain algorithm to optimize \(Q\). For a tutorial on using this method, see Esfahlani _et al._[115].

For our specific implementation, we used a variant of the modularity function (\(Q\)) that treats the contribution of positive and negative edges asymmetrically [116]. For an implementation of this method, see the **community_louvain** function in the brain connectivity toolbox: [https://sites.google.com/site/bctnet/](https://sites.google.com/site/bctnet/)

### Output lesions

In our **Results** section on the _Function_ of FC modules in recurrent neural networks, we use a lesioning method that we refer to as "output lesions" to test if the RNN is using information from its modules to perform the task. These output lesions are lesions to the weights of the output layer of the RNN. The output layer performs a linear transformation of activity in the recurrent layer such that the \(M\)-dimensional activity of the recurrent layer is translated into the activity of \(O\) output neurons. The output lesions were performed in the following way:

We implemented the modularity maximization method optimized using the Louvain algorithm to partition the functional connectivity of the recurrent activity into \(modules\). We then use these modules to lesion the connections from a given module to the output nodes:

\[\text{for selected module }m,\in modules,\,W_{out}[:,m]=0.\]

here, \(W_{out}\) be the weight matrix for connections from the recurrent layer to the output layer

### Generative models

In the **Supplemental Section** on the origin of functional connectivity modules in RNNs (**Section** 3), we developed three generative models of functional connectivity modules. Each of these generative models involved creating a new input weight matrix, and the parameters of the models allowed us to explore how modularity (\(Q\)) changed as we varied parameters in the model.

The first generative model was the sign-based model. This model creates an input weight matrix \(W\) of size \(N\times M\) out of two values \((+1,-1)\). That is:

\[W_{ij}=\begin{cases}-1&\text{if }i\in S\\ 1&\text{otherwise}\end{cases}\]

\[S=\{i_{1},i_{2},\dots,i_{k}\}.\]

\[k=\left\lceil N\times p\right\rceil.\]

where \(S\) is a random subset of row indices up to \(k\), and \(k\) is defined based on the percentage \(p\) of negative weights to be added to the weight matrix \(W\). All columns of this matrix \(M\) are given the same value.

The second generative model was the difference-based generative model. This model creates an input weight matrix \(W\) of size \(N\times 2\) by randomly assigning one columns value from each row to \(\alpha\) and assigning the other columns value to \(\alpha+\gamma\):

\[W_{i1},W_{i2}=\begin{cases}(\alpha_{i},\alpha_{i}+\gamma)&\text{with probability }0.5\\ (\alpha_{i}+\gamma,\alpha_{i})&\text{with probability }0.5\end{cases}\]
\[\alpha_{i}=\left|x_{i}\right|,\text{ where }x_{i}\sim\mathcal{N}(0,1).\]

\(\alpha\) is the absolute value of a random variable drawn from a Gaussian distribution with a mean of zero and a variance of 1.

Our third and final generative model included both sign and difference elements and is likewise refered to as the sign-difference generative model. The only difference between the previous model and this model is the definition of \(\alpha\).

\[\alpha_{i}\sim\mathcal{N}(0,1).\]

That is, whereas the previous model took the absolute value of the randomly drawn variable, this model does not, allowing approximately half of these values to be negative.

### Fixed point approximation

In our **Results** section on _Dynamics_, we illustrate that the perceptual decision-making task results in an approximate line attractor (Fig. 4b). In order to approximate the attractors that form this line, we used a simple gradient based method. For the implementation of the fixed-point approximation method that we used (also showing an approximate line attractor with the perceptual decision-making task) see: [https://github.com/gyyang/nn-brain/blob/master/RNNY2BDynamicalSystemAnalysis.ipynb](https://github.com/gyyang/nn-brain/blob/master/RNNY2BDynamicalSystemAnalysis.ipynb).

Briefly, this process involves optimizing the hidden/recurrent activity of the RNN such that the mean squared error (MSE) between that activity and the hidden/recurrent activity one step forward in time is minimized. Attractors occur where the derivative of a system is equal to zero. Intuitively then, if the difference between a current state and the next state is brought to zero, the system is in an attractor state. Alternatively, the loss could be very low but not zero, a case that is sometimes referred to as a slow-point [75]. This optimization is implemented by randomly initializing the recurrent state of the RNN in many different states, and then using back-propagation to minimize the MSE between current activity and the future activity.

### Input projection lesions

In our **Results** section on _Dynamics_, we explored the effects of lesions to the weights of the recurrent layer of RNNs trained on the perceptual decision-making tasks. We started by using sign-based partitons or difference-based partitions to define the input projection modules (Fig. S4d and i respectively). Briefly, sign-based partitions were determined by taking the sign of weights from one of the input neurons. Neurons receiving positive weights were in one module. Neurons receiving negative weights were in another. Difference-based partitions were determined by comparing the weights of the two input neurons onto each recurrent neuron. If \(inputweight1>inputweight2\) then the recurrent neuron was placed in module 1, and if \(inputweight1<inputweight2\) then the recurrent neuron was placed in module 2.

Then, we tested the effects of lesions to RNNs that were"perturbed" with a large continual input value of 4. More specifically, we lesioned positively weighted connections in the recurrent layer that were associated with the same input projection module:

\[W^{\prime}_{ij}=\begin{cases}0&\text{if }i,j\in M_{k}\text{ and }W_{ij}=\max_{a,b\in M_{k}}(W_{ab}),\\ W_{ij}&\text{otherwise}.\end{cases}\]

Here, \(W\) is the weight matrix of the recurrent layer, \(W^{\prime}\) is the modified weight matrix after the lesioning, \(M_{k}\) is a module in the network, and \(i\) and \(j\) are neurons in the recurrent layer. If neurons \(i\) and \(j\) belong to the same module \(M_{k}\) and the weight \(W_{ij}\) is the maximum weight among all the connections within the module, then the weight is set to zero. Otherwise, the weight remains the same. We iterated this equation such that \(W^{\prime}\) became \(W\) on future iterations. In this way, we lesioned an increasing number positively weighted connections within the same module and on each iteration we tested the distance of the perturbed activity from the decision boundary. Given that these RNNs were trained on cross entropy loss, the output node with the greatest activity corresponded to the RNNs "choice". For this reason, we defined the decision boundary as the difference between the output activity of the two output neurons that correspond to decision 1 and decision 2.

### Lesion similarity

In our **Results** section on _Dynamics_, we explored the effects of lesions to neurons in the recurrent layer of the RNN by quantifying how they changed the phase space of the system. First, we defined the phase space of the RNN by performing principal component analysis (PCA) on the activity of its recurrent layer across many task trials. Then, we defined a grid of \(25\times 25\) points (for a total of 625) in the 2-dimensional space defined by the first two principal components. These points are spread equally so that they span the entire activity space explored during task trials. Then, we use each of these points as an initial recurrent state for the RNN (projecting it back into the original activity space) and step it forward in time by one time step. This results in a new set of 625 points corresponding to the location of the recurrent state after one time step from each point on the original grid. We stored these points as the original points \(Orig\) (importantly, we
stored each of these points in the original \(M\)-dimensional activity space).

For the lesioning analysis, we would lesion a node and then use the same grid of initial recurrent states to initiate the system and run it for a single time step. This resulted in a new set of forward stepped points \(L\) for the lesioned model. We created \(M\) separate lesioned models where \(M\) is the number of neurons in the recurrent layer. For each model, we lesioned a different neuron, and then took the Euclidean distance between each point in \(L_{i}\) and the corresponding point in \(Orig\), resulting in a 625 length array of lesion distance values for every lesioned model.

Next, in order to quantify how similar the effects of different lesions were, we calculated the similarity of lesions by correlating the lesion distances from each lesioned model, resulting in an \(M\times M\) matrix of similarity values that we refer to as the _lesion similarity_ matrix.

### Mouse resting state fMRI data

All in vivo experiments were conducted in accordance with the Italian law (DL 2006/2014, EU 63/2010, Ministero della Sanita, Roma) and the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. Animal research protocols were reviewed and consented by the animal care committee of the Italian Institute of Technology and Italian Ministry of Health. The rsfMRI dataset used in this work consists of \(n=19\) scans in adult male C57BL/6J mice that are publicly available [117; 118]. Animal preparation, image data acquisition, and image data preprocessing for rsfMRI data have been described in greater detail elsewhere [119]. Briefly, mice were anesthetized with isoflurane (5% induction), intubated and artificially ventilated (2%, surgery). The left femoral artery was cannulated for continuous blood pressure monitoring and terminal arterial blood sampling. At the end of surgery, isoflurane was discontinued and substituted with halothane (0.75%). Functional data acquisition commenced 45 minutes after isoflurane cessation. Mean arterial blood pressure was recorded throughout imaging sessions. Arterial blood gasses (\(paCO_{2}\) and \(paO_{2}\)) were measured at the end of the functional time series to exclude non-physiological conditions. rsfMRI data were acquired on a 7.0-T scanner (Bruker BioSpin, Ettlingen) equipped with BGA-9 gradient set, using a 72-mm birdcage transmit coil, and a four-channel solenoid coil for signal reception. Single-shot BOLD echo planar imaging time series were acquired using an echo planar imaging sequence with the following parameters: repetition time/echo time, 1200/15 ms; flip angle, 30\({}^{\circ}\); matrix, \(100\times 100\); field of view, \(2\times 2cm^{2}\); 18 coronal slices; slice thickness, 0.50 mm; 1500 volumes; and a total rsfMRI acquisition time of 30 minutes. Timeseries were despicked, motion corrected, skull stripped and spatially registered to an in-house EPI-based mouse brain template. Denoising and motion correction strategies involved the regression of mean ventricular signal plus 6 motion parameters. The resulting time series were band-pass filtered (0.01-0.1 Hz band) and then spatially smoothed with a Gaussian kernel of 0.5 mm full width at half maximum. After preprocessing, mean regional time-series were extracted for 15314 regions of interest (ROIs) derived from a voxelwise version of mouse structural connectome [120; 68; 121].

### Mouse Anatomical Connectivity Data

The mouse anatomical connectivity data used in this work were derived from a voxel-scale model of the mouse connectome made available by the Allen Brain Insitute [120; 68] and recently made computationally tractable [121]. Briefly, the structural connectome was obtained from imaging enhanced green fluorescent protein (eGFP)-labeled axonal projections derived from 428 viral microinjection experiments, and registered to a common coordinate space [122]. Under the assumption that structural connectivity varies smoothly across major brain divisions, the connectivity at each voxel was modeled as a radial basis kernel-weighted average of the projection patterns of nearby injections [120]. Leveraging the smoothness induced by the interpolation, neighboring voxels were aggregated according to a Voronoi diagram based on Euclidean distance, resulting in a \(15314\times 15314\) whole brain weighted and directed connectivity matrix [121].

### Human connectomic data

Structural, diffusion, and functional human brain magnetic resonance imaging (MRI) data was sourced from the Human Connectome Project's (HCP) young adult cohort (S1200 release). This contained structural MRI (T1w), resting-state functional MRI (rs-fMRI), and diffusion weighted imaging (DWI) data from 1000 adult participants (53.7% female, mean age = 28.75, standard deviation of age = 3.7, age range = 22-37). A comprehensive report of imaging acquisition and preprocessing is available elsewhere [123]. In brief, imaging was acquired with a Siemens 3T Skyra scanner with a 32-channel head coil. The rs-fMRI data was collected with a gradient-echo echo-planar imaging (EPI) sequence (run duration = 14:33 min, TR = 720 ms, TE = 33.1 ms, flip angle = 52\({}^{\circ}\), 2-mm isotropic voxel resolution, multi-band factor = 8) with eyes open and instructions to fixate on a cross [124]. DWI data was acquired using a spin-echo planar imaging sequence (TR= 5520 ms, TE = 89.5 ms, flip angle = 78\({}^{\circ}\), 1.25 mm iso- tropic voxel resolution, b-values = 1000, 2000, 3000 s/mm2, 90 diffusion weighed volumes for each shell, 18 b = 0 volumes) [125]. HCP minimal processing pipeline was used to preprocess the functional and diffusion imaging data [123]. In particular, rs-fMRI underwent gradient distortion correction, motion correction, registration to template space, inten
sity normalization and ICA-FIX noise removal [123; 126]. The diffusion preprocessing pipeline consisted of b0 intensity normalization, EPI distortion correction, eddy-current-induced distortion correction, registration to native structural space, and skull stripping [123].

#### Diffusion tractography

A probabilistic streamline tractography pipeline tailored for the computation of high-resolution human connectomes [73] was implemented in MRtrix3 [127] which also adopted recent recommendations detailed elsewhere [73]. In particular, an unsupervised heuristic was used to estimate macroscopic tissue response functions for white matter (WM), gray matter (GM), and cerebrospinal fluid (CSF) [128]. Multi-shell, multi-tissue constrained spherical deconvolution was used to estimate fiber orientation distributions (FODs) [129]. This information was used to apply combined intensity normalization and bias field correction [130]. Liberal and conservative brain masks were respectively utilized in the last two steps to mitigate the detrimental effects of an imperfect mask on the respective procedures [73]. The normalized FODs were used to perform anatomically constrained tractography (ACT) [131]. A tissue-type segmentation was used to create a mask at the GM-WM boundary to enable brain-wide streamline seeding. Whole-brain tractography was conducted using 2nd-order integration over FODs (iFOD2) [132]. A total of five million streamlines (per participant) were generated that satisfied length (minimum length = 4mm) and ACT constraints.

#### Human structural connectome

The high-resolution structural connectivity network was constructed from the whole-brain tractograms [73]. Notably, streamline endpoints were mapped onto HCP's CIFTI space comprising two surface meshes of the cortex, in addition to volumetric delineations of several subcortical structures. The surface meshes included a subset of the fs-LR template mesh after removal of the medial wall with respectively 29696 and 29716 vertices for the left and right cortex. In addition, 31870 voxels from the MNI template space were included for subcortical brain regions, resulting in a total of 91282 high-resolution network nodes. Euclidean distance was used to assign each streamline to its closest node pair via nearest neighbor mapping. Streamlines with endpoints falling more than 2mm away from all nodes were discarded. Connectome spatial smoothing with a 6mm FWHM kernel was performed to account for accumulated errors in streamline endpoint location and to increase intersubject reliability of connectomes [72]. A group-level consensus connectome was constructed by aggregating individual connectivity matrices across all participants. To minimize the pipeline's computational complexity, group aggregation was applied before smoothing, considering that both aggregation and connectome spatial smoothing entailed linear operations, and the order of implementation does not impact the resulting connectomes.

#### Human thalamic projections

We next sought to estimate the thalamocortical structural connectivity projections from the smoothed high-resolution group-level human connectome model. To this end, the cifi structures for the left and right thalamus were combined to form a binary mask of the thalamus. This mask was subsequently used to aggregate group level thalamic projections to the left cortex.

#### Human functional network

The dense group-level functional connectivity network was provided by HCP. Specifically, this network was constructed from a pipeline combining high-resolution individual rs-fMRI data across all participants. First, the minimally preprocessed individual rs-fMRI data were aligned by a multimodal surface matching algorithm (MSMAll) [133]. Next, all timeseries were temporally demeaned followed by a variance normalization [134] and were passed to a MELODIC's Incremental Group-PCA [135]. The group-PCA outputs were renormalized, eigenvalue reweighted, and correlated to create a dense functional connectivity matrix (\(91282\times 91282\)). Finally, a subset of this dense connectome was extracted to denote left cortical functional connectivity (\(29696\times 29696\)).
