Compound Stimuli Reveal the Structure of Visual Motion Selectivity in Macaque MT Neurons

Anesthetized recording procedures

We recorded from seven anesthetized, paralyzed, adult male macaque monkeys (Macaca fascicularis) and one adult female macaque (Macaca mulatta) using standard procedures for surgical preparation and single-unit recording, as described previously (Cavanaugh et al., 2002). We maintained anesthesia and paralysis by intravenously infusing sufentanil citrate (6–30 μg kg⁻¹ h⁻¹), and vecuronium bromide (Norcuron, 0.1 mg kg⁻¹ h⁻¹), respectively, in isotonic dextrose-Normosol solution (4–10 ml kg⁻¹ h⁻¹). Vital signs [heart rate, lung pressure, electroencephalogram (EEG), electrocardiogram (ECG), body temperature, urine flow and osmolarity, and end-tidal CO₂ partial pressure (pCO₂)] were continuously monitored and maintained within appropriate physiological ranges. Atropine was applied topically to dilate the pupils. Gas-permeable contact lenses protected the eyes, which were refracted with supplementary lenses chosen by direct ophthalmoscopy. Experiments typically lasted 5–7 d at the end of which the monkey was killed with an overdose of sodium pentobarbital. We conducted all experiments in compliance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals and with the approval of the New York University Animal Welfare Committee.

The monkey was positioned so their eyes were 57–114 cm from the display. Grating and plaid stimuli each lasted for 1000 ms and were presented in randomly interleaved blocks. We used 0.5–3 MΩ impedance quartz-platinum-tungsten microelectrodes (Thomas Recording) to make extracellular recordings in the brain through a craniotomy and small durotomy. Electrolytic lesions were made at the end of each recording track for histological confirmation of MT recording sites. For each isolated unit, we determined eye dominance and occluded the non-preferred eye. While isolating neurons in V1 for recording, we selected those with strong direction tuning.

Awake recording procedures

To verify that the observations we made in the anesthetized preparation were not affected by anesthesia, we also recorded from two awake, actively fixating, adult male macaques (one M. mulatta and one Macaca nemestrina). No differences were observed between awake and anesthetized data. A headpost was surgically implanted for head stabilization using the design and methods described previously (Adams et al., 2007). In a second surgical procedure, a chamber was implanted for chronic electrode recording over the superior temporal sulcus (STS) of the left hemisphere, using the techniques and a variant of the design described in (Adams et al., 2011). Before surgery, we used structural magnetic resonance imaging (MRI) and Brainsight software (Rogue Research) to design a chamber with legs matched to the curvature of the monkey’s skull (Johnston et al., 2016) above the STS.

We acclimated each monkey to his recording chair and experimental surroundings. After this initial period, he was head-restrained and rewarded for looking at the fixation target with dilute juice or water. Meanwhile, we used an infrared eye tracker (EyeLink 1000; SR Research) to monitor eye position at 1000 Hz via reflections of infrared light on the cornea and pupil. The monkey sat 57 cm from the display.

The monkey initiated a trial by fixating on a small white spot (diameter 0.1°), after which he was required to maintain fixation for a random time interval between 2350 and 4350 ms. A grating or plaid stimulus would appear 100 ms after fixation began and last for 250 ms. Stimulus conditions were presented in randomly interleaved blocks. The monkey was rewarded if he maintained fixation within 1–1.75° from the fixation point for the entire duration of the stimulus. No stimuli were presented during the 300–600 ms in which the reward was being delivered. If the monkey broke fixation prematurely, the trial was aborted, a timeout of 2000 ms occurred, and no reward was given.

For extracellular recordings, we advanced 0.5–3 MΩ impedance tungsten microelectrodes with epoxylite insulation (FHC) through a 23 Gauge stainless steel dura-penetrating guide tube. We identified area MT from gray matter-white matter transitions and isolated neurons’ brisk, direction-selective responses.

Analysis of neuronal response

