A model for the origin and development of visual orientation selectivity

Orientation selectivity is a key property of primary visual cortex that contributes, downstream, to object recognition. The origin of orientation selectivity, however, has been debated for decades. It is known that on- and off-centre subcortical pathways converge onto single neurons in primary visual cortex, and that the spatial offset between these pathways gives rise to orientation selectivity. On- and off-centre pathways are intermingled, however, so it is unclear how their inputs to cortex come to be spatially segregated. We here describe a model in which the segregation occurs through Hebbian strengthening and weakening of geniculocortical synapses during the development of the visual system. Our findings include the following. 1. Neighbouring on- and off-inputs to cortex largely cancelled each other at the start of development. At each receptive field location, the Hebbian process increased the strength of one input sign at the expense of the other sign, producing a spatial segregation of on- and off-inputs. 2. The resulting orientation selectivity was precise in that the bandwidths of the orientation tuning functions fell within empirical estimates. 3. The model produced maps of preferred orientation–complete with iso-orientation domains and pinwheels–similar to those found in real cortex. 4. These maps did not originate in cortical processes, but from clustering of off-centre subcortical pathways and the relative location of neighbouring on-centre clusters. We conclude that a model with intermingled on- and off-pathways shaped by Hebbian synaptic plasticity can explain both the origin and development of orientation selectivity.


Introduction
neuron are co-extensive; a Hebbian development process functionally segregates the two signs of input. Intracortical inhibition, which is assumed to derive from a widespread, slow-acting network that indiscriminately reduces membrane potential, contributes to orientation selectivity through the iceberg effect [17]. Our aim is to show how developmental refinement of geniculocortical connections can lead to cortical response characteristics-time courses, precise orientation tuning, and orientation preference maps-that match well with those seen in the laboratory.
To make the modelling manageable, the scope of the model is limited in several ways. First, the model is designed to describe the cat's visual pathway because the visual literature for this species is particularly rich (including almost all the animal studies cited above). Second, the subcortical pathway is chosen to pass through the X-type retinal ganglion cell because of its relatively high acuity. Last, the model is restricted to monochromatic, monocular stimuli, and the input layers of primary visual cortex. The model builds on a previous one [18] by adding a retinal ganglion cell array, a development process and dynamic intracortical inhibition. Earlier accounts of this work have appeared in abstract form [19], [20].

Model structure
A flow diagram of the model is shown in Fig 1A. There are multiple subcortical channels, each of which passes through either on-centre or off-centre neurons. Each channel consists of four neurons-photoreceptor, bipolar cell, retinal ganglion cell and relay cell in the dorsal lateral geniculate nucleus-in series. The input to each channel is a dot product of the stimulus with a Gaussian weighting function representing subcortical spatial spread. Apart from shared input, signals in each channel are assumed to be independent. Thus, the membrane potential p in neurons 1, 2, 3 and 4 for channel j is given by: t dp j1 ðtÞ dt ¼ À g js ðx; yÞ � sðt; x; yÞ À p j1 ðtÞ t j dp j2 ðtÞ dt ¼ À n j p j1 ðtÞ À p j2 ðtÞ t j dp j3 ðtÞ dt ¼ p j2 ðtÞ þ p s À p j3 ðtÞ t j dp j4 ðtÞ dt ¼ hðp j3 ðtÞÞ À p j4 ðtÞ The derivation of these equations is provided in the Methods, along with definitions and values of the variables. Subcortical channels converge onto layer 4 and 6 neurons in primary visual cortex via a Gaussian convergence function. Cortical neurons are of two types, excitatory and inhibitory. They receive the same subcortical input, and inhibitory neurons therefore have receptive fields similar to those of excitatory neurons [8,21]. Inhibitory neurons differ from excitatory neurons in that they are split into two compartments. The first compartment, consisting of dendrites and soma, has fast dynamics corresponding to fast-spiking neurons [21]. Simple cells stimulated with flashed stimuli, however, have long-lasting inhibitory tails [22], so we have included a second compartment-axon and terminals-with much slower dynamics. The axons converge on excitatory neurons with another Gaussian convergence function. The effect of this slower component is to provide a widespread, slow-acting inhibition to excitatory neurons [14,22]. Defining inhibitory somas, inhibitory axon terminals and excitatory neurons as stages 5, 6 and 7, the membrane potential in neuron k for each of these stages is given by: t dp k5 ðtÞ dt ¼ X j g kc ðx j ; y j Þw jk hðp j4 ðtÞÞ À p k5 ðtÞ t inh dp k6 ðtÞ dt ¼ hðp k5 ðtÞÞ À p k6 ðtÞ t dp k7 ðtÞ dt ¼ X j g kc ðx j ; y j Þw jk h p j4 ðtÞ � � À X l g ke ðx l ; y l Þp l6 ðtÞ À p k7 ðtÞ Fig 1B illustrates signal processing in a typical neuron. The sum of the weighted synaptic inputs is integrated over time to produce a generator potential. This potential is rectified to produce action potential rate; the exceptions are the photoreceptors and bipolar cells, which do not produce action potentials. Fig 1C shows the signal processing in an inhibitory neuron. Subcortical inputs undergo temporal integration as they pass through the dendrites and soma. The rectification function is then applied to the sum. The resulting action potential is once again integrated within the axon and inhibitory network to produce the inhibitory input to excitatory cortical cells.

