## Abstract

In V1, local circuitry depends on the position in the orientation map: close to pinwheel centers, recurrent inputs show variable orientation preferences; within iso-orientation domains, inputs are relatively uniformly tuned. Physiological properties such as cell's membrane potentials, spike outputs, and temporal characteristics change systematically with map location. We investigate in a firing rate and a Hodgkin–Huxley network model what constraints these tuning characteristics of V1 neurons impose on the cortical operating regime. Systematically varying the strength of both recurrent excitation and inhibition, we test a wide range of model classes and find the likely models to account for the experimental observations. We show that recent intracellular and extracellular recordings from cat V1 provide the strongest evidence for a regime where excitatory and inhibitory recurrent inputs are balanced and dominate the feed-forward input. Our results are robust against changes in model assumptions such as spatial extent and strength of lateral inhibition. Intriguingly, the most likely recurrent regime is in a region of parameter space where small changes have large effects on the network dynamics, and it is close to a regime of “runaway excitation,” where the network shows strong self-sustained activity. This could make the cortical response particularly sensitive to modulation.

## Introduction

Neurons in the sensory cortices compute representations of the environment relevant for perception and action. These computations are modulated by the spatial (e.g., Levitt and Lund 1997; Seriès et al. 2003; Schwartz et al. 2007), temporal (e.g., Dragoi et al. 2000, 2001; Clifford et al. 2007; Kohn 2007), and behavioral (e.g., Lamme et al. 1998; Sharma et al. 2003; Bar 2004) context and are based on the integration of signals received via afferent and local recurrent connections. However, the relative contributions of afferent and recurrent inputs to cortical computation as well as their contribution to the emergence of tuned neuronal responses have long been a matter of debate. Orientation selectivity in the primary visual cortex (V1) is a paradigmatic example of such a computation, and the last 4 decades witnessed a vivid and polarized discussion about the underlying cortical mechanisms (Reid and Alonso 1996; Sompolinsky and Shapley 1997; Ferster and Miller 2000; Ringach et al. 2003; Finn et al. 2007)—in particular about the influence of the local excitatory and inhibitory recurrent synaptic connections and their interplay with the afferent drive (Martin 2002). A range of studies have provided computational models in support of different cortical mechanisms (for a review, see e.g., Teich and Qian 2006). Modeling studies (Somers et al. 1995; Hansel and Sompolinsky 1996; Troyer et al. 1998; McLaughlin et al. 2000; Kang et al. 2003) have also demonstrated that the properties of cortical networks, and hence the way information is processed in cortex, can change dramatically with the cortical operating regime, characterized by the relative strengths of the afferent and the recurrent inputs.

The recurrent circuitry of cortical networks can vary widely within an area, depending on its functional architecture. For example, the distribution of orientation preferences of neurons neighboring a particular V1 neuron depends on that neuron's location in the orientation map (Fig. 1A). We can distinguish 2 different extremes of regions in this map, those close to the singularities (pinwheel centers), where neurons with most or all the preferred orientations are represented in a small neighborhood, and those regions, where one particular preferred orientation dominates and only varies slowly with location (orientation domains). The location dependence of neuronal properties has been demonstrated experimentally for neurons lying along the continuum between these extremes: Neurons close to pinwheel centers have a more broadly tuned membrane potential (*V*_{m}) when compared with neurons in orientation domains (Fig. 1B, see also Schummers et al. 2002) and more broadly tuned excitatory (*g*_{e}) and inhibitory (*g*_{i}) total conductances (Fig. 1B, see also Mariño et al. 2005). However, the firing rate (*f*) of V1 neurons is highly selective near pinwheel centers and in orientation domains (Fig. 1B, see also Maldonado et al. 1997; Ohki et al. 2006). When the temporal response dynamics of cat V1 neurons is analyzed using the “reverse correlation” technique (Schummers et al. 2007), it is also found to be similar in pinwheel centers and orientation domains. The variability of the response time course, however, is significantly larger in pinwheel center neurons (Fig. 1C). Here, we show that all these physiological observations strongly constrain the regime in which visual cortical networks are likely to operate.

We have demonstrated previously that the experimentally measured response properties of cat V1 neurons are consistent with the predictions of a Hodgkin–Huxley network model dominated by recurrent interactions and with balanced contributions from excitation and inhibition (Mariño et al. 2005; Schummers et al. 2007). However, this cannot rule out alternative cortical operating regimes. Here, we turn to model-based data analysis in order to assess a continuum of network models that encompasses the full range from feed-forward via inhibition- and excitation-dominated models to models with excitation and inhibition in balance. Thus, in contrast to traditional modeling approaches, we determine the complete space of models able to account for the data, rather than demonstrating one model to be compatible with the data set.

We employ 2 generic model classes of different complexity: an analytically tractable firing rate network model (cf. Kang et al. 2003) and a physiologically more realistic Hodgkin–Huxley-based network model. Using a Bayesian approach, we calculate how much evidence the experimental data provide for different network parameters, in particular for different strengths of the afferent versus the lateral excitatory and inhibitory connections. Using data on the tuning of the neurons’ spike output, their membrane potential, their total excitatory, and their total inhibitory input conductance, we find that the experimental data provide a strong support for only one regime. This regime is characterized by significant excitatory and inhibitory inputs, dominating the afferent input. We then investigate how the strength of lateral excitatory and inhibitory connections in the Hodgkin–Huxley network model affects the response dynamics of the model neurons. Comparing simulated reverse correlation experiments with in vivo data, we find that only a strongly recurrent regime can replicate both the response dynamics and the observed variability of V1 neurons. Thus, also the analysis of the response dynamics points to the same recurrent operating regime.

## Methods

### Orientation Selectivity Index

Orientation tuning was analyzed using the orientation selectivity index (OSI; Swindale 1998), which is given by

*R*(ϕ

*) is the value of the quantity whose tuning is to be analyzed (e.g., the spiking activity) in response to a grating stimulus of orientation ϕ*

_{i}*. For all measurements, the stimulus orientations ϕ*

_{i}*,*

_{i}*i*= 1, …,

