Cerebral cortex is a complex dynamical system dominated by feedback circuits, but we limit our exploration to the feed-forward component of this system, which endows neurons with their identity-defining so-called “classical” RFs and feature-tuning properties. We further confine our exploration to the central pathway in the feed-forward elaboration of cortical neurons’ properties, which proceeds through a repeating sequence of two cortical layers. Cortical layer 4 (L4) is the principal initial recipient of the feed-forward afferent input to a cortical area. L4 converts that input into a new form and sends it, in particular, to layer 3 (L3) of the same cortical area for its feature-extraction operation. The product of that L3 operation is then sent to L4 of the next cortical area, where the same two-stage feature-extracting operation is repeated, but on a higher level, building on the advances made by the preceding cortical area (Rockland and Pandya, 1979; Felleman and Van Essen, 1991; Callaway, 2004).

In addition to their afferent input from L4, L3 neurons receive extensive contextual input via long-range horizontal connections from surrounding columns up to several millimeters away from their resident column (Gilbert and Wiesel, 1983; DeFelipe et al., 1986; Lund et al., 1993; Burton and Fabri, 1995). These contextual inputs are expected to guide L3 neurons to the sources of mutual information in these two, afferent and contextual, input sets (Phillips and Singer, 1997). Two principally different kinds of such sources are possible. First, the two input sets might have partially overlapping RFs; in other words, they might share some neurons along their afferent pathways in common. Such an internal source of mutual information in the two sets is trivial and has to be avoided. Second, the afferent and contextual input sets might be impacted by the same environmental agent. Such an external source of mutual information has the potential of being behaviorally significant and therefore worthy of recognition. One way to ensure that mutual information in the two input sets comes from external sources is to use inputs with non-overlapping RFs. Indeed, long-range horizontal connections come from far enough to have non-overlapping RFs but are close enough to reflect the same distal variables in the engaged environment.

The value of the feature (φ) extracted by the i^th L3 cell is computed from the afferent inputs to its resident cortical column m, as in Equation 1:

where f mi is the chosen feature-specific transfer function, and I m → is a vector of activities of all the afferent axons innervating column m. According to the contextual guidance proposal, the function f mi is chosen so as to maximize correlation of φ mi with the “best” function g mi over features extracted in other, surrounding columns and delivered via long-range horizontal connections, as in Equation 2:

where corr is Pearson’s correlation coefficient, and φ → context is a vector of all available L3 contextual features combined.

Since different neurons in a column should extract different features, we should elaborate the choice of the feature-extracting transfer function f mi for a single neuron i in column m: the choice is to maximize correlation of φ mi with the “best” function g mi over features in other columns, subject to the constraint that correlation of this feature with features computed by other neurons in the same layer in the same column should not be excessive.

It should be noted that features extracted by cortical neurons – especially in high-level cortical areas – are highly nonlinear functions of peripheral patterns of receptor activations. No recursive application of linear transform functions would be able on its own to produce such features. Cortical neurons must be able to use nonlinear transfer functions. However, experience-driven learning of nonlinear transfer functions in neural networks can be highly problematic, unless a kernel-based strategy is used. Kernel methods, popular in machine learning, offer a highly effective strategy for dealing with nonlinear problems by transforming the input space into a new “feature” space where a nonlinear problem becomes linear and thus more tractable with efficient linear techniques (Schölkopf and Smola, 2002). According to Favorov and Kursun (2011), such a kernel-based function linearization strategy happens to be used in the neocortex in its principal input layer, layer 4. This insight suggests that cortical columns first perform a nonlinear function-linearization transform of their afferent inputs in L4 and then learn linear transform functions in L3.

An important feature of L4 functional architecture is the presence of untuned feed-forward inhibition, which reflects the overall strength of the stimulus activating a local L4 network but is insensitive – invariant – to spatial details of the stimulus patterns (Kyriazi et al., 1996; Bruno and Simons, 2002; Swadlow, 2003; Hirsch et al., 2003; Sun et al., 2006; Cruikshank et al., 2007). Favorov and Kursun (2011) showed that the presence of such untuned feed-forward inhibition converts a conventional neural network into a functional analog of Radial Basis Function (RBF) networks (Lowe, 2003), which are well known for their universal function approximation and linearization capabilities (Park and Sandberg, 1991; Kurková, 2003). Input transforms performed by such networks automatically linearize a broad repertoire of nonlinear functions over the afferent inputs. This capacity for pluripotent function linearization suggests that L4 can contribute importantly to cortical feature extraction by performing such a transform of afferent inputs to a cortical column that makes possible for neurons in the other layers of the column, including L3, to extract nonlinear features of afferent inputs using mostly linear operations.

A biologically realistic and highly effective pluripotent function linearizer has the following ingredients (Equation 3): (1) activity of each excitatory L4 cell is computed, in part, as a weighted sum of its afferent inputs, which are Hebbian; (2) lateral interconnections among L4 cells are used to diversify the afferent connectional patterns among L4 cells in a cortical column and give them a rich variety of RF properties; and (3) feed-forward inhibition makes L4 cells behave similarly to RBF units and is principally responsible for function linearization capabilities. Following Favorov and Kursun (2011), we describe L4 operation as:

