Population Dynamics of Early Visual Cortex During Working Memory

Subjects:

Using procedures approved by New York University, five healthy subjects (3 female; right-handed; 25–42 years old) gave informed consent and participated in the study. Each subject participated in six scanning sessions: one to obtain three high-resolution anatomies, one to measure retinotopy, and four to measure responses during memory-guided prosaccade and antisaccade tasks. Portions of the data were previously published in Saber et al. (2015).

Memory-Guided Saccade Tasks:

Subjects performed two memory-guided saccade tasks (Figure 1). Subjects were instructed to remember the location of a briefly presented visual stimulus and immediately plan a saccade to, in the case of memory-guided prosaccades, or away from the stimulus, in the case of memory-guided antisaccades. The stimulus, a high contrast flickering visual checkerboard (1° diameter, spatial frequency 2 cpd, 8 Hz flicker), appeared at one of eight locations separated by 30° angles at 10° eccentricity, avoiding vertical and horizontal meridians. Each trial began with 1.25s of central fixation followed by a brief presentation of the visual stimulus. Subjects maintained the planned prosaccade or antisaccade over a long and variable delay period (3–15s). Disappearance of the fixation along with a brief auditory tone prompted the memory-guided saccade. Subjects held their gaze until the correct location was presented as feedback and corrected any error in their gaze. Subjects then returned to fixation for the inter-trial interval (ITI, 10.5–15s). Partial trials were also included to help break the dependency of the delay on the target. These trials were aborted after the visual stimulus, signaled by a change in the color of the fixation spot. In four sessions, each with eight 322s scans, subjects performed 192 prosaccade and 192 antisaccade trials, blocked by condition.

Oculomotor Methods:

Eye position was recorded in the scanner at 1000 Hz (EyeLink 1000, SR Research, Ontario, Canada). We computed error in degrees of angle of the primary saccade as well as the final error after any quick corrections were made prior to the target feedback. All trials were inspected for noncompliance. Trials in which subjects made an eye movement to the wrong saccade target were rare (0.3% and 0.6% of all pro- and antisaccade trials, respectively). We excluded trials where gaze deviated from fixation greater than 2° during the delay-period or immediately following visual target (4.6% and 2.7% of pro- and anti-saccade trials, respectively).

MRI Acquisition:

All MRI data were acquired on a 3T Allegra head-only scanner (Siemens, Erlangen, Germany) using a head transmit coil (NM-011) and two surface receive coils (NOVA Medical, Wakefield, MA): (1) a four-channel phased array receive surface coil (NMSC-021) for retinotopic mapping, (2) a four-element phased parallel array (NMSC-011) for memory-guided saccade tasks. We acquired T2*-sensitive echo-planar images (repetition time (TR) 1.5s; echo time 30ms; flip angle 75°; 26 slices; 3×3×3mm voxels; 192×192mm field of view). Three T1-weighted MPRAGE scans were averaged and used for gray matter segmentation, cortical flattening, registration, and visualization for creating ROIs (see below). Functional scans were corrected for head motion and aligned across sessions, detrended and high-pass filtered with a cutoff frequency of 0.0167 Hz, and converted to percentage signal modulation.

Retinotopic Mapping:

We defined retinotopic visual areas (V1, V2, V3, V3A/B, IPS0) using a standard phase encoded approach (Wandell, Dumoulin, & Brewer, 2007). Polar angle components of the retinotopic maps were estimated with a high contrast flickering wedge that rotated (clockwise or counter-clockwise) around central fixation (subtended 90° of polar angle). Radial components were estimated with a high contrast flickering annulus expanding or contracting from fixation (subtended 4.2° radius). Each scan consisted of 10.5 cycles with a period of 24s, and lasted 252s. Each subject completed six to eight runs of the rotating wedge aperture and four runs of expanding/contracting annulus aperture. We used flattened cortical surface representations of each subject’s occipital and parietal cortices to visualize amplitude, coherence and phase maps. Visual area boundaries were drawn on retinotopic maps based on standard conventions (Larsson & Heeger, 2006; Wandell et al., 2007). We used voxels from these ROIs in our IEM described below.

Estimating Delay Period Activity:

