Mechanisms for Stable, Robust, and Adaptive Development of Orientation Maps in the Primary Visual Cortex

Measuring stability and selectivity.

The stability of orientation map development is illustrated in Figure 1 with data from chronic optical imaging in ferrets, reprinted from Chapman et al. (1996). Stability and selectivity measures for maps in eight ferrets (of either sex; T. Bonhoeffer, personal communication) are shown in Figure 2 and described below.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Stable development of orientation maps in ferret striate cortex. Orientation maps recorded in the primary visual cortex of one ferret (animal 1-3-3630) at the postnatal ages indicated are shown. In these polar maps, pixel color indicates orientation preference, and pixel brightness indicates the strength of orientation tuning. The selectivity of each map is normalized independently, making blood vessels visible at P31 (e.g., the blue streak in the top left corner), but not at P42, once orientation selectivity has developed. As the maps mature, iso-orientation domains become visible as colored regions that become more strongly responsive over time, without changing the overall map pattern. Data are reproduced from the study by Chapman et al. (1996).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Development of orientation maps measured by chronic optical imaging in ferrets. A, Recorded selectivity (red square markers) and stability values (blue round markers) for all eight ferrets, as a function of postnatal day [replotted data from the study by Chapman et al. (1996), their Figs. 4, 7]. Stability is quantified by the average difference in preference of each orientation map with the final map, as defined by Equation 1. A common selectivity and stability scale is used to allow comparison between ferrets. All values are in arbitrary normalized units, using the lowest recorded value as the zero reference point and the highest recorded value as the maximum (see Materials and Methods). B, The mean selectivity and stability across ferrets as a function of postnatal day, with the ±95% confidence intervals for each day indicated by the shaded area around the mean line. Selectivity and stability increase steadily and simultaneously over development so that once neurons are selective, they are organized into a map with the same form as the final measured map. These results show that selectivity does not precede increasing stability, indicating stable map development. C, As a reference for later modeling work, the final orientation preference map (without selectivity) for ferret 1-3-360 on day P42 (Fig. 1, rightmost plot) is shown. The map is organized into regularly repeating hypercolumns in all directions, as seen by the ringness in the Fourier power spectrum, plotted with the highest amplitude component in black. The center value is the DC component, and the midpoint of each edge represents half of the highest possible spatial frequency in the cardinal directions (i.e., the Nyquist frequency). This Fourier spectrum is used to calculate the average periodicity of the map, which is then plotted as a hypercolumn area Λ², covering one period in the cardinal directions (white boxed area). Figure 3 illustrates how to use these calculations to determine whether a model map resembles this animal data.

Qualitatively, map development is considered stable if the map pattern remains constant once neurons have become highly selective; unstable development would be characterized by maps that reorganize substantially after orientation-selective patches first emerge (Chapman et al., 1996). In electrophysiological measurements from ferret V1, average selectivity (orientation tuning) is low from about postnatal day 23 (P23; when visual responses can first be measured) until P28–P34, reaching adult levels by P42–P49 (Chapman and Stryker, 1993). To examine stability of maps using optical imaging, Chapman et al. (1996) thus focused on the period from P31 to P45. This period also corresponds to the critical period for visual experience—disrupting activity before P50 eliminates selectivity even with normal visual experience after this age (Chapman and Gödecke, 2000). For this modeling study, we will also focus on this period of normal orientation selectivity emergence; stability after the critical period can be modeled trivially by disabling or greatly reducing plasticity, and will not be investigated further here.

The red lines with square markers in Figure 2A show the average selectivity of the optical imaging signal from P31–P45 for all eight ferrets (Chapman et al., 1996); one additional data point from ferret 1-1-4479 at P55 is not shown. Although the selectivity values are on an arbitrary scale due to the measurement technique, the values for each postnatal age correspond closely to those found electrophysiologically (Chapman et al., 1996), and thus cover the range from highly unselective to adult-like selectivity.