where F L 4 i is the activity of L4 neuron i; a_j is the activity of afferent input neuron j; w_i,j is the weight, or efficacy, of the excitatory synaptic connection from afferent neuron j to L4 neuron i; λ is a lateral connection scaling constant; F L 4 k is the activity of a neighboring L4 neuron k; ρ i , k is the correlation coefficient between activities of L4 neurons i and k; θ is a feed-forward inhibition scaling constant; and [·]⁺ indicates that if the quantity in the brackets is negative, the value is to be taken as zero.

In the presence of both feed-forward inhibition and plastic lateral connections, which are unique to L4 in that higher correlation in firings of the pre- and post-synaptic cells leads to decrease – rather than increase – of synaptic strength (Egger et al., 1999; Sáez and Friedlander, 2009), a modeled network of L4 neurons, trained on visual inputs, develops biologically accurate diversity of multi-subfield RFs and acquires orientation tuning matching in sharpness that of real L4 neurons, as well as a host of other real L4 functional properties (Favorov and Kursun, 2011).

In the above L4 model, neurons accomplish their pluripotent function linearization by acting in local groups. To explain, consider a local group of L4 neurons that are innervated by a common set of afferent neurons. Together, such a set of N convergent afferent neurons can be viewed as defining an abstract N-dimensional afferent input state space, each dimension corresponding to one of the constituent afferent neurons. A targeted L4 neuron i takes a particular direction in this afferent space (defined by the vector of its afferent connection weights w → i in Equation 3) as its conic RBF center. Neighboring L4 neurons innervated by the same set of afferent neurons chose different RBF centers (i.e., different w → ) in their common afferent space, influenced in their choices by lateral interactions with each other. Together, a local group of L4 neurons will spread their RBF centers evenly throughout their common afferent space so as to map it most efficiently, with more active regions of that space mapped at higher resolution (Deco and Obradovic, 1995). A wide range of nonlinear functions defined over this space – including transfer functions that can extract contextually predictable features – can be then approximated by weighted sums of the activities of the mapping L4 neurons.

What is the size of such RBF-network-like local groups of L4 neurons that together can perform pluripotent function linearizations? Although it is limited, anatomical evidence based on typical sizes of afferent axon arborizations in L4 and the lateral spread of dendrites and axon collaterals of L4 cells within the confines of L4 suggests that such groups should be larger than a minicolumn (i.e., >50 μm in diameter) but smaller than a macrocolumn (i.e., <500 μm). We will refer to such intermediate-size function-linearizing columns in this paper as mesocolumns.

Structurally, minicolumns are the radially oriented cords of neuronal cell bodies evident in Nissl-stained sections of the cerebral cortex (Mountcastle, 1978; Buxhoeveden and Casanova, 2002). They are the narrowest (~50 μm diameter) columnar aggregates of neurons in the neocortex (Favorov and Diamond, 1990; Tommerdahl et al., 1993) and thus can be viewed as the smallest building blocks of cortical columnar organization (Mountcastle, 1978). Published estimates of L4 cell densities in visual and somatosensory cortical areas (Beaulieu and Colonnier, 1983; Budd, 2000; Markram et al., 2015; Meyer et al., 2010) suggest that a single minicolumn has between 30 and 60 excitatory L4 neurons. Such a number is clearly not enough for a functionally useful RBF mapping of a minicolumn’s afferent space.

Minicolumns are packed together in the cortex in an essentially hexagonal pattern. From a geometric perspective, the next larger-size columnar entity to consider is a group of 7 minicolumns, one surrounded by 6 others. Such columns will be 3 minicolumns wide and thus ~150 μm in diameter. They will have between 200 and 400 excitatory L4 neurons. In Section 3 we will show that such numbers of L4 neurons are sufficient for the purposes of contextually predictable feature extraction in L3. We propose that, based on the available evidence, such groupings of 7 minicolumns are the most plausible candidates for the role of function-linearizing mesocolumns (Figure 1). A group of 7 such mesocolumns, in turn, make the next larger-size columnar entity ~450 μm in diameter, corresponding in size to the well-known macrocolumns (Mountcastle, 1997).

The lateral spread of axonal projections of L4 cells in L3 and the lateral spread of basal dendrites of pyramidal cells in L3 indicate that L4 neurons send their output to L3 neurons not only in their own mesocolumn but also in the 6 surrounding mesocolumns (Lubke et al., 2003; da Costa and Martin, 2010). Consequently, the “classical” RFs and feature-tuning properties of L3 neurons in a given mesocolumn are, essentially, the product of weighted summation of output activities of L4 neurons of the same and 6 surrounding mesocolumns (Figure 1):

where φ i is the feature-expressing activity of L3 neuron i in the central mesocolumn; N_L4mc is the number of L4 cells in a mesocolumn; F L 4 mj is the activity of an L4 neuron j in mesocolumn m; and u_i,mj is the strength of their connection.

In L3, the afferent inputs from L4 target basal dendrites of pyramidal cells, whereas the contextual inputs from surrounding cortical territories target their apical dendrites (Gilbert and Wiesel, 1983; Kisvardy et al., 1986; Lubke et al., 2003; Petreanu et al., 2009). Synaptic inputs to basal dendrites are integrated in the soma, leading to spike generation in the initial axon segment. But the apical dendrite, including its terminal tuft extension in layer 1, has its own site of synaptic input integration and is able to generate its own spikes (Bernander et al., 1994; Cauller and Connors, 1994; Schiller et al., 1997; Stuart and Spruston, 1998; Larkum et al., 1999, 2007). Output activity of the apical dendrite in the i^th L3 cell is, essentially, the product of weighted summation of output activities of L3 neurons of surrounding columns:

