{'fraction_non_alphanumeric': 0.031365653178481605, 'fraction_numerical': 0.013576788598769093, 'mean_word_length': 4.512623985572588, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Laminar model for the joint development of ocular dominance columns and CO blobs in V1

In this paper, we present a multi--layer, activity--dependent model for the joint development of ocular dominance (OD) columns and cytochrome oxidase (CO) blobs in primate primary visual cortex (V1). For simplicity, we focus on layers 4C and 2/3 with both layers receiving direct thalamic inputs and layer 4C sending vertical projections to layer 2/3. Both the thalamic and the vertical connections are taken to be modifiable by activity. Using a correlation--based Hebbian learning rule with subtractive normalization, we show how the formation of an OD map in layer 4C is inherited by layer 2/3 via the vertical projections. Competition between these feedforward projections and the direct thalamic input to layer 2/3 then results in the formation of CO blobs superimposed upon the ocular dominance map. The spacing of the OD columns is determined by the spatial profile of the intralaminar connections within layer 4, while the spacing of CO blobs depends both on the width of the OD columns inherited from layer 4 and the spatial distribution of intralaminar connections within the superficial layer. The resulting CO blob distribution is shown to be consistent with experimental data. In addition, we numerically simulate monocular deprivation and find that while the CO blob distribution is unaltered, the OD pattern undergoes modification. The OD stripes of the deprived eye narrow, whereas the OD stripes for the remaining open eye widen.


Introduction
Most neurons in the primary visual cortex (V1) respond preferentially to a variety of features associated with visual stimuli, including orientation, ocular dominance (OD), spatial frequency, and direction selectivity. As one progresses tangentially across V1, the response properties vary in a nearly continuous fashion. Seminal studies by Hubel & Wiesel [1,2] revealed the existence of ocularity and orientation preference feature maps. Furthermore, they found that neurons through the cortical layers have similar response properties suggesting that the visual cortex is organized in a columnar fashion. It was further conjectured that the feature preferences of cortical neurons are generated by the convergence of thalamic afferents on to input layer 4, and then passed on to other layers through vertical interlaminar projections.
It follows from this that the formation of feature preference maps can be understood in terms of the development of feedforward connections from thalamus to layer 4. The latter has been the focus of most models of cortical development [3,4].
Although the columnar hypothesis proposed by Hubel and Wiesel may be approximately true in cat, it is clearly an oversimplification in primates such as the macaque, where feature preferences are not uniform throughout the depth of a column [5,6]. For example, there are two major classes of thalamic cells projecting to input layer 4C of the macaque monkey, magnocellular (M) and parvocellular (P). The M and P pathways connect to layer 4C in a graded fashion with the M pathway dominating upper 4Cα and the P pathway dominating lower 4Cβ [7,8].
Interestingly, orientation-selective cells within layer 4C also appear to emerge in a graded fashion, being found predominantly in mid and upper layer 4Cα [9,10], and coinciding with the emergence of horizontal axon projections from excitatory spiny stellate cells [11,12]. It is possible that such projections, rather than direct thalamic inputs, provide the anatomical substrate for the generation of orientation preference [13]. A second example arises from the observation that direction-selective cells first emerge in layer 4B but are not found in superficial layers 2/3 [14].
Yet another important example of nonuniformity through the layers is the occurrence of cytochrome oxidase (CO) blobs in superficial layers of macaque V1 [15,16].
These regions of higher metabolic activity receive direct inputs from a third class of thalamic cells, namely those forming the koniocellular (K) pathway [17,18]. The spatial distribution of CO blobs within cortex is correlated with a number of feature preference maps. For example, the blobs are found at evenly spaced intervals along the center of OD columns [16], and neurons within the blobs tend to be less binocular and less orientation selective [19]. The blobs are also linked with low spatial frequency domains [20] and appear to be coincident with singularities in the orientation preference map [21]. The arrangement of CO blobs is reflected anatomically by the distribution of intrinsic horizontal connections, which tend to link blobsto-blobs and interblobs-to-interblobs [19,22,23], and by extrinsic corticocortical connections linking blobs to specific compartments in V2 and other extrastriate areas [19,24]. Taken together these observations suggest that the CO blobs are sites of functionally and anatomically distinct channels of visual processing. Note that CO blobs have also been found in cat, but tend to exhibit weaker functional properties [25,26]. For example, although blobs avoid OD borders in cat visual cortex, they are not strictly aligned along the center of the columns.
From a developmental perspective, the fact that cells in different layers of primate V1 may have different sets of feature preferences and receive distinct classes of thalamic input, suggests that one should consider a multi-layer cortical model. Such a model would first need to explain how direct thalamic inputs to different layers develop certain common features. For example, how do the OD maps of the M,P and K pathways align with each other? Second, it would need to account for the emergence of layer specific features as exemplified by the CO blobs. In this paper we construct such a laminar network model in order to investigate the joint development of OD columns and CO blobs in primate V1. For simplicity, we focus on layers 4C and 2/3. Both layers receive direct thalamic inputs and layer 4C sends vertical projections to layer 2/3. We assume that the OD map first develops in layer 4C and that this is then inherited by the thalamic inputs to layer 2/3 via the vertical projections from layer 4C. We also assume that both the thalamic and the vertical connections are modifiable by activity during the early development of OD columns. Using a correlation-based Hebbian learning rule with subtractive normalization [27], we show how the feedforward vertical projections from layer 4C to layer 2/3 compete with the direct thalamic input to layer 2/3, resulting in the formation of CO blobs superimposed upon the OD map. As in previous activity-based models [28,29,27], the spacing of the OD columns is determined by the spatial profile of the intralaminar connections within layer 4. This is consistent with recent experimental findings by Hensch and Stryker [30], who show that altering the lateral inhibitory circuits in kitten visual cortex during development results in a widening or narrowing of the ocular dominance columns. The spacing of CO blobs, on the other hand, depends both on the width of the OD columns inherited from layer 4 and the spatial distribution of intralaminar connections within superficial layer 2/3. We analyze our multi-layer developmental model by carrying out a perturbation expansion with respect to the strength of coupling of the intralaminar connections.
Such a perturbation expansion is implicit in standard single-layer correlation-based Hebbian models [27,3]), in which the steady-state activity within the layer is obtained by inverting a nonlocal operator that depends on the distribution of intralaminar connections. The resulting inverse operator consists of a convolution with a modified interaction kernel, which is usually taken to be a Mexican hat function.
Such an interaction kernel generates a spatially correlated competition between left and right eye inputs that ultimately leads to the formation of OD columns. One of the problems with taking a Mexican hat function, however, is that the convolution kernel is not invertible. On the other hand preserving invertibility means that the OD pattern becomes more sensitive to initial conditions. As we show in this paper, such issues become even more significant within the context of a multi-layer net-work model. Another feature that emerges from our analysis is that it is necessary to use an adaptive subtractive normalization scheme in order to preserve locally the total synaptic density throughout the developmental process. That is, the form of the subtractive normalization changes whenever one of the weight components vanishes.
It should be noted that a developmental model of the laminar structure of V1 has also been constructed by Grossberg and Williamson [31]. However, their model is concerned with the development of the feedforward, horizontal and feedback circuits within V1 that are thought to be involved in perceptual grouping, attention and learning in adult cortex. They do not address the issue of how cortical feature maps develop across different layers. The joint development of OD columns and CO blobs has previously been studied by Nakagama and Tanaka [32]. They consider an abstract thermodynamic model for the joint development of CO blobs and OD columns, in which layers 2/3 and 4C are collapsed into a single effective layer. Although their model successfully generates patterns that are consistent with experimental data, it is difficult to interpret the model directly in terms of the various synaptic pathways from LGN to V1 and within V1. Moreover, the single layer model assumes that the development of the CO blobs mutually interacts with the development of the OD columns, whereas in our more biologically realistic twolayer model the development of the OD columns influences the initial formation of the CO blobs but not vice-versa.

