Functional specialization in the middle temporal area for smooth pursuit initiation

Experimental procedure

We used a paradigm similar to that of Hohl et al. (2013) to test the influence of stimulus size on smooth pursuit initiation. Stimuli were generated in MATLAB with the Psychophysics Toolbox (Pelli, 1997) and presented on an OLED TV (LG 55EA9800, 122.7 cm × 79.86 cm) at a frame rate of 120 Hz. Each pixel subtended 0.067° of visual angle at a viewing distance of 57 cm, and the screen covered 47° and 35° of horizontal and vertical visual angle. The experiments were carried out in total darkness.

A total of seven subjects, five females and two males, mean age 26.8 years (range 20–34), with normal vision were recruited for this experiment. Subjects, except one of the authors (SB), were recruited via online announcements posted on Student’s Society of McGill University (SSMU) Market Place. Except for one of the subjects (author SB), subjects were naïve and unexperienced in psychophysics experiments. Chin and forehead rests were used to stabilize the subjects’ heads and to maximize the accuracy of eye tracking. All subjects gave written informed consent to participate in the study, which was approved by the Ethics Committee of the Montreal Neurological Institute and Hospital (NEU-06-033).

Figure 1A contains a schematic representation of the task paradigm. Each trial started with a red fixation point subtending 0.2° (luminance 55.2 cd/m²), positioned in the center of the screen. Simultaneous with the fixation point, a red circle subtending 1°, with the same luminance as the fixation point, appeared at 10° in the periphery. The peripheral circle was the main target of smooth pursuit. After a fixation period chosen randomly between 500 and 1500 ms, the main target became dim (luminance 1.5 cd/m²) and a patch of small white dots (each dot 0.1° of visual angle, 0.3% density, 100% coherence, luminance 75 cd/m²) appeared around the main target (Hohl et al., 2013). The low luminance of the target during this time period was chosen to maximize the influence of random dots on smooth pursuit initiation. The patch of random dots and the main target had the same direction, chosen randomly to be either leftward or rightward, and the same speed, which was chosen on each trial randomly between 5° and 20° s^-1.

Figure 1. Schematic representation of the task paradigm and eye tracking preprocessing.

(a) Schematic representation of the task paradigm. Each trial started with a red fixation point, positioned in the center of the screen (red circle at the center). Simultaneous with the fixation point, another red circle appeared at 10° in the periphery. The peripheral red circle was the main target of smooth pursuit. After a fixation period chosen randomly between 500 and 1500 ms, the main target became dim and a patch of small white dots appeared around the main target (small white circles around the peripheral pursuit target). The dashed circle around the white dots marks the border of the random dots patch and was not shown in the experiments. The patch of random dots and the main target had the same direction, chosen randomly to be either leftward or rightward, and the same speed, which was chosen on each trial randomly between 5° and 20° s^-1. 150 ms after the onset of target motion, the patch of random dots disappeared, while the main target continued with the same velocity for another 750 milliseconds. The subject had to saccade toward the moving target (red, low contrast circle) as soon as he/she detected the motion, and follow it with post-saccadic smooth pursuit. At the end of each trial, the subject should go back to the central fixation point. The task had three conditions, with the only difference between the conditions being the size of the patch of random dots (2, 6, or 20 degrees). All the stimuli were presented on a dark background. (b) Top: samples of the horizontal eye displacement. Time zero on the x axis corresponds to the initiation of the target motion. Bottom: The same eye displacement samples aligned at 20 ms after the saccade. The highlighted area (20 ms to 80 ms after the saccade) shows the part of the eye movements that was used as the open loop smooth pursuit for further analysis. (c) Samples of estimated eye velocity. The gray curves show the average of eye velocity for 10 trials for the eye position traces shown in Figure 1B. The black curve shows the average eye velocity for all the trials in one condition (190 trials). The highlighted area (20 ms to 80 ms after the saccade) shows the part of the eye velocity that was used for further analysis.

150 ms after the onset of target motion, the patch of random dots disappeared, while the main target continued with the same velocity for another 750 milliseconds. The subject had to saccade toward the moving target (red, low contrast circle) as soon as he/she detected the motion, and follow it with post-saccadic smooth pursuit. The task had three conditions, with the only difference between the conditions being the size of the patch of random dots (2, 6, or 20 degrees). All the stimuli were presented on a dark background.