where C is the set of all the L3 neurons in surrounding columns that contribute contextual input to cell i; v i , c is the strength of connection to cell i from contextual cell c; and φ c is basal dendrite output of cell c (Equation 4). With such separate contextual input integration, the apical dendrite can guide basal dendrites in their selection of afferent connectional patterns (and vice versa) so that they will maximize covariance of the cell’s apical and basal outputs A_i and φ i , as was proposed and successfully demonstrated in a basic model by Kording and Konig (2000).

The number of excitatory cells in L3 is approximately the same as in L4 (Markram et al., 2015; Meyer et al., 2010), suggesting that the L3 compartment of a mesocolumn contains between 200 and 400 pyramidal neurons. Under mutual competitive pressure to diversify their RF tuning properties, similar to plastic local lateral connections among neighboring L4 cells driving them to select different features, these 200–400 neurons in a mesocolumn will compete in their search for contextually predictable features. Together, they will find and tune to all the different contextually predictable features present in their shared afferent input from L4 (Equation 4).

In geometric terms, together the 1,400–2,800 L4 cells in 7 mesocolumns that provide afferent input to the L3 compartment of a central mesocolumn create that mesocolumn’s high-dimensional state space. Since we are concerned with L3 features that are computed linearly in that state space (Equation 4), such features correspond to particular directions in the mesocolumn’s L3 state space: for a given L3 neuron i, its afferent connectional vector u → i determines that neuron’s preferred direction in its mesocolumn’s L3 state space and thus its preferred feature.

Any arbitrary direction in the L3 state space will express some feature of input patterns taking place in the mesocolumn’s RF. However, most of such arbitrarily chosen features will not be perceptually significant. Contextually predictable features occupy a lower-dimensional subspace of the L3 state space. The size and contents of this contextually predictable subspace have never been revealed before, even for the primary sensory cortical areas. In the next section (Section 2.4), we derive a computational algorithm for estimating the principal axes (basis vectors) of this subspace, and in Section 3 we apply this algorithm to natural images in an attempt to reveal the contextually predictable subspace of a mesocolumn in the primary visual area, V1. Also in Section 3 we investigate what features (i.e., directions in the mesocolumn’s L3 state space) individual L3 cells will choose if they are modeled as comprising two dendritic compartments – basal and apical, each receiving Hebbian connections from either L4 of its own macrocolumn or from L3 of surrounding macrocolumns, respectively – and are trained on natural images. We show that such modeled L3 cells do indeed select features in the mesocolumn’s contextually predictable subspace.

We begin by formalizing terminology to be used in the rest of this paper:

We formulate our approach based on the following considerations. We define the 1st axis of the canonical feature subspace (i.e., the 1st canonical variate) to be the basis vector with the maximal correlation with the contextual input, the second axis (i.e., the 2nd canonical variate) to be the basis vector with the second largest correlation with the contextual input, and so on until the last axis. In the cortex, different mesocolumns develop their own sets of afferent, lateral, and contextual connections based on their particular histories of sensory experiences. However, since neighboring mesocolumns will end up being exposed to and being shaped by the same regularities in their sensory experiences, any emergent differences among them will not be functionally significant. Thus, in deriving our algorithm, we can make an assumption that all the neighboring mesocolumns involved in contextual guidance will have the same matrices of L3 afferent [u_i,mj] and contextual [ v i , c ] connections (Equations 4, 5) and also identical sets of canonical variates.

To quantify the correlation of a canonical variate with the contextual input, we use the mean of pairwise correlations of that variate in the central mesocolumn and the same variate in each of the neighboring mesocolumns that contribute the contextual input. If we label the direction of a canonical variate in the L3 state space as b → (b stands for “basis”), then we define the contextual correlation of this variate as in Equation 6:

where N_m is the number of mesocolumns contributing contextual input. F → af f 0 and F → aff m are the afferent inputs to the L3 compartment of the central (0 ^th ) mesocolumn and the m^th neighboring mesocolumn, respectively, from their flattened 7 × N L 4 mc dimensional vectors of the outputs of L4 neurons of the same and 6 immediately surrounding mesocolumns:

Our objective function for the i^th canonical variate is to find such a direction b → i in the L3 state space that will maximize its contextual correlation r i (subject to the constraint that b → i ⊥ b → j for all j < i ) . Our objective function is designed to maximize the canonical correlation of the i^th component, concurrently ensuring orthogonality with all previously computed components. This methodology, which constructs orthogonal vectors sequentially, beginning with the first, systematically generates a series of orthogonal vectors. Each vector maximizes the variance subject to the orthogonality constraints imposed by its predecessors.

Continuing with our assumption that different mesocolumns in a contextually related cortical territory have the same internal connectivities, we also assume that all N_m mesocolumns in our model have the same means and covariance matrices of the afferent inputs to their L3 compartments:

Then, we can write our objective function as:

Thus, the objective function in Lagrangian formulation is given by:

The optimization task specified in Equation 10 has a well-structured generalized eigenproblem and can be efficiently solved using established numerical algorithms (Hotelling, 1936; Hardoon et al., 2004; Kursun et al., 2011; Golub and Van Loan, 2013; Alpaydin, 2014), in which the eigenvector having the largest eigenvalue giving us the first canonical variate b → 1 , the eigenvector having the second largest eigenvalue giving us the second variate b → 2 , and so on:

where K cross is the cross-covariance matrix (Hotelling, 1936; Kursun et al., 2011):

According to our proposed division of tasks between L4 and L3 – with the L4 mesocolumnar network linearizing feature-extracting functions that will be computed by L3 cells – the first step in estimating L3 canonical variates is to develop RF and functional properties of cells in the L4 compartment of mesocolumns. This is done by repeatedly exposing L4 cells to images and adjusting the weights of their Hebbian input and intrinsic connections, gradually driving them into stable connectional patterns. The emergent functional RF properties of the model L4 cells, which come to closely resemble those of simple cells in cat V1, are comprehensively described in Favorov and Kursun (2011), and for brevity we omit their description here.

On their own, the trained L4 cells have very low pairwise contextual correlations with L4 cells in surrounding macrocolumns. This is shown in Figure 4A by plotting the distribution of maximal correlations of coincident activities of individual L4 cells in the central macrocolumn and L4 cells in the first and second rings of surrounding macrocolumns. However, our expectation is that optimally chosen weighted sums of multiple L4 cells will have much higher contextual correlations with surrounding macrocolumns. In Section 2.4 we introduced a particular algorithm for finding such optimal weighted sums, giving us the axes of the canonical feature subspace, i.e., canonical variates. We apply this algorithm to the outputs of L4 cells in the modeled field of 19 macrocolumns to obtain canonical variates. We test the strength of their contextual correlations by introducing a third ring of 12 mesocolumns, chosen to be at a such distance from the central mesocolumn that any mutual information they might have in their RFs will have to come from the environmental sources rather than from sharing any pixels in common. To compute their contextual correlations, we used responses Φ (Equation 20) of the canonical variates in the central and these 12 distant surrounding mesocolumns to 1,000 image patches taken at random in the 100 dataset images. For each canonical variate, its contextual correlation is expressed by Pearson correlation coefficient computed between the 1,000 responses of that variate in the central mesocolumn and the 1,000 averages of responses of that variate in the 12 distant mesocolumns. The magnitudes of the computed contextual correlations are plotted in Figure 4B (white bars) for the first 20 canonical variates. The first 7 variates have particularly high correlations (r² ≥ 0.1). Correlations of the 8th to 15th variates, although low, are nevertheless statistically significant (at α = 0.05 with Bonferroni correction), suggesting that even these canonical variates might reflect some causally significant factors in the environment.

To demonstrate the necessity of the L4 function-linearization operation for maximizing contextual correlations, we also developed canonical variates directly from the LGN afferent inputs to the L4 compartments of the central and 6 surrounding mesocolumns (together constituting a macrocolumn) rather than using L4 outputs of these 7 mesocolumns. Unlike the L4-based variates, all but the first of the LGN-based variates showed no statistically significant contextual correlations (black bars in Figure 4B).

To test generalizable nature of the canonical variates extracted from the 100 IARP TC-12 images, canonical variates were also extracted from the 5,000 COCO images. The magnitudes of their contextual correlations, shown for the first 40 variates, are plotted as black bars in Figure 4C, superimposed on the first 40 canonical variates extracted from IARP images (plotted as white bars). As Figure 4C shows, although they come from different sources, magnitudes of the two sets of canonical variates are very similar, with the first 15 variates having statistically significant contextual correlations. But how similar are the features extracted from the two image sources? As basis vectors, the first 15 contextually predictable canonical variates enclose the canonical feature subspace of the mesocolumn’s entire feature space (as defined in Section 2.4). To determine how much the IARP and COCO canonical feature subspaces overlap, we performed Canonical Correlation Analysis (CCA; Hotelling, 1936), in which we treated the first 15 IARP and first 15 COCO canonical variates as 2 sets of input variables and used 5,000 training image patches, taken at random from the COCO dataset, to compute their loadings. Next, we used these loadings to compute 15 canonical correlations of the two sets of variables over a different set of randomly picked 1,000 COCO image patches. If the two feature subspaces, enclosed by the 15 IARP and 15 COCO canonical variates, match closely, the 15 canonical correlations would all be close to 1. On the other hand, if the two subspaces do not overlap at all, the 15 canonical correlations would all be close to zero. The actual computed correlations are plotted in Figure 4D, revealing that the first 6 CCA variates had very high correlations, whereas the last 5 CCA variates had very low correlations. Thus, we can conclude that canonical feature subspaces extracted from the IARP and COCO datasets mostly overlap, albeit not completely.