Laminar architecture of V1
We begin by briefly describing the layered structure of macaque V1. Further details and references can be found in the reviews by Levitt, Lund and Yoshioka [7] and Callaway [8]. There are three major classes of LGN neurons projecting to the cortex: magnocellular (M), parvocellular (P), and koniocellular (K). Each class of LGN neuron receives input from a different type of retinal ganglion cell and sends axons to distinct regions of V1.
(i) The M pathway originates from retinal ganglion cells that innervate the two most ventral M layers of the LGN. They have large receptive fields, are color insensitive, respond transiently to visual stimuli, prefer low spatial frequencies and are relatively sensitive to luminance contrasts. Thus neurons in the M pathway are excellent for detecting subtle changes in luminance or rapidly moving stimuli, but are poorly suited for the analysis of fine shape or color. The corresponding M-type LGN neurons project predominantly to layer 4Cα with a weaker projection to layer 6.
(ii) The P pathway originates from retinal ganglion cells that innervate the four most dorsal P layers of the LGN. They tend to be more numerous in the central retina and have smaller receptive fields than those in the M pathway, allowing them to convey information on a finer spatial scale. They are also able to detect color contrast. On the other hand, P cells have a slower, more sustained response to visual stimuli, making them less useful for the detection of rapid movement. The P-type LGN neurons project predominantly to layer 4Cβ with weaker projections to layers 4A and 6.
(iii) The K pathway originates from the small retinal ganglion cells that innervate the thin intercalcated regions between the M and P layers of macaque LGN. Until recently, these cells were thought to play a minor role in visual perception. However, there is a growing recognition of the importance of the K pathway, in large part due to extensive studies of new world galago monkeys, whose LGN has more distinct K layers than old world primates such as the macaque, thus making them more amenable to experimental analysis [17]. Across various primate species one finds that K-type LGN cells have distinct neurochemical and receptive field properties, and specifically target the cytochrome oxidase (CO) blobs in layer 2/3 of cortex [18]. There are also some K cells that target layer 1. Although the functional role of the K pathway is still unclear, one likely contribution in diurnal primates is color processing at short wavelengths.  Figure 1 Schematic diagram illustrating the main input pathways from LGN to V1 of macaque, and the major types of projections from layer 4C to superficial layer 2/3. See text for details. Adapted from [8].
Spiny stellate neurons in layer 4C of macaque V1, which are the recipients of thalamocortical inputs, can be subdivided into at least four groups according to the location of their dendritic fields within the layer and the targets of their projecting neurons in superficial layers [22,23]. This is illustrated in figure 1. The first group consists of neurons located in upper 4Cα whose dendrites are confined to 4Cα and whose axons arborize in layers 4B and in the CO blobs of layer 2/3. The second group has narrowly stratified dendrites that are confined to lower 4Cα and axons that target the interblobs in layer 2/3. Both of these neuron types are expected to belong to the M pathway. The third group straddles the border between 4Cα and 4Cβ and is thus likely to receive both M and P afferents from the LGN, whereas the fourth group is found in layer 4Cβ and has dendrites confined to layer 4Cβ. The axons of cells belonging to the last two groups arborize in layers 4A and layer 2/3 blobs and interblobs. Note that there are some differences across primate species.
For example, the interblob regions of diurnal primates such as the squirrel monkey receive inputs only from P-type neurons in layer 4Cβ, whereas the interblob regions of nocturnal primates such as the galagos receive inputs only from M-type neurons in layer 4Cα [33]. In cats, W and Y cells of the LGN project directly to CO blobs whereas X cells project indirectly via layer 4 to both blobs and interblobs [34]. Another important aspect of the laminar architecture of primary visual cortex is that each cortical layer sends its primary interlaminar output to one other layer. This leads to a simplified two-level model of information processing by local circuits of V1 as illustrated in figure 2 [8]. Each level consists of one feedforward module and one feedback module. The feedforward modules (layer 4C at level 1 and layers 2-4B at level 2) receive strong excitatory inputs from the previous level and make strong connections to the next higher level. The corresponding feedback modules (layer 6 at level 1 and layer 5 at level 2) receive weaker input from the feedforward modules at the previous level and at the same level, and provide heavy feedback to the latter. It follows that the most direct path for information flow from the LGN to extrastriate cortex is from LGN to layer 4C to layers 2-4B to layer 4 of extrastriate cortex.
Note that the feedforward connections tend to be focused whereas the feedback connections tend to be more diffuse in nature and play a modulatory role. A similar type of circuit exists in cat V1 if we consider layer 4C to be analogous to cat's layer 4 and layers 2-4B analogous to layer 2/3.

