{'fraction_non_alphanumeric': 0.02838842461986125, 'fraction_numerical': 0.010911800713259168, 'mean_word_length': 4.918513073611257, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Recurrent networks improve neural response prediction and provide insights into underlying cortical circuits

Feedforward CNN models have proven themselves in recent years as state-of-the-art models for predicting single-neuron responses to natural images in early visual cortical neurons. In this paper, we extend these models with recurrent convolutional layers, reflecting the well-known massive recurrence in the cortex, and show robust increases in predictive performance over feedforward models across thousands of hyperparameter combinations in three datasets of macaque V1 and V2 single-neuron responses. We propose the recurrent circuit can be conceptualized as a form of ensemble computing, with each iteration generating more effective feedforward paths of various path lengths to allow a combination of solutions in the final approximation. The statistics of the paths in the ensemble provide insights to the differential performance increases among our recurrent models. We also assess whether the recurrent circuits learned for neural response prediction can be related to cortical circuits. We find that the hidden units in the recurrent circuits of the appropriate models, when trained on long-duration wide-field image presentations, exhibit similar temporal response dynamics and classical contextual modulations as observed in V1 neurons. This work provides insights to the computational rationale of recurrent circuits and suggests that neural response prediction could be useful for characterizing the recurrent neural circuits in the visual cortex.


Introduction
Feed-forward deep neural networks have been shown to be effective models for modeling and predicting neural responses of early visual areas in the brain [1][2][3][4][5][6]. However, abundant recurrent connections exist within and between visual areas and are crucial to generating neural responses [7,8]. While biologically-inspired recurrent circuits incorporated in convolution neural networks have been shown to be more efficient than purely feedforward models, accomplishing comparable performance in image classification tasks with less layers and less parameters [9][10][11], their underlying computational rationales remains to be understood. In this paper, we investigated the use of recurrent circuits in neural system identification with the following motivating questions: (1) Since neural data in the visual cortex are generated by recurrent circuits, would the use of a similar architecture provide a better model for predicting neural response? (2) What are the computational rationales and network features responsible for the performance advantages of recurrent circuits? (3) Finally, could neural response prediction using a network model with recurrent circuits provide a useful approach for characterizing and understanding the underlying recurrent neural circuits in the cortex?
Our task is to perform neural system identification by training neural network models with recurrent circuits to predict the mean firing rates of recorded macaque early visual cortical neurons (V1 and V2) in response to natural input images, when the animals were performing a fixation task. We separate the model into two stages, one stage including feedforward convolution layers and recurrent circuits to generate the inputs for the target neuron to be modeled, and a second stage, called readout, that models processing by the target neurons and temporal average of their responses to generating the target signals to be predicted. Three sets of data are tested, with distinct stimulus sizes and presentation duration (see Section A and Figure 1 of Appendix for details), to characterize the conditions under which recurrent circuits could exhibit performance gain over feedforward circuits. Neural response prediction models are sensitive to hyperparameters. To establish the computational advantages of the recurrent circuits over feedforward network in general, rather than specific to some optimized hyperparameters and architectural choices, we compare each network with recurrent circuits against its hyper-parameter and parameter matched feedforward model, across a large number of hyperparameters for model size (width and depth), activation functions (ReLU or softplus), model initialization seeds, and the amount of training data. In addition, we investigate four readout paradigms designed to model how the targeted neurons process the temporal responses of the hidden units from the recurrent circuits as well as the computation involved in creating the targeted signals for the network to fit.
Our contributions are three-fold. First, our experimental results conclusively show that recurrent networks perform better than hyperparameter-matched feedforward networks for predicting the responses of neurons in the early visual cortex. This performance gain is general, and not specific to a variety of hyperparameter choices, but more significant for more limited data size. This performance gain could arise from at least two possible sources: recurrent network (1) being a better or more efficient function approximator or (2) being a better match of the underlying circuitry in the early visual cortex, or both. Our second contribution addresses the first possibility by proposing that recurrent circuits can be conceptualized as a multi-path ensemble, reaping the benefits of ensemble computing [12][13][14]. We demonstrated empirically multi-path ensemble model approximates well the performance of recurrent circuits. This supports the conceptualization but also allows us to dissect the recurrent circuits from the multi-path ensemble perspective. Our third contribution addresses the second possibility by showing that the hidden units of the recurrent circuits exhibit response temporal dynamics and contextual modulation observed in early cortical neurons even when the training stimuli and targeted signals were not designed to create these dynamics and contextual modulation effects. Our findings provide fresh insights to the rationales underlying the computational advantages of recurrent circuits in neural networks and possibly in the brain, and suggest that the response prediction paradigm potentially could be useful for characterizing cortical recurrent circuits in the brain.

Related work
Recent work has demonstrated that adding biologically-inspired recurrent circuits to feedforward CNNs can lead to performance comparable to that of the state-of-the-art in image classification, but with significant less model parameters [9][10][11]15]. This suggests that recurrent circuits may provide a better or at least more efficient functional approximation. Here, we consider this approach for tackling the problem of neural response prediction. Rather than seeking state-of-the-art models or performance, we seek to assess whether recurrent networks are better than hyperparameter-matched feedforward models across many hyperparameters -where each recurrent layer is equated to two feedforward layers, to study the performance gain of recurrent computation in a more pure form.
In earlier work that incorporated recurrent circuits into CNNs for image classification, it was also found that the recurrent hidden units can exhibit temporal response dynamics reminiscent of that IT neurons [10] and the learned weight matrices of the recurrent kernels exhibit patterns that might give rise to known physiological responses [16]. These works however did not explicitly perform "neurophysiological experiments" on the hidden units to see if the neurons actually exhibit classical contextual modulation that presumably emerge from computation in recurrent circuits. Inspired by the works of [17] on retina, we performed some neurophysiological tests on the hidden units in the recurrent circuits trained to fit V1 and V2 responses. Our work differs from [17] in that we fit network to predict mean response to static images, rather than temporal response to dynamic video. While it is possible to use a feedforward network to predict temporal responses to static images, it inherently puts the feedforward network to a competitive disadvantage relative to the matched recurrent network which has intrinsic dynamics. It is more compelling to let recurrent network compete in a problem in which feedforward network has been shown to be able to achieve the state-of-the-art performance, than in a problem where a feedforward network is not known to do well. Interestingly, even when we trained the network to fit mean responses, we found the hidden units exhibit response dynamics that are similar to what we observed in V1 and V2 neurons to static images.
Although recurrent networks have been shown to produce performance improvements in a number of goal-driven tasks, the theoretical reasons underlying their success are not well understood. Deep equilibrium theory Bai et al. [18], showing weight-tied deep recurrent networks converge to a fixed point of the recurrent dynamics, provides an important perspective, but it has not helped us understand performance of the recurrent circuits in our setting. Here, we propose a different perspective, in the form of multi-path ensemble computing, to reason about recurrent computation. Treating recurrent computation as a multi-path ensemble is novel even though multi-path ensembles have been evoked to explain deep networks with residual connections [19,20], [21], [22,23], [24], and DenseNet [19]. Liao and Poggio [25] did suggest the equivalence between a specific type of recurrent networks with weight-tied ResNets, and Chen et al. [26] formulated deep weight-tied ResNets to approximate the dynamics of continuous recurrent models specified in ordinary differential equations. Our work explicitly connects recurrent models to multi-path ensemble models and explores the use of this reformulation to characterize recurrent networks in terms of path length statistics to understand the differential performance of the different models.