Going back to Figure 4B, as it shows, only the first canonical variate does not depend on L4 function-linearization preprocessing. The reason is that it reflects the overall magnitude of activity evoked in the macrocolumn’s L4 compartment (Figure 5A) and thus the overall stimulation intensity of the macrocolumn’s RF. Since the other canonical variates depend on L4 function-linearization preprocessing, they must be tuned to various structural features of the image patterns occurring in the mesocolumn’s RF. What these features are, either in our model canonical variates or in real L3 neurons, is not obvious but some insight is traditionally gained in V1 studies by characterizing responses of V1 neurons to moving grating images of various orientations and spatial frequencies. Figure 5B shows orientation and positional tuning of the statistically significant first 15 canonical variates, revealing that variates 8 through 11 are sensitive to both orientation and position while others are sensitive to grating orientation but not its position in the RF (translational invariance), thus falling into the categories of the simple and complex cells, respectively (Hubel and Wiesel, 1962). With real V1 neurons exhibiting diversity in the degrees of their orientation and grating phase tuning, the standard metric used to place any given V1 cell on the simple vs. complex cell spectrum is the F1/F0 ratio, which is the ratio of the 1st and 0th Fourier harmonics of a neuron’s activity during stimulation of its RF with an optimal sinewave moving grating (Skottun et al., 1991; Ringach et al., 2002). V1 cells with F1/F0 > 1 are classified as simple and cells with F1/F0 < 1 are classified as complex. Figure 5C shows F1/F0 scores of the statistically significant canonical variates 2–15, showing that 70 and 30% of variates fall into the complex cell and simple cell categories, respectively.

In principle, function-linearization capabilities of mesocolumns’ L4 compartment depend on the number of cells they employ (Favorov and Kursun, 2011): the larger the number of L4 cells in a mesocolumn, the broader the repertoire of nonlinear functions it can linearize. This is shown in Figure 6, in which the total contextual correlation of the first 20 canonical variates, computed as the sum of squared contextual correlations of individual variates (Watanabe, 1960), is plotted as a function of the number of cells in each mesocolumn’s L4 compartment. Significantly, there is little further gain in total correlation after the number of L4 cells in mesocolumns is increased beyond 150–200, which suggests that they linearize all the contextually predictable features available for extraction in the mesocolumn’s RF.

L3 cells are expected to be driven by their apical dendrites to tune to contextually predictable – canonical, according to our terminology – features. Such features occupy a particular subspace in the mesocolumn’s state/feature space, and the extracted canonical variates give us the principal axes of this canonical feature subspace. In choosing their features, L3 cells should be attracted to the canonical variates according to their contextual predictability. We explored this feature-selecting mechanism in our extended field of 19 macrocolumns (Figure 3A) by providing the L3 compartment of the central mesocolumn (red-shaded in Figure 3A) with 150 cells, each modeled as a pair of basal and apical dendrites. In each L3 cell, its basal dendrite was given afferent input from 1,050 L4 cells of its own and 6 immediately adjacent mesocolumns (150 L4 cells per mesocolumn), which together make up a macrocolumn (blue-shaded in Figure 3A). The vector of the weights of these L4 connections to the cell’s basal dendrite ( u → in Equation 4) determines that dendrite’s preferred direction in the mesocolumn’s L3 state space and thus that cell’s preferred feature. The apical dendrite of each L3 cell was given contextual input from the 18 more distant mesocolumns (black-shaded in Figure 3A) in the form of their 15 statistically significant canonical variates (for further details, see Methods section 3.1.6). The apical dendrite learns to produce output that best matches the output of the basal dendrite and vice versa.

Feature selectivities of L3 cells were developed by repeatedly exposing the model to randomly picked dataset images and adjusting the weights of the Hebbian basal and apical connections of L3 cells, gradually driving them into stable connectional patterns (Figure 7). After such training, the apical and basal dendrites of L3 cells developed prominent correlations in their responses to image patches (Figure 8A), demonstrating that all 150 L3 cells succeeded in tuning to contextually predictable features. Furthermore, Figure 8B shows that cross-correlations between basal outputs of different L3 cells residing in the same mesocolumn are low, indicating that these cells tuned to diverse canonical features.

Figure 8C provides some insight into the nature of the features chosen by the 150 L3 cells. It plots correlations of the basal dendrite of each L3 cell in the central mesocolumn with each of the first 15 canonical variates computed for that mesocolumn. The plot shows that each L3 cell developed either positive or negative sensitivity to each of the first 11 canonical variates, declining gradually from the first to the last variate. Combined with information in Figure 8B, this indicates that L3 cells picked different mixes of positive and negative sensitivities to the 11 variates. If we view the canonical feature space defined by the first 11 variates as an 11-dimensional hypercube, Figures 8B,C indicate that L3 cells picked different corners of this hypercube. Unlike L3 cells, L4 cells do not have any preferential sensitivity to the canonical variates (compare white- and black-shaded bars in Figure 8D). The prominent preferential sensitivity of L3 cells is the product of L3 self-organization. As our no-context L3 model shows, L3 self-organization without contextual guidance from surrounding columns also can to some degree enhance cells’ sensitivity to canonical variates (compare black- and gray-shaded bars in Figure 8D), but much less than under contextual guidance.