Eye movement data analysis

Eye position was recorded using a desk-mounted EyeLink 1000 eye tracker (SR Research) at 1000 Hz. Recordings were monocular (right eye only), and the data were analyzed offline. Saccades were detected from the recorded eye positions using an automated acceleration-based algorithm (Groh et al., 1997), followed by visual inspection to correct false detections. On each trial, the period of 20–80 ms after the saccade was chosen as the open-loop phase of smooth pursuit (Born et al., 2000). After low-pass filtering (0.1–25 Hz), the average velocity of the eye during this period was calculated and passed to subsequent stages for further analysis.

Statistical analysis

The relationship between eye velocity and target velocity was studied with a linear model. Separate models were fitted to the data from each individual subject. The models had the form:

where y_j and x_j represent the eye and target velocities, respectively; α_j is the intercept, β_j is the slope of the fitted linear model, and ∊_j is the error term with distribution 𝒩(0,σ_j) for the j^th subject. In this formulation, the slope and the intercept can vary across subjects. The average of the normalized estimated parameters across the subjects $(\widehat{\alpha },\widehat{\beta })$ were used for group-level analysis.

We also fitted a linear mixed effect model to the data from all subjects pooled together. The model had the form:

where y and x represent the eye and target velocities from all the subjects, respectively, with the target velocity being the fixed effect. As with the individual-subject models, this model can have varying slope and intercept; the subject identity is the random effect. ∊ is the error term with distribution 𝒩(0, σ).

The separately fitted linear models and the mixed effect model produced similar results. Therefore, here we only report the results of the separately fitted linear models, unless stated otherwise.

We used the Statistics and Machine Learning Toolbox^TM in MATLAB 2016 for fitting the linear and the mixed-effect models.

Computational model

To interpret our behavioral findings, we extended an existing model of MT (Hohl et al., 2013) to include recent observations on the relationship between surround suppression and noise correlations (Liu et al., 2016). Except where noted below, all parameters were constrained either by the existing model (Hohl et al., 2013) or by MT recordings (Liu et al., 2016).

Simulation of MT population. We simulated a population of MT neurons with parameters related to direction, speed, and size tuning (Hohl et al., 2013). For each cell, the preferred velocity (𝜗_pref) and direction (θ_pref) were sampled randomly from 0.5-512 degrees/second and 0-180 degrees, respectively. Our method for characterizing the size tuning curves for MT neurons are detailed in (Liu et al., 2016). Briefly, for each stimulus size, we first calculate the neuronal d’ as:

where r_pref and r_null are the mean responses to the preferred and null directions, and ${\sigma }_{pref}^2$ and ${\sigma }_{null}^2$ are the response variances. The resulting size tuning curves are then fitted with the Difference of Error Functions (DoE) function (Equation (4) below) (Liu et al., 2016).

The following equations determine the mean firing rate of each cell:

where the speed, direction, and size tuning are respectively given by f_s(𝜗), f_d(θ), and f_r(p). R₀ indicates the baseline mean firing rate, and G represents the neuron’s gain. (𝜗,θ,p) are the stimulus parameters (speed, direction, and size.) σ_s and σ_τ also determine the tuning curve widths, and σ_e and σ_i are excitatory and inhibitory sizes. A_e and A_i also determine the excitatory and inhibitory heights, respectively. Except for the size tuning parameters, which have been sampled from our MT electrophysiology recordings, all the other parameters are chosen based on the values used in (Hohl et al., 2013). On each trial, the firing rates of the simulated MT neurons are sampled from a multivariate Gaussian distribution with the mean vector of R(𝜗,θ,p) and a covariance matrix calculated based on the method in Shadlen et al. (1996). The correlation structure between the neurons ($r_{noise}^{ij}$) is:

where $\Delta v_{\text{random}}^{ij}$ and $\Delta {\theta }_{\text{random}}^{ij}$ are calculated based on (7) and (8), but the preferred directions and speeds are replaced with random samples from the range of possible values.

