Skip to main content
Advertisement
  • Loading metrics

Laminar Neural Field Model of Laterally Propagating Waves of Orientation Selectivity

  • Paul C. Bressloff ,

    Contributed equally to this work with: Paul C. Bressloff, Samuel R. Carroll

    bressloff@math.utah.edu

    Affiliation Department of Mathematics, University of Utah, Salt Lake City, Utah, United States of America

  • Samuel R. Carroll

    Contributed equally to this work with: Paul C. Bressloff, Samuel R. Carroll

    Affiliation Department of Mathematics, University of Utah, Salt Lake City, Utah, United States of America

Abstract

We construct a laminar neural-field model of primary visual cortex (V1) consisting of a superficial layer of neurons that encode the spatial location and orientation of a local visual stimulus coupled to a deep layer of neurons that only encode spatial location. The spatially-structured connections in the deep layer support the propagation of a traveling front, which then drives propagating orientation-dependent activity in the superficial layer. Using a combination of mathematical analysis and numerical simulations, we establish that the existence of a coherent orientation-selective wave relies on the presence of weak, long-range connections in the superficial layer that couple cells of similar orientation preference. Moreover, the wave persists in the presence of feedback from the superficial layer to the deep layer. Our results are consistent with recent experimental studies that indicate that deep and superficial layers work in tandem to determine the patterns of cortical activity observed in vivo.

Author Summary

A major challenge in neurobiology is understanding the mechanisms underlying the formation and propagation of cortical waves and the underlying neural circuitry that supports them. A variety of neurological disorders such as epilepsy and spreading depression during migraine episodes are characterized by spatially localized oscillations and waves propagating across the surface of the brain. Moreover, traveling waves can be induced in vitro by electrically stimulating disinhibited cortical slices, and they are also observed in vivo. Physiological traveling waves are also correlated with perceptual waves such as those observed during binocular rivalry. Spatially structured neural network models (neural fields) are becoming increasingly popular as a method for exploring, both analytically and numerically, the conditions under which cortical traveling waves can occur. However, such models tend to ignore the laminar structure of cortex and the fact that sensory cortical cells are often tuned to particular features of a stimulus. In this study we develop a a laminar neural field model of propagating waves in primary visual cortex (V1), in which the orientation selectivity of V1 cells is taken into account. We demonstrate the important role of vertical inter-laminar connections and horizontal intra-laminar connections in supporting the propagation of orientation selective waves.

Introduction

A major challenge in neurobiology is understanding the mechanisms underlying the formation and propagation of cortical waves and the underlying neural circuitry that supports them [1, 2]. A variety of neurological disorders such as epilepsy [35] and spreading depression [68] are characterized by spatially localized oscillations and waves propagating across the surface of the brain. Moreover, traveling waves can be induced in vitro by electrically stimulating disinhibited cortical slices [912] and are observed in vivo. Physiological traveling waves are also correlated with perceptual waves such as those observed during binocular rivalry [1316]. Stimulus-induced and spontaneous laterally propagating waves in vivo have been observed in the primary visual cortex (V1) of anesthetized rodents [1719], ferrets [20], cats [2123], and primates [23, 24]. These observations have been obtained using various experimental methods, including optical imaging with voltage-sensitive dyes [17, 18, 22, 24], measurements of local field potentials (LFPs) [21, 23, 25], and calcium imaging [19, 20]. Two particular features of propagating activity in V1 motivate the modeling study of this paper. The first concerns the fact that most V1 cells respond preferentially to local stimuli with certain preferred properties such as orientation and left/right eye preference (ocular dominance). This means that propagation in cortical space is correlated with both retinotopy and stimulus feature preferences. Indeed, one can observe the lateral spread of orientation selectivity in carnivore V1 based on voltage-sensitive dye imaging [22], LFPs [21, 23] and epifluorescent imaging of calcium waves [20]. There is also indirect evidence for the propagation of orientation-dependent activity in V1 from experimental studies of binocular rivalry waves [13, 14]. The second feature concerns the laminar structure of V1, in particular, growing evidence that propagating activity in cortex is initially generated by local recurrent connections in deep (infragranular) layers, which then spreads vertically to superficial (supragranular) layers. This has been observed both in mouse V1 [19] and other cortical areas [2628].

In this paper we develop a continuum neural field model of propagating waves in V1 that takes into account both the orientation-dependence of V1 neurons and the laminar structure of cortex. We focus on animals that have structured orientation preference maps such as ferrets, cats and primates rather than the “salt-and-pepper” organization found in rodents. That is, we take superficial layers of cortex to have a hypercolumnar structure consisting of orientation columns organized around a set of pinwheels [2931], with strong local recurrent connections and weaker (modulatory) long-range horizontal connections that link neurons in different hypercolumns with similar orientation preferences [3235]. Following several modeling studies [3638], we assume that within each hypercolumn, neurons with sufficiently similar orientations tend to excite each other whereas those with sufficiently different orientations inhibit each other, and this serves to sharpen a particular neuron’s orientation preference. (Note, however, that the precise role of local recurrent connections in orientation tuning is still controversial, given the lack of direct evidence for antagonistic inhibition within cat V1 [39, 40].) Such a tuning mechanism suggests that local connections are structured with respect to orientation preference rather than retinotopy, and thus cannot provide a substrate for laterally propagating waves—this is also consistent with the observation of standing waves of orientation-dependent activity observed by Benucci et al. [21]. The weakness of the horizontal connections means that they also cannot support wave propagation on their own. In conclusion, the functional anatomy of superficial layers is consistent with experimental studies indicating that wave propagation is initiated in deep layers. Moreover, there is growing evidence that neurons in deep layers are more poorly tuned for orientation than those in superficial layers. For example layer 5 neurons in mouse V1 exhibit very little selectivity [41] and are weakly tuned in tree shrew [42]. Although orientation selectivity is observed in all layers within macaque V1, it appears to be weaker in deep layers [43, 44]. In the case of macaque, direct thalamic inputs from the lateral geniculate nucleus (LGN) are not orientation tuned, so that one possible source of weak orientation in deep layers is vertical feedback connections from input layer 4 or superficial layers. Note, however, that any orientation tuning in cat deep layers persists when feedback from superficial layers is suppressed [45]. If we assume that deep layer neurons are not strongly tuned for orientation, then it is possible that strong recurrent connections within the deep layers are spatially structured with respect to retinotopy rather than orientation preference, and thus support wave propagation.

The above experimental observations and hypotheses suggest that to a first approximation, propagating activity is initiated in orientation-independent deep layers, which subsequently induces the propagation of orientation-dependent waves in superficial layers via a combination of vertical inter-laminar connections and weak long-range horizontal connections within superficial layers. In this paper, we establish the validity of this emerging picture using a combination of mathematical analysis and numerical simulations. For simplicity, we consider a bilayer neural field model, which could loosely be interpreted in terms of an orientation-dependent superficial layer coupled to an orientation-independent deep layer. In our mathematical analysis, we construct traveling wave solutions in a simpler one-dimensional (1D) version of our model, and establish that inclusion of feedback from superficial to deep layer has little effect on wave propagation. We then show numerically that orientation selective traveling waves also occur in a two-dimensional (2D) model.

Finally, it should be noted that the laminar structure of V1 has been modeled extensively by Grossberg and collaborators [4648]. These models are highly detailed and have been used to investigate a wide variety of phenomena in early visual processing. In contrast, our much simpler model is used to investigate large-scale cortical dynamics rather than visual information processing. There has been relatively little work on analyzing the dynamics of layered neural fields, one exception being recent work on the effects of inter-laminar coupling on regularizing wave propagation in stochastic neural fields [4951], see also [52].

Laminar Neural Field Model

One mathematical approach to studying cortical waves involves the analysis of continuum neural field models, in which the large—scale dynamics of spatially structured networks of neurons is described in terms of nonlinear integro-differential equations, whose associated integral kernels represent the spatial distribution of neuronal synaptic connections [5358]. In this section we construct a neural field model of V1 that takes into account (i) the orientation selectivity of V1 neurons and (ii) the laminar structure of cortex.

