## Abstract

Disparity tuning measured in the primary visual cortex (V1) is described well by the disparity energy model, but not all aspects of disparity tuning are fully explained by the model. Such deviations from the disparity energy model provide us with insight into how network interactions may play a role in disparity processing and help to solve the stereo correspondence problem. Here, we propose a neuronal circuit model with recurrent connections that provides a simple account of the observed deviations. The model is based on recurrent connections inferred from neurophysiological observations on spike timing correlations, and is in good accord with existing data on disparity tuning dynamics. We further performed two additional experiments to test predictions of the model. First, we increased the size of stimuli to drive more neurons and provide a stronger recurrent input. Our model predicted sharper disparity tuning for larger stimuli. Second, we displayed anticorrelated stereograms, where dots of opposite luminance polarity are matched between the left- and right-eye images and result in inverted disparity tuning in the disparity energy model. In this case, our model predicted reduced sharpening and strength of inverted disparity tuning. For both experiments, the dynamics of disparity tuning observed from the neurophysiological recordings in macaque V1 matched model simulation predictions. Overall, the results of this study support the notion that, while the disparity energy model provides a primary account of disparity tuning in V1 neurons, neural disparity processing in V1 neurons is refined by recurrent interactions among elements in the neural circuit.

## Introduction

Most research on the neurophysiology of binocular vision in primary visual cortex (V1) has focused on hypotheses generated by the feedforward disparity energy model (Ohzawa et al., 1990). The disparity energy model, however, is unable to fully explain disparity tuning in V1 (Cumming and Parker, 1997; Samonds et al., 2009; Tanabe et al., 2011), and local solutions can fail to find the correct solution of disparity (Chen and Qian, 2004; Read and Cumming, 2007). More recent theoretical and experimental research suggests models that include neuronal interactions could provide a more accurate description of disparity tuning and such interactions could facilitate disparity processing to find more reliable estimates of disparity (Menz and Freeman, 2003; Chen and Qian, 2004; Read and Cumming, 2007; Samonds et al., 2009; Tanabe et al., 2011).

There are two primary characteristics of disparity tuning in V1 that are inconsistent with the disparity energy model that we will address in this article. First, we previously found that disparity tuning curves evolved over time causing the preferred disparity to be more prominent with respect to nonpreferred disparities (Samonds et al., 2009). The sharpened peak and broadened valleys of disparity tuning curves over time are inconsistent with the Gabor function of disparity predicted by the disparity energy model (Ohzawa et al., 1990). Second, stimulation with anticorrelated stereograms (Julesz and Tyler, 1976) generates inverted disparity tuning curves that have weaker modulation amplitudes than disparity tuning curves that result from standard correlated stereogram (Julesz, 1964) stimulation, although the disparity energy model predicts that the modulation amplitudes should be equal (Cumming and Parker, 1997).

In the present study, we developed a simple single-layer neuronal network model with feedforward inputs based on the disparity energy model (Ohzawa et al., 1990) and recurrent inputs from neighboring neurons constrained by neurophysiological data (Menz and Freeman, 2004; Samonds et al., 2009). We examined whether or not our model could explain the aforementioned variations of V1 disparity tuning from the disparity energy model and performed two experiments to test model predictions. First, the recurrent model replicated sharpening over time, but it also predicted more pronounced sharpening with larger stereograms because more neurons were driven and therefore a larger number of recurrent inputs were activated. When we increased the aperture size of stereograms, progressively sharper tuning was also observed in neurophysiological recordings. Second, the recurrent model produced the observed result that inverted disparity tuning curves from anticorrelated stereogram stimulation had weaker modulation amplitudes than disparity tuning measured from correlated stereogram stimulation. The model additionally predicted much weaker and negligible sharpening of disparity tuning during anticorrelated stereogram stimulation compared with disparity tuning during correlated stereogram stimulation. This was also confirmed with neurophysiological recordings. Both the reduced modulation amplitudes and reduced sharpening occurred in the model because recurrent inputs were weaker when the peaks of anticorrelated disparity tuning curves for recurrently connected neurons did not line up across spatial scale. Overall, these results provide stronger support that facilitative and suppressive interactions among V1 neurons indeed contribute to disparity processing.

## Materials and Methods

##### Neurophysiological recordings.

The data for this study were collected simultaneously with data reported in two previous articles where the details about the specific methods can be found (Samonds et al., 2009, 2012). In brief, two different recording procedures were used on three rhesus monkeys (*Macaca mulatta*) that were approved by the Institutional Animal Care and Use Committee of Carnegie Mellon University and are in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals. The first recording procedure used for two monkeys (male and female) used two to eight tungsten-in-epoxy and tungsten-in-glass microelectrodes in a chamber overlying the operculum of V1 (Samonds et al., 2009). In the second procedure used on the third monkey (male), we recorded from neurons using a chronically implanted 10 × 10 Utah Intracortical Array (400 μm spacing) inserted to a depth of 1 mm in V1 (Samonds et al., 2012). Spike sorting was used to isolate single units (Samonds et al., 2009, 2012).

Dynamic random dot stereograms (DRDS) with 25% density of black (<0.1 cd/m^{2}) and white (50.7 cd/m^{2}) dots on a mean gray background (25.3 cd/m^{2}) and a 12 Hz refresh rate were centered on the mean position of the receptive fields for the population of neurons determined by both minimum response fields based on bar stimuli (Samonds et al., 2009) and spike-triggered receptive fields based on reverse correlation with white noise stimuli (Kelly et al., 2007). Shutter goggles were used to present images to the left and right eyes separately. Because of the small size of the receptive fields (<1 degree) and their tight clustering (highly overlapping), all receptive fields were well within the DRDS stimuli. For the varying aperture experiment, the DRDS was presented in a 2-, 3-, or 4-degree diameter aperture. No dots were presented outside of the aperture and the aperture size was constant across all disparities. For the correlated versus anticorrelated DRDS experiment, the aperture had a 3.5-degree diameter. For correlated DRDS, there was 100% correspondence between black and white dots in the left- and right-eye images (Julesz, 1964). For anticorrelated DRDS, there was 100% correspondence of black dots to white dots, and white dots to black dots, in the left- and right-eye images, respectively (Julesz and Tyler, 1976). Eleven disparities between corresponding dots were tested for both experiments: ±0.94, ±0.658, ±0.282, ±0.188, ±0.094, and 0 degrees. Because the emphasis in this study was quantifying how binocular disparity tuning evolved over time, we only examined the most robust data from our recordings (*n* = 184 neurons). The responses of the neurons had to have highly significant disparity tuning (one-way ANOVA, *p* < 0.01) and a disparity discrimination index (DDI) >0.4 (Prince et al., 2002a; Samonds et al., 2009).

##### Model.