Comparison between the recurrent circuit models and feedforward models in neural response prediction
Our first objective is to evaluate whether and what recurrent network models perform better in predicting the mean firing rate of individual neurons than feedforward networks of matched model parameter sizes and hyperparameters.
Datasets: In this neural system identification problem, We train each model on three datasets independently. All three sets contains V1 and V2 neurons' responses to natural images recorded using multi-electrode arrays from macaque monkeys performing a fixation task. We call them the RF-S, WF-S, and WF-L datasets. RF-S's stimuli were presented in an aperture close to the size of the classical receptive field (RF) and in short (S) duration (60 ms) each in rapid succession. The RF-S data are public-domain data under a Creative Commons license, from [5] contains 115 V1 neurons' responses to 7250 images. The WF-S and WF-L data were collected by us, containing 79 and 34 V1/V2 neurons responding to 8000 and 2250 images respectively in response to stimuli presented in 8 o and 10 o visual angle aperture (wide-field (WF)) respectively. The presentation duration of each stimulus is 47 ms for the WF-S dataset, using the same paradigm as in [5]), and 500 ms for the WF-L dataset. Each stimulus was repeated about 10 times in both of these sets. These datasets are characterized by different extents of spatial and temporal contextual modulation. WF-L involves a large spatial context and long temporal processing for each image, thus significant contextual modulation is expected. WF-S involves extensive spatial context but does not allow much temporal processing. RF-S, with minimal spatial context and temporal integration, is expected to involve the least amount of contextual modulation. We expect the recurrent circuits learned will strongly depend on the degree of contextual modulation. Note that the WF-L and WF-S datasets include both V1 and V2 neurons.
Feedforward Baselines: Our baseline feedforward models are based on those in Klindt et al. [2], Cadena et al. [5] (Fig. 1a) which are state of the art in predicting neural responses of early visual areas. These models consist of a "shared core" of convolutional layers used for the whole dataset, followed by a factorized readout for each individual neuron. We trained feedforward models of many different model sizes and hyperparameters; each baseline model is a feedforward model (Fig. 1a) that contains a batch normalization [27] layer, a set of convolutional processing blocks (CPBs), as well as pooling, fully connected, and activation layers. The hyperparameters varied include amount of training data (25%, 50%, and 100%), activation function (ReLU or Softplus), model size (number of layers/channels), and others, resulting in at least 288 feedforward models for each dataset.
Recurrent Models: For each FF model, we will compare it with recurrent networks (with the same hyperparameters, and parameters) with the simplest recurrent interaction, in which a weighted sum (by a lateral interaction kernel) of all the neurons in a neuron's 3 x 3 neighborhood at each iteration is added to its response in the previous iteration. Since the CPB's kernel is 3x3xk where k is the number of feature channels, the lateral interaction kernel is thus also 3x3xk. From the perspective of model size, each Recurrent Convolutional Processing Block (RCPB) in the recurrent model is equivalent in the number of parameters to two convolution processing blocks (CPB) in the FF baseline model as shown in (Fig. 1b,d). A recurrent convolutional processing block (RCPB) extends regular CPB by having two inputs: a bottom-up input from the previous stage of the model in the current iteration and a lateral input from the output of the RCPB in the previous iteration (Fig. 1d). The batch normalization parameters are learned separately for each inference iteration, following the design of Spoerer et al. [15]. We experimented with networks with one RCPB and with two RCPBs. Two-RCPB models are parameter-matched to feedforward models with four CPBs. During model inference, information flows over each CPB only once, whereas it flows over each RCPB for a certain number of iterations. The number of iterations allowed is a model hyperparameter; when the number of model iterations is 1, a recurrent model is effectively a feedforward one. For each FF model and for each type of recurrent model architecture type, we trained six recurrent models, one for each number of iterations, with the number of iterations range from 2 to 7. Note that the model with 1 iteration is effectively a feedforward model but with less parameters than the feedforward model being compared, and thus will not be included in the comparison.  Variations in Architecture and Readout Modes: For both feedforward and recurrent models, the model inference can be divided into two phases: convolutional and readout phases (Fig. 1). The convolutional phase contains a BN layer followed by (recurrent) convolutional processing blocks, each of which contains convolution, BN, and activation layers; the readout phase takes the output of the final convolutional block through average pooling, factorized fully-connected [2], and activation layers to get the final model output. The readout phase represents the operation performed by the neuron to be modeled on the bottom-up input signals provided by the convolution phase, as well as the operation performed by the experimenter, averaging the temporal response in our case, to produce the teaching signals to be fitted by the model.
There are four types of readout paradigms in the readout phase for each of the recurrent models (Fig 2), and the learnable parameters are the same for the four readout modes, as well as that of the feedforward counterpart. For a model with T iteration, its RCPB running for T iterations will generate T intermediate outputs. The target neurons integrate and process these outputs differently in the four readout modes. no-avg, simply uses the output after the T iteration, as typical for convolutional RNNS [11]. early-avg averages the RCPB's outputs across time and feeds the mean to the readout phase after the T iteration, as was done in some bio-inspired RNNs for object recognition [15]. These two readout modes do not generate a temporal response at the ACT stage of the readout phase. In both late-avg and 2-avg modes, the target neurons process the output of the RCPB units in "real-time" at every iteration and produce an output at the ACT stage at every iteration, hence generating a temporal response for the targeted neuron. The temporal responses of the targeted neurons are then averaged to produce the mean response to match the teaching signals. Thus the Ave operation at the end models the temporal averaging that had been performed on the neural data. The main difference between them is that the target neurons process the output of the RCPB independently at each iteration in the late-avg case, but process the cumulative average of the output of the RCPB up to that iterations in the 2-avg case. Thus, target neurons in 2-avg remembers the history of the RCPB's responses for each stimulus up to that iteration when generating their responses at that time point. The four readout modes represent our hypotheses or assumptions on the computations underling the generation of the experimental neural data. Neural response prediction performance can serve as a metric to evaluate these assumptions.
In summary, for each feedforward model, there are 24 hyperparameter-match recurrent counterparts ( six numbers of iterations (2 to 7), and four readout modes). Thus, for the 288 feedforward models tested for each dataset, there would be 6912 recurrent models to be compared against for each dataset.  Model training and evaluation: For each data set, roughly 64%, 16%,and 20% of images were allocated for training, validation, and testing, respectively. The validation set was used to control early stopping of the models, since we did not select hyperparameters. We constructed and trained all models using PyTorch [28] on an in-house computing cluster with 4-GPU (NVIDIA GeForce GTX 1080 Ti or similar) nodes. and followed Klindt et al. [2] on data preprocessing and model optimization details. Given a trained model, we use average CC 2 norm [4,29] over all neurons and all testing images (20 % of all samples) to quantify the model's performance on each data set. Conceptually, CC 2 norm measures the fraction of variance explained by the model, with trial-to-trial variance in neural responses discounted. To get the CC 2 norm for a neuron, we first compute the raw Pearson correlation CC raw between the ground truth time-and trial-averaged neural responses and the model responses, then we divide CC raw by CC max , which estimates the maximal Pearson correlation coefficient an ideal model can achieve given the noise in the neural data [29,30], followed by squaring.