In order to motivate our model, we first briefly recall the original ice cube model of V1 due to Hubel and Wiesel [59], see Fig 1a. The latter treats superficial layers of V1 as a collection of horizontally arranged orientation hypercolumns. Each hypercolumn consists of cells with similar spatial (retinotopic) coordinates x, which can be further partitioned into columns consisting of cells with similar orientation preference θ (highlighted by different colors in Fig 1). This structure can be modeled by collapsing each hypercolumn into a single point (through some form of spatial coarse-graining) and treating V1 as a continuum of hypercolumns [38, 60, 61]. Thus neurons are lndexed by the pair {x, θ} with x ∈ ℝ2 labeling the hypercolumn at (coarse-grained) position x and θ, −π/2 < θπ/2, labeling the orientation preferences of neurons within the hypercolumn This so-called R2 × S1 structure of the superficial layers can be made more explicit by representing orientation hypercolumns as vertically drawn cylinders, which collectively form a fibration of the retinotopic plane (Fig 1b), see [62, 63]. (These abstract vertical columns should not be confused with the actual vertical columns passing through the various cortical layers.) If the synaptic weights within each orientation hyperolumn are taken to be π-periodic (topology of a circle), then one can treat each hypercolumn as a ring network that can amplify weakly tuned inputs from the thalamus or deeper cortical layers to form population based tuning curves [36, 37]. On the other hand, patchy horizontal connections between hypercolumns can be represented as links between similar orientation columns within different hypercolumns.

thumbnail
Fig 1. Schematic illustration of R2 × S1 model of V1.

(a) Ice cube model of Hubel and Wiesel [59] showing cortical layers, orientation hypercolumns, and left(L)/right(R) ocular dominance columns. (b) Abstraction of superficial layers 2/3 as a collection of vertical cylinders (hypercolumns) that form a fibration of the plane. Patchy horizontal connections whoz are represented as links between the vertical cylinders. There are also local recurrent connections wloc within each hypercolumn.

https://doi.org/10.1371/journal.pcbi.1004545.g001

Another feature of the original ice cube model is that the vertical distribution of stimulus features across different cortical layers is taken to be uniform so that to a first approximation the layered structure of cortex can be ignored. However, it has subsequently been shown that in non-human primates, for example, different cortical layers exhibit different properties [6466]. Roughy speaking, the layers can be grouped into input layer 4, which is the main input layer from the thalamus, superficial layers 2, 3 that project to higher cortical regions, and deep layers 5, 6 that project back to thalamus and the brain stem. There is also extensive communication between the layers via so-called vertical connections. This is encapsulated by the modified canonical microcircuit shown in Fig 2. It has traditionally been assumed that the major feedforward pathway from thalamus is LGNL4 → L2/3 → L5/6 with L5/6 sending feedback projections to L2/3. However recent experimental work on rat barrel cortex has established that there is a strong direct input from thalamus to L5/6 [67]. This suggests that there could be a second feedforward pathway LGNL5 → L2/3 with respect to which L2/3 → L5 is viewed as feedback. Although the situation in the visual cortex of primates is less clear, we will explore the potential consequences of this second feedforward pathway for the propagation of orientation-selective waves.

thumbnail
Fig 2. Schematic diagram illustrating an updated version of the classical microcircuit introduced by Douglas and Martin [68].

This minimal circuitry comprises superficial (layers 2 and 3) and deep (layers 5 and 6) pyramidal cells and a population of smooth cells in layer 4C. Feedforward thalamic inputs target all cell populations. Weaker (modulatory) connections are shown in gray. In the original microcircuit model it was assumed that direct thalamic inputs to deep layers were weak. However, a recent experimental study indicates that these inputs are comparable to those innervating layer 4C [67]. Thus the classical direct pathway ThL4 → L2/3 → L5/6 (red) is supplemented by an additional direct pathway ThL5/6 → L2/3 (blue).

https://doi.org/10.1371/journal.pcbi.1004545.g002

In this paper, we consider a bilaminar neural field model representing an orientation-dependent superficial layer coupled to an orientation-independent deep layer. Note that neurons in deep layers of V1 tend to be more weakly tuned for orientation than neurons in superficial layers [41, 42, 44]. This tuning could be even weaker at the population level if the distribution of orientation preferences is more disordered in deep layers. As a first approximation, we treat the deep layer as non-orientation selective. However, the main point is not that deep layer neurons are independent of orientation, but that recurrent connections within the deep layer are untuned to orientation. A schematic diagram of the basic network architecture is shown in Fig 3, which combines the ℝ2 × S2 product structure of Fig 1 with the canonical microcircuit shown in Fig 2. The corresponding neural field equations take the form (1a) (1b) Here u(x, t) denotes orientation-independent activity in the deep layer and v(x, θ, t) is the orientation-dependent activity in the superficial layer. The recurrent weight distributions in the two layers are denoted by wd and w, respectively, and the corresponding relaxation time constants are given by τd and τs. We fix the units of time by setting τd = 1 (typical physical values of τd are 1–10 ms). We also include vertical feedforward input from deep to superficial layers with strength γd and feedback from the superficial to the deep layer with strength γs; in the latter case we average with respect to the orientation preference of neurons in the superficial layer. This is partly motivated by the observation that removing feedback does not appear to affect orientation tuning in cat layer 5 [45]. If this observation is combined with the assumption that recurrent connections in deep layers are not tuned for orientation, then the main source of orientation selectivity in deep layers is due to an external drive from an oriented stimulus. In the case of self-sustained propagating activity, such an input is absent so we can neglect any orientation dependence in the deep layer. Finally, the firing rate functions fj are taken to be sigmoidal, (2) where ηj is the gain and κj, κj > 0, is the threshold of the jth layer. In the high gain limit, ηj → ∞, we have fj(u) → H(uκj), where H is the Heaviside function.

thumbnail
Fig 3. Schematic illustration of the laminar neural field model.

Thalamic inputs feed into both the deep (d) and superficial (s) layers. At retinotopic position x, there is orientation independent activity u(x, t) in the deep layer and orientation-dependent activity v(x, θ, t) in the superficial layer. Intracortical connections in the deep layer are denoted by wd, whereas the intracortical connections w in the superficial can be decomposed into local wloc and long-range parts whoz. Vertical feedforward and feedback connections between the layers are also shown. Since we are interested in the self-sustained propagation of activity, we drop any external inputs Id and Is from thalamus or layer 4.

https://doi.org/10.1371/journal.pcbi.1004545.g003

In this paper, we primarily analyze traveling wave dynamics with a 1D network which has isotropic horizontal connections. However, we briefly present results for the more relevant 2D network with anisotropic horizontal connections. We therefore first decompose the recurrent weights of the superficial layer into short-range connections that depend only on the neurons’ orientation preferences, and long-range patchy horizontal connections between compatible neurons: (3) with wloc and whoz even, π-periodic functions of θ, and δ(x) the Dirac delta function. The θ-dependent part of the horizontal connections, whoz, is taken to be positive and narrowly peaked around zero so that only neurons with similar orientation preferences excite each other. Finally, ws(xx′), determines the strength of connections between two populations at x, x′ and is assumed to depend on the spatial separation ∣xx′∣. (For simplicity, we ignore any anisotropy in the horizontal connections.) Substituting the particular decomposition of w into Eqs (1a) and (1b) yields (4a) (4b) where x ∈ ℝ and One final aspect of our neural field model is that we do not explicitly model distinct excitatory and inhibitory populations. This is a common simplification of neural fields, in which the combined effects of excitation and inhibition are incorporated using, for example, Mexican hat functions [56, 58]. Typically, wave phenomena are modeled using excitatory neural fields, motivated by the fact that waves in vitro are observed in disinhibited cortical slices [912] and certain types of epilepsy appear to involve reductions in inhibition. On the other hand, the formation of population orientation tuning curves (stationary orientation bumps) in the ring model is modeled using a combination of short-range excitation and longer-range inhibition [36, 37]. Long-range horizontal connections are mediated by the axons of excitatory pyramidal neurons. However, they innervate both pyramidal neurons and feedforward interneurons so that they can effectively be treated as either excitatory or inhibitory, depending on stimulus conditions [35, 69, 70]. As a final level of complexity, there is a diversity of interneurons with different afferents and efferents distributed across different cortical layers [71]. In this paper we make the following simplifying assumptions with regards excitation and inhibition:

  1. Intracortical connections in the deep layer are excitatory so that wd(x) is a positive function of x. (One could use a Mexican hat function, provided that the inhibition is not too strong to prevent wave propagation.)
  2. Long-range horizontal connections in the superficial layer are excitatory so that ws(x) is a positive function of x.
  3. The local connections wloc(θ) in the superficial layer are excitatory for small θ (similar orientations) and inhibitory for larger θ (Mexican hat function).
  4. The vertical connections between the deep and superficial layers are excitatory so γs, γd > 0.