*N*, are uniformly distributed over 0°–180°. Then the OSI is a measure of tuning sharpness ranging from 0 (unselective) to 1 (perfectly selective). In addition, the OSI was used to characterize the sharpness of the recurrent input a cell receives based on the orientation preference map (cf. Figs. 1A and 2). To calculate this “map OSI,” we estimate the local orientation preference distribution by binning the orientation preference of all pixels within a radius of 250 μm around a cell into bins of 10° size; the number of cells in each bin replaced the quantity

*R*(ϕ

*).*

_{i}### Calculating the Bayesian Posterior

For the experimental data (Fig. 1B), the relationship between the OSI of the tuning for any of the quantities *g*_{e}, *g*_{i}, *V*_{m}, or *f* and the map OSI is fitted by a linear regression line with slope *sl*_{data} and intercept *ic*_{data}. For the analysis of our model results, we therefore accordingly assumed a linear relationship and described it with the 2 free parameters *sl* (slope) and *ic* (intercept): OSI_{tuning}=*sl* · OSI_{input}+*ic*.

For every set of model parameters, we then calculate the likelihood of the specific tuning of the *n* = 18 experimental data points given the model's slope and intercept under the assumption of Gaussian additive noise

The standard deviation (SD) of the noise was set to the SD of the experimental data around the linear regression lines for the different quantities *g*_{e}, *g*_{i}, *V*_{m}, and *f* (σ = 0.106, 0.097, 0.115, and 0.126, respectively; cf. Fig. 1B). The posterior distribution of *sl* and *ic* is then calculated according to Bayes’ rule, yielding the Bayesian posterior (BP)

In the following, we assume a noninformative (flat) prior *p*(*sl*, *ic*) = constant.

The BP(*sl*) for different slopes *sl* irrespective of the intercept is obtained by marginalizing the BP over all intercepts *ic*:

The posterior is then normalized to its maximum, that is, to its value for *sl* = *sl*_{data}, yielding the normalized BP

Most of our analysis is restricted to NBP(*sl*) because our main interest is the variation of the tuning properties with the map location (captured by the slope of the regression line) and not factors affecting the absolute tuning (captured by the intercept of the regression line), which would, for example, be affected by shifts in the baseline response. By multiplying the individual normalized BPs for *V*_{m}, *g*_{e}, *g*_{i}, and *f*, we can combine the evidence from all measurements.

In order to find this combined BPs of the OSI values of *V*_{m}, *g*_{e}, *g*_{i}, and *f* for the populations of pinwheel and orientation domain cells (cf. Fig. 6C), we first calculated the average model prediction for the OSI of *V*_{m}, *g*_{e}, *g*_{i}, and *f* by pooling only cells in pinwheel areas (map OSI < 0.3) or orientation domains (0.6 < map OSI < 0.9); larger values were not selected as these extreme local map OSIs occur in the corners of the orientation map and are thus not available for all orientations. We then repeated the analysis for the experimental data and assessed the likelihood of the experimental data points under the assumption of Gaussian additive noise in each group, where the SD was again estimated from the measured data. The normalized BP was then calculated as above by assuming again a noninformative (flat) prior.

### Significance Testing for Reverse Correlation Simulation

To assess the significance of the difference in the variance between the responses of cells located close to pinwheel centers (input OSI < 0.3) and in orientation domains (0.6 < input OSI < 0.9), we employed a pointwise bootstrap method: 30 pinwheel and 30 orientation domain cells were drawn at random from the set of all model cells 1000 times. For every draw, mean and variance of the temporal response was calculated separately for pinwheel and orientation domain cells. For each time lag τ, it was then counted in which proportion of draws the variance of the pinwheel cell responses was larger than that for the orientation domain cells. If this proportion was above 95%, we regarded the difference as significant. For the experimental data, we applied the same method, sampling a subset of 15 pinwheel and 15 orientation domain neurons at each repetition.

### The Firing Rate Model

The firing rate model (cf. Kang et al. 2003) consists of 3 populations of threshold-linear neurons, each arranged in a 2-dimensional grid of 64 × 64 cells. The populations represent fast excitatory (τ_{E1} = 5 ms), slow excitatory (τ_{E2} = 50 ms), and fast inhibitory (τ_{I} = 5 ms) neurons, where τ* _{j}* is the synaptic conductance time constant of the population

*j*. Lateral connection strengths between the neurons are weighted according to a Gaussian distribution with the same spatial extent (σ = 125 μm) for excitation and inhibition (see the calibration of the lateral connection extent in the Supplementary Information); periodic boundary conditions were used. Recurrent excitatory input is contributed as a weighted mixture by the fast (40%) and the slow excitatory population (60%). In addition, both excitatory and inhibitory cells receive identically tuned feed-forward input, determined by an artificial orientation preference map consisting of 4 pinwheels (Fig. 2). All firing rates and excitatory and inhibitory inputs are measured after the network settles in a steady state for a given stimulus. A complete description of the firing rate model can be found in the Supplementary Information.

### The Hodgkin–Huxley Network Model

The detailed model is similar to Mariño et al. (2005) and consists of Hodgkin–Huxley type point neurons (Destexhe and Paré 1999; Destexhe et al. 2001). Synaptic conductances were modeled as originating from γ-aminobutyric acid_{A} (GABA_{A})-, α-amino-3-hydroxyl-5-methyl-4-isoxazole-propionate (AMPA)-, and N-methyl-D-aspartic acid (NMDA)-like receptors (Destexhe et al. 1998); additional conductances represent background activity (Ornstein–Uhlenbeck conductance noise). Orientation preferences were assigned according to the location in the calibrated artificial orientation map (Fig. 2). The network was composed of 50 × 50 excitatory and 1/3 × (50 × 50) inhibitory neurons and corresponds to a patch of cortex 1.56 × 1.56 mm^{2} in size. In order to avoid boundary effects, we used periodic boundary conditions. As in the firing rate model, spatially isotropic synaptic connections with a calibrated radial profile (σ_{E} = σ_{I} = 125 μm) were used both for excitatory and inhibitory cells. Afferent inputs to excitatory and inhibitory cortical cells were moderately tuned (circular Gaussian tuning function with σ_{Aff} = 27.5°) and modeled as Poisson spike trains. In order to calculate the OSI–OSI relationship for *V*_{m}, *g*_{e}, *g*_{i}, and *f*, an input spike train with a constant rate was applied and the network was simulated for 1.5 s with 0.25 ms resolution (usually, the network settled into a steady state after a few hundred milliseconds). For analyzing the tuning properties, we calculated the firing rate, the average membrane potential, and the average total conductances for every cell, in an interval between 0.5 and 1.5 s.