Segregation of on-and off-channels
Our aim in this paper is to describe a physiologically plausible model that reproduces key aspects of orientation selectivity. It has become increasingly clear over recent years that cortical properties depend heavily on the response characteristics of subcortical channels [10,11] and the way in which the cortex combines its subcortical inputs [12]. We therefore start by describing the spatial distribution of subcortical channels in the model and the process by which visual development weights the cortical inputs. Fig 2A shows all subcortical channels in a 6˚×6˚patch of visual field. Each channel is represented by a circle at the centre of its receptive field, red for on-centre channels and blue for offcentre. The calculation of the locations ensured that closest neighbours were almost always of opposite polarity, as required by anatomical measurements [9]. The weighting of each channel's synapse with a cortical neuron is represented by circle diameter, and all weights are assumed equal at the start of the development process.
Response-dependent development of the visual system is driven by intrinsic connections [23] and at least two stimulus sources: waves of activity traversing the retina [24], and moving visual stimuli encountered after eye opening [25]. Neuronal responses in the model were therefore stimulated with a drifting sinusoidal grating over the full range of orientations. Each cycle in the development process then consisted of increasing the weight of all geniculocortical synapses for one randomly chosen subcortical channel. If the action potential rate of a cortical neuron increased as a result, the synapse between the channel and the cortical target remained strengthened. Otherwise, the synaptic weight was weakened to less than its original value. Fig 2B and 2C show the weights at intermediate steps of visual development for the excitatory cortical neuron at the middle of the visual field patch, and Fig 2D shows the final result. On-centre channels dominate off-centre channels in oval areas of visual field and off-centre channels dominate in other areas. The source of the on/off segregation is easy to understand. Neighbouring channels tend to have opposite sign so that their signals typically cancel each other when they sum at a cortical neuron. A strengthening of one channel and weakening of its opposite-sign neighbour will increase the cortical response, and the Hebbian development process preserves this bias. The channel sign that dominates a visual area will depend on the randomised process by which channel locations were assigned. Fig 3A and 3B show the start and end of development for the neuron in Fig 2. Fig 3C shows synapse strength for a cortical neuron located at the position marked by the yellow dot. The on-off segregation is in much the same direction as that in part B of the figure, but the pattern is more odd-than even-symmetric. Part D of the figure shows synaptic strengths for a cortical neuron in a third location. In this case the orientation of on-off segregation differs from the other two locations, indicating that each cortical neuron has its own spatial pattern of inputs. All of these on-off segregations are suggestive of the receptive field patterns to which they could contribute: we will show receptive fields at the end of the results.

Intracortical inhibition
Hebbian changes could potentially increase a synapse's strength beyond the physiological limit. But the model includes intracortical inhibition which is driven by the same geniculocortical input as are the excitatory neurons. As synapses increase in strength so does inhibition: this limits excitatory responses in the cortex, preventing further synaptic strengthening. We now demonstrate this growth in inhibition during development and its effect on response time courses. Fig 4A shows the time course of two representative cortical neurons at the start of visual development. Responses to a 2 Hz drifting grating are shown. The blue curve indicates the somal generator potential of a geniculorecipient inhibitory neuron, the orange curve shows potential in the axonal terminal in the same neuron, and the red curve gives potential in an excitatory neuron with which the axon makes an inhibitory synapse. All response amplitudes  The generator potential time course is shown for two cortical neurons-inhibitory and excitatory-responding to a drifting grating. Both neurons are located at the centre of the visual field patch, and both compartments of the inhibitory neuron are represented. The near-cancellation of activity in neighbouring on-and offchannels results in low-amplitude responses in both neurons at the start of development (A), and the functional segregation of on-and off-channels produce larger response amplitudes post-development (B). C. Action potential rate in the excitatory neuron is shown for four steps in the development process, demonstrating the iceberg effect. D. Comparison of the shape of the time course with that recorded in real cortex by Dean and Tolhurst [26]. The horizontal axis shows the fundamental Fourier amplitude of impulse rate divided by mean rate, and the vertical axis gives the fraction of excitatory cells in the model (blue) and for simple cells (red). are low due to the near-cancellation of neighbouring on-and off-channels. Fig 4B, obtained at the end of development, shows much larger amplitudes due to the functional segregation of on-and off-channels. The figure also shows that the mean potential of the inhibitory terminal has risen due to rectification at the axon initial segment, and that the excitatory cell's mean potential has fallen as a consequence.
The result is the iceberg effect [17] illustrated in Fig 4C. This shows action potential rate in the excitatory neuron for four steps in the development process. As the excitatory drive to this neuron rises so does inhibition, and only the peak of the underlying response is seen. We used Fourier analysis to compare the shape of this time course with that found in the published literature. The horizontal axis in Fig 4D shows relative modulation, defined as the ratio of the Fourier fundamental component of the impulse rate divided by mean rate. Dean and Tolhurst [26] showed that as the response becomes more peaked, and the peristimulus time histogram is reduced to a single bin, this ratio approaches 2. A frequency histogram of the ratio is shown in Fig 4D for