Multi-path reformulation of recurrent circuit models
After demonstrating networks with recurrent circuits in the appropriate readout modes consistently outperform their hyperparameter-matched feedforward counterparts, we develop an approach to characterize and understand the computational rationales and network properties that contribute to the performance gain provided by the recurrent circuits in this task.
We propose that a recurrent circuit can be conceptualized as a multi-path ensemble, with each iteration generating additional and longer paths. Ensemble computing is known to have the so-called diversity bonus, allowing a network to combine many weak classifiers or approximators to build a stronger classifier or approximators [12][13][14]. We proceed to reformulate the recurrent model's recurrent processing blocks (RCPBs) into an approximately equivalent multi-path ensemble with multiple feedforward paths of different lengths. The weights of the kernels across the paths are shared in order to keep the model parameters constant. Fig. 4 shows a one-RCPB network's information flow path at the first three iterations. With each iteration, the network adds a longer path. Fig. 4b approximates information flow as summations of simpler individual feedforward paths with shared components. Decomposing the operation of the recurrent circuit into into a multi-path ensemble allow us to dissect the recurrent circuits by eliminating paths as well as compare the path statistics of the different models or readout modes to understand their performance difference. We found two path statistics particularly illuminating. The first is the weighted average lengths of the paths of the multi-path ensemble, and the second is the diversity of the distributions of the weighted path lengths. Each component box in Fig. 4 can be associated with a certain weight, which is essentially a gain factor, representing the ratio of the magnitude of each box' output over that of its input. Parameters learned along different paths have different magnitudes which should be related to the relative contribution of the different paths to the final output. The product of the weights of all the box components along a particular path gives the total weight of that path. See Appendix F for details and an example of the weight computation. With these weights, we can compute the effective path length and path distribution diversity of the multi-path ensemble for the different model architectures.