For the reverse correlation simulations, the stimulus was modeled as a time series consisting of 20 ms frames of 1 of 16 orientations or a “blank” frame. For each neuron, the input rate was again calculated according to a circular Gaussian orientation tuning curve (σ_{Aff} = 27.5°) with the peak at the preferred orientation of the neuron, but the time series was subsequently convolved with different temporal filters, before Poisson-distributed spikes were generated. Each model neuron received 1 out of 4 different temporal kernels, in order to capture the variability present in LGN and V1 simple cell responses (Alonso et al. 2001). The network was simulated for at least 100 s with 0.25 ms resolution. The spike output of the excitatory cells was then analyzed according to the reverse correlation paradigm as described in previous studies (Ringach et al. 1997; Schummers et al. 2007). All analyses of the reverse correlation kernels calculated were restricted to the period from 0 to 100 ms. Later responses cannot fit the data from the reverse correlation experiments because the model does not simulate the effect of long-range connections (omitted for simplicity) on the response time course. A detailed description of the Hodgkin–Huxley network model can be found in the Supplementary Information.

## Results

In order to constrain the cortical operating point, we evaluate experimental data from 2 studies on orientation selectivity in cat V1 (Mariño et al. 2005; Schummers et al. 2007): The first data set (Mariño et al. 2005) consists of intracellular measurements quantifying the tuning of the membrane potential (*V*_{m}), the spike response (*f*), the total excitatory (*g*_{e}), and the total inhibitory (*g*_{i}) input conductance as a function of the position in the orientation map in response to full-field drifting grating stimulation (see Fig. 1B). The tuning of each of these properties was quantified for each cell using the OSI and so was the distribution of orientation selective cells in the neighborhoods of the neurons. This allows assessing the tuning of the response properties depending on the position of the neurons in the orientation map. The OSI of the membrane potential as well as the OSI of the total excitatory and inhibitory conductances vary strongly with map location, whereas the OSI of the firing rate does not.

The second data set (Schummers et al. 2007) consists of extracellular measurements of V1 neuronal responses in a reverse correlation paradigm (cf. Ringach et al. 1997). Here, a comparison between the response time course at preferred orientation for cells close to a pinwheel center and cells in the orientation domain shows that the mean response averaged across all neurons in both populations is similar (Fig. 1C1), whereas the corresponding variance across neurons is higher for the population close to pinwheel centers (Fig. 1C2).

In a first set of analyses, we assessed orientation tuning in 2 network models of V1 and quantitatively compared the model responses with the intracellular measurements obtained with the drifting grating stimulation. For every parameterization of our model, OSI–OSI plots similar to those of Figure 1B can be generated from the simulation data, and the best linear fit is determined by linear regression. Comparison with the experimental data then allows a quantitative evaluation of the confidence that a given model is correct given the experimental evidence using Bayesian analysis. The firing rate model is governed by fewer free parameters and is less prone to overfitting. The Hodgkin–Huxley network on the other hand allows taking into account the full set of experimental data including the membrane potential tuning, while ensuring through the comparison with the simpler firing rate model that the results are plausible.

Beyond considering the steady state, the Hodgkin–Huxley network also allows investigating the temporal dynamics, which we do in a second set of analyses. For every parameterization of the model, we obtained temporal response kernels using the reverse correlation technique. The comparison of the mean response kernels and their variance across neurons in the pinwheel and the orientation domain regions to the experimental data (cf. Fig. 1C) imposes further constraints on the potential operating regime.

### Dependence of Orientation Tuning on Map Location in a Firing Rate Model

We set up a firing rate (mean field) model of a 2-dimensional cortical orientation map that consists of 4 pinwheels (see Fig. 2). All excitatory and inhibitory neurons in the model receive orientationally tuned afferent inputs with similar tuning widths and with their preferred orientations being assigned according to the orientation map. We then systematically varied the recurrent excitatory synaptic strength to excitatory (*S*_{EE}) and inhibitory (*S*_{IE}) postsynaptic cells while keeping the recurrent inhibitory synaptic strengths (inhibition to excitation, *S*_{EI}, and inhibition to inhibition, *S*_{II}) fixed. For every parameter combination, we determined the firing rate (*f*), the total excitatory (*g*_{e}), and the total inhibitory (*g*_{i}) inputs as a function of stimulus orientation. Excitatory and inhibitory inputs directly correspond to mean conductances in conductance-based models (see Supplementary Information). Therefore, the input tuning of model cells will be compared with the conductance tuning of the recorded neurons. We then determine the best-fitting regression lines for the OSIs of *f*, *g*_{e}, and *g*_{i} as a function of the local map OSI for each model parameterization. The normalized BP for the slopes of the regression lines is used as a measure of how well the model with the particular combination of parameters explains the experimental data.

Figure 3A shows the BP (gray value; see scale bar) as a function of the synaptic weights of recurrent excitation (*S*_{EE}) and inhibition (*S*_{EI} × *S*_{IE}). Following Kang et al. (2003), we distinguish 4 parameter regimes based on the properties of the intracortical feedback kernel: “FF,” feed-forward, where the afferent input dominates; “EXC” (corresponding to regime I in Kang et al. 2003), excitatory dominated; “INH” (corresponding to regime IV in Kang et al. 2003), inhibitory dominated; “REC,” recurrent (corresponding to regimes II and III in Kang et al. 2003), which is characterized by strong recurrent excitation and inhibition. The recurrent regime is located along the border to a parameter regime where the network is either “unstable,” that is, where the output rates diverge or where it enters the so-called marginal phase (“MP”), whose key feature is the sharpening of a broadly tuned feed-forward input (Ben-Yishai et al. 1995; Hansel and Sompolinsky 1996). In the MP, there exist certain combinations of model parameters for which the preferred orientation of the output of some model neurons is different from the preferred orientation of their feed-forward drive. This was not observed for combinations of parameters in the other regimes.