An important extension to the model of Hohl and colleagues is the addition of a term α^ij, which models the relationship between surround suppression and noise correlations observed by (Liu et al., 2016); its value is determined by the surround suppression indices of the i^th and j^th neurons. The surround suppression index for the i^th neuron is calculated as $S{I}_i=(d_m^{\prime }-d_L^{\prime })/d_m^{\prime },$, where $d_m^{\prime }$ is the maximum response across the responses to different stimulus sizes, and $d_L^{\prime }$ is the response to the largest stimulus size (Liu et al., 2016). For two neurons with strong surround suppression, α^ij will be close to one, and $r_{noise}^{ij}$ in (9) will not depend on the tuning similarity of the two neurons. The opposite is true for two cells with no surround suppression. This captures the observation (Liu et al., 2016) that pairs of surround-suppressed neurons have noise correlations that are largely independent of the signal correlations, which has important implications for neural coding (Averbeck et al., 2006; see Discussion below).

The surround suppression dependency (α^ij) can be defined in different forms, and indeed its definition affects the structure of noise correlation. We incorporated noise correlations measured from pairs of real MT neurons to obtain a more realistic modeling of surround suppression dependency in our model. We assumed four possible structures for α^ij: 1) α^ij is directly proportional to the sum of the surround suppression indices of the two neurons,

2) α^ij is equal to 0, 0.5, or 1, based on the surround suppression indices of the two neurons,

3) α^ij is drawn randomly from Gaussian distributions with mean values of 0, 0.5, and 1, and standard deviations of 1. The mean value of the Gaussian distribution is chosen based on the surround suppression indices of the two neurons,

and 4) α^ij is equal to zero for all pairs of neurons which represents the structure in which noise correlation does not depend on surround suppression at all.

To choose between these four models, we compared their marginal likelihoods $p(\widehat{r}|M)$. Since our noise correlation models are parameterized by τ_d and τ_s, the Bayesian model evidence is estimated by integrating over τ_d and τ_s:

where τ = (τ_d, τ_s), M_k represents the k^th suggested model of surround suppression dependency as explained above, and k = 1,2,3,4. $\widehat{r}=({\widehat{r}}_1,\cdots ,{\widehat{r}}_N)$ are the measured noise correlations from MT population. For all M_k, p(τ|M_k) is assumed to be a uniform distribution U(0,100). The integral is evaluated on a grid of τ = (τ_d, τ_s) sampled from (0,100] with a resolution of 4.5 for τ_d and τ_s. In this comparison, the third model, described by (12), produced the largest marginal likelihood (see Results and Figure 6). The values of τ_d and τ_s were then chosen to maximize the likelihood

where N is the number of noise correlation samples measured from area MT. The likelihood $p({\widehat{r}}_n|\mathbf{\tau },M_{\text{3}})$ was estimated via repeated simulation of the noise correlation under different values of τ_d and τ_s. In our analysis, τ_d = 45, and τ_s = 4.5 maximized the likelihood function, and were used in our simulations.

Given the noise correlation $r_{noise}^{ij}$, the covariance matrix of the population is then defined based on the method suggested in Shadlen et al. (1996).

In our model, the number of MT neurons in the population that respond to a specific stimulus is determined by the size of the stimulus. The cortical magnification in area MT, Magnification Factor = 6 × eccentricity^-0.9 (Erickson et al., 1989; Van Essen et al., 1981), is used to map visual space in degrees to cortical space in millimeters. Integrating over the cortical space based on this equation yields the cortical area (in square millimeters) that is being activated by a stimulus with a specific diameter. The absolute number of activated neurons can then be obtained by multiplying the activated cortical area by the number of neurons per square millimeter in area MT. In our simulations, we used 20 $\frac{neurons}{m{m}^2}$ to obtain population sizes that are consistent with previous studies (Liu et al., 2016).

Decoding stimulus velocity. Given the simulated MT population responses, vector averaging (Hohl et al., 2013; Lisberger & Ferrera, 1997) was used to decode the stimulus speed from the population activity:

where,

c_i and SI_i are the center of the receptive field and the suppression index of the i^th neuron, respectively. W_D(.) and W_SI(.) determine the contribution of each neuron based on their receptive field location and surround suppression index. In other words, using W_D(.), we limited the read-out to the cells with receptive fields centered within D degrees of the center of the stimulus. Moreover, for different values of the parameter a in W_SI(.), the contribution of neurons changed based on their SI values. For example, a = 0 yields an equal contribution of all the neurons, regardless of their suppression index, in vector averaging, while for a → ∞, only the neurons with suppression indices larger than median(SI) (neurons with surround suppression) take part in estimating motion direction and velocity.