When evaluated by their responses to moving images of sinewave gratings, 80% of L3 cells fall in the complex-cell category, whereas 20% of L3 cells fall in the simple-cell category (Figure 9A). This suggests that as a group, L3 cells in a mesocolumn should be able to represent both the orientation and position of grating images in their output activity vector φ → (Equation 4). To show how well they can do it, we expose the RF of the central mesocolumn to a grating pattern of randomly selected orientation, spatial frequency, and position in the RF, and compute the φ → 0 response it evokes in the 150 L3 cells in the central mesocolumn. We then rotate the grating pattern by a randomly chosen angle α, compute the new φ → α response of L3 cells, and measure the angle between the two L3 output vectors φ 0 → and φ α → . Figure 9B plots the average L3 angle as a function of the angle between the two gratings (gray curve). For a comparison, we also plot the average angle computed for the output vectors of the 150 L4 cells in the central mesocolumn (black curve). In Figure 9C, instead of rotating, we translate the grating pattern by a randomly chosen fraction (phase) of the grating’s period, compute the new response of L4 and L3 cells, and measure the angle between the two L3 output vectors and between the two L4 output vectors. Figure 9C plots the average L3 (gray) and L4 (black) angles as a function of the phase shift between the two gratings. The plots show that both L4 and L3 output vectors can discriminate even small differences in gratings’ orientation or position in the RF. It is interesting to note that even at the maximal orientation (90°) or spatial phase (180°) differences between two gratings, L3 output vectors show less than maximal (90°) separation, reflecting the fact that other than for their orientation or phase, the two gratings are the same.

Figure 10 demonstrates the importance of contextual guidance for the development of biologically realistic feature properties in L3 cells. To test orientation tuning of L4 and L3 cells in the model, each cell was stimulated with moving grating patterns of the optimal spatial frequency and the full 180° range of orientations. The tightness of the cell’s orientation tuning was expressed by the standard half-width and half-height (HWHH) of the orientation tuning curve. Figure 10 plots the F1/F0 ratio determined for each cell against its HWHH. The model’s L4 cells are shown as blue circles, L3 cells as red dots, and L3 cells of the no-context model as green dots, revealing that all L4 cells are most tuned to orientation (average HWHH = 18°, matching real cat V1) and belong to the simple-cell category, L3 cells are also well-tuned to orientation and have biologically accurate proportion of simple- and complex cell categories, whereas L3 cells in the no-context model fail do develop translational invariance and have poor orientation tuning.

In the real cortex, long-range horizontal connections link cortical columns separated by up to several millimeters in a cortical area. They preferentially link cortical sites that share similar functional properties but have non-overlapping RFs (Gilbert and Wiesel, 1983; DeFelipe et al., 1986; Lund et al., 1993; Burton and Fabri, 1995; Bosking et al., 1997). The fact that input patterns encountered by mesocolumns in their RFs possess contextually predictable features makes it possible for such long-range horizontal connections to establish Hebbian links between distant cortical columns even though they have non-overlapping RFs. When L3 cells tune to contextually predictable features, they become correlated with similar L3 cells in surrounding columns in their stimulation-evoked activities. Figure 11 shows the magnitude of such correlations between L3 cells in the central mesocolumn and functionally identical L3 cells in the first and second ring of surrounding mesocolumns. For a comparison, Figure 11 also shows that even in the first ring of mesocolumns, functionally identical L4 cells have very low correlations, which means that they would not be able to establish lateral Hebbian connections.

Canonical features learned by L3 cells characterize input patterns in many different ways that reflect the orderly aspects of the sensed outside world. This makes it possible for different input patterns, which are in some way objectively related, to be preferentially clustered in the L3 output space. We demonstrate this clustering tendency on an example of 4 different 512 × 512 pixel texture images shown in Figure 12A. We exposed the field of 19 mesocolumns to 500 randomly picked locations in each of the 4 texture images, and for each location we averaged the responses of L3 cells tuned to the same feature across 19 mesocolumns, resulting in a 150-dimensional feature vector representation of the imaged texture field. We next did Principal Component Analysis (PCA) on the 500 × 4 feature vectors. In Figure 12B, we plot the computed scores of the first 3 principal components, color-coding them according to the texture images from which they originated. This 3-D plot reveals that L3 responses to the 4 different textures occupy non-overlapping regions in the principal components space. For a comparison, Figure 12C shows the same plot for responses of thalamic LGN cells, which provided the input to the 19 cortical mesocolumns. As expected, LGN responses to the 4 different textures show no sign of preferential clustering.

We expressed the similarity of L3 output vectors evoked in response to different randomly picked locations in the same texture image by computing correlations between pairs of such vectors. Figure 13 plots average correlations (squared, keeping the sign) for each of the 13 texture images in the Brodatz (1966) database. Also plotted are average correlations computed for the L4 response vectors and vectors of 15 significant canonical variates. As the plot shows, unlike L4 vectors, both L3 cell and canonical variate vectors evoked by different views of the same texture show substantial similarity. Also noticeable is that L3 cells consistently show greater similarity than canonical variates (compare red and green bars), demonstrating the advantages of the overcomplete representation of the mesocolumn’s 15-D canonical feature space by 150 L3 cells.

The results of our cortical model simulations support biological plausibility of the proposed mechanism of cortical feature tuning and offer new insights into the nature of information extraction by neocortical networks. While potential usefulness of contextual guidance for feature tuning has long been recognized (see Introduction), so far it has only been explored at an abstract level or using greatly reduced “toy” models. In this paper we explore actual biological means by which neocortex tunes its neurons to contextually predictable features. Such means are likely to be the backbone of functional organization of cortical columns.

The starting point in our biologically grounded exploration was consideration that the tuned features have to be nonlinear. Our proposed two-stage solution is that, similar to artificial neural networks, feature tuning relies on hidden layer-like preprocessing, performed by the afferent input layer L4. That is, a local group of L4 neurons together perform a nonlinear transform of their thalamic inputs, which is akin to a basic RBF transform. Such transform accomplishes pluripotent function linearization, thus allowing L3 cells in the second stage to extract their features by simple linear summation of their L4 inputs.