Results

One of the particularly useful aspect of continuum neural field models is that one can combine numerical simulations with analytical approaches adapted from the theory of partial differential equations [2]. In particular, the mathematical analysis allows us to identify parameter regimes where we expect traveling waves to occur, and to determine how the wave speed depends on these parameters. However, in order to carry out the mathematical analysis, one has to make some additional simplifying assumptions. Therefore, we proceed by initially considering a 1D version of our laminar model given by Eqs (4a) and (4b), in which feedback from the superficial to the deep layer is neglected (γs = 0). As a final simplification, we follow Amari [55] by taking the high gain limit, η → ∞, so that fj(u) → Hj(u) ≡ H(uκj), where H is the Heaviside function. Following our mathematical construction of traveling wave solutions, we show numerically that the resulting waves persist in the presence of feedback and for high gain smooth firing rate functions. Technical aspects of the analysis are presented in Materials and Methods, including a proof that feedback has no affect on the existence and speed of the traveling wave in the Heaviside case. Having developed our theory for 1D waves, we then show numerically how analogous waves occur in the full 2D network model. Note that throughout this section we ignore any external inputs (Is = Id = 0 in Fig 3), since we are interested in the self-sustained propagation of activity due to recurrent connections rather than activity driven explicitly by a propagating external input. Thus, we assume that there are only transient inputs, which initiate the wave by determining appropriate initial conditions.

thumbnail
Fig 4. Schematic illustration of two types of traveling wave.

(a) A traveling front. (b) A traveling pulse.

https://doi.org/10.1371/journal.pcbi.1004545.g004

We proceed with our analysis with feedback connections removed, γs = 0. However, we show in Materials and Methods that our results are unchanged when including feedback. In the absence of feedback, the deep layer is independent of activity in the superficial layer and Eq (4a) reduces to the classical scalar neural field equation (5)

The standard analysis of Amari [55] can then be used to establish the existence of a traveling wave solution in the deep layer. Formally speaking, a traveling wave is a solution that travels at constant velocity and fixed shape. For one-dimensional systems, one can distinguish two types of solitary traveling wave: a traveling front linking a stable down state to a stable up state, and a traveling pulse that begins and ends at the down state, see Fig 4. In this paper, we focus on traveling fronts. (In order to obtain traveling pulses, it would be necessary to include some form of local negative feedback such as spike frequency adaptation or synaptic depression [72].) Suppose that the excitatory weight distribution wd(x) is an even function of x and is monotonically decreasing for x ≥ 0. A common choice is the exponential weight distribution (6) where σd determines the range of synaptic connections. The latter tends to range from 100 μm to 1 mm. Mathematically speaking, a traveling front solution of Eq (5) takes the form u(x, t) = U(z) for some fixed profile function U, where z = xct is a traveling wave coordinate and c denotes the wavespeed. We also require that such that U(z) only crosses the threshold κd once. Since Eq (5) is translation invariant in space, we are free to take the threshold crossing point to be at the origin, U(0) = κd, so that U(z) < κd for z > 0 and U(z) > κd for z < 0. Substituting the traveling front solution into Eq (5) and requiring the resulting solution to be bounded then yields an implicit formula for the wave speed given by (7)

One of the useful aspects of the above analysis that it allows us to derive an explicit expression for the wavespeed as a function of physiological parameters such as firing threshold and range of synaptic connections. In the case of the exponential weight distribution (6), the relationship between wavespeed c and threshold κd is (8a) (8b) This establishes the existence of a unique front solution for fixed κd, which travels to the right (c > 0) when and travels to the left (c < 0) when . If the threshold κd is too large or the recurrent connections are too weak so that , then no traveling wave exists. It can be proven that the traveling front is stable [73]. Moreover, given the existence of a traveling front solution for a Heaviside rate function, it is possible to extend the analysis to a smooth sigmoid nonlinearity using a continuation method [74]. Finally, reintroducing the time constant τd, we see from Eq (8b) that the wavespeed has the order of magnitude cσd/τd = 10 cm/s for σd = 1 mm and τd = 10 msec. In the following we fix the units of space and time by setting σd = 1 = τd.

No spatially coherent orientation selective waves in the absence of horizontal connections (ws = 0)

Now suppose that the traveling front solution of the deep layer drives activity in the superficial layer via the vertical feedforward connections. For the moment, we ignore the effects of horizontal connections by setting ws = 0. (From our analysis of waves in the deep layer, we know that weak horizontal connections cannot on their own sustain propagating activity in the superficial layer.) Eq (4b) then reduces to (9) where we have used Hd(u(x, t)) = H(U(xct) − κd) = H(ctx) with c > 0. It is important to note that although Eq (9) and its generalization for ws > 0 suggests one could eliminate the deep layer and simply drive the superficial layer explicitly with a propagating external input, two important features are lost. First, one needs to understand the intrinsic cortical mechanisms that support wave propagation (in this case recurrent excitation in the deep layer) and determine the corresponding wavespeed and second, we will later include the effects of feedback from the superficial to the deep layer. (In Materials and Methods, we show that the inclusion of feedback does not affect the wavespeed c in the deep layer so that Eq (9) still holds.) The neural populations labeled by spatial position x act independently of one another so that we effectively have a continuum of independent orientation tuning equations (ring models) with time-dependent input γd H(ctx) for fixed x. Again, we can analyze each ring model using the standard analysis of Amari [55]. In spatial regions sufficiently far from the front at x = ct, we can assume that each ring network has reached a quasi-stationary state V(θ), which satisfies an equation of the form (10) with χ = 1 if xct and χ = 0 if xct. As shown by Amari [55], for certain choices of local recurrent connections, there exist stationary bump solutions, which in the present context can be interpreted as population-level orientation tuning curves. Conditions for the existence and stability of bumps can be determined in terms of the function In particular, the existence condition for a symmetric bump centered about θ = 0 is where Δ is the half bump width satisfying W(2Δ) = κsγd, and the stability condition is W′(2Δ) < 0. The existence and stability of bumps can then be analyzed graphically as illustrated in Fig 5 for wloc given by a cosine function: (11)

thumbnail
Fig 5. Existence of orientation bumps in the ring model.

(a) Plot of Mexican hat local connections, wloc = (-2 + 8 cos(2θ))/π. (b) Plot of W(2θ). The intercepts of W(2θ) with the horizontal line y = κsγ determine the bump solutions in the presence of a constant input γ.

https://doi.org/10.1371/journal.pcbi.1004545.g005

The latter realizes a basic assumption of the ring model [36, 37], namely, that neurons with similar orientation preferences excite each other, whereas those with sufficiently different orientation preferences inhibit each other.

In our model we assume that the superficial layer operates under the following conditions:

  1. Wmax < κs—condition for the down state V(θ) ≡ 0 to be the only stable steady state when χ = 0.
  2. There exists a unique stable bump solution Vbump(θ) of half-width Δ0 when χ = 1. This holds provided that Wmin < ksγd < Wmax.
  3. κsγd < 0—condition for the down-state not coexisting with the bump state when χ = 1.