We selected the parameters of the weight functions (i.e. D in W_D and a in W_SI) by comparing our behavioral data with the simulated pursuit velocity for a predetermined range of parameters. Specifically, as a measure of behavioral suppression, we used the observed ratio of pursuit variability at 20 degrees to the variability at 6 degrees (see Experimental Procedure),

where, ${\sigma }_{\text{6}}^j$ and ${\sigma }_{\text{20}}^j$ are smooth pursuit variabilities at 6 and 20 degrees of stimulus size, respectively for the j^th subject. In the model, we simulated ${\widehat{SI}}_{behav}$ on a grid of a and D values (D ∊ [4,8,16,32], a ∊[0: 2: 20]. By projecting the $S{I}_{behav}^j$s on the simulated ${\widehat{SI}}_{behav}$ values, we chose the parameter values that generated the closest behavioral suppression to the average behavioral data. To this end, the 3 nearest neighbors for each $S{I}_{behav}^j$ among the simulated ${\widehat{SI}}_{behav}$s were detected, and the average of their corresponding D and a values were selected for the simulations.

Analysis of neural data. We re-analyzed previously published data recorded from monkey area MT (Liu et al., 2016). In these experiments, the motion stimulus was presented for 50 ms to the monkey while electrophysiology recording was being done simultaneously from area MT (for details see (Liu et al., 2016)). Using linear microelectrode arrays, in each session of experiment, four to seven single units were recorded simultaneously. We used the population of recorded neurons in each session to decode the direction of stimulus motion from the population activity in each trial. Specifically, we would like to see how the decoder’s certainty changes as a function of stimulus size.

To this end, we used two-class linear discriminant analysis (LDA) (Bishop, 2007) to decode the direction of motion from the activity of the neurons in each trial. Given the LDA prediction ($\widehat{\theta }$) and the true stimulus motion direction (θ) for each trial, we calculated the uncertainty coefficient of the decoder as:

where, I(θ, $\widehat{\theta }$) is the mutual information of the decoded and real stimulus direction, and H(θ) is the entropy of the real stimulus direction (Cover & Thomas, 2012). The uncertainty coefficient quantifies how precisely the decoded direction predicts the real direction of the stimulus. Importantly, high uncertainty is assigned to a decoder which is inconsistent in decision making; wrong but consistent decisions do not result in high uncertainty. When sampling from a trial to another trial from the decisions of an inconsistent decoder (with high uncertainty coefficient) for a constant stimulus direction, one obtains more variable decisions, compared to those of a consistent decoder. In other words, the output of an uncertain decoder has high trial-to-trial variability. Therefore, as a proxy to behavioral variability, we measured the decoder uncertainty for different stimulus sizes. We then calculated the weighted average of uncertainties over the recording sessions. Once the weights assigned to each session was proportional to the number of surround suppressed neurons (SI_i ≥ median(SI)) red curve in Figure 2A) and once proportional to the number of non-surround suppressed neurons (SI_i < median(SI)) (blue curve in Figure 2A). Comparison of these two weighted averages in Figure 2A shows the effect of surround suppression on the decoder uncertainty.

Figure 2. Direction decoding uncertainty from middle temporal (MT) population.

(a) Direction decoding uncertainty. We decoded the direction of moving stimuli based on the responses of groups of simultaneously recorded MT neurons. The uncertainty of the decoder in assigning directions to the stimuli is plotted as a function of different stimulus sizes. The uncertainty of trained decoders are averaged once weighted based on the number of non-surround suppressed neurons in the sets (blue curve) and once based on the number of surround suppressed neurons (red curve). Therefore, the red curve shows the average decoder uncertainty when the majority of neurons in the population have surround suppression, while the blue curve shows the uncertainty for a population mainly including surround supressed neurons. (b) A schematic representation of our hypothesis. The first column represents the stimulus in three different sizes. The second column schematically demonstrates the activity that each stimulus evokes in area MT. The small circles and the triangles show the neurons with strong and weak surround suppression, respectively. The dotted lines within the boxes separate the surround-suppressed and non-surround-suppressed neurons. The stimulus induced activity level is color coded and explained by the size tuning curves at top row. Blue and red colors represent low and high firing rates, respectively. The largest stimulus evokes the highest average level of activity in neurons without surround suppression (red triangles) while it induces a lower firing rate in surround suppressed cells (blue circles). The third column represents the decoder. f(MT) estimates the velocity of the stimulus in the first column given the activity of MT population in each trial. According to one hypothesis, the neurons with strong surround suppression have a stronger contribution (solid gray arrow) in the read-out compared to the neurons with weak surround suppression (dashed white arrow.) The trial-to-trial variability of the estimated velocity is depicted in the fourth row. The least variability, among these three stimulus sizes, is given by the medium sized stimulus (middle row) where the signal to noise ratio is higher than the two other conditions for the cells with surround suppression.