To assess stability quantitatively over this same time period, Chapman et al. (1996) computed the correlation of orientation preference at each developmental age with the organized preference map observed in the final recording for that animal. For convenience, we have linearly rescaled the “orientation similarity index” values from that study to vary from 0.0 (uncorrelated preference) to 1.0 (identical preference), which we compute as a stability index (SI): where i iterates over all n data points (pixels in the imaging frame), F_i is the final orientation preference value for the neural unit corresponding to pixel i, and O_i is the corresponding preference value at an earlier age. The blue lines with round markers in Figure 2A show the results of this calculation for eight ferrets from Chapman et al. (1996) measured at different postnatal days. Note that even identical underlying maps would give values below 1.0 on this measure, if there is any noise or variation in measuring the map, and so the values shown have been normalized by the highest and lowest stability values across all ferrets studied.

Given this data, stable development is defined as a map having high SI values whenever the average selectivity value is high. When the selectivity is very low, presumably before the neurons have developed at all, the map is arbitrary and dominated by measurement noise, but if development is stable, the SI value should increase as soon as the selectivity value increases. Unstable development would be visible as a selectivity value that increases well before stability does, representing neurons that have achieved selectivity before assuming the final preference value observed in the last map measured. As can be seen in Figure 2A, the SI value and the average selectivity are highly correlated in every ferret shown, with both values typically increasing over time rather than selectivity predating stability (though ferret 1-1-4479 could be seen as relatively unstable). An overall correlation between selectivity and stability is a hallmark of stable map development, and is visible as a general trend in the average data from all ferrets, plotted in Figure 2B.

Assessing the organization of orientation maps.

Recording the selectivity and stability index is sufficient for evaluating development in experimental recordings, where all maps are necessarily biologically realistic. In simulated models of map development, there is no guarantee that simulated maps have an orientation structure matching those observed in real animals. An unrealistic map may be both stable and selective, and so we require that a good model of map formation also results in realistic maps, similar to those found in animals. Ideally, map quality would be assessed with an automated metric that reports how close the maps are to animal data, e.g., on a scale from 0 to 1.0 as for stability.

There have been several previously published measures that could be used to quantify the degree of map organization. Biological maps tend to be smooth, which can be measured by the local homogeneity index (LHI; Nauhaus et al., 2008) or the earlier “local input orientation selectivity index” (Schummers et al., 2004), which both convey the local homogeneity in orientation preference at a particular cortical location. Biological maps typically have non-uniform orientation histograms (Coppola et al., 1998; Müller et al., 2000; Tanaka et al., 2009), which have been used to characterize model maps as well (Bednar and Miikkulainen, 2004). Biological maps tend to have a ring-shaped Fourier power spectrum (Blasdel, 1992a,b; Erwin et al., 1995; Fig. 2C), due to the regular repetition of orientation patches across the cortical surface. Orientation maps also contain pinwheels that may be identified, classified by polarity (whether orientation preference increases clockwise or counterclockwise) and analyzed; e.g., the average distance between nearest-neighbor pinwheels tends to be greater between those of matching polarity than those of opposite polarity (Müller et al., 2000).

Although these measures can each help quantify the realism of simulated maps, they do not offer a simple, unique reference value that can characterize all biological maps. A useful metric should assign biological maps to a narrow possible range of values, ideally with theoretical justification for what values should be expected in the biological system and why. Furthermore, a reliable metric should be consistently defined between species, if it is to avoid uncertainty in fitting species-dependent free parameters.