Under these operating conditions, we then have the following scenario with regards to the response of the superficial layer to the propagating drive from the deep layer. Ahead of the front (x > ct) the superficial layer is in the down state. However, as soon as the front passes a particular location x, the down state disappears and the ring network at x evolves to a stationary bump solution at a relaxation rate determined by the time constant τs. One subtle feature is that there has to be some orientation-dependent perturbations in the superficial layer, otherwise the activated network evolves to an orientation-independent state that fluctuates around the threshold κs[75]. Such perturbations occur in the presence of horizontal connections, see below. However, since there are no horizontal connections linking distinct hypercolumns, each activated ring network forms a bump whose peak is uncorrelated with other bumps. Hence as activity propagates along the deep layer, connections to the superficial layer drive the activity to form orientation tuning curves with a random tuned orientation at each location x.

Spatially coherent orientation selective waves in the presence of weak horizontal connections (ws > 0)

We now establish that weak long-range horizontal connections serve to correlate the phase of the bump along the spatial direction leading to the propagation of a coherent orientation tuning curve. In Fig 6 we show numerical simulations of the full system given by Eqs (4a) and (4b), which includes feedback from the superficial layer to the deep layer, with the initial conditions where Vbump(θ) is the stationary bump solution of half-width Δ. We take the various weight distributions to have the particular form with parameter values specified in the caption of Fig 6. Note, however, that our results are insensitive to the precise choice of weight functions, provided they satisfy the various conditions highlighted in the text. In Fig 6, we plot the deep layer activity u(x, t) along side a density plot of the superficial layer activity v(x, θ, t) in the (x, θ) plane at successive times t for both a Heaviside and smooth firing rate function. Each snapshot clearly shows a sharp front separating an orientation bump centered about θ = 0 to the left of the front and an orientation-independent low activity state (down state) to the right of the front. Moreover, the front itself propagates to the right in successive frames. Additionally, we see that the results are qualitatively the same in the smooth firing rate case, except that the wave speed is larger. Now suppose that we simulate Eqs (4a) and (4b) as in Fig 6, and allow the solution to evolve for one time unit, t = 1. At this point, we remove feedforward connections from the deep to the superficial layer and allow the solution to evolve for another time unit. In Fig 7 we show a density plot of both the deep and superficial layer solutions in the (x, t) plane (we take maxθ v(x, θ, t)). We see that after the feedfoward connections are removed, the activity in the superficial layer decays towards the rest state, while the activity in the deep layer persists as a traveling wave. This implies that input from the deep layer to the superficial layer is necessary for wave propagation, which is consistent with the various experimental studies highlighted in the Introduction [19, 2628]. Furthermore, we see that the wave speed for the deep layer is approximately unchanged indicating that feedback connections to the deep layer are both unnecessary and negligible for wave propagation in the deep layer, which agrees with our analysis for Heaviside firing rates in Materials and Methods.

thumbnail
Fig 6. Numerical simulation of orientation-dependent traveling wave solution in the superficial layer.

(A) Plot of deep layer traveling wave solution for both Heaviside firing rate (grey) and smooth firing rate (black). (B) Orientation-dependent traveling wave solution in superficial layer for smooth firing rate function. We superimpose the threshold crossing contour for both Heaviside firing rate (white dashed) and smooth firing rate (white solid). In each time frame we show a contour plot of activity v(x, θ, t) in the (x, θ)-plane with high (superthreshold) activity indicated by orange (light) and low (subthreshold) activity indicated by brown (dark). Parameter values: wd,0 = 2, ws,0 = 0.1, w0 = -1, w2 = 1, kd = ks = 0.5, γd = γs = 1, σd = 0.1, σs = 1, τd = τs = 0.1, η = 10.

https://doi.org/10.1371/journal.pcbi.1004545.g006

thumbnail
Fig 7. Propagation failure in the superficial layer in the absence of feedfoward connections.

Density plot of traveling wave solutions in the (x, t) plane for deep and superficial layers before and after feedforward connections are removed at time t = 1 (indicated by the red line). For t > 1 the solution in the superficial layer decays to the rest state while the solution in the deep layer persists with the wave speed approximately unchanged indicating that the feedback connections play a negligible role in wave propogation. For comparison we superimpose a threshold crossing contour for both solutions with the Heaviside (dashed white line) and the smooth (solid white line) firing rate functions. Parameter values are the same as Fig 6.

https://doi.org/10.1371/journal.pcbi.1004545.g007

Inclusion of horizontal connections yields the following generalization of Eq (9): (12) with c independent of the feedback and determined solely by the deep layer. Unfortunately, it is not possible to obtain an explicit traveling wave solution of Eq (12). However, we observe that there are four distinct regions of interest: quiescent, queuing, integrating, bump. This is illustrated in Fig 8, where we superimpose these regions on a density plot of activity in the moving frame for the traveling wave solution. That is, we set v(x, θ, t) = V(z, θ) with z = xct and plot V(z, θ) in the (z, θ)-plane. This graphical construction can be used to qualitatively describe the traveling wave solution, as detailed in Materials and Methods. In order to define the different regions we introduce two quantities: the range L of long-range connections and the point z0 at which activity in the superficial layer first reaches the threshold κs, given that the front in the deep layer is located at z = 0. We find that if ws,0 is sufficiently small, so that the superficial layer does not support wave propagation on its own, then −L < z0 < 0 necessarily. The quiescent region is where neurons receive no synaptic inputs from either the deep layer or from active areas of the superficial layer, since they are out of range of the corresponding horizontal connections. Since the zero state is the only stable state when there is no input from the deep layer we find that V(z, θ) = 0 for z > L−∣z0∣ > 0. In the queuing region z0 < z < L−∣z0∣ neural populations have not yet received input from the deep layer, but have received a weak orientation dependent input from the active region via the horizontal connections. This leads to the formation of sub-threshold orientation bumps, which cannot cross threshold since the horizontal connections are weak. On the other hand, in the integrating region z0 < z < 0, neurons also receive input from the deep layer. Following our discussion on the existence of bumps, the input from the deep layer places the network in a parameter regime where the only stable solution is a bump. Hence the populations start integrating towards threshold as sub threshold bumps. Finally for z < z0, the populations continue to integrate towards the steady state bump solution Vbump(θ).

thumbnail
Fig 8. Graphical construction of a traveling wave solution.

The four regions quiescent, queuing, integrating, bump are superimposed on a traveling wave solution of the neural field Eq (12). Here L denotes the range of long-range connections and z0, z0 < 0 is the point at which activity in the superficial layer first reaches threshold. (The location of the front in the deep layer is at z = 0.) In the quiescent region (z > L − ∣z0∣) there is no activity. In the queuing region (0 < z < L − ∣z0∣) neural populations have not received input from the deep layer, but receive a weak orientation-dependent input from the active region via the horizontal connections. This leads to the formation of sub-threshold orientation bumps, which grow in the integrating region (z0 < z < 0) where populations now also receive orientation-independent input from the deep layer. Finally, the growing bumps cross threshold in the bump region (z < z0) evolving to a bump of half-width Δ as z → −∞. (We use a different color code than Figs 6 and 7 in order to highlight the different regions. Here red indicates high activity and blue low activity.)

https://doi.org/10.1371/journal.pcbi.1004545.g008

Traveling waves in 2D

As a starting point to understanding the dynamics of traveling waves of orientation selectivity in 2D, we numerically solve the full neural field Eqs (1a) and (1b). Analogous with the results in 1D we find that feedback connections do not affect the wave speed or the threshold crossings. Therefore, suppose that u(x, t) is a solution to Eq (1a) such that u(r, t) > κd for r = (x, y) ∈ U(t) ⊂ ℝ2 and u(r, t) < κd for rU(t), where U(t) is some set which changes in time. Hence where is the indicator function for the set U. The dynamics for the the superficial layer then evolve according to (13) Here we decompose the synaptic weight kernel as This is similar to the 1D case, except that we have included the function A(arg(r′−r)−θ) to account for anisotropies in the horizontal connections. Here arg(r) denotes the angle from the x-axis to the vector r and A is some positive, symmetric function, narrowly peaked about θ = 0, describing the strength of connections between neurons at locations r and r′ with common preferred orientation θ. That is, populations with similar orientation preferences in different hyper columns are only connected if they lie along the line with angle similar to their common orientation preference.