Analysis of neural data from MT

To generate predictions based on the properties of different neuronal populations, we re-analyzed previously published data recorded from monkey area MT (Liu et al., 2016). Specifically, we decoded the direction of stimulus motion from groups of four to seven simultaneously recorded MT neurons. We then used an optimal decoder (see Methods) (Powers, 2011) to estimate the uncertainty associated with the neural representation of velocity in the MT population. As shown in Figure 2, the uncertainty varies with stimulus size, and this relationship is different for surround-suppressed and non-suppressed subpopulations. For a population of neurons with surround suppression, there is an increase in the decoder’s uncertainty for large stimuli (Figure 2, red curve), while the uncertainty decreases and then plateaus with increasing stimulus size for the neurons with weaker surround suppression (Figure 2, blue curve).

Psychophysics experiments

In our psychophysics experiments, human subjects pursued a moving target that was surrounded by coherently moving dots (Figure 1A). For each subject and each trial, we quantified the first 60 ms of the post-saccadic smooth pursuit velocity. Figure 1B shows the raw eye position data for one subject, as well as the segment of the pursuit eye movement that was used for further pursuit velocity analysis (pink). Figure 1C shows examples of measured post-saccadic eye velocity.

In each trial of the task, a random velocity between 5 and 20 degrees/sec was assigned to the pursuit target. The scatter plot in Figure 3A shows the eye velocity as a function of target velocity on every trial in one subject.

Figure 3. Sample of regression analysis on one subject.

Top: The gray circles represent the single trials for one of the conditions (target size = 6 degrees.) The horizontal and the vertical axes show the target and the eye velocities, respectively. The green dashed line demonstrates the linear model fitted to the measured target/eye velocity pairs. Bottom: The Root Mean Square Error (RMSE) of the fitted lines in (a), for the three conditions. The error bars show the standard error on samples of RMSE values estimated by bootstrapping the data.

In Figure 3A, the small gray circles represent the eye and target velocity in each trial. To measure the accuracy and precision of the smooth pursuit system, we then fitted a linear model to the pairs, as explained in the Methods section. This process was repeated across the subjects and across task conditions (i.e. different sizes of random dots patches.)

For an ideal smooth pursuit system, we expect to obtain a linear fit with a slope of one, zero intercept, and zero noise variance. This hypothetical system is able to follow the target velocity exactly (unity gain, no bias) and with perfect reproducibility across trials (no variance). The parameter that is of interest in this study is the one that quantifies the variability of the real smooth pursuit system. In the linear model shown in Figure 3A, the Root Mean Square Error (RMSE) captures the variability of smooth pursuit.

In order to observe the effect of target size on the pursuit variability, we have plotted (Figures 3B) the estimated RMSE of the linear fits for the three target sizes in a representative subject (the same subject as in Figure 3A). As can be seen in the figure, increasing the diameter of the patch of random dots from 2° to 6° decreases the variability of pursuit initiation (Cohen’s d effect size = 1.59, Wilcoxon signed rank test: z = 26.92, p < 0.0001), while increasing the patch size from 6° to 20° increases the variability (Cohen’s d effect size = 1.89, Wilcoxon signed rank test: z = 27.39, p < 0.0001). Thus for this subject the relationship between eye velocity and stimulus size is consistent with the idea that pursuit initiation is driven by the population of surround-suppressed MT neurons (Figure 2B).

