Abstract
While our understanding of the way single neurons process chromatic stimuli in the early visual pathway has advanced significantly in recent years, we do not yet know how these cells interact to form stable representations of hue. Drawing on physiological studies, we offer a dynamical model of how the primary visual cortex tunes for color, hinged on intracortical interactions and emergent network effects. After detailing the evolution of network activity through analytical and numerical approaches, we discuss the effects of the model’s cortical parameters on the selectivity of the tuning curves. In particular, we explore the role of the model’s thresholding nonlinearity in enhancing hue selectivity by expanding the region of stability, allowing for the precise encoding of chromatic stimuli in early vision. Finally, in the absence of a stimulus, the model is capable of explaining hallucinatory color perception via a Turing-like mechanism of biological pattern formation.
Author Summary
We present a model of color processing in which intracortical neuronal dynamics within the visual cortex serve as the substrate for hue perception. Our analytical and numerical treatments of the emergent behavior seek to characterize the population dynamics underlying chromatic processing within the visual cortex, as well the roles of the various cortical parameters in determining the selectivity of the steady-state network response. We show that the system is self-organizing, capable of encoding stable representations of hue regardless of the stimulus strength, and generating spontaneous color hallucinations in the absence of any input.
INTRODUCTION
Our experience of color begins in the early visual pathway, where, from the moment light strikes the retina, cone-specific neuronal responses set off the mechanisms by which the photons’ chromatic information is converted to the hues we ultimately see. While color vision scientists agree that the single-cell processing of chromatic stimuli occurs along the two independent cone-opponent L–M and S–(L+M) pathways (Conway, Eskew, Martin, & Stockman, 2018; Kaiser & Boynton, 1996), there is yet no consensus as to how and where the divergent signals are synthesized to encode hue. To complicate matters, cone-opponency, observed in electrophysiological recordings of single neurons (Shapley & Hawken, 2011), is often confounded with hue-opponency, a phenomenon first theorized by Ewald Hering in the 19th century and later mapped out in clinical studies by Jameson and Hurvich (De Valois, Cottaris, Elfar, Mahon, & Wilson, 2000; Jameson & Hurvich, 1955; Shevell & Martin, 2017).
Best depicted in the Derrington-Krauskopf-Lennie (DKL) stimulus space (Figure 1), cone-opponency predicts that neurons tuned to either the L–M or S–(L+M) pathway will not respond to light whose wavelengths isolate the other (Derrington, Krauskopf, & Lennie, 1984). It is tempting to equate these null responses to the four unique hues of color-opponent theory, in which unique blue, for example, is observed when the “redness” and “greenness” of a perceived color exactly cancel. But the wavelengths of the unique hues specified by perceptual studies (Jameson & Hurvich, 1955) only roughly match the wavelengths isolating either cone-opponent pathway (Wool et al., 2015; Wuerger, Atkinson, & Cropper, 2005; Xiao, 2014), and, more fundamentally, we do not yet understand the mechanisms behind the processing that the analogy implies (Mollon & Jordan, 1997; Stoughton & Conway, 2008; Valberg, 2001). That is, how do we get from the single neurons’ chromatic responses to our perception of color?
The DKL space maps chromatic stimuli onto a circle with two “cardinal” axes representing the L–M and S–(L+M) pathways. The excitatory or inhibitory effect of a stimulus on cone-opponent cells tuned to either pathway can be thought of as a projection of its location in DKL space onto the relevant axis. Stimuli isolating the two pathways correspond roughly to wavelengths associated with the red, green, blue, and yellow unique hues of color-opponent theory, leading to the common, but mistaken, conflation of chromatic stimulus and color perception.
The DKL space maps chromatic stimuli onto a circle with two “cardinal” axes representing the L–M and S–(L+M) pathways. The excitatory or inhibitory effect of a stimulus on cone-opponent cells tuned to either pathway can be thought of as a projection of its location in DKL space onto the relevant axis. Stimuli isolating the two pathways correspond roughly to wavelengths associated with the red, green, blue, and yellow unique hues of color-opponent theory, leading to the common, but mistaken, conflation of chromatic stimulus and color perception.
The necessary processing has often been attributed to higher-level brain function (De Valois & De Valois, 1993; Lennie, Krauskopf, & Sclar, 1990; M. Li, Liu, Juusola, & Tang, 2014; Mehrani, Mouraviev, & Tsotsos, 2020; Zaidi & Conway, 2019) or yet unidentified higher order mechanisms (Valberg, 2001; Wuerger et al., 2005). A central question of color vision research is whether these mechanisms rely on parallel or modular processing to encode stimulus chromaticity (Conway, 2009; Garg, Li, Rashid, & Callaway, 2019; Liu et al., 2020; Nauhaus, Nielsen, Disney, & Callaway, 2012; Schluppeck & Engel, 2002; Shapley & Hawken, 2011). If signaling about chromaticity is transmitted with information about other visual features, such as brightness, orientation, and spatial frequency, how do these features get teased apart? If not, where is the purported color center of the brain (Conway, Moeller, & Tsao, 2007; Gegenfurtner, 2003)?
Several authors have addressed these questions through combinatorial models that parameterize the weights of the L, M, and S cones contributing to successive stages of processing (De Valois & De Valois, 1993; Gegenfurtner & Ennis, 2015; Judd, 1949; Mehrani et al., 2020; Stockman & Brainard, 2010). Though differing in their assumptions of modularity, the theories share a mechanistic framework for the transition of single-cell receptive field properties (Brown, 2014). Starting with cells in the retina and lateral geniculate nucleus (LGN) known to be tuned broadly to the cone-opponent axes, these proposed mechanisms build up to cells in various cortical areas more narrowly tuned to divergent (and debated) chromatic directions in DKL space. While parsimonious, this formalism comes at the cost of tuning the cone weights arbitrarily, disregarding specific properties of real neurons’ receptive fields (Eskew, 2009; Kaiser & Boynton, 1996; Stockman & Brainard, 2010). Furthermore, the linear combinatorial mechanism is not, on its own, able to account for the variety of color cells observed in the visual cortex (Garg et al., 2019; Johnson, Hawken, & Shapley, 2001; Shapley & Hawken, 2011). In addition to the forward flow of chromatic information through the successive stages of processing, the encoding of color reflects the neuronal dynamics within each. Modelers agree that the next forays into a mechanistic theory of color vision should consider these intracortical circuits, but disagree about where such interactions first become important (De Valois & De Valois, 1993; Hanazawa, Komatsu, & Murakami, 2000; Liu et al., 2020; Wachtler, Sejnowski, & Albright, 2003).
Electrophysiological studies of macaque visual cortex have shed some light on this question, showing that the processing of individual hues previously associated with higher level mechanisms has its origins in the primary visual cortex (V1) (Garg et al., 2019; Gegenfurtner, 2003; Hanazawa et al., 2000; Li et al., 2022; Wachtler et al., 2003; Xiao, 2014; Xiao, Casti, Xiao, & Kaplan, 2007). These experiments have identified the emergence of neurons in V1 tuned to the gamut of hues in DKL space, as well as to the role of processing nonlinearities in determining their tuning curves (De Valois et al., 2000; Hanazawa et al., 2000; Lennie et al., 1990; Wachtler et al., 2003). Puzzlingly, these cells mainly inhabit the so-called CO “blobs,” patchy regions rich in cytochrome oxidase that display a sensitivity to stimuli modulating either of the cone-opponent axes rather than the full set of hues (Landisman & Ts’o, 2002b; Li et al., 2022; Livingstone & Hubel, 1984; Salzmann, Bartels, Logothetis, & Schüz, 2012). Some have speculated that this colocalization stems from a mixing of cell populations encoding the two cardinal pathways (Li et al., 2022; Xiao, 2014) while others indicate a distinct population of hue-sensitive neurons in the “interblob” regions, more conclusively associated with orientation tuning (Garg et al., 2019; Landisman & Ts’o, 2002a). As a whole, however, these studies point to the need for a population theory of chromatic processing remarkably early in the visual pathway.
In this article, we present a model of color processing in which intracortical neuronal dynamics within V1 serve as the substrate for hue perception. Drawing on the canonical Wilson-Cowan neural field equations and the ring model of orientation tuning, we show that this population approach allows us to account for cells responsive to the full range of DKL directions without the need to fine-tune input parameters (Ben-Yishai, Bar-Or, & Sompolinsky, 1995; Burstein, 2022; Hansel & Sompolinksy, 1998; Wilson & Cowan, 1972, 1973). The threshholding we employ bears in mind the input-response nonlinearities of previous combinatorial models, but zooms out of the single-cell, feedforward interpretation of input as the stimulus-driven LGN afferents to individual neurons. Rather, we model input as the total synaptic current into a population of cells, taking into account both the cone-opponent LGN afferents as well as the hue-dependent connectivity between distinct neuronal populations.
The resulting demarcation between the cone-opponent and the hue-selective mechanisms in the same population of cells points to the importance of V1 in the transition from chromatic stimulus to color perception. To characterize this role, we study the effects of the model’s connectivity parameters and processing nonlinearities on the narrowness and stability of the hue tuning curves. In the final part of the paper, we show that the model is able to explain color responses in the absence of LGN input, evoking color hallucinations via a Turing-like mechanism of spontaneous pattern formation in DKL space.
MODEL
In light of the patchy distribution of color-sensitive cells reported in Landisman and Ts’o (2002b), Li et al. (2022), Livingstone and Hubel (1984), and Salzmann et al. (2012), we model the color map of V1 as a set of neuronal networks, each encoding the chromaticity of its corresponding region of the visual field. This organization brings to mind the hypercolumnar structure of orientation preference within V1 (Hubel & Wiesel, 1974), which, on the basis of its feature-based connectivity properties, allows for the representation of network activity as a function of a localized feature space. Here, we assume a mean hue-dependent activity a(θ, t) where θ represents a direction in the DKL stimulus space, a strictly physiological conception of “hue” from the hues categorizing color perception, as explained above. In drawing this distinction, and in agreement with Wool et al. (2015) and Li et al. (2022), we give no special status to V1 cells tuned to the DKL directions associated with the unique hues of color-opponent theory, while simultaneously emphasizing the cone-opponent nature of feedforward afferents from the LGN.
We express the input as the average effect of the net synaptic current on the membrane potential of a cell, following the conventions of Ermentrout (1998) and Carandini and Ringach (1997). Thus, h(θ, t) has units of mV and can take on both positive and negative values, chosen here so that a(θ, t) typically ranges from 0 to 60 spikes/sec, consistent with electrophysiological experiments penetrating individual color-responsive cells (Conway, 2001; Johnson et al., 2001; Landisman & Ts’o, 2002a; Wachtler et al., 2003).
Feedforward Input
To parameterize the input, prior work has relied on the direct relationship between cortical feature preferences and properties of the visual stimulus (Ben-Yishai et al., 1995; Bressloff & Cowan, 2003b). Cells in the cortex labeled, for instance, by their spatial frequency preferences can be mapped directly onto a visual space parameterized by the same variable. Thus, the activity of each neuronal population is no longer labeled purely by its position on the cortical sheet, but also by its preferred stimulus in an analogous feature space.
The corresponding network topology may be modeled on the cortical histology, such as the orientation map of Bosking, Zhang, Schofield, and Fitzpatrick (1997) or spatial frequency maps addressed in Bressloff and Cowan (2002), Bressloff and Cowan (2003a), and Bressloff and Cowan (2003b). Conversely, it may be based entirely on functional considerations, as, for instance, in the orientation tuning model of Sompolinksy et al. (Ben-Yishai et al., 1995; Hansel & Sompolinksy, 1998), also known as the “ring model,” which posits a topology based on the experimentally motivated assumption that populations with similar orientation preferences are maximally connected (Ben-Yishai et al., 1995) and on the argument that the important features of such a connectivity are captured by its first-order Fourier components (Hansel & Sompolinksy, 1998).
Recurrent Interactions
Next, we map the weighting function onto hue space, drawing from the hue tuning micro-architecture revealed by the imaging studies of Liu et al. (2020) and Xiao et al. (2007). These studies point to a linear relationship between distance and hue angle difference, which minimizes the wiring length of cells tuned to similar hues (Liu et al., 2020). The hue-preferring cells inhabit the so-called “color regions,” defined as such for their activation by red-green grating stimuli (Liu et al., 2020). These regions predominantly overlap with the V1 CO blobs (Landisman & Ts’o, 2002b; Li et al., 2022) and are responsive to the full range of hues, much like the patchy distribution of orientation maps within the V1 hypercolumns. Thus, in a similar manner to the local feature processing models of Bressloff and Cowan (2003b) and Ben-Yishai et al. (1995), we model the CO blob as a single color-processing unit consisting of N neurons labeled by the continuous hue preference variable θ ∈ [−π, π] (Bressloff & Cowan, 2003b).
Figure 2 shows the distribution of hue-responsive neurons within a typical color region (Figure 2A) as well as a more coarse-grained demarcation of peak activity within several of these regions (Figure 2B). To describe the spatial organization of their hue preference data, Xiao et al. (2007) and Liu et al. (2020) applied a linear fit to the cortical distance between two cell populations as a function of the difference in their preferred hue stimuli Δθ ≡ |θ − θ′| apart in DKL space. Note, this implies a discontinuity between θ = 0 and θ = 2π, allowing for the 2π periodicity of the hue preference label. Liu et al. (2020) report that the linear fit was able to capture the micro-organization of 42% of their tested hue maps, and a regression performed by Xiao et al. (2007) on an individual hue map gave a squared correlation coefficient of R2 = 0.91.
(A) Hue map of individual hue-selective cells obtained by 2-photon calcium imaging of neuronal responsiveness to seven test hues. Scale bar: 200 μm. (B) Regions of peak response to test hues (solid contours). The dashed white lines demarcate the color-preferring regions, colocalized with the CO blobs. Scale bar: 0.5 mm. Modified with permission from Liu et al. (2020).
(A) Hue map of individual hue-selective cells obtained by 2-photon calcium imaging of neuronal responsiveness to seven test hues. Scale bar: 200 μm. (B) Regions of peak response to test hues (solid contours). The dashed white lines demarcate the color-preferring regions, colocalized with the CO blobs. Scale bar: 0.5 mm. Modified with permission from Liu et al. (2020).
Cortical connectivity functions. (A) A difference of two Gaussians, one characterizing the excitatory interactions (here with σE = 40°) and the other the inhibitory interactions (σI = 90°). This is the connectivity typically assumed in mean field models of cortical processing. (B) The difference of cosines formulation (Equation 9), with J0 = −2 and J1 = 3, captures the local excitation and lateral inhibition assumed in panel A.
Cortical connectivity functions. (A) A difference of two Gaussians, one characterizing the excitatory interactions (here with σE = 40°) and the other the inhibitory interactions (σI = 90°). This is the connectivity typically assumed in mean field models of cortical processing. (B) The difference of cosines formulation (Equation 9), with J0 = −2 and J1 = 3, captures the local excitation and lateral inhibition assumed in panel A.
Here, it is also important to distinguish between the connectivity function and the center-surround receptive fields of single- and double-opponent color cells (Shapley & Hawken, 2011). While the structures of both can be approximated by the same functional form, the resemblance is superficial: the former characterizes the interactions between different neuronal populations, and the latter is a property of single cells, often adapted for computer vision algorithms (Somers, Nelson, & Sur, 1995; Turner, Schwartz, & Rieke, 2018).
RESULTS
Evolution of Network Activity
The limits of integration δ1 and δ2 in Equation 19 are the angles corresponding to the intersection of α (in gray) and cos(θ + γ) (in black). Here, c = 1, β = 1, and T = −1. = π/8. J0 and J1 are as in Figure 3.
The limits of integration δ1 and δ2 in Equation 19 are the angles corresponding to the intersection of α (in gray) and cos(θ + γ) (in black). Here, c = 1, β = 1, and T = −1. = π/8. J0 and J1 are as in Figure 3.
Separating Equation 25 into its real and imaginary parts, and noting that a real-valued activity profile a(θ, t) requires c0 ∈ ℝ and c1 = −, reduces the system to a set of equations for c0(t), (t), and (t).
Written in this form, the system provides a representation of the time evolution of a(θ, t) in terms of the coupled evolution of the constants c|ν|≤1. It is important to note that these equations are nonlinear due to the implicit Heaviside in our determination of δ1(t) and δ2(t). While our reformulation of the right-hand side of Equation 24 allows for the explicit representation of the coupling of cν via the nonlinearity, it is also this coupling that proves the analytical solution of the trajectories intractable. Thus, to describe the behavior of the time-dependent solution, we turn next to a numerical analysis of the system’s phase portrait—that is, to an exploration of the features and stability of the system’s emergent steady states.
Steady-State Solution
We approach the solution to Equation 1 with a Forward Euler method, propagating the activity from a random array of spontaneous initial values between 0 and 0.2 spikes/sec to its steady-state value. Within each time step (typically chosen to be 1 msec), we coarse-grain the network into n = 501 populations with hue preferences separated evenly across the DKL angle domain [−π, π]. The choice of an odd n allows us to numerically integrate Equation 10 using the Composite Simpson’s Rule, whereupon we rectify {h(θ, t) − T} and evaluate the right-hand side of Equation 1. Below, we use the term tuning curve only in reference to the emergent steady-state activity profiles.
Steady-state activity profile for a neuronal network encoding stimulus = π/8. Parameters are as in Figure 4.
Steady-state activity profile for a neuronal network encoding stimulus = π/8. Parameters are as in Figure 4.
We emphasize that the values of the cortical parameters J0, J1, c, T, and β are bounded by the physiological properties of V1. Varying these parameters in the subsequent analysis is therefore an investigation of their relative effects on hue processing, and we are not fine-tuning their weights to obtain specific hue tuning curves.
Here, we explore a range of values for the cortical and stimulus parameters under the constraint that the network activity remains between 0 and 60 spikes/sec, as motivated above. We further restrict J1 > 0 and J0 < 0 to elicit the local excitation and global inhibition connectivity ansatz of previous neural field models. Our main aim is to graphically characterize the relative effects of the parameters on the width, Δc = − , and peak height, a∞(), of the network tuning curves. Together, these two properties reflect the network selectivity and emergent signal strength, respectively. Note that these effects are robust to small additive white noise and may also be gleaned from the net input, expressed as in Equation 20 and evaluated at the steady-state values of the coefficients.
It is also important to note here the difference between a network tuning curve and a single-neuron tuning curve. The former is a coarse-grained representation of the CO blob response, with the horizontal axis representing the gamut of hue preferences within a single network. A relatively large tuning width would therefore indicate considerable responses from a wide range of hue tuning cells and poor network selectivity. The single-neuron tuning curve, on the other hand, is an electrophysiological recording of an individual cell’s response to a set of hue stimuli, with the horizontal axis representing the range of stimulus hue angles used in the experiment. The peak location of the single-neuron tuning curve would therefore indicate the hue preference of the individual neuron, while the width would represent its selectivity for that specific hue. Thus, though the two types of tuning curves are labeled and shaped similarly, the latter is only useful to characterize our network’s constituent neurons and not the emergent properties of the population as a whole (Bressloff & Cowan, 2003b).
Roles of the stimulus strength and cortical threshold.
We begin by considering the role of the stimulus signal strength c on the hue tuning width and peak height. Figure 6 shows typical tuning curves for two values of J1. We find that the stimulus strength has a quickly saturating effect on Δc for all J1 > 0, which is more pronounced at lower values of c as J1 → 0. Above saturation, the main contribution of the chromatic signal is to increase the network response, that is, to increase a∞().
Effect of c on the tuning curve properties. The tuning role of c quickly saturates, while its effect on the network response rate grows without bound. For = 0, β = 1, T = −1, and J0 = −1. (A) J1 = 0.2. (B) J1 = 0.7.
Effect of c on the tuning curve properties. The tuning role of c quickly saturates, while its effect on the network response rate grows without bound. For = 0, β = 1, T = −1, and J0 = −1. (A) J1 = 0.2. (B) J1 = 0.7.
We also note that at T = 0, the trend reverses, such that increasing c has no effect on the tuning width at T = 0 and a widening effect for T > 0. Figure 7 illustrates this reversal with four tuning curves of matched parameters and varying values of T. The coupling of c and T must be considered because some neural field models (see Amari, 1977; Carandini & Ringach, 1997; Dayan et al., 2001) take T = 0 for mathematical simplicity. Indeed, we might expect that there is no more physiological significance to choosing a threshold potential of T = 0 mV than any other value, beyond their relative magnitudes to h(θ, t). However, the independence of c and Δc at T = 0 and the significance of the relative signs of c and T elsewhere suggest quite the opposite. The effect of the chromatic input on tuning the network hue selectivity weakens not only once the anisotropic strength parameter, J1, is large enough to predominate, but also as T → 0.
Effect of c on the tuning curve for varying values of T with β = 1, J0 = −1, J1 = 0.2, and = 0. Note that the small network response rates are due to the low values of c chosen here. (A) T = −5. (B) T = −1. (C) T = 0. (D) T = 0.5.
Effect of c on the tuning curve for varying values of T with β = 1, J0 = −1, J1 = 0.2, and = 0. Note that the small network response rates are due to the low values of c chosen here. (A) T = −5. (B) T = −1. (C) T = 0. (D) T = 0.5.
The coupling of c and T is equally significant to the effects of T on the tuning curve properties. Figure 8 shows that below a certain value, T primarily modulates a∞(). However, for comparable magnitudes of the stimulus strength and threshold, |c| ∼ |T|, we see a transition in which T also begins to sharpen the tuning curve and continues to do so until the threshold surpasses h(θ, t) for all θ (i.e., for = = 0). Accordingly, for higher stimulus strengths, the thresholding nonlinearity plays a greater role in modulating the network selectivity at lower and a wider range of T values.
Effect of T on the tuning curve properties. = 0, β = 1, J0 = −3, J1 = 2, and c = 10.
Effect of T on the tuning curve properties. = 0, β = 1, J0 = −3, J1 = 2, and c = 10.
Roles of the cortical weights.
The anisotropic connectivity strength J1 exhibits similar relationships to the tuning curve properties to those of c. That is, for T < 0, a∞() grows and Δc narrows with increasing J1 (see Figure 9A). The trend with respect to Δc reverses for T > 0 (Figure 9B), whereas the trend with respect to a∞() remains unaffected.
Effect of J1 on the tuning curve properties for varying values of T. β = 1, c = 0.3, J0 = −10, and = 0. (A) T = −1. (B) T = 0.2.
Effect of J1 on the tuning curve properties for varying values of T. β = 1, c = 0.3, J0 = −10, and = 0. (A) T = −1. (B) T = 0.2.
These similarities are a mark of the competition between the external input and the cortical parameters in driving the network selectivity and reflect the fact that both parameters modulate the anisotropic terms of the model. This means that the role of J1 in driving network selectivity becomes more significant with decreasing stimulus strength (see Figure 10). However, a large external input does not suppress the contribution of J1 to the overall network activity. That is, increasing J1 results in raising a∞(), regardless of the strength of the stimulus. Similarly, a relatively large value of J1 does not restrict the growth of the network response with increasing stimulus strength. Thus, the anisotropic tuning introduced by the external input and the recurrent interactions act cooperatively to raise the network’s response to the stimulus hue, and competitively to tune its selectivity.
Effect of J1 on the tuning curve properties for different stimulus strengths. β = 1, T = −5, J0 = −9, and = 0. (A) c = 1. (B) c = 3.
Effect of J1 on the tuning curve properties for different stimulus strengths. β = 1, T = −5, J0 = −9, and = 0. (A) c = 1. (B) c = 3.
In contrast, J0 acts cooperatively with the external stimulus to sharpen the curves. As shown in Figure 11, the tuning curves narrow with decreasing values of J0, that is, with an increase in the relative strength of global inhibition to global excitation, a trend which is conserved for various stimulus strengths. Furthermore, there is no trend reversal at T = 0. Rather, for much of the parameter space, J0 acts with the thresholding to sharpen the tuning curves, as is illustrated in Figure 12. This could be expected from the fact that at each point throughout the dynamics, both T and J0 act isotropically on all hue preferences, lowering or raising the input for all contributing neurons. However, this commonality also means that for |T| >> |c| (where the effect of T on Δc saturates, as explained above), the thresholding suppresses the role of J0, analogous to the competition between c and J1. Finally, Figures 11 and 12 also show that increasing the global inhibition acts to reduce the value of a∞() for all c and T.
Effect of J0 on the tuning curve properties for varying stimulus strengths. β = 1, T = −2, J1 = 2, and = 0. (A) c = 1. (B) c = 6.
Effect of J0 on the tuning curve properties for varying stimulus strengths. β = 1, T = −2, J1 = 2, and = 0. (A) c = 1. (B) c = 6.
Effect of J0 on the tuning curve properties for varying values of T. β = 1, c = 1, J1 = 2, and = 0. (A) T = −5. (B) T = 0.2.
Effect of J0 on the tuning curve properties for varying values of T. β = 1, c = 1, J1 = 2, and = 0. (A) T = −5. (B) T = 0.2.
We thus conclude that the emergent hue curves in V1 are both inherited from the LGN and built on the recurrent interactions. The competition between J1 and c points to a continuum of regimes in which either hlgn or hctx dominates. However, in all regimes, J0 works cooperatively with c to narrow the curves, and all the parameters work together to raise the network response. Likewise, the competition between J0 and T (both cortical parameters) is modulated by the value of c, and the location of the peak is always completely determined by the LGN signal, regardless of the relative magnitudes of the cortical and stimulus strength parameters (see Methods: Evolution of Peak Angle).
Comparisons with the orientation tuning ring model.
Finally, we seek to compare the emergent properties of the hue tuning model with those of the orientation tuning ring model (Ben-Yishai et al., 1995; Hansel & Sompolinksy, 1998). This leads us to separate the analysis into two regions: one corresponding to the analytical regime with J0 < and J1 < , and the other to the extended regime with J1 ≥ and J0 constrained as described in the section Stability Analysis. As detailed in Methods: Linear Solution, the former regime defines the (J0, J1) parameter space wherein the model permits a closed-from stable solution for cases in which the input into all cells is above threshold. By contrast, the dynamics in the latter, extended regime always implement thresholding and thus do not permit the linear closed-form solution. For comparison purposes, note that these parameter regimes are analogous to the orientation model’s homogeneous and marginal regimes, respectively, labels which refer to the system’s responses to unoriented stimuli.
The differing assumptions underlying the formulation of h(θ, t) have important implications for the subsequent parameter analyses adopted by our two models. In the orientation tuning model, the authors detail the pronounced shift in the relative roles of the cortical and stimulus parameters in narrowing the tuning curve. In this setup, for ϵ → 0.5, an increase in c widens the tuning curve, whereas for ϵ → 0, the tuning curve selectivity is completely determined by the cortical parameters. The latter scenario constrains the value of the analogous anisotropic cortical parameter, J2, to the marginal regime.
In contrast, our model does not apportion separate regions of the parameter space to external and recurrent mechanisms. Rather, in both the analytical and extended regimes, the roles of c and J1 exist on a spectrum, where the effect of each parameter is suppressed by larger values of the other. Of course, this suppression is more stark in the extended regime because it covers larger values of J1. In this sense, our color model draws a similar conclusion to that of the orientation model: when the anisotropic tuning provided by the recurrent interactions is large, the tuning from the stimulus is negligible, and vice versa. However, we emphasize that the transition is not sharp and that c does have an effect on the tuning curve selectivity in the extended regime (see Figure 6B), as does J1 in the analytical regime.
In this regard, the two models are more consistent in their interpretations of J0’s contribution to the selectivity of the tuning curves. That is, in the two regimes of each model, the inhibition acts cooperatively with the thresholding to sharpen the tuning curves. Here again, the orientation model makes a distinction between the marginal phase (i.e., ϵ = 0 and J2 ∈ marginal regime), wherein the tuning curve width is completely determined by the cortical anisotropy, and all other cases, where the isotropic inhibition and stimulus come into play. For these cases, the authors argue, J0 does not act alone to narrow the curve: though J0 may sharpen the tuning curves, it is the anisotropy from the input or cortical interactions which acts as the source of the orientation selectivity. Although our color model’s tuning mechanism, too, requires a source of anisotropy, we have emphasized above that there is no single source of hue selectivity. When J1 is small, in either regime, both the stimulus and the uniform inhibition are significant to the hue tuning mechanism.
Ultimately, the orientation model sets up a dichotomy between two specific regions of parameter space. In the nonmarginal case, c is the primary player in the tuning mechanism, and in the marginal case, this role belongs to J2. The uniform inhibition is thus given a secondary “sharpening” role. By contrast, in choosing a fully anisotropic hlgn, our color model does not encompass an analogous marginal phase with an always dominating J1. Rather, even for large J1, the uniform inhibition is at least equally important to the modulation of the tuning width. In fact, as we have shown above, for larger values of c, J0 is more effective than J1 in narrowing the tuning curves, for both the analytical and extended regimes.
We thus stress that the two regimes, though analogous to those of the orientation model, do not constitute a division in the hue processing mechanism. Rather, we define the boundary between the analytical and extended regimes solely by whether or not the linear case exists. It is therefore determined by the values of J0 and J1 for which the linear solution applies, given that the values of c, T, and β keep the input above threshold throughout the dynamics (Methods: Linear Solution). We note that for each combination of J0 and J1 within the analytical regime there exists also a nonlinear case, in which h(θ, t) is cut off by the thresholding nonlinearity and, thereby, the linear solution does not apply. Our definition differs from that of the orientation model, which demarcates the boundary between the homogeneous and marginal phases based on the emergent steady-state tuning curves alone. For more on this approach, see the discussion of the broad and narrow profiles in Hansel and Sompolinksy (1998). As we will show next, the boundary is integral to the corresponding stability analysis of the steady-state tuning curves.
Stability Analysis
To analyze the stability of the emergent tuning curves, we turn once more to our separable activity ansatz assumed in the eigenfunction decomposition of Equation 15. This means that we are faced again with a nonlinearity induced coupling of the time-dependent coefficients and, consequently, the analytical intractability of the associated stability analysis. We therefore set up the Jacobian matrix for a numerical analysis of the local stability.
We note that the existence of a steady state is a function of the cortical strengths J0 and J1. By fixing the values of β, c, , and T, we are left with a two-parameter family of differential equations, allowing us to analyze this dependence numerically in the associated (J0, J1) parameter space.
Carrying out a parameter sweep across this space, we find that the system features a bifurcation curve, below which the model permits steady-state solutions and above which no equilibrium exists. To determine stability within the former region, we compute 𝕁 at the emergent steady-state tuning curves of various points in the parameter space. Solving the associated characteristic equations, we observe that the eigenvalues are always real and negative, and thus conclude that all emergent steady-state tuning curves are stable.
Figure 13 shows the bifurcation diagrams for two families of equations, distinguished by their values for T. Most notable is the extended regime, which permits stable steady-state solutions beyond the boundary set by the linear case (Methods: Linear Solution). As this parameter regime is not accessible to the linear solution, the tuning curves in this regime are necessarily a product of the thresholding nonlinearity and are thus always cut off below |θ| = π. The thresholding nonlinearity therefore not only expands the region of stability, but also ensures that the tuning curves emerging within the extended regime are selective for hue. As we have seen, this expansion is pivotal when the external input is weak and the anisotropic cortical strength plays the larger role in narrowing the tuning curves. Furthermore, regardless of input strength, it allows for a larger overall network response, as the peak activity, a∞(), grows with increasing J1. Finally, as we will see in the following section, in the absence of any stimulus (i.e., for c = 0), the extended regime features the spontaneous generation of stable tuning curves and may thus serve as the bedrock for color hallucinations.
Bifurcation diagram for β = 1 and c = 1 for two values of threshold (shown in the legend). The gray and white regions correspond to the analytical and extended regimes, respectively. The black dashed line is the bifurcation curve, above which the tuning curves grow without bound. The overlaid symbols correspond to points tested in a parameter sweep over the extended regime. Notably, the parameter sweep produces the same bifurcation curve for both values of T. Here, we must note that for critical values of T, for which the input is not large enough to generate activity, the model permits the trivial a(θ) = 0 steady-state solution in both the analytical and extended regimes. This solution, however, is unstable to perturbations large enough to make the input cross the threshold. For more on bifurcation theory in the context of neural fields, see Gross (2021) and Hesse and Gross (2014). See also Hansel and Sompolinksy (1998) for an analogous “phase diagram” analysis of the orientation ring model.
Bifurcation diagram for β = 1 and c = 1 for two values of threshold (shown in the legend). The gray and white regions correspond to the analytical and extended regimes, respectively. The black dashed line is the bifurcation curve, above which the tuning curves grow without bound. The overlaid symbols correspond to points tested in a parameter sweep over the extended regime. Notably, the parameter sweep produces the same bifurcation curve for both values of T. Here, we must note that for critical values of T, for which the input is not large enough to generate activity, the model permits the trivial a(θ) = 0 steady-state solution in both the analytical and extended regimes. This solution, however, is unstable to perturbations large enough to make the input cross the threshold. For more on bifurcation theory in the context of neural fields, see Gross (2021) and Hesse and Gross (2014). See also Hansel and Sompolinksy (1998) for an analogous “phase diagram” analysis of the orientation ring model.
However, looking back at Figure 13, perhaps most striking is the horizontal portion of the bifurcation curve at J0 = for J1 < , which sets the same stability conditions on J0 and J1 as in the linear case. This is despite the fact that many of the points in the analytical regimes of the two featured families correspond to solutions that implement thresholding, thus signifying that the analytical regime is not an exclusively linear one.
The key to understanding the shape of this region lies in noticing that the bifurcation diagram does not change for varying values of c, T, and , as shown in Figure 13 for the two values of T. The stability conditions on J0 and J1 are thus uniquely determined by β alone. Furthermore, for the general diagram (i.e., with β fixed and c, T unfixed), each point of the analytical regime permits linear solutions, in addition to the ones that implement thresholding. Accordingly, the uniqueness of the bifurcation diagram implies that at each point of the analytical (J0, J1) subspace, the stability of the latter, nonlinear solutions is equivalent to that of the linear solutions. This means that the boundary at J0 = set by the linear case (Methods: Linear Solution) applies to the full, nonlinear model as well.
A Turing Mechanism for Color Hallucinations
Biological Turing patterns.
Underpinning our hue tuning model is the mathematics of reaction-diffusion systems, for which, in particular, Alan Turing’s treatment of biological pattern formation offers many valuable insights (Turing, 1952). The general Turing mechanism assumes a system of two interacting chemicals, whose local reaction and long-range diffusion properties govern the dynamics of their relative concentrations. In the original framework these chemicals are termed “morphogens” to elicit their form-producing capabilities within a developing embryo, whose anatomical structure emerges as a result of their underlying concentration dynamics. This, for instance, may be attributed to the morphogens’ catalysis of organ formation in different parts of the developing organism.
In the absence of diffusion (i.e., with DX = DY = 0), the system has a homogeneous steady-state solution, (X, Y) = 0, whose stability is determined by a Jacobian composed of the reaction rates, , and hence by the system’s local chemical properties alone (Hoyle, 2006). Note that at this point the system is circularly symmetric with respect to interchanging any two cells on the ring.
From here, a set of additional conditions may be placed on DX and DY so that the system undergoes a Turing bifurcation, wherein at least one λμ becomes positive and the homogeneous steady state loses its stability. With the addition of a small random perturbation, the instability results in the growth of the corresponding eigenmodes eiμθ, such that, over time, Equation 35 is dominated by the eigenmodes with largest λμ. These represent stationary waves whose wavelengths are set by the circumference of the ring (i.e., by the spatial properties of the medium) and whose growth is bounded by the higher order terms that had been initially ignored in the near-equilibrium formulation (Murray, 2003; Scholz, 2009; Turing, 1952). The underlying circular symmetry is thus broken and a spatial pattern is formed.
In his seminal paper, Turing extrapolated this mechanism to explain various biological phenomena, such as the development of petals on a flower, spotted color patterns, and the growth of an embryo along various directions from an original spherical state. A hallmark of each of these examples is that there is no input into the system, so the emergent patterns reflect a mechanism of spontaneous symmetry breaking, onset by a perturbation of “some influences unspecified” (Turing, 1952). In light of this, we ask, can the visual cortex self-generate the perception of hue?
Spontaneous symmetry breaking and color hallucinations.
To assess our model’s ability to self-organize in the absence of visual input, we set c = 0 and seek to establish the presence of a Turing mechanism marked by the following three features:
A system comprised of local excitation and long-range inhibition.
Spontaneous symmetry breaking in the absence of input within a region of a parameter space defined by the relevant bifurcation parameter(s).
The emergence of patterns that are bounded by the system’s nonlinearities.
The onset of color hallucinations in the (J1, T) parameter space. The model generates spontaneous hue tuning curves beyond J1 = ≡ and below T = 0.
The onset of color hallucinations in the (J1, T) parameter space. The model generates spontaneous hue tuning curves beyond J1 = ≡ and below T = 0.
Spontaneous pattern formation in the absence of input (c = 0). β = 1. (A–B): J0 = −2, J1 = 0.1 (A) T = 0 (B) For T < 0, the homogeneous steady-state value increases. Here, T = −10. (C–D): Pattern formation in the extended regime for J0 = −2, J1 = 0.4. (C) No hue tuning curve emerges for T ≥ 0. Here, T = 0. (D) T = −10. A hue tuning curve is generated in the absence of external input. (E–F): T = −10, J0 = −7, J1 = 6. The emergent tuning curve is more selective for larger values of J1. For each run, the activity is peaked about a different angle, set by the random initial conditions. The peak value and tuning width are consistent between trials.
Spontaneous pattern formation in the absence of input (c = 0). β = 1. (A–B): J0 = −2, J1 = 0.1 (A) T = 0 (B) For T < 0, the homogeneous steady-state value increases. Here, T = −10. (C–D): Pattern formation in the extended regime for J0 = −2, J1 = 0.4. (C) No hue tuning curve emerges for T ≥ 0. Here, T = 0. (D) T = −10. A hue tuning curve is generated in the absence of external input. (E–F): T = −10, J0 = −7, J1 = 6. The emergent tuning curve is more selective for larger values of J1. For each run, the activity is peaked about a different angle, set by the random initial conditions. The peak value and tuning width are consistent between trials.
We further observe that beyond J1 = , a stable homogeneous steady-state solution remains at a∞(θ) = 0 for T ≥ 0 (Figure 15C). However, at T = 0, this radial symmetry is broken, and the cortex generates spontaneous tuning curves with peak locations determined by the random initial conditions (Figure 15D–F). Thus, the system bifurcates when J1 = and T = 0, permitting the onset of color hallucinations in a region defined by these two values (Figure 14). Note that the unimodal tuning curves predict stationary, single-hued phosphenes. Extension to other CO blob networks would therefore indicate a hallucination comprised of multiple phosphenes of varied hues, each determined by the local cortical activity at hallucinatory onset.
Bearing these predictions in mind, we point to a recent functional MRI study of blind patients experiencing visual hallucinations (Hahamy, Wilf, Rosin, Behrmann, & Malach, 2021). The study attributes these visions to the activation of the neural networks underlying normal vision, precipitated by the hyper-excitability of the cortex to spontaneous resting-state activity fluctuations when it is deprived of external input. This is suggestive of the required lowering of neuronal threshold at the onset of color hallucinations predicted here. Notably, a reduction in membrane potential threshold has also been attributed to the action of hallucinogens (De Vos, Mason, & Kuypers, 2021; Varley, Carhart-Harris, Roseman, Menon, & Stamatakis, 2020).
Finally, we note that the stability of the emergent tuning curves is determined by the bifurcation diagram of Figure 13. This means that, in addition to expanding the region of stability in the presence of chromatic stimuli, the model’s nonlinearity allows for stable, spontaneous color hallucinations in their absence.
Having thus established a Turing-like mechanism for our model’s self-organization, we end with an analogy to Turing’s original diffusion-driven formulation. In his concluding example, Turing applies the mechanism to explain the growth of an embryo along various axes of its original spherical state. This growth is driven by diffusion, directed by the “disturbing influences,” shaped by the system’s chemical and physical properties, and bounded by the system’s nonlinearities. It is all too clear to see the parallels with our hue tuning model, wherein a hallucination is driven by the anisotropy of the cortical interactions, its hue determined by the initial conditions, its selectivity shaped by the cortical parameters, and its stability ensured by the thresholding nonlinearity.
DISCUSSION
This paper presents a neural field model of color vision processing that reconceptualizes the link between chromatic stimuli and our perception of hue. It does so guided by the premise that the visual cortex initiates the mixing of the cardinal L–M and S–(L+M) pathways and thereby transforms the discrete cone-opponent signals to a continuous representation of chromatic information. Such mixing mechanisms have been implemented by previous combinatorial models of color processing, though through a largely feedforward approach or at the level of the single neuron.
Our theory bears in mind the mixing mechanism, but reframes the stage-wise combinatorial scheme to one based on the nonlinear population dynamics within the visual cortex. Accordingly, we propose a hue-based cortical connectivity, built upon the cortical hue map micro-architecture revealed by recent optical imaging studies of V1. By considering the intracortical network interactions, we have accounted for V1 cells responsive to the gamut of DKL directions without the need to fine-tune the cortical parameters. We do so without restricting to a particular category of V1 neuron, as both single-opponent and double-opponent, and altogether novel types of cells, have been suggested as the primary messengers of chromatic information. Rather, we zoom out from the individual neuron’s receptive field to model the aggregate, population-level properties and, in particular, the stable representation of hue. We thereby offer that chromatic processing in the visual cortex is, in its essence, a self-organizing system of neuronal-activity pattern formation, capable of encoding chromatic information in the presence of visual stimuli and generating information in their absence.
Further, in assuming modularity for chromatic processing, we have not ruled out a mechanism for joint feature processing. Our choice to focus on the unoriented color cells within the CO blob regions allowed us to parse out the chromatic pathway for an independent analysis, but should not be interpreted as a claim about its functional independence. We leave open the question of the functional and anatomical separation of the various visual streams.
Equally unsettled is the question of how much S cone input contributes to the mixing of the cone-opponent channels, with some studies showing a relatively weak S cone input into the neurons of V1, compared to its L and M cone counterparts (Li et al., 2022; Xiao, 2014). The variations across these experiments may stem, in part, from differences in optical imaging and electrode penetration techniques, including the particulars of the chromatic stimulus used (Li et al., 2022; Liu et al., 2020; Polimeni, Granquist-Fraser, Wood, & Schwartz, 2005; Salzmann et al., 2012). On the whole, however, single-cell recordings have identified two main types of color-responsive regions: color patches that contain neurons tuned exclusively to stimuli modulating either of the cone-opponent pathways, and patches with neurons exhibiting a mixed sensitivity to a combination of the two (Landisman & Ts’o, 2002a; Li et al., 2022; Livingstone & Hubel, 1984). Further experiments on the connectivity between these regions, and among the single- and double-opponent color cell populations of which they consist, may point to added micro-architectures for the hue maps, along the lines of the geometric orientation models of Bressloff and Cowan (2003a) and Bressloff et al. (2001).
Finally, we emphasize that the mechanism we offer departs from previous combinatorial color models that predict hue sensation at the final stage of processing (De Valois & De Valois, 1993; Mehrani et al., 2020), as well as neural field models that conflate cone- and color-opponency in their interpretations (Faugeras, Song, & Veltz, 2022; Smirnova, Chizhkova, & Chizhov, 2015; Song, Faugeras, & Veltz, 2019). The emergent hue tuning curves we have characterized are a network property reflective of the physiological neuronal responses, and should not be confounded with our perception of hue. A photon of wavelength 700 nm striking a retina is no more “red” than any other particle—color is a perceptual phenomenon not yet represented in these first stages of vision. By recognizing that the hue tuning mechanism of the visual cortex is an early stop in the visual pathway, we point to the need for further field theory approaches to our understanding of color perception.
METHODS
Linear Solution
Evolution of Peak Angle
Linear Stability Analysis
This section presents the mathematical details for obtaining Equation 30.
ACKNOWLEDGMENTS
The authors acknowledge the fruitful and stimulating discussions with Wim van Drongelen and Graham Smith, and research support from the Oberlin College libraries.
AUTHOR CONTRIBUTIONS
Zily Burstein: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Project administration; Resources; Software; Supervision; Validation; Visualization; Writing – original draft; Writing – review & editing. David D. Reid: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Project administration; Software; Supervision; Validation; Visualization; Writing – review & editing. Peter J. Thomas: Conceptualization; Methodology; Resources; Writing – review & editing. Jack D. Cowan: Conceptualization; Methodology; Project administration; Resources; Writing – review & editing.
FUNDING INFORMATION
Peter J. Thomas, Directorate for Mathematical and Physical Sciences (https://dx.doi.org/10.13039/100000086), Award ID: DMS-2052109. Peter J. Thomas, Office of Research Infrastructure Programs, National Institutes of Health (https://dx.doi.org/10.13039/100016958), Award ID: R01 NS118606.
TECHNICAL TERMS
- Cone-opponency:
Receptive field property of neurons in the early visual pathway, by which chromatic stimuli are processed through the comparison of the relative L, M, and S cone responses.
- Stimulus space:
A geometrical construct in which chromatic stimuli are represented by the relative cone responses they yield.
- Color-opponent theory:
Theory of color appearance that postulates that the four unique hues—red, green, blue, and yellow—are perceived antagonistically. That is, there is no such thing as a reddish green or a bluish yellow.
- Cardinal pathways:
In DKL space, the two orthogonal axes representing stimuli isolating the L–M and S–(L+M) cone-opponent pathways.
- Wilson-Cowan neural field equations:
Coupled set of partial differential equations describing the network dynamics of excitatory and inhibitory neural populations.
- Spontaneous pattern formation:
A system’s ability to self-generate new symmetries in the absence of external input.
- Activation function:
A function mapping the afferent input into a population of neurons (expressed as a current or membrane potential) to the population’s firing rate or probability of firing.
- Bifurcation curve:
A curve in parameter space defining a transition in the system’s dynamics and stability.
- Reaction-diffusion system:
Mathematical model of two or more interacting substances consisting of local dynamics (via the reaction terms) and spatial dynamics (via the diffusion terms).
- Globally asymptotic stable steady state:
A system’s equilibrium solution which attracts all other solutions regardless of initial conditions.
REFERENCES
Author notes
Competing Interests: The authors have declared that no competing interests exist.
Handling Editor: Gustavo Deco