Evaluation of response dynamics and contextual modulation in the hidden units
Our third objective is to determine the extent to which the recurrent circuits learned by our model produce response properties known to belong to recurrent circuits in visual cortex. We evaluated three features of the recurrent circuits learned by our recurrent models: their average temporal dynamics, the extent of end-stopping they show, and the extent of surround suppression they show. Neurons in early visual cortex typically show a rapid, high responses to the onset of a static stimulus, which then decays and plateaus at a lower level [31]. V1 neurons additionally are known to exhibit end-stopping (suppression by oriented bars longer than the preferred length) [32] and surround suppression (suppression by sine-wave gratings larger than the preferred size) [33]. We present oriented bars and gratings of varying sizes to our trained recurrent models and extracted the responses of centrally-positioned hidden units in the RCPBs, whose place in the recurrent circuit corresponds to superficial neurons in V1 or V2. We then examine the temporal dynamics of these responses, and their tuning to the length of bars and size of gratings (at optimal spatial frequencies and orientations), averaged over 150-200 units across several hyperparameter values. See Appendix C for more details.

Recurrent models outperformed similarly-sized feed-forward models
For models trained for each of the three datasets, we compare the feedforward models with their 24 types of parameter-matched recurrent counterparts (4 readout modes, 6 iteration numbers). Figs. 3B shows three plots. Each plot summarizes the performance of the 24 model types for each dataset, trained with 100% training data. The averaged performance of all the tested feedforward models is indicated by the horizontal solid line with the dotted line indicating 1 standard error. The readout modes are indicated by different colors. The mean performance of 48 recurrent models of this type of models with 7-iteration and 2-ave read-out mode is circled in the WF-L plot, with the vertical bar indicating standard error. The performances of the 48 recurrent models correspond to the 48 green dots in Fig. 3A. Each green dot compares the performance of one recurrent model of this type, trained on 100% of the training set, against its parameter and hyperparameter-matched feedforward model, i.e. the recurrent model with 1 RCPB is matched with a feedforward model with 2 CPBs, with the same activation functions, loss, etc. Fig. 3A also shows the prediction performances of these models when only 50% and 25% of the training sets were used, revealing the percentage improvement is more significant when training data is more limited. We perform pairwise T-test to evaluate the significance of prediction performance gain for each model type for different amount of the training data, and found the performance gain is significantly greater for the WF-L-trained models (performance gain is 11%, 6% and 5% for 25%, 50% and 100% training data) than the WF-S-trained models (performance gain is 6.5%, 3.5% and 2% for 25%, 50% and 100% training data). The performance gain for RF-S trained model is significant but small. In short, everything being equal, the 1-RCPB models and the 2-RCPB models performs consistently better than their 2-CPB and 4-CPB feedforward counterparts respectively, except for the 1-RCPB models trained for the RF-S data. For completeness, we showed all 24 panels for each of the three datasets in Appendix L.
Across all three datasets, we observe that the 2-avg (red curve) and late-avg (green curve) readout mode models consistently increase for models of more iterations, and outperform the other models for higher iteration numbers. We suggest this is likely because the temporal averaging of the neural response is explicitly incorporated in these readout modes.

Multi-path approximations work and reveal varying effects of paths of different lengths
To verify the approximation of the multi-path ensemble to the original recurrent models, we compare their neural response prediction performance. Figure 4c shows that their performance metrics are highly correlated across hyperparameter values (r=0.88). Furthermore, the effective average path lengths of the recurrent models and that of the multi-path reformulation are also highly correlated (r=0.96). This gives us confidence that the approximation is good and so we can draw conclusions about the corresponding recurrent models by analyzing the multi-path ensembles.  The multi-path ensemble reformulation allows us to analyze and dissect the recurrent circuits from the perspective of path-lengths and path distributions. Figure 5 panels A and B recast the performance graphs of Figure 3B as a function of effective path lengths instead of iteration numbers for the WF-L and RF-S datasets (WF-S results are somewhere in between and shown in Appendix I. The graph showed that at approximately the same averaged lengths, as indicated by the circled data points, the two readout modes that explicitly models the temporal averaging of the responses at the end, fare better than the early-avg, and far better than the no-avg models.  A and B show average performance as a function of path length for 2-layer recurrent models on the WF-L and RF-S datasets, respectively. C shows the results of ablating paths of one, two, or three lengths for models trained on WF-L (top) and RF-S (bottom), with shorter paths mattering more for the RF-S data (additional results in Appendix K). Panels D and E show distributions of path lengths for the best-performing iteration numbers of each readout mode on the WF-L and RF-S datasets respectively. The distributions are more biased towards shorter paths for the RF-S models, though in all cases no-avg has a strong bias towards the longest path.
Given that the notion of diversity is important for ensemble computing, we investigate whether the distribution of the contribution of paths of different lengths, beyond just the average length, plays a role in the differential performance of the four readout modes. Figure 5 panels D and E show the distribution of contribution by the paths of different lengths for each readout mode at similar averaged path lengths (the circled points in Figure 5A and B). We found that the better performing model architecture tend to be more evenly distributed in path-lengths, with a stronger emphasis on the earlier iterations particularly visible for the RF-S dataset. The contribution of the longer paths becomes more significant with the increase in the spatial context of the stimuli and in the duration of the stimulation, which would allow a greater amount of contextual modulation. This finding is also supported by the ablation studies in panel C, the performance of multi-path models retrained with paths of one, two or three consecutive lengths removed (since removing just one length has a very small effect). For the RF-S dataset, ablating the very shortest paths produces the largest performance decrease, while slightly longer paths (i.e. lengths 2 and 3) are more important for the WF-L trained models. Results for models trained on subsets of the WF datasets containing only V1 neurons, and the full WF-S data, in Appendix I and J, lie somewhat in-between, suggesting that spatial and temporal context are important contributors to the models of neural response.   Figure 6A shows the averaged temporal response to sine-wave grating stimuli; the initial burst followed by a decay to a steady state of the 2-avg matches visual cortical neurons [31]. Two well-observed physiological effects that have been thought to arise from recurrent circuits are end-stopping and surround suppression. End-stopping is the reduction of neural responses to a bar of the optimal orientation of the neuron when the bar's length exceeds the size of the receptive field or is longer than a preferred length [32]. Surround suppression is the reduction of neural responses when the size of the sine-wave grating of a neuron's preferred orientation is larger than the receptive field, encroaching on the receptive field's near and/or far surround [33]. Figure 6B and C display the average tuning of model units to bars and gratings as a function of length/size, showing the expected rise-and-fall for 2-avg units, but not for early-avg and late-avg respectively. In addition, the average RF size of about 1.2 degrees of visual angle (9 pixels in the models' input) is roughly the expected RF size of V1 neurons at our eccentricity, and the average level of surround modulation, compared to that of real neurons in [34], is quite close (Figure 6D and E).