Accordingly, we developed a new map-quality metric based on the empirical observation that pinwheel count in biological orientation maps scales linearly with hypercolumn size across many different species (Kaschube et al., 2010). This linear relation gives a consistent number of pinwheels per hypercolumn area (Λ²), implying a constant pinwheel density when averaged over sufficiently large cortical areas. This dimensionless, statistical measure of pinwheel distribution is thought to reflect a universal constant of map organization, converging to π across carnivorans, primates, cats, and tree shrews (Kaschube et al., 2010; Keil et al., 2012). This value was predicted by a theoretical model of map organization and has strong empirical evidence, with a mean pinwheel density across four species (tree shrew, galago, cat, and ferret) statistically indistinguishable from π (Kaschube et al., 2010); see data from three species in Figure 3D.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Evaluating map quality as the deviation from π pinwheel density (ρ). A, High-quality, realistic, orientation preference and selectivity map with approximately π pinwheel density (3.146), from the final model to be discussed in this paper (GCAL). Corresponding preference-only map is shown in Ca. Pinwheel density is defined as the average number of pinwheels (white circles) per hypercolumn area (Λ²) indicated by the white boxed area (A, B). Pinwheels are identified at the intersection of the zero contours of the real and imaginary components in polar representation (white and black contours respectively). The periodicity of hypercolumns is estimated from the radius of the ring in the Fourier power spectrum (FFT) using the fitting method described in Kaschube et al. (2010). B, A lower-quality map generated from the first model introduced in this paper (L), which has visible discontinuities in OR preference seen in Cb. The greater pinwheel density (5.917) is due to a higher pinwheel count and a larger hypercolumn area due to more widely spaced orientation blobs. The histogram (FFT plot inset) indicates mean spectral power as a function of radius, the red line indicates the least-squares fit (Kaschube et al. 2010, their supporting online material, Eq. 7), and the blue arrow indicates the estimated peak spectral power radius (A, B). C, A selection of model maps ordered by pinwheel density with pinwheel count to hypercolumn area ratio, shown in parentheses. Lower-quality maps usually have higher estimated pinwheel counts and correspondingly higher pinwheel densities, with pinwheel counts so large for very poor maps as to be effectively undefined. D, Pinwheel density of three species (diamonds) and simulated maps (circles) as a function of hypercolumn size [data replotted from the study by Kaschube et al. (2010)]. Horizontal lines indicate median values of each cluster with the medians of high quality model maps also clustered around π. E, Normalized, heavily tailed gamma distribution used to transform pinwheel density into a suitable metric on a unit interval with values for maps a, b, and c.

To establish the pinwheel density for any given map, the total pinwheel count must be determined. Pinwheel centers are located at the intersection of the zero contours of the real and imaginary components in the polar representation of orientation preference (Löwel et al., 1998). These contours are shown in Figure 3A as black and white lines for a sample simulated preference and selectivity map with approximately π pinwheel density. The preference-only channel of this high-quality map is shown together with corresponding pinwheel density and pinwheel count in Figure 3Ca. All the maps shown in Figure 3 are for illustration only and are derived from simulations to be described later; simulated data are necessary to demonstrate the behavior of the metric for nonbiological maps.

For simulated maps using uniform random input statistics, a single number is sufficient to characterize the Fourier plots: the radius of the isotropic ring, computed using the methods described in the study by Kaschube et al. (2010). To eliminate disruptions in hypercolumn size due to border effects at the sheet edges, only the central 1.0 × 1.0 area from a simulated V1 sheet of area 1.5 × 1.5 is analyzed for all model maps presented. The first step is to integrate the spectral power as a function of radius, as shown by the histograms in the FFT plots. The peak values of these histograms are estimated using a least-squares fit of a Gaussian curve with additional linear and quadratic terms (Kaschube et al., 2010, their Equation 7, supplementary materials). The computed fits are shown by the red curves overlaying the histograms with the fitted peaks indicated by the blue arrows, marking the estimated ring radius. This value reflects the periodicity with which hypercolumns repeat across the map and is equal to the hypercolumn area in units of Λ² when squared, shown for a ferret map in Figure 2C and for simulated maps in Figure 3, A and B.

An identical analysis is shown for a lower-quality simulated orientation map in Figure 3B, where the map was generated by running a simulation outside its optimal operating range. Disruptions to the smooth organization of preference result in clumps of pinwheels, explaining the higher pinwheel density seen in Figure 3Cb. These clumps of pinwheels could in principle be eliminated arbitrarily as duplicates, but because they represent genuine areas of poorly organized preferences in these noise-free simulations, they have been retained as indicators of poor map structure. This trend is illustrated in Figure 3C by a selection of maps of increasingly poor quality that illustrate how the pinwheel count and density tend to increase as map quality decreases.

The validity of using π pinwheel density as a reference value is demonstrated in Figure 3D. The experimental values of pinwheel density for tree shrews, galagos, and ferrets are plotted as diamonds, using data reproduced from the study by Kaschube et al. (2010). The values for simulated maps are plotted as circles for the high- and low-quality map simulations (n = 40 and n = 33 respectively). The high-quality maps are obtained from the final model developed in this paper, while the lower-quality maps are computed using the initial, simplified model introduced in the next section. The horizontal lines through each cluster indicate the median value of the cluster points. The median is indistinguishable from the mean for all clusters except the low-quality map cluster, which has outliers of relatively high pinwheel density not visible in the plot.