We used Expo software (http://corevision.cns.nyu.edu) on an Apple Mac Pro to process signals and sort single units. Signals were amplified 1000×, bandpass filtered (300 Hz to 10 kHz), fed into a multiple-window time-amplitude discriminator, and time-stamped with 100 μs resolution. During each single-unit isolation attempt, discriminator windows and thresholds were manually set, online, to most unambiguously and stably distinguish a single spike wave form from noise and other units. If spike waveforms from multiple neurons overlapped to the extent to which they could not be separated, the discriminator allowed us to detect this based on subrefractory period interspike intervals. If this type of multiunit activity was detected, we would then either manually refine the discrimination windows to effectively isolate the unit, move the electrode a small distance, or abandon that neuron for experimentation and move on to the next unit that could be fully isolated. During the experiments, each unit's waveform was continuously monitored for isolation stability. If necessary, minor, manual adjustments detailed above were made to preserve isolation quality. Spike waveform consistency and interspike interval distributions over the entire duration of the experiments were verified post hoc. While running experiments, Expo automatically generated and displayed tuning curves, from which half-maximum spiking response rates, and their corresponding stimulus values, could be directly computed and displayed through the graphical interface.

Visual stimulation

We presented visual stimuli on a gamma-corrected CRT monitor (an Eizo T966 during anesthetized experiments, and an HP P1230 during awake experiments; mean luminance, 33 cd/m²) at a resolution of 1280 × 960 with a refresh rate of 120 Hz. Stimuli were generated and presented on the same Mac Pro using Expo.

For each isolated unit, we presented windowed sinusoidal grating stimuli to determine, by hand, initial estimates of each cell’s receptive field and preferred size. We used a standard sequence of tuning experiments to make precise estimates of the cell’s tuning preferences. Each tuning curve measured in this sequence featured 100% contrast single gratings varying along a single stimulus dimension, beginning with size tuning and followed by direction, spatial frequency, and temporal frequency tuning. After each of these individual tuning experiments finished, we determined the preferred stimulus value of the dimension tested and used it in subsequent experiments. Next, we measured pattern direction selectivity, at optimal spatial and temporal frequencies, with interleaved drifting gratings and plaids. In the rare cases in which this experiment yielded a different grating direction preference from the previously determined value, we repeated the full sequence of tuning curve measurements, to make sure optimal values would be used for the single component and planar plaid experiments which followed. The receptive fields of all recorded neurons were centered between 2° and 30° from the fovea.

Next, we ran the planar plaid study, which required no further stimulus optimization. In the planar plaid study, stimuli were chosen to span four different direction tuning curves at the optimal spatial frequency (Fig. 7A,B). The first two were based on single gratings at 50% contrast, one with temporal frequency held constant at the optimal value (“constant frequency”), presented in all directions in 30° intervals, and the other with constant, optimal velocity (“constant velocity”) from –90° to 90° relative to the preferred direction, in 15° intervals. Since a given velocity corresponds to spectral content lying on a tilted plane in frequency space, constant velocity gratings had a temporal frequency that varied with the cosine of their direction.

The last two tuning curves consisted of 120° “plaids” (sums of two gratings with orientations 120° apart). The component gratings had the same temporal frequency, and were presented at 50% contrast each. The “pattern direction” of motion (direction consistent with rigid translation, equal to the average direction of the two gratings) was sampled the same set of directions used for single grating tuning curves.

Following the planar plaid study, we ran the single component study. It included presentation of 225 drifting grating stimuli, each at 100% contrast. Stimuli were arranged to widely sample the three dimensions of spatiotemporal frequencies near a given neuron’s tuning preferences, using multiple tuning curves, each varying along a single dimension: direction, spatial frequency, or temporal frequency. These tuning curves were measured at optimal and suboptimal values, the latter of which were determined by reading out (or if necessary, linearly interpolating) the stimulus values which elicited a response at half the neuron’s maximum spike rate (computed directly in Expo, see above, Analysis of neuronal response) in the preceding standard tuning curve experiments (for details, see Extended Data Table 4-3). By sampling spatiotemporal frequencies in this way, we could efficiently concentrate stimuli to reveal subtle changes of each neuron’s selectivity in a manner that does not assume a particular shape of selectivity or manner of tuning specific to either V1 or MT.

Note that the first two direction tuning curves of the single component study differ from the two grating tuning curves in the planar plaid study in that in the planar plaid study: (1) gratings were at 50% contrast instead of 100%, and (2) constant frequency gratings spanned the whole range of directions rather than just the semicircle of directions centered at the preferred one. Although gratings were only presented at 50% contrast in the planar plaid study, its implications for 3D frequency selectivity shape should generalize to the 100% contrast case, since differences in response strength for these contrast levels are negligible (Sclar et al., 1990; Carandini et al., 1997).

Frequency-separable and velocity-separable models

The MT linear weighting functions for both the frequency-based and velocity-based models are defined as a separable product of tuning functions over direction w_d(d), spatial frequency w_s(s), and temporal frequency w_t(t).

Specifically, the frequency-separable linear weighting for a grating is defined as follows: (1)

Direction tuning above is represented by a von Mises function: (2)where μ_d and σ_d represent the direction preference and bandwidth, respectively, and is the modified Bessel function of order 0 (which normalizes the integral of the numerator). Spatial frequency is represented by a log-Normal function, parameterized by spatial frequency preference μ_s and bandwidth σ_s: (3)

Finally, temporal frequency is represented by a Gaussian in coordinates which are linear at low frequencies and logarithmic at higher ones: (4)where g(t) is (5)

Using this functional form for temporal frequency tuning allows to be logarithmic at high temporal frequencies, but also be zero-valued and continuous at zero temporal frequency. The parameter τ determines the temporal frequency at which the function transitions from linear to logarithmic, and μ_t and σ_t are the temporal frequency preference and bandwidth, respectively.

The velocity-separable linear weighting function is defined as follows: (6)where the velocity-separable temporal frequency function, v_t, is defined as a Gaussian, again linear at low frequencies and logarithmic at higher ones: (7)

The only difference between w_t(t) (Eq. 4) and (Eq. 7) is that in the latter, temporal frequency tuning is centered on the preferred speed plane P(d, s): (8)

Note that a consequence of Equations 6–8 is that the velocity-separable model “shears” vertically in the direction of the temporal frequency axis, rather than along the vector orthogonal to the preferred velocity plane (Fig. 1). This was done deliberately, to account for the broad temporal frequency tuning MT neurons exhibit near the preferred direction and spatial frequency (Figs. 4, 8).

Previous models included a V1 normalization stage, either explicitly or implicitly simulated, at this part of the computation (Simoncelli and Heeger, 1998; Rust et al., 2006; Nishimoto and Gallant, 2011). Normalization at the V1 stage has been previously shown to be an important contributor to MT tuning properties. While we could have made it an explicit piece of the model here, Rust et al. (2006) showed that, in some cases, it can be combined with the MT normalization stage to yield a single normalization computation. Moreover, V1 contrast normalization is engaged only when contrast varies widely, and cross-orientation suppression is strongest for components close in orientation. Since the grating components in the plaid stimuli in our experiments are always 50% contrast and 120° apart, we assume such cross-orientation and contrast normalization effects in V1 are negligible. Thus, we assume the next model stage consists of MT neurons summing the responses to each plaid component.

Finally, the full MT model response is computed by raising the linear responses to a power β, and then normalizing them: (9)where the sums are over the components of the stimulus (n = 2 for plaids, n = 1 for gratings), and the α₀ and α₁ parameters represent the spontaneous and maximum discharge rates of the cell. The relative gains of responses to grating and plaid are controlled by the term in the numerator in Equation 9.

The normalization signal, , is meant to approximate the effects of tuned normalization. In the original Simoncelli and Heeger cascade model, MT normalization signals were computed by summing over a simulated population of MT neurons, but this construction would be computationally prohibitive in the context of fitting the model to spiking data. We parameterize the tuning as follows: (10)

This function is maximally active at zero temporal frequency (with a value of 1) and minimally active with a value of γ₀ at t_max. t_max is the highest temporal frequency simulated and experimentally presented. We used this form of normalization for fitting tractability and because it has useful properties, namely, it: (1) ensures there is no suppression at the preferred temporal frequency, (2) can be completely disabled by setting γ₀ = 1, and (3) can be sublinear, linear, or super-linear.

The rationale for this particular approximation of MT tuned normalization is based on the following thought experiment. We start with the assumption that neurons in the MT normalization pool fall into component-, intermediate-, and pattern-selective subpopulations (Fig. 1C). Next, we assume all neurons’ selectivities in each of these subpopulations evenly tile the space of spatiotemporal frequencies. Now we will consider the consequences of summing together the spatiotemporal frequency selectivities of all the neurons in each subpopulation. Specifically, the manner in which these selectivities sum together and overlap each other in frequency space determine if and how normalization is tuned. Any systematic biases in tuning overlap, as a function of spatiotemporal frequency, yield tuned normalization.

Component neurons from both model variants have narrow tuning, overlapping only between neurons with adjacent spatiotemporal tuning preferences. As a consequence, the summed responses of the population of component neurons evenly tile frequency space, producing an untuned normalization signal.

Frequency-separable pattern-selective neurons have broad direction tuning, so their overlap will occur most strongly in direction. The overlap, however, will be separable in spatial and temporal frequency, so for any subpopulation with the same spatial and temporal frequency tuning at all directions, the overlap will be confined to a donut-shaped region centered on those spatial and temporal frequencies. Since we assume the population of frequency-separable pattern selective neurons are evenly distributed across all preferred spatial and temporal frequencies, the tuning overlap will also be evenly distributed, yielding an untuned normalization signal.

Finally, velocity-separable pattern-selective neurons, which are organized along tilted planes which pass through the origin, will have strong overlap at zero and low temporal frequencies regardless of their preferred direction. As such, a pool of velocity-separable pattern neurons tiling frequency space generate a tuned normalization signal which strongly emphasizes low/zero temporal frequency. The same conclusions can be drawn for intermediate neurons in each separable model, although their selectivities will overlap less, yielding a similar, but weaker, tuned normalization signal.

Estimating model parameters for individual cells

In total, the model has nine free parameters for the single-grating study and 10 for the planar plaid study. For the former, they are: the direction preference and bandwidth (μ_d and σ_d), spatial frequency preference and bandwidth (μ_s and σ_s), temporal frequency preference, bandwidth, and log-linear transition (μ_t, σ_t, and τ), and the spontaneous and maximum firing rates (α₀ and α₁). For the latter experiment, μ_s, σ_s, and μ_t are unconstrained by the data and are therefore held fixed at experimentally determined values, but the exponent (β), semi-saturation constant (α₂), and normalization parameters (γ₀ and γ₁) are free. To avoid model fits producing spuriously wide temporal frequency tuning, we included temporal frequency tuning data collected immediately prior in the fitting of the planar plaid dataset. That temporal frequency tuning data, along with the planar plaid stimuli which sample different directions, constrain μ_d, σ_d, and σ_t. In each study, the frequency-separable and velocity-separable models have the same parameters, and only differ in the parameterization of their temporal frequency linear weighting functions, w_t and v_t (Eqs. 4, 7).

For each cell, we optimized the model parameters by minimizing the negative log-likelihood (NLL) over the observed data, assuming spike counts arise from a modulated Poisson model. An additional parameter, σ_G, describes across-trial fluctuations in neural response gain (Goris et al., 2014) and was optimized to the data independently from the frequency-separable and velocity-separable models and held constant during model fitting. We performed the optimization in successive steps, using optimal values from one step as initialization values for the next. First, we fit τ, then added the rest of the MT linear weighting parameters, and then in the case of the planar plaid experiment, the MT parameters controlling the MT nonlinearity.

Experimental design and statistical analysis

In the single component study, we recorded single-unit responses of 13 V1 neurons and 39 MT neurons from seven anesthetized, paralyzed, adult male macaque monkeys (M. fascicularis) and one adult female macaque (M. mulatta). From those same eight monkeys, we recorded 21 V1 neurons and 53 MT neurons for the planar plaid study. These monkeys were also used for other visual system experiments not reported here, 12 of the 53 MT neurons reported had another experiment run after the single component and planar plaid studies were run. Of the 53 MT neurons in the planar plaid study, 23 were also in the single component study. All of the 13 V1 neurons from the single component study are in the set of 21 V1 neurons in the planar plaid study. We additionally recorded 58 MT neurons from two awake, actively fixating, adult male macaques for the planar plaid study (31 from M. mulatta “Albert” and 27 from M. nemestrina “LW”). In all studies, neurons were only excluded from analysis if spike isolation degraded during the experiment, or if spike rates were too low (e.g., always below 10 spikes/s) or variable to reliably predict direction tuning. No additional neurons were rejected from any subsequent analyses.

Following stimulus onset, we counted spikes within a 1000 ms window (anesthetized experiments) or a 250 ms window (awake experiments). For each cell, latency of these windows (relative to stimulus onset) were chosen by maximizing the sum of response variances computed for each stimulus condition (Smith et al., 2005). Error bars on tuning curve responses indicate ±1 SD.

We used standard methods to compute each cell’s pattern index (Movshon et al., 1985; Smith et al., 2005). First, we computed partial correlations between the actual (constant temporal frequency) plaid responses and idealized predictions of pattern and component direction selectivity (r_p and r_c, respectively). We then converted these values to Z-scores to stabilize the variances of the correlations (Z_p and Z_c). Finally, the pattern index is the difference of these two quantities: Z_p – Z_c. Cells were classified as pattern selective if , or component-selective if . Both thresholds correspond to a significance of p = 0.90. Confidence intervals on pattern index were computed from the 95th percentile of 100 bootstrapped estimates (Efron and Tibshirani, 1993; Rust et al., 2006).

For optimal and non-optimal spatial and temporal frequency tuning curves in the single component study in Figure 2B–D, we fit a difference of log2-Gaussians (Hawken et al., 1996). For each neuron, the stimulus value corresponding to the peak of this fitted difference of log2-Gaussians function was used as the fitted preferred stimulus in the figure. To test the robustness of these tuning curve fits, we ran a bootstrap analysis in which trials from each tuning curve were pseudorandomly resampled 1000 times, with replacement, with the restriction that no stimulus condition had zero trials sampled. The error bars in Figure 2B–D represent the 95% confidence intervals of these bootstrapped fitted peak stimulus values. Some tuning curves had flat tops, yielding unreliable tuning preference estimates. We therefore excluded neurons (7 MT, 0 V1) from all analyses in Figure 2B–D, which had a confidence interval exceeding 1.5 decades in any of the three. The conclusions are the same with or without these neurons.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

At the preferred direction, V1 is frequency-separable and MT is velocity-separable. A, Spatial frequency tuning curve data from an example MT cell, measured at three temporal frequencies (error bars denote ±1 SD). The light gray shaded area denotes the spontaneous firing rate, ±1 SD. The fitted SF tuning preferences for the three curves are shown above as triangles. B, The fitted SF preferences from each cell are plotted against the TFs at which they were presented. Both axes are on a normalized scale, representing the ratio of the non-optimal frequencies, relative to the optimal frequency. Each line is the best fit line to the data for one cell. The data along the ordinate axis are aligned to the offset of each best fit line. Lines and points are shaded by the pattern index corresponding to each individual cell. Red corresponds to MT neurons, with darker shades corresponding to higher pattern index, and blue corresponds to V1 neurons, with darker shades corresponding to lower pattern index. The blue and red triangles indicate the mean slopes for all V1 and MT neurons, respectively. Error bars indicate the 95% confidence intervals of 1000 bootstrapped fitted peak stimulus values (for details, see Materials and Methods). C, Same as B but based on TF tuning curves measured at optimal and suboptimal SFs. D, Same as C, but based on TF tuning curves at optimal and suboptimal directions. Here, the points are aligned to the origin, which represents the preferred TF at the preferred direction. Mean slopes in D are computed separately for the different suboptimal directions.

For the constrained parameter search during model fitting, we used a simplex algorithm (the MATLAB function ‘fmincon’). To avoid overfitting and obtain estimates of parameter stability (i.e., the error bars in Figs. 9A,B, 10), we fit the model on 100 bootstrap resamplings of the data. Bootstrapping was done on a per stimulus-condition basis; that is, trials within each stimulus condition were sampled with replacement, ensuring that there were no stimulus conditions without data. Error bars on model fits indicate ±1 SE.

To compare model fits to a given neuron, we computed “velocity superiority”, the difference of the normalized NLLs of the velocity-separable and frequency-separable models. The NLLs were normalized by their corresponding “null” and “oracle” models, which serve as lower and upper bounds, at 0 and 1, respectively. The null model assumes the cell has two possible response rates: one when a stimulus is present and another when there is no stimulus. These are fixed to the measured mean spontaneous and maximal stimulus-driven response rates, respectively. The oracle model serves as an upper bound for the models’ performance, computed by using the measured mean responses to each stimulus condition to predict the neuron’s response to any individual presentation of that stimulus. We used the Wilcoxon signed-rank test to test velocity superiority significance and Pearson’s r to assess correlation between velocity superiority and other quantities, such as pattern index.