We next examined the relationship between pursuit variability and stimulus size for our cohort of seven subjects. To summarize the data from all subjects, we fitted linear models to all the data we acquired from individual subjects (see Methods for the details). The lines shown in Figure 4A are the averages of the estimated models across the subjects for the three conditions. For a simpler presentation of the data, and to make the plot less cluttered by individual data points, we grouped single pairs of eye/target velocity to quantiles. Each circle in Figure 4A represents one of these estimated target/eye velocity quantiles. The differences in the RMSE of the fitted models, as can be seen in this figure, capture the effect of stimulus size on smooth pursuit variability. The changes in the average pursuit variability across subjects is plotted in Figure 4B. The estimated RMSE values are normalized for each subject before averaging. The error bars are the standard error across subjects. The results indicate that, although there was substantial variability, the multi-subject analysis yielded effects similar to those shown for the example subject in Figure 3 (see Figure S1 for the estimated RMSE for the individual subjects). For the pursuit variability, we see a decreased RMSE from 2° to 6° (Cohen’s d effect size = 0.4376, Wilcoxon signed rank test: z = 1.5978, p = 0.1101), and a large increase from 6° to 20° (Cohen’s d effect size = 1.3112, Wilcoxon signed rank test: z = 2.0863, p = 0.0372). As reported previously (Born et al., 2000), there was also a decrease in pursuit gain when stimulus size increased beyond 6 degree (Cohen’s d effect size = 1.7473, Wilcoxon signed rank test: z = 2.5893, p = 0.0096).

Figure 4. Linear models averaged across subjects.

(a) The pooled data of all the subjects for the three conditions are shown here (Black = 2°, Green = 6°, Orange = 20°.) The horizontal and vertical axes show the target and the eye velocities, respectively. Each circle represents one of the estimated target/eye velocity quantiles. On the horizontal and vertical axes, the marginal values of the target and eye velocity data points are shown. The lines overlaid on the samples represent the averaged fitted mixed-effect models obtained through averaging the estimated β_js and α_js. (b) The average Root Mean Square Error (RMSE) of the fitted linear models is shown for the three conditions. The error bars show the standard errors.

In general, the results presented in this section showed that increasing the size of the pursuit target can decrease and then increase the pursuit initiation variability. As briefly explained in the Introduction, this pattern of results could arise from the sensitivity of the neurons involved in the sensory representation of target motion. To explore this point in greater detail, we performed simulations of a detailed computational model based on the properties of MT neurons.

Computational model

To suggest a plausible neuronal mechanism for the behavioral effects reported in the previous section, we devised a model that combines well-established aspects of motion coding (Hohl et al., 2013) with recent data on noise correlations in MT (Liu et al., 2016). We then tested different versions of this model on the same task as in our psychophysics experiments.

As detailed in the Methods, the model has two main stages. The first transforms the stimulus into a pattern of activation across a simulated population of MT neurons, with standard functions representing tuning for direction, speed, and size. We also included noise correlations among neurons and the dependency of these correlations on other tuning parameters (Liu et al., 2016). This provided an estimate of the neuronal response across MT for each stimulus used in our experiments.

The second stage involves a readout of the MT population, in which vector averaging is used to map MT activity into smooth pursuit velocity (Lisberger, 2010). Our formulation of the vector averaging operation provided the possibility of assigning different weights to different neuron types (surround-suppressed versus non-suppressed), which allowed us to consider different hypotheses about the readout of the MT population. In principle, the readout could make use of the entire MT population, as is often assumed. Or it could rely primarily on the MT subpopulation with strong surround suppression, as suggested by the results of Born et al. (2000). To evaluate these possibilities, we fit the model weights to the behavioral data reported in the previous section. The parameters of the simulated MT population (e.g. noise correlation structure) were chosen to maximize the similarity of the simulated and real MT population (see the Methods for details).

In the model, the vector averaging operation involves parameters corresponding to the integration area in visual space and the weights assigned to each neuron depending on its suppression index. The integration area limits the vector averaging only to the neurons with receptive field centered within a specific distance (D) from the stimulus center, as suggested by previous data (Mukherjee et al., 2016). The suppression index-related parameter (a) determines the contribution of each neuron in the vector averaging based on its surround suppression index.