Orientation tuning
To measure orientation tuning in the model we drifted a sinusoidal grating across the visual field at a variety of orientations. Fig 5A gives the result for the cortical neuron whose geniculocortical strengths are shown in Fig 3B. Early in the development process, response amplitudes are low and orientation tuning is weak. Later in development two tuning peaks appear, one for motion in one direction and the other for motion in the opposite direction. The preferred orientation matches that expected from the synaptic strengths. To check tuning precision, we fitted a curve to the response at each cortical location: Fig 5B shows the result of fitting the curve to the post-development data in part A of the figure. Tuning bandwidth was measured as the half width at half height of the fitted curve, as shown. Fig 5C shows tuning bandwidths as a function of peak response amplitude at four steps of development, where each circle represents one excitatory cortical neuron. Not surprisingly, bandwidths are large and scattered at the start of development. Bandwidths at the end of development, however, are tightly clustered and as low as 15˚. To compare this result with real cortex we calculated the frequency histogram of bandwidth, as shown in Fig 5D. The histogram falls between two recent empirical estimates [27,28], indicating that bandwidth estimates in the model are realistic.

Orientation map
Previous work, particularly optical imaging, has shown that preferred orientation forms characteristic patterns across the cortical surface [29]. We next wished to determine whether our model reproduces such patterns. Preferred orientation, taken from the maximum of the motion direction tuning curve, is shown as a function of visual field location in Fig 6. Part A of the figure shows that an orientation map exists even before development starts. This is as expected from the orientation tuning curve in Fig 5A, development cycle 0, which shows a mature preferred orientation even though tuning is poor. Like the empirical work, the map at the end of development ( Fig 6B) displays regions in which orientation changes little and other regions containing pinwheels, around which orientation varies across its whole range. The consistency of the map during development matches findings in ferret primary visual cortex [16]. Orientation tuning curves for the excitatory neuron whose weights are shown in Fig 3B. Four cycles in the development process are given. At the beginning of development, orientation tuning is poor. Later in development two tuning peaks arise, one due to grating motion in one direction and the other for the opposite direction. B. A tuning curve was fitted to the post-development data in part A. The model comprised the sum of two von Mises functions separated by 180˚. C. Tuning bandwidth, measured as half width at half height, was obtained from the fitted curve in part B. Bandwidth is shown as a function of maximum response amplitude for each neuron in the central 2˚×2˚of the visual field patch; each circle represents one neuron. Bandwidths range from small to large at the start of development but are uniformly low post-development. D. Frequency histogram of bandwidths comparing values in the model with those found in the cortex by Koch et al. [27] and Nauhaus et al. [28]. Bandwidth estimates in the model fall within these empirical findings. The model provides us with a unique opportunity for exploring the source of these patterns. The first clue comes from comparing the pattern at the start and end of development (Fig 6A  and 6B). The patterns are very similar (correlation coefficient = 0.97, n = 6559, p<round-off error, null hypothesis: maps are uncorrelated), indicating that the basic pattern is set prior to any developmental changes in the cortex. This suggests that the orientation maps are determined by the spatial layout of subcortical channels. We tested this idea by subtracting the map of off-centre channels from that for the on-channels. The result, in Fig 6C, shows areas clearly dominated by one or the other contrast polarity. To see whether this inhomogeneity could produce the orientation map we took the dot product of the map in C with Gabor functions varying in both spatial phase and orientation. The details of this calculation are provided in the Methods section. The result, in part D of the figure, shows clear similarities with the orientation map (correlation coefficient = 0.88, n = 6559, p<round-off error). In a further six simulations, obtained by varying the randomisation seed (and therefore the retinal ganglion cell array), the correlation coefficient was never less than 0.87. It seems, therefore, that clumping of on-centre channels in one region and of off-channels in a nearby region can produce orientation preference maps that look much like those recorded in the laboratory.

Orientation map periodicity
We compared the statistics of the map with published work by calculating its periodicity. Each orientation in Fig 6B was converted to a unit vector and a Fourier transform was calculated for the real and imaginary parts of the vector map, as shown in Fig 7. The peak magnitudes of both transforms were at 0.63 cycles/deg from the origin, yielding a periodicity of 1.6 degrees of visual angle. Converting to distance, using the cortical magnification factor calculated in the Methods section, yielded a periodicity of 1.1 mm. This compares well with the estimates of Löwel et al. [30] and Diao et al. [31] who found values of 1-1.1 and 1.1 mm, respectively. We also counted pinwheels, as defined in the Methods. The result was 2.8 pinwheels per orientation hypercolumn, close to the mean value (3.1) found in three species by Kaschube et al. [32] and in four species, including the cat, by Schottdorf et al. [33]. We are therefore confident that the orientation preference map obtained from the model faithfully reproduces empirical data.
What is the source of the periodicity in the orientation map? Given that the map can be predicted from subcortical arrays, we looked at subcortical sources. Fig 8A shows an analysis like that in the previous figure except that here the analysed map is subcortical. Pulses representing off-centre ganglion cell locations were subtracted from on-centre pulses, the map was smoothed with the spatial profile of the geniculate centre mechanism, and the result was Fourier transformed. The maximum magnitude in this map is found at 0.60±0.05 cycles/deg from the origin (mean±standard error, obtained from the standard and six other maps), representing a periodicity of 1.7˚. This value is close to that found for the orientation map (1.6˚), indicating that we should look subcortically for the source of the periodicity.
One possibility is aliasing [34]: the on-and off-centre ganglion cells in our model have a density of 24.4 and 26.6 cells/deg 2 , respectively. As described in the Methods, this results in a periodicity of 6.5˚which is about four times the required value. Another possibility is shown in Fig 8B, the upper part of which shows the one-dimensional centre mechanism profiles of neighbouring on-and off-centre cells. The lower part of Fig 8B shows the difference between impulses representing off-centre subcortical channels from those representing on-centre channels. The map has been smoothed with a Gaussian profile with standard deviation 0.1˚. There are areas that are clearly dominated by one or the other sign of response. D. Orientation map produced from taking the dot product of the map in part C with Gabor functions varying in both spatial phase and orientation. The Gabor functions had a standard deviation of 0.7˚. There are strong similarities with the map in part B, indicating that the orientation preference map can be largely predicted by the spatial layout of subcortical channels.
https://doi.org/10.1371/journal.pcbi.1007254.g006 these two profiles to mimic the results of cortical convergence. The Fourier transform of this difference curve, shown in Fig 8C, has a peak 0.56 cycles/deg from the origin, which translates to a periodicity of 1.8˚. This value is close to the periodicity of 1.6˚found in the orientation map, strongly suggesting that the latter value is set by the differencing of on-and off-centre receptive field profiles.