It would be interesting to evaluate pinwheel density over development using the chronically recorded ferret data, illustrated in Figure 1. This measurement is not possible without access to the raw data, as preprocessing may have introduced biases to the orientation structure that need to be carefully controlled for when assessing pinwheel density (Kaschube et al., 2010). The pinwheel density at the end of maturation in ferrets has been established using 82 maps recorded from adult animals, shown by a cluster of green diamonds centered around π in Figure 3D.

For the poorest-quality maps, the pinwheel density is unbounded, as illustrated by Figure 3C. In these highly disorganized maps (to the right of Fig. 3C), an unrealistically large number of locations match the automated criterion for a pinwheel; few or any of these would satisfy any subjective definition of a pinwheel. Even so, the “pinwheel count” is still a useful measure in these cases, because such high numbers reliably indicate low-quality maps.

To turn pinwheel density into a useful metric between unity (high-quality maps) and zero (low-quality maps), a heavily tailed squashing function is needed, mapping all the maps with unrealistically high pinwheel counts to zero. We chose a normalized gamma distribution, shown in Figure 3E. The labels a, b, and c show how the final metric value is computed for the three example maps shown. The gamma distribution Γ(k, θ) is characterized by k = 1.8 and θ = π/(k − 1). This latter constraint on θ ensures that the mode of the curve occurs at π. This kernel is then normalized by the value evaluated at π to ensure that maps with exactly π pinwheel density have a metric value of 1. The only free parameter k is set to the value k = 1.8 to allow detailed discrimination of high-quality maps without discarding all the information available for poor maps with a pinwheel density up to 30; other values will place more or less penalty for very poor maps with unrealistically high pinwheel densities, but will otherwise lead to similar results.

Of course, a value of unity on the pinwheel density metric (or any metric) is not sufficient by itself to guarantee that a given simulated map is indistinguishable from a biological recording. For instance, it should be possible to construct a synthetic map designed specifically to attain π pinwheel density that nonetheless appears highly unnatural. Such deliberate “gaming” of the map metric is not possible using the simulations presented in this paper, as the maps emerge gradually from a developmental process, without modeler control over the specific placement of pinwheels in the final organization. Even so, it is important to verify that simulated maps with a high metric value satisfy the types of subjective criteria used to evaluate biological results, such as those introduced by Blasdel (1992b). For each result where maps are evaluated by metric, we have thus also performed subjective evaluations of map quality, such as assessment of whether the FFT is ring shaped, and whether there are saddle points, linear zones, 180° pinwheels, fractures, and so on. We have found no examples of simulated maps that achieve a realistic pinwheel density but do not also look realistic with respect to these other properties. The automated pinwheel density metric value thus correlates well with our overall subjective assessment of map quality. We have also ensured that all the plots presented in this paper are representative of those simulations that are not shown, for the same conditions.

Model architecture.

The mechanistic models we will consider are variants of the original SOM (self-organizing map) algorithm introduced by von der Malsburg (1973). In this type of model, the inputs are actual visual images or patterns of spontaneous activity, which can be directly related to the visual environment or to the imaging of retinal activity patterns. Biologically plausible properties of single neurons can be integrated into the network to explain many of the observed phenomena, as opposed to just the geometric properties of the map pattern.

Unlike correlation-based learning models that use Hebbian learning over large batches of inputs (Linsker, 1986; Miller, 1994), self-organizing map models operate using incremental Hebbian learning rules. Although the mathematical structure of the final model is often less amenable to analysis than those of linear, feedforward networks, the incremental nature of self-organizing maps make them more suitable for studying map development. Incremental learning allows gradual changes in network organization to be tracked as a stream of inputs drive the development of the network forward.

