) Distribution of variance ratio after one-way ANOVA analysis of the CNN L1 and L2 responses across samples and families from randomly selected 103 units. L1 responses are more variable across samples in the same family. L2 responses are more variable across families than L1. Geometric mean variance ratios are shown with triangle marker which are 0.18 in L1 and 0.29 in L2 and the difference is significant (, test on the log variance ratio). This trend of L2 versus L1 units is qualitatively similar to the properties of V2 neurons shown in , though the variance ratios are higher for the cortical data.
, a popular deep CNN model widely applied in computer vision and neuroscience. Here we refer to the network as AlexNet. We chose AlexNet as our base model, because it includes computations that are loosely matched to visual cortex, such as pooling and local (divisive) response normalization. In addition, the receptive field size ratio can be controlled roughly to match the V1 to V2 ratio. Our approach was to keep a single base model, and then manipulate the architecture in various ways. We later also considered some other hierarchical architectures in the literature. The CNN model was trained on the ILSVRC2012 ImageNet dataset, a popular large-scale image database. We then presented each of the texture images to the resulting model.
We took the layer outputs after pooling and normalization, referring to them as L1, L2 and so forth. To get a better handle on where exactly in the neural network compatibility with V2 first emerges, we also considered layer outputs at all other intermediate points of L2 in the network (Fig. 0(b)). This allowed us to better understand how the computational building blocks (e.g., convolution, pooling, normalization) in the CNN may give rise to the differences observed in texture selectivity between V1 and V2; in other words, at what point in the CNN there is a transition from V1-like behavior to V2-like behavior.
I-a Texture generation and neurophysiology data
The neurophysiology data for V1 and V2 is described in  and , with recordings from macaque monkeys. We used the synthetic textures of , which were generated from a set of 15 real texture images. Each synthetic texture imagewas generated using the approach of Portilla and Simoncelli . We refer to each set of textures generated from the same source image as family and all the images within the same family as samples. Naturalistic textures for a given family were generated each with a different random seed. Spectrally matched noise images (which we denote noise images) were generated by randomizing the phase of the synthetic images. The noise images have the same spatial frequency distribution of energy as the original ones, but lack the differences in higher order statistics.
Overall, the image set included 15 samples from each family, resulting in 225 texture and 225 noise images. We downsampled the textures so that the effective portion of the image that the CNN units are sensitive to is equated to the receptive field size of the cortical neurons. For more detail on this process and the texture generation, see Methods.
For the data comparison, we largely focused on a population metric of modulation index which is indicative of selectivity to the textures versus the noise. We initially made a qualitative comparison between the biology and model, and then developed a way to quantitatively compare the macaque data to the CNN model units. We also considered other qualitative comparisons between the CNN and the cortical data.
I-B Qualitative correspondence of the CNN to the neurophysiology cortical data
We equated the number of CNN neural units in our simulations to the number of units in the neural population as in  and  (102 V1 and 103 V2 neurons). For the CNN, we considered the total number of units as the number of channels in a given layer, times a center spatial neighborhood (see Methods). We then randomly selected 103 neural units as shown in Fig. 0(c).
We first asked whether there is qualitative correspondence between the CNN model and the experimental neural cortical data, and compared these with a number of metrics. This then prompted the quantitative approach in the next section.
I-B1 Visualization of the CNN population selectivity
To first gain intuition that L1 and L2 differ in their texture representation, we visualized the unit population activity. We transformed the responses of CNN layers from a high-dimensional space (where dimensionality is the number of units in the given CNN layer) to a 2-dimensional space to visualize the unit response properties of L1 and L2. We used the Barnes-Hut t-distributed stochastic neighbor embedding (t-SNE) ,  algorithm to achieve this visualization. t-SNE is a technique for dimensionality reduction that tries to model small pairwise distances to capture local data structures in the low-dimensional space.
each point results from embedding an image represented by a high-dimensional response vector into two dimensions. Therefore, we have a total of 225 points, that come from the same number of images from 15 texture categories. Each point represents the population response to a texture sample, and samples belonging to a same family share the same color. L1 responses are more scattered and do not group the same texture family together. This is apparent both when randomly choosing 103 units (Fig.1(a)); and when considering all units in an spatial neighborhood (Fig. 1(b)), amounting to a total of 3072 units. In the L2 response space, samples from the same texture family group more tightly together. This is less apparent when using 103 random units in L2 (Fig. 1(d)), but becomes striking when considering all units in a spatial neighborhood (Fig. 1(e)), amounting to a total of 2048 units. L2 therefore better separates the texture responses than L1, qualitatively similar to what has been shown for V2 versus V1 in the biological data . However, a larger number of units form the CNN are required to match the texture discrimination capabilities of the V2 population.
I-B2 Tolerance across versus within texture family
We also found that, with the same number of units, L2 responses are more tolerant than L1 to the variations in image features across samples within a texture family, qualitatively similar to the V2 versus V1 data (Fig. 1(c), 1(f)).
The variance ratio is calculated by taking the geometric mean of the ratio of variance across-family to the variance across-samples, with a large value indicating the high tolerance of neurons to the statistical variation of samples in the same family. More specifically, we ran a one-way ANOVA analysis of the CNN L1 and L2 unit responses across samples and families, for a set of randomly selected 103 units. We found that the L2 variance ratio (0.29) is significantly greater (, -test on the log variance ratio) than L1 (0.18). This gives an indication that L2 unit responses show stronger variability toward different families and higher tolerance to the samples within a same family. These observations are qualitatively consistent with the neurophysiology results reported in  for V2 versus V1, although the variance ratios are higher in the cortical data than in the CNN (0.63 in V2 versus 0.4 in V1). This indicates that the cortical neurons show stronger variability towards families than the CNN. A note of caution in comparing the results is that the cortical data analysis included a nested ANOVA with consideration of the noise that occurs from stimulus repeats. In the CNN, there was no source of noise for stimulus repeats and so we did not consider this component in the ANOVA.
I-B3 Differentiating L2 units from L1 in CNNs
To further show the distinction between L1 and L2 for texture sensitivity, we followed the approach in  for V1 and V2, to compute a modulation index metric. The modulation index captures the differential response between textures and noise. As for the data, the mean modulation index was computed for each of the 15 texture families, resulting in 15 mean modulation index values. This was done for each of the network layers, L1 and L2. We computed the modulation index from the responses of all samples from each family, both natural and noise, and averaged over the number of neural units in the respective layer.
The modulation index is defined as the difference of responses from the textures to noise for each neural unit, divided by their sum, according to (1):
where, and are the responses to naturalistic textures and noise respectively. Fig. 2(a) shows the average modulation index for all texture categories in the CNN, both for L1 (red) and L2 (blue).
Averages are obtained from 10000 repeats, where at each iteration, we randomly select 103 neural units and compute the modulation index.
High modulation index of a population of neural units towards a family means that this group of units are highly sensitive to this specific family; hence they show high differential response. Since first and second order statistics are matched for both natural and noise images, a differential response also means that units latch onto higher order statistical properties of the stimuli.
We found that the average modulation index of L1 neural units is near zero and the modulation index of L2 is substantially higher than L1. The diversity in modulation index for the different texture families is shown in Fig. 3, for both the CNN (Fig. 2(a)) and the biological data (Fig. 2(b)). The average modulation index of L2 (0.18) is higher than L1 (-0.04). The difference between the indexes of L1 and L2 is significant (, t-test considering signs; , t-test ignoring signs and considering only the magnitudes) and is qualitatively comparable, but stronger than the biological data (V1: and V2: ). More specifically, Fig. 2(a) shows that the L2 modulation is more drastic in some texture families than others, as also observed for the V2 data . However, the rank order of the textures was different between the CNN and the biological data (prompting our quantitative subset selection approaches below).
As a control to examine the influence of the CNN architecture alone, instead of using the model weights that resulted from training on the ImageNet database, we generated random weights (in the interval ) for the L1 and L2 layers and averaged over 100 iterations. While keeping the L1 weights as trained, randomization only in the L2 units decreased the average modulation index to 0.05 (from 0.18). The difference between the L1 and L2 modulation index still remained significant (, t-test). We found that randomizing both L1 and L2 significantly decreased the sensitivity of the L2 units to the textures versus the noise and yielded an average modulation index of -0.02 and 0.01 respectively for L1 and L2. This lead to no significant difference between the L1 and L2 modulation index (, t-test on the magnitude). We can see that, irrespective of texture types, L2 loses its distinctive behavior in the complete absence of the trained weights and shows similar modulation to L1. This indicates that the selectivity we report in L2 is non-trivial and the CNN model with trained weights corresponds better to the biological data than the CNN architecture itself.
I-C Systematic quantification of the CNN to the visual brain data
We have thus far shown some qualitative correspondence between the CNN and the biological cortical data. Fig. 4 illustrates that for a random selection of CNN units, indeed this correspondence is only qualitative. The mean modulation indices for the various textures in the V2 cortical data versus the L2 in the CNN are correlated but clearly do not lie on a straight line. We wondered whether there exists a set of 103 CNN units that can well fit the cortical data. That is, rather than considering all units in the CNN or some random selection of units, we posited that perhaps a subset of the units could better explain the subset of experimentally recorded V2 neurons. For the remainder of the paper, we focus on this question with respect to the modulation index metric, capturing sensitivity to the textures versus the noise.
We note that the subset selection aspect of this problem makes it different from a standard regression and from approaches we are aware of for fitting neural data, which often take a weighted average of the units , . Why search for a subset of CNN units that can fit the cortical modulation index data? Our rationale was that finding such a subset would suggest that at the population level, there is some overlap between the CNN units and the cortical neurons in their representation for the texture versus the noise stimuli, as explained below.
To show a systematic quantification, we probed the CNN to select a subset of 103 neural units that is most consistent with the cortical neurophysiology experiments (Fig. 0(c)). We considered several approaches for subset selection. First, we employed a greedy technique, which we call subset greedy, to choose a subset of 103 units that best match the data from the brain. Briefly, from the set of all possible neural units, the greedy approach chooses the first unit with the closest euclidean distance to the V2 mean modulation index data; then the second unit is added to this subset so as to minimize the euclidean distance; and so on until a total of 103 units are chosen (see Methods).
For comparison to our greedy fitting approach, we also applied an optimal weighted average or full population approach. The full population approach finds a weighted sum of the neural units (under the constraint that the weights sum to 1) that is the closest in squared Euclidean distance to the experimental data. Notice that the weighted average may include all available units and weight units differently. The greedy approach is, in contrast, an approximation that finds a subset of 103 units with equal weights that best matches the neural data.
Our rationale is that the full population approach shows the best fit one can obtain with units from a given layer. However, it does not show an actual population representation that matches the data; it only reveals a linear transform of the representation. The suboptimal method chooses a subset of 103 units and thus uses the actual CNN representation. This subset selection approach is therefore more comparable to the analysis in the cortical experiments, in which the modulation index is computed as an equally weighted average of the neural units. In addition to the greedy approach, we also applied another suboptimal model selection technique, which is a regularized version of the full population fit that selects 103 units and we termed it assubset regularized. See Methods section for more details about the model fitting techniques.
In the next sections, we show results for fitting the CNN neural population to the V2 texture data with these approaches. We find that the L2 population can well fit the V2 data, but that the L1 population provides a poor fit. We then proceed to cross-validated fits, showing that this main result holds when we train the subset selection on one set of texture and noise images and test on the left out images. Finally, we show cross-validation results for a wider range of neural network manipulations.
I-C1 L2 population fits are well matched to the V2 data compared to L1
In this section, we discuss the fitting results for the full data set. In the next section, we then discuss the equivalent results for the cross-validation.
We found that a subset of 103 L2 neural units exist that provide a good fit to the V2 neurophysiology data (Fig. 5; second row). In contrast, all three fitting approaches showed that for the L1 units, no such set exists that can fit the V2 data well (Fig. 5; first row). We found that even the full population optimally weighted L1 fit (Fig. 4(b)) could not fit the V2 biological data well. This indicated that the second layer, but not the first layer of the CNN, is better matched to the V2 data in terms of the sensitivity to textures versus spectrally matched noise.
We quantified the fits with the Euclidean error distance between the mean modulation indices in the neurophysiological data versus the modulation indices obtained from the CNN for each family. A smaller Euclidean distance indicates a better fit to the V2 data and higher correspondence to the brain (see Methods section for details). The rationale behind using the Euclidean distance as a measure of correspondence is that it is directly related to the root mean squared error (MSE) up to a normalization constant. Our optimal weighted and subset regularized fits are done in terms of squared Euclidean distance, which for the optimal fitting method makes the error and regularization terms work at similar scales. In the subset greedy approach, MSE and Euclidean (and even squared Euclidean) distances indicate the same outcome.
We obtained Euclidean errors of 0.03, 0.01 and 0.08 for the L2 subset greedy, full population and subset regularized approaches respectively. In contrast, the L1 fits yielded Euclidean errors of 0.61, 0.30 and 0.61 for the three fitting approaches. This poor L1 fit did not change if we took the outputs from any other stage in the first layer. We therefore found that the correspondence of the CNN with V2 for texture sensitivity emerges in the second layer of the CNN after the (ReLU) rectification, but that no stage in the first layer (even after pooling and normalization) could account for the V2 texture sensitivity data.
As we have seen in Fig. 5, the L1 units do not fit the V2 data, even using the full population. We also explored the L1 fit of the CNN unit responses to the V1 data. Both the subset greedy and full population selection could find a set of units from L1 that fits the V1 data well. In addition, we considered the size of the spatial neighborhood from which we chose the units. In the greedy approach, selecting L1 units from a center spatial area fit V1 better (error 0.02) than choosing them from a neighborhood (error 0.14). At the same time, selection from a larger pool of L1 units did not rescue the V2 data fit (error 0.44). This indicated that L1 can fit the V1 data much better than the V2 data. The converse was not true: That is, both L1 and L2 could fit the V1 modulation index data well, possibly since L2 inherits some properties from L1.
We wondered if the main result was peculiar to the  network that was pre-trained on the ImageNet database. We therefore trained the CNN on another popular large scale image database known as the Places365 database . We found a similar trend, namely that the L2 units were more selective than the L1 units. We obtained L1 and L2 Euclidean errors of (greedy selection: 0.63 vs 0.10; full population: 0.49 vs 0.01; and subset regularized selection: 0.62 vs 0.11).
I-C2 Cross-validating the L2 population fits
We have thus far fit the full set of texture data with the CNN. We found it interesting that the L2 but not the L1 units could well fit the V2 data. Indeed, L1 could not well fit the V2 data, even though the training and test data sets were the same. To test the generalization capability of our method, we ran cross-validation fits.
We obtained better cross validation results by extending the image dataset to a total of 225 texture and noise images for each family. We learned the population (e.g., of 103 units) using 224 texture and noise images from each family for the training, and made a prediction of the mean modulation index for the left-one-out set of 15 images (see Methods). Fig. 6 shows the exact same fits as for Fig. 5, but with the cross-validation test results. This reveals the same main trend as we found for the fitting of the full texture data without cross-validation: L2 (second row) but not L1 (first row) could provide a good fit to the V2 data. Euclidean errors for the subset greedy method were 0.20 and 0.22 for the training and test predictions, respectively. Considering the whole population, we obtained train and test errors of 0.15 and 0.19; and for the regularized fits we obtained errors of 0.20 and 0.23 respectively. Note that the train and test errors were rather close. The train errors were higher than the fitting errors for the original texture images (Fig. 5
), probably because of the added texture images to which we were fitting.
These fitting errors were all lower than the random selection of 103 units in the population that we examined earlier (compare to Fig. 4; Euclidean error 0.37). In addition to the Euclidean error, we also considered the explained variance (): this was 0.60 for the subset greedy, 0.54 for the subset regularized, and 0.70 for the full population. In contrast, the explained variance for the random population of 103 units was 0.40. Note that this represents a lower bound, since we are not considering the variability due to samples in a family, nor are we taking into account variability in the experiments due to stimulus repeats.
In the remainder of the paper, we report cross-validated results for all conditions. Fitting without cross-validation resulted in comparable results, but in some cases lead to potential over-fitting (for instance, the almost perfect fit of L2 to the V2 data). This was particularly the case for the full population fit, which has the added luxury of using any number of (weighted) units in the population. For this reason, we report all remaining results for the cross-validated fits.
I-C3 CNN population fits are reduced for the architecture alone, with random rather than learned weights
Given the good correspondence of the L2 neural units to the V2 data, we wondered at what point this fit breaks or can be reduced. In particular, we wanted to address the following: Is training on natural images important for the CNN model to develop texture selectivity? What is the effect of the CNN architecture itself? Can the CNN architecture already develop texture selectivity, just by stacking up multiple layers?
We therefore considered a control of taking random weights for the CNN layers, instead of weights learned from the supervised model trained on images. We generated random weights in the interval of for the L1 and L2 CNN layers. We then asked how well the L2 layer (with the random rather than learned weights) can fit the V2 data. We repeated this process and took the average of 10 iterations from the layer responses. We found that taking random weights in the CNN resulted in a larger error and therefore a reduced fit, compared to the trained network. This can be seen in Fig. 7. The Euclidean errors were 0.52, 0.50, and 0.53, respectively for the greedy subset, full population, and regularized subset techniques.
Fig. 8 shows a more comprehensive summary of the cross-validated fitting errors across a range of random weight controls. The figure first summarizes the main results on the trained weights, before showing the results for the various randomized conditions. As described in the previous subsection, the trained L1 neural units (layer 1) exhibited the largest fitting errors among all layers and controls, meaning they resulted in a poor fit (hence little correspondence) to the neurophysiology V2 data (see also Fig. 5(a), 5(b), 5(c)). The trained L2 neural units exhibited the lowest fitting errors in all fitting techniques (greedy subset 0.22, full population 0.19, and regularized subset 0.23), meaning that L2 achieved better correspondence to the V2 data (5(d), 5(e), 5(f)).
We exhaustively explored a range of controls for randomizing the CNN layer weights (Fig. 8), and therefore considering the influence of the architecture alone. Overall, assigning random weights to the network layers increased the cross-validated fitting error. Randomizing only layer 2 weights and keeping the trained L1 weights (rand 2) lead to comparatively much better fits than randomizing both layer 1 and layer 2 weights (rand 12). The rand 2 fits were much closer to the L2 trained model, indicating that training the first layer alone went a long way in obtaining a good fit; but was still significantly worse than the L2 trained ( in greedy; in optimal; in regularized; one sample -test).
We wondered if deeper random architectures could lead to better correspondence with the brain data. We therefore considered randomizing layer 1, 2 and 3 weights (rand 123) and randomizing layer 1, 2, 3 and 4 weights (rand 1234). For these conditions, we fit the outputs of layer 3 and 4 respectively, to the data. The goal here was to see if stacking more random layers helped in obtaining a better fit to the data. However, the error remained high even when we stacked together 4 layers (compare rand 12 with rand 123 and rand 1234). Therefore, a deeper random network did not rescue the fit.
We also asked what happens if the trained weights within each filter are shuffled to destroy any spatial correlations. This maintains the distribution of the trained weights in each of the filters. To test this, we spatially shuffled the trained weights for each of the filters in both layers (shuffle 12). We found that this resulted in a better fit than the randomized counterpart, but still remained rather poor (compare rand 12 vs shuffle 12). In this case, the fitting errors stayed slightly lower than the random because the network might benefit from having the weights come from the same distribution as the trained (albeit that the weights are scrambled). Nevertheless, the error remained high, even in the shuffled version.
These factors give a clear indication that training the model on image classification leads to a better correspondence with the brain texture selectivity data. Even one trained layer is able to influence the subsequent layer(s) to gain some cortical correspondence.
On one hand, this reveals the necessity of training the deep learning models on the natural image dataset beforehand to achieve a better match to the V2 texture sensitivity data (and also high recognition accuracies for that matter). Other studies have also indicated the necessity of model training. Indeed, the primate brain is also ”trained” on natural scene data (albeit not necessarily in the same supervised manner). Nevertheless, it is interesting that the architecture alone can partly account for the data (see also , ).
I-C4 Effect of normalization and different CNN computations in texture selectivity
An important question regarding the CNN is how the various computations influence the compatibility of the model to the data.
AlexNet includes a local response normalization in L1 and L2, whereby each unit response is divisively normalized by the responses of 5 units (including the self) that spatially overlap. This loosely mimics the divisive cross-orientation suppression in cortical V1 neurons . We therefore wondered: What happens to the selectivity if we ignore the local normalization of L1 (norm1) and L2 (norm2) layers altogether?
We first considered removing normalization from the CNN with random L2 weights, following the approach of the previous subsection. We found that the local normalization had mild impact on the compatibility with the biological data. This can be seen in Fig. 8 (rand 2 vs rand 2 noN; in greedy, in full population, which was high and did not pass significance, and in subset regularized; independent sample -test).
A more direct way to tease apart the different computations (conv, ReLU, pool, norm) involved in L2, is to consider the intermediate outputs of the CNN trained on the ImageNet database as in Fig. 0(b), and to quantify the impact of each of these on the compatibility with the V2 data. This gives us a sense of how much each of the computations contribute to capturing the high-order statistics. We therefore generated outputs from each of the points in L2.
Outputs from conv2 (i.e., after only the convolution in the second layer) had high fitting errors. This is because the response from the conv layers can be negative before the ReLU. We found that compatibility to the V2 data starts to develop already after the rectification (i.e., ReLU2) stage (with Euclidean errors of subset greedy 0.33, full population 0.31, subset regularized 0.30). The fitting errors after pooling were (subset greedy 0.23; full population 0.20; subset regularized 0.24). After the local normalization (i.e., the point in L2 that we initially referred to in all our measurements), the fitting errors were (subset greedy 0.22; full population 0.19; subset regularized 0.23). The main improvement in the fit appeared to be at the L2 pooling stage.
As a third way to gauge the importance of local normalization, we retrained AlexNet on the ImageNet database, but with the local normalization layers removed. From an object recognition perspective, removing the normalization layers in the CNN model decreased the accuracy with a small margin (from 57.0% to 55.71%), echoing previous observations. Fitting errors with and without normalization were: subset greedy: 0.22 vs 0.30; full population: 0.19 vs 0.27; and subset regularized: 0.24 vs 0.30. This indicated that the trained model with local normalization resulted in a better fit than the model trained without normalization. Taken together, our results suggest that normalization had only a mild role in improving compatibility.
I-C5 Fitting higher layers of the CNN
So far, our main focus has been when selectivity to textures first develops, strongly differentiating L1 from L2. Nevertheless, it is interesting to consider how this changes across higher levels of the deep network. The results for Euclidean error are summarized in Table I. An increase in error might be expected at higher layers, similar to the observation that area V4 better discriminates synthesized “jumbled image” textures than area IT , . However, for higher layers of the CNN beyond L2, the error remained lower. Overall, the main change in error was between L1 and L2, compared to the smaller changes for the higher layers.
I-D Selectivity in other hierarchical model architectures
We have thus far focused on the AlexNet CNN, and obtained a good fit to V2 with L2 (but not L1) units. Our main approach was to take a single popular base model and its computational building blocks, and examine how manipulating certain pieces impacts the compatibility of the CNN model with the brain V2 data.
In this section, we consider several other popular models in the literature. We note that we are taking these models ”off the shelf” as is usually done in the literature, but that the various models differ from each other in many ways (e.g., the exact architecture and Linear Nonlinear operations, the receptive field size in each layer, and the number of units). Nevertheless, this exercise is informative for considering changes that develop across the first and second layers.
I-D1 Analysis of the VGG network layers
We have focused in detail on one particular popular deep CNN, AlexNet . We also applied our methodology, using the same input images as in AlexNet, to the VGG16  network (which we denote VGG), trained on the same ImageNet dataset. We found similar trends in the VGG layers as we have seen in AlexNet.
However, we found that the VGG develops selectivity more gradually than AlexNet. Starting from L1, the fit kept improving, with the largest reduction between L3 and L2 (similar to the L2 to L1 comparison in AlexNet). The Euclidean errors quantify this trend (Table I). In particular, between L2 and L3, the errors reduced from 0.42 to 0.15 in greedy, from 0.33 to 0.12 with the full population, and from 0.43 to 0.17 in the subset regularized technique. The fit for L4 remained similar to L3, with some increase in error for L5 with the subset selection methods. Since the VGG units have a smaller receptive size than AlexNet, we wondered if further downsampling the textures to obtain a more realistic RF size match between L2 and V2, would improve the result. However, downsampling further showed unreasonably high fitting errors, perhaps because in this case the images become too blurry and partially lose their texture properties.
In sum, for the VGG network, compatibility with the texture data first emerged in L3, and the change between L2 and L3 was the most striking. This means that L3 in VGG becomes more comparable with L2 in AlexNet (and to V2 in the brain). This is probably because the increase of the RF size in AlexNet is more rapid than in the VGG. L2 RF size in AlexNet is
(with stride = 2 in L1), but the same in the VGG net is only; in contrast, the RF size in L3 of the VGG net is . This is also consistent with the observation that the modulation index for the texture images computed over the Portilla and Simoncelli parameters  is weaker and more variable at small sizes (Corey Ziemba, personal communication; ). The VGG also does not include normalization, but based on our manipulations with AlexNet, we believe that its effect is minimal for the trained network in terms of compatibility with the V2 data. Overall, we found that middle layers in the VGG showed better compatibility with the biology.
I-D2 Analysis of HMAX and ScatNet
We have thus far focused on the CNN class of models. We wanted to further ask if differences develop across the first two layers for some related hierarchical visual models, such as HMAX ,  and ScatNet [15, 16], and how their selectivity compares to the CNN model.
HMAX is a popular model , , (see also ) motivated by the hierarchical organization of visual cortex. HMAX builds an increasingly complex and invariant feature representations in the upper layers and has been shown to be robust in recognizing both shape and texture-based objects. We presented the same texture stimulus ensembles and generated the response from both layers, named C1 and C2 (see Methods). We found that the HMAX model layer 2 fit the V2 data better than the HMAX layer 1. However, compared to the AlexNet CNN, the Euclidean errors for HMAX were considerably higher (greedy: 0.53 vs 0.22; full population: 0.44 vs 0.19; subset regularized: 0.52 vs 0.24) as shown in Fig. 8(a).
Can we further pinpoint what is different in the HMAX selectivity to the textures versus the noise? Our subset selection approach finds a population of units that can best fit the data. However, it doesn’t give an indication of the selectivity of the population on average to the textures versus the noise. We therefore considered as a complementary approach, randomly selecting 103 units 10000 times from the HMAX model, similar to what we did for AlexNet in our initial analysis. We found that the mean modulation index for the first two layers of HMAX respectively were (C1: -0.04 and C2: -0.01; see Fig. 8(b)). This negative mean modulation index means that on average, a random selection of 103 units in the HMAX was more selective to the noise than to the textures in the second layer. This was quite different from the V2 cortical data (e.g., -0.05 in the HMAX versus 0.12 in the V2 data). This is also in contrast to the CNN, which on average was more selective to the textures than the noise in L2, more similar to the V2 data (0.18 versus 0.12 in the V2 data). The CNN model was therefore more compatible to the cortical data.
, that computes invariant representations by wavelet transforms and pooling (see Methods). This network is able to discriminate textures by incorporating higher order moments. As with the other models, we found that the second layer of ScatNet fits the V2 data better than the first layer of ScatNet (Fig.9). However, in comparison, the CNN L2 units achieved better compatibility to the biological cortical data than the ScatNet layer 2 (Fig. 4(b) and 4(e)). The Euclidean error distance in ScatNet layer 2 was higher than in the CNN L2 (greedy: 0.39 vs 0.22; full population: 0.39 vs 0.19; subset regularized: 0.39 vs 0.24). We also found that the average modulation index of randomly selected 103 units was (M1 -0.14 and M2 0.05; see Fig. 8(b)). This indicates the reduced sensitivity of ScatNet units towards textures versus spectrally matched noise.
Our observation is that the CNN overall shows better correspondence with the V2 data compared to HMAX and ScatNet across the first two layers. Here the receptive field size is not a limiting factor as in the VGG, since the size develops more rapidly. It is important to realize that the CNN model differs from HMAX and ScatNet across many factors. First, these models differ in the approach to training. The CNN is trained on natural images in a supervised discriminative manner. The universal filters in the second layer of the HMAX model are obtained by extracting prototypes from a set of natural images. ScatNet uses an energy preservation property. One possibility is that the discriminative training of the CNN (as opposed to image prototyping or preservation) encourages faster development of sensitivity to the higher order statistics that distinguish the texture images from the noise. In addition, both the HMAX and ScatNet models we used consisted of two layers and the learning is not influenced by later layers, but AlexNet is deeper and has eight layers that are included in the supervised training; this may pose an advantage in the learning.
However, the amount of features learned is significantly outnumbered in the other models compared to the CNN. Whereas for the AlexNet CNN model we used 512 units ( neighborhood, 128 filters) in L2, the HMAX has 3200 units (8 patch sizes, 400 prototypes), and the ScatNet has 1536 units ( neighborhood, 384 filters). This may pose an advantage for the other models in selecting a population subset that best fits the data. There are also differences in the exact computational building blocks. However, all models include a pooling operation, with the max pooling in both HMAX and the CNN. All models also include a linear stage, with convolution in both the ScatNet and the CNN.
A number of related hierarchical models in the machine learning and neuroscience literature include similar basic computational building blocks stacked together, namely, convolution, spatial pooling, and sometimes local response normalization. This paper is motivated by the intriguing question: What makes CNNs, which are only very crudely matched to the brain architecture, able to capture some aspects of cortical processing surprisingly well?
We specifically focused on texture processing in early areas of visual cortex, for which  have compellingly shown develops in V2 but not in V1. This provided a rich data set to quantitatively compare the CNN and other related classes of hierarchical models. We found that L2 (but not L1) of AlexNet could well fit the V2 data. We also showed some qualitative similarities for texture selectivity between the first two layers of the CNN and the first two cortical visual areas.
We were interested in the question of when this compatibility first emerged, since in cortex there is a clear difference between the first two cortical neural areas for the texture stimuli. For a number of models we examined (AlexNet , VGG16 Net , HMAX , ScatNet ), we found that the details of the architecture and training made a difference regarding the compatibility of the model with the V2 data, and how quickly this developed across layers of the hierarchical network.
We presented initial versions of this work and the ability of deep neural networks to qualitatively capture some of the V2 versus V1 texture data in abstract form . Ziemba et al. have shown that a descriptive model of V2 can capture some of the qualitative results on the increased sensitivity to naturalistic textures in abstract and thesis form [47, 48]. Zhuang et al. showed increased sensitivity to textures versus noise in higher layers of deep neural networks, and related this to sparsity . The work described here, in contrast, asks at what point texture selectivity first emerges in the CNN layers, and focuses on both qualitatively comparing and quantitatively fitting the experimental data to changes that develop across the first two cortical areas.
What factors are important for better compatibility between the CNN and the cortical V2 data already in the second layer? We found that training on natural images was necessary for the model to develop compatibility with the cortical data. Various manipulations of random or shuffled weights could partly account for the modulation index data, but lead to a reduced fit between the CNN model and the V2 data relative to the learned weights. This indicated that the architecture by itself was not sufficient to obtain good compatibility. However, interestingly, having a trained first layer but random weights in the second layer, resulted in a good fit. The importance of retaining the trained weights (rather than random weights) in the first layer may be because the filters need to be matched to the frequencies and orientations that appear more often in the natural images in order to pick out the higher order structure of the textures in the subsequent layer.
In addition, we found that the receptive field size needed to develop sufficiently fast (comparing AlexNet and the VGG), and that the local normalization in Alexnet only had a limited role in obtaining good fits to the texture data. Both the HMAX and ScatNet followed a similar trend of better compatibility in the second versus the first layer, but overall did not fit the V2 data as well, and showed less sensitivity to the textures versus the noise in the second layer. This may be due to the discriminative learning in the CNN, versus the image prototyping or preservation of energy in the other models.
Although the AlexNet CNN showed compatibility with the V2 data, it too had some deviations from the cortical data. For instance, in the qualitative results in Fig. 2 it is intriguing to see that, with the same amount of neural units, the brain V2 outperforms the CNN L2 at grouping together different texture categories. In addition, as indicated by the modulation index, the CNN on average was more selective to the textures versus the noise than the V2 population. Its rank order of the selectivity to the textures on average also deviated from the data. This suggests that the CNN still has room for improvement in terms of capturing the cortical data. One possibility for improvement in future work is incorporating more realistic models of surround normalization (see, e.g., the range of surround normalization models used for modeling V1 data in ).
Local response normalization is a computation prevalent in visual cortex . It’s possible that the limited role of normalization in obtaining compatibility with the V2 data is due to the homogeneous nature of the textures. It is possible that divisive normalization will play a more important role in capturing data for non homogeneous images. Therefore, future work should test compatibility with V2 data over a broader range of natural stimuli and tasks. Future work should also incorporate other well known forms of cortical divisive normalization (e,g., surround normalization) into CNNs and consider their role in capturing cortical data across the hierarchy.
In addition, adding more units per layer may also play a role, and we found that choosing from a larger spatial neighborhood could improve the CNN fit. However, larger spatial neighborhoods did not reveal a good fit for the L1 units in the AlexNet, and the first layer of all models was not compatible with the V2 data. This also resonates with the original texture model of  that was actually used to generate the experimental stimuli; although the model does not have an explicit V2 unit representation, the textures are generated by joint statistics between V1 model units, i.e. by a 2-layer model. Though we have not exhausted all the possible hierarchical models, our method is pragmatic enough to be applied to any hierarchical models to test and find correspondence with the brain neurophysiology data.
While we have focused on texture selectivity in V2, there is room to explore compatibility of CNNs with changes across the early cortical hierarchy for other stimulus properties. For instance, biological studies have found selectivity in V2 to conjunctions of orientations, and to figure ground , , , with some aspects addressed in computational models of V2 , , , ,  . There may not be one unique CNN architecture that explains the neural data, but we believe that testing across fairly early visual areas (e.g., V1 and V2) with less stacked computations, and for a wider range of stimuli and tasks, can facilitate understanding of the critical factors and computations. Beyond area V2, studies have also shown that CNN units are compatible with shape tuning properties in visual area V4 .
In the quantitative comparisons between the modeling and data, we developed approaches for subset selection. These were more appropriate for the given cortical data than a weighted average of the neural population, which is typically used in fitting data. This is because the subset approach more faithfully represented the data analysis, which included an equally weighted average modulation index. This approach also allowed us to ask the question about whether there exists a population of units in the CNN that can well represent the experimental data. We therefore chose a neural population from the representation itself rather than a linear transform of the representation. By finding a subset of neural units that are most compatible with the data, it may be possible in the future to drive new experiments in which stimuli are generated from this population of model units and tested on the data. This may be applied more generally in the future to modeling other data sets and neural areas.
On one hand, our results add to the intriguing findings that CNNs trained on natural images have some compatibility with biological data, and moreover we found that this holds across low levels of the cortical hierarchy. But we believe that our approach goes beyond showing compatibility, by providing a direction for manipulating these early layers and teasing apart what aspects of CNNs (and other related hierarchical models), the training, and computational building blocks are most critical. This creates the opportunity for more discussion and systematic study of the various building blocks of deep networks, and opens the door to answer the long standing research question about correspondence between primate and machine vision.
We are grateful to Adam Kohn, Ruben Coen Cagli, and Corey Ziemba for discussions and very helpful comments on the manuscript; to Eero Simoncelli and Corey Ziemba for discussions and providing us with the original texture dataset used in the experiments; and to David Grossman and Ariel Lavi for discussions during their REU research. This work was supported by a Google Faculty Research Award, the National Science Foundation (grant 1715475), and a hardware donation from NVIDIA.
Iv-a Texture generation
and generated from an original set of 15 texture images. From each original texture, multiple synthetic texture images that matched the statistics of the original image were generated. Naturalistic textures for a given family were generated each with a different random seed, using an iterative process of constraining Gaussian white noise images to have similar marginal and joint statistical properties of the original textures
. Spectrally matched noise images were generated by randomizing the phase, i.e. computing the Fourier transform and inverse Fourier transform after phase randomization. From the 15 original textures, we have 15 different samples from each family, resulting in a total of 225 images of naturalistic textures, and 225 images of spectrally matched noise, as used in and . For cross-validation, we generated extra images per texture family from the model of .
Iv-B Matching receptive fields with the physiology data
We wanted to match as much as possible the spatial extent of the images that the model receptive field (versus the typical experimental neuron) is sensitive to. The input images were size , and the average receptive field size for V2 was approximately (with the V1 receptive field approximately half that size). The receptive field size of units in AlexNet is for L2 and for L1. We downsampled the input images by a factor of to obtain images of size , so that the effective size of the L2 receptive field was closer to the neurons recorded from in the experiment. We could not get an exact match, due to the constraint of downsampling by factors of 2. For HMAX, we downsampled only by a factor of to obtain images of size , since the C1 receptive fields in HMAX can be up to pixels. For ScatNet, the filter sizes are determined relative to the input images, so there was no need to downsample. We also originally ran the whole set of simulations without downsampling the images at all, and the results remained qualitatively similar, except that there were light improvements in the compatibility to the data with the appropriate downsampling.
Following the experiments, we contrast normalized the images before feeding them to the networks (CNN, ScatNet, Hmax). The luminance () is given by the mean pixel intensity of the downsampled image (). The contrast () is given by the standard deviation. The contrast normalized images () are then defined as follows:
where the desired contrast defines the range of the input pixels and the desired luminance defines the intensity centered around the range. We use 0.22 as the value of and since the desired luminance is gray, we use 0.5 as the value for .
Iv-C Deep CNN models for texture simulations
In our simulations, we mostly used the pre-trained AlexNet model, trained on natural images and specifically on the ILSVRC 2012 ImageNet  dataset. We also re-trained the network on ImageNet, or on the Places365 database , yielding similar results. We also contrasted this with an equivalent model architecture that included random weights (in the interval ) rather than pre-trained weights. AlexNet consists of five convolution followed by three fully connected layers. The first and second convolution layers are followed by local response normalizations and max-pooling. We used CaffeNet which is a variant of AlexNet where normalization follows the pooling. We refer to this as AlexNet for convenience. We examined the outputs from the first and second normalization layers (along with a more exhaustive examination of other layers), and compared them to the experimental data for V1 and V2 neuron outputs.
We used a modified version of AlexNet by changing the stride at L1 from 4 to 2. This allowed us to significantly reduce the receptive field size in L2 (from to ; with L1 of size ), making it more comparable with the biological ratio of V2 to V1 RF size. This modification also matched the experimental data better in our simulations. In addition, we simulated the response of the first 48 (instead of 96) L1 filters as they are the ones that show orientation selectivity; the remaining are more color selective. These 48 filters are the input to the first 128 (out of 256) channels in layer 2, so we considered these first 128 filters from L2. We focused mostly on the center four () spatial positions from each of these selected channels. We tested our method on larger neighborhoods and obtained qualitatively similar results.
We later considered the VGG16 network , trained on the same ImageNet dataset as AlexNet. VGG16 is a 16 layer network stacked with multiple (usually 2 or 3) convolution layers with filters and then followed by a pooling layer. We examined outputs from those five pooling layers, which we refer to as L1 through L5.
Iv-D Other hierarchical models for texture simulations
We also fit the texture data to the HMAX model . We generated responses from both layers, named C1 and C2. The C1 layer includes a range of receptive field sizes, with the largest approximately 49 by 49 pixels. The original HMAX model pools over the full extent of the image for L2; we reduced this to a 128 by 128 pixel region so that the receptive field size ratio between a C2 and C1 unit is more comparable to the receptive field size ratio between V2 and V1. To compute the C2 responses, we employed the universal set of units provided with the HMAX author’s code, available at . Furthermore, since the original code outputs the squared norm of the distance to the prototypes, we added the exponential function and selected a scaling factor (denoted as in ) of 1 to resemble the Gaussian like tuning of the neural responses.
Iv-E Local response normalization in CNN
Local response normalization plays an important role in hierarchical object recognition models. Local response normalization is used in both layer 1 and layer 2 of the AlexNet and has been shown to improve the recognition accuracy. It is loosely related to cross-orientation suppression in the brain, by normalizing the responses of groups of units with spatially overlapping receptive fields (5 in the case of AlexNet). If is the rectified linear activation at the position in each -th channel (or unit), then the normalized response is defined by 
where is the size of the normalization neighborhood, is the total number of channels in the layer. Constants are the hyper-parameters with the default values of 2, 5, , and 0.75 respectively.
Normalization is done across the spatially overlapping unit activations across channels. Each response is selected and normalized with corresponding values of all the channels across the channel dimension.
From a machine learning perspective, local response normalization is specifically useful to normalize the unbounded activations coming from the ReLU (Rectified Linear) non-linearity. It detects high-frequency features with large responses and penalizes the responses which are uniformly large in a local neighborhood. It is a type of regularization that encourages competition amongst units in the network.
Iv-F CNN population fitting and subset selection approaches
We considered the total number of units in the CNN as the number of units in a given layer (e.g., 48 for layer 1 and 128 for layer 2), times a center 2-by-2 spatial neighborhood. The rationale was that experimental data can be collected for receptive fields at different spatial positions. We chose a 2-by-2 spatial neighborhood, and did not find a significant difference when exploring larger spatial neighborhoods. We selected 103 units from L2 and 102 units from L1, to match the population numbers in the neurophysiology experiments of . Before starting the subset selection procedure, we removed from consideration the CNN neural units that had zero response to any family. This amounted to 432 units out of a total of 512 units from which we selected the subset of size 103.
Iv-F1 Subset greedy approach
We consider the subset greedy technique known as forward selection to choose a subset of 103 units that best match the data from the cortical neurons. In the greedy approach, the goal is to each time select another neural unit so as to minimize the Euclidean distance between the neural data and the CNN model modulation indices, until we have chosen 103 units from the available CNN layer population. The approach is greedy, because it optimizes the selection of the next unit as best it can given the current set of units. However, it does not guarantee a globally optimal solution.
Formally, let be a 15-dimensional vector containing the average modulation indices per texture family from the 103 recorded neurons in the physiological experiments, be the set of CNN neural units, and the average modulation indices per texture family of the th simulated unit in . Starting from , the greedy algorithm adds a unit to the current set of selected units that minimizes the squared Euclidean distance between the biological data and model average modulation indices:
The above procedure is repeated until the desired number of units is obtained. In particular, we stop at units.
Iv-F2 Optimal weighted average or full population
In this approach, we find a weighted average of the set of simulated units that best fit the biological data, by solving a constrained optimization problem. The constraint guarantees that the sum of the weights add to 1. This approach does obtain an optimal solution, and therefore shows the best one can do. However, note that it does not as faithfully capture the analysis of the neural data, since for the data analysis an equally weighted average of 103 neural responses give rise to the modulation index data.
Formally, we have:
where is the matrix of average modulation index values for all families computed according to (1), and is the vector of weights for all the simulated neurons. In terms of Euclidean distance, this is the best fit that could be attained by considering a weighted average on the simulated units. Nevertheless, note that the solution need not be sparse since there is no mechanism forcing the weights to become zero, and the disparity of the weight values can be hard to interpret.
Iv-F3 Subset regularized average followed by threshold
As noted above, for the optimal weighted average, weights can be very disparate. Since we seek to select a subset of the simulated units whose regular average (all weights are equal) follows closely the physiological experiments, we relax the selection problem by solving a regularized version of the optimization problem (5), as follows:
where is the trade-off parameter that promotes weight equalization. For , which is equivalent to solving (4), we found that only of the simulated neural units have the weights with only a handful of them containing large values that account for . As increases, the regularization term pushes the weights towards the center of the simplex. For instance, for , we found that approximately of the neural units have weights
. The subset of units is selected by applying a threshold to the estimated weights, as proposed in and then choosing the 103 units with the highest weights. However, a main difference from  is that our two-stage procedure is applied to the solution of (5) instead of (4). This approach also yields an excellent fit to the V2 data for L2 units, as shown in Fig. 5 (third row, middle column). We used for all fits; lower increased the fitting error but did not alter the trends (and vice versa).
For the cross-validation, we extended the image dataset. In the original data, there were only 15 images generated in each family. We therefore extended the set by generating 210 additional images (texture and noise) from each of the texture families. We optimized the learning, assuming that each group of 15 images, out of the 225 in each family, should yield an average modulation index that is as close as possible to the mean modulation index for that family in V2. Therefore, for 225 images in each family, we randomly divided the images into groups of 15. This yielded 15 data points per family, and a total of 225 data points for all 15 families. We then applied a 225 fold leave-one-out technique, training on 224 points and leaving out one point (corresponding to leaving out one set of 15 images). We thus learned the population (e.g., of 103 units in the greedy subset selection method) with the 224 training data points, and made a prediction of the mean modulation index for the left-one-out set of 15 images.
Iv-H Euclidean distance as the measurement of correspondence
To quantify the error between the CNN model and the neurophysiology data, we use the Euclidean distance metric. We calculate the Euclidean distance between the values computed from the CNNs and the neurophysiology V2 data. For the CNNs, the computed values include the modulation indices obtained from the various fitting and subset selection techniques, or predicted modulation index in the case of the cross validation.
Given two vectors and are two points in Euclidean -space, the Euclidean distance is computed using the 2-norm, as follows:
In all our experiments, is usually 15, the number of texture categories, as we take the average over samples and/or units. Lower Euclidean distances indicate a better fit of the model to the V2 data, and therefore higher correspondence of the model to the brain.
-  A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in Advances in neural information processing systems (NIPS), pp. 1097–1105, 2012.
-  Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” Nature, vol. 521, pp. 436–444, 2015.
-  Y. LeCun, B. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. Hubbard, and L. D. Jackel, “Handwritten digit recognition with a back-propagation network,” in Advances in neural information processing systems (NIPS), 1989.
-  Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” Proceedings of the IEEE, vol. 86, no. 11, pp. 2278–2324, 1998.
-  M. Zeiler and R. Fergus, “Visualizing and understanding convolutional networks,” in European Conference on Computer Vision (ECCV), 2014.
-  S.-M. Khaligh-Razavi and N. Kriegeskorte, “Deep supervised, but not unsupervised, models may explain IT cortical representation,” PLoS Comput Biol, vol. 10, no. 11, p. e1003915, 2014.
-  D. L. Yamins, H. Hong, C. F. Cadieu, E. A. Solomon, D. Seibert, and J. J. DiCarlo, “Performance-optimized hierarchical models predict neural responses in higher visual cortex,” Proceedings of the National Academy of Sciences, vol. 111, no. 23, pp. 8619–8624, 2014.
-  N. Kriegeskorte, “Deep neural networks: A new framework for modeling biological vision and brain information processing,” Annual Review of Vision Science, vol. 1, pp. 417–446, 2015.
-  D. L. Yamins and J. J. DiCarlo, “Using goal-driven deep learning models to understand sensory cortex,” Nature Neuroscience, vol. 19, no. 3, pp. 356–365, 2016.
-  R. M. Cichy, A. Khosla, D. Pantazis, A. Torralba, and A. Oliva, “Deep neural networks predict hierarchical spatio-temporal cortical dynamics of human visual object recognition reveals hierarchical correspondence,” arXiv preprint arXiv:1601.02970, 2016.
-  G. Umut and v. G. Marcel A. J., “Deep neural networks reveal a gradient in the complexity of neural representations across the ventral stream,” The Journal of Neuroscience, vol. 35, no. 27, pp. 10005–10014, 2015.
-  D. Pospisil, A. Pasupathy, and W. Bair, “Comparing the brain’s representation of shape to that of a deep convolutional neural network,” in The First International Workshop on Computational Models of the Visual Cortex: Hierarchies, Layers, Sparsity, Saliency and Attention, ACM, 5 2016.
-  C. F. Cadieu, H. Hong, D. L. K. Yamins, N. Pinto, D. Ardila, E. A. Solomon, N. J. Majaj, and J. J. DiCarlo, “Deep neural networks rival the representation of primate IT cortex for core visual object recognition,” PLoS Compt Biol., vol. 10, no. e1003963, 2014.
-  S. A. Cadena, G. H. Denfield, E. Y. Walker, L. A. Gatys, A. S. Tolias, M. Bethge, and A. S. Ecker, “Deep convolutional models improve predictions of macaque V1 responses to natural images,” bioRxiv, p. 201764, 2017.
-  J. Bruna and S. Mallat, “Invariant scattering convolution networks,” Transactions on Pattern Analysis and Machine Intelligence (PAMI), vol. 35, no. 8, pp. 1872–1886, 2013.
L. Sifre and S. Mallat, “Rotation, scaling and deformation invariant
scattering for texture discrimination,” in
Computer Vision and Pattern Recognition (CVPR), pp. 1233–1240, 2013.
-  T. Serre, L. Wolf, S. Bileschi, M. Riesenhuber, and T. Poggio, “Robust object recognition with cortex-like mechanisms,” IEEE Trans. on Pattern Analysis and Machine Intelligence, vol. 29, no. 3, pp. 411–426, 2007.
-  M. Riesenhuber and T. Poggio, “Hierarchical models of object recognition in cortex,” Nature Neuroscience, vol. 2, no. 11, pp. 1019–1025, 1999.
-  J. Freeman, C. M. Ziemba, D. J. Heeger, E. P. Simoncelli, and J. A. Movshon, “A functional and perceptual signature of the second visual area in primates,” Nature Neuroscience, vol. 16, no. 7, pp. 974–981, 2013.
-  C. M. Ziemba, J. Freeman, A. Movshon, and E. P. Simoncelli, “Selectivity and tolerance for visual texture in macaque V2,” Proc of Natl Academy of Science (PNAS), vol. 113, no. 22, 2016.
-  W. H. Merigan, “Cortical area V4 is critical for certain texture discrimination, but this effect is not dependent on attention,” Visual Neuroscience, vol. 31, no. 23, pp. 8543–8555, 2000.
-  A. Hanazawa and H. Komatsu, “Influence of the direction of elemental luminance gradients on the responses of V4 cells to textured surfaces,” Journal Neuroscience, vol. 21, pp. 4490–4497, 2001.
-  F. Arcizet, C. Jouffrais, and P. Girard, “Natural textures classification in area V4 of the macaque monkey,” Exp Brain Research, vol. 189, pp. 109–120, 2008.
-  A. S. Nandy, T. O. Sharpee, J. H. Reynolds, and J. F. Mitchell, “The fine structure of shape tuning in area V4,” Neuron, vol. 78, pp. 10780–10793, 2013.
-  G. Okazawa, S. Tajima, and H. Komatsu, “Image statistics underlying natural texture selectivity of neurons in macaque V4,” Proceedigs of National Academy of Science (PNAS), vol. 112, no. 4, 2015.
-  G. Okazawa, S. Tajima, and H. Komatsu, “Gradual development of visual texture-selective properties between macaque areas V2 and V4,” Cerebral Cortex, vol. 27, no. 10, pp. 4867–4880, 2016.
-  B. Julesz, “Experiments in the visual perception of texture,” Scientific American, vol. 232, no. 4, pp. 34–43, 1975.
-  H. Tamura, S. Mori, and T. Yamawaki, “Textural features corresponding to visual perception,” IEEE Transactions on Systems, Man, and Cybernetics, vol. 8, no. 6, pp. 460–473, 1978.
-  B. Julesz, “Textons, the elements of texture perception and their interactions,” Nature, vol. 290, pp. 91–97, 1981.
-  J. R. Bergen and E. H. Adelson, “Early vision and texture perception,” Nature, vol. 333, no. 6171, pp. 363–364, 1988.
-  J. Malik and P. Perona, “Preattentive texture discrimination with early vision mechanisms,” J. of the Optical Society of America. A, Optics and image science, vol. 7, no. 5, pp. 923–932, 1990.
-  T. S. A. Wallis, C. M. Funke, A. S. Ecker, L. A. Gatys, F. A. Wichmann, and M. Bethge, “A parametric texture model based on deep convolutional features closely matches texture appearance for humans,” Journal of Vision, vol. 17-5, no. 12, 2017.
-  D. J. Heeger, “Normalization of cell responses in cat striate cortex,” Visual neuroscience, vol. 9, no. 02, pp. 181–197, 1992.
-  M. Carandini and D. J. Heeger, “Normalization as a canonical neural computation,” Nature Reviews Neuroscience, vol. 13, no. 1, pp. 51–62, 2012.
-  O. Russakovsky, J. Deng, H. Su, J. Krause, S. Satheesh, S. Ma, Z. Huang, A. Karpathy, A. Khosla, M. Bernstein, A. C. Berg, and L. Fei-Fei, “ImageNet Large Scale Visual Recognition Challenge,” International Journal of Computer Vision (IJCV), vol. 115, no. 3, pp. 211–252, 2015.
-  J. Portilla and E. P. Simoncelli, “A parametric texture model based on joint statistics of complex wavelet coefficients,” International Journal of Computer Vision (IJCV), vol. 40, no. 1, pp. 49–71, 2000.
-  L. van der Maaten, “Accelerating t-SNE using tree-based algorithms,” Journal of Machine Learning Research (JMLR), vol. 15, pp. 3221–3245, 2014.
-  L. van der Maaten and G. E. Hinton, “Visualizing data using t-SNE,” Journal of Machine Learning Research (JMLR), vol. 9, pp. 2579–2605, 2008.
-  B. Zhou, A. Lapedriza, A. Khosla, A. Oliva, and A. Torralba, “Places: A 10 million image database for scene recognition,” IEEE Transactions on Pattern Analysis and Machine Intelligence, 2017.
-  K. Jarrett, K. Kavukcuoglu, M. Ranzato, and Y. LeCun, “What is the best multi-stage architecture for object recognition?,” in Computer Vision, 2009 IEEE 12th International Conference on, pp. 2146–2153, IEEE, 2009.
-  A. Saxe, P. W. Koh, Z. Chen, M. Bhand, B. Suresh, and A. Y. Ng, “On random weights and unsupervised feature learning,” in Proceedings of the 28th international conference on machine learning (ICML-11), pp. 1089–1096, 2011.
-  N. C. Rust and J. J. DiCarlo, “Selectivity and tolerance (”invariance”) both increase as visual information propagates from cortical area V4 to IT,” Journal of Neuroscience, vol. 30, no. 39, pp. 12978–12995, 2010.
-  K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” in International Conference on Learning Representations (ICLR), 2015.
-  C. Ziemba, R. Perez, E. Simoncelli, and J. Movshon, “Selectivity of contextual modulation in macaque V1 and V2,” in Annual Meeting, Neuroscience, Nov 2017.
-  K. Fukushima, “Neocognitron: A self-organizing neural network model for a mechanism of pattern recognition unaffected by shift in position,” Biological cybernetics, vol. 36, no. 4, pp. 193–202, 1980.
-  M. N. U. Laskar, L. G. S. Giraldo, and O. Schwartz, “Deep learning captures V2 selectivity for natural textures,” in Computational and Systems Neuroscience (Cosyne), Abstract, 2017.
-  C. M. Ziemba, R. L. Goris, J. A. Movshon, and E. P. Simoncelli, “Characterizing receptive field selectivity in area v2,” in MODVIS Workshop, VSS abstract, 2015.
-  C. M. Ziemba, Neural representation and perception of naturalistic image structure. PhD thesis, Center for Neural Science, New York University, New York, NY, May 2016.
-  C. Zhuang, Y. Wang, D. Yamins, and X. Hu, “Deep learning predicts correlation between a functional signature of higher visual areas and sparse firing of neurons,” Front. Comput. Neuroscience, vol. 11, no. 100, 2017.
-  R. Coen-Cagli, A. Kohn, and O. Schwartz, “Flexible gating of contextual modulation during natural vision,” Nature Neuroscience, vol. 18, pp. 1648–1655, 2015.
-  M. Ito and H. Komatsu, “Representation of angles embedded within contour stimuli in area V2 of macaque monkeys,” The Journal of Neuroscience, vol. 24, pp. 3313–3324, 2004.
-  Y. El-Shamayleh and J. A. Movshon, “Neuronal responses to texture-defined form in macaque visual area V2,” The Journal of Neuroscience, vol. 31, no. 23, pp. 8543–8555, 2011.
-  H. Zhou, H. S. Friedman, and R. Von Der Heydt, “Coding of border ownership in monkey visual cortex,” The Journal of Neuroscience, vol. 20, no. 17, pp. 6594–6611, 2000.
-  D. J. Heeger, E. P. Simoncelli, and J. A. Movshon, “Computational models of cortical visual processing,” Proceedigs of National Academy of Science (PNAS), vol. 93, no. 4, pp. 623–627, 1996.
-  H. Lee, C. Ekanadham, and A. Ng, “Sparse deep belief net model for visual area V2,” in Advances in neural information processing systems, pp. 873–880, 2008.
-  R. Coen-Cagli and O. Schwartz, “The impact on midlevel vision of statistically optimal divisive normalization in V1,” Journal of Vision, vol. 13, no. 8, 2013.
-  H. Hosoya and A. Hyvarinen, “A hierarchical statistical model of natural images. explains tuning properties in V2,” The Journal of Neuroscience, vol. 35, no. 29, pp. 10412–10428, 2015.
-  L. Zhaoping, “Border ownership from intracortical interactions in visual area V2,” Neuron, vol. 47, no. 1, pp. 143–153, 2005.
-  J. Hegde and D. C. V. Essen, “Selectivity for complex shapes in primate visual area V2,” Journal of Neuroscience, vol. 20, no. RC61, 2000.
-  http://maxlab.neuro.georgetown.edu/hmax.html#code, 2007. [Online; accessed 23-April-2018].
-  http://www.di.ens.fr/data/software/scatnet/download/, 2013. [Online; accessed 23-April-2018].
-  P. Li, S. Sundar Rangapuram, and M. Slawski, “Methods for Sparse and Low-Rank Recovery under Simplex Constraints,” ArXiv e-prints, May 2016.