All units in the model were complex cells with feedforward inputs determined by the energy model (Ohzawa et al., 1990; Cumming and DeAngelis, 2001) (Fig. 1, Input). Each cell was given a preferred phase disparity φ_{d}, spatial frequency σ, and receptive field center location *x*_{0}, *y*_{0}. We denote the set of these neuronal tuning parameters by **θ** = {φ_{d}, σ, *x*_{0}, *y*_{0}}. The feedforward input given left and right stimuli *X*_{L} and *X*_{R} was then given by a sum of *N* simple cell responses, each with quadratic nonlinearities, as follows:
where *r*_{S} is the response of a simple cell of phase φ_{i} with left and right receptive fields *R*_{L} and *R*_{R}. *R*_{L} is a Gabor filter centered at *x*_{0}, *y*_{0} with spatial frequency σ and phase φ_{i} − (φ_{d}/2). The summand ranges over two or more quadrature pairs of phase φ_{i}. Thus, the feedforward input into model complex cells follows the energy model (Ohzawa et al., 1990) and was given by the sum of squared binocular Gabor simple cell responses.

When the input stimulus was a DRDS of uniform disparity *d*, the expected value of the feedforward input can be shown to be as follows:
Using Parseval's theorem, we see that the terms under the summand in Equation 3 do not depend on φ_{i} or *N*, and so we can write the following:
This feedforward disparity tuning curve can be shown to be approximately Gabor as a function of *d* (although the model does not make this approximation). The response to anticorrelated DRDS was similar, except the first term of Equation 4 was negative.

In each simulation, we modeled 32 different preferred disparities φ_{d} (even increments from −π to π), eight different preferred spatial frequencies σ (subtending four octaves), and a 21 × 21 grid of spatial locations *x*_{0} and *y*_{0}, for a total of 112,896 simulated neurons. Because preferred disparity is not uniformly distributed in V1, we sampled neurons according to empirically observed distributions (Prince et al., 2002b; Liu et al., 2008; Poole et al., 2010).

The neural response at time *t* is given by *r*(*d, t*), using a standard dynamic neural field model:
where ** r**(

*d, t*) is the vector of all current neural activity levels,

*u*is the membrane potential, and

**(**

*w***θ**) denotes the lateral connections into neuron

**θ**.

The neural static nonlinearity *g* was given by the sigmoid function:
with *u*_{0} chosen to produce a baseline firing rate of 20 spikes per second (sps), and *M* chosen to produce a maximum firing rate of 200 sps. In our model simulations, no neural firing rate exceeded 100 sps. Therefore, the nonlinearity *g* was effectively monotonically increasing and always expansive (∂*g*/∂*u* > 0; Fig. 1, small grid in each neuron box). Other expansive nonlinearities, such as *g*(*u*) = *u ^{2}*, produced similar results.

Synaptic weights between neurons were chosen to match our conclusions from spike time correlation studies (Samonds et al., 2009): facilitative connections were drawn between neurons of similar tuning at the same and nearby spatial locations, and inhibitory connections were drawn between neurons of differing tuning within the same spatial location. Within a spatial location *x*_{0}, *y*_{0}, the weight *W*(**θ**_{i},**θ**_{j}) between neurons *i* and *j* was chosen to be proportional to the Pearson correlation between the feedforward tuning curves of those two neurons (Fig. 1, Local). This choice was designed to approximate Hebbian learning between neurons: connections between neurons strengthen when their firing patterns correlate. Note, however, that the neural response to a natural stimulus is generally substantially reduced and sparser, with fewer neurons firing, in comparison with the response to DRDS. To emulate this, feedforward tuning curves were thresholded before computing neural correlations. Thus, *W*(**θ**_{i},**θ**_{j}) was determined by the correlation between *max*(*T,I*_{F}(*d* | **θ**_{i})) and *max*(*T,I*_{F}(*d* | **θ**_{j})), where *T* was set to the median neural input value. Note that all neurons within a spatial location were interconnected. For example, neurons with differing spatial frequencies σ may have facilitative or inhibitory connections, depending on whether their tuning curves were positively or negatively correlated. The resulting weight distribution was sharp (excitatory connection strength tapered rapidly as two neurons differed in disparity tuning), and was primarily inhibitory (>75% of all lateral connections were inhibitory).

Across spatial locations, neurons were connected only if they had matching disparity and spatial frequency preferences (Fig. 1, Long-range). All cross-spatial connections were positive, and set by a Gaussian, as follows:
with σ_{G} set to 1.0.

##### Quantifying sharpening.

Unlike previous studies that used reverse correlation analysis to measure tuning curves from responses to rapidly changing stimuli (Menz and Freeman, 2003; Chen et al., 2005; Xing et al., 2005), we measured tuning curves at various delays from stimulus onset to stimuli presented continuously for one second. A lot of consideration and testing went into our choice of how to quantify the changes (e.g., sharpening) of disparity tuning over time (Samonds et al., 2013). We define a sharp tuning curve as one with mean firing rates that are very informative about the most likely value of the stimulus. Among tuning curves with fixed mean firing rate and amplitude, a rapidly firing cell with a sharp tuning curve conveys more information and describes the input stimulus more precisely than a rapidly firing cell with a dull tuning curve. We examined fitting a Gabor function to the data, fitting a difference of Gaussians function to the data, calculating the Fourier transform, and calculating sample skewness. All methods including our final selection required that tuning curves were reliably measured and robust over time because we were computing the mean firing rate in 100 ms windows, using 10–60 trials for each disparity. That is why we chose strict selection criteria for the data that we analyzed (see above). This is because flat or noisy tuning curves can lead to fits with outlier parameters or outlier results with all the potential methods.

The primary problem with a Gabor function or a Fourier transform is that both methods would be chosen to characterize sharpening assuming that a single frequency component was increasing over time in the disparity tuning function (Samonds et al., 2013). That assumption was not the case based on our observations of disparity tuning sharpening over time. The primary peak in the disparity tuning function was increasing in frequency, but the valleys were decreasing in frequency (Samonds et al., 2009). Prince et al. (2002a) also previously noted that their Gabor fits deviated from their data with side flanks that were wider than the peak.

Although a Gabor is the standard function for describing disparity tuning over a diverse population of neurons, an alternative function that is not constrained to a single frequency or bandwidth component, and is therefore ideal for capturing the dynamics of disparity tuning described in the study by Samonds et al. (2009), is the difference of Gaussians. A difference of Gaussians function does not always fit well with particular disparity tuning curves and even in a simplified form still requires us to fit 6 parameters to 11 data points, which like a Gabor function, leaves it highly sensitive to the same problems that we encountered with Gabor fits with parameter initialization and outlier results (Samonds et al., 2013). The difference of Gaussians function, however, did provide good fits for some of our robust examples. Although outliers and noise in parameter estimates from difference of Gaussians fits made it difficult to reliably characterize trends over a population of neurons, at least the trends in our robust examples were consistently reflecting some of the properties of sharpening (Samonds et al., 2013).