The models presented here self-organize using the same principles as the LISSOM (Laterally Interconnected Synergetically Self-Organizing Map) algorithm (Miikkulainen et al., 2005), which was inspired by earlier SOM models (Kohonen, 1982; Obermayer et al., 1990). The simplest model presented here, model L (laterally connected), may be considered a simplified version of the LISSOM model, omitting ad hoc modeler-determined changes to the Hebbian learning rate, activation thresholds, and lateral excitatory radii over time. The fundamental operation of all these self-organizing map models has been explained in terms of dimensionality reduction, specifically a discretized approximation of the principal surface of the input (Ritter et al., 1992).

The architecture of the four models evaluated in this paper is shown in Figure 4. The models each consist of four sheets of neural units representing the input (retinal photoreceptors), the ON-center and OFF-center pathways from the photoreceptors to V1 via the retinal ganglion cells and the lateral geniculate nucleus, and V1. A sheet is a two-dimensional array of firing-rate point neurons, with activation and plasticity equations as described below. The simplest variant, model L (laterally connected) consists of four afferent connections (one from the photoreceptor sheet to each ON/OFF sheets and one from each ON/OFF sheet to V1) with lateral connectivity only in the V1 sheet.

The AL (adaptation, laterally connected) model has the same set of connections as L but includes homeostatic adaptation in V1, described by the equations below. The GCL (gain control, laterally connected) model has the same V1 sheet structure as L but includes lateral, divisive inhibition in the ON/OFF sheets that implement contrast-gain control. The final GCAL model includes both the homeostatic adaptation of the AL model as well as the contrast-gain control of the GCL model. Apart from these specific differences, each model shares the same parameters and mechanisms, and can thus be compared directly against the others.

As model L has the same basic architecture as the LISSOM model, the spatial extents of connections and weight patterns are taken from and described in the study by Miikkulainen et al. (2005). Models GCL and GCAL introduce one new spatial parameter, determining the lateral extent of the contrast-gain control mechanism in the LGN layers. No parameters need to be specified with very high precision for our conclusions to hold; in the GCAL model, two significant digits of precision are sufficient to develop qualitatively indistinguishable maps. We estimate that all parameters may be changed by ∼10% without affecting the overall behavior (except for the homeostatic smoothing parameter, which may not be >1.0).

The models are implemented in the Topographica simulator, freely available at www.topographica.org. Topographica allows simulation parameters to be specified in measurement units that are independent of the level of detail used in any particular run of the simulation. To achieve this, Topographica provides multiple spatial coordinate systems, called sheet and matrix coordinates. Parameters (such as sheet dimensions) are expressed in continuous-valued sheet coordinates. In practice, of course, sheets are discretized using some finite matrix of units. Each sheet has a density, which specifies how many units (matrix elements) in the matrix correspond to a unit length in sheet coordinates. Each of these simulations used a density of 98 neural units per sheet coordinate in V1, and 24 units per sheet coordinate in the retinal photoreceptor and ON/OFF sheets. Results are independent of the density used, except at very low densities where discretization artifacts can become prominent.

In all simulations shown, the area of the V1 sheet plotted and analyzed has sheet-coordinate dimensions 1.0 × 1.0, thus consisting of 98 × 98 neural units. The underlying simulated V1 area is actually 1.5 × 1.5 (147 × 147 neural units), which is then cropped to 1.0 × 1.0 for analysis and display, to eliminate border effects due to partial patterns of lateral connectivity for neurons near the edge of the cortical sheet. The ON/OFF channels and photoreceptors are similarly extended for simulation to 3.0 × 3.0 and 3.75 × 3.75 in sheet coordinates, respectively, to avoid having any afferent receptive field cropped off, and to avoid edge effects in the ON and OFF channels due to gain control connections when present. In Figure 4, all sheets are shown at the same size for visibility, but the actual area mapped topographically between sheets was the central 1.5 × 1.5 region of each sheet that corresponds to V1.

Temporal properties.

As a simplification, we have ignored the detailed temporal properties of the subcortical neural responses and of signal propagation along the various types of connections. Instead, the model ON/OFF units have a constant, sustained output, and all connections have a constant delay, independent of the physical length of that connection. One training iteration in the model represents one visual fixation (for natural images) or a snapshot of the relatively slowly changing spatial pattern of spontaneous activity (for retinal waves); i.e., an iteration consists of a constant retinal activation, followed by processing at the ON/OFF and cortical levels.

