Abstract
A new computational model of the primary visual cortex (V1) of the macaque monkey was constructed to reconcile the visual functions of V1 with anatomical data on its LGN input, the extreme sparseness of which presented serious challenges to theoretically sound explanations of cortical function. We demonstrate that, even with such sparse input, it is possible to produce robust orientation selectivity, as well as continuity in the orientation map. We went beyond that to find plausible dynamic regimes of our new model that emulate simultaneously experimental data for a wide range of V1 phenomena, beginning with orientation selectivity but also including diversity in neuronal responses, bimodal distributions of the modulation ratio (the simple/complex classification), and dynamic signatures, such as gamma-band oscillations. Intracortical interactions play a major role in all aspects of the visual functions of the model.
SIGNIFICANCE STATEMENT We present the first realistic model that has captured the sparseness of magnocellular LGN inputs to the macaque primary visual cortex and successfully derived orientation selectivity from them. Three implications are (1) even in input layers to the visual cortex, the system is less feedforward and more dominated by intracortical signals than previously thought, (2) interactions among cortical neurons in local populations produce dynamics not explained by single neurons, and (3) such dynamics are important for function. Our model also shows that a comprehensive picture is necessary to explain function, because different visual properties are related. This study points to the need for paradigm shifts in neuroscience modeling: greater emphasis on population dynamics and, where possible, a move toward data-driven, comprehensive models.
Introduction
Orientation selectivity (OS) is one of the most salient characteristics of the primary visual cortex (V1) of the macaque monkey. Hubel and Wiesel (1962) proposed a theory for OS framed in terms of feedforward convergence of LGN cells with receptive fields (RFs) aligned in visual space. We were compelled to revisit their ideas in the light of experimental findings about the extreme sparseness of LGN input to V1.
We focus on V1 layer 4Cα because orientation selectivity and spatial frequency (SF) selectivity are already observable in this cortical input layer (Livingstone and Hubel, 1984; Ringach et al., 2002; Zhu et al., 2010). Layer 4Cα receives input from the magnocellular LGN layers (Chatterjee and Callaway, 2003). According to anatomical evidence (see Results), the magnocellular LGN input to V1 is very sparse (Connolly and Van Essen, 1984; Silveira and Perry, 1991; Angelucci and Sainsbury, 2006). Layer 4Cα in each V1 hypercolumn (Lund et al., 2003) receives only ∼10 direct magnocellular inputs, and a total of 30 LGN inputs altogether counting up direct and ramifying axons from neighboring hypercolumns. These 30 LGN inputs are driving ∼3000 excitatory (E) and 1000 inhibitory (I) cortical neurons.
The anatomical data cited above lead to the following geometric/combinatorial questions. In a patch of 10–20 LGN cells, one-half ON, one-half OFF, is it possible to fit two or three approximately parallel receptive field subregions, as in the Hubel–Wiesel picture, each driven by ON (or OFF) LGN cells, with distances between the subregions that would support the correct spatial frequency preference of the V1 neuron? And, would model cortical neurons acquire the ability to discern many orientations around the clock, as observed in real V1 (Obermayer and Blasdel, 1993), from the same set of sparse inputs?
The new computational model of layer 4Cα that we have constructed answers both of these questions in the affirmative. Orientation selectivity is passed to cortex, with a great deal of help from recurrent cortical interactions. As shown in Results, our model exhibits the full range of orientation and spatial frequency selectivity found in V1. It takes weakly tuned inputs from sparse LGN afferents, and, via recurrent inhibitory and excitatory synaptic interaction, sharpens them and smoothes out the discrete set of orientations from LGN into a continuous map.
Our model goes far beyond accounting for OS and orientation maps. It shows that many V1 functions are linked. These include, in addition to OS, the bimodal distribution of modulation ratio (MR) that is an index of the simple/complex classification (see Fig. 4; Skottun et al., 1991), the wide diversity in selectivity and neuronal response magnitudes (see Figs. 6, 8; Ringach et al., 2002), as well as the emergence of gamma rhythms when driven by optimal gratings (see Fig. 9; Volgushev et al., 2002; Kayser and König, 2004).
Our analysis shows that dynamic interactions among cortical neurons are the single most important factor that binds together these different functions of V1. Recurrent excitation, indirect suppression, and the competition and balancing between the excitatory and inhibitory populations shape neuronal responses to stimuli. These interactions give rise to many emergent phenomena, i.e., phenomena that cannot be predicted from the behaviors of single neurons alone.
We propose that the model presented here can serve as a starting point in the creation of a larger dynamic picture, one that connects anatomical structures, dynamics, and cortical functions, and extends to the rest of cerebral cortex.
Materials and Methods
We will discuss in separate sections the general layout, or architecture, of our network model, the equations of neuronal dynamics, and parameter determination. For readers who wish to bypass model details, we recommend proceeding directly to Results.
Model layout
The main component of our model is that of a small piece of layer 4Cα of the macaque V1 cortex located at ∼5° eccentricity. Also modeled are the two sheets of LGN cells that project to this region of 4Cα, as well as outputs of layer 6, the feedback layer to 4Cα.
Magnocellular LGN layers.
The input to layer 4Cα in our model consists of regions of the magnocellular layers of the LGN. Data from the studies by Connolly and Van Essen (1984) and Silveira and Perry (1991) indicate that there are on average ∼10 LGN cells (5 ON and 5 OFF) that project to each hypercolumn (HC); see Results, Input patterns of LGN to V1. In the model, we started modeling the LGN with a (regular) triangular lattice on the plane, which we identified with visual space. Assuming that the centers of the receptive fields of ON LGN cells correspond to points in this lattice, we placed the OFF cells at the barycenters of the triangles defined by the ON lattice. Lattice spacing was chosen so that neighboring ON (or, respectively, OFF) cells are ∼0.125° apart, resulting in 10 cells/HC. We then perturbed randomly and independently each lattice point to match the retinal mosaics of M retinal ganglion cells.
Figure 1A shows the distribution of receptive field centers of ON and OFF LGN cells in the lattice, together with the scale of one HC in visual space. Determining plausible input patterns of LGN cells to cortex based on anatomical and physiological data is a major part of the present modeling effort. Our choices of LGN input patterns are reported in Results.
Layer 4Cα.
We modeled nine HCs of 4Cα arranged in three ocular dominance columns (Fig. 1B). Each HC measures 0.5 × 0.5 mm2 and is divided into regions where cells are intended to have similar orientation preferences, following the pinwheel pattern of orientation selectivity seen in experiments (Obermayer and Blasdel, 1993). Following Beaulieu et al. (1992), we used cell densities of ∼4000 neurons/HC, three-quarters of which are E cells and the rest are I cells. The I neurons are assumed to be a homogeneous population of local-circuit basket cells, a reasonable approximation for layer 4Cα (Defelipe et al., 1999). E and I neurons are uniformly distributed in the model cortex.
The probability of connection among model cells is dependent on distance and cell types (E or I), while the strength of connection is independent of distance. This approach is based on intracellular data (Oswald and Reyes, 2008) and is more realistic than all-to-all coupling with synaptic strength decreasing with distance, as was done in some previous large-scale models of V1 (McLaughlin et al., 2000; Tao et al., 2004, 2006). As discussed in Results (see Fig. 6), the observed diversity in OS and other cortical properties depended on the probabilistic connectivity used here. For presynaptic E neurons, connection probabilities are given by Gaussians with SD = 200/√2 μm; and for presynaptic I cells, SD = 125/√2 μm (Fitzpatrick et al., 1985; Yoshioka et al., 1994). Peak connection probability among E cells is ∼15% on average (i.e., we scaled the Gaussian describing the spatial dependence of connection probability to have a peak of 0.15). According to Holmgren et al. (2003) and Oswald and Reyes (2008), I cell connections are much denser (from 50% to 100%); we set the peak connection probabilities for E → I, I → E, and I → I to be 60% (Fig. 1C,D,E). The consequences of all these choices can be seen in Figure 1C–E, which shows examples of the spatial extent of E–E, E–I, I–E, and I–I couplings in the model. As can be observed in Figure 1, the spatial reach of E neurons is larger than that of I neurons, while the E–E coupling is the sparsest.
The numbers of connections and cell densities imply that, on average, an E cell has ∼200 E cells and ∼100 I cells presynaptic to it, while an I cell has, on average, ∼750 E cells and ∼100 I cells presynaptic to it (Fig. 1C–E). These are mean values; the actual numbers of presynaptic inputs fluctuate from cell to cell. Connectivity between each specific pair of neurons was determined first by a coin flip (with a suitable probability of “yes”). Outliers in the distributions of presynaptic neurons were then discarded to prevent, for example, domination by cells with too many presynaptic E neurons (Fig. 1C–E).
Observe that the numbers of presynaptic neurons computed from cell densities and connection probabilities are in the hundreds, not tens of thousands (as in the continuum limit of van Vreeswijk and Sompolinsky, 1996); each presynaptic cortical cell makes multiple synaptic connections on its target neuron (Martin, 1988).
Outputs from layer 6.
Layer 4Cα receives substantial feedback primarily from layer 6 (L6) (Callaway, 1998), the effects of which are incorporated in the model. L6 receives input from multiple layers of V1 as well as LGN (Callaway, 1998; Sincich and Horton, 2005), and its neuronal population is known to be very diverse (Wiser and Callaway, 1996); constructing a detailed model of its dynamics would present formidable challenges. For our purposes, however, it is sufficient to model only the outputs of the subpopulation of L6 that project to 4Cα. We modeled these outputs as a collection of spike trains, one for each L6 neuron.
From the cell density of L6 (Beaulieu et al., 1992) and the small fraction of E cells that project to layer 4Cα (Wiser and Callaway, 1996), we estimated that ∼300 L6 neurons/HC have ascending axons that terminate in 4Cα. The axonal spreads are large (Wiser and Callaway, 1996); in the model, we assumed that most (approximately five-sixths) of the synaptic contacts of a layer 6 cell occur in a disk of radius ∼180 μm, while the rest can extend as far as ∼360 μm. In the absence of experimental guidance on connection probabilities between L6 and 4Cα neurons, we assumed that they are similar to those within 4Cα; this gives, for example, an average of ∼50 presynaptic L6 neurons, represented by a corresponding number of spike trains, for each E cell in 4Cα.
Spontaneous spiking of L6 neurons was modeled as Poissonian, at 0.5–10 spikes/s (Ringach et al., 2002). L6 spike trains were tuned for orientation (Ringach et al., 2002). Specifically, if f6,max represents the mean firing rate of an L6 neuron when driven by its optimally oriented grating pattern at full contrast, we scaled this number down to approximately 0.25 * f6,max for the response to an orthogonal grating. Also incorporated in the model are the facts that (1) a majority of L6 neurons are complex (Ringach et al., 2002), and (2) there is partial synchronization in the gamma band of frequencies for nearby L6 neurons (Xing et al., 2012).
Equations of neuronal dynamics
Dynamics of LGN cells.
Following the study by Lin et al. (2012), we modeled the dynamics of individual LGN cells by the integrate-and-fire equation, as follows:
where V is membrane potential; c = 100 is the leak rate (t is in seconds); I+ and I− are deterministic currents entering ON and OFF LGN cells, respectively; Snoise is the coupling coefficient of a Poisson-timed noise term; si is ± with equal probability; and ti are the arrival times of the Poissonian noise inputs. When V reaches a threshold Vthresh ∼1, the potential is reset to 0, and a spike is sent to all postsynaptic cells in 4Cα. For a sinusoidal grating drifting in the direction u (unit vector) with temporal frequency w, the current I(t, x) into an LGN cell, the center of whose RF corresponds to location x in visual space is given by the following:
where I0 = 100 (as in the study by Lin et al., 2012), ε is contrast, C(|u|) is a function encoding the dependence of firing rate on the spatial frequency of the grating (modeled as a standard difference-of-Gaussians, as in the study by Zhu et al., 2010), and p(t,x) = wt + <x,wu> gives the phase at location x (where <x,wu> is the inner product). Constants are chosen so that the background (bg) firing rate of an LGN cell is ∼20 spikes/s, and, when driven by a grating at full contrast, the peak firing rate is ∼100 spikes/s (Hawken et al., 1996). When driven, this LGN model produces firing patterns significantly closer to those of real LGN cells (Reich et al., 1997; Lin et al., 2012) than Poisson models. The more realistic LGN model excited V1 cells more effectively than Poisson drive (Lin et al., 2012).
Dynamics of cortical cells in layer 4Cα.
Consider a model neuron, n, of type σ = E or I. To model the time evolution of its membrane potential, we used normalized voltage units where resting potential Vrest = 0, and spike threshold Vth = 1. The membrane potential of the nth neuron vn is driven toward the normalized spike firing threshold Vth = 1 by the following leaky integrate-and-fire (LIF) equation:
where vn is in normalized units, time t is in seconds, the leakage conductance gR,σ is set to 50/s for σ =E and 1.33 * 50/s for σ =I (Beierlein et al., 2003). Here gEn(t) and gIn(t) are the E and I conductances of neuron n at time t; their time evolutions are described below. Finally VE and VI are the E and I reversal potentials, which in normalized coordinates are 14/3 and −2/3 respectively. When vn reaches Vth, a spike is fired, and vn is reset to 0, where it is held for a refractory period = 2 ms. The biophysical constants above are textbook (Koch, 1999).
The evolution of gIn(t), the I conductance of neuron n, is given by the following:
where SσI is the synaptic weight, or the constant describing the change in I conductance in neuron n upon receiving synaptic input from an inhibitory neuron. Synaptic weights vary by ∼10% from neuron to neuron; for simplicity, we write only their mean, which we assume depends only on the type σ of neuron n, hence the notation SσI. The first summation in the equation for gIn(t) is over P4C,I (n), defined to be the set of all I neurons in layer 4Cα that synapse on neuron n. The second summation is over the spikes fired by neuron i. Specifically, tik is the time of the kth spike in neuron i. GGABA (s) describes the time course of I conductance for a neuron when a spike is fired by a presynaptic I neuron at time s = 0.
The E conductance gEn (t) is the sum of four synaptic conductances: (I), coming from LGN; (II), from layer 4Cα; (III), layer 6; and (IV), neuromodulatory influences from the rest of the brain or body. We discuss each of these terms separately. First:
As above, Sσ,LGN is the synaptic weight from an LGN spike for a postsynaptic cell of type σ, PLGN (n) is the set of LGN cells presynaptic to neuron n, and GAMPA (s) is the conductance time course for AMPA. Likewise:
The meanings of the terms are analogous to that given before: P4C,E (n) denotes the set of E neurons in layer 4Cα that synapse on neuron n, and ρσA and ρσN denote the fractions of synaptic input received by AMPA and NMDA receptors in a postsynaptic neuron of type σ. The feedback term (III) is identical to (II), except that P4C,E(n) is replaced by P6,E(n), the set of E neurons in layer 6 that synapse on neuron n, and SσE is replaced by S6σE, a different number. Finally, the term (IV), in which we have lumped the influence of multiple modulatory substances, to be thought of as mostly but not limited to acetylcholine, is modeled approximately (in the absence of more precise information) as follows:
Here “amb” is shorthand for “ambient.” The summation is over Poisson pulses, of size Samb and at rate rσ, occurring at times tk. The times of occurrence of ambient pulses are independent from neuron to neuron.
The functions GGABA, GAMPA, and GNMDA are similar to those used by Tao et al. (2004), with rise times of a few milliseconds; and decay times of ∼5 ms for GGABA, ∼3 ms for GAMPA, and ∼80 ms for GNMDA. Following the study by Stratford et al. (1996), we further assumed that each E-to-E spike carries a synaptic failure rate of up to 20%.
We remark that even after all the parameters have been fixed, one cannot solve for the quantities tik (= the arrival times of spikes) in the equations above explicitly. A neuron spikes when its membrane potential reaches a certain threshold that is largely determined by the net effect of the spiking activities of its presynaptic E and I neurons, which themselves are influenced by neurons presynaptic to them. In particular, firing rates cannot be expressed in terms of the parameters above. We choose parameters aiming to hit a target distribution of firing rates; that is the best we can do. As simple as the form of the LIF equations may appear, the interactions among neurons are complicated, and it is these interactions that shape collective behavior.
Parameter determination
Parameters were set systematically based as much as possible on experimental data. First, we fixed those constants for which there is some biological guidance, such as the ratios of AMPA to NMDA in E synapses. We set ρAE = 0.8, ρNE = 0.2, ρAI = 0.67, and ρNI = 0.33, as NMDA is known to be more prevalent in I neurons (Lisman et al., 2008). We set SI,LGN = 2.1 * SEE, SI,LGN = 3 * SEE, S6EE = 0.5 * SEE (Stratford et al., 1996; Beierlein et al., 2003), and assumed a synaptic failure rate up to 50% for L6 neurons (Stratford et al., 1996). As the mean firing rate of the particular subpopulation of L6 that project to 4Cα is not known, we set f6,max at ∼40 spikes/s, well within the range of the entire L6 population.
The choice of Samb = 0.01 was arbitrary, the only known constraint being that Samb should be small (acetylcholine is known to be delivered in very small pulses; Disney et al., 2012). There was less guidance for the ambient neuromodulatory rates rEA and rIA, except that the neuromodulation alone should not cause spiking in layer 4Cα (Soma et al., 2012). We used rEA, rIA ∼ 250/s. The choices of the neuromodulatory rates have little effect on orientation selectivity
Determining SEE, SEI, SIE, and SII.
We set SEE = 0.028, a justification being that starting from a reasonable membrane potential and conductance, it takes 10–20 excitatory pulses in quick succession to drive a cell across threshold if we use an SEE value between 0.02 and 0.03 (Stratford et al., 1996).
Once SEE was fixed, we used the following experimental guidance to set SEI. We may express the ratio of I conductance to E conductance as follows:
where W = (firing rate of I neurons in 4Cα)/(firing rate of E neurons in 4Cα); X = (no. of presynaptic I neurons in 4Cα)/(no. of presynaptic E neurons in 4Cα); Y = correction of W * X to include sources of E currents from outside of 4Cα; and Z = SEI/SEE.
Substituting gIE/gEE ∼ 3 based on experimental data (Borg-Graham et al., 1998), W ∼ 4 (Swadlow, 1988), X ∼ 1/2 (see Materials and Methods, Model Layout, above), and Y ∼ 2/3 (a guess based on the results of the study by Binzegger et al., 2009), we obtain Z ∼ 2.25. In the model, we used SEI = 2 * SEE.
Then we fixed the relation between SII and SEI, for example, by setting SII to be uniformly distributed between 0.65 * SEI and 0.85 * SEI. Our rationale for taking SII to be smaller than SEI is that, in addition to the usual chemical synapses, some pairs of I neurons are known to be electrically coupled (Gibson et al., 2005), and electrical coupling means that the I neurons excite rather than suppress one another.
Of the four synaptic weight constants SXY, the most sensitive is SIE; more than the other SXY weight constants, a small change in SIE can have a large effect on the behavior of the system. This is in part due to the large numbers of presynaptic E neurons to I cells (Fig. 1D). For this reason, we located a suitable value of SIE through parameter exploration. With all SXY weight constants fixed, except for SIE, we determined SIE guided by the bg firing rates for E neurons. More precisely, we found that, without exception, the mean bg E firing rate decreases monotonically as SIE/SEE is increased (Fig. 2). This dependence is expected, and it allows one to choose the value of SIE/SEE that gives the mean bg firing rates of E neurons that one requires to emulate experimental mean bg firing rate data at 3–4 spikes/s.
We do not claim to have performed an exhaustive search for all parameters (which is not feasible given the dimensionality of parameter space), but have found, through numerical exploration, that there is considerable flexibility in the sense of producing viable model outputs provided that excitation and inhibition in the model are well balanced. Below, we give three examples of how this E/I balancing is accomplished in the procedure above, as well as the properties that are affected. None of the parameter variations below has significant effects on OS, diversity, and dynamic linkage to function, as discussed in Results.
Example 1: varying SEI.
We have explored varying SEI = c * SEE for c in the range (1.6, 2.8). Using the parameter determination procedure (Fig. 2) for SIE/SEE, we find a larger SIE value for smaller c values and a smaller SIE value for larger c values, so that the total inhibitory capability of the system remains approximately constant. What changes is the I cell firing rate: a larger c value leads to lowered I cell firing, as can be expected since each I cell spike is more powerful.
Example 2: increasing SII.
If SII increases, our procedure for determining SIE (Fig. 2) gives a larger SIE value to compensate for the extra inhibition to I neurons. A side effect of increased SII is that I cell voltages sit lower on average (i.e., the membrane is more hyperpolarized), and I cell spiking becomes a little less synchronized when driven.
Example 3: varying the amount of ambient drive.
Varying ambient input spike rates rEA and rIA by similar amounts has little effect, except that too low an ambient drive causes the system to become more synchronized than is realistic in background. Favoring ambient E or I over the other is balanced by the procedure for selecting SIE, as indicated before.
Trimming variation and introducing “egalitarian” redistributions.
E neurons that receive significantly more excitatory input may fire more, silencing other E neurons via recurrent inhibition. Such competition is part of the dynamic picture, but in too extreme a form, it makes the model unrealistic. For this reason, we had to trim down the variability in the number of presynaptic E neurons for E cells, as explained in Materials and Methods, Model Layout. It turns out that we will also have to introduce variability for a different reason: to compensate neurons that do not receive adequate feedforward input. This is discussed in Results, Input patterns of LGN to V1.
Monocular versus binocular simulations.
In model simulations, both eyes were assumed to be open when the visual stimuli were presented. Because much of the experimental data used for validation was obtained with one eye occluded, we also simulated monocular stimulus presentation, and the results (data not shown) are very similar to the binocular case.
Summary of parameter values used in the simulations in this article
Leakage conductance, excitatory neurons, gRE = 50/s; leakage conductance, inhibitory neurons, gRI = 66/s; coupling weight LGN → E cells, SE,LGN = 0.059; coupling weight LGN → I cells, SI,LGN = 0.084; coupling weight cortical E → E, SEE = 0.028; coupling weight cortical E → I, SIE = 0.0095; coupling weight cortical I → E, SEI = 0.056; coupling weight cortical I → I, SII = uniformly distributed in (0.036, 0.048); coupling weight layer 6 → E, S6EE = 0.014; coupling weight layer 6 → I, S6IE = 0.0047; coupling weight ambient, Samb = 0.01; ambient firing rate, rEA = rIA = 250/s; AMPA fraction of E conductance, ρAE = 0.8; NMDA fraction of E conductance, ρNE = 0.2; AMPA fraction of I conductance, ρAI = 0.67; NMDA fraction of I conductance, ρNI = 0.33; and peak-driven firing rate in layer 6, f6,max ∼ 40/s.
Results
Input patterns of LGN to V1
Calculation of sparse LGN magnocellular input to V1 cortex layer 4Cα
An important new feature of our model is that it takes into account the very sparse LGN input to macaque V1. That the LGN input is extremely sparse can be gleaned from known anatomical data. According to Silveira and Perry (1991), the M retinal ganglion cell density at 5° eccentricity is 3500/mm2. If we take the conversion factor of [200 μm (on the retina)]/[degree of visual angle] for the Java monkey Macaca fascicularis (Silveira and Perry, 1991), then M cell density = 140 M ganglion cells/deg2 at 5° eccentricity. The cortical magnification factor in M. fascicularis V1 at 5° eccentricity is ∼2 mm/deg (Xing et al., 2009). Therefore, one hypercolumn in V1 cortex (0.5 × 0.5 mm) at 5° eccentricity corresponds to approximately (0.25°)2 = 1/16 deg2. Then the number of M ganglion cells that provide direct input to one hypercolumn would be 1/16 * 140, or approximately nine cells. We round that number off to 10. A similar calculation can be performed with the macaque LGN data from the study by Connolly and Van Essen (1984). They provide an estimate of the magnocellular LGN density at 5° of ∼160/deg2, which is in close agreement with the M retinal ganglion cell density estimate by Silveira and Perry (1991) of 140/deg2. This yields the estimate of 10 magnocellular LGN cells providing direct LGN input to each hypercolumn. The data from the study by Connelly and Van Essen (1984) indicate that a similar estimate would apply also at 10° eccentricity. The concordance of data between the studies by Silveira and Perry (1991) and Connelly and Van Essen (1984) suggests that within the central 10° of the visual field the same number of M retinal ganglion cells and magnocellular cells represent the same region of visual space.
With approximately 10 M retinal ganglion cells and a corresponding number of magnocellular LGN cells providing direct input to each HC, branching of magnocellular axons from neighboring HCs (Blasdel and Lund, 1983; Lund et al., 1993; Angelucci and Sainsbury, 2006) provides another 20 or so LGN inputs to each HC, making a total of no more than 30 thalamic inputs that drive the 4000 neurons in the HC of the model. Furthermore, Angelucci and Sainsbury (2006) reported that in layer 4Cα, a region 300 μm in diameter is contacted by no more than 11 LGN cells, and they proposed that the numbers of LGN inputs to cortical cells are likely much smaller. Their estimate is consistent with the numbers above, as well as the RF profiles obtained by Yeh et al. (2009) (and C. I. Yeh, D. Xing, P. E. Williams, and R. M. Shapley, unpublished observations).
Given the sparseness of LGN input to V1, it is a challenge for any model to generate the full repertoire of the visual functions of V1. Our model met this challenge with a specific set of LGN input patterns, which were designed to be as consistent as possible with anatomical and physiological data. This is discussed in the Numbers and templates section below. We also found that the distribution of the MR and SF preferences of V1 neurons are strong constraints on the choice of LGN input patterns. This is discussed in the Modulation ratios and nLGN section below.
Numbers and templates
Following the optical imaging data from the study by Obermayer and Blasdel (1993), we sought to assign to each neuron in a designated orientation domain in layer 4Cα (Fig. 1B) a configuration of LGN cells that has the potential to produce two or three ON/OFF subregions aligned in the direction specified. This task has the following constraints: (1) each ON (or OFF) subregion cannot be represented by more than one row of LGN cells, or the number of LGN inputs would be too large, contradicting the findings of the study by Angelucci and Sainsbury (2006); (2) for the same reason, each row cannot exceed two to three LGN cells; (3) because the peak response of V1 is to drifting gratings with SF at 2.5–3 c/d (cycles/deg) (De Valois et al., 1982), rows of LGN cells corresponding to adjacent subregions should be separated by three-sixteenths to one-quarter degree [in particular, choosing ON/OFF pairs that are nearest neighbors, as proposed by Ringach (2007) and Paik and Ringach (2011), would lead to a quite different SF preference]; and (4) because of the limited reach of LGN axonal terminal trees (Angelucci and Sainsbury, 2006), a V1 neuron cannot be connected to LGN cells that have direct projections that are >450 μm away.
LGN input and corticocortical coupling of E and I cells. A, Receptive field centers of M retinal ganglion cells or magnocellular LGN cells in the locations of their projections to V1 cortex. ON cells are indicated by white, and OFF cells by black, against a gray background. The ganglion cell density is from the study by Silveira and Perry (1991). The square drawn in the center corresponds to the retinal region that projects to a single hypercolumn (HC). B, The nine HCs of the model with the intended orientation map drawn in. C–E, Connectivity examples: presynaptic E input (red) to an E cell (black dot; C); presynaptic E inputs (red) to an I cell (D); and presynaptic I inputs (blue) to an E (or I) cell (E). The square in each plot represents an HC.
Permitting each cortical E cell to have 0–6 LGN inputs, arranged in two or three rows of one to three cells each, we found that constraints (1)–(4) above can be met by the LGN lattice in Figure 1A, but that there is only a small number of viable configurations. Some example LGN templates are shown in Figure 3. As shown in Figure 1B, we divided each HC into six distinct orientation domains, and stipulated that within each domain neurons received LGN inputs favoring one of the orientations of 0°, 30°, 60°, 90°, 120°, or 150°, with 0° taken to be vertical. Because the lattice we used has a threefold rotational symmetry, it is sufficient to enumerate two sets of admissible templates, one preferring, for example, horizontal and the other 30° from horizontal. Rotating these two sets of templates by 60° and 120°, one obtains templates for the other four orientations.
Determination of SIE/SEE. A, Dependence of background firing rates on SIE/SEE. The target mean background firing rate for E cells is 2.5–4 spikes/s (dashed lines), so the SIE/SEE ratio needs to be 0.3–0.32. Note that in this window, I firing rates are three to four times the E firing rates. B, Distribution of the bg firing rates of the population of E neurons with SIE/SEE = 0.3.
LGN → V1 connectivity. A, Templates of LGN–V1 cell connections for excitatory cells in a vertical-preferring patch. All viable configurations of ON (white) and OFF (black) LGN cells are enumerated, ordered by increasing numbers of cells in the configuration (nLGN). B, Templates for cells in a patch that prefers Horizontal orientation, ordered similarly. The model used templates for 0°, 30°, 60°, 90°, 120°, and 150°, obtained by rotating the two sets of templates shown.
We used the following algorithm for assigning to each cortical neuron a set of afferent LGN cells: first pick the number of LGN inputs (nLGN) randomly according to a distribution the determination of which we postpone to the Modulation ratios and nLGN section below. Then, from the list of templates in the intended orientation of the cortical neuron and nLGN (Fig. 3), pick a template at random, with equal probability of picking templates with opposite polarity. Finally, go to the LGN sheet for the corresponding eye, and select a configuration of LGN cells matching this template, making sure that it projects to a location in cortex within 400–500 μm of the cortical neuron in question.
Cells near the border of two orientation domains have positive probabilities of being assigned templates with either orientation. These probabilities are determined by a narrow Gaussian function: at 10 μm from the boundary, the probability of picking the “wrong” template is 0.38; at 30 μm away, the probability drops to 0.01.
As for I cells, it does not matter a great deal whether one uses oriented or unoriented (or scrambled) LGN inputs. We have tried both, and the results of model simulations were similar. In both cases, the I cells showed a significantly higher circular variance (CircVar) than E cells (see Results, Tuning and orientation-specific responses of model neurons), and the rest of the model was not substantially affected. In the simulations shown in this article, we used unoriented I cells, with nLGN ranging from zero to eight, its distribution given by a Gaussian with mean = 3, SD = 2.
Modulation ratios and nLGN
We found that getting the model to match the distribution of the MR across the V1 population was a strong constraint on the number of LGN inputs (nLGN) and, therefore, on the allocation of LGN templates. The MR of the responses of a cortical cell has been taken to indicate the degree to which the cell sums its inputs linearly (De Valois et al., 1982; Skottun et al., 1991). Many experimentalists have measured MR and reported its diversity in V1 (Ringach et al., 2002). Here MR is defined to be F1/F0, where F1 is the amplitude of the fundamental frequency of the response of the cell when it is driven periodically, and F0 is the mean. Specifically, MR = 0 corresponds to a cycle average that is constant over time; a sine wave function with amplitude equal to its average yields an MR = 1, and MR = 2 corresponds in practice to all spikes being concentrated in a single bin in each cycle average histogram. Following the studies by De Valois et al. (1982) and Skottun et al. (1991), we call a model cell “complex” if its MR is <1 when driven by its optimal grating at full contrast, and “simple” when MR is >1.
It is commonly thought—and receptive field mappings appear to support the hypothesis—that complex cells correspond to those V1 neurons with either no or few LGN inputs. In our model, such cells tend to be suppressed by cells with larger numbers of LGN afferents if nothing is done to boost their firing rates, yet it is an empirical fact that when driven, complex cells have firing rates that are higher than average [Ringach et al., 2002; http://www.ringachlab.net/lab/Data.html (the Ringach laboratory web site)]. This suggests that some boosting is necessary. We will return to this point later, preferring to first proceed to our main observation about the dependence of MR on nLGN.
With the firing rates of low nLGN cells adequately boosted (through increased connectivity or synaptic weights), we found that there is a strong correlation between the nLGN of a V1 cell and its MR. To illustrate this point, we presented the model with a battery of 64 gratings in eight different orientations, each with eight different SF values. In these model simulations of cortical responses, as for all of the other results of simulations in this article, the responses of the model were obtained for 20 s of visual stimulation. The MR of a cell was taken to be its MR for the grating to which it had the strongest response. When cells in the center HC of our model were divided according to nLGN and histograms of MR in each category were plotted, we observed a steady upward drift in MR as nLGN was increased. This phenomenon appears to be robust: it is independent of specific choices of parameters that affect mean firing rates (within a reasonable range) or of the allocation of nLGN among cells in the model.
Experimental results (Ringach et al., 2002; M. J. Hawken, unpublished observations) show that in layer 4Cα of the macaque V1 cortex, the distribution of the MR is bimodal, with ∼30% of the cells being complex. To carve out a bimodal distribution emulating data on relative incidence of simple/complex cells, we found it necessary to use nLGN distributions that approximately reflect these percentages, that is, to give approximately one-third of cortical cells zero to two LGN inputs, and approximately two-thirds of them, four to six LGN inputs. Using this as a guide, we located a range of nLGN distributions that emulate data, with some emulating data better than others.
The allocation of nLGN used in the model simulations presented is shown in Figure 4. Here complex cells correspond mostly to V1 cells with two LGN cells with ON–OFF pairs randomly oriented and positioned as close to one another as possible. We chose to give two LGN inputs to cells that would likely turn out to be complex because of the results of Hoffman and Stone (1971), Tanaka (1983), and Hirsch et al. (2003), who suggested that some complex cells have LGN inputs. However, others have reported no direct input of LGN to complex cells (Alonso et al., 2001). To find out whether or not it mattered that complex cells received nLGN = 2 or 0, we also ran model simulations in which complex cells corresponded primarily to cells with nLGN = 0, and those simulations worked as well as the nLGN allocation scenario reported here.
Number of LGN inputs (nLGN) and modulation ratio (MR) distributions. Allocation of nLGN inputs used in the model, chosen to produce a bimodal MR distribution that emulates data. Top seven panels, Histograms of MR in the central HC plotted in a vertical column separately for neurons with nLGN = 0, 1, …, 6. Bottom, Sum of the seven histograms yields the bimodal distribution of MR for all cells in the central HC.
We now return to the issue of firing rates. Other things being equal, driven firing rates in our model increase with increasing nLGN. This is reasonable, as V1 is a high-gain system and sensitivity to additional external drive is a hallmark of such systems. From the empirical fact that complex cells fire at a higher rate, we infer that some corrective mechanism is present. One possibility is that, via some form of homeostasis, neurons with less LGN input receive more excitation from other sources. We implemented this in the model: E neurons with fewer LGN inputs were given higher cortical input, from within 4Cα and L6. The idea of compensating a cell for low LGN inputs with higher cortical input was borrowed from the study by Tao et al. (2004), but, unlike the model in Tao et al. (2004), all of the neurons in our model receive most of their excitatory synaptic current from other cortical cells; therefore, the boosts of low nLGN cells are fractionally much smaller.
Specifically, the following enhancements in cortical inputs were used. Within layer 4Cα, peak connectivity (as expressed in the number of presynaptic E neurons) ranged from 13% to 18%, with lower corticocortical connectivity for cells with higher nLGN. Furthermore, V1 cells with zero to two LGN inputs had about twice as many presynaptic inputs from layer 6 as cells with five to six LGN inputs. Through exploration of many simulations, we found that the cortical enhancement of complex cells could be redistributed between 4Cα and L6, or between the amount of connectivity and synaptic weights, SEE and S6EE, as long as they produced the same combined enhancement of the firing rate.
The collection of templates in Figure 3 and the nLGN distribution in Figure 4 complete our model description.
Tuning and orientation-specific responses of model neurons
Now comes the test: does our model perform like a real V1 cortex? To test the model, we used as visual stimuli sine gratings of high contrast, drifting at a rate of 4 Hz, which is similar to stimuli used in many experimental studies of V1 function. The model data shown in Figures 5, 6, 7, 8, 9, and 10 are collected from the central HC (Fig. 1B, shading) to avoid boundary effects.
Orientation and SF tuning in the model. A, Five E cells from the horizontal (90°) patch of the central HC, driven at 4 Hz by a battery of gratings at 100% contrast. Left, <orientation, spatial frequency> matrix of peak responses. Middle, Orientation tuning curves with CircVar values displayed. Right, Cycle average firing rates with MRs. B, Five I cells from the same patch as in A, and their peak responses, orientation, and spatial frequency tuning as illustrated in A. The dashed horizontal lines indicate the mean background firing rates.
Histograms of firing rates, orientation, and spatial frequency tuning in the central HC. The top three rows show histograms for excitatory complex cells, excitatory simple cells, and inhibitory cells. Bottom row shows distributions of CircVar and of the SF preference of the synaptic current from LGN to all E cells. (No meaningful notion of firing rate distribution for LGN current.)
Orientation tuning curves of two representative model E cells and their synaptic inputs. The orientation tuning of the peak of the averaged response is plotted vs orientation. Top row, A simple E cell; bottom row, a complex E cell. Both cells are drawn from the same vertical-preferring patch of model 4Cα. Tuning curve, Orientation tuning of the spike rates of the cells (the E cells have CircVar ∼0.5); LGN, LGN input current as a function of orientation angle; nLGN = 5 for the simple cell and nLGN = 0 for the complex cell (i.e., the complex cell illustrated here receives zero LGN input current). The three panels on the right show input currents from various sources as functions of orientation. E, Recurrent cortical excitatory current; FB, cortical feedback excitatory current (from layer 6); I, cortical inhibitory current. Note that the inhibitory current is plotted as positive in magnitude, but its sign is opposite to that of the excitatory currents.
Continuity of orientation maps. A, Patch-averaged peak firing rates for E cells as a function of spatial frequency and orientation. The first five panels are for patches preferring 0°, 22.5°, 45°, 67.5°, and 90°; the sixth panel is located at the pinwheel center. The sizes and locations of these patches in the central HC are indicated (right). The grid square of each color map denotes a grating with specified orientation and SF; its color corresponds to peak rates averaged over all E cells in the patch. B, The seven panels show distributions of orientation preferences for patch neurons in seven wedges in the central HC corresponding, in an ideal orientation map, to regions preferring 30°, 35°, 40°, …, 60°, respectively (see main text and inset for detail). The histograms depict responses to 12 gratings 5° apart between 0° and 180°. The red vertical line in each plot denotes the hoped for average orientation preference. C, Plot of average visual angle preference as a function of angle around the pinwheel center.
Dynamics of V1 population impulse trains and mean firing rates. A, Rasters of spikes of a target V1 E cell and its four LGN inputs. These are responses to an optimal drifting grating, drifting at 4 Hz (250 ms period). White dots, ON center LGN spikes; black dots, OFF center LGN spikes; red dots, V1 spikes. B–D, Rasters of background spike firings in the cortical patch that prefers vertical gratings in the central HC (B); rasters when the same patch is driven at the optimal (vertical) orientation of the patch (C), driven by the orthogonal orientation (D). B–D, x-axis is time; y-axis is neuron index within the patch, sorted by number of LGN inputs; lower indices correspond to fewer LGN inputs. E cell spikes are red dots; I cell spikes are blue dots. E, Further information for the regime depicted in C. Top, Summed spike activity of the E population as a function of time, shown in 5 ms windows. The vertical scale of the plot is the fraction of the total number of cells, illustrating that, in what may appear to be synchronous spiking events in C, only a small fraction of all E cells are active. Bottom panels, Histograms of mean firing rates in impulses per second for E and I neurons in the same patch, driven in the optimal orientation.
Intracellular voltages and conductances of driven E and I cells. A, Intracellular voltage [in normalized units (0,1); see Materials and Methods, Equations of neuronal dynamics], E conductance, and I conductance for one representative simple E cell in the vertical patch (as in Fig. 9). B, Same for one representative complex E cell in the patch. C, Same for one representative I cell in the same patch. Occurrences of spikes indicated by red dots in the voltage plots.
Effectiveness of sparse LGN coding
In the Introduction, we raised the issue of whether so few LGN inputs to a V1 HC can confer OS on V1 neurons. Indeed, there is reason for concern: from Results, Input patterns of LGN to V1 (Fig. 4), the number of LGN inputs to a cortical neuron is, on average, ∼4 compared with ∼200 presynaptic E neurons from layer 4Cα, and 30–50 presynaptic E neurons from L6. Even with SE,LGN ∼ 2 * SEE and driven LGN firing rates approximately three times that of cortical neurons, the fraction of input currents to a 4Cα neuron that is feedforward is small—in fact the simulations indicate that the excitatory drive from LGN is <20% of the total in simple cells, and <10% in complex cells.
Nevertheless, the model produces model E cells with the range of properties expected from experimental data. Figure 5 shows the OS properties of example neurons from the model. The color maps in Figure 5A show the responses of five representative model E cells from a horizontal-preferring patch. Their responses, measured in terms of peak firing rates, are to an orientation × SF matrix of gratings (De Valois et al., 1982), with orientation plotted on the vertical axis and SF on the horizontal axis. The eight orientations range between 0° and 180° evenly spaced, and the eight SF values range from 0.5 to 16 c/d. In the middle and right columns are the tuning curves and cycle average firing rate distributions of the same neurons.
The five E neurons shown in Figure 5 have peak firing rates ranging from 20 to 80 spikes/s, and they have quite varied tuning curves and profiles of cycle averages. If one compares these tuning curves to those of real V1 neurons (Ringach et al., 2002), one finds that model E cells have as much OS as in real cortex, with the same diversity: some of the model neurons are highly selective, while others have broad tuning curves.
Figure 5B demonstrates that model I cells are less selective for orientation than E cells. Though they show some selectivity, the five representative model I cells produce <ori, SF> response surfaces that are markedly less peaked than those of E cells (Nowak et al., 2008).
The model allows us to analyze OS across the population. Next, we present model population statistics, collected from the central HC of our nine-HC model cortex (Fig. 1B).
Firing rates
To demonstrate that the model emulates V1, we present in Figure 6A the peak firing rates of model neurons when driven by their preferred gratings (from among the 64 used for the simulations depicted in Fig. 5). The distributions in Figure 6A, which show complex E cells firing more than simple E cells (D. Ringach, R. M. Shapley, and M. J. Hawken, unpublished observations), and I firing rates three to five times those of E firing rates (Swadlow, 1988; Cardin et al., 2007), are consistent with data.
We also call attention to the wide range of firing rates seen in Figure 6A. The diversity of neuronal responses resembles what is found in V1 (Ringach et al., 2002) and is an emergent property in the model, meaning it is a property that was not programmed into the equations governing the dynamics of individual neurons (see Materials and Methods, Equations of neuronal dynamics) but occurs as a result of interactions among neurons. Analysis revealed the source of the observed diversity as caused by the fact that small random biases in connectivities and synaptic weights are magnified through neuronal interactions, exaggerating—through recurrent inhibition—the dominance of neurons that receive slightly more E input or slightly less I input. Competition among neurons drives much of the dynamics in this model, and we propose that the same may be true in the real cortex, as indicated by the data from the study by Monier et al. (2003).
Circular variance
Circular variance (CircVar) is a widely used measure of OS (Mardia, 1972; Batschelet, 1981; Ringach et al., 2002). CircVar = 1 for nonselective neurons, and CircVar = 0 for perfectly selective neurons. As one can see in Figure 6B, the mean CircVar for E neurons is between 0.6 and 0.7, which is in agreement with data. Equally if not more striking are the spreads of the distributions of CircVar in the model. They show great diversity: some neurons are very broadly tuned, while others are quite sharply tuned. Such diversity is also characteristic of real V1: both mean and variance are consistent with experimental data (Ringach et al., 2002; Zhu et al., 2010). As with firing rates, diversity in CircVar in the model is due in part to the probabilistic connectivity between cortical cells (Fig. 1C,D).
Model I cells have some but significantly less OS (higher CircVar) than E cells, which is also consistent with available data (Nowak et al., 2008). Among the reasons why I cells in the model are less tuned than E cells even though they receive excitatory input from a comparable region of the HC (Fig. 1C,D) may be (1) that I cells receive unoriented LGN inputs, and (2) suppression of responses to the orthogonal-to-preferred orientation has a stronger effect on E cells than on I cells because SII is smaller than SEI.
Least tuned of all are LGN input currents to cortical cells, with a CircVar of ∼0.8–0.9; this is discussed in more detail in conjunction with Figure 7.
Bandwidth is another measure of OS that is widely used (De Valois et al., 1982). We calculated bandwidth distributions for model neurons (simulation results not shown) and found very good agreement quantitatively, both in mean values and in the variance, with experimental data (De Valois et al., 1982; Ringach et al., 2002).
SF preferences
Histograms of SF preferences are shown in Figure 6C. As is consistent with data (De Valois et al., 1982), the overall preferred SF is 2–3 c/d. This property is effectively built into the model by the spatial frequency tuning function for LGN cells (Materials and Methods, Equations of neuronal dynamics; Zhu et al., 2010) and the distances between ON and OFF subregions in the LGN templates (Figs. 1, 3). Model complex cells prefer a slightly higher SF due to the geometry of their LGN afferents: tightly packed two LGN cells respond better to slightly higher SF (see Results, Input patterns of LGN to V1). This model result also is consistent with data from the study by De Valois et al. (1982).
Sources of input currents to a cortical cell
Having established that our model neurons enjoy OS, we now analyze the model to understand how it came about. Figure 7 displays the tuning curves of spike rates and currents that go into two typical cortical neurons, one simple E cell and one complex E cell that were chosen to have approximately the same moderate amount of OS (CircVar ∼0.5). The top row of Figure 7 shows information for the simple cell. Its orientation-tuning curve for peak spike rate is plotted in the leftmost panel of Figure 7. The graph labeled “LGN” shows its summed LGN input current at the peak of the cycle-averaged response. LGN current is only weakly OS, which is consistent with the statistics in Figure 6. Weak OS of LGN input in our model is not surprising: given that the RF subregions in our LGN templates are so short and so far apart (Fig. 3), an LGN configuration of, say, 3-ON, 2-OFF responds to gratings of many different orientations. The improvement in OS is caused by a combination of the strong effect of tuned recurrent excitation from neighboring E cells (Fig. 7, panel labeled “E”) that have spike rate-tuning curves like the tuning curve displayed on the left, and the effect of recurrent inhibition (Fig. 7, panel labeled “I”), which tends to suppress responses at nonpreferred orientations (Zhu et al., 2010; Xing et al., 2011). OS is further enhanced by cortical dynamics, by the temporal synchronization of recurrent excitation (Fig. 9). In addition, feedback from L6 (Fig. 7, “FB”) slightly augments OS.
The bottom row in Figure 7 shows the corresponding plots for a complex cortical neuron with nLGN = 0 and low MR. We have chosen such a cell to illustrate the point that, even with no direct LGN input, a cell can acquire OS from neighboring cells via corticocortical interactions. One can observe in this case the strong, tuned, recurrent excitatory drive from neighboring E cells (Fig. 7, panel labeled “E”). We point out that ours is the only realistic model of the visual cortex that provides a robust mechanism for OS in complex cells. Complex cells in the Tao et al. (2004, 2006) and Zhu et al. (2010) models did not have much OS, and in those of McLaughlin et al. (2000) and Wielaard et al. (2001) there were no complex cells to speak of.
To summarize, our model shows that to cultivate a property like OS it is necessary only to “seed” it, that is, to introduce a weak form of OS to a fraction of the population (e.g., those cortical cells with four or more LGN inputs), and the effects will be amplified through corticocortical interactions. Similar views have been expressed before (Ben-Yishai et al., 1995; Hansel and Sompolinsky, 1997). Here we demonstrate the idea in a realistic model, with robust OS in both simple and complex cells emerging as a result of corticocortical interactions.
Continuity of orientation maps
An important question raised in the Introduction is whether or not the model can generate orientation preferences that are intermediate between the small number of orientations that can be specified by the known sparse LGN connectivity. Simulation results, using the same 8 × 8 matrix of <orientation × SF> gratings as in the results shown in Figure 5, answered this question in the affirmative (Fig. 8). In the <orientation × SF> matrix, we deliberately used grating orientations that were not commensurate with those in the LGN templates (Results, Input patterns of LGN to V1). The templates have a threefold rotational symmetry, while the gratings used in the experimental simulations represent eight orientations at equal angles apart. Figure 8A shows the peak firing rates in response to the <orientation × SF> matrix averaged over several hundred neurons in six different patches, the precise locations of which in our model cortex are also shown. Except for the pinwheel-center patch, each of the other patches shows a clear preference for an orientation and a range of SF values. The orientation preference coincides with the location of the patch within the HC, and the SF preference peaks at 2–3 c/d, consistent with macaque V1 data (De Valois et al., 1982). Since no cortical cell has LGN afferents aligned with grating stimuli at 45° or 135°, the only way that some neurons can acquire orientation preferences in these directions is through the effects of corticocortical interactions.
To take a closer look at the ability of our model to produce orientation preferences not represented by LGN templates, we focus on two adjacent regions of the HC, one containing neurons assigned LGN templates aligned with 30°, and the other templates aligned with 60°; the two regions are separated by a narrow corridor in which some neurons receive the 30° and some the 60° templates (see Results, Input patterns of LGN to V1). The question of interest here is, how well can the model “see” the visual angles in between?
Figure 8B shows statistics for orientation preferences when the model is presented with 12 gratings, aligned with orientations evenly spaced between 0° and 180°. The seven panels show histograms of preferences for neurons located in seven different (narrow) wedges of the HC, chosen so that, with an ideal orientation map, they would have orientation preferences for 25–35°, 30–40°, 35–45°, …, 55–65°, respectively. In the first and last panels, it is not surprising that 30° and 60° are preferred, given that they coincide with LGN template assignments. But observe how, in increments of 5 visual degrees, the peaks of the histograms shift gradually yet definitively to the right. In the third panel of Figure 8B, for example, all of the neurons receive 30° templates, yet, due to their relative proximity to the 60° region there is an unmistakable bias in that direction. Neurons in the fourth panel of Figure 8B have the strongest preference for 45°, even though there are no 45° LGN templates. The only way that these graded preferences could come about is through corticocortical interactions.
Notice also the diversity in preferences within each region in Figure 8B. While the neurons collectively have a clear orientation preference, individually they show great diversity in their preferences. Many cells deviate significantly from the mean of the local population. The diversity of orientation preference is also an emergent property; it comes about as a consequence of the biases in cortical input coupled with the dominance of recurrent and feedback excitation in the model. These features depend on the probabilistic connectivity between cortical cells described in Figure 1, C and D.
Further evidence of smoothing of the orientation map is displayed in Figure 8C, which plots vector means of the most preferred orientations for 36 wedges organized around the pinwheel center of the HC. Each wedge spans 10° in visual angle with overlapping windows, as in Figure 8B. Evidently our model network was able to smooth out the very sparse information passed along by LGN to create a continuous orientation map without abrupt jumps in orientation preference that might be expected from the sparse LGN input alone.
Robustness of OS
The simulation results in this article assume that SE,LGN = 2.1 * SEE following the study by Stratford et al. (1996); see Materials and Methods, Parameter determination. Other articles have reported different ratios of thalamocortical-to-corticocortical synaptic strength; e.g., it was found to be ∼1 for cat S1 (Schoonover et al., 2014) and ∼4 for mouse S1 (Gil et al., 1999). Therefore, we thought it was necessary to explore whether or not varying the SE,LGN/SEE ratio would affect OS adversely, and found that it does not. With other parameters set to be similar to those in the present model, we found that setting SE,LGN/SEE = 1 causes firing rates to decrease and the fraction of complex cells to increase, and vice versa when SE,LGN/SEE = 3. However, when suitably rescaled, the orientation tuning and orientation maps, as shown in the plots in Figure 8A, remain substantially unchanged. This proves that OS is extremely robust in this kind of model with approximately balanced excitation and inhibition.
Temporal dynamics in background and under drive
The dynamics of spike firing in local populations is a crucial emergent property of the model that makes a significant contribution to OS. Raster plots of spike firing in 500 ms windows are shown in Figure 9B–D, for a patch of ∼400 E and I neurons in an orientation domain. Figure 9B depicts spike firing in background, and Figure 9C when the population is driven by a grating drifting at the optimal orientation, while Figure 9D shows how the same population responds to a grating drifting at the orthogonal-to-preferred direction. Notice the gamma rhythm (25–90 Hz) due to partial synchronization in the spike firing in Figure 9C. The periodicity of the visual drive (4 Hz) can be deduced from Figure 9A, which shows the spike firing of a target E cell and its four LGN inputs. Observe by comparing panels A and C in Figure 9 that gamma rhythms are tied neither to the spiking of LGN cells, which can fire at considerably higher rates (∼100 spikes/s), nor to the much slower temporal frequency of the drifting grating; they are generated within cortex with their own autonomous resonant frequencies. This is an emergent property of population activity.
In our model, synchronization of spike firing in the gamma band is graded and depends on stimulus orientation. These results are consistent with experimental data (Eckhorn et al., 1988; Gray and Singer, 1989; Volgushev et al., 2002; Kayser and König, 2004; Henrie and Shapley, 2005; Womelsdorf et al., 2012). The Fourier transform of the summed spike density of the neurons (data not shown) has a gamma peak that was high and narrow for the optimal stimulus, weaker and broader for the nonpreferred stimulus, and low and very broad under background conditions. Not only is the greater synchronization at the preferred orientation an orientation-specific response, it contributes to OS by amplifying the effect of recurrent excitation locally.
It is important to emphasize that, unlike the PING model (Whittington et al., 2000), synchronization events in our model involve only fractions of the population, and participation varies widely from neuron to neuron. The wide distribution of model firing rates shown in Figure 9E does not support synchronous spiking involving the full population. Notice in fact that most model neurons fire at rates too low to follow the gamma-band fluctuations. In this respect, the firing rate distributions of the model resemble many experimental datasets (for review, see Roxin et al., 2011; Buzsáki and Mizuseki, 2014; M. J. Hawken and R. Shapley, unpublished V1 observations). Even though it may appear from the raster plot in Figure 9C that the population firing is highly synchronous, only a fraction of the population participates in any firing event. Figure 9E, which gives the spike counts for E neurons in Figure 9C, shows clearly that no more than 25% or so of the E population fire spikes in any of the 5 ms windows.
For a general discussion of the dynamic mechanism behind such partial synchronization in neuronal networks, see the study by Chariker and Young (2015). Gamma rhythms are very important signatures of cortical dynamics, but in the interest of limiting the length of this article, we have elected to postpone a more detailed analysis to a forthcoming article.
To elucidate how gamma-band rhythms are related to the dynamics of individual neurons, we show in Figure 10 the voltage and conductance traces of three model neurons: two E cells, one simple and one complex (Fig. 10A,B); and an I cell (Fig. 10C). These traces strongly resemble data from intracellular recordings (Hasenstaub et al., 2005). Figure 10A–C show that fluctuations in conductance are entrained to gamma-band fluctuations, even though most of the rising peaks do not cause spikes to fire. Indeed, how often a neuron spikes, and when, is to some degree accidental, depending on whether a threshold is crossed. And that in turn has to do with the interplay between the E and I populations, even the arrival times of the spikes of presynaptic E/I neurons. Figure 10 demonstrates clearly that firing rates are determined by dynamics, and cannot be computed from the LIF equations and parameters alone, i.e., they are also emergent properties of the model.
Model predictions
A fraction of the simulation results we presented was for purposes of validating the model by comparison with existing data. Among the many model outputs that compare well with existing data are background and driven firing rates and their distributions (Figs. 2, 6, 9), modulation ratio (Fig. 4B), orientation and spatial frequency tuning of individual E and I cells (Fig. 5), population statistics on peak firing rates, CircVar and SF distributions (Fig. 6), and gamma-band synchronization that is somewhat tuned for orientation (Fig. 9).
A number of our model results go beyond the data that are currently available. We propose them as model predictions to be tested in future experiments, as follows.
Configurations of LGN cells presynaptic to cortical E cells, as proposed in Results, Input patterns of LGN to V1 (Fig. 3).
Weak OS of I cells: there are some data about weak OS of I cells in cat cortex (Nowak et al., 2008) and also in mouse visual cortex (Isaacson and Scanziani, 2011; Liu et al., 2011), but quantitative data about the distribution of OS in I cells are lacking.
The even weaker OS of LGN input currents to a cortical cell: data about mouse V1 indicate very weak OS of LGN input to V1 cells (Lien and Scanziani, 2013); LGN input to V1 in macaque is likely very different from that in mouse, but the model predicts that its OS also is very weak, as shown in Figure 6.
The sharpening of OS for 4Cα cells through interaction with other cortical neurons, with or without direct LGN input, as shown in Figure 7.
Properties of the HC orientation preference map: robust population averages vs the diversity of responses among individual neurons revealed in Figure 8, A and B.
The continuity of orientation preference maps, as depicted in Figure 8C (the present resolution of optical imaging is not sufficient to test predictions 5 and 6, and high-resolution measurements with two-photon imaging or another high-resolution technique will be needed).
The partial participation of individual cells in gamma rhythms, as evidenced by the irregular spiking and entrainment of the conductances shown in Figures 9 and 10. While some intracellular recordings in V1 in vivo (Volgushev et al., 2002) suggest the intermittent participation in the gamma rhythm that the model predicts, more data are needed to test this specific prediction.
What distinguishes our model from feedforward (or modified feedforward) models (e.g., Priebe and Ferster, 2008; Persi et al., 2011; Rubin et al., 2015; Goris et al., 2015) are the model predictions we have to offer. Predictions such as the sharpening of OS, higher OS in E cells than I cells, smooth orientation maps on fine spatial scales, orientation-tuned gamma synchronization, and diversity of OS and orientation preferences within local populations are based on model phenomena that emerge as a result of dynamic cortical interactions. Feedforward models, which downplay dynamic interactions, generally do not enable such predictions.
Discussion
The ability of the model cortex to (1) sharpen weakly tuned orientation inputs (Figs. 5, 6, 7) and (2) fill in the gaps of orientations not represented by LGN (Fig. 10) is very strong evidence in support of the hypothesis that, although visual signals in the brain originate from the retina and LGN, dynamic interactions of cortical neurons play a major role in the processing of this information. In our model, OS is extremely robust, and does not require any fine-tuning of parameters as long as excitation and inhibition are roughly balanced.
How does the model work and what insight does this give us into cortical function in general?
Much of the visual processing leading to OS occurs within the cortex as a result of the dynamic interactions among cortical neurons. Therefore, these dynamic interactions are a focal point of our investigation. We have seen in Results (Tuning and orientation-specific responses of model neurons) how, through recurrent excitation, cortex produces a coherent and robust orientation preference in local populations, derived from “seeds” of orientation bias transmitted through weak feedforward currents. Another useful function of corticocortical interactions is that they modulate the effects of sparse LGN inputs by filling in potential gaps in the orientation map. No less important is the constant competition between the E and I populations. The precarious balancing act between these two populations puts the system in a state of high gain, making it very sensitive to small changes in input, as we saw in the analysis of the MR distribution. Furthermore, recurrent excitation coupled with recurrent inhibition helps to produce diversity in neuronal responses enabling cortex to discriminate among a wider range of stimuli (Kang et al., 2004).
Through dynamic interactions among neurons, our network model exhibits behaviors far more complex than can be expected from the LIF equations (Materials and Methods, Equations of neuronal dynamics) that govern the dynamics of individual cells. For a start, there is a great deal of balancing and regulation within the cortical network. When stimulated, E neurons magnify the drive at the same time that I neurons put things in check. Properties (e.g., OS) conferred upon a fraction of the population are shared with the rest of the population, yet the same system that has an averaging effect on its constituent cells (as in the smoothing of the orientation map) also produces great diversities in firing rates, in circular variance, in spatial frequency preferences, and in modulation ratios. The diversity of functional properties is emergent, i.e., it cannot be deduced from the behaviors of single neurons or even a handful of neurons but arises only through the interaction of the many constituent parts of the network.
Emergent behaviors are the hallmark of complex dynamical systems. As our model outputs closely resemble many sets of experimental data, we propose that the visual cortex is, like our model, a complex dynamical system. As such, it is driven largely by population activity, and to understand it requires that we study not only individual neurons but also population dynamics, in particular dynamic interactions among neurons connected by local and medium range circuitries.
Orientation-specific dynamic responses
Our model generates partially synchronized firing in the gamma band (30–90 Hz) when driven by drifting gratings (Fig. 9A–C). The response is stimulus dependent and orientation specific: there is very little synchronized firing in the background, but strong synchrony was evident at high contrast in domains preferring an orientation aligned with the grating. Synchronized spiking is more effective in producing recurrent excitation locally, and indirect inhibition in a slightly larger area through the widespread excitation of I neurons by E neurons (Fig. 1D). Therefore, gamma-band synchronization at the preferred orientation is a very effective way for the cortical network to enhance OS.
As mentioned in Results, others have conjectured that gamma-band rhythms contribute to OS, based on cortical data (Volgushev et al., 2002; Womelsdorf et al., 2012) and on an abstract model (de Almeida et al., 2009). Now we have shown in our realistic, data-driven model that gamma-band dynamics are closely linked to OS. This is an entirely emergent property of the model. The dynamical characteristics observed in the model are very similar to data on visually induced gamma-band activity in real V1 (Eckhorn et al., 1988; Gray and Singer, 1989; Volgushev et al., 2002; Kayser and König, 2004). Notice also that the synchronization of the model network is asymmetric; the I cells are much more synchronized to the population than are the E cells (Fig. 9), which is another emergent model property that resembles V1 data (Hasenstaub et al., 2005).
Unlike the earlier PING model (Whittington et al., 2000), our model predicts that, in gamma-band fluctuations, only small fractions of the local population are involved in each firing event, and that individual neurons exhibit quite irregular behaviors despite the quasi-oscillatory activity seen on the population level. These predictions echo those of Rangan and Young (2013) and Chariker and Young (2015).
Comparison with earlier work
To our knowledge, the model presented here is the first realistic model that has modeled accurately the sparseness of magnocellular LGN inputs to the macaque V1 cortex and also successfully derived OS directly from these inputs. Earlier, more abstract models proposed that recurrent cortical interactions could produce high levels of OS even when the OS of the LGN input was weak (Douglas et al., 1995; Ben-Yishai et al., 1995; Somers et al., 1995; Hansel and Sompolinsky, 1997). Our model supports these ideas but goes beyond them in its realistic depiction of the anatomy, network dynamics, and the linkage between network dynamics and cortical function.
The present model is very different and much improved compared with earlier large-scale models, such as those of McLaughlin et al. (2000) and Tao et al. (2004, 2006). First, our focus on dynamics and on population activity is novel; this represents a paradigm shift that we believe will have implications for future directions of the field (see below). Second, the present model incorporates much more realistic neuroanatomy. In addition to sparse LGN → V1 connectivity, our modeling of the sparse E–E connections and denser I connections is data driven, as is our modeling of the feedback from layer 6; we have used biological guidance in virtually all aspects of our modeling when it is available. Third, our parameter determination process is more systematic, and it has been made transparent (Materials and Methods, Parameter determination). Last but not least, the present model is more comprehensive than all previous ones that we know of. We sought not just to replicate one or two aspects of V1 behavior but to reproduce, in a single model and from a single set of parameters, many V1 phenomena, such as orientation and spatial frequency selectivity, spontaneous and driven firing rates, gamma rhythms, simple versus complex cells, and neuronal diversity in all aspects of the above.
The present model differs also from rate models of V1 function, such as those of Priebe and Ferster (2008), Persi et al. (2011), Rubin et al. (2015), and Goris et al. (2015), in its emphasis on temporal dynamics as opposed to the use of fixed transducer functions as a means to account for cortical transformation of visual inputs into the population distributions of spike-firing rates. We appreciate the mathematical appeal of filters and simplified rate models, but approaches of this type downplay dynamic interactions among cortical neurons. As we have shown, such interactions are an integral part of cortical processing.
Future directions in computational neuroscience
A comprehensive, large-scale model of the cortex is not without cost. In particular, it cannot be as simple as phenomenological models that focus on one or two sets of behaviors at a time, or on the properties of single neurons. Comprehensive and phenomenological models serve very different purposes: phenomenological models typically seek to suggest analogies or to offer simplified mathematical descriptions. Comprehensive population models, such as the one we have presented here, seek to link cellular properties and network structure to dynamics and function, the ultimate goal being to use these models to test hypotheses and to suggest future experiments. We propose that for areas of the brain about which there are sufficient data, such as the visual cortex, it is time to move to next-generation models that are more comprehensive, data driven, and dynamic. Such a move would constitute a paradigm shift in computational neuroscience, and the present model is a step in that direction.
With regard to other regions of the cerebral cortex, much more anatomical and functional data than are currently available are needed for detailed modeling. However, the linkage between dynamics and function is not limited to V1 cortex. In many cortical areas, cortical activation is associated with increased gamma-band power in the local field potential (motor cortex, Murthy and Fetz, 1996; auditory cortex, Brosch et al., 2002; parietal cortex, Pesaran et al., 2002; prefrontal cortex, Gonzalez-Burgos and Lewis, 2008), as in V1 (Eckhorn et al., 1988; Gray and Singer, 1989; Volgushev et al., 2002; Henrie and Shapley, 2005; Ray and Maunsell, 2010). We propose that analysis of the connections among anatomy, function, and dynamics in the visual cortex may reveal principles that generalize to the entire cerebral cortex.
Footnotes
This work was supported in part by National Science Foundation Grant DMS-1363161 to L.-S.Y and also by a grant from the Swartz Foundation. This research was performed on the High Performance Computing resources at New York University. We thank Kevin Lin for preliminary studies; and Luiz Carlos Silveira, Michael Hawken, Dario Ringach, I-Chun Lin, Alessandra Angelucci, and Luke Hallum for helpful discussions.
The authors declare no competing financial interests.
- Correspondence should be addressed to Lai-Sang Young, Courant Institute of Mathematical Sciences, 251 Mercer Street, New York University, New York NY 10012. lsy{at}cims.nyu.edu