To simplify our analysis, we chose a method that required no fits, no parameter initialization, and no interpolation. The statistical measurement of the sample skewness is the third standardized moment of a distribution and can be computed directly from the mean firing rates *f*(*d*) for each disparity *d* tested:
Skewness is invariant with respect to the mean and variance of the tuning curve so changes in skewness cannot be attributed to changes in the baseline firing rate or the amplitude of the tuning curve over time. The skewness of the distribution of firing rates over disparity captures all the features of sharpening that we observe with disparity tuning over time (Samonds et al., 2009, 2013). When a small amount of disparities have firing rates far above the mean firing rate across the entire disparity tuning curve, skewness has a high positive value. This happens with narrow positive peaks and broad negative peaks, and skewness will increase if positive peaks become narrower and/or negative peaks become broader. When a small amount of disparities have firing rates far below the mean firing rate across the entire disparity tuning curve, skewness has a high negative value. This happens with narrow negative peaks and broad positive peaks and skewness will decrease if negative peaks become narrower and/or positive peaks become broader. If there are an equal amount of disparities with responses equally above and below the mean (e.g., sinusoidal function), skewness will be equal to zero. Finally, skewness increases if the response to secondary peaks is reduced. Overall, skewness increases for all the features that we observed in disparity tuning over time (Samonds et al., 2009, 2013).

The difference of Gaussians positive Gaussian peak width produced the most consistent results from all the other methods we considered and although both the positive Gaussian peak width and skewness capture the primary characteristic of disparity tuning over time (narrowing peak), skewness is less noisy (especially with lower firing rates and over a larger population of neurons) and outlier results are less extreme because the computation is much simpler and more robust to noisy tuning curves. Additionally, skewness produces a stronger result because it is also capturing the additional characteristics of disparity tuning such as broadening of negative peaks and suppression of secondary positive peaks (Samonds et al., 2013).

All of the methods described above describe the shape of disparity tuning invariant of the strength of the response. However, none of the methods alone can distinguish between what potential mechanisms actually caused the shape of the tuning curves to change over time (e.g., sharpening). In this study, we used data produced from a recurrent neuronal network model that we compared with neural recordings to test whether recurrent connectivity could be the underlying mechanism of sharpening. An alternative source of sharpening could be an expansive output nonlinearity (Ohzawa et al., 1990). An expansive output nonlinearity would sharpen a tuning curve over time if the mean firing rate was increasing. We attempted to avoid this confound by measuring skewness starting at the peak of the population response (100 ms) over an interval where the mean firing rate of the population decreases.

## Results

We developed a neural network model with the straightforward organization of a single layer and recurrent connections among binocular disparity-tuned neurons that represent what has been inferred based on cross-correlation results (Menz and Freeman, 2004; Samonds et al., 2009). The organization includes local positive connections among neurons with similar disparity tuning, different spatial scales, and overlapping receptive fields (Fig. 1; see also Materials and Methods). There are also local negative connections among neurons with different disparity tuning and overlapping receptive fields. And finally, there are distant positive connections among neurons with the same disparity tuning, same spatial scale, and with neighboring receptive fields across the visual field. For inputs (time step = 1), we used tuning curves based on the disparity energy model (Ohzawa et al., 1990). We then compared the dynamics of model disparity tuning curves after applying several iterations of recurrent inputs to the dynamics of disparity tuning curves measured from recordings in the primary visual cortex of awake, behaving macaques while they fixated on DRDS.

Eight spatial frequencies (subtending four octaves) and 32 disparity increments were included in the model. The mean firing rate, disparity units assigned to the model, and the range of disparities and frequencies tested were all matched to what was observed from recordings. Only one model simulation was used, and no additional adjustments were made after running the model. Because the model and recordings provided neurons with a variety of preferred disparities and spatial frequencies, examples were chosen for comparison from the model and recordings that had similar spatial frequencies and preferred disparities.

### Model captures sharpening of disparity tuning

Model tuning curves and tuning curves measured from recordings evolved in a similar manner. Peaks became narrower while valleys became wider (Fig. 2*A*,*B*) as we reported previously from neurophysiological recordings (Samonds et al., 2009). We quantified this behavior with the statistical measurement of the sample skewness of the distribution of mean firing rates over disparity (Eq. 9; see also Materials and Methods) (Samonds et al., 2013). Skewness was measured from recordings using the mean firing rates from 100 ms sliding windows every millisecond. Again, the skewness measured from model data and data from recordings were consistent with each other (Fig. 2*C*,*D*) and consistent with our previous observations (Samonds et al., 2009). The skewness increased more strongly in the earliest iterations and during the earliest portion of the neuronal responses soon after the peak of the response onset. The skewness continued to increase over iterations or time, but at a progressively slower rate. This also happened in the model because the behavior converged to a steady state where the tuning curves no longer changed with a greater number of iterations. Therefore, throughout this article, and as we have done previously (Samonds et al., 2009), we will present all skewness measurements versus log-time-steps (log-iterations) or log-time, and we computed tuning curves from recordings over time using progressively larger windows of time (50, 100, 200, and 400 ms) starting soon after the peak of the response onset in the population (100 ms). For the model, we will compare the input (disparity energy model) and steady-state tuning curves (result of the final iteration). For the recordings, we will compare the tuning curve measured from the initial window of time to subsequent windows of time from response onset. Our initial window is still delayed and therefore would likely include some sharpening from recurrent interactions so it is possible that some changes in actual neurons are happening faster than we can measure them.

Skewness can vary substantially depending on the preferred disparity and spatial frequency of a neuron. For model neurons, the skewness varied from −1.5 to 3.0 for our sample of tuning curves, which is consistent with the variation that we observed in the actual neurons as well. The model does also produce variation in temporal dynamics (some neurons converge to a steady-state faster than other neurons) depending on the preferred disparity and spatial frequency. For actual neurons, the skewness measurements tended to be too noisy in individual cases, making it impractical to systematically compare their skewness convergence with the variations in convergence behaviors in model neurons. The direction of change of skewness over time, however, did not vary in the model. Overall, the examples of different preferred disparities and spatial frequencies that we present throughout this article are representative of the variety of behavior observed in the model and recordings.

### Disparity tuning dynamics for tuned inhibitory neurons

In models, binocular disparity-tuned neurons that are classified as tuned inhibitory (Poggio and Fischer, 1977; Poggio et al., 1988) do not behave in the same manner as tuned excitatory neurons, especially with respect to an expansive output nonlinearity (Read et al., 2002; Haefner and Cumming, 2008). We initially separated our disparity tuned neurons into tuned inhibitory and tuned excitatory categories to test for differences in behavior in the model and neurophysiological recordings. Neurons were classified as tuned inhibitory when they had one primary negative peak and two positive peaks that were not significantly different (Samonds et al., 2012), and our population of 184 disparity-tuned neurons described in this article included 39 (21%) tuned inhibitory neurons. The primary difference we observed in the model, and for the neurophysiological data, is that features of sharpening were inverted for tuned inhibitory neurons with respect to tuned excitatory neurons: the primary negative peak narrowed and became more prominent, and skewness decreased over time.