For one iteration, assume that input is drawn on the photoreceptors at time t, and the connection delay is defined as δt = 0.05. Then, at t + 0.05, the ON/OFF cells compute their initial activation; at time t + 0.10, the ON/OFF cells incorporate lateral inhibition for gain control (if present in this model); and at t + 0.15, V1 begins computing. At t + 0.20 and every 0.05 iterations until t + 1.0, V1 continues to compute in a series of settling steps, a fixed number sufficient to ensure that cortical activity is no longer changing significantly (data not shown). During this period, the retinal activity is assumed to be constant, to reduce computational requirements.

Stimuli.

Images are presented to the model at each iteration by activating the retinal photoreceptor units. The activation value ψ_i of unit i in the photoreceptor sheet (P) is given by the gray-scale value in the chosen image at that point.

The images used here are either oriented and elongated two-dimensional Gaussian patterns, disk-shaped patterns of noisy activation, or natural image patches, which together cover a range of different input types so that we can evaluate the robustness of map development. The center coordinates and orientation of each Gaussian are chosen from a uniform random distribution and cover an area of the photoreceptor sheet that is a third wider and taller than V1, so that even V1 units near the borders will see an approximately uniform distribution of each part of the Gaussian blobs. For Gaussian patterns, the contrast is defined as the percentage of the input range 0.0 to 1.0.

Noisy disk patterns create a circular region of activity with a radius that is a substantial portion of the size of the photoreceptor sheet to generate distinct edges; the activation is smoothed at the edges of the disk with a Gaussian blur (Bednar and Miikkulainen, 2004). Disk centers are chosen from a uniform random distribution larger than the retinal photoreceptor sheet size to allow for V1 units to see a uniform distribution of disk areas. Uniform zero-mean random noise is then added to the disk pattern.

Natural image patterns are photoreceptor-sheet-sized patches from images of natural objects and landscapes, from a data set by Shouval et al. (1996). The pictures were taken at Lincoln Woods State Park in Rhode Island, and scanned into a 256 × 256 pixel image. No corrections were used for the optical distortions of the instruments, and there was no preprocessing of the images. The data set used to approximate vertical goggle rearing shown in Figure 11 consists of the same set of image patches, but blurred by convolution with an anisotropic Gaussian kernel. The convolution kernel is a 128 × 128 matrix, consisting of a centered, vertical Gaussian pattern with an aspect ratio of 10 (σ_x = 0.25, σ_y = 0.025 in unit sheet coordinates) and total weight of unity.

The ON/OFF sheets.

At each iteration, a new retinal input is presented, and the activation of each unit in each sheet is updated. Neurons in all sheets are firing-rate point neurons, with a state characterized by an activation level. For all models, the activation level η for a unit at position j in an ON/OFF sheet O at time t + δt is defined as follows: The constant γ_O = 14.0 is an arbitrary multiplier for the overall strength of connections from the photoreceptor sheet to the ON/OFF sheets, chosen to give typical activations in the range 0.0 to 1.0, whereas γ_S is the strength of the feedforward contrast-gain control; ψ_i is the activation of unit i in the two-dimensional array of neurons on the photoreceptor sheet from which ON/OFF unit j receives input (its afferent connection field F_j,P), and η_i,O(t) is the activation of other ON/OFF units on the previous time step (received over the suppressive connection field F_j,S). The activation function f is a half-wave rectifying function that ensures the activation of ON/OFF units is always positive.

In the L and AL models where contrast-gain control is not applied, k = 1 and γ_S = 0, whereas in the GCL and GCAL models, k = 0.11 and γ_S = 0.6. The constant k ensures that the output is always well defined for weak inputs, and γ_S is chosen to rescale activation values so that the numerical results are comparable with and without gain control.

The weights ω_ij represent the fixed connection weights from photoreceptor i to the ON or OFF unit j defined with a standard difference-of-Gaussians kernel. The connection fields for ON units have a positive center and negative surround, and vice versa for OFF units. More precisely, the weight ω_ij from an ON-center cell at location (0, 0) in the ON sheet and a photoreceptor sheet in location (x, y) on the photoreceptor sheet is given by the following: The width of the central Gaussian is defined by σ_c = 0.037, and σ_s = 0.15 determines the width of the surround Gaussian, where Z_c and Z_s denote the normalization constants that ensure the center and surround weights each always sum to 1.0. The weights for an OFF-center cell are the negative of the ON-center weights (i.e., surround minus center). The center of the connection field of each ON/OFF unit is mapped to the location in the photoreceptor sheet corresponding to the location of that unit in sheet coordinates, making the projection retinotopic.