Receptive fields
The geniculocortical weight maps in Fig 3 are reminiscent of simple cell receptive fields. We calculated receptive fields in the model with a sparse noise stimulus: squares of light or dark, as shown in Fig 9A, were briefly presented at a variety of visual field locations. Responses in excitatory cortical neurons were calculated and the peak impulse rate was recorded at each location. Contour plots of the responses to light and dark are shown in red and blue, respectively. Fig  9B, 9C and 9D show the resulting receptive fields for the neurons whose synaptic weights are given in Fig 3B, 3C and 3D, respectively.
Stimulus locations far from the neuron produce a generator potential that fails to reach spike threshold, and the receptive fields are therefore smaller than the corresponding weight maps. In particular, the weight map in Fig 3B is nearly even-symmetric so that only the central subfield survives in the receptive field. Otherwise the synaptic weight patterns are faithfully represented in their receptive fields. Another feature of interest is the peak impulse rate shown next to each subfield. Responses to dark stimuli are clearly greater than to light stimuli with the same contrast magnitude, reflecting the dark dominance seen in real cortex [35,36]. Dark dominance in the model arises from two sources: off-centre ganglion cells outnumber on-  The periodicity of the subcortical channel map was calculated by subtracting pulses representing off-centre cell locations from those representing oncells, smoothing with the geniculate centre mechanism spatial profile, and Fourier transforming the result. B. The periodicity was also calculated by comparing on-and off-centre pairs. The upper panel shows the centre mechanism profile for neighbouring on-and off-centre geniculate cells, and the lower panel shows centre cells, and off-centre geniculate responses are faster than their on-centre counterparts. Both of these asymmetries reflect empirical findings [9,37].

Discussion
In this paper we have described a visual system model supporting the following conclusions.
• A Hebbian process is sufficient to functionally segregate on-centre and off-centre inputs to primary visual cortex.
• This segregation can produce orientation-selective neurons.
• The resulting selectivity has a precision mirroring that of real cortical neurons.
• The cortical map of orientation preference arises not from cortical sources but from local clustering of on-centre and off-centre neurons in the retina.
This last conclusion leads to two predictions for future experiments. First, measurement of the locations of on-and off-centre β ganglion cells in the cat retina will allow the calculation of the preferred orientations in the corresponding region of primary visual cortex. An experiment testing this prediction would be difficult to perform because it would require simultaneous measurement of retinal arrays and preferred orientation in the cortex. We therefore offer a second prediction, which should be rather more straightforward to test: that the periodicity of the orientation map depends on the spatial profile of the geniculate relay cell's centre mechanism. This result follows from the analysis illustrated in Fig 8. We now discuss some of the issues arising from our results.

Previous models
A number of previous models have addressed issues such as orientation selectivity and cortical mapping. How do our results fit in with this previous work? One of the earliest studies was by von der Malsburg [38], whose model consisted of a retinal layer connected directly to a cortical layer. Retinocortical connections varied in strength through Hebbian plasticity, producing a cortical map of preferred orientation. The retina, however, consisted entirely of on-cells, limiting the orientation selectivity of cortical cells. Miller [39] described a model for the development of orientation selectivity, and showed that it could produce both receptive fields like those in simple cells, and cortical maps of orientation preference resembling laboratory observations. To do this Miller assumed that nearby like-sign subcortical neurons have positively correlated activity but that the correlation between more distant neurons is between those of opposite sign. We have not made any assumption about the statistics of subcortical activity. Instead, we have used the known statistics of subcortical receptive field locations [9] to calculate orientation selectivity: simple cell-like receptive fields result from spatial clustering of same-sign subcortical channels rather than from response correlations.
Somers, Nelson and Sur [40] described a model in which subcortical inputs provide a mild orientation bias to cortical neurons, but the authors do not provide a mechanism by which such a bias might develop. Instead, they showed that the orientation tuning can be sharpened by connections between excitatory cortical cells, leading to positive feedback. These connections correspond to the prolific excitatory-to-excitatory synapses found in anatomical studies [41]. A more recent study [27] has also shown that mutual excitation can improve tuning their difference. Profile maxima are normalised. C. The on-off difference was Fourier transformed: the normalised magnitude of the transform is shown. The peak of the transform is at 0.56 cycles/deg which converts to a periodicity of 1.8˚in the subcortical map.
https://doi.org/10.1371/journal.pcbi.1007254.g008 Fig 9. A. Receptive fields were calculated using a sparse noise stimulus. The stimulus was a 1˚×1˚square, shown to scale, presented at the nodes of a square grid with node spacings of 0.25˚. The patch had a duration of 50 ms and a contrast of ±1. B, C, D. The maximum impulse rate at each stimulus location was compiled into contours. Responses to lights and darks are shown in red and blue, respectively. Peak impulse rate is shown next to each subfield and the contours give 0.05, 0.5 and 0.95 of this peak rate. The cortical cells represented in B, C and D are the same as those whose synaptic weights are shown in Fig 3B, 3C and 3D Origin and development of orientation selectivity sharpness. We have chosen to avoid recurrent connections from excitatory cells for simplicity and ease of interpretation.
One of the advantages of the Somers et al. [40] model is that the units it uses, such as micrometres of cortical surface, correspond with those used in empirical studies. We have taken their approach even further by specifying spatiotemporal and response parameters with units of degrees of visual angle, seconds, millivolts and impulses/s. This aids comparison of the model with experimental observations. Another advantage of our model is that its stimulus can be any spatiotemporal pattern, provided that it is monocular and monochromatic. The model can therefore be tested with stimuli not described here.