We used a voxel-wise generalized linear model (GLM) to estimate the responses of each voxel to the stimulus encoding, memory retention interval (i.e., delay period), and memory-guided saccade intervals of each trial (Srimal & Curtis, 2008). Each epoch was modeled by the convolution of a canonical model of the hemodynamic impulse response function with a square wave (box-car regressor) whose duration was equal to the duration of the corresponding interval, detrended and high-pass filtered. Importantly, we estimated beta coefficients for each of the three trial epochs (stimulus, delay, response) for every trial independently. Modeling every event in each trial guaranteed that the beta coefficient corresponding to each epoch estimated the magnitude of BOLD activity specific to each epoch. As described below, we used the beta coefficients from the stimulus encoding epoch of the prosaccade task for model training and the beta coefficients from the delay periods of the prosaccade and antisaccade tasks for model testing.

Inverted encoding model (IEM):

To reconstruct the spatial representation of the retinotopic maps, we used an IEM (Sprague & Serences, 2013). The IEM assumes that the activity of each neuron is a linear combination of the outputs of a set of basis functions, or information channels, tuned to different stimulus values. Similarly, it assumes that the activity of each voxel, modeled as the GLM beta coefficients, is a weighted sum of the activity of all neurons in that voxel. As a result, the activity of each voxel can be equation ted as a linear combination of the output of these basis functions for a given stimulus. This can be shown mathematically as:

Where B is a matrix that represents the GLM beta coefficients corresponding to an event of interest (i.e., stimulus encoding, delay period) for all voxels in a given ROI and for all trials. W is a matrix of regression weights and C contains the coefficients of the basis functions. The spatial representation of visual space is a linear combination of all basis functions, each determined by its corresponding coefficient in C. To calculate the channel coefficients in a given trial using equation (1), we use regression weights, W, that link the outputs of the basis functions to each voxel’s activity. Thus, we can rewrite equation (1) as:

Where B_train in (2) is a subset of the data used for training. However, during training we do not have access to the real contribution of each basis function to the population activity and thus we need to use a set of hypothetical channel coefficients, C_train. For each trial, we calculate these hypothetical coefficients as the overlap between each basis function and that part of the visual field covered by the stimulus in that trial. Using the regression weights, $\widehat{W}$, and the population activity, B_test, in a given trial, we can calculate each channel’s coefficients as:

To avoid circularity in our analysis, we cross-validated our model by using the GLM derived beta coefficients corresponding to the stimulus encoding epoch of the pro-saccade task to calculate the regression weights, $\widehat{W}$. We then used these regression weights to estimate the channel coefficients, C_test in (3), given the population activity, B_test, during the delay period in the pro-saccade and both encoding and delay periods in the anti-saccade task. For the stimulus encoding epoch of the pro-saccade task, we used a leave-one-out method to train and test the model. For each trial in the dataset, we took that trial out and used the rest of trials to calculate the regression weights, $\widehat{W}$. We then used these regression weights to estimate the channel coefficients, C_test in (3), given the population activity, B_test,, corresponding to the excluded trial. Thus, adding all basis functions, each multiplied by its corresponding channel coefficient in C_test, gives us the reconstructed spatial representation of visual space in that trial. Figure 1C shows a schematic of how our model works.

Based on the spatial selectivity of voxels in each ROI in early visual cortex, we constructed a set of basis functions consisting of 2-D Gaussian filters with equal FWHMs distributed uniformly along a circle corresponding to a 10° eccentricity (i.e., eccentricity for the visual targets and saccade goals). The number and width of these basis functions were chosen in a way to both avoid rank-deficiency (caused by too much overlap between channels) and obtain smooth reconstructions. To make sure that using basis functions restricted to the stimulated ring does not bias the reconstructed spatial representation, we also tested a set of basis functions tiling the whole visual field. Since we did not find any meaningful differences in the reconstructed visuospatial maps, here we only report the results from the basis functions restricted to the actual eccentricity used in the experiment.

Statistical analysis:

To evaluate the differences in behavioral performance between pro- and anti-saccade tasks, we measured three parameters for each trial; the primary and final saccade accuracies and the reaction time. We then, for each parameter, performed a permutation test between the measurements in each task, and a random subset of concatenated measurements in both tasks.