The weights ω_ij,S in the denominator of Equation 2 specify the spatial profile of the lateral inhibition received from other ON/OFF units when contrast-gain control is active. The weights of these connections have a fixed, circular Gaussian profile so that a neuron located at (0, 0) in either the ON or OFF sheet will have the following weights: where (x, y) is the location of the presynaptic neuron, σ_S = 0.125 determines the width of the Gaussian, and Z_S is a normalizing constant that ensures that the total of all the lateral inhibitory weights ω_ij to neuron j sum to 1.0. When contrast-gain control is enabled, these recurrent lateral connections are activated once per iteration, before activity is sent to the V1 sheet.

The V1 sheet.

Each V1 neuron in each model receives connections from three different connection types or “projections” (p), i.e., the afferent projection from the ON/OFF sheets (both channels concatenated into one input vector; p = A), the recurrent lateral excitatory projection (p = E), and the recurrent lateral inhibitory projection (p = I) from other V1 neurons.

The contribution C_j,p to the activation of unit j from each projection type (p = A, E, I) is calculated as follows: where η_i,p is the activation of unit i taken from the set of neurons in V1 to which unit j is connected (its connection field F_j), and w_ij,p is the connection weight from unit i in V1 to unit j in V1 for the projection p. Afferent activity (p = A) remains constant after the first update from the retina, but the other contributions change over 16 settling steps, depending on the activity in V1.

The contributions from all three projections to V1 (afferent and lateral) described above are combined using Equation 6 to calculate the activation of a neuron j in V1 at time t: The projection strength scaling factors for each projection type p are γ_A = 1.5, γ_E = 1.7, and γ_I = −1.4 for all models, set to provide a balance between excitation and inhibition, and between afferent and lateral influences, to provide robust formation of activity bubbles that allows smooth maps to form. For models L and GCL, f is a half-wave rectifying function that ensures positive activation values. In the case where single-neuron automatic adaptation is included (models AL and GCAL), f has a variable threshold point (θ) dependent on the average activity of the unit as described in the next subsection, but in all cases the gain is fixed at unity.

Once all 16 settling steps are complete, the settled V1 activation pattern is deemed to be the V1 response to the presented pattern. At this point we use the V1 response to update the threshold point (θ) of V1 neurons (using the adaptation process described below) and to update the afferent and lateral inhibitory weights via Hebbian learning. V1 activity is then reset to zero, and a new pattern is presented. Note that both adaptation and learning could instead be performed at every settling step, but this would greatly decrease computational efficiency.

Adaptation.

To set the threshold for activation, each neuron unit j in V1 calculates a smoothed exponential average of its settled activity patterns (η̄_j): The smoothing parameter (β = 0.991) determines the degree of smoothing in the calculation of the average. (η̄_j) is initialized to the target average V1 unit activity (μ), which for all simulations is η̄_jA (0) = μ = 0.024. The threshold is updated using the following: where λ = 0.01 is the homeostatic learning rate. The effect of this scaling mechanism is to bring the average activity of each V1 unit closer to the specified target. If the activity in a V1 unit moves away from the target during training, the threshold for activation is thus automatically raised or lowered to bring it closer to the target. Note that an alternative rule with only a single smoothing parameter (rather than β and λ) could be formulated, but the rule as presented here makes it simple for the modeler to set a desired target activity μ.

Learning.

Initial connection field weights are isotropic 2D Gaussians for the lateral excitatory projection and uniformly random within a Gaussian envelope for afferent and lateral inhibitory projections. Specifically, a neuron located at (i, j) will have the following weights: where (x, y) is the sheet-coordinate location of the presynaptic neuron, u = 1 for the lateral excitatory projection (p = E), and u is a scalar value drawn from a uniform random distribution for the afferent and lateral inhibitory projections (p = A, I), σ_p determines the width of the Gaussian in sheet coordinates (σ_A = 0.27, σ_E = 0.025, σ_I = 0.075), and Z_p is a constant normalizing term that ensures that the total of all weights ω_ij to neuron j in projection p is 1.0. Weights for each projection are only defined within a specific maximum circular radius r_p (r_A = 0.27, r_E = 0.1, r_I = 0.23).