Each of the parameter regimes displays a characteristic relation between the OSI of the output rate, the excitatory and the inhibitory conductance, and the tuning of the local input area (map OSI). Typical examples of these relationships are shown in Figure 3B. For small values of the recurrent synaptic strength, the excitatory cells are mainly driven by the feed-forward input that leads to a map invariant spike tuning and only a weak dependence of the total excitatory input tuning on the position in the orientation map (Fig. 3B, FF). The tuning of the total inhibitory input (*g*_{i}) depends more strongly on map location. This is due to the afferent input to the inhibitory population. Even for low values of *S*_{IE}, this population therefore inhibits the excitatory population via the connection *S*_{EI}, whose strength was kept fixed. Further increasing *S*_{IE} while keeping *S*_{EE} small (Fig. 3B, INH) leads to a sharpening of the firing rate tuning via the iceberg effect. Because the effect is stronger close to pinwheel centers, it leads to a sharper firing rate tuning when compared with orientation domains. In this regime, the tuning of the total excitatory input remains invariant across the map and the slope of *g*_{e} is therefore close to zero. When *S*_{EE} is increased for small *S*_{IE} (Fig. 3B, EXC), on the other hand, the slope of *g*_{e} increases and the firing rate tuning is broadened near pinwheels due to the higher excitatory recurrent input at orthogonal orientations. Only when both the excitatory and the inhibitory strength are large (Fig. 3B, REC), we obtain slopes for the total excitatory and inhibitory input tuning that are as steep as in the measured data, whereas the spike rate tuning remains approximately independent of the map location. For this reason, the BP is maximal in a regime with considerable contribution of the excitatory and inhibitory recurrent inputs. The OSI–OSI relationship in the MP (Fig. 3B, MP) is similar to the excitatory dominated regime. It shows a strong dependence of tuning selectivity on map location for total input and firing rate.

Thus, the Bayesian analysis indicates that of the parameterizations tested, neither purely feed-forward regimes nor regimes dominated by either inhibitory or excitatory connections can satisfactorily account for the relationship between the orientation tuning of neurons and their location in the orientation map—quantified through their OSI–OSI relationships. Models with significant recurrent excitation and inhibition on the other hand are able to reproduce the experimentally observed relationship reasonably well.

### Dependence of Orientation Tuning on Map Location in a Hodgkin–Huxley Network Model

In order to validate the above finding in a more biologically plausible network, we set up a network of Hodgkin–Huxley point neurons with conductance-based synapses, allowing us to additionally evaluate the tuning properties of the membrane potential. We used the same 4 pinwheel map as for the firing rate model. Every excitatory and inhibitory neuron received well-tuned time-invariant afferent input, with a reduced afferent input strength for inhibitory neurons. Again, the synaptic strengths (maximum conductance values) of excitatory connections to excitatory ($g-EE$) and inhibitory cells ($g-IE$) were varied systematically, keeping all other parameters constant. For each such model parameterization, we determined the best-fitting regression lines for the OSI of *V*_{m}, *f*, *g*_{e}, and *g*_{i} as a function of the local map OSI. Every OSI–OSI plot is then quantitatively compared with the experimental data using the normalized BP for the slopes of the regression lines. Figure 4 shows the posteriors for the slopes of *V*_{m}, *f*, *g*_{e}, and *g*_{i}. Where the average firing rate exceeds 100 Hz (above the thick solid line), we refer to the regime as unstable because the network shows self-sustained activity, that is, the network activity remains at high firing rates if the afferent input is turned off. We do not evaluate the posteriors in this region. For parameter combinations approaching this “regime of instability” (above the thin solid line), the model neurons’ responses increasingly lose their orientation tuning (no model neuron has a firing rate tuning with an OSI above 0.3). This untuned network response is reflected in a drastically decreased location dependence of the conductance tuning (see contour lines for *g*_{e} and *g*_{i} in Fig. 4). The average firing rate in this regime is very high (>50 Hz) and increases for higher values of $g-EE$. As can be seen in the figure, below the solid lines, the normalized posteriors of *V*_{m}, *g*_{e}, and *g*_{i} support values for $g-EE$ and $g-IE$ that are close to the untuned and the unstable regime. The normalized posterior for *f* is consistent with such parameter values, too. However, it is much less informative than the normalized posteriors for *g*_{e}, *g*_{i}, and *V*_{m} as many other parameter combinations of $g-EE$ and $g-IE$ lead to a high posterior, as well. The very weak location dependence of spike tuning—as quantified via the OSI–OSI relationship—is consistent with the data for a large range of model parameters.

To assess their combined evidence, we calculate the product of the individual posterior values from Figure 4, as a function of $g-EE$ and $g-IE$ (Fig. 5A), and show the OSI–OSI relationship for 4 representative points (Fig. 5B). The OSI–OSI relationships of *f*, *g*_{e}, and *g*_{i} are qualitatively similar to the corresponding plots of the mean field model (cf. Fig. 3B) and support the hypothesis that the “most likely” operating regime is recurrent. The OSI–OSI plot of *V*_{m} is flat in the feed-forward regime (Fig. 5B, FF), shows a negative slope in the inhibitory (Fig. 5B, INH) and a positive slope—similar to what is observed experimentally—in the excitatory dominated (Fig. 5B, EXC) and recurrent (Fig. 5B, REC) regimes. The recurrent regime (Fig. 5B, REC) is close to the instability line. In fact, increasing $g-EE$ by just 10% is sufficient for the network model to show self-sustained activity with firing rates of more than 150 Hz.

Figure 6A shows the effective recurrent input currents into excitatory (Fig. 6A1) and inhibitory (Fig. 6A2) neurons. The effective input current is given by the current through recurrent synapses, normalized by the current through afferent synapses. In the most likely operating regime, the recurrent excitatory current dominates the afferent current by a factor of 1.7. The current through inhibitory synapses exceeds the afferent current by a factor of 1.4, suggesting that excitation and inhibition are relatively balanced. For the most likely regime, we also observe a precise covariation of the tuning of the excitatory and inhibitory conductances, that is, the ratio of their slopes of the corresponding OSI–OSI plots are close to 1 in the area of the highest posterior (see Fig. 6B).