Simplifying the task of feature tuning to that of linear operation over the L4 transform allowed us to determine – using an objective-function optimizing algorithm derived in Section 2.4 – that an RF of a representative V1 cortical column is likely to possess around 15 independent contextually predictable features, which we called canonical variates (Figure 4B). Furthermore, we determined that the transform-performing group of L4 cells will need at least 150 members in order to maximize the contextual predictability of all the canonical variates in its RF (Figure 6). Such a number is much greater than the 30–60 excitatory cells found in the L4 compartment of a single 0.05 mm-diameter cortical minicolumn (which is the narrowest columnar entity in the neocortex), but much smaller than the total number of excitatory L4 cells in the 0.5 mm-diameter macrocolumn. This finding leads us to propose a new class of cortical columns, the 0.15 mm-diameter mesocolumn. It is estimated to comprise a local group of 7 minicolumns, all innervated by the same bundle of afferent axons, and to have 200–400 L4 excitatory cells, which together perform the pluripotent function linearizing transform of the mesocolumn’s afferent input.

Moving to the second stage of cortical feature extraction, which takes place in the upper layers, both anatomical evidence and model studies indicate that an L3 pyramidal cell extracts its feature from the L4 input it receives from not just its own mesocolumn, but from a macrocolumn-size group of 7 mesocolumns. Anatomical segregation of the ascending L4 inputs to the basal dendrites and the long-range horizontal inputs to the apical dendrites of L3 pyramidal cells, combined with their separate spike-generating synaptic input integrating centers allows basal and apical dendrites to have separate – classical and extra-classical – RFs and develop RF properties that will maximize covariance of the cell’s apical and basal outputs (Kording and Konig, 2000; Kay and Phillips, 2011). Our model simulations, based on Hebbian plasticity of apical and basal synapses, show that contextual inputs to the apical dendrites readily drive basal dendrites to select contextually predictable (i.e., canonical) features in their classical RFs. Similarly to real V1 cortex, 80% of model L3 cells acquire complex-cell RF properties while 20% acquire simple-cell properties (Gilbert, 1977; Kim and Freeman, 2016). Overall, the design of the model and its emergent properties are fully consistent with the known properties of cortical organization.

If cells in the mesocolumn’s L3 compartment did not push each other to select different features, they would all tune exclusively to the first – most predictable and thus most attractive – canonical variate. However, diversification pressures drive L3 cells to choose the second best solution. Rather than tuning to one or a mix of few of the canonical variates, all L3 cells in the mesocolumn become sensitive to all first 11 variates. This sensitivity declines gradually from the 1st to the 11th variate in all cells (Figure 8C). For each variate, all L3 cells develop approximately the same correlation with it but they differ in the sign of that correlation. Thus, each L3 cell carries maximal information it can about all first 11 variates (rather than emphasizing a subset of variates), with each variate contributing either positively or negatively in the pattern unique to that cell. As a result, L3 can be considered as approximating Hadamard-like domain transform of the first 11 canonical variates, decomposing them into a set of constituent functions (canonical features) over all variates. The L3 transform differs from Hadamard transform in that its transform functions are not orthogonal, and their number (200–400) greatly exceeds the number of variates (11). That is, L3 generates an overcomplete representation of canonical variates.

Considered geometrically, diversification pressures among cells in the mesocolumn’s L3 compartment drive them to choose different preferred directions in the mesocolumn’s canonical feature subspace. If we view the canonical feature space defined by the first 11 variates as an 11-dimensional hypercube, we find that L3 cells pick different corners of this hypercube. Such a hypercube will have 2048 corners to choose from. It is intriguing that while 200–400 L3 cells in a mesocolumn will not be able to pick all these corners, the larger columnar entity comprising a group of 7 mesocolumns – together making up a macrocolumn and sharing L4 input – will have just the right number of L3 cells for such a task.

As pointed out in Introduction, the orderly – as evidenced by their contextual predictability – nature of canonical features reflects the orderly structures in the environment. In tuning selectively to canonical features, L3 performs selective filtering of the information it receives from L4, emphasizing information about orderly aspects of the sensed environment and downplaying local, likely to be insignificant or distracting, information. Despite selective filtering and overcomplete representation of the canonical feature subspace, L3 output preserves excellent discrimination capabilities (Figure 9) while acquiring novel categorization/abstraction ability to preferentially cluster in the L3 output space different input patterns that are in some way objectively related (Figure 12). Furthermore, reduced sensitivity of L3 output to distracting irrelevant details should help the L4 in the next cortical area to minimize the Curse of Dimensionality and to succeed in the next round of pluripotent function linearization, and for the next L3 to find higher-order canonical features.