Figure 2 (top right panels) demonstrates the disparity tuning dynamics for an example tuned inhibitory model and recorded neuron. Over time, the primary negative peaks became narrower and more prominent with respect to other disparities (Figs. 2*E*, black vs gray; *F*, black vs progressively lighter curves). The change in shape over time was confirmed by showing that the skewness was decreasing at an approximately linear rate versus log-time-steps and log-time (Fig. 2*G*,*H*, respectively).

We summarize the inverted behavior of tuned inhibitory neurons (primary negative peak, *n* = 39 neurons) with respect to tuned excitatory neurons (primary positive peak, *n* = 145 neurons) in the bottom row of Figure 2. We sorted tuning curves from the most- to the least-preferred ranked disparity based on the primary positive peak or negative peak, respectively, before averaging. For tuned excitatory neurons, the responses to nonpreferred disparities became relatively weaker over time compared with the response to the preferred disparity in the population average tuning curve (Fig. 2*I*). Additionally, the average skewness increased at an approximately linear rate versus log-time (Fig. 2*J*). For tuned inhibitory neurons, responses to nonpreferred disparities (based on the primary negative peak) became relatively stronger (farther away in mean firing rate difference) over time compared with the response to the preferred disparity in the population average tuning curve (Fig. 2*K*). As in the examples, for tuned inhibitory neurons, the average skewness decreased at an approximately linear rate versus log-time (Fig. 2*L*). Overall, the behavior was very similar between the two subpopulations, but inverted with respect to each other when based on tuning shape, sharpening, and skewness. Results were very similar for alternative criteria in dividing neurons into the tuned excitatory and tuned inhibitory classes, such as the ratio of positive peak height (maximum − mean) compared with negative peak height (mean − minimum). Since the average behavior of tuned inhibitory neurons was consistently inverted for both the aperture size and anticorrelated versus correlated DRDS experiments, we inverted their tuning curves and skewness measurements for any population analysis described in subsequent sections, and we confirm the consistency of the inversion in the final section. Because these neurons only represent 21% of the population, even if we did not invert these tuning curves, the general observation was that during DRDS stimulation with correlated DRDS using a DRDS with a large aperture (>3 degrees), disparity tuning sharpened and skewness versus log-time had a significantly average positive slope (*p* < 0.02).

### Model predictions

To conduct a stricter test of whether or not recurrent interactions like those incorporated into our model could predict the dynamics of disparity tuning observed in our recordings, we examined the dynamics of disparity tuning in the model and from recordings while applying more complex manipulations to DRDS stimuli. First, we examined the disparity tuning dynamics while increasing the size of the DRDS, therefore covering a greater number of receptive fields and exciting a greater number of neurons. Then we compared the disparity tuning dynamics between traditional stereograms (Julesz, 1964) and anticorrelated stereograms (Julesz and Tyler, 1976), where the input tuning curve ends up inverted in the disparity energy model (Cumming and Parker, 1997).

### Disparity tuning dynamics depend on DRDS aperture size

When neurons in primary visual cortex are driven by progressively larger drifting sinusoidal luminance gratings with varying orientation, the orientation tuning curves exhibit progressively more sharpening (Chen et al., 2005; Xing et al., 2005). Because the size of these gratings in these studies extended well beyond the classical receptive field and the orientation tuning sharpened over time, this result suggests that the larger gratings were recruiting a larger number of recurrent inputs that had a delayed and increasingly stronger contribution to orientation tuning. We tested for whether similar behavior occurred for binocular disparity tuning in our model and recordings with DRDS with progressively larger apertures.

We first inspected the steady-state tuning curves in the model and tuning curves of individual neurons in the latest part of the stimulation period analyzed and compared the curves computed from responses to DRDS with varying aperture size. Figure 3, *A* and *B*, shows disparity tuning curves measured from the responses to DRDS with varying aperture size for two example model neurons and two example neurons from recordings, respectively. As the aperture size increased (increasingly lighter curves), tuning curves had narrower peaks and the responses to nonpreferred disparities were relatively suppressed compared with the response to the preferred disparity. Overall, the tuning curves were sharper for the responses to DRDS with larger apertures (Fig. 3, light gray curves).

Next, we looked at how the disparity tuning curves sharpened over time for varying DRDS aperture size for the model and recordings. Figure 4, *A* and *F*, shows the input and initially measured (100–150 ms after stimulus onset) normalized tuning curves for an example neuron in the model and recordings, respectively. For the model, the input tuning curves for the three DRDS aperture sizes were exactly the same (Fig. 4*A*), and for the recordings, the tuning curves for the three DRDS aperture sizes were very similar (Fig. 4*F*). Figure 4, *B* and *G*, shows the steady state and latest measured (450–850 ms after stimulus onset) tuning curves for the same neuron in the model and recordings, respectively. They reveal that over multiple iterations (or time), the response was relatively weaker to nonpreferred disparities compared with the preferred disparity, and the peaks became narrower for larger DRDS aperture sizes (increasingly lighter curves). This change in shape can be more clearly illustrated by plotting the skewness of the distribution of mean firing rates in the tuning curve over multiple iterations or over time. The skewness for both the model disparity tuning and disparity tuning measured from the recorded neuron increased more for larger DRDS aperture sizes (Fig. 4*E*,*J*) compared with smaller DRDS aperture sizes (Fig. 4*C*,*H*).

We summarized the behavior observed in Figure 4 by computing a population average of normalized disparity tuning and skewness over time steps or time for model neurons and recorded neurons for varying DRDS aperture sizes (Fig. 5). The responses were sorted with respect to ranked disparity from the most preferred to least preferred before averaging. As in the examples, for the population average of normalized disparity tuning over time, the responses were relatively weaker to nonpreferred disparities compared with the preferred disparity for DRDS stimulation with the largest aperture size, revealing a clear change in shape (Fig. 5*C*,*K*, black vs light gray). As aperture size increased, you can clearly see that the response to the least preferred disparity became relatively weaker or closer to the dashed line. This change in shape was confirmed by observing a greater increase in average skewness versus log-time for DRDS stimulation with the largest aperture size (Fig. 5*G*,*O*) compared with DRDS stimulation with the smallest aperture size (Fig. 5*E*,*M*). For each recorded neuron (*n* = 81), we performed linear regression on the normalized firing rate versus log-disparity rank for the disparity tuning measured in the latest response window (i.e., a log-fit for the black curves in Fig. 5*I–K*). The average fall-off rate in relative mean firing rate with progressively more nonpreferred disparities (Fig. 6*A*) was significantly larger for the 4-degree DRDS aperture size versus the 2- and 3-degree DRDS aperture size (*p* < 0.001), and the 3-degree DRDS aperture size versus the 2-degree DRDS aperture size (*p* < 0.05). Therefore, the response to the preferred disparity became significantly more prominent with increasing size of a DRDS aperture. Additionally, for each recorded neuron, we performed linear regression on skewness versus log-time. The average slope of the fit increased with DRDS aperture size (Fig. 6*B*) and was significantly positive (*p* < 0.01) and greater for a DRDS with a 4-degree aperture size versus a DRDS with a 2-degree aperture size (*p* < 0.05). Therefore, disparity tuning sharpened more with increasing size of a DRDS aperture.