Reduced two-layer cortical model
The laminar architecture outlined above is already present in newly born primates, as are ocular dominance columns and cytochrome oxidase blobs. The prenatal development of OD columns is usually modeled in terms of a Hebbian-like competitive mechanism for the modification of left/right eye afferents from the LGN to layer 4C [28,29,27,3,35].. In the case of primates, it is assumed that prenatal neural activity is generated spontaneously either by retinal waves or by activity arising from the cortico-geniculate feedback loop [36]. Most developmental models do not distinguish between sublayers 4Cα and 4Cβ nor the corresponding M and P pathways highlighted in figure 1. However, the observation that a subclass of M-type neurons in upper layer 4Cα innervate the interblobs, while LGN neurons from the K layers and a second subclass of M-type neurons in lower 4Cα innervate the blobs, suggests that there is some Hebbian-like competition between these two classes of afferents within layer 2/3. This motivates us to consider a two-layer cortical model for the joint development of OD columns and CO blobs, which is consistent with the two-level architecture shown in figure 2. The basic assumptions of our model are as follows: 1.) Ocular dominance columns in layer 4C are generated by an activity-dependent synaptic modification of feedforward afferents from the M and P layers of the LGN, involving Hebbian-like competition between left and right eye afferents.
The resulting pattern of OD columns is passed on to superficial layer 2/3 by vertical feedforward connections from layer 4C to layer 2/3.
2.) Cytochrome oxidase blobs in layer 2/3 are generated by an activity-dependent synaptic modification of feedforward afferents from the K layers of the LGN and M-type feedforward afferents from layer 4C. This leads to an alternating pattern of blobs that are innervated by K-type inputs and one class of Mtype inputs, and interblobs that are innervated by a second class of M-type inputs. Since the drive from layer 4C is approximately monocular, the resulting pattern of CO blobs aligns with the OD columns through the underlying Hebbian mechanism. The absence of feedback from layer 2/3 to layer 4C implies that the development of the CO blobs does not affect the formation of the OD columns.
Our model assumes that the vertical connections between layer 4C and layer 2/3 are modifiable by neural activity during the period when OD columns and CO blobs form. There is some evidence for this in studies of cat visual cortex, where projections from layer 4C to layer 2/3 continue to mature until around postnatal week 6 whilst OD columns and CO blobs begin to appear around postnatal week 3 [37]. Interestingly, vertical connections between other layers tend to mature much earlier. Indeed, various in vivo and in vitro studies have shown that layerspecific intracortical connections in mammalian visual cortex tend to develop with a high degree of specificity from the outset and are not the result of subsequent remodeling (see the reviews by Katz and Callaway [38] and Bolz et al. [39]). This suggests that growing axons are able to distinguish specific cortical layers during development using activity-independent cues such as molecular markers. However, this picture is likely to be oversimplified, since a number of recent studies indicate that activity-dependent mechanisms may regulate layer-specific axonal targeting for certain types of pyramidal neuron targeting certain layers [40,41,42]. Plasticity in vertical connections has also been found in experimental studies of development in the rodent somatosensory cortex [43]. Early in development (first postnatal week) whisker manipulations induce rapid changes in the whisker barrel map in layer 4C, consistent with plasticity of thalamic afferents. On the other hand, in older animals plasticity tends to occur first in layer 2/3 (and layer 5) and only later (if at all) in layer 4. A major cellular component of plasticity in layer 2/3 is deprivationinduced weakening of the layer 4 to layer 2/3 vertical projections via long-term depression [44]. Plasticity in layer 2/3 has also been found during the critical period of cat V1, and appears to precede modifications in the direct thalamic input to layer 4 during monocular deprivation experiments [45]). To what extent this involves changes in vertical rather than intralaminar horizontal connections has not yet been determined.
Another assumption of our model is that the CO blobs and OD columns emerge at roughly the same time, with the OD columns occurring slightly earlier. This is consistent with experimental findings, even if the precise temporal ordering is not yet clear. In macaque both CO blobs and OD columns emerge prenatally, so that at birth the pattern of OD columns and their spatial relationship with blobs is adultlike. In contrast, the cat's visual cortex is quite immature at birth. For example, supragranular layers of cat V1 differentiate postnatally and the blobs in these layers are normally first visible around 2 weeks of age (about one week after eye opening and one week before the critical period). This is approximately coincident with the earliest observation of OD columns in cat [46].