To this point, all models assumed uniform values for $g-EE$ and $g-IE$ across the entire network. However, the strength of the recurrent connections may well differ between pinwheel regions and orientation domains. To assess whether the data provide evidence for this hypothesis, we average the OSI of *V*_{m}, *f*, *g*_{e}, and *g*_{i} of all model cells in pinwheel regions (map OSI < 0.3) and in orientation domains (0.6 < map OSI < 0.9). We then compare the average OSIs in pinwheel and orientation domain cells separately to the OSIs of the experimental measurements from cells in the corresponding map regions and quantify the deviation between model prediction and data by a BP, whose value is high if prediction and data match well and low otherwise (see Methods). The regimes of highest posterior for the pinwheel and for the orientation domain data overlap strongly (Fig. 6C); hence, there is no evidence for different values of the recurrent synaptic strengths $g-EE$ and $g-IE$ in pinwheel areas and orientation domains. The regions also overlap with the values of highest posterior for the slopes of the OSI–OSI relationship shown in Figure 5A.

In summary, we find that the analysis of the firing rate model and the Hodgkin–Huxley network model leads to similar conclusions with regard to the most likely operating regime given the experimental data. The most likely regime is one with a significant contribution of the recurrent excitatory and inhibitory connections, and it is located close to a regime of instability in parameter space.

### Influence of the Spatial Extent and the Strength of the Lateral Inhibitory Connections

Up to now, we assumed equal spatial extent of lateral excitatory and inhibitory connections in line with Mariño et al. (2005). Earlier studies, however, suggested that inhibitory connections are more (Fitzpatrick et al. 1985; Callaway 1998) or less (Buzás et al. 2001) spatially restricted than excitatory ones. Does the most likely operating regime of the Hodgkin–Huxley network model critically depend on the spatial extent of the lateral inhibitory connections? We first reduced the range of synaptic connections from inhibitory cells to 50% of the range of the excitatory cells, while leaving all other parameters fixed. Figure 7A1 shows the product of the individual posteriors, summarizing the comparison of OSI–OSI relationships between the model and the experimental data. The area of high posterior remains similar in extent and in its location relative to the instability line compared with the case of equal spatial extent of excitatory and inhibitory connections shown in Figure 5A. If we increase the range of inhibitory synaptic connections to 150% of the original value, the area of high posterior again remains similar in extent and in location relative to the instability line (Fig. 7A2). Thus, the experimental data support a recurrent operating regime, independent of assumptions about the relative spatial range of lateral inhibitory and excitatory interactions.

Other parameters that were fixed in all previous simulations, but might crucially influence the most likely operating regime, are the strengths (peak conductance values) of recurrent inhibitory synapses. Therefore, we again varied $g-EE$ and $g-IE$ but now set $g-II$to 50% (Fig. 7B1) or 150% (Fig. 7B2) of its standard value, without changing any of the other parameters. When $g-II$ is decreased, the network remains stable for higher values of $g-EE$ (Fig. 7B1). The reason is that inhibitory cells increase their firing rate which leads to strong inhibitory input to excitatory neurons. Accordingly, when $g-II$ is increased, the network becomes unstable for lower values of $g-EE$ (Fig. 7B2). However, in both cases, the area with high posterior values shifts accordingly and thus remains in a recurrent regime, close to the “line of instability.” Finally, we set $g-EI$ to 50% (Fig. 7C1) or 150% (Fig. 7C2). This again leads to a shift of the stability line, but the area of high posterior values shifts accordingly.

In summary, Figure 7A–C shows that the main conclusion drawn from the data, namely that the cortical operating point is characterized by strong recurrent inputs, is robust against changes in several basic model assumptions. Each of the models can account for the location dependence of orientation tuning in the experimental data: The maximum posterior values in each parameterization of the model are of similar magnitude in Figure 7A–C (between 0.22 and 0.41). This implies, on the other hand, that the physiological data considered here do not impose any hard constraints on model parameters such as the range of excitatory versus inhibitory interactions. Additionally, although the location dependence of orientation tuning restricts the most likely operating regime to high values of recurrent excitation and inhibition, it does not allow determining the absolute values of the strength of the afferent and the recurrent synapses.

### Dependence of the Orientation Tuning Dynamics on Map Location in a Hodgkin–Huxley Network Model

To this point, we used a time-invariant input and considered the steady-state behavior only, abstracting from the dynamics of the V1 responses. Including temporal filters, which describe the dynamics of the afferent input to each model cell, in the Hodgkin–Huxley network model, we simulated reverse correlation experiments. The stimulus is a sequence of gratings with randomly chosen orientations, and a new grating is presented every 20 ms. Single-unit recordings from cat primary visual cortex under this paradigm (Schummers et al. 2007) showed that neurons close to pinwheel centers and neurons in orientation domains exhibit a similar time course in their averaged responses but differences in their intercell variability. The mean responses of orientation domain cells are more similar to one another than those of pinwheel cells (cf. Fig. 1C). We have shown previously (Schummers et al. 2007) that these characteristics of the temporal dynamics of V1 can be reproduced in a Hodgkin–Huxley network model and that the differential variability can be attributed to the variability present in the afferent input provided by neurons with different temporal characteristics. The recurrent connectivity removes some of this variability, but the degree of “smoothing” differs between orientation domains and pinwheel centers, causing the observed differences.

Here, we investigate how the temporal characteristics in pinwheel and orientation domain neurons vary with different parameterization of the Hodgkin–Huxley network. The temporal responses of 4 different models (one for each previously identified operating regime; cf. dots in Fig. 8A) to reverse correlation stimulation are shown in Figure 8B, using the same format as for the experimental results in Figure 1C1. Note that the timescale is shifted for the simulated data: here, zero corresponds to the input onset for V1, and in the experiment, zero refers to stimulus onset. Because the network model does not include long-range connectivity, it cannot account for the late part of the response. Therefore, all analyses presented here are restricted to the period from 0 to 100 ms (simulation timescale).