Development
Previous work has shown that orientation preference maps are firmly established early in visual life, even though orientation tuning is weak. Chapman et al. [16] used optical imaging to measure maps in ferret cortex around the time of eye opening. They found that maps in each animal were geometrically stable over many weeks even though responses were weak in the earliest recordings. Our model provides a basis for this finding by proposing that the stability originates in the retinal ganglion cell array. Crair et al. [42] recorded orientation preference maps in both normal and binocularly deprived cats and found that the maps were present in animals without visual experience. This suggests that the stimuli we have used to drive Hebbian development are due to moving waves of subcortical activity rather than motion in the visual world.
Later work, however, has introduced a new factor into this parsimonious picture. Smith et al. [23] used calcium imaging of spontaneous activity in ferret visual cortex to show correlated activity between widespread cortical patches. They found that correlated patches had similar orientation preferences when stimulated. The patchy correlations were evident before eye opening, at a time when orientation preference maps were still forming. How the subcortically driven maps that we have simulated fit in with the long-range cortical networks found by Smith et al. remains to be determined.

Orientation selectivity
Previous models have assumed that orientation selectivity depends on discrete lines of on-and off-centre cortical inputs [5,13] or on small numbers of nearest-neighbour retinal ganglion cells of opposite sign [10,11]. Both of these models now look impractical. Consider, for example, a visual field point that is 11˚from the cat's central area, as in the model we have described. According to the estimates in the Methods section, the density of X-type retinal ganglion cells at the corresponding retinal location is about 51 neurons/deg 2 [43]. Given that each X-cell drives one or very few relay neurons in the lateral geniculate nucleus [44,45], the densities of geniculocortical inputs will be similar. Those inputs produce a simple cell receptive field with a median area of 2.7 deg 2 [46]. We can therefore expect that about 51×2.7 = 140 subcortical inputs converge on a simple cell. Reid and Alonso [7] estimated that a third of those inputs are functionally effective, while our estimate (Fig 3B, which shows that one or other input sign dominates at each point) is closer to one half. We therefore estimate that a simple cell is driven by 47 to 70 subcortical neurons, substantially higher than in the earlier models.
A recent paper [47] modelled the contribution of retinal on/off pairs to cortical orientation preference maps. The authors concluded that paired interactions do not lead to realistic orientation preference maps. Our results provide an alternative model, namely that orientation preference arises from localised high densities of same-sign retinal neurons. In particular, orientation preference at a specific visual field location is determined by a retinal area in which on-centre neurons are more densely packed than are off-centre neurons, and a nearby area in which off-centre neurons are more densely packed (Fig 6C). Whether our model reflects the seeding of real orientation preference maps remains to be tested in the laboratory.

Random variation
There are two sources of randomness in the model, both structural. The first is the random variation of a channel location about its node on a rectangular grid. This variation seeds the orientation map. The second source of random variation is the sequence in which channels are tested with an increase in geniculocortical strength. Changing this sequence has little effect on the results. We have chosen not to add noise to the membrane potential: all of the differential equations defining the model are deterministic. This choice comes with two disadvantages. First, two studies [48,49] have shown that membrane potential noise helps to preserve the precision of orientation tuning as contrast increases. Contrast invariance in our model may suffer from the lack of membrane potential noise. Second, the frequency histograms generated by the model are substantially narrower than those recorded in the laboratory (Figs 4D and 5D). Membrane potential noise would broaden the model histograms. We have chosen, however, to exclude such noise from the model in order to make it simpler to interpret.
We calculated impulse rate by thresholding membrane potential without considering random fluctuations in potential. There is strong empirical support for this approach. Carandini and Ferster [50] recorded membrane potential in primary visual cortical cells, low-pass filtered the potential and then applied a half-wave rectifier to predict impulse rate. The prediction worked well over a wide variety of conditions, including adaptation and a range of contrasts. Our model uses the relationship they measured between membrane potential and impulse rate, as shown in Eq 2.