Mathematical formulation of model
In order to develop the mathematical formulation of our model we make two additional simplifications. First, we ignore the spread of synaptic inputs between the LGN and layer 4C, and assume that there is a one-to-one retinotopic map between each LGN layer and each cortical layer. (This assumption will be relaxed later on, see section 3.2). Thus, after an appropriate rescaling, we can represent each layer by Σ ⊂ R 2 and uniquely label each point in an LGN layer by the corresponding point in cortex that it is connected to. Second, we consider only a single M-type feedforward pathway from layer 4C to layer 2/3, which we assume innervates the interblobs in the mature cortex. Inclusion of a second M pathway together with a pair of P pathways as described in section 2.1 complicates the model without altering the basic mechanism identified in our study. A schematic diagram of the simplified model is shown in figure 3. Let w L (r) and w R (r) denote the densities of left and right M-type LGN connections to a point r = (x, y) of layer 4C in cortex at time t, and denote the corresponding K-type densities from LGN to layer 2/3 by k L (r) and k R (r). Suppose that in each layer there are weak cortical-cortical interactions between neurons at r and r given by the synaptic density J i (r − r ) where layers 4C and 2/3 are indexed by i = 1 and i = 2 respectively. There are also vertical feedforward connections from M -type neurons in layer 4C to layer 2/3 with synaptic density m(r). Assuming a linear model for the cortical activity V i (r, t) in the two layers, leads to the pair of coupled linear equations where I L,R (r) and I L,R (r) denote the left and right inputs to the M and K layers of the LGN respectively. Note that although magnocellular and koniocellular geniculate cells have different receptive field properties, we expect I L (r), I L (r) to be strongly correlated, and similarly for the pair I R (r), I R (r). The constant coupling constant ε determines the relative strength of the lateral connections.
LGN Since development takes place on a much slower time-scale than the dynamics of cortical activity, we can take V i to be given by its steady-state value. However, calculating the steady-state explicitly requires inverting the linear operator In the case of weak intracortical interactions, this inversion can be carried out by performing a perturbation expansion in ε. The first-order approximation is thus and and δ is the Dirac delta function. Note that we explicitly include the δ-function component of the inverted operator L −1 i . Usually this term is ignored and G i is treated as a smooth weight function [27,3].
As we show in the next section, the δ-function can have a significant effect on the weight dynamics. Given the steady-state postsynaptic response to an input in the case of fixed weights, we now allow the feedforward weights to vary slowly in time according to a linear Hebb rule with subtractive normalization [47]: where . . . denotes averaging over the distribution of M-type inputs I L,R and Ktype inputs I L,R . The decay terms γ(w) and η(k, m) enforce various conservation constraints whose specific form will be derived in section 3. Equations (2.5)-(2.7) are supplemented by additional constraints that ensure the weights remain positive and bounded. That is, The matrix C determines the correlations between the two eyes. We will take the eyes to be symmetric and to drive development equally. That is, we take C LL = C RR = C S and C LR = C RL = C D so that C and is symmetric such that the diagonal terms are the same eye correlations C S and the off-diagonal terms are the different eye correlations C D . This leads to the vector equation where H 1 (r) = G 1 (r)Q(r), a = (1, 1) T and w = (w L , w R ) T . The subtractive term γ(ω)a is chosen so that the total synaptic density w L (r) + w R (r) is conserved at each point in cortex. Thus, we take Exploiting the fact that the input correlation matrix C has eigenvalues µ ± = C S ± C D with corresponding eigenvectors e ± = (1, ±1), it is straightforward to show that the vector equation (3.2) decomposes into the pair of decoupled equations with w ± = w L ± w R . Conservation of total synaptic density is thus achieved by setting µ = µ + /2 so that dw + (r, t)/dt = 0 for all r and t. For a more general choice of correlation matrix C, the subtractive normalization term takes the form The remaining equation for dw − /dt describes the growth of ocular dominance columns and can be analyzed using Fourier transforms. For the sake of illustration, suppose that Σ is the unbounded two-dimensional plane. Then A p e ip·r e λ(p)t/τw dp, where p = |p| and λ(p) The complex amplitude A p is determined by the initial conditions: We expect correlations between the same eye to be greater than correlations between opposite eyes, that is, which is a real symmetric function of the wavenumber p, has a unique positive maximum at p = p c . The standard way of implementing this condition is to take where H consists of short-range excitation and long-range inhibition, as realized by either a difference-of-Gaussians or a difference-of exponentials. The former takes the form with σ E and σ I determining the range of excitation and inhibition respectively. The ultimate growth of the modes is restricted by the constraints 0 ≤ w L,R ≤ W .
Numerically one finds alternating bands of left/right eye ocular dominance columns whose mean spacing is determined by the critical wavelength (see also section 4).
Layer 4C can thus be partitioned into left and right OD regions Σ L,R such that lim t→∞ w(r, t) = W(r) with W(r) = (W, 0) T for all r ∈ Σ L and W(r) = (0, W ) T for all r ∈ Σ R . A similar analysis holds when Σ is a bounded domain with the integral in equation (3.7) replaced by a discrete sum.