We find that in an excitation-dominated regime (Fig. 8B, EXC), the responses of orientation domain cells—for their preferred stimulus and averaged over all cells—are longer than those of cells close to pinwheel centers. This is intuitive: For an orientation domain cell, many neighboring cells have similar preferred orientations (cf. Fig. 1A2) and, therefore, the excitation-dominated recurrent input will amplify the response to the cell's preferred orientation most strongly. A pinwheel center cell, on the other hand, receives recurrent input from cells with a broader range of preferred orientations and, therefore, its recurrent input in response to the preferred orientation will not affect the cells response as much. It is thus the stronger recurrent drive in orientation domains that prolongs the response to a phasic input compared with pinwheel centers, where recurrent input is weak. However, such a prolonged response for cells in the orientation domain is inconsistent with the observations in vivo (Fig. 1C1). For the other example regimes, on the other hand, the phasic responses of orientation domain cells are similar (Fig. 8B, FF and REC) or even slightly shorter than those of the pinwheel cells (INH).

The difference between the response curves of pinwheel and orientation domain cells is quantified in Figure 9A. The figure shows the mean difference between the averaged response curves of pinwheel and orientation domain cells in the interval 0–100 ms (i.e., a measure for the size of the shaded areas in Fig. 8B) as a function of $g-IE$ and $g-EE$. Averaged responses are relatively similar, that is, differences between response curves are small to moderate (darker colors) for a relatively wide range of feed-forward, recurrent, and moderately inhibition-dominated regimes. The values are compatible with what has been observed in vivo (0.06 for the interval 50–150 ms). For the example models shown in Figure 8A (dots), we obtain values of 0.02 (FF) and 0.06 (REC and INH), respectively. Strongly inhibitory as well as excitatory regimes show markedly larger differences (up to 0.15 for the more extreme cases and 0.08 for the EXC regime shown).

The difference in the variance of the temporal responses between pinwheel and orientation domain cells observed in vivo provides a second important constraint. Figure 8C shows the response variance as a function of time for the selected example models. Only the excitatory dominated (Fig. 8C, EXC) and the recurrent regimes (Fig. 8C, REC) show the observed higher response variance in cells close to pinwheel centers. To quantify the response variance as a function of $g-IE$ and $g-EE$, we compute, for every combination of parameters, the proportion of the interval 0–100 ms, for which the variance between pinwheel cells is significantly higher than the variance between orientation domain cells (see Methods). Figure 9B shows that this fraction is highest along the line of instability (EXC and REC regimes), where values similar to the experimental data are achieved (Fig. 1C2). Experimental data show significantly increased variance for 49% of the interval from 50 to 150 ms, though the variance difference is also partly in the late response, which the model cannot reproduce. Nonetheless, the value compares well with the 30% and more for the excitatory dominated and recurrent regimes (47% and 30%, respectively, for the 2 examples shown) and is incompatible with the 0% found for both the feed-forward and inhibitory dominated regimes.

In order to understand the reason for the dependence of the variance differences on $g-IE$and $g-EE$, one needs to consider the smoothing influence of the recurrent input. In Schummers et al. (2007), we demonstrated that in order to reproduce the differential response variability, the model requires that afferent input is provided by neurons with different temporal characteristics. As has been argued above, by local averaging through recurrent excitatory connections, variability is reduced. In the highly recurrent regime close to the line of instability, the recurrent excitatory input into orientation domain cells can dominate their afferent input, whereas for pinwheel center neurons, the relative contribution of the recurrent input is much weaker, leading to a particularly large difference of variability in this regime. This is confirmed by Figure 9C, which compares this relative strength of the afferent drive around pinwheel centers with that of orientation domain cells. As can be seen, only in regimes with a sufficiently strong excitatory drive, that is, the regimes close to the line of instability, is the relative influence of the afferent drive stronger around pinwheel centers.

In sum, the temporal dynamics of V1 responses and their dependence on location in the orientation map also constrain the likely operating regime: The response shape as well as the similarity in the mean response are not compatible with either excitatory dominated or strongly inhibitory regimes; the response variability found in vivo can, on the other hand, only be reproduced in regimes close to the line of instability, that is, excitatory dominated or balanced recurrent regimes. Taken together, only the strongly recurrent regime with balanced excitation and inhibition can account for both aspects of the in vivo data. This is fully consistent with the model-based analysis of stationary responses, that is, the previously considered orientation tuning properties of the membrane potential and the input conductances of V1 neurons and their dependence on location in the orientation map. Both results constitute converging evidence for a cortical network operating in a regime close to the “border of instability,” beyond which the network settles into a state of high, self-sustained activity.

## Discussion

In this paper, we assessed which cortical operating regime is best suited to explain the dependence of orientation tuning properties on cortical map location as observed in recent experiments (Mariño et al. 2005; Schummers et al. 2007). Intracellular measurements (Mariño et al. 2005) quantify for the same neurons: 1) input (the total excitatory and the total inhibitory conductance), 2) state (subthreshold tuning of the membrane potential), and 3) output (spike response) as a function of local context (position in the orientation map). The second data set (Schummers et al. 2007) complements this by quantifying the time course of the responses. We selected these particular data sets because access to both sub- and superthreshold tunings, recorded for the same neurons at specific map locations, is key to constraining the model. Spike output alone, for example, is compatible with a wide range of model parameters (cf. Fig. 4B).

We systematically varied the strength of recurrent excitation and inhibition in a firing rate network model and in a biologically more plausible Hodgkin–Huxley network model. We investigated which combination of strengths of recurrent excitation and inhibition could account best for the location dependence of tuning properties and for the measured response dynamics. For both types of models and both data sets, our results indicate that significant recurrent contributions are important for reproducing the experimental results. These contributions need to consist of excitatory and inhibitory components balancing one another and dominating the afferent input. Due to the strong recurrent drive, the most likely regime is located close to a region in parameter space where the network settles into a state of strong, self-sustained activity.

### Robustness of the Simulation Results

Although not supported by our recent data (Mariño et al. 2005), several other modeling studies used a different spatial range for lateral excitatory and inhibitory connections (Ben-Yishai et al. 1995; Somers et al. 1995; McLaughlin et al. 2000; Kang et al. 2003). We therefore analyzed the response of the Hodgkin–Huxley network model for different corticocortical synaptic radii. Specifically, we varied the spatial range of the lateral inhibitory connections keeping the range of the lateral excitatory connections constant. Varying the ratio of the range by ±50% did not affect the position of the high posterior region in parameter space being in the recurrent regime close to the instability line (cf. Fig. 7A).