Manipulation of these two parameters led to a range of models that varied in terms of the behavioral surround suppression (i.e. the ratio of the variability at 20 degrees to the variability at 6 degrees of stimulus diameter; Equation (18)). As detailed in the Methods, we used an optimization procedure to select the values of D and a (Figure 5). In short, we evaluated the behavioral surround suppression index of our model (Equation 18) for different values of D and a on a two-dimensional grid, and the parameter values that yielded the closest fit to our experimental data were selected. The optimal parameters were D = 4° and a = 15, which correspond to a modest integration radius (Liu et al., 2016; Mukherjee et al., 2017) and a strong reliance on surround-suppressed MT neurons (Born et al., 2000).

Figure 5. The variability of smooth pursuit initiation simulated with our model.

The trial-to-trial variability of the simulated pursuit velocity for different stimulus sizes of 0.4–20.4 degrees. The red curve illustrates the variability for the simulations based on a = 15 (optimal value; only surround suppressed neurons contribute in the readout). The blue curve shows the simulations corresponding to a = 0 (equal contribution of all neuron types in the readout). Both simulations are based on D = 4 (optimal value). The error bars are calculated by iterating the simulations 1000 times.

Figure 6. The average firing rate (left) and signal-to-noise ratio (SNR, right) for the population of simulated MT neurons.

The insets are the enlarged version of the marked parts in the main graphs. The circled points A and B correspond to the firing rate (on left) and SNR (on right) at 2.2° and 4.4° of stimulus size, respectively. The red and blue curves correspond to the neurons with and without surround suppression.

Once the parameters were established, we calculated the variance of the model output by setting the pursuit target motion parameters (direction and velocity) to specific values, and iterating the simulations 1000 times. In the model, variability arises from the noise distributions assigned to each neuron and the noise correlation structure across the population (see Methods, Equation 7–Equation 9). Figure 5 shows the variance of the model output as a function of stimulus size for the optimal parameter values (red line); for comparison we also show the model output for the case in which equal weight is assigned to all neurons, regardless of surround suppression strength (a = 0; blue line).

As can be seen in Figure 5 (red curve), increasing the size of the target initially (from 0.4° to 4.4°), decreases the trial-to-trial variability of the decoded velocity, while further increases in stimulus size increase the variability. This pattern is similar to our psychophysical observations. In contrast, when surround-suppressed MT neurons are not accorded greater weight (Figure 5, blue line), the decoding velocity shows little dependence on size beyond a certain point.

Two factors can explain the initial increase in the precision of smooth pursuit, which was observed in both our psychophysics and simulated data. Relative to the smallest stimuli, larger stimuli induce a larger population response, both by filling the classical receptive fields of the neurons and by activating more neurons; the same two factors improve the signal-to-noise ratio (SNR) of the MT population (Figure 6). The more interesting phenomenon is the decrease in pursuit precision with further increases in stimulus size. In this case, a model that gives equal weight to all MT neurons predicts a saturation in pursuit precision with increasing stimulus size (Figure 5, blue lines), while a model that gives greater weight to surround-suppressed neurons predicts the size-dependent effects we have seen in our data (Figure 5, red lines). The decreased precision in the model arises from the decreased sensitivity of surround-suppressed MT neurons for increasing stimulus size. This decreased sensitivity for large stimuli is evident in the average firing rate of the surround-suppressed MT neurons (Figure 6, left, red line) compared to that of the non-surround-suppressed neurons (Figure 6, left, blue line).

Sources of motor variability

Our results have important implications for understanding the neuronal basis of behavioral variability in smooth pursuit initiation (Osborne et al., 2005). Previous studies attributed a large portion of this variability to the noisy responses of MT neurons. Specifically, by showing that noise correlations in area MT exhibit a specific structure that is detrimental to motion coding, Huang & Lisberger (2009) claimed that the neuronal trial-to-trial variability propagates from area MT all the way through to behavior. This hypothesis, of a sensory source for motor variability, is based on the assumption that area MT is a homogenous population of neurons with identical noise correlation structure and equal contributions to the encoding of stimulus motion during smooth pursuit initiation. However, Liu and colleagues showed recently that the noise correlation between MT neurons with surround suppression is less detrimental for motion coding than other MT subpopulations (Liu et al., 2016). Specifically, because the noise correlations are largely unrelated to the signal correlations, the noise in the responses can be averaged out, to a large extent, over the encoding population. Hence the variability in the output of MT could be considerably lower than previously suspected, such that the limiting factors for pursuit precision are the sensitivity of the MT population and the noise introduced as part of the sensorimotor transformation. This interpretation leaves open the possibility that other areas downstream of area MT also contribute in the variability of smooth pursuit initiation (Lee et al., 2016).