Inclusion of δ component
One subtle aspect of the above analysis that is often overlooked concerns the inversion of the linear operator L 1 , under the assumption that the lateral interactions within layer 4C are relatively weak compared to the feedforward afferents from the LGN. In particular, G 1 (r) = δ(r) + εJ 1 (r) so that H 1 (r) = Q(0)δ(r) + εQ(r)J 1 (r).
The contribution from the term Q(0)δ(r) is usually ignored and Q(r)J 1 (r) is identified with the Mexican hat function H(r). At first sight such a simplification does not appear to be significant, since the Fourier transform of the term Q(0)δ(r) simply introduces a p-independent shift in H(p) so that the critical wavenumber p c remains the same, see figure 5(a). However, such a shift destabilizes the uniform mode p = 0, which then amplifies the homogeneous component of the initial state.
The amplitude A 0 of this component will depend upon the initial balance between left and right eye inputs, since A 0 = Σ {w L (r, 0)−w R (r, 0)}dr. We define a measure, ξ, for the total fraction of left eye input innervating layer 4C as (3.11) The measure ξ varies from 0 to 1, where 0 corresponds to total right eye dominance, 1 is total left eye dominance, and 1/2 is a balanced L/R state. A comparison of the dynamics with and without the δ-function term is illustrated in figure 5 induces a shift to the right of the sensitivity plot shown in figure 5(b). In order to understand this result, let us return to the Fourier decomposition (3.7) of w − (r, t).
The growth of the OD pattern can be viewed as a race between all modes p for which λ(p) > 0. Classically, one considers the mode that corresponds to the critical wavenumber p c , i.e., the "fastest" growing mode, to dominate the final pattern that emerges once the saturating nonlinearities of equation (2.8) have been imposed.
With the inclusion of the δ-component, the zeroth mode has also entered the race so that an approximately uniform initial state means that |A 0 | |A p | for all p = 0 which, in effect, gives the zeroth mode a large head start in the race. Thus, the resulting pattern becomes dependent on A 0 and hence the initial balance between left and right eye inputs. Let T and T represent the times for the critical mode p c and the zeroth mode, respectively, to reach the upper bound W . In order that segregation of the L and R visual pathways occurs, the parameter regime must be such that T < T . That is This inequality is satisfied provided that where α c = λ(pc)−λ(0) . Note that if |A pc | |A 0 |, then the zeroth mode dominates and the cortex would develop into a wholly left or right eye driven state. Inclusion of the noise term ξ(r) in the initial state increases the ratio |A pc |/|A 0 | so that the inequality (3.14) is satisfied.

Development of OD columns and CO blobs in layer 2/3
As the ocular dominance columns form in layer 4C, they influence the development of the feedforward afferents m(r) connecting layer 4C to layer 2/3 as well as the koniocellular afferents k L,R (r) projecting from the LGN to layer 2/3. Motivated by experimental data in the cat [37], we will assume that the OD columns in layer 4C form before the OD columns and CO blobs in layer 2/3 so that we can take w(r, t) → W(r) in equations (2.6) and (2.7). (Note that even if they were to develop simultaneously, the analysis would still be valid to a lowest order approximation, since there is no direct feedback from layers 2/3 to layer 4C). Following along similar lines to the analysis of equation (2.5), we introduce the input correlation matrices C and B according to