Discussion
In this study, we showed that recurrent models with the appropriate computational architecture significantly outperform feedforward models of matched size and hyperparameters in scenarios where neural responses are significantly modulated by the spatial and temporal context. The performance gain is significant and consistent, particularly when the readout mode takes into account the experimenter's temporal average of the neural responses in constructing the target responses. Furthermore, it is reassuring that the hidden units in 2-avg models exhibit similar temporal dynamics and contextual modulation as observed in early visual cortical neurons in neurophysiological experiments, even though the models are trained only to predict mean responses with no constraints imposed explicitly to bring about these effects. Taken together, these results suggest that neural response prediction can provide a useful paradigm for characterizing the underlying recurrent circuits of the visual cortex.
By conceptualizing and reformulating recurrent networks as multi-path ensembles, we suggest recurrent circuits might derive their computational advantages from the notion of ensemble computing. More importantly, the reformulation allows us to perform analysis of the contribution of different paths based on path-length statistics and ablation analysis. The results of our analysis suggest models with a more uniform distribution of path-lengths and paths of longer lengths are needed to model greater amount of contextual modulation. The recognition that recurrent circuits can be conceptualized as multi-path ensembles might suggest a new direction of research to investigate how more powerful multi-path ensemble can be constructed and optimized by using recurrent network, both for machine learning and for understanding the brain.
There are several limitations to our work. While using small models is standard for neural response prediction, and allows us to explore a wide variety of hyperparameter combinations, the method of multi-path reformulation developed for such shallow models might be difficult to transfer to the much deeper models generally used in computer vision, or to other domains in which recurrent networks are applied. This means that the positive societal benefits of our work resulting from greater model interpretability may be limited, but we do not foresee any negative impacts from this formulation or better characterization of early visual cortical circuits. The hyperparameters we sweep over may not include those that result in optimal performance, and architectural improvements like the Gaussian readout of [35] may reduce the improvement from recurrent connections that we find. In addition, our recurrent circuit models are far simpler than the true circuits in visual cortex, which include a detailed laminar structure, extensive feedback, dynamically varying connectivity, and more. Thus, there is a limit to how well neural circuits can be characterized by our current models. [

A Details of neural data sets
The WF-S and WF-L datasets were collected in our laboratory from an awake behaving macaque monkey performing fixation task using an SC96 array, a 96-electrode multi-electrode array (Gray Matter Research, MT), implanted over the V1 operculum. Each set was obtained from 5 days of recording, with 4000 trials per day.
The RF-S is the dataset from Cadena et al. 2019 Plos Computational Biology, collected using 32-channel linear multi-electrode array. The data are available publicly (https://doi.gin.gnode.org/10.12751/g-node.2e31e3/) under a Creative Commons license, obtained and used in this paper with permission from the owners.  were presented in rapid sequence with each stimulus presented for 47 or 40 ms respectively for WF-S and RF-S. Each stimulus was presented for 500 ms per trial for WF-L, and responses to 2250 images (8-10 repeats each) were collected across five days. Only neurons that could be tracked and exhibit consistent tuning to daily calibration images across those days were kept for analysis.

B Details of the explored models B.1 Baseline models
During inference, a baseline model takes an image x ∈ R H×W of height H and width W as input and generates an prediction r ∈ R N of the neural responses r ∈ R N to the image. Mathematically, the model inference process is defined by Eqs. (1). r = Act(FC(Pool(y (M ) ))) (1a) For a model of M convolutional processing blocks (CPBs), the model inference starts with Eq. (1c) to normalize the input image and obtain the initial CPB's output y (1) , followed by a few applications of Eq. (1b) to get outputs of later CPBs y (2) , . . . y (M ) . Finally, Eq. (1c) is applied to obtain the predicted neural responses r. Conv, BatchNorm, Act, Pool, FC represent convolution, batch normalization, nonlinear activation, average pooling, and factorized fully connected layers [2], respectively. Different convolutional processing blocks have different hyper parameters and learned parameters, as denoted by superscripts (m) in Eq. (1d).

B.2 Recurrent models
The model inference process is defined by Eqs. (2). and the detailed information flow for a recurrent model across iterations are illustrated in Figure 8.

C Neurophysiological experiments
We tested the hidden units in the recurrent layer block in our models using oriented bars and sinewave gratings, as commonly used in visual neuroscience, to evaluate whether these units exhibit end-stopping and surround suppression effects, even though the models were not designed or trained for these effects. We extracted the responses of units in the center position of the RCPB, so their receptive fields were located near the center of the image.

C.1 End-stopping
First, we determined the precise RF location and preferred orientation of each unit with oriented bars of an intermediate length (10 pixels) presented in a 5 × 5 grid of locations around the center of the image, spaced 5 pixels apart. The bars spanned the range of 0 to 360 degrees of orientation, equally spaced at 22.5 degree intervals, for a total of 16 orientations, and were all 2 pixels wide, black against a white background. Using the orientation and position of bar that maximized the response of each unit in the first iteration, we presented bars varying in length from 5 to 45 pixels, at intervals of 2 pixels. For the analysis in Figure 6 of the main paper, we included only units that had a stable response of at least 1.0 across all time bins to at least one of the lengths, since some units were observed to drop to 0 or very low values. About 60-80% of units across models had these stable responses.

C.2 Surround suppression
We determined the orientation tuning of units to sine-wave gratings separately from that to oriented bars, using gratings that covered the whole image. However, we used the receptive field locations determined by responses to oriented bars for this experiment. Gratings varied in frequency from 0 to 15 cycles per image (roughly the maximum achievable with the low-resolution images the models were trained on) at steps of 1 Hz, and the optimal frequency was selected for each unit. Then, gratings of increasing diameter, from 5 to 45 pixels stepping by 2, were presented at each unit's preferred location and frequency. The gratings were presented on a gray background, with slight Gaussian blurring at the end of the aperture to reduce edge effects, and at 30% contrast, since surround suppression is generally observed at lower contrasts, and this value resulted in the closest average suppression index to previously reported neurophysiological data. Like for end-stopping, we only included units with a stable response to at least one of the grating sizes, which included about 60-80% of units.

C.3 Temporal dynamics
The average temporal dynamics were computed over all units' responses to all orientations and frequencies of the full-image sine-wave gratings used to determine orientation tuning for the surround suppression experiment. Dynamics in response to oriented bars looked the same. Figure 9 shows the results of 2-avg models trained on the RF-S data on our experiments (panels A and B), as well as the feedforward responses (first iteration) of models trained on the WF-L data. In all of these cases, end-stopping and surround suppression are not shown. The lack of appropriate contextual modulation in the first iteration of the response shows that the effects arise as a result of the learned recurrent circuit, and are not just a property of the feedforward part of the RCPBs. The lack of modulation in the RF-S-trained models suggests that neural responses exhibiting surround modulation (like the long-presentation, wide-field WF-L data) are necessary to learn the appropriate circuits via neural response prediction.

C.4 Feedforward responses and RF-S trained 2-avg models do not exhibit V1-like responses
C.5 2-recurrent-layer WF-L trained models exhibit similar responses as 1-layer ones Figure 10 shows the responses of hidden units from the second recurrent layer of two-recurrent-layer models on our experiments. The temporal dynamics, end-stopping, and surround suppression are still present for WF-L-trained models, in their last iteration after recurrent processsing. Figure 9: Results on the end-stopping and surround-suppression neurophysiological experiments for the 2-avg models in circumstances where they fail to reproduce V1-like effects. Panels A and B show the results for 2-avg models trained on the RF-S data instead of the WF-L data of the main paper's Figure 6, using the last iteration as in the main paper, for end-stopping and surround suppression respectively. Panels C and D show the results for 2-avg models trained onthe WF-L data, but using the very first iteration of the response, before recurrent connections have come into play. The structure of the plots is as described in the caption of Figure 5 in the main paper. In all cases, end-stopping and surround suppression are not present. Following practices in the literature [2,5], for each data set, we used cropped and downsampled images as the input and normalized average neural responses as the output for modeling purposes.

D.2 Model configurations
To comprehensively evaluate recurrent models vs. feed-forward models, we trained and tested multiple variants of recurrent and feed-forward models with different model configurations.
Kernel sizes Following previous studies [2,5] and as well as based on our pilot experiments, in either of a feed-forward model or a recurrent model, the first convolutional layer has a relatively large kernel size (9) and subsequent layers have a small kernel size of 3. The pooling layer always has a kernel size of 3 and a stride of 3.
Feed-forward models For each data set, we trained a large set of feed-forward models with different configurations differing on aspects like amount of training data, model size, loss functions, randomness parameters, etc., as listed in Table 1. Each aspect of configuration was explored for certain reasons as follows.
• Training data and model size-related model hyperparameters. we explored model configurations in these two aspects because we feel training data amount and model size affect the relative performance difference between feed-forward and recurrent models. In particular, we hypothesized that the additional circuit priors from recurrent models would be more useful when there was less training data and the model size was larger. • Model size-independent model hyperparameters and loss functions We explored model configurations in this aspect because we wanted to know if the advantage of recurrent models over feed-forward models is limited to certain choices of model size-independent model hyperparameters and loss functions, or if such advantage is universal across choices of model size-independent hyperparameters and loss functions. • Randomness hyperparameters We explored different model initialization seeds because we wanted to obtain robust metric numbers for our models by averaging results over different seeds. Ideally we should also explore over other randomness hyperparameters, such as the seed used to split data into training, validation, and testing sets. We did not do so due to computational resource constraints. Recurrent models For each feed-forward model, we trained 24 corresponding recurrent models of the same size, to comprehensively compare feed-forward and recurrent models. To keep the comparison as fair as possible, each of the 24 recurrent models have the same hyperparameters as the feed-forward one on those hyperparameters listed in Table 1, except with a decreased number of convolutional layers to match the model size because a recurrent convolutional processing block has roughly twice as many parameters as a convolutional processing block with the same number of channels. In particular, we replace every two CPBs of kernel size 3 in a feed-forward model with one RCPB of kernel size 3 in a recurrent model. Note that a recurrent model does not involve recurrent computation in its first CPB, which has a large kernel size of 9, and only has recurrent computation at subsequent CPBs of kernel size 3; we made this decision based our pilot experiments and for simplicity in matching the model size between recurrent and feed-forward models. In addition, the 24 recurrent models differ in their recurrence-related model hyperparameters listed as follows.
• Number of iterations, with six possible values 2 to 7. This hyperparameter affects the amount of recurrent computation during model training and inference. We wanted to explore this hyperparameter because obviously amount of recurrence may affect model performance. Note that we always use the same number of iterations during training and testing; e.g. a recurrent model with T iterations during training is always evaluated with T iterations during testing.
• Readout mode, with four possible values (no-avg, early-avg, late-avg, 2-avg. We want to explore this hyperparameter because we wanted to know if some readout modes perform better than others and if there exists some optimal way to make use of intermediate outputs generated by a recurrent model.

D.3 Model optimization
We trained all models using PyTorch [28] and followed Klindt et al. [2] for model optimization.
Objective function For each model, the objective function to be minimized was the sum of two parts.
Neural prediction loss Depending on the choice of hyperparameters, the loss can be either mean squared loss or mean Possion loss between between the predicted neural response and recorded ground truth averaged across neurons and images. Note that for recurrent models with late-avg or 2-avg readout modes, following the practice in Spoerer et al. [15], we computed the neural prediction loss by averaging the losses over individual iterations. Regularization terms We applied L1 sparsity penalties as regularization; in addition, for the first convolutional layer of the model, we applied smoothness regularization as in Section 5.2 of Klindt et al. [2].
Optimization algorithm Given the objective function, we used the same optimization algorithm as in the public code of Klindt et al. [2]; in particular model parameters were optimized sequentially in three phases, by three Adam optimizers with decreasing learning rates. Early stopping was applied in each phase, guided by the neural prediction loss evaluated on the validation set. For each data set, roughly 64 %, 16 %,and 20 % of images were used for training, validation, and testing, respectively. Due to computation resource constraints, we only created one particular split of training, validation, and testing sets for each dataset.

E Model performance evaluation in detail
Given a trained model, we use average CC 2 norm over all neurons to quantify its performance on a data set. For each neuron, we compute its CC 2 norm based on Eqs. (3).
Concretely, we first compute the raw Pearson correlation CC raw between the ground truth trialaveraged neural responses r as defined by Eqs. (3d),(3e) and the model responsesˆ r using Eq. (3b). We then divide CC raw by CC max , which is defined in Eq. (3c) and estimates the maximal Pearson correlation coefficient an ideal model can achieve given the noise in the neural data [29,30]; inside the square root of Eq. (3c), the numerator is the difference between variance of response sums per stimulus k r m,k and sum of response variances per trial Var({r m,k }), and the denominator is the variance of the trial-average neural responses scaled by a factor related to the number of trials K. Finally, we get the square of the normalized Pearson correlation coefficient using Eq. (3a). As squared CC raw gives the fraction of variance in neural responses explained by the model in a simple linear regression, squared CC norm gives the normalized explained variance that accounts for noise in the neural data. Note that we used neural responses on all stimuli instead of on testing set stimuli to compute CC max in Eq. (3c) for more accurate estimation. In addition, because the WF-L data set has a variable number of trials per image (8 to 10), we only used the first eight trials for simplicity.
While the derivation of CC max is technically complicated, it's positively related with the Pearson correlation between mean responses computed over half of trials and those computed over the other half [29,30]. We use CC 2 norm instead of raw CC 2 raw to quantify model performance because the former is more robust to trial to trial response fluctuations. Our main results in Section 4.1 still hold if we measure the performance in CC 2 raw .
F Path length calculations for the multi-path models The relative contributions of each path of a multi-path model can be quantified; summarizing these contributions by path length allows us to compute the average path length of each model, as well as its distribution, or diversity, as presented in Figure 5 of the main paper.   Figure 11: Computation of average path length and path diversity for an example recurrent model. The model has one recurrent block (RCPB) and uses the no-avg readout mode with three iterations. In (a), the relevant information flow is shown, with a strength assigned to each component. In (b), the original information flow is reformulated as a simple multi-path ensemble with three paths, and the unnormalized strength of each paths 1 ,s 2 ,s 3 is computed as the product of strengths along that path. In (c), the average path length of the model is defined as the average path length weighted by the normalized strengths s 1 , s 2 , s 3 , and the path diversity of the model is the discrete probability distribution implied by normalized strengths s 1 , s 2 , s 3 . The length of each path is defined as the number of convolutional layers in that path. Activation layers are omitted in (a) and (b) for brevity and they do not affect the calculation due to having a strength of 1.
G Performance versus iterations divided by number of recurrent layers Figure 12 shows performance of recurrent models as a function of number of iterations, with parameter-matched feedforward models' performance shown as a horizontal line. It is analogous to Figure 3 of the main paper, but divides the results up between models with one and two recurrent layers. Two-recurrent-layer models tend to further outperform their parameter-matched feedforward models, particularly for the RF-S data, where one-recurrent-layer models are no better than feedforward. H Path length distributions for one-recurrent-layer models Figure 13 shows the performance versus path length and path length distributions for models with one recurrent layer, analogous to the results for two-recurrent layer models shown in Figure 5 of the main paper. Ablation studies were not performed on these models. Figure 14 shows the performance as a function of path length and path length distributions for each 2-layer model readout mode trained on the WF-S data. Notably, the path length distributions fall somewhere in between those observed in the RF-S and WF-L datasets in the main paper's figure, with shorter paths having a higher impact than for WF-L models, but a flatter distribution than the RF-S models. The WF-S data has wide-field images that allow for surround modulation (like the WF-L data), but short presentations that limit the amount of modulation possible (like the RF-S).

I Path length distributions for WF-S models
J Multi-path analysis for V1-only WF-L/WF-S models The WF-L dataset (as well as the WF-S dataset) consist of both V1 and V2 neurons, while the RF-S dataset contains only V1 neurons. To see whether the differences between models trained on the WF-L and RF-S in Figure 5 of the main paper depend more on the stimulation paradigm Figure 13: Performance plotted against average path length, and path length distributions for onelayer recurrent models trained on the RF-S and WF-L datasets. For these models as well as the two-recurrent-layer ones, training on the RF-S dataset results in higher weighting of shorter paths.
(WF-S) Figure 14: Path length and diversity for models trained on the WF-S data. Analogous to panels A, B, D, E of Figure 5 in the main paper, but for the other dataset.
or the cortical area of the modeled neurons, we trained models on a subset of the WF-L dataset containing only neurons in V1, 6 out of the 34 neurons. Figure 15 shows the analyses of Figure  5 for both 1-recurrent-layer and 2-recurrent-layer models trained on this V1-only WF-L subset. The results generally resemble the results for the full WF-L dataset in Figure 5 of the main paper (2-recurrent-layer models) and Figure 13 (1-layer recurrent models), with path length distributions notably preferring longer paths than those for models trained on the RF-S data in Figure 5. However, the ablation study results are in-between, with the removal of mid-length paths having the largest effect when 3 paths are removed (like the full WF-L), but the shortest paths mattering most when only 1 or 2 paths are removed (more like the RF-S). This suggests that the differences we see between models trained on the WF-L and RF-S are substantially due to the difference in presentation paradigms, but in part driven by the different cortical areas modeled.
In addition, while we did not include WF-S-trained models in Figure 5 of the main paper, we replicated those results (included in Figure 14) on models trained on the subset of the WF-S data containing only V1 cells. The results, displayed in Figure 16, are fairly similar to those of models  Figure 5 and elsewhere, ablation studies were only performed on the 2-recurrent-layer models.
trained on all the WF-S neurons, including the V2 cells. This supports our claims that the difference we see in the importance of paths of different lengths between datasets is driven mostly, but not entirely, by the presentation paradigm of the stimuli than by the cortical area modeled. Figure 16: Performance versus average path length (A), path length distributions (B), and ablation results (C) for models trained on only the V1 neurons in the WF-S dataset. As in Figure 5 and elsewhere, ablation studies were only performed on the 2-recurrent-layer models. Comparable to Figure 14, but only training on V1 neurons.

K Ablation study
To compute the baseline shown as horizontal dashed lines in the main paper's Figure 5C and 17, we trained the multi-path models corresponding to the original recurrent models with 2 layers, 7 iterations, and 16 or 32 channels using 100% of training data. We used all the training data because the reformulation of recurrent models as multi-path models worked better with more training data, and we did not train recurrent models with other numbers of layers or channels because models with 3 layers or more than 32 channels exceeded GPU memory limits during training and 8-channel recurrent models did not outperform feed-forward models much. We also restricted this analysis to the 2-avg models, since the neurophysiological experiments and performance values highlighted those as the most relevant.
To compute the performance metrics of ablated models in the figures, we trained additional multi-path models like those for getting the baseline, but with paths of certain lengths removed. In each of the figures, the left panel shows the performance metric change when paths of a particular length (1,2,3,4,5,6) were removed; the middle panel shows the performance metrics change when paths of two adjacent lengths (1-2, 2-3, 3-4, 4-5, 5-6) were removed; the right panel shows the performance metric change when paths of three adjacent lengths (1-3, 2-4, 3-5, 4-6, 5-7) were removed. L Extensive comparison of recurrent models against feed-forward models

WF-S
The scatterplots in this section show the performance of 24 types of recurrent models against hyperparameter-matched feedforward models for the three datasets, split up by readout mode and number of iterations. Each point is the average performance value of two random initializations for a model with a particular set of hyperparameters. For the RF-S dataset, the scatter plots include both 1R-layer and 2R-layer models, like they do for the other datasets. Upon scrutiny, we found 1R recurrent models do not yield a statistically significant improvement on this dataset, but 2R recurrent models do, as shown in the histograms of Figure 23, where much of the mass is concentrated above zero. Results are also shown for models trained subsets of the WF-S and WF-L datasets containing only V1 neurons, for models trained on 100% of the data only. They continue to show the improvement of recurrent networks.       : Performance improvement of 2-recurrent-layer recurrent models against feedforward models trained for the RF-S dataset. Histograms of (R-F)/F performance metric for the 24 varieties of recurrent models. R is the average CC norm for the recurrent models and F is the average CC norm for the matched-size feedforward model.