Surround suppression in the gain of smooth pursuit

Previous studies have largely focused on the gain of smooth pursuit and its relationship to stimulus parameters. For example, low-contrast stimuli shift the population of MT neurons toward an overestimation of stimulus speed, as neurons that encode high speeds become active for slow stimuli (Pack et al., 2005). At the same time, pursuit gain actually decreases with decreasing contrast (Priebe & Lisberger, 2004). To reconcile these results, Priebe & Lisberger (2004) made use of a modified version of vector averaging, in which an additional term was introduced in the denominator to decrease the estimated velocity when the average firing rate in MT is low.

This approach can be accommodated within the framework of Bayesian motion coding (Weiss et al., 2002), wherein, low contrasts are associated with less precise stimulus coding, which leads to a greater reliance on prior assumptions. In models that include a prior belief that objects are stationary (i.e. zero mean velocity prior), decreasing the precision of image motion leads to a lower perceived velocity; this effect has been confirmed psychophysically (Ascher & Grzywacz, 2000; Hürlimann et al., 2002; Stocker & Simoncelli, 2006).

The relationship between the current results and the Bayesian coding scheme is not entirely clear, as the larger stimulus can in principle be integrated so as to reduce the precision of the velocity estimate. However, it does not appear that the brain performs this integration, given the observation that pursuit gain decreases for large stimuli (Born et al., 2000), while the variability increases (Figure 4). These findings highlight the importance of considering the neural representation of the stimulus, which in this case appears to be determined primarily by surround-suppressed MT neurons.

It should also be noted that the Bayesian model depends heavily on the form of the prior assumption, and it is not difficult to find experimental situations in which the low-velocity prior is violated (Thompson et al., 2006). Thus, the relationship between MT activity and the gain of pursuit initiation is likely more complex than what has been captured by existing models (Brooks et al., 2011). We have largely avoided this issue by focusing on the variability of pursuit initiation.

Relation to ocular following response

A similar behavior to smooth pursuit eye movement is the ocular following response (OFR), which is a short-latency involuntary eye movement elicited by sudden movement of large-field visual stimuli (Gellman et al., 1990; Miles et al., 1986). These movements are important for stabilizing the eyes on the visual scene. The OFR and smooth pursuit eye movements share part of their neuronal pathways, including brain areas MT and MST, as well as other subcortical regions (e.g. cerebellar nuclei) (Kawano et al., 1994; Optican et al., 1986). Given their commonalities, one might suspect that the OFR system partially interferes with the smooth pursuit system in our psychophysics experiments, and could be responsible for at least part of our behavioral effects. Particularly, given the fact that OFR is triggered by large-field moving stimuli, the largest stimulus in our psychophysics task (20 degrees) could potentially engage the OFR system. However, certain aspects of our results challenge this possibility. First: as mentioned in the Results and reported elsewhere (Born et al., 2000), the gain of eye movements in our task decreases for the largest stimulus compared to the smaller one. This contradicts our expectations of OFR gain which should increase monotonically with stimulus size (Masson & Perrinet, 2012; Miles et al., 1986). Second, in our paradigm, subjects are instructed to follow a stimulus whose motion is consciously anticipated. Unlike our task, in conditions where OFR has been reported, the subjects were unaware of stimulus motion and their eye movements were involuntary (Glasser & Tadin, 2014). Third, our moving texture (random dots) appears on the screen at the same time that the main target starts to move in the periphery (see Experimental Procedure). In the OFR paradigms, however, the stationary texture on the screen starts drifting at an unexpected time. Moreover, the absence of foveal motion signal decreases the OFR gain significantly (Miles et al., 1986). Together, these aspects of our task paradigm and behavioral results reduce the chance of OFR system involvement in our observations, albeit do not entirely dismiss its possibility.

In summary, we showed that the neurons with surround suppression in area MT are the main contributors to smooth pursuit initiation. This causes a specific pattern of changes in the trial-to-trial variability of pursuit initiation that was validated by our psychophysics experiments. In addition, our proposed model was shown capable of reproducing this behavioral pattern.