It is clear from equations (3.17) and(3.18) that when considering correlation-based
Hebbian learning in a multi-layer network, the resulting dynamical equations for the weights rapidly become complicated due to the presence of multiple convolutions. This makes it difficult to gain insights into the basic mechanisms underlying feature map formation in the higher layers, and also makes numerical simulations more difficult. However, considerable simplification can be obtained by carrying out a perturbation expansion to first order in the coupling parameter ε. Such an approximation is already implicit in the inversion of the linear operators L i . In order to proceed we will need to include the Dirac delta function in L −1 i . Indeed, considerable insight into the underlying linear algebraic structure of the Hebbian dynamics can be obtained by setting ε = 0 in equations (3.17) and (3.18) so that G i (r − r ) → δ(r − r ) . This corresponds to ignoring the effects of lateral interactions; such interactions will subsequently be incorporated by carrying out a first-order perturbation expansion in ε.
The binary-valued functions χ L,R (r) arise from the simplification that there is a one-to-one vertical interlaminar projection from layer 4C to layer 2/3. We now relax this assumption by allowing the vertical connections to sample over a localized region of layer 4C. The sampling region reflects the spatial extent of the dendritic field of the forward projecting neurons, which we assume is a fixed arbor function.
We can incorporate the effect of a dendritic field into the above model by taking the functions χ L,R (r) to represent the degree of left and right eye drive from layer 4C projected to point r in layer 2/3 such that 0 ≤ χ L,R (r) ≤ 1. More specifically, introducing the normalized Gaussian J f (r) with space constant σ f and setting W = 1, we consider Under this formulation, the synaptic input to a cortical location r due to the left and right eye streams is sampled from the footprint and is subject to a transfer function Φ, where we take Φ(ν) = √ ν. With this choice of Φ, χ 2 L (r) + χ 2 R (r) = 1 for all r. Naturally the smooth variation in χ L,R sampled from the above layer naturally results in the emergence of binocular regions along the OD borders in layer 2/3 (discussed below). As there is likely a similar integration of left and right eye drive occurring in layer 4C due to input from the LGN and accounting for cortical-thalamic feedback (reviewed in [48]), we use this formulation χ L,R (r) in Near the center of a left-eye OD column Υ L (r) ≈ 1 and Υ R ≈ 0, near the center of a right-eye OD column Υ L (r) ≈ 0 and Υ R ≈ 1, and near the OD border Υ L (r) ≈ In the case of primates, the development of OD columns and CO blobs occurs prenatally so that the Hebbian process is driven by spontaneous activity in the left and right eyes. Therefore, we expect opposite eye correlations to be small and we can set C D = C D = B D = 0 for simplicity. As a further simplification we also set C S = C S > B S , although our results generalize to the case C S = C S . We then rewrite equations (3.21)-(3.23) in the matrix form where we have set τ k = τ m = 1 for convenience and dropped the explicit dependence on r. Following along analogous lines to the analysis of OD column formation in layer 4C, we choose the subtractive normalization term η 0 in order to conserve the total synaptic density k L + k R + m at each point in layer 2/3. However, there is an additional subtlety that arises when considering three rather than two weight components: the normalization term must be redefined in the event of one compo-nent vanishing. As we show below, this occurs in the center of OD columns that are inherited from layer 4C, where either k R or k L vanishes and the remaining component subsequently competes with the m input. In order to conserve the total synaptic density throughout the developmental process, we define the subtractive normalization in a dynamic way: It can be seen that the k − dynamics decouples, whereas the components k + and k R compete against one another. We find that for a wide range of parameter values the right-eye population (k R ) vanishes leaving the left-eye population (k + ) intact.
The adaptive normalization scheme then leads to the new dynamics (3.30) in which the remaining left-eye K and M inputs now compete with one another.
We simulate the dynamics for a representative choice of parameters and display the results in figure 6(a). Case II: OD border (Υ L (r) = 1/2 : χ L = χ R = χ) Given the constraint that k L + k R + m is locally conserved, we can reduce the dynamics at the given point r to the following pair of equations where k tot = k L + k R . The associated matrix is degenerate and has the following non-zero eigenvalue Provided 3C S − 4χB S > 0, the K and M pathways will compete as the associated For the given initial conditions, we find that k L , k R → 0 at the OD border leaving a purely interlaminar, binocular magnocellular drive to the cortical point r, see figure 6(c). Similar behavior occurs in a region surrounding the OD border where Υ L (r) ≈ 0.5, see figure 6(b).
The above examples show that, given an approximately balanced uniform initial state, layer 2/3 rapidly inherits the OD columns of layer 4, with k R (r) = 0 when r ∈ Σ L and k L (r) = 0 when r ∈ Σ R . Moreover, with the inclusion of the dendritic sampling from layer 4C, the koniocellular weights k L,R both decay to zero at the OD borders where Υ L , Υ R ≈ 0.5 (see figure 6(c)), ultimately sequestering the CO blobs away from the OD borders. When we include spatial interactions, the competition within an OD column between the K-type and M -type inputs will result in the segregation of the K and M pathways, i.e., the formation of the blob/interblob areas.