We would like to understand how anisotropic horizontal connections correlate the phase of the bump such that spatially structured patterns of orientation tuning emerge in the region U(t). Recently, we have studied this problem when the region U is fixed in time [75]. In this case we find that the solution inside the region U is approximately given by where V(θ) is the bump solution and ϕ(r, t) is the phase of the bump. We then derived an integro-differential equation describing the spatiotemporal dynamics of the phase of the bumps ϕ(r, τ) on the slow time scale τ = ɛt. We found that when the horizontal connections are taken to be isotropic, A ≡ 1, that the synchronous state ϕ(r) = ϕ0 is a steady state solution for any region U. However this is no longer necessarily the case when the connections are anisotropic, rather the geometry of the region and the structure of the anisotropy determine the solutions.

Our first example is the 2D analog of the 1D traveling wave solution namely, a plane wave. Note that the 2D neural field equation for the deep layer is and thus if (an example of such a function is a Gaussian) we can find a solution u = ue(x, t) which is constant in the y-direction, given by which supports the 1D traveling front solution constructed above. By rotational symmetry of the system, it then follows that u(r, t) = ue(Rϕ r, t) for any rotation by ϕ, Rϕ, is also a solution. This corresponds to a front which is constant along the vector k = (cosϕ, −sinϕ) and traveling in the direction k. One can then perform a similar simplification for the activity in the superficial layer. The simplest case is when the connections are isotropic, A ≡ 1, so that v(r, θ, t) = ve(x, θ, t) where ve is a solution to Eq (9b). In this case any phase in the bump yields a valid solution. However this is no longer the case with anisotropic connections due to the fact that integrating out one direction will not get rid of the dependence of wA on the orientation θ. Deriving an effective 1D equation for anisotropic connections is a non-trivial task so we leave this for future work. For now we numerically simulate the network with the full system in Eqs (4a) and (4b) with initial conditions We choose to set v = 0 in order to see what the phase of the bumps naturally settle to.

In Fig 9 we see that with anisotropic connections the phase is parallel to the boundary of the front. In other words, the orientation is orthogonal to the direction of wave propagation. However with isotropic connections the phase tends towards some synchronous solution entirely unrelated to the boundary or direction of propagation, which is consistent with our analysis of 1D waves. Simulating the network repeatedly with noise we find that this synchronous phase is random. Although isotropic horizontal connections cause the network to have a synchronous phase, without anisotropy there is no preference for what the phase of synchrony will be. Our numerical results suggest that, in the case of anisotropic connections, for U(t) = {r ∈ ℝ2 | x ≤ ct} we have that v(r, θ, t) = ve(x, θ, t) where ve is the 1D traveling bump solution, except that the center of the bump is at θ = π/2. It then follows from shift-twist invariance [60, 75] that for the rotated region the corresponding solution is given by v(r, θ, t) = ve(Rϕ x, θϕ, t), which has bump center at θ = π/2 + ϕ. Hence for waves in the deep layer traveling in the direction k = (sinϕ, cosϕ) the corresponding activity in the superficial layer is a traveling bump with orientation coinciding with the direction k.

thumbnail
Fig 9. Two-dimensional plane waves.

Numerical simulation of Eqs (4a) and (4b) for (A) anisotropic connection A(θ) = 1 + cos(2θ) and (B) isotropic connections A ≡ 1. All parameter values are the same as in Fig 6 and we take x0 = −0.5. All weight functions are the 2D analogs of the functions used for 1D.

https://doi.org/10.1371/journal.pcbi.1004545.g009

Our second and final example demonstrates the dynamics when the deep layer supports a traveling wave solution with a circular front. Such a solution is given by a target wave [7679]. In this case the region is given U(t) = {r ∈ ℝ2|0 ≤ |r| ≤ ct}. In the stationary case, U(t) = {r ∈ ℝ2|0 ≤ |r| ≤ R} for any R, we found that the angular solution ϕs(r) = arg(r) is a solution for both isotropic and anisotropic connections, whereas the synchronous solution only exists for isotropic connections. We find that the results are analogous in the traveling wave scenario. In Fig 10, we show a numerical simulation with initial conditions We see that with anisotropic connections the phase of the bumps tend towards the angular solution ϕ(r) = arg(r), suggesting that the anisotropic connections along with the geometry of the region determine the phase structure of the solution. On the other hand, with isotropic connections we see that there is no such structure in the phase. Rather, it appears that the solution tends towards a synchronous solution.

thumbnail
Fig 10. Two-dimensional target patterns.

Numerical simulation of Eqs (4a) and (4b) for (A) anisotropic connection A(θ) = 1 + cos(2θ) and (B) isotropic connections A ≡ 1. All parameter values are the same as in Fig 6 and we take r0 = 0.1. All weight functions are the 2D analogs of the functions used for 1D.

https://doi.org/10.1371/journal.pcbi.1004545.g010

Discussion

In this paper we have considered two important features of visual cortex that are usually neglected when developing neural field models of cortical waves, namely, the laminar structure of cortex and the orientation preferences of V1 cells in superficial layers. The latter naturally leads to the introduction of neural fields on product spaces [60, 61, 80]. It turns out that obtaining an orientation-dependent traveling front solution of the laminar model is non-trivial, requiring a delicate geometric construction. Two major predictions of our model are (i) long-range horizontal connections in superficial layers play a crucial role in supporting orientation-dependent propagating waves, and (ii) the waves originate in deep layers and are not strongly affected by feedback from superficial layers.

It is important to emphasize that we made a number of assumptions regarding the functional architecture of V1, which are still not settled experimentally and are likely to depend on the particular animal model. These assumptions include the following:

  1. Neurons in deep layers are weakly tuned for orientation. In particular, orientation tuning only occurs in response to oriented visual stimuli, but does not arise spontaneously due to a lack of orientation-dependent recurrent connections.
  2. Feedback from superficial to deep layers is orientation-independent.
  3. Traveling waves can occur in an isolated deep layer but not in an isolated superficial layer.
  4. There is a direct feedforward drive from the deep layer to the superficial layer, which is necessary in order to observe propagating activity in the superficial layer.
  5. Sharpening of orientation tuning occurs in superficial layers due to local recurrent excitation and lateral inhibition. Long-range connections are assumed to be excitatory.

It should be noted that advances in the use of viruses [81] means that it should eventually be possible to more accurately test these various assumptions by inactivating specific neuron populations in specific layers, and observing the effects on propgating activity. One possible application of the theory developed in this paper is to the study of binocular rivalry waves in V1 [13, 14, 82]. Binocular rivalry is the phenomenon where perception switches back and forth between different images presented to the two eyes; this switch does not occur simultaneously across the visual field but spreads like a propagating front. Previously, we developed a model of binocular rivalry waves consisting of two mutually-coupled 1D neural fields that were driven by left-eye and right-eye stimuli, respectively. Recurrent connections within each 1D network were assumed to be excitatory, whereas connections between the two networks were inhibitory (cross-inhibition). We showed that incorporating some form of slow adaptation such as synaptic depression into the model broke the symmetry between the left and right eye neural fields, thus allowing a front to propagate [15, 16]. Such a front represents the neural correlate of the propagating perceptual switch from one eye’s image to the other. One simplification of our previous work was to ignore the orientation preference of the neurons. However, it has been observed experimentally that the speed of binocular rivalry waves depends on the orientation of the left and right eye stimuli [13]. This suggests replacing the single layer left-eye and right-eye neural fields by a pair of laminar neural fields similar to the one considered in this paper.

Another possible extension of our work would be to consider other types of spatiotemporal patterns arising in a laminar neural field. For example, suppose that the deep-layer exhibited a Turing-like instability due to an increase in excitability, resulting in the spontaneous formation of a spatially periodic pattern. What would be the corresponding pattern of orientation-dependent in the superficial layer? This could provide a laminar-based extension of the theory of geometric visual hallucinations [60, 83].

Finally, it would be interesting to explore other aspects of how the functional architecture of cortex influences wave propagation. For example, by incorporating other feature preference maps such as ocular dominance (left/right eye preference) and spatial frequency [84] or variables associated with texture processing [85].) In doing so it might be necessary to distinguish between pinwheels and linear zones of the orientation preference map. It is important to note that the coupled ring model does not explicitly incorporate pinwheels. Hence, it cannot determine possible fine-scale differences between activity propagating near pinwheels and linear zones. On the other hand, the ring model itself can be interpreted as describing changes in orientation preference in linear zones as one rotates around a pinwheel.