To verify that the IEM can reconstruct the WM content uniformly well across different saccade target locations, we plotted the model’s reconstructed locations versus the true target locations. Within each topographic map, we calculated the peak of the reconstructed heat map averaged across all trials for each saccade target location and plotted this against the true target locations. Then, we estimated the slope of lines fit to the these points, where the closer the slope of a line was to one indicated a more precise readout across different locations. Moreover, to evaluate the model performance on a trial-by-trial basis, we defined the readout precision at each trial as $1-\frac{\mathit{\text{error}}}{180}$ where error is the difference between the angular position of the peak in the reconstructed retinotopic map and the actual location in the corresponding trial. We then, measured the skewness of the distribution of precisions. The precision varies between 0 (maximum error) and 1 (minimum error). Thus, the higher the trial-by-trial model performance, the more skewed the precisions distribution towards 1. We measured the significance of the skewness towards 1, by comparing the skewnes corresponding to the original and randomized order of data through a permutation test. Moreover, we identified correctly decoded trials based on the distance between the IEM readout peak and the actual target location in each trial. We then, through a permutation test, compared the proportion of correctly decoded trials from the original and randomized order of data in each ROI. To do so, we shuffled the trial orders 1000 times, for each subject and each ROI, and reconstructed the visual-spatial maps each time. Then, we ran a student t-test between the proportion of correctly decoded trials corresponding to the original order and shuffled.

To measure the spatiotemporal interaction of the reconstructed representation of population activity under different conditions, we fit a 2-D Gaussian to the reconstructed visual-spatial maps and measured the gain, width and precision of the fit. We then used a 2-way ANOVA to test the significance of the interaction between each of these parameters at two time epochs, the encoding and delay periods, and two locations, the stimulus location and saccade goal, which is the same on pro-saccade, but different on anti-saccade trials. We measured the Pearson’s correlation for each voxel between the delay period GLM beta coefficient across trials and the trial-by-trial success of our model’s decoding of memorized space. To identify all voxel clusters across the whole brain that significantly predict the model performance in an ROI, we used, independently for each voxel, a permutation test to compare the z-transformation of above correlations corresponding to the original and randomized order of trials. We then used these clusters to compare between subjects the brain areas the activity of which can significantly predict the model precision in each of the five ROIs. We identified two clusters of voxels from two subjects as the same if 1) the distance between the boundaries of two clusters was not larger than two voxels in the Talaraich TT_N27 brain space, and 2) both clusters were in the same brain area according to Freesurfer’s cortical labels (aparc.a2009s.annot) from the automatic cortical parcellation.

Performance equated on memory guided pro- and anti-saccades

Figure 1B depicts the behavioral results for both pro and ant-saccade trials. Overall, 85% and 91% of all memory-guided pro-saccade and anti-saccade response times (SRTs), respectively, fell within the normal range of visually guided saccades (i.e., 200–500 ms), and SRT did not significantly differ between the conditions (permutation test, P=0.17). The accuracy of memory-guided pro-saccades and anti-saccades was also indistinguishable (permutation test, P=0.56). Moreover, the final eye position after small corrective saccades (Mackey, Devinsky, Doyle, Golfinos, & Curtis, 2016; Mackey, Devinsky, Doyle, Meager, & Curtis, 2016), but before feedback were not significantly different between the two conditions (permutation test, P=0.53). Thus, we were able to compare neural differences unconfounded by behavioral differences between the pro-saccade and anti-saccade conditions.

Modeled population activity in visual cortex reconstructs a representation of WM

We collected BOLD fMRI signal while the subject performed memory guided pro- and anit-saccade tasks (Figure 1A). In both tasks, each trial can be divided into four epochs; central fixation, stimulation (stimulus encoding), delay period, and response. We mainly focused on stimulus encoding and delay period as these are the intervals that the WM content is being encoded and maintained. We used an IEM (Sprague & Serences, 2013), described in detail in Methods, to reconstruct the spatial representation of neural population activity within each retinotopic map during the encoding and delay periods independently for both pro-saccade and anti-saccade trials (Figure 1C).