O(ε) analysis
In order to incorporate the leading order effects of the lateral interactions, we now determine the O(ε) contributions to equations (3.17) and (3.18): and We have introduced the convolution operator * defined according to for any function f (r) where K i (r) = J i (r)Q(r) for i = 1, 2. As in the analysis of the zeroth order dynamics, we assume for simplicity that B D = C D = 0 and C S = C S . Substituting the approximations V 2 I = V 2 I 0 + ε V 2 I 1 and V 2 V 1 = V 2 V 1 0 +ε V 2 V 1 1 into equations (2.6) and (2.7), and using equations (3.19), (3.20), (3.34), (3.34) then gives Here we use the notation m L (r, t) = χ L (r)m(r, t) and m R (r, t) = χ R (r)m(r, t) (3.38) and set ψ L (r) = K 1 * Υ L (r) and ψ R (r) = K 1 * Υ R (r). (3.39) As in the O(1) case, we have to consider an adaptive normalization scheme in order locally to conserve total synaptic density throughout the developmental process.
Thus η is of the form where η j i denotes the subtractive normalization term for the ith order contribution at the jth stage. The terms η j 1 are chosen so as to conserve total synaptic density at O(ε). The O(1) analysis showed how the layer 2/3 system can rapidly inherit the ocularity preference of layer 4C. We will assume that the OD segregation is maintained at O(ε) on taking the subtractive normalization η as defined above. The preservation of the OD map in layer 2/3 will be confirmed numerically in section 4. The occurrence of such a segregation, which is confirmed numerically in section 4, generates the CO blobs. To further illustrate this, consider a point r at the center of a left-eye dominated region where χ L ≈ 1 and χ R ≈ 0. Equations (3.35)- (3.37) reduce to the simpler form where k L ± = k L ± m. A similar equation holds for r at the center of a right-eye dominated region under the change L → R. The fact that the convolution with respect to the interaction kernel H 2 in equation (3.42) is restricted to same-eye ocular dominance strips (r, r ∈ Σ L ) rather than the whole two-dimensional layer, suggests that the pattern formation process is quasi-one-dimensional. That is the center of each OD column can be idealized as a narrow stripe that is partitioned into alternating regions of K and M afferents, that is, blobs and interblobs. To a first approximation, the spacing of the pattern within a single stripe will be determined by the Fourier transform of H 2 along the stripe. Weak interactions between sameeye stripes will then determine the alignment of patterns across the cortex. This basic picture is consistent with the numerical simulations presented in section 4.
As in layer 4C, we take the effective interaction function to be H 2 = δ + εH with H given by the Mexican hat function (3.10). Hence, as in the layer 4C model, the presence of the Dirac delta function means that the resulting CO blob pattern will be sensitive to the initial balance of M and K afferents. This is handled in our numerical simulations (see below) by adding to the uniform initial state a modified white noise term analogous to equation (3.12).

Numerical results
In this section, we present numerical results obtained by solving the reduced set of We considered a wide array of weight functions and found that an approximate balance of excitation and inhibition performs most robustly. For example, when excitation strongly dominates the system's sensitivity to homogeneous perturbations is exacerbated, and in the regime where inhibition dominates the growth terms for high wavenumbers do not significantly differ from the growth term corresponding to the critical wavenumber, hence there is a mixing of modes resulting in a highly disorder pattern. With lateral circuitry that has an approximate balance of excitation and inhibition, our model captures several qualitative aspects of CO blob and OD column formation in superficial layers of primate V1. In particular, it reproduces the experimental observation that the CO blobs are localized around the centers of the columns. In our numerical simulations, we neglected the inhomogeneous terms involving ψ L , ψ R . This approximation would be valid, for example, under the assumption that the intracortical interactions in the upper laminae are stronger than in the input layer 4C, that is, |K 1 | |K 2 |. It can be checked numerically that the model performs similarly under the inclusion of the ψ L , ψ R terms, although it does become more sensitive to the choice of initial conditions.
As a further comparison with experimental data, we calculate some quantitative aspects of a typical CO blob distribution generated by our model as shown in figure   9(a). Following Murphy et al. [49], we first determine the mean CO blob spacing of the pattern using a Delauney triangulation. We find that the mean edge width of the Delaunay triangulation is approximately twice the mean width of an OD column. Taking the latter to be approximately 385 µm [50], yields a mean blob spacing of approximately 716 µm. This is a slight overestimate of the mean value of 590 µm found experimentally by Murphy et al. [49], which could be an artifact of the heightened disorder in the simulated system compared to experimentally found patterns. However, the Gaussian-like shape of the histogram of edge lengths (figure 9(b)) is consistent with experimental data. Another characteristic feature of CO blob distributions in macaque is that the spacing between blobs across OD columns is smaller than between blobs along the same OD column. The former is typically 353 µm [51,52] whereas the latter is approximately 550 µm [16]. An interesting feature of our model is that the profile of the local connectivity in layer 4C determines the periodicity of the OD pattern, whereas the local connectivity in  appear roughly as skewed Gaussians, which is consistent with experimental findings.
The ratio of the CO blob spacing along the OD band and across the border is found to be approximately 1.56 in monkeys [51,16], which is in very good agreement with our simulated distribution that has a ratio of 1.50. Furthermore, Horton [16] noted that there is also an asymmetry in the shape of the blobs themselves. The blobs are oval shaped of dimension 250 x 150 µm, which suggests that the O(1) dynamics determine the length of the minor axis, whereas the spatial interactions within an OD stripe determine the length of the major axis of the CO blob shape. We do not track statistics for the blob shapes, but note that they have an oval-like profile. The differences in our model results and experimental observations could be rectified by an alternate choice of the layer 2/3 circuitry and perhaps an alternate initial balance of the K and M pathways. The experimental data does not suggest a precise choice for these parameters and hence we are only suggesting the underlying structure of the model and representative parameter choices.