Although skewness is invariant with respect to changes in the mean and variance of the tuning curve (see Materials and Methods), it cannot distinguish sharpening caused by recurrent interactions versus sharpening caused by an increase in mean firing rate over time coupled with an expansive output nonlinearity, which is part of the disparity energy model (e.g., squaring the response) (Ohzawa et al., 1990). We attempted to avoid this confound by measuring skewness starting at the peak of the population response (100 ms) over an interval where the mean firing rate of the population decreases. To confirm whether or not this was true, we performed linear regression on the mean firing rate versus log-time (Fig. 5*P*), and the mean firing rate significantly decreased over the interval that we measured skewness for DRDS stimulation with all three aperture sizes (*p* < 0.001). Finally, we also examined the mean firing rate with aperture size (Figs. 5*L*, 6*C*) since skewness increased with aperture size. The average mean firing rate (Fig. 6*C*) significantly decreased with increasing aperture size (*p* < 0.001) so the sharpening we observed cannot be explained by an expansive output nonlinearity alone and is incompatible with the disparity energy model. The decrease in mean firing rate with increasing aperture size additionally supports that the increased DRDS aperture size is recruiting a greater number of recurrent inputs outside of the classical receptive field rather than simply increasing the excitatory input to the classical receptive field.

### Disparity tuning dynamics for correlated versus anticorrelated DRDS

In a traditional DRDS, there is 100% correspondence or correlation between the left and right-eye images (Julesz, 1964). Each black and white dot in the left-eye image has a matching black and white dot, respectively, in the right-eye image, but all shifted at the same horizontal disparity. Neurons with binocular disparity tuning in V1 also respond selectively to anticorrelated DRDS (Cumming and Parker, 1997), where each black (and white) dot in the left-eye image has a matching white (and black) dot, respectively, in the right-eye image that are at the same horizontal disparity (inverted polarity) (Julesz and Tyler, 1976). However, the observed modulation amplitude of the inverted disparity tuning curve for anticorrelated DRDS is weaker than the disparity tuning curve for correlated DRDS (Cumming and Parker, 1997; Ohzawa et al., 1997; Nieder and Wagner, 2001). The disparity energy model predicts that the tuning between the two stimuli will be inverted, but have equal strength in modulation (Eq. 4). To try to explain this discrepancy, we examined the tuning curves measured from the responses to correlated and anticorrelated DRDS for both model and neurophysiological data (*n* = 103 neurons). We used the measurement of skewness to reveal and compare the dynamics of disparity tuning for correlated and anticorrelated DRDS stimuli for our model and neurophysiological data.

Examples of the model data (Fig. 7*A*) and data from recordings (Fig. 7*B*) illustrate that the tuning curves based on anticorrelated DRDS (gray) were inverted with respect to correlated DRDS (black) as predicted by the disparity energy model and as reported previously (Cumming and Parker, 1997). Both the model results and data from recordings had tuning curves for anticorrelated DRDS with smaller modulation amplitudes than tuning curves for correlated DRDS. Population averages of tuning curves were generated by sorting the data from the most preferred disparity to the least preferred ranked disparity (both based on correlated DRDS) before averaging. These population averages show that the inverted tuning for anticorrelated DRDS stimulation was consistent across the populations (Fig. 7*C*,*D*). The reduced modulation amplitude for our population of tuning curves based on anticorrelated DRDS stimulation (Fig. 7*F*, *n* = 103, μ = 0.60) was also consistent with previous reports (Cumming and Parker, 1997; Ohzawa et al., 1997; Nieder and Wagner, 2001). However, our model with recurrent interactions also replicated the reduced modulation amplitude for anticorrelated DRDS stimulation (μ = 0.63, although with a narrower and more skewed distribution; Fig. 7*E*), so our model could explain this phenomenon that is not predicted if disparity tuning is modeled with the disparity energy model alone (Cumming and Parker, 1997).

Applying a threshold (Eq. 7) to a disparity energy model neuron can produce reduced modulation amplitude for anticorrelated DRDS stimulation without any recurrent interactions among disparity-tuned neurons (Lippert and Wagner, 2001). However, the threshold only reduces the modulation amplitude for tuned excitatory neurons and actually increases the modulation amplitude for tuned inhibitory neurons (Read et al., 2002). Nonetheless, because there are more tuned excitatory neurons (*n* = 82) than tuned inhibitory neurons (*n* = 21) (Prince et al., 2002b; Liu et al., 2008; Poole et al., 2010), a threshold alone could still result in an average amplitude modulation ratio of less than one. Our model also replicated this bias (see Materials and Methods) so the threshold in our model indeed produced an initial (before any recurrent interactions) average amplitude modulation ratio that is less than one (μ = 0.84). However, previous studies have observed no systematic relationship between amplitude modulation ratio and phase disparity (Nieder and Wagner, 2001) and there are clear examples of amplitude modulation ratios of less than one for tuned inhibitory neurons (Cumming and Parker, 1997). For our neurophysiological data, we did observe that the amplitude modulation ratio is slightly higher for tuned inhibitory neurons (*n* = 21, μ = 0.66) versus tuned excitatory neurons (*n* = 82, μ = 0.58), but the amplitude modulation ratio was still far below one for tuned inhibitory neurons and the difference between the populations was not statistically significant (*p* = 0.27). Although the recurrent interactions in our model did substantially reduce the amplitude modulation ratio from μ = 0.84 to μ = 0.63, the amplitude modulation ratio of tuned inhibitory neurons was still higher than the amplitude modulation ratio of tuned excitatory neurons. However, we could reduce and produce amplitude modulation ratios significantly below one for tuned inhibitory neurons by adjusting the balance between the threshold and the recurrent interactions in the model. Because the threshold increases amplitude modulation and the recurrent interactions reduce amplitude modulation, this was accomplished by relatively weakening the threshold (a more gradual increase in rate or smaller exponent) and/or strengthening the recurrent input.