In the memory-guided pro-saccade task, the subject is instructed to make a saccade towards the stimulus location after a variable delay period. We hypothesized that the activity patterns of neural population ensembles encode and maintain the stimulus location. Utilizing an IEM, we reconstructed the spatial representation of population activity in the five identified retinotopic areas. To average over all possible stimulus locations, we rotated the estimated channel coefficients, aligning all trials to one location, 45° in the upper right quadrant. Figure 2 shows the reconstruction results for a sample subject (Figure 2A) and for the average of all subjects (Figure 2B) during both encoding and delay epochs. The lines on the polar plot depict the relative channel coefficients across the circular space of the model. The heat maps in the background show the reconstructed spatial representation of the neural population activity, which is the linear combination of the information channels each multiplied by its corresponding coefficients. In all retinotopic areas we identified, the IEM model yielded a high fidelity reconstruction of the stimulus location during both time periods. This suggests that the population activity in these retinotopically organized cortical areas encodes and maintain the stimulus location.

Next, we quantified the accuracy and precision of the model in each retinotopic area by estimating the trial-by-trial performance (Figure 3). For each trial, we estimated the accuracy by comparing the peak of the reconstructed visual space within each retinotopic area with the true stimulus location, separately at the encoding and delay epochs. To do this over all stimulus locations, we regressed the decoded locations, averaged over all trials with the same stimulus location, against the true locations. In all five retinotopic areas, we found a tight coupling between reconstructed location and true location as demonstrated by regression slopes that were near 1 (range across ROIs [0.94, 1.05] with largest P = 1.6×10⁻⁶). The top row of Figure 3 depicts the relation between the decoded location and true location in a sample retinotopic area (V3) during both encoding and delay epochs. We defined the precision of the model as the inverse of the decoded errors across trials, where smaller decoding error approached a precision of 1. The distribution of precision over trials was highly skewed toward 1 for each of the retinotopic areas (range across ROIs [−1.36, −0.51] with largest P = 1.2×10⁻¹³; note negative numbers indicate skewness towards 1). The middle row of Figure 3 depicts the distribution of precision in a sample retinotopic area (V3) during the encoding and delay epochs. Moreover, using various thresholds we asked how often did the model correctly decoded the stimulus location within a confidence interval of 15°, 30°, or 45° of the true location. The bottom row of Figure 3 depicts the percentage of correctly decoded trials at different confidence intervals. For all retinotopic areas, the percentage of correctly decoded trials was significantly above chance even at the lowest interval of 15° (largest P = 2×10⁻¹¹⁸ for all confidence intervals). Overall, our model of the population activity in these retinotopically organized cortical areas accurately and precisely represented the location of the encoded and maintained stimulus location.

Dissociating WM representations of past visual stimulation from future saccade goals

In the memory-guided pro-saccade task, the subject was instructed to plan and make a saccade towards the stimulus location. Therefore, we do not know if the population response encodes the visual target or the saccade goal because they are the same. To dissociate these two possible causes, we used a memory guided anti-saccade task in which the saccade goal is dissociated from the stimulus location (Figure 1A). We used the same IEM as in the pro-saccade task to reconstruct the spatial representation of neural populations in each of the retinotopic areas during encoding and delay time epochs. To avoid any circularity in the reconstruction, we used the data corresponding to the encoding phase of the pro-saccade task to calculate regression weights for the reconstruction at encoding and delay periods of the anti-saccade task.

Figure 4 depicts the reconstructed retinotopic maps for both stimulus encoding and delay epochs for a sample subject and the average over all subjects. As with the prosaccade data, the lines on the polar plot depict the relative channel coefficients across the circular space of the model. Here, we aligned all trials such that the location of the visual stimulus was at 45° and the goal of the antisaccade was at 135°. In all retinotopic areas we identified, the IEM model reconstruction yielded two peaks, one at the location of the stimulus and another at the location of the saccade goal. Remarkably, this suggests that the population activity in these retinotopically organized cortical areas encode and maintain locations that were never visually stimulated, and therefore, must be the result of top-down stimulation.