Monocular Deprivation
In order to consider the effects of monocular deprivation during development, we relax the assumption that our correlation matrices are symmetric. In particular, we consider the correlation matrices for the magnocellular(M), koniocellular(K) and correlations between the K and M pathways are respectively given by The evolution of the synaptic weights is then given by where η 0 is the corresponding subtractive normalization term and χ L,R are defined in equation (3.24).
We consider normal development until MD onset, then reduce the correlation coefficients associated with the deprived eye, similar to the approach in [27] to simulate MD. At the initiation of MD, there is an approximate balance of left/right eye dominated regions with CO blobs centered within the OD stripes. Below in figure   10, we display the simulated effects of MD on the OD pattern and distribution of the CO blobs. Without any loss of generality, we take the right eye to be deprived (black regions) while the left eye corresponds to the non-deprived eye (white regions), and CO blobs are colored in gray. At the onset of MD, the correlation coefficients corresponding the right eye activity driving cortical activity are reduced (as in [27].
MD results in the broadening of the OD stripes associated with the non-deprived eye broaden and a reduction of the OD stripes associated with the deprived eye, while there are minimal alterations of the CO blob positions. This numerical study is inline with the experimental finding of Murphy et al. [53].

Discussion
In this paper, we constructed a multi-layer, activity-dependent model for the joint development of ocular dominance columns and cytochrome oxidase blobs in primary visual cortex of primates. We assumed that the OD map first forms in layer 4C and is then inherited by the supragranular layers via feedforward vertical projections from layer 4C. Competition between these feedforward projections and the direct thalamic input to layer 2/3 then resulted in the formation of CO blobs superimposed a) b) before MD after MD Figure 10 Effects of monocular deprivation on the ocular dominance pattern and CO blob distribution. Here black corresponds to the deprived eye, whereas white indicates the remaining open eye, while CO blob regions are colored in gray. In (a), the OD pattern and CO blob distribution immediately before the onset of MD is shown, whereas in (b) the aftereffects of MD is pictured. Before deprivation the non-deprived eye area constituted 48.2% of region, while post deprivation the non-deprived controlled 69.3% of the cortical region. Note that while the OD pattern is significantly modified, the CO blob distribution remains preserved. The domain size is 6 × 12 and is meant to simulate a 2.3 mm by 4.6 mm region of cortex. The parameters used in the simulation prior to MD are identical as those in figure 7 with the exception that = 0.06 and C d = 0. upon the OD map. The CO blob distribution obtained from numerical simulations was consistent with experimental data. First, the CO blobs were located in the centers of OD columns. In our model this was related to the fact that the feedforward projections sampled over a local region of layer 4C resulting in binocular regions along the OD borders. Second the CO blob lattice exhibited a characteristic asymmetry, in which the spacing of a CO blob from the nearest neighbor blob driven by the same-eye was on average greater than the spacing from the nearest neighbor blob driven by the alternate-eye. In our model this was a consequence of the fact that the periodicity of the OD map that develops in layer 4C was determined via the lateral circuitry in layer 4C, whereas the periodicity of the CO blob lattice within an OD stripe was determined via the circuitry in layer 2/3. For simplicity, we used a linear correlation-based Hebbian rule with subtractive normalization. The advantage of a linear model is that one can explicitly determine the propagation of weight modifications through the layers. However, one of the possible limitations of such a model is the sensitivity to initial conditions. It might be possible to construct a more robust model by considering a nonlinear competitive Hebbian learning rule [54], although the analysis would be considerably more involved.
As a first approximation we considered a purely feedforward model in this paper.
This was partly motivated by the observation that there are no direct feedback connections from layer 2/3 to layer 4C. However, connections between layer 5 and layer 6 could provide a more indirect feedback pathway for activity in layer 2/3 to modulate the activity in layer 4C. Furthermore, it has been demonstrated in cat that at the onset of monocular deprivation (MD), plasticity first occurs at the OD borders in layer 2/3 and only after a delay of several days does any modification in the direct thalamic input to layer 4 occur, suggesting that layer 2/3 influences the development of layer 4 during the critical period [45]. Studies in mice indicate that layer 2/3 plasticity induced by monocular deprivation involves a rapid restructuring of dendritic spines that is regulated by GABA inhibition (reviewed in [55]). Hence, just as in early development when OD columns and CO blobs first form, it is likely that the laminar structure of cortex plays an important role during the critical period. In particular, one cannot model monocular deprivation simply as a change in input correlations to a single layer network. We hope in future work to extend our laminar network model to analyze development during the critical period, and to understand the delayed occurrence of plasticity in layer 4 during monocular deprivation by considering the role of feedback from layer 2/3 to layer 4. Finally, note that the principles outlined in this paper could be extended to include orientation selectivity and spatial frequency in the visual system, and could also be employed to study the development of laminar structures in other cortical areas.