Furthermore, we varied the strengths of recurrent inhibitory synapses, which led to a change in the range of parameters, for which the network is stable (Fig. 7B,C), that is, does not show strong, self-sustained activity. However, the area with high posterior values always changes accordingly and remains in a recurrent regime, close to the line of instability. Thus, our main conclusion concerning the most likely operating regime remains valid under such changes.

In addition, we assumed that the local connectivity pattern is independent of the location in the orientation preference map, as indicated by experimental data (Mariño et al. 2005). Also the peak conductances for recurrent synapses in the model were identical for neurons close to pinwheel centers and in orientation domains. This seems reasonable because we found no indication that the optimal parameters for the peak conductances are different in orientation domain and pinwheel center neurons (cf. Fig. 6C).

An important parameter that is not well characterized experimentally is the ratio of AMPA versus NMDA receptors of excitatory synapses in visual cortices (Myme et al. 2003). In the Hodgkin–Huxley network model, 70% of the synapses were fast (AMPA-like) and 30% were slow (NMDA-like). If the fraction of slow synapses is increased, the network remains stable for higher values of recurrent excitation but the position of the region of high posterior shifts too and remains close to the instability line in a dominantly recurrent regime (data not shown). However, the number of fast versus slow synapses is crucial for obtaining a higher response variance in pinwheel regions compared with orientation domains. We tested several ratios of fast versus slow synapses and found that more than about 50% of the synapses have to be fast in order to obtain the experimentally observed higher variance close to pinwheel centers.

Some very recent results (Nauhaus et al. 2008) question the observation that the firing rate tuning is entirely location independent, arguing that it, too, should be correlated with the tuning of the local recurrent inputs of a cell. Such correlations can be observed in the models closer to the instability line: Lowering the value of the inhibitory peak conductance from the most likely operating point by just 10% increases the slope of the location dependence of the firing rate tuning to 0.25 without affecting the tuning of the other properties strongly (see contour lines in Fig. 4). The product of the posteriors for such a dependence of firing rate tuning on the local input OSI (having all other OSI–OSI relationships fixed) would thus shift slightly but be nonetheless in a strongly recurrent regime close to the line of instability as the most likely regime is strongly determined by the location dependence of *V*_{m} and *g*_{e}.

The conclusion that significant recurrent excitation and inhibition are necessary to account for the experimentally observed tuning properties was drawn from computational models of different complexity and 2 independently recorded data sets. With the comparatively simple firing rate model, we demonstrated that a strongly recurrent operating regime accounts best for the experimental results. The more complex Hodgkin–Huxley network model provided confirmation of this result, taking also into account the tuning of the membrane potential. Adding realistic temporal dynamics to the afferent input of the Hodgkin–Huxley network model, we found that also the map location dependence of the temporal response properties of a second set of recordings from cat V1 neurons could only be accounted for in a regime with significant recurrent contributions. Thus, similar results were obtained with 3 independent approaches of differing complexity, providing an additional evidence for the robustness of the results.

### Biological Plausibility of the Models

In our models, we assumed a relatively sharply tuned feed-forward excitatory drive to all model cells. Earlier studies support the assumption that the afferent input to layer 4 simple cells in cat V1 is sharp and well selective for orientation (Reid and Alonso 1995; Alonso et al. 2001). Furthermore, inactivating the recurrent input into visual cortex by cooling or electrical shocks does not remove the selectivity of the cell responses (Ferster et al. 1996; Chung and Ferster 1998). Looking at the response time course, both intracellular data (Gillespie et al. 2001) and optical imaging (Sharon and Grinvald 2002) show that the orientation selectivity is well-tuned right from stimulus onset and does not sharpen over time in cat V1. Moreover, also neurons in superficial layers receive their strongest drive from well-tuned layer 4 simple cells (Alonso and Martinez 1998; Martinez and Alonso 2001). Thus, our assumption of a fairly selective feed-forward excitatory drive is likely to be valid also for neurons in layer 2/3.

In our model networks, the orientation of the peak of the conductance tuning for both *g*_{e} and *g*_{i} is always at the preferred orientation as determined by the firing rate of the cell, consistent with an earlier study (Anderson et al. 2000) and our own observations (Mariño et al. 2005). In a detailed recent intracellular work, inhibition was found to be preferentially tuned to the cross-orientation in 40% of all cells (Monier et al. 2003). However, this variation can possibly be explained by the laminar position of the cell in the network: Whereas cells in layer 4 and layer 2/3 show matching preferred orientation for their input conductance and their firing rate, layer 5 cells characteristically display a net hyperpolarization that is stronger at the orientation orthogonal to the preferred orientation of the firing rate (Martinez et al. 2002; Hirsch and Martinez 2006).

Other recent physiological studies emphasized the role of both a tuned and an untuned inhibitory components in the generation of orientation selectivity (Ringach et al. 2003; Shapley et al. 2003) that could be realized in layer 4 by 2 functionally different groups of inhibitory interneurons (e.g., Hirsch et al. 2003; Nowak et al. 2008). However, another recent study did not find evidence for untuned inhibitory neurons in any layer of V1 (Cardin et al. 2007). Our model networks have only a tuned inhibitory component: The firing rate of inhibitory cells displayed a tuning similar to excitatory cells. However, we find only a small parameter regime where the slope for *g*_{i} is steeper than what we observed in the data. An additional untuned inhibitory baseline would reduce the slope of *g*_{i} significantly, and model predictions would no longer match the data well. Therefore, the observed map dependence of *g*_{i} tuning does not allow for dominant untuned recurrent inhibition.

A characteristic of the physiological data and, therefore, also of the most likely regime in our model networks is the strong covariation of the tuning of the excitatory and the inhibitory conductance with the local input region. This covariation indicates a strong overall contribution of the local cortical synaptic network to the input conductance (Mariño et al. 2005). In line with these findings, other recent studies stressed the importance of well-balanced recurrent excitation and inhibition, for example, as a mechanism that determines the spike firing regularity in ventral cochlear nucleus chopper neurons of the rat (Paolini et al. 2005) or as a way to generate multiple states of activity in local cortical circuits (Shu et al. 2003).