Materials and Methods

We present the details of our analysis of an orientation-dependent traveling wave solution of the neural field Eq (12), based on the graphical construction shown in Fig 8. We then show how feedback does not modify the speed of a traveling wave in the deep layer. Finally, we briefly summarize our numerical methods.

Graphical construction of traveling wave in the superficial layer

We seek a solution of Eq (12) of the form v(x, θ, t) = V(z, θ) with z = xct and (14) where Vbump satisfies for with wloc(θ) + ws,0 whoz(θ) satisfying conditions (i)-(iii) below Eq (11). Substituting the wave solution into Eq (12) gives (15) where we have absorbed the τs term into c for notational clarity. If such a solution exists then there must be a point z0 < 0 such that V(z, θ) < κs for all (z, θ) ∈ (z0, ∞) × [−π/2, π/2]. That is there must be a point where the solution first reaches threshold. Furthermore this point must be negative since γH(−z) = 0 for z > 0 and hence the dynamics remains below threshold in this region, as we show more explicitly below. We then observe that for z > z0 (16) and for z < z0 (17) where (18) Therefore, the region z > z0 is only dependent on the activity in region z < z0 while the active region receives no input from the inactive region and is thus not dependent on the solution in the inactive region. In the following we let V+(z, θ) denote the solution for z < z0 and V(z, θ) denote the solution for z > z0.

Using the existence of the point z0, we now determine the quiescent, queuing, integrating, and bump regions shown in Fig 8. In order to facilitate the analysis, we assume that ws has compact support, say such that ∣z0∣ < L. Then for z > L−∣z0∣ we have that and thus (19)

Hence, the requirement that V(z, θ) → 0 as z → ∞ translates to V(z, θ) ≡ 0 for all z > L − ∣z0∣ and we call this the quiescent region.

For 0 < z < L − ∣z0

In this queuing region, the populations have not received input from the deep layer but receives a weak synaptic input from the active region. If V+(z, θ) have the same bump center for all z, then V(z, θ) will also have the same center and will remain below threshold. More precisely, suppose that there is a unique θ0 such that for all z < 0. Then since V+(z, θ) must be symmetric about θ0 we have that V + (z, θ) > κs for θ ∈ (θ0 − Δ(z), θ0+Δ(z)) for each z > 0 and some function, Δ(z) describing the z-dependent bump half width. Hence multiplying both sides by and integrating we obtain, using the initial condition V(L−∣z0∣, θ) = 0, where (20) with Thus the location of the max of V(z, θ), with respect to θ, is determined by the location of the max of Ωhoz(θθ0, Δ(z)), which is θ0. Furthermore, in this region the activity remains below threshold for sufficiently small ws. To see this it is simple to show that for all 0 < z < L−∣z0∣. Hence, as long as ws,0 whoz,0 < ks then the solution remains below threshold. Thus, before the populations receive input, they are queued up as sub-threshold bumps to start integrating to super-threshold bumps with the appropriate center. This is the most important region with respect to propagation and is where horizontal connections play a crucial role. Without horizontal connections there is no queuing, and the solution remains at V(z, θ) = 0 until z = 0. Once the populations receive input, they start integrating towards the threshold state (unless perturbed by orientation-dependent noise), as described for the case ws ≡ 0. However, with horizontal connections the super-threshold activity forms an orientation-dependent sub-threshold bump which provides a mechanism for both the destabilization of the threshold state as well as correlating the bumps in a way that the bump center remains the same as it propagates.

In the integrating region z0 < z < 0, we have

The populations now receive orientation-independent input from the deep layer and a weak orientation-dependent input from horizontal connections. Since the zero state no longer exists the resulting solution evolves as a subthreshold bump state until it reaches threshold at z = z0. Finally, in the bump region z < z0, the remaining dynamics are governed by the equation (21)

Here, the solution integrates beyond threshold and continues to evolve towards the stationary bump.

Feedback has no affect on wavespeed in the deep layer

We show that, in the case of Heaviside firing rate functions, the feedback connections do not affect the wave speed, rather only the shape of the wave. Traveling wave solutions in the deep layer with feedback satisfy Following the analysis of bump solutions in Results, see also Fig 5, we look for traveling wave solutions U(z, θ) and V(z, θ) such that and Note that the asymptotic condition for the up state in the deep layer now includes a feedback term 2Δγs where Δ is the half-width of the bump solution Vbump. Assuming that V(z, θ) is a traveling wave, there must exist some point z0 such that the solution first reaches threshold. Moreover, if ws,0 is sufficiently small, so as to not support propagation on its own, then z0 < 0 necessarily (see Fig 11). That is, the threshold crossing for the superficial layer lags behind the deep layer, and beyond this point the solution integrates towards a bump solution. Thus at each z < z0 there is a function describing the instantaneous bump width, Δ(z) with the conditions

thumbnail
Fig 11. Schematic illustration of superficial layer traveling front solution V(z, θ) in the z-θ plane.

The curve determined by ±Δ(z) separates the plane into a region of super threshold activity (grey) and sub threshold activity.

https://doi.org/10.1371/journal.pcbi.1004545.g011

This function describes the threshold crossing contour, i.e. V(±Δ(z), z) = κs. We can then write Hs(V(z, θ)) = H(Δ(z)−∣θ∣)H(z0z). Plugging this form into the integral yields

Therefore, when carrying out the computation for the positive wave speed as above we obtain where the integration with the Heaviside term vanishes since it is identically zero for z0 < 0. Therefore the wave speed remains unchanged when adding feedback connections. This can be also understood by viewing neural fields in terms of interface dynamics as in [86], where they show that with a Heaviside firing rate function, the evolution of solutions only depend on the activity at the interface (or the threshold crossings). Calculating the solution yields where U0(z) is the traveling wave solution in the absence of feedback. Since Δ(z) ≥ 0 and γs > 0 the second integral term is positive for z < z0 and identically zero for zz0. Therefore the threshold crossing for U0(z) and U(z, θ) coincide and thus Hd(U(z)) = Hd(U0(z)) = H(−z). Hence, the dynamics in the superficial layer become which is precisely what one obtains when ignoring feedback, i.e. setting γs = 0. A mathematically convenient consequence of this fact is that we can first prove existence of traveling wave solutions in the deep layer, in the absence of feedback, independently of the dynamics in the superficial layer. We can then use this solution to provide an effective time-dependent input to the superficial layer, and analyze its dynamics independently of the deep layer. Note, however, that this is purely an artifact of the Heaviside firing rate and in general may not apply to arbitrary firing rate functions. However, in the Results, we show numerically that the wave properties persist when using a smooth firing rate with sufficiently large gain, η.

Numerical methods

All numerical simulations were performed in Matlab. One dimensional numerical simulations were performed using a forward Euler method scheme in time and a trapezoidal rule for integration in x and θ. Time steps were taken to be Δt = 0.001, spatial steps Δx = 0.01 and orientation steps Δθ = 0.01π. The numerical methods employed for 2D simulations were similar. However, rather than using the trapezoidal rule for integration in x, we used Fast Fourier Transforms to compute the convolutions. Recall that the convolution theorem states where f*g is the convolution of f and g and is the Fourier transform of f. We can therefore write the neural field equation as We therefore compute integration in θ using a trapezoidal rule for each fixed x. We then take the result and perform fast Fourier transforms to compute the remaining convolution in x. Time steps were taken as Δt = 0.001, spatial steps Δx = 0.02 and orientation steps Δθ = 0.01π. To plot the vector fields we find a ϕ(x, t) such that at each x, t, which denotes the tuned orientation of the population. We then define the vector field as

Author Contributions

Wrote the paper: PCB SRC. Model: PCB. Analysis of Model: SRC PCB. Numerical Simulations: SRC.