In the model, as images are presented to the photoreceptors, V1 afferent connection weights ω_ij,A from the ON/OFF sheets are adjusted once per iteration (after V1 settling is completed) using a simple Hebbian learning rule. This rule results in connections that reflect correlations between the presynaptic ON/OFF unit activities and the postsynaptic V1 response. Hebbian connection weight adjustment at each iteration is dependent on the presynaptic activity, the postsynaptic response, and the Hebbian learning rate: where for unit j, α is the Hebbian learning rate for the afferent connection field F_j. Unless it is constrained, Hebbian learning will lead to ever-increasing (and thus unstable) values of the weights (Rochester et al., 1956). In all of the models, the weights are constrained using divisive postsynaptic weight normalization (Eq. 10), which is a simple and well-understood mechanism. Afferent connection weights from ON and OFF units are normalized together in the model. We expect that a more biologically motivated homeostatic mechanism for normalization such as multiplicative synaptic scaling (Turrigiano, 1999; Turrigiano and Nelson, 2004; Sullivan and de Sa, 2006) or a sliding threshold for plasticity (Bienenstock et al., 1982) would achieve similar results, but we have not tested these.

The learning rates for the afferent projection, lateral excitatory projection, and lateral inhibitory projections are α_A = 0.1, α_E = 0.0, and α_I = 0.3, respectively. The density-specific value used in the equation above is then calculated as α = α_A/τ_A, where τ_A is the number of connections per connection field in the afferent projection. To increase computational efficiency, lateral excitatory connections do not learn during development in the simulations presented here. The effect of lateral learning has been explored in detail in previous similar models (Miikkulainen et al., 2005), and the maps generated by the GCAL model are visually indistinguishable regardless of whether or not lateral excitatory connections are plastic (data not shown).

Analysis of model maps.

Model orientation maps are calculated based on the vector average method (Blasdel and Salama, 1986; Miikkulainen et al., 2005). Sine grating inputs that cover the full range of parameter values (combinations of all orientations, frequencies, and phases) are presented, and for each orientation, the peak response of the neuron is recorded. The orientation preference is calculated by constructing a vector for each orientation θ (between 0 and 180°), with the peak response as the length and 2θ as its orientation. These vectors are summed and the preferred orientation is calculated as half of the orientation of the summed vector. The selectivity is given by the magnitude of the summed vector. Average orientation selectivity values are the mean value of orientation selectivity across all units in the map and are normalized by dividing by the maximum average selectivity measured across all simulations reported in this paper.

Orientation tuning curves are measured by presenting sine gratings with 20 different orientations, each at eight different phases and for a range of contrasts. For each contrast, the responses of each neuron unit are measured, and the maximum response at each orientation over all phases is recorded. The tuning curve is constructed from these maximum responses.

To ensure consistent map measurement across all conditions, orientation maps in Figures 3⇓⇓⇓⇓⇓–9 are measured before lateral interactions and the V1 activation function take affect (i.e., on the afferent input activity only). Maps for models with gain control change little if these steps are included. However, for models without gain control, if activation thresholds and lateral interactions were simulated during map measurement, maps would need to be measured at many different contrasts, because orientation tuning in these models is not contrast invariant, and therefore orientation selectivity is difficult to determine. Instead, we use the orientation selectivity of the afferent connections as a measure of the map selectivity in all cases. This simplification allows fair and consistent comparison of selectivity values between the simulations of all four models. Similarly, connection field plots show the weights of connections between the ON sheet and V1 at different times during development. These weights are used as an approximation to the receptive fields that can be measured using a computationally intensive reverse correlation process that takes into account the full connectivity, which again must control for contrast in the models without gain control. The true orientation selectivity of neurons in the full recurrent network is reflected in the orientation tuning curves for the GCAL model (as shown in Fig. 12).