### Computational Contributions of the Cortical Network

In the most likely operating regime, intracortical processing is not necessary to establish sharp orientation selectivity. In the Hodgkin–Huxley network, the recurrent processing sharpens a moderately tuned afferent input (circular Gaussian with *σ*_{Aff} = 27.5°) only a little, to on average σ_{Exc} = 23.3° (excitatory cells). Nevertheless, a significant level of local recurrent input is important to amplify the firing rate response. In the recurrent regime (REC in Fig. 5A), neurons responded to their preferred stimulus with a 1.8-fold increased firing rate compared with the feed-forward regime (FF in Fig. 5A). Previous modeling studies have already attributed such an amplification to the cortical network (Douglas et al. 1995; Somers et al. 1995), and so has experimental work (Sharon and Grinvald 2002). Sharon and Grinvald used optical imaging with voltage sensitive dyes to investigate the development of the tuning width and the modulation depth (difference between the responses to the preferred and the orthogonal orientation) over time. Although they found no significant change in tuning width over time, the modulation depth peaked tens of milliseconds after the cortical response began, rather than at stimulus onset. They attributed this peak in the modulation depth to an intracortical synaptic mechanism likely to be of excitatory and inhibitory synaptic origin. It is thus possible that recurrent excitatory amplification plays an important role in the temporal evolution of the orientation selective response. Similar mechanisms have been suggested to operate in auditory cortex (Wehr and Zador 2003).

But why would it be efficient for a cortical network to operate in a regime where strong recurrent activity seems to be “wasted” because a tuned afferent input is mapped onto an almost similarly tuned cortical output? Our analysis demonstrates that due to the strong recurrent contribution, the most likely operating point is close to the line of instability across which the cortical network exhibits self-sustained activity. In this region, small changes in the network parameters lead to relatively large changes in the neuronal activity. An indication for this sensitivity is the change in the mean firing rate and response length when increasing the recurrent excitatory peak conductance $g-EE$. In the regimes close to the line of instability (EXC and REC in Fig. 5A), an increase of 5% leads to an increase in firing rate of 41% (EXC) and 22% (REC) and a prolongation of the response kernel (length at half height) by 10% (EXC) and 20% (REC), respectively. In the other regimes (FF and INH in Fig. 5A), it only changes the firing rate by 3% (FF), respectively 2% (INH), and has a negligible effect on the length of the response kernel. In the identified most likely operating regime, a 10% change in firing rate, which is of similar magnitude as observed firing rate changes under attention in macaque V1 (McAdams and Maunsell 1999), is easily achieved by increasing $g-EE$ by 2%. It is thus conceivable that contextual modulations, adaptation, or attentional effects temporarily shift the operating point, which means that the cortical response is particularly sensitive to small changes in afferent or feedback signals. Experimental data, for example detailing spatial and temporal contextual interactions (Schwartz et al. 2007), may allow to investigate the functional implications of the cortical operating regime. For example, Dragoi et al. (2001) found pronounced adaptation-induced short-term plasticity close to pinwheel centers. It would be interesting to investigate if those findings can be explained by the network model operating in a strongly recurrent, balanced regime.

The peculiar location of the cortical operating point close to the line of instability is reminiscent of the theory of “computing at the edge of chaos” (Langton 1990), which has attracted some attention in the past decade. Theoretical studies have demonstrated that abstract networks can accomplish more challenging computational tasks near the transition from ordered to chaotic dynamics (Bertschinger and Natschlager 2004) and that such a system is responsive to inputs over a wider dynamic range (Kinouchi and Copelli 2006). Such theoretical considerations are complemented by in vivo data showing recurrent excitation and inhibition in balance (Haider et al. 2006), which has also been interpreted as a facilitator for easy transition between different cortical states. Thus, our findings here add to the evidence that the recurrent network contributes to maintaining an operating point that allows easy switching between different cortical states.

## Conclusions

In our study, we use physiological measurements in a quantitative way to constrain the potential operating regime for the computation of orientation in cat primary visual cortex. A range of different model-based data analyses lets us conclude that to represent orientation, visual cortical neurons in cat process a moderately tuned feed-forward input by means of a relatively strong, well-balanced corticocortical recurrency. The tuning of the excitatory conductances provides the most compelling evidence for this conclusion, whereas the spike tuning does not put strong constraints.

The robustness of the results was demonstrated by varying some of the certain model parameters (such as connection range and synaptic strength) as well as by the mere fact that similar conclusions could be drawn from 3 essentially independent modeling approaches. Paraphrasing the robustness to some model parameters, we can also conclude that the available data cannot constrain a few other open issues, such as the range of the inhibitory connections with respect to the excitatory ones.

Speculating about the actual role of the recurrent connectivity, 2 potentially related hypotheses are offered: First, balanced recurrency may serve to amplify a moderately tuned but relatively weak feed-forward input. Second, the corticocortical recurrency establishes a cortical operating point close to a line of instability, across which the network activity increases dramatically. This could make the cortical response particularly sensitive to modulatory effects or feedback signals. An exploration of these hypotheses requires a detailed investigation in a paradigm likely to activate such modulatory connections, for example through long-range interactions or attentional modulations.

Given the similarities of cortical architecture between areas related to different functions, one can expect that these results apply to cortical areas and functions beyond visual cortex and the interpretation of sensory signals. It, therefore, remains an interesting question, whether a balanced recurrent regime is realized in other brain areas.

## Supplementary Material

Supplementary information can be found at: http://www.cercor.oxfordjournals.org/.

## Funding

German Federal Ministry of Education and Research (BMBF) (01GQ0414 to K.O.); and National Institutes of Health (J.M., J.S., D.L., and M.Sur).

We thank Gidi Farhi and Thorsten Dietzsch for helping with the simulations and the analytical treatment of the firing rate model. We thank Oliver Beck and Peter Wiesing for providing the code of the Hodgkin–Huxley network model used in Mariño et al. (2005). *Conflict of Interest*: None declared.