References

  1. 1. Muller L, Destexhe A. Propagating waves in thalamus, cortex and the thalamocortical system: experiments and models. J Physiol (Paris). 2012;106:222–238.
  2. 2. Bressloff PC. Waves in Neural Media: from single neurons to neural fields. Springer; 2014.
  3. 3. Milton J, Jung P. Epilepsy as a dynamic disease. Springer; 2003.
  4. 4. Soltesz I, Staley K. Computational neuroscience in epilepsy. San Diego: Academic Press; 2008.
  5. 5. Kramer MA, Cash SS. Epilepsy as a disorder of cortical network organization. Neuroscientist. 2012;18:360–372. pmid:22235060
  6. 6. Martins-Ferreira H, Nedergaard M, Nicholson C. Perspectives on spreading depression. Brain Res Rev. 2000;32:215–234. pmid:10751672
  7. 7. Somjen GG. Mechanisms of spreading depression and hypoxic spreading depression-like depolarization. Physiol Rev. 2001;81:1065–1096. pmid:11427692
  8. 8. Dahlem MA, Chronicle EP. A computational perspective on migraine aura. Prog Neurobiol. 2004;74:351–361. pmid:15649581
  9. 9. Chervin RD, Pierce PA, Connors BW. Periodicity and directionality in the propagation of epileptiform discharges across neocortex. J Neurophysiol. 1988;60:1695–1713. pmid:3143812
  10. 10. Golomb D, Amitai Y. Propagating neuronal discharges in neocortical slices: Computational and experimental study. J Neurophysiol. 1997;78:1199–1211. pmid:9310412
  11. 11. Pinto D, Patrick SL, Huang WC, Connors BW. Initiation, propagation, and termination of epileptiform activity in rodent neocortex in vitro involve distinct mechanisms. J Neurosci. 2005;25:8131–8140. pmid:16148221
  12. 12. Richardson KA, Schiff SJ, Gluckman BJ. Control of traveling waves in the mammalian cortex. Phys Rev Lett. 2005;94:028103. pmid:15698234
  13. 13. Wilson HR, Blake R, Lee SH. Dynamics of traveling waves in visual perception. Nature. 2001;412:907–910. pmid:11528478
  14. 14. Lee SH, Blake R, Heeger DJ. Traveling waves of activity in primary visual cortex during binocular rivalry. Nat Neurosci. 2005;8:22–23. pmid:15580269
  15. 15. Bressloff PC, Webber M. Neural field model of binocular rivalry waves. J Comput Neurosci. 2012;32:233–252. pmid:21748526
  16. 16. Carroll SR, Bressloff PC. Binocular rivalry waves in directionally selective neural field models. Physica D. 2014;285:8–17.
  17. 17. Xu W, Huang X, Takagaki K, Wu JY. Compression and reflection of visually evoked cortical waves. Neuron. 2007;55:119–129. pmid:17610821
  18. 18. Han F, Caporale N, Dan Y. Reverberation of recent visual experience in spontaneous cortical waves. Neuron. 2008;60:321–327. pmid:18957223
  19. 19. Stroh A, Adelsberger H, Groh A, Rhülmann C, Fischer S, Schierloh A, et al. Making Waves: Initiation and Propagation of Corticothalamic Ca2+ Waves In Vivo. Neuron. 2013;77:1136–1150. pmid:23522048
  20. 20. Smith GB, Sederberg A, Elyada YM, Hooser SDV, Kaschube M, Fitzpatrick D. The development of cortical circuits for motion discrimination. Nat Neurosci. 2015;18:252–261. pmid:25599224
  21. 21. Benucci A, Frazor R, Carandini M. Standing waves and traveling waves distinguish two circuits in visual cortex. Neuron. 2007;55:103–17. pmid:17610820
  22. 22. Chavane F, Sharon D, Jancke D, Marre O, Fregnac Y, Grimvald A. Lateral spread of orientation slectivity in V1 is controlled by intracortical cooperativity. Front Syst Neurosci. 2011;5(4):1–26.
  23. 23. Nauhaus I, Busse L, Ringach DL, Carandini M. Robustness of traveling waves in ongoing activity of visual cortex. J Neurosci. 2012;32:3088–3094. pmid:22378881
  24. 24. Grinvald A, Lieke EE, Frostig RD, Hildesheim R. Cortical point-spread function and long-range lateral interactions revealed by real-time optical imaging of macaque monkey primary visual cortex. J Neurosci. 1994;14:2545–2568. pmid:8182427
  25. 25. Sato TK, Nauhaus I, Carandini M. Traveling waves in visual cortex. Neuron. 2012;75:218–229. pmid:22841308
  26. 26. Wester JC, Contreras D. Columnar interactions determine horizontal propagation of recurrent network activity in neocortex. J Neurosci. 2012;32:5454–5471. pmid:22514308
  27. 27. Plomp G, Quairiaux C, Kiss JZ, Astolfi L, Michel CM. Dynamic connectivity among cortical layers in local and large-scale sensory processing. Eur J Neurosci. 2014;40:3215–3223. pmid:25145779
  28. 28. Krause BM, Raz A, Uhlrich DJ, Smith PH, Banks MI. Spiking in auditory cortex following thalamic simulation is dominated by cortical network activity. Front Syst Neurosci. 2014;8:(170).
  29. 29. Blasdel GG, Salama G. Voltage-sensitive dyes reveal a modular organization in monkey striate cortex. Nature. 1986;321:579–585. pmid:3713842
  30. 30. Bonhoeffer T, Grinvald A. Orientation columns in cat are organized in pinwheel like patterns. Nature. 1991;364:166–146.
  31. 31. Blasdel GG. Orientation selectivity, preference, and continuity in monkey striate cortex. J Neurosci. 1992;12:3139–3161. pmid:1322982
  32. 32. Malach R, Harel YAM, Grinvald A. Relationship between intrinsic connections and functional architecture revealed by optical imaging and in vivo targeted biocytin injections in primate striate cortex. Proc Natl Acad Sci USA. 1993;90:0469–10473.
  33. 33. Yoshioka T, Blasdel GG, Levitt JB, Lund JS. Relation between patterns of intrinsic lateral connectivity, ocular dominance, and cytochrome oxidase—reactive regions in macaque monkey striate cortex. Cerebral Cortex. 1996;6:297–310. pmid:8670658
  34. 34. Bosking WH, Zhang Y, Schofield B, Fitzpatrick D. Orientation selectivity and the arrangement of horizontal connections in tree shrew striate cortex. J Neurosci. 1997;17:2112–2127. pmid:9045738
  35. 35. Angelucci A, Levitt JB, Walton EJS, Hupe JM, Bullier J, Lund JS. Circuits for local and global signal integration in primary visual cortex. J Neurosci. 2002;22:8633–8646. pmid:12351737
  36. 36. Somers DC, Nelson S, Sur M. An Emergent Model of Orientation Selectivity in Cat Visual Cortical Simple Cells. J Neurosci. 1995;15:5448–5465. pmid:7643194
  37. 37. Ben-Yishai R, Bar-Or RL, Sompolinsky H. Theory of orientation tuning in visual cortex. Proc Nat Acad Sci. 1995;92:3844–3848. pmid:7731993
  38. 38. Bressloff PC, Cowan JD. Amplitude equation approach to contextual effects in visual cortex. Neural Comput. 2002;14:493–525. pmid:11860680
  39. 39. Priebe NJ, Ferster D. Direction selectivity of excitation and inhibition in simple cells of the cat primary visual cortex. Neuron. 2005;45:1133–1145.
  40. 40. Priebe NJ, Ferster D. Mechanisms underlying cross-orientation suppression in cat visual cortex. Nature Neuroscience. 2006;9:552–561. pmid:16520737
  41. 41. Niell CM, Stryker MP. Highly Selective Receptive Fields in Mouse Visual Cortex. J Neurosci. 2008;28:7520–7536. pmid:18650330
  42. 42. Hooser SDV, Roy A, Rhodes HJ, Culp JH, Fitzpatrick D. Transformation of Receptive Field Properties from Lateral Geniculate Nucleus to Superficial V1 in the Tree Shrew. J Neurosci. 2013;33:11494–11505. pmid:23843520
  43. 43. Ringach DL, Shapley RM, Hawken MJ. Orientation Selectivity in Macaque V1: Diversity and Laminar Dependence. J Neurosci. 2002;22:5639–5651. pmid:12097515
  44. 44. Shushruth S, Nurminen L, Bijanzadeh M, Ichida JM, Vanni S, Angelucci A. Different Orientation Tuning of Near- and Far-Surround Suppression in Macaque Primary Visual Cortex Mirrors Their Tuning in Human Perception. J Neurosci. 2013;33:106–119. pmid:23283326
  45. 45. Schwark HD, Malpeli JG, Weyand TG, Lee C. Cat area 17. II. Response properties of infragranular layer neurons in the absence of supragranular layer activity. J Neurophysiol. 1986;56:1074–1087. pmid:3783230
  46. 46. Raizada RDS, Grossberg S. Towards a theory of the laminar architecture of cerebral cortex: computational clues from the visual system. Cereb Cortex. 2003;13:100–113. pmid:12466221
  47. 47. Grossberg S. How does the cerebral cortex work? Development, learning, attention, and 3D vision by laminar circuits of visual cortex. Behav Cogn Neur Rev. 2003;2:47–76.
  48. 48. Grossberg S, Versace M. Spikes, synchrony, and attentive learning by laminar thalamocortical circuits. Brain research. 2008;1218:278–312. pmid:18533136
  49. 49. Kilpatrick ZP. Interareal coupling reduces encoding variability in multi-area models of spatial working memory. Frontiers in computational neuroscience. 2013;7.
  50. 50. Kilpatrick ZP. Coupling layers regularizes wave propagation in stochastic neural fields. Physical Review E. 2014;89:022706.
  51. 51. Bressloff PC, Kilpatrick ZP. Nonlinear Langevin equations for wandering patterns in stochastic neural fields. SIAM J Appl Dyn Syst. 2015;14:305–334.
  52. 52. Folias SE, Ermentrout GB. Bifurcations of stationary solutions in an interacting pair of E-I neural fields. SIAM J Appl Dyn Syst. 2012;11:895–938.
  53. 53. Wilson HR, Cowan JD. Excitatory and inhibitory interactions in localized populations of model neurons. Biophys J. 1972;12:1–23. pmid:4332108
  54. 54. Wilson HR, Cowan JD. A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik. 1973;13:55–80. pmid:4767470
  55. 55. Amari S. Dynamics of pattern formation in lateral inhibition type neural fields. Biol Cybern. 1977;27:77–87. pmid:911931
  56. 56. Ermentrout GB. Neural networks as spatio-temporal pattern-forming systems. Rep Prog Phy. 1998;61:353–430.
  57. 57. Coombes S. Waves, bumps and patterns in neural field theories. Biol Cybern. 2005;93:91–108. pmid:16059785
  58. 58. Bressloff PC. Spatiotemporal dynamics of continuum neural fields: Invited Topical review. J Phys A. 2012;45:033001 (109 pp.).
  59. 59. Hubel D, Wiesel T. Functional architecture of macaque monkey visual cortex. Proc R Soc London Ser B. 1977;198:1–59.
  60. 60. Bressloff PC, Cowan JD, Golubitsky M, Thomas PJ, Wiener M. Geometric Visual Hallucinations, Euclidean Symmetry and the Functional Architecture of Striate Cortex. Phil Trans Roy Soc Lond B. 2001;356:299–330.
  61. 61. Bressloff PC, Cowan JD, Golubitsky M, Thomas PJ. Scalar and pseudoscalar bifurcations: pattern formation on the visual cortex. Nonlinearity. 2001;14:739–775.
  62. 62. Petitot J. The neurogeometry of pinwheels as a sub-Riemannian contact structure. J Physiol Paris. 2003;97:265–309. pmid:14766146
  63. 63. Ben-Shahar O, Zucker SW. Geometrical Computations Explain Projection Patterns of Long Range Horizontal Connections in Visual Cortex. Neural Computation. 2004;16:445–476. pmid:15006089
  64. 64. Callaway EM. Local circuits in primary visual cortex of the macaque monkey. Ann Rev Neurosci. 1998;21:47–74. pmid:9530491
  65. 65. Hirsch JA, Martinez LM. Laminar processing in the visual cortical column. Curr Opin Neurobiol. 2006;16:377–384. pmid:16842989
  66. 66. Casagrande VA, Yazar F, Jones KD, Ding Y. The morphology of the koniocellular axon pathway in the macaque monkey. Cereb Cortex. 2007;17:2334–2345. pmid:17215477
  67. 67. Constantinople CM, Bruno RM. Deep cortical layers are activated directly by thalamus. Science. 2013;340:1591–1594. pmid:23812718
  68. 68. Douglas RJ, Martin KAC. A functional microcircuit for cat visual cortex. J Physiol. 1991;440:735–769. pmid:1666655
  69. 69. Lund JS, Angelucci A, Bressloff PC. Anatomical substrates for functional columns in macaque monkey primary visual cortex. Cerebral Cortex. 2003;12:15–24.
  70. 70. Schwabe L, Obermayer K, Angelucci A, Bressloff PC. The role of feedback in shaping the extra-classical receptive field of cortical neurons: a recurrent network model. J Neurosci. 2006;26:9117–9129. pmid:16957068
  71. 71. Neske GT, Patrick SL, Connors BW. Contributions of Diverse Excitatory and Inhibitory Neurons to Recurrent Network Activity in Cerebral Cortex. J Neurosci. 2015;35:1089–1105. pmid:25609625
  72. 72. Pinto D, Ermentrout GB. Spatially structured activity in synaptically coupled neuronal networks: I. Traveling fronts and pulses. SIAM J Appl Math. 2001;62:206–225.
  73. 73. Coombes S, Owen MR. Evans functions for integral neural field equations with Heaviside firing rate function. SIAM J Appl Dyn Syst. 2004;4:574–600.
  74. 74. Ermentrout GB, McLeod JB. Existence and uniqueness of travelling waves for a neural network. Proc Roy Soc Edin A. 1993;123:461–478.
  75. 75. Carroll SR, Bressloff PC. Phase equation for patterns of orientation selectivity in a neural field model of visual cortex. SIAM J Appl Dyn Syst. 2015; Accepted.
  76. 76. Laing CR. Spiral waves in nonlocal equations. SIAM J Appl Dyn. 2005;4:588–606.
  77. 77. Troy WC, Shusterman V. Patterns and Features of Families of Traveling Waves in Large-Scale Neuronal Networks. SIAM J Appl Dyn Syst. 2007;6:263–292.
  78. 78. Shusterman V, Troy WC. From baseline to epileptiform activity: a path to synchronized rhythmicity in large-scale neural networks. Phys Rev E. 2008;77:061911.
  79. 79. Kilpatrick ZP, Bressloff PC. Spatially structured oscillations in a two-dimensional neuronal network with synaptic depression. J Comp Neurosci. 2010;28:193–209.
  80. 80. Bressloff PC, Carroll SR. Spatiotemporal dynamics of neural fields on product spaces. SIAM J Appl Dyns Syst. 2014;13:1620–1653.
  81. 81. Nielsen KJ, Callaway EM, Krauzlis RJ. Viral vector-based reversible neuronal inactivation and behavioral manipulation in the macaque monkey. Frontiers in Systems Neuroscience. 2012;6:(48). pmid:22723770
  82. 82. Kang M, Heeger DJ, Blake R. Periodic perturbations producing phase-locked fluctuations in visual perception. J Vision. 2009;9(2):8.1–12.
  83. 83. Ermentrout GB, Cowan J. A mathematical theory of visual hallucination patterns. Bio Cybern. 1979;34:137–150.
  84. 84. Bressloff PC. Spatially periodic modulation of cortical patterns by long-range horizontal connections. Physica D. 2003;185:131–157.
  85. 85. Chossat P, Faugeras O. Hyperbolic planforms in relation to visual edges and textures perception PLoS Comput. Biol. 2009;5:e1000625.
  86. 86. Coombes S, Schmidt H, Bojak I. Interface dynamics in planar neural field models. J Math Neurosci. 2012;2(1):1–27.