The general idea of using spatiotemporal coherence to discover useful regularities in inputs was introduced by Becker and Hinton (1992) and later elaborated by Becker (1996). Their IMAX learning procedure discovers regularities in multiview inputs by maximizing mutual information between outputs of two nonlinear multilayer network modules that receive nonoverlapping, but spatially or temporally related, input samples, thus tuning to higher-order input features reflecting common distal causes in the external world. Details of IMAX design, however, make it unsuitable for implementation in the cerebral cortex (Becker, 1996). Phillips and Singer (1997) suggested a way of making computation of mutual information biologically more plausible, and it is one of the cornerstones of their Coherent Infomax theory. They consider abstract local processors, loosely analogous to unspecified local cortical circuits, that receive both the afferent input from their RFs and lateral (contextual field) input from other such local processors. The contextual field input guides local processors to tune to those stimulus features in their RFs that are predictably related to the context in which they occur. According to Coherent Infomax, contextual inputs can be used not only to guide learning but, importantly, also modulate short-term processing of sensory information. Phillips and Singer (1997) derived a particular mechanism for how contextually-guided learning might be accomplished. Unfortunately, that mechanism is limited in its practical utility due to its inability to search for nonlinear correlations. In its later development, Kay and Phillips (2011) showed that Coherent Infomax is consistent with a particular Bayesian interpretation for the contextual guidance of learning and processing and suggested learning rules that are more computationally feasible within systems composed of very many local processors.

Rather than invoking abstract local processors, Kording and Konig (2000) proposed that contextual guidance of feature tuning is implemented in individual pyramidal neurons, in which the apical dendrite acts – in addition to the soma – as a second site of integration capable of generating action potentials. Synaptic inputs to the soma site, coming from the cell’s RF, mainly determine the output activity of the post-synaptic neuron. Contextual inputs to the apical site gate synaptic plasticity. This separation makes it possible for contextual information to avoid confounding the effects of processing and learning. In “toy” simulations of such 2-site neurons receiving nonoverlapping but correlated inputs to their somata while sending their “teaching” outputs to each other’s apical site, cells learned to represent only the coherent part of the input, which would be expected to be relevant to the processing at higher stages. Kording and Konig termed their design Relevant Infomax.

To explain how 2-site pyramidal neurons might be able to tune to nonlinear features in their inputs, the challenge which was not addressed by the Kording and Konig model, Favorov and Ryder (2004) proposed that since dendritic trees are fundamentally nonlinear integrators, they might be able to operate functionally as error backpropagating multilayer perceptrons (MLP). In their SINBAD (acronym for Set of INteracting BAckpropagating Dendrites) neuron model (Ryder and Favorov, 2001), the apical dendrite in each pyramidal cell functions as one MLP and the basal dendrites function as the second MLP, using each other’s output activities as their reciprocal backpropagating teaching signals. While SINBAD cells are very powerful in discovering high-order nonlinear regularities hidden in multiview sensory inputs, effectively approximating Gebelein’s maximal correlation (Kursun and Favorov, 2010), it has become clear that they are not biologically feasible because, while action potentials do backpropagate from the initial axon segment up the apical dendrite, their experimentally observed amplitude modulation is not consistent with what would be required in the error signal. Furthermore, this design depends on a complete separation of the inputs to the apical and basal dendrites, which is not observed in the real cortex. Instead, a much more biologically appealing solution for the necessity of tuning cells to nonlinear features is to make use of pluripotent function linearization in L4 (Favorov and Kursun, 2011), followed by linear learning in L3, as is explored in this paper.

The model of contextual guidance of feature selection explored in this paper is not complete. In addition to spatial context, which was investigated here, contextual guidance can come from temporal context in which orderly features occur, as well as from higher-level understanding of the overall situation. In this paper, we only used static images and thus confined ourselves to spatial features of orderly structures, leaving temporal features of orderly processes for later studies. We anticipate that studies of feature acquisition under temporal contextual guidance and feedback from higher-level cortical areas will make it necessary to expand our current L4-L3 model by adding deep layers and layer 2, resulting in a cortical column model incorporating all cortical layers.

The biological realism of neurons modeled in this paper is not complete. Unlike real neurons, which have binary outputs and are either excitatory or inhibitory, but not both, the modeled cells have outputs that are continuous variables in a negative–positive range and have connections that in the process of learning can change their sign. Adding this degree of biological realism to the model will be insightful, but we do not expect it to negate the lessons learned using the current model. Also, some of the mathematical techniques used in the model, such as normalization of the connection weights in Equations 17, 23, 25, might only approximate the true homeostatic mechanisms in the cortex (e.g., Turrigiano et al., 1998) and should be investigated further.

Sensory cortical columns are engaged not only in feature extraction and sensory information transmission to higher cortical areas, but also in other tasks, such as across-column binding by selective spike synchronization (Uhlhaas et al., 2009; Singer and Lazar, 2016), dynamic contrast enhancement and focused attention (Schummers et al., 2005; Tommerdahl et al., 2010; Tallon-Baudry, 2012), predictive computation (Bubic et al., 2010; Favorov et al., 2015; Marvan and Phillips, 2024; George et al., 2025), etc. Correspondingly, output of real pyramidal cells in L3 is determined not only by synaptic integration of L4 inputs by the basal dendrites, as was done in the current paper, but also by local excitatory and inhibitory inputs, input from the apical dendrite, and other sources (Angelucci and Bressloff, 2006). Our current model lacks all this machinery since its sole purpose was to investigate mechanisms determining classical RF and feature tuning properties of cortical neurons. However, assuming that our proposed mesocolumn-based mechanism of 2-stage feature extraction is biologically realistic, our current model provides a starting point, constraints, and guidance in building a progressively more comprehensive model of cortical functional organization.