Not only was the modulation amplitude generally weaker for anticorrelated compared with correlated DRDS, but the shape of the disparity tuning curves differed between the stimuli as well. We inspected the steady-state tuning curves in the model and tuning curves of individual neurons in the latest part of the stimulation period analyzed and compared the curves computed from responses to correlated DRDS to the curves computed from responses to anticorrelated DRDS. Figure 8, *A* and *C*, shows disparity tuning curves measured from the responses to correlated and anticorrelated DRDS for two example model neurons and two recorded neurons, respectively. The tuning curves measured from anticorrelated DRDS stimulation were not simply inverted tuning curves measured from correlated DRDS with weaker modulation amplitudes. To illustrate this more clearly, we normalized the tuning curves by the peak response and inverted the tuning curve measured from the response to anticorrelated DRDS (Fig. 8*B*,*D*, dashed line and open circles, respectively). Tuning curves measured from correlated DRDS had narrower peaks and broader valleys, and secondary peaks were relatively suppressed compared with what we observed in anticorrelated DRDS tuning curves. The responses to nonpreferred disparities were relatively much weaker compared with the preferred disparities for tuning curves measured from the responses to correlated versus anticorrelated DRDS. Overall, the tuning curves were sharper for the responses to correlated DRDS compared with anticorrelated DRDS.

Next, we looked at how model tuning curves and tuning curves measured from responses to correlated and anticorrelated DRDS sharpened over time (Fig. 9). Figure 9, *A* and *C*, shows how example normalized tuning curves for a neuron in the model and recordings, respectively, change over time when presented a standard correlated DRDS. Over multiple iterations or time, the response was relatively weaker to nonpreferred disparities compared with the preferred disparity. Peaks became narrower and valleys became wider. This change in shape can be more clearly illustrated by plotting the skewness of the tuning curve over multiple iterations or over time. The skewness for both the model disparity tuning curve and disparity tuning curve measured from the recorded neuron increased at an approximately linear rate versus log-time-step (iteration) and log-time, respectively (Fig. 9*B*,*D*). Although the tuning curves measured from anticorrelated DRDS stimulation of the same model neuron and example recorded neuron (Fig. 9*E*,*G*) also changed over time in relative magnitude at different disparities, the shape did not appear to change as much and the tuning curve did not sharpen in the same consistent manner as during correlated DRDS stimulation. This qualitative observation was confirmed with the skewness measurement over iterations or time (Fig. 9*F*,*H*) by revealing no clear increase and consistent change in the skewness of disparity tuning during anticorrelated DRDS stimulation.

We summarized the behavior observed in Figure 9 by computing a population average of normalized disparity tuning and skewness over time steps or time for model neurons and recorded neurons during correlated and anticorrelated DRDS stimulation (Fig. 10). The responses were sorted with respect to ranked disparity from the most preferred to least preferred (both based on DRDS stimulation) before averaging. As in the examples, for the population average of disparity tuning over time, the responses were relatively weaker to nonpreferred disparities compared with the preferred disparity for DRDS stimulation with a clear change in shape (Fig. 10*A*,*C*, black vs light gray), which was confirmed by observing that the average skewness increased at an approximately linear rate versus log-time (Fig. 10*B*,*D*). During anticorrelated DRDS stimulation, there was some relatively weakened response for nonpreferred disparities compared with the preferred disparity for model neurons in the population average (Fig. 10*E*), but less than what was observed during DRDS stimulation (Fig. 10*A*). Also, there was almost no noticeable change in shape of the population average of disparity tuning during anticorrelated DRDS stimulation, which was confirmed by observing little change in skewness over time steps (Fig. 10*F*). Any changes in disparity tuning measured from anticorrelated DRDS were even less clear in the population averages of the responses to recorded neurons (Fig. 10*G*,*H*). For each recorded neuron (*n* = 103), we performed linear regression on skewness versus log-time. The average slope of the fit was significantly positive (*p* < 0.001) for correlated DRDS stimulation and was significantly greater for correlated DRDS versus anticorrelated DRDS stimulation (*p* < 0.01).

To again rule out the possibility that increasing skewness was solely the result of sharpening caused by an expansive output nonlinearity, we also performed linear regression on the mean firing rate versus log-time. For the model and this experiment, the overall mean firing rate was stronger for disparity tuning measured from the responses to correlated versus anticorrelated DRDS stimuli (Fig. 7*C*,*D*), so with an expansive output nonlinearity, we predicted that disparity tuning would be overall sharper for correlated versus anticorrelated DRDS stimulation. Indeed, the initial skewness measurements are higher for the responses to correlated versus anticorrelated DRDS (Fig. 10*B*, vs *F*, *D* vs *H*). However, the expansive output nonlinearity does not predict that the skewness would increase for either condition over time because the mean firing rate significantly decreased over the interval that we measured skewness for both correlated (*p* < 0.001) and anticorrelated (*p* < 0.05) DRDS stimulation. Additionally, if we sample tuning curves for correlated and anticorrelated DRDS stimulation so that they have equal distribution of tuning strength (based on DDI), the slope measurements are nearly identical to the measurements based on all neurons: *n* = 29 neurons, 0.10 (*p* = 0.10) versus −0.03 (*p* = 0.58) skewness/log-time for correlated versus anticorrelated DRDS stimulation. This supports that the tuning curves measured from anticorrelated DRDS stimulation are not sharpening over time even when they have the same tuning strength as curves measured during correlated DRDS stimulation. Overall, the sharpening we observed over time for the responses to correlated DRDS cannot be explained by stronger mean firing rates, stronger tuning, or an expansive output nonlinearity.

Even though the statistical tests of the population averages reveal a significant decrease in mean firing rate over time during the interval where we measured skewness, there is diversity in how the mean firing rate evolves over time for individual neurons and for some neurons, the mean firing rate continually increases over time (Samonds et al., 2009). Therefore, we examined the slopes of skewness and mean firing rate versus log-time to make sure that particular extreme examples did not disproportionately contribute to any particular significant trends of skewness over time that we observed. The vast majority of the strong increases of skewness over time coincided with decreases in mean firing rate over time and for the small number of examples where both skewness and mean firing rate increased over time, the mean firing rate increased proportionally less than when mean firing rates decreased over time. We note that in the study by Samonds et al. (2009), even in those examples where skewness was increasing and mean firing rate was continually increasing, there were still features of sharpening over time that could not be explained by an expansive output nonlinearity, such as suppressed secondary peaks. Overall, there was no significant correlation between how skewness or mean firing rate varied over time (*r* = 0.00, *p* = 0.97).

### Disparity tuning dynamics for both experiments for tuned inhibitory neurons