As with the prosaccade data, we quantified the accuracy and precision of the model in each retinotopic area by estimating the trial-by-trial performance (Figure 5). Focusing solely on the delay period epoch, we compared separately both the decoded versus true location of the encoded stimulus and the decoded versus true location of the antisaccade goal. As depicted in the top row of Figure 5 for sample area V3, the true locations of the stimulus and antisaccade target were well predicted by the reconstructed locations as evidenced by regression slopes that were near 1 (range across ROIs [0.93, 1.12], largest P = 1.0×10⁻⁴). In each of the retinotopic areas, except V2, the distribution of precision over trials was highly skewed toward 1 when decoding the stimulus (range across ROIs [−0.5, −0.24] with largest P = 0.003) and antisaccade target (range across ROIs [−1.0, −0.4] with largest P = 0.009) locations. For V2, the average skewness was −0.4 (P = 0.17). In a sample retinotopic area (V3), the middle row of Figure 5 depicts the distribution of precision of decoding the stimulus and antisaccade target locations based on the pattern of activity in the delay period. Next, we calculated the percentage of correctly decoded stimulus locations and correctly decoded antisaccade target locations within a boundary of 15°, 30°, or 45° of the true location. The bottom row of Figure 5 depicts the percentage of correctly decoded trials at different boundaries. For all retinotopic areas, the percentage of correctly decoded stimulus and antisaccade target locations was significantly above chance even at the lowest interval of 15° (largest P = 8.7×10⁻⁵⁶ for all boundaries). Overall, our model of the population activity in these retinotopically organized cortical areas accurately and precisely represented the location of both the past stimulus location and future antisaccade goal.

Population dynamics track the evolution of the shifting prioritized location

In order to perform the antisaccade task, one must first encode the location of the visual stimulus, then compute, plan and maintain the location of the antisaccade target. As the trial evolves, the prioritized location therefore shifts from the stimulus location to the saccade goal. Consistent with this idea, and as can be seen in Figure 4, it appears as if the reconstruction of the visual stimulus is greater during the encoding epoch than the delay period epoch. It also appears as if the reconstruction of the antisaccade target is greater during the delay than then encoding epoch. Building on this observation, we quantified the gain of the stimulus and saccade target locations in the reconstructed visual space during the encoding and delay period epochs. To this end, we fit 2-D Gaussian functions to the peak locations in the reconstructions at both the stimulus and antisaccade target locations. In Figure 6, we plot the average gain at the stimulus and target locations during the encoding and delay epochs for each retinotopic area. As can be seen, there is a strongly significant reconstruction gain interaction between epoch (encoding/delay) and location (stimulus/saccade target) in all retinotopic areas (two-way ANOVA for V1, V2, V3, V3AB, and IPS0: smallest F(1,99) = 597, P = 10⁻⁸¹). This pattern supports the hypothesis that the spatiotemporal pattern of the population activity tracks the prioritized location as it evolves from the visual stimulus to the antisaccade goal.

Delay period activity in frontal and parietal cortex sculpts the population activity in early retinotopic areas improving the precision of readout

We demonstrated that the pattern of population activity in visual cortex represents the antisaccade goal. This representation is not the result of bottom-up visual stimulation, and therefore, must be instantiated by top-down feedback signals. In order to identify the source of these top-down signals, we searched for brain areas whose trial-by-trial fluctuations in the magnitude of delay period activity predicted the precision with which our model could accurately reconstruct the remembered location. We first, for each retinotopic area, created a vector of values that represented the model accuracy, inversely proportional to the error of the reconstructed WM location. Second, using a generalized linear model (GLM) we estimated, for all voxels in the brain, the delay period activity for each trial. Then, we calculated the correlation between decoding accuracy and delay period activity of each voxel. In individual subjects, we found clusters of voxels whose delay period activity significantly predicted the decoding accuracy of the antisaccade goal from the pattern of activity in area V3 (Figure 7A). In the sample subject shown, we find significant clusters in dorsolateral prefrontal cortex, superior precentral sulcus, dorsal parietal cortex, and lateral occipital cortex. Despite some individual differences in the exact locations of significant voxel clusters within brain areas, we found consistent overlap in these four brain areas (Figure 7B). The histograms in Figure 7 depict the distribution of voxel correlations within each of the four brain areas, as well as their significance thresholds determined by permutation testing (Methods). In addition, Table 1 contains a list of all clusters of voxels in the brain, containing more than four contiguous voxels with overlap in at least four out of five subjects, whose activity is significantly correlated with the decoding accuracy in areas V2 and V3. We focus here on areas V2 and V3 because the correlations were particularly strong and consistent, compared to the correlations with decoding from areas V1, V3AB and IPS0, which was less consistent. The patterns of correlations identify both attentional and motor planning networks (Corbetta & Shulman, 2002) as potential sources for the modulatory top-down control signals prioritizing the saccade target location in the population activity in early visual cortex.