Resting activity
Resting activity is a key component in our model. In keeping with subcortical measurements [51] we assume that X-cells have a substantial resting impulse rate, which helps to linearise their responses to stimuli of moderate contrast. This resting activity translates to an even higher impulse rate in the cortical inhibitory neurons, which corresponds with empirical findings [21]. In turn, the high firing rate in inhibitory neurons produces a resting hyperpolarisation in the cortical excitatory neurons. This matches the absolute contrast threshold found in simple cells [52] and the hyperpolarisation in intracellularly recorded simple cells [53]. All of these resting activities in the model result from a single parameter, the constant depolarisation added to the ganglion cell generator potential (Eq 8). While resting activity in the real visual system may derive from multiple sources, it is interesting that we have been able to construct an internally consistent model in which resting activities result from a single subcortical source.

Model equations in the time domain
Here we derive the equations describing the model. Each neuron is represented by a single nonlinear differential equation, and time courses are obtained by simultaneous numerical integration of the equations for all neurons. The difference between membrane potential at the initial segment of a neuron's axon and action potential threshold determines the action potential rate. This difference is therefore called the generator potential, denoted by p(t) where t is time. The input to the neuron is a set of synaptic potentials, v i (t), each of which is weighted by a gain, g i . The model assumes that the neuron is a low-pass filter that integrates the difference between the driving potential and generator potential: where τ is the time constant. The neuron's action potential rate, a(t), is obtained by rectifying the generator potential: ( To complete Eq 1 we need to know how the activity in one processing stage, z, depends on that in the previous stage, z−1. We assume that the postsynaptic potential, v i , is proportional to presynaptic impulse rate. The general equation for a model neuron is then: where proportionality constants have been absorbed into gain g i . The general equation requires modification for each stage of the model. The stage numbers are 1 to 7 representing, in order, photoreceptors, bipolar cells, ganglion cells, geniculate neurons, inhibitory neuron somas, inhibitory neuron axons, and excitatory cells. Fig 1 shows the signal-processing sequence. The photoreceptors receive their input from the visual stimulus rather than a presynaptic neuron, and do not produce action potentials: t dp j1 ðtÞ dt ¼ À g js ðx; yÞ � sðt; x; yÞ À p j1 ðtÞ ð4Þ where p represents the difference between membrane potential and resting potential rather than generator potential, j is channel number, the subscript s indicates subcortex, x and y give visual field location and s(t,x,y) is the stimulus. We use the (�) symbol to represent both vector and integral dot products. The gain is a Gaussian function of location representing subcortical spatial spread due to optical blurring and neural convergence: where g s and r s are the contrast sensitivity and radius of the centre mechanism (the model does not include a surround mechanism) and (x j , y j ) is the spatial location of channel j. The dot symbol in Eq 4 represents the dot product: The minus sign preceding the dot product in Eq 4 corresponds to photoreceptor hyperpolarisation by light. Bipolar cells do not produce action potentials and can be on-or off-centre: t j dp j2 ðtÞ dt ¼ À n j p j1 ðtÞ À p j2 ðtÞ ð7Þ where τ j and n j are the time constant and sign for channel j. Ganglion cells produce action potentials so p represents generator potential for these and subsequent neurons. Ganglion cells also have a resting impulse rate, which is implemented in the model by adding a constant depolarisation, p s , to the driving potential: t j dp j3 ðtÞ dt ¼ p j2 ðtÞ þ p s À p j3 ðtÞ ð8Þ Relay cells in the dorsal lateral geniculate nucleus inherit the constant depolarisation from, and rectify, their input: t j dp j4 ðtÞ dt ¼ hðp j3 ðtÞÞ À p j4 ðtÞ ð9Þ The input to inhibitory cortical neuron k is obtained by taking the dot product of the subcortical input with a Gaussian convergence function: t dp k5 ðtÞ dt ¼ X j g kc ðx j ; y j Þw jk hðp j4 ðtÞÞ À p k5 ðtÞ ð10Þ where the convergence function g kc is obtained by subscript substitution into Eq 5, c stands for cortex, and w jk is the strength of the synapse from subcortical input j. The generator potential in the inhibitory neuron initial segment is rectified and integrated by its axon and connection into the inhibitory network: Finally, excitatory cortical neuron k is driven by the sum of its subcortical and inhibitory inputs: t dp k7 ðtÞ dt ¼ X j g kc ðx j ; y j Þw jk h p j4 ðtÞ � � À X l g ke ðx l ; y l Þp l6 ðtÞ À p k7 ðtÞ ð12Þ where e stands for excitatory neuron. Eqs 4-12 together define the model.

Model equations in the frequency domain
We transformed the model into the frequency domain for two reasons. First, we reduced the possibility of mathematical and computational errors by ensuring that the solutions in the temporal and frequency domains agreed to within round-off error. Second, we reduced computation time by performing most of the calculations in the frequency domain. The Fourier transform of Eq 4 is: itoP j1 ðoÞ ¼ À g js ðx; yÞ � Sðo; x; yÞ À P j1 ðoÞ ð13Þ where i ¼ ffi ffi ffi ffi ffi ffiffi À 1 p , ω is temporal frequency, and Fourier transforms are shown in upper case. Thus: Similarly, the transforms for the following subcortical stages are: where δ is the Dirac delta function and F is the Fourier transform. At the cortical level:

Solution for drifting grating
Most of the simulations in the paper use a drifting grating as stimulus. We made these simulations faster by using the analytical solution for the dot product, Eq 6; the solution follows. There is no surround antagonism in the basic model. To compensate, stimuli are defined in terms of contrast rather than luminance. Contrast is obtained by finding the difference between local and background luminance, and dividing the difference by background luminance. The equation for a drifting grating is: where c is the contrast, ψ stim is the spatial frequency, ω stim is the temporal frequency, and u = cos(θ)x+sin(θ)y is the distance in the grating's direction of motion, θ, with u j being the location of the jth channel. The dot product is then: (where v is distance perpendicular to u) which has transform g js ðx; yÞ � Sðo; x; yÞ ¼ pg s c exp À r 2 s c stim

Neuronal location
A neuron's location is defined in the model by the centre of the convergence function that weights its inputs. Off-centre subcortical channels were located at the nodes of a square grid aligned with the visual field patch, with a node at the centre of the patch. Each location was then perturbed with a Gaussian deviate in both the horizontal and vertical directions. Jang and Paik [54] provided evidence for developmental repulsion between on-and off-centre ganglion cells. Accordingly, on-centre channels were distributed similarly to off-centre channels except that the grid was offset: the four nodes closest to the centre of the visual field patch were equidistant from the centre. Cortical neurons were located on an unperturbed grid aligned with the patch, with the central node at the centre of the patch. Both an inhibitory and excitatory neuron were located at each grid node. The model used a retinal ganglion cell map based on the work of Wässle et al. [9]. Fig 10  illustrates two analyses designed to test whether the map we used accurately reproduces the statistics of measured maps. Part A shows the frequency histogram for the distance between nearest on-and off-centre neighbours on the left and right, respectively. Wässle et al.'s data are shown in red and the model in blue. The widths of the measured and modelled histograms are similar, as required. Part B provides a test for whether the modelled on-and off-centre maps are statistically independent, as defined by Rodieck [55]. Each histogram shows the density of on-centre cells in annuli of the stated distance from reference off-centre cells. Histograms are provided for two random seeds used to generate the maps. For each seed, a linear regression performed on the data showed that there was no linear trend (seed 1: F(1,18) = 1.12,p = .30; seed 2: F(1,18) = 0.23,p = .64). Our modelled maps are therefore consistent with the conclusion of Eglen et al. [56] that on-and off-centre maps are independent outside a small inner area.

Development
The development process adjusted the strength of the synapse of each geniculate neuron onto each of its cortical targets. At the start of development all of these synapses were assigned a weight of 1. For each development cycle a geniculate neuron was selected, with all neurons equally likely to be chosen. All synaptic weights for this neuron were increased by 0.2 and the model was stimulated with gratings drifting in 16 directions evenly distributed across the whole range. Each excitatory cortical neuron's impulse rate was calculated, and if the maximum response increased relative to the previous cycle, the weight increase was retained for that neuron. Otherwise, the weight was reduced by 0.2 relative to its value on the previous cycle. Weights were restricted to lie between 0 and 2.
The number of development cycles was determined as follows. Each geniculocortical synapse needed five cycles to change from its starting value to its minimum or maximum. The number of cycles was therefore set at five times the number of geniculate neurons, 5×3281 = 16405, and rounded to the nearest thousand, 16000. Figs 2 and 4C show that synaptic strength and response amplitude, respectively, change very little between cycles 10666 and 16000. Indeed, response amplitude changed by at most 2% between these two cycles. We therefore refer to cycle 16000 as "After development". Previous work has shown that inhibitory connections mature during visual development [57]. We therefore increased the strength of the synapse of inhibitory neurons onto excitatory neurons during development. All of these synapses had equal strength, which increased linearly with development cycle from 1 to g e . Table 1 provides a glossary of model parameters and their values. The following text describes these parameters and explains how they were set.

Spatial parameters
Location and size of visual field patch. We simulated a visual field patch centred on the horizontal meridian, and 11˚from the central area. The size, 8˚×8˚, is substantially larger than a typical cortical receptive field. Concentration of retinal ganglion cells. The mean density of β ganglion cells at 11e ccentricity is 1275 cells/mm 2 [43]. Given that β cells are the morphological correlates of Xtype ganglion cells, the mean concentration of X-cells is then 1275×(0.2) 2 = 51 cells/deg 2 . Wässle et al. [9] counted 65 on-and 71 off-centre cells in their analysed sample (their Fig 6). The concentration of X-type on-centre cells is therefore (65/(65+71))×51 = 24.4 cells/deg 2 and, similarly, 26.6 for off-centre cells.
Subcortical channel location. Wässle et al. found that the packing of same-sign β cells ranged between square and hexagonal arrays; we used a square grid for simplicity. They measured the distance between same-sign nearest neighbours and found that the standard deviation divided by distance was 0.189. We therefore placed on-channels on a square grid with 1= ffi ffi ffi ffi ffi ffi ffi ffi ffi 24:4 p ¼ 0:20 � spacing and then perturbed the locations with a 0.20×0.189 = 0.038˚Gaussian deviate. Similarly, off-channels were placed on a 0.19˚±0.037˚grid. Centre mechanism size. Saul and Humphrey [59] measured the centre size of X-type geniculate cells at a mean eccentricity of 11˚, which is the reason we chose this eccentricity for the visual field patch. From their mean radius of the centre mechanism of (non-lagged) Xcells, r s = 0.40˚.
Cortical magnification factor. We use the value measured by Tusa et al. [60] at 11e ccentricity along the horizontal meridian to obtain the cortical magnification factor: 0.45 mm 2 /deg 2 .
Radius of geniculocortical convergence. This radius is derived from the work of Jones and Palmer [46], who modelled simple cell receptive fields as Gabor functions. Their Table 1 provides the width of the Gaussian functions they used to compute the Gabors. Converting from standard deviation to radius and taking the median across their sample yields the value 1.03˚. This value represents the dot product of the subcortical centre mechanism with geniculocortical convergence. Deconvolution gives r c = 0.95˚.
Location of cortical neurons. The critical issue in choosing cortical cell density is that it be substantially less than the radius of the geniculocortical convergence. We chose to situate both excitatory and inhibitory neurons on a square grid with 0.1˚spacing.
Radius of inhibitory-excitatory convergence. We have not found a measurement of inhibitory receptive field radius as rigorous as the Jones and Palmer analysis of simple cells. We have therefore set the radius equal to r c .

Temporal parameters
Time constants. Each subcortical channel consists of a cascade of four first-order lowpass filters. The impulse response of a series of z low-pass filters with time constant τ peaks at (z−1)τ. A cortical cell receiving this input is the fifth cell in the cascade and therefore has a peak time of 4τ. Since simple cell impulse responses peak at values as low as 40 ms [22], τ = 40/ 4 = 10 ms.
Subcortical time constants. It has recently been shown that off-centre X-type geniculate cells lead their on-centre counterparts. The leading edge of the impulse response in off-cells precedes that in on-cells by a mean of 3 ms when measured at 40% of maximum response [37]. We replicated this result by setting time constants in stages 2, 3 and 4 of the subcortical channels as follows: τ j = 11 ms, on-channel, τ j = 9 ms, off-channel.
Inhibitory time constant. Inhibitory strength is an important variable in our model, but setting this strength involves a trade-off between the inhibitory time constant and the inhibitory-to-excitatory gain. We have chosen to set the inhibitory time constant to a relatively large value and use empirical evidence to set inhibitory-excitatory gain as described below. The time constant, 200 ms, is large enough to match the long-lasting inhibitory tail seen in simple cells responding to flashed stimuli [22].

Intensive parameters
Generator gain. The form of the generator function and its gradient, g rect = 7.2 Hz/mV, are taken directly from the work of Carandini and Ferster [50].
Geniculate contrast sensitivity. This parameter is calculated by integrating the centre mechanism's spatial profile over both dimensions: Z Z g s pr 2 We set this equal to the contrast sensitivity of the X-type ganglion cell centre mechanism, 620 Hz/contrast-unit (from the 2 Hz data in Fig 12 of Frishman et al. [61]), multiplied by the attenuation between retina and geniculate, 0.73 (from Fig 5A of Kaplan et al. [51]). Finally, converting from Hz to mV, g s is given by: g s g rect ¼ 620 � 0:73 ¼ 450 Hz=contrastÀ unit Subcortical resting depolarisation. The resting impulse rate of geniculate cells averages 14 Hz [51]. Thus: p s ¼ 14 Hz=g rect ¼ 1:9 mV Cortical contrast sensitivity. The contrast sensitivity of cortical neurons is best calculated from empirical measurements of cortical membrane potential, which avoid the complications of action potential thresholding. Carandini and Ferster [50] measured the membrane potential of simple cells stimulated with gratings of optimal orientation and spatial frequency. For the three cells in their Fig 13, response amplitude divided by contrast has a maximum gradient that averaged 70 mV/contrast-unit. The geniculocortical gain g c was set so that model contrast sensitivity was close to this value.
Inhibitory-excitatory gain. This gain sets the resting hyperpolarisation in excitatory cells. Anderson et al. [53] measured this hyperpolarisation in nine simple cells and found a median difference of 9 mV between the threshold and resting potential. Gain g e was set so that the resting hyperpolarisation in excitatory cells approximated this value. Fig 6D shows the orientation preference map calculated from the spatial map of subcortical channels. The calculation used the following steps.

Predicting the orientation map
1. The subcortical map was represented by a grid with elements fine enough (0.005˚×0.005˚) that on-and off-channels did not coincide.
2. On-and off-channels locations were assigned the grid values 1 and -1, respectively, and all other locations zero. Fig 6C shows a smoothed version of this location map.
3. Gabor functions were constructed with a standard deviation of 0.7˚, the standard spatial frequency (0.5 cycles/deg), 8 orientations uniformly distributed across the range 0 to 180˚, and 8 phases uniformly distributed across the range 0 to 360˚.
4. The dot product of the channel grid and each of the (stationary) Gabors was calculated.
5. For each grid location the maximum value of the product across orientations and phases was determined.
6. The orientation that yielded the maximum value at each grid point was used as the preferred orientation.

Analysing the orientation map
Periodicity. The periodicity of the orientation preference map shown in Fig 6B was calculated in several steps. First, the orientation preference at each location was converted to a vector with unit length. Second, each vector was replaced by its real part and the resulting map was cropped to 5˚×5˚to avoid edge effects. Third, a two-dimensional Fourier transform was computed; its magnitude is shown in Fig 7A, with the component at zero spatial frequency set to zero. In the same way, a transform was also calculated for the vector imaginary part and is