To verify that the behavior was consistently inverted for tuned inhibitory neurons with respect to tuned excitatory neurons, we examined the results for each experiment on the subpopulations separately. As we increased the size of the DRDS aperture, there were greater increases in skewness for tuned excitatory neurons (Fig. 11*A*,*C*, top row, left-to-right) and greater decreases in skewness for tuned inhibitory neurons (Fig. 11*A*,*C*, bottom row, left-to-right). We also compared the changes in skewness measured for excitatory and inhibitory tuned neurons while changing from correlated to anticorrelated DRDS stimulation. During correlated DRDS stimulation, the skewness increased for tuned excitatory neurons (Fig. 11*B*,*D*, top row) and decreased for tuned inhibitory neurons (Fig. 11*B*,*D*, bottom row). During anticorrelated DRDS stimulation there was no clear change in skewness over time for both tuned excitatory neurons (Fig. 11*B*,*D*, top row) and tuned inhibitory neurons (Fig. 11*B*,*D*, bottom row). Overall, these results support that the behavior of tuned excitatory neurons and tuned inhibitory neurons are consistent with each other for the two experiments, but inverted with respect to the direction of change for skewness.

### Recurrent network versus a feedforward network

To deconstruct the contributions of the expansive output nonlinearity, the disparity tuning-dependent connectivity, and recurrence in our network, we measured disparity tuning dynamics for two simple feedforward networks and compared these results to our recurrent network results.

For the first feedforward model, recurrent connections were removed from the original model except for self connections to simulate neural dynamics. This resulted in a feedforward model where only the expansive output nonlinearity could produce sharpening of disparity tuning over time. Because the mean firing rate increased slightly over time in this feedforward model, there was only a small amount of sharpening (<1% increase in skewness) over time and a small reduction in the amplitude modulation ratio from correlated to anticorrelated DRDS stimulation only because of the greater ratio of tuned excitatory versus tuned inhibitory neurons (Fig. 7*E*).

For the second feedforward model, we retained the lateral connections in our original model so that they were still consistent with disparity tuning-dependent connectivity between neurons (Menz and Freeman, 2004; Samonds et al., 2009). However, recurrence was removed so that the remaining lateral connections occurred only in one direction. In other words, we removed the connections at all but one neuron, measured the response of that neuron, and then repeated this procedure for each of the neurons in our population. This effectively modified our recurrent network into a feedforward model with an additional layer of neurons. For the original layer, we had disparity energy neurons representing all possible preferred disparities and spatial frequencies in our model. These neurons then provided inputs to an equal amount of neurons (representing all possible preferred disparities and spatial frequencies) with weighting equivalent to our original lateral connections. One way to visualize this would be to consider that the center neuron in the local circuit in Figure 1 represents an example neuron in the new layer, while all neurons connected to this center neuron represent example neurons in the original layer. The difference in the multilayer feedforward version of our model with respect to the recurrent model in Figure 1 was that all connections from the center neuron in the new layer back to the surrounding neurons in the original layer were no longer present.

In the original full model with both the expansive output nonlinearity and recurrent connections, sharpening due to the nonlinearity and sharpening due to recurrent connections interact with each other significantly. If we remove either the nonlinearity or the recurrent connections, we end up with less sharpening. Without recurrence, the multilayer feedforward model still had 56% of the increase in skewness compared with the recurrent model suggesting that our lateral connections implemented in a feedforward manner alone could account partly for the observation of sharpened disparity tuning. Indeed, previous feedforward models of disparity tuning in V1 have replicated sharpening behavior such as suppressed secondary peaks (Tanabe et al., 2011). However, the mean firing rate increased more substantially over time for our multilayer feedforward model while the mean firing rate decreased over time for the recurrent model, so we cannot rule out that some of the remaining sharpening might have been caused by the expansive output nonlinearity in this feedforward model. Additionally, although the weighted inputs or the expansive output nonlinearity in the multilayer feedforward model alone can produce some proportion of the overall sharpening observed, the recurrent model provides the simplest explanation of the slowly increasing skewness over time, especially considering that the observed mean firing rate of recordings was decreasing. Furthermore, the weighted inputs in the multilayer feedforward model did not replicate the aperture size experiment results and only produced negligible differences (<1%) in the modulation amplitude reduction between correlated and anticorrelated DRDS stimulation beyond what resulted from applying the expansive output nonlinearity alone (Fig. 7*E*). Overall, our model required the weighted inputs between disparity-tuned neurons to circulate through recurrent connections to gain any significant power. An alternative feedforward model than the one we tested could potentially replicate more of the steady-state behavior of our recurrent model given enough freedom of complexity and number of layers. What makes the recurrent model a more appealing explanation, however, is the simplicity of requiring only a single recurrent layer and that the model captures both the steady-state and dynamic behavior that we observed in neurophysiological recordings.

## Discussion

We introduced a simple neural network model with local feedforward responses based on the disparity energy model (Ohzawa et al., 1990) and recurrent connectivity based on observations made from neurophysiological recordings from small populations of V1 neurons (Menz and Freeman, 2004; Samonds et al., 2009). Our model allowed us to produce a rich dataset of dynamic disparity tuning curves. Because the underlying neural architecture is known in the model, we can understand what features of the network caused specific changes in disparity tuning over time. This insight then allows us to make more confident interpretations and predictions about what features of the V1 network are causing similar changes in disparity tuning over time in neurophysiological recordings. We used the statistical measurement of skewness to quantify changes in tuning curves over time, which allowed us to robustly measure features of tuning curve sharpening such as a narrowing peak and suppressed secondary peaks.

The DRDS aperture experiment provides convincing evidence of the role of recurrent inputs in sharpening disparity tuning. As greater numbers of recurrently connected neurons in our model were excited by their mutually preferred disparity with larger DRDS stimulation, there was stronger sharpening of disparity tuning. This was also true for disparity tuning curves measured from our recordings, and similarly, orientation tuning curves are sharper when increasing the size of drifting sinusoidal gratings (Chen et al., 2005; Xing et al., 2005). However, the sharpening that we observed occurred over 100s of milliseconds, while the sharpening observed in orientation tuning-based studies occurred over 10s of milliseconds (Menz and Freeman, 2003). Xing et al.'s (2005) model suggested that the behavior was the result of increased tuned suppressive recurrent inputs. In our model, the behavior was a result of increased facilitative recurrent inputs. However, both models are relatively simple and do not capture the full scale of network interactions in V1. For example, the recurrent inputs outside of the classical receptive field in our model were based on the simplest interpretation of spike timing cross-correlation results (Samonds et al., 2009), but cross-correlation histograms can have multiple potential interpretations with respect to the underlying circuitry (Moore et al., 1970). So even though we did not include tuned suppressive inputs from beyond the classical receptive field in our model, our neurophysiological data does not eliminate the possibility that they might exist between disparity-tuned neurons and might be involved in sharpening disparity tuning. Overall, all three studies, including this study, provide convincing evidence that increasing stimulation that is well outside the most liberal estimates of the classical receptive field of V1 neurons increases the proportion of recurrent inputs that sharpen tuning for a particular feature, such as orientation or disparity, regardless of whether the inputs are facilitative and/or suppressive. If V1 responses are used to infer disparity or orientation, then increasing the aperture size can be interpreted as increasing the amount of evidence about that feature and a sharper tuning curve can lead to a more confident estimate of that feature.

Next, we examined the difference in the disparity tuning dynamics between correlated and anticorrelated DRDS stimulation. The disparity energy model predicted, and previous studies have shown, that V1 neurons have inverted disparity tuning for anticorrelated compared with correlated DRDS stimulation (Cumming and Parker, 1997). The disparity energy model fails to predict the reduced modulation amplitude of disparity tuning for anticorrelated compared with correlated DRDS stimulation (Cumming and Parker, 1997). Our neurophysiological data were consistent with the previous observations and our model was able to capture the reduced modulation amplitude. Additionally, our model predicted reduced firing rates during anticorrelated DRDS stimulation, as well as more complex differences that we observed in the tuning dynamics between correlated and anticorrelated DRDS stimulation. Clear sharpening was measured qualitatively and quantitatively with skewness during correlated DRDS stimulation, while similar behavior was not clearly observed or significant during anticorrelated DRDS stimulation.

In our model, when a correlated DRDS was shown to tuned excitatory neurons, the neurons that responded most were those with positive peaks aligned with the DRDS disparity, and those neurons mutually facilitated each other (Fig. 12*A*). When an anticorrelated DRDS was shown to these same neurons, the neurons that responded most were those with negative peaks (based on correlated DRDS stimulation) aligned with the anticorrelated DRDS, and those neurons had different preferred disparities and spatial scales and therefore, misaligned positive peaks (based on anticorrelated DRDS stimulation; Fig. 12*B*). This misalignment with respect to anticorrelated DRDS-based positive peaks means they had less positive recurrent input, which lead to weaker facilitative interactions. For tuned inhibitory neurons, the neurons that were activated most were those with positive peaks aligned with the DRDS disparity, and those neurons suppressed the tuned inhibitory neuron at the negative peak (Fig. 12*C*). When an anticorrelated DRDS was shown to these same neurons, the neurons that responded most were those with negative peaks (based on correlated DRDS stimulation) aligned with the anticorrelated DRDS, and those neurons had different preferred disparities and spatial scales and therefore, misaligned positive peaks with the tuned inhibitory negative peak (based on anticorrelated DRDS stimulation; Fig. 12*D*). This misalignment with respect to anticorrelated DRDS-based positive peaks means they had less negative recurrent input, which lead to weaker suppressive interactions. From a functional perspective, correlated DRDS represent more natural visual inputs compared with anticorrelated DRDS, which are perceptually confusing and different depths are not perceived (Cumming and Parker, 1997). If the purpose of the recurrent inputs is to perform cooperative stereo computations (Samonds et al., 2009), then it makes sense that they would be organized to deal with the more natural stimuli and sharpen disparity tuning in that case, while they would fail to function during the unnatural and unexpected anticorrelated stimuli.

Our model was reasonably robust to parameter selection. As long as neurons with similar disparity tuning and multiple spatial scales facilitated each other, and there was an expansive nonlinearity, we could produce the primary results reported in this study regardless of our choice of connection weights: (1) disparity tuning sharpened over time (increasing skewness), (2) there was increased sharpening with increasing DRDS aperture size, and (3) there was more sharpening for correlated versus anticorrelated DRDS.

There are, however, two results of our model where parameter selection was not as robust. First, a careful choice in connection weights was necessary to produce model data with a decrease in firing rate over time. To achieve decreasing firing rates over time, the local negative recurrent inputs had to be strong enough where they did play at least some role in sharpening disparity tuning for all neurons in the model. Additionally, stronger negative recurrent inputs produced a greater amount of inverted sharpening (sharpened negative peak) for tuned inhibitory neurons (Fig. 2*E*,*C*). Although Tanabe et al.'s (2011) model was not based on recurrent connectivity, their results also suggest an important general role of suppressive inputs in sharpening disparity tuning by reducing the response to secondary peaks. Second, we had to carefully adjust the balance between the threshold parameters and the overall strength of the recurrent interactions in the model to produce amplitude modulation ratios of less than one for tuned inhibitory neurons. Our network would probably be more flexible about reproducing this result if we included a more realistic input (Haefner and Cumming, 2008). For simplicity, our input was limited to a population of phase-shifted disparity tuned neurons (Ohzawa et al., 1990) and the properties of the tuning curves in V1 suggest that disparity tuning is more complex involving a hybrid of phase-shifted and position-shifted receptive fields (Anzai et al., 1997; Livingstone and Tsao, 1999; Prince et al., 2002b), as well as a combination of positive and negative inputs (Livingstone and Tsao, 1999; Haefner and Cumming, 2008; Tanabe et al., 2011). Our model also does not capture all the potential network, and even additional local, behavior of V1 neurons such as adaptation and feedback (Teich and Qian, 2003, 2006; Schwabe et al., 2006). The motivation of the model was to provide us with more confidence in our original interpretation (Samonds et al., 2009) that there is a link between organized circuitry among disparity-tuned neurons and sharpening of disparity tuning over time. Future experiments, a more complex model, and more complex methods will be required to more definitively decipher the specific contributions of facilitative and suppressive interactions.

The disparity energy model captures a substantial amount of the observed disparity tuning behavior in primary visual cortex. However, the feedforward model fails to capture more complex behavior when we introduce stimulus modifications that encourage interactions among disparity-tuned neurons, which are revealed when we examine the disparity tuning over time. Although there is still much to be explored and debated about the specific underlying computational goals, the evidence of interactions among disparity-tuned neurons and the sharpening of disparity tuning suggests that V1 is playing at least some role in a neural computation that helps to solve the stereo correspondence problem (Samonds and Lee, 2011). Our recurrent model is consistent with the original concept of Julesz (1970) that this solution is a long-range process serving the role of a “search dense surfaces” through the array of local disparities of the image features. Psychophysical evidence supports that such long-range processes are operating to integrate across surfaces in stereoscopic space (Tyler and Kontsevich, 1995; Tyler and Likova, 2011). The present results provide converging evidence about the mechanism tuning and extent of spatial integration underlying these stereoscopic surface interactions.

## Footnotes

This work was supported by NIH F32 EY017770, NSF CISE IIS 0713206, NIH R01 EY022247, Air Force Office of Scientific Research (AFOSR) FA9550-09-1-0678, NIH P41 EB001977, and a grant from Pennsylvania Department of Health through the Commonwealth Universal Research Enhancement Program. We appreciate the technical assistance provided by Karen McCracken, Ryan Poplin, Matt Smith, Ryan Kelly, and Nicholas Hatsopoulos.

- Correspondence should be addressed to Jason M. Samonds, 4400 Fifth Avenue, 115 Mellon Institute, Pittsburgh, PA 15213. samondjm{at}cnbc.cmu.edu