Comparison of randomized multifocal mapping and temporal phase mapping of visual cortex for clinical use.

fMRI is becoming an important clinical tool for planning and guidance of surgery to treat brain tumors, arteriovenous malformations, and epileptic foci. For visual cortex mapping, the most popular paradigm by far is temporal phase mapping, although random multifocal stimulation paradigms have drawn increased attention due to their ability to identify complex response fields and their random properties. In this study we directly compared temporal phase and multifocal vision mapping paradigms with respect to clinically relevant factors including: time efficiency, mapping completeness, and the effects of noise. Randomized, multifocal mapping accurately decomposed the response of single voxels to multiple stimulus locations and made correct retinotopic assignments as noise levels increased despite decreasing sensitivity. Also, multifocal mapping became less efficient as the number of stimulus segments (locations) increased from 13 to 25 to 49 and when duty cycle was increased from 25% to 50%. Phase mapping, on the other hand, activated more extrastriate visual areas, was more time efficient in achieving statistically significant responses, and had better sensitivity as noise increased, though with an increase in systematic retinotopic mis-assignments. Overall, temporal phase mapping is likely to be a better choice for routine clinical applications though random multifocal mapping may offer some unique advantages for selected applications.

Introduction

Functional magnetic resonance imaging (fMRI) is being used increasingly for mapping key brain structures prior to surgical treatment of tumors and other types of focal pathology. In such applications, the goal is to identify viable areas of the brain that might be at risk of damage due to resection, radiation or other invasive treatment. For tumors of the occipital lobe or adjacent portions of the parietal or temporal lobes, damage to visual cortex can cause partial or complete blindness or other disruptions of visual perception (Martin et al., 2012).

Temporal phase mapping

Conventionally, cortical maps of the visual field are charted using a temporal phase mapping technique that consists of a rotating checkered wedge, or an expanding checkered ring stimulus (Fig. 1A, B) (DeYoe et al., 1996; Engel et al., 1994, 1997; Sereno et al., 1995). Temporal phase mapping can evoke robust responses in primary visual cortex (V1) and a number of extrastriate visual areas. Moreover, the checkered ring stimulus can identify cortical locations that support central vision, which if damaged can impair reading and other critical visual functions. The ability to identify the representation of central vision makes temporal phase mapping particularly useful for pre-surgical planning and generally superior to simple flashed checkerboards or pulsed lights (DeYoe et al., 2011). It is also time efficient in that all eccentricities or polar angles throughout the visual field can be mapped in less than 4 min. Temporal phase mapping has been used in a variety of clinical disorders, including inherited photoreceptor abnormalities (Baseler et al., 2002), amblyopia (Conner et al., 2007), glaucoma (Duncan et al., 2007a), albinism (Hoffmann et al., 2003), achiasma (Hoffmann et al., 2012; Sinha and Meng, 2012), scotoma (Sunness et al., 2004), long-period deprivation of visual input (Levin et al., 2010) and developmental reorganization of cortical visual field maps (Muckli et al., 2009). It has also played an important role in the ongoing debate over cortical plasticity (Baseler et al., 2011).

Fig. 1

Visual stimuli. A. Phase-based eccentricity mapping. Left to right: outlines of elementary ring steps; sample stimulus image composed of 4 elementary rings at one time point. B. Phase-based polar angle mapping. Left to right: outlines of 20 elementary wedge steps; sample stimulus image composed of 5 elementary wedges at one time point; stimulation sequences for three visual field locations indicated by red dots. C–E. Random multifocal stimuli for 13, 25 and 49 segments. D. (right) Three stimulation sequences associated with three visual field locations indicated by red dots. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

As illustrated in Fig. 1A and B, the fMRI signals produced by temporal phase mapping are periodic waveforms that are distinguished from each other by their temporal phases. As the wedge/ring sweeps through each visual field location, fMRI activation sweeps through retinotopically corresponding locations in visual cortex. The timing of the activation at a particular brain voxel is determined by the distance of a voxel's “population receptive field” (pRF) (Dumoulin and Wandell, 2008) from the stimulus onset location plus a delay due to the sluggish hemodynamic mechanism responsible for the fMRI signal. Consequently, the ability to precisely locate the pRF of a given voxel is limited by the accuracy of the estimation of the temporal properties of the signal and the variability of the local hemodynamics. Computing the temporal phase of the fMRI signal is typically accomplished by cross-correlation with a reference waveform (e.g. sinusoid) (Bandettini et al., 1993; Saad et al., 2003) or through Fourier analysis (Boynton et al., 1996; Engel et al., 1997). Removing the additional phase delay caused by the BOLD hemodynamics can be achieved with fMRI scans using two stimuli moving in opposite directions, such as clockwise versus counterclockwise rotating wedges or expanding versus contracting rings (Sereno et al., 1995). For each voxel, the phases obtained with stimuli moving in opposite directions can be averaged to cancel out the phase shift caused by the hemodynamic delay. The resulting corrected time series can then be averaged across repetitions for improved signal-to-noise ratio (SNR).

The dependence of temporal phase mapping on precise timing of the fMRI signal makes it susceptible to signal distortions. Though correlation methods often used with temporal phase mapping can provide good immunity to some types of pulse or burst noise, other types of noise can blur the small phase difference of responses evoked by adjacent visual field locations, thus introducing errors in preferred stimulus location. In addition, if a voxel contains a mixture of neurons with spatially distinct receptive fields, as can occur for a voxel straddling a sulcus, then the voxel's response will be a sum of multiple (approximate) sinusoids having the same period but different phases. Such a sum, if non-zero, results in a single sinusoid with an erroneous phase that is intermediate between those of the true individual components. As a result, the estimated preferred stimulus location for that voxel will also be in error.

Randomized multifocal mapping

A potential solution to the problem of temporal distortion is to use a code-based mapping paradigm such as randomized multifocal stimulation. Fig. 1C–E illustrates multifocal visual stimuli consisting of multiple checkered segments that are each presented in a unique randomized temporal pattern. Multifocal stimulus paradigms have been used to explore visual cortex using a variety of neurophysiological techniques including visual evoked potentials (VEP) (Baseler et al., 1994; Slotnick et al., 1999), magnetoencephalography (MEG) (Tabuchi et al., 2002), and fMRI (Hansen et al., 2004; Vanni et al., 2005). In this approach, the stimulus segments to which a voxel responds can be identified by the unique ON/OFF pattern of the fMRI response. In other words, each stimulus sequence can be viewed as a unique temporal code for a specific region of the visual field. To identify which stimulus segments activate a particular voxel, a conventional multiple regression analysis can be employed with the unique time series of each stimulus segment used as a regressor (Ward, 2006).

Multifocal mapping has been tested in clinical applications including post-surgical mapping of primary visual cortex (Vuori et al., 2012), and the measurement of training-induced changes in the cortical representation of a hemianopic field (Henriksson et al., 2007). It has also been suggested that multifocal mapping might save time compared to the individual presentation of multiple static stimuli which could be beneficial for scanning elderly glaucoma patients (Duncan et al., 2007b).

One potentially important advantage of multifocal mapping is its ability to correctly identify voxels that respond to multiple separate locations in the visual field. Moreover, the unique ON/OFF ‘digital code’ associated with each stimulus segment potentially makes multifocal mapping more tolerant of temporal distortions. Indeed, computational simulations indicate that a random multifocal stimulus paradigm can provide more accurate retinotopic assignments in the face of noise and hemodynamic variations (Ward et al., 2012).

In this paper we describe a randomized multifocal mapping technique that permits more flexible presentation sequences than popular maximum-length-shift-register sequences (m-sequences) (Buracas and Boynton, 2002; Sutter, 2001), and compare it with conventional temporal phase mapping on a number of clinically relevant measures including effectiveness, efficiency and reaction to increased noise level. The use of more flexible presentation sequences permits greater control over scan duration, which can be advantageous for clinical applications. This includes real-time brain mapping paradigms in which scan length is adjusted adaptively based on recursive estimation of data quality (Ward et al., 2012).

Material and methods

Subjects

Seven subjects (three female), 18–35 years old, participated in this study. All subjects had normal or corrected-to-normal visual acuity and provided informed, written consent in accordance with methods approved by the Institutional Review Boards of both the Medical College of Wisconsin and Marquette University.

Visual stimulation

Visual stimulation was achieved using a custom back-projection system providing a 20° radius field of view. Custom computer graphical images were generated using a Cambridge Research System ViSaGe visual stimulus generator and displayed via a Sharp XG-C330X color video projector. All stimuli consisted of 8 Hz flickering checkerboard patterns with black/white contrast of 96% superimposed on a uniform gray background of 265 cd/m2. Luminance levels for black and white checks were 23 and 1120 cd/m2 respectively. Subjects were asked to passively fixate on a green marker of 0.1° radius of the center of the display throughout the stimulus presentation.

The randomized multifocal stimulus (Fig. 1C–E) consisted of an array of segments arranged in concentric rings so that each segment covered a predetermined range of eccentricity and polar angle thereby allowing both dimensions to be mapped simultaneously. As shown in Fig. 1C–E, we tested three randomized stimulus patterns consisting of 13 (3 rings × 4 wedges + 1 center segment), 25 (3 rings × 8 wedges + 1 center segment), and 49 segments (3 rings × 16 wedges + 1 center segment). The 3 rings of segments extended radially from 0.32° to 20° visual angle, with borders at eccentricities of 3.2° and 8.55° visual angle. The rings were scaled in size based on the retino-cortical mapping function reported by Schira et al. (2007, 2010), so as to stimulate roughly comparable-sized patches of primary visual cortex.

Fig. 1 also shows single, 2 s stimulus frames and stimulus timing waveforms for 3 locations in the visual field. For the random stimuli, individual segments were turned ON and OFF during the fMRI scan using unique, 600 s random presentation sequences, three of which are illustrated at the right of Fig. 1. Typically, one unique random sequence was repeated twice with each subject and the resulting fMRI time series were concatenated to provide a single 1200 s series. Within each stimulus sequence, ON block durations varied from a minimum of 10 s to 40 s and intervening OFF blocks varied from 2 to 88 s while maintaining either a 25% or 50% duty cycle overall depending on the experiment. Specific sets of sequences were selected from 1000 randomly generated sets so as to minimize the error in estimating the fMRI response attributable to each stimulus segment, in other words, to maximize estimation accuracy. Specifically, a covariance matrix was calculated for each set of random stimulation sequences after convolution with a standard hemodynamic response function (Cohen, 1997). The covariance matrices were then compared to identify one stimulus set with the smallest, maximum covariance.

For comparison, we used conventional temporal phase mapping stimuli consisting of an expanding ring or rotating wedge (Fig. 1A, B). To map visual field eccentricity from 0.4° to 20° radius, this range was divided into 16 elementary rings scaled according to the formula Rmin = Rmax ∗ f, where Rmin/Rmax was the innermost/outermost radius — 0.4°/20° (vis. ang.) and N was the number of the elementary rings. At the beginning of each run, the stimulus consisted of just the center segment (area within 0.4°) and then recruited one more elementary ring each 2 s until the center segment plus 3 elementary rings were included. Then the combined ring shifted outward by one elementary ring every 2 s until it slowly disappeared in steps at the outer margin of the display. The display was uniform gray for 2 s then the sequence repeated from the center. Each full cycle contained 20 steps (40 s), and each scan had 20 cycles, thus spanning a total of 800 s. To map polar angle, a checkered wedge spanning a quarter of the visual field rotated counter-clockwise in 20 steps (40 s) beginning centered on the lower vertical meridian. The rotation of the wedge was also repeated 20 times per scan for a total of 800 s.

Image acquisition

Brain images were obtained with a 3 T General Electric Signa Excite 3.0 MRI system equipped with an 8 channel High Resolution coil. A 96 × 96 voxel matrix covering a 24 × 24 cm field of view and 25 slices was used to obtain 2.5 mm3 voxels. A vendor-supplied, gradient-recalled EPI pulse sequence was used for fMRI imaging. The MR parameters were as follows: 77° flip angle, TR of 2000 ms, and TE of 30 ms. The first 4 images of each scan, acquired while a uniform gray stimulus was displayed, were discarded to allow brain tissue magnetization to achieve steady state. A high resolution, T1-weighted, spoiled GRASS (gradient-recalled at steady state) anatomical image was also collected. The anatomical data set covered the whole brain with a voxel size of 0.9375 mm × 1.0714 mm × 1.0 mm (flip angle = 12°, TR = 8.2 ms, TE = 3.2 ms). Raw MRI data were converted into image format using GE Signa software and assembled into a time series of volumetric imaging data using the AFNI analysis package (Cox, 1996). Head motion correction was applied using AFNI's 3dvolreg routine, which registered each volume of images at each time point to the base volume. The base volume was selected to be the 50th volume of the functional data set that was collected closest to the collection of SPGR (Detailed description of each AFNI function used in this paper is available at AFNI website http://afni.nimh.nih.gov/pub/dist/doc/program_help/). For the image in the bottom half of Fig. 3, the fMRI time courses were spatially filtered with a 3.5 mm averaging kernel. fMRI data used in all other analyses were not filtered spatially or temporally.

Fig. 3

Retinotopic maps from one subject on flat brain surfaces. Black outlines: Visual area boundaries as defined in Fig. 2. White curve: border between V3 and higher visual areas. Top: Maps without spatial or temporal smoothing. Bottom: Maps after a spatial smoothing with a 3.5 mm spherical kernel applied to volumetric fMRI time series data. Note that VO1/VO2 on the left hemisphere, denoted by dashed outlines, was drawn with uncertainty based on the limited information provided by phase-encoded retinotopic maps in that region.

Regression analysis of random mapping data

For the random mapping data, linear regression was used to detect activated voxels and determine the stimulus segments to which they responded. This was accomplished using AFNI's 3dDeconvolve (Cox, 1996; Ward, 2006), which performs multiple linear regression using multiple input stimulus time sequences. In the multiple linear regression model for random mapping data, regressors consisted of the individual time sequence for each stimulus segment convolved with the default AFNI hemodynamic response function (HRF). The default HRF was a gamma variate function HRF(t) = t8.6 * e− t/0.547 (Cohen, 1997), where t represents time in seconds. A nuisance regressor consisting of a 3rd order polynomial was also included, to model baseline and very low frequency drifts. The regression analysis generated a least squares estimate of the fit coefficients and a partial F-statistic for each voxel with respect to each random segment's regressor. The goodness of fit for the whole model with all regressors was quantified by a full F-statistic. For each voxel, the random stimulus segment with the largest partial F-statistic was recorded as the ‘best segment’ for that voxel and the ring and wedge comprising this segment were recorded as ‘best ring’ and ‘best wedge’. Both the full and largest partial F-statistics for each voxel were converted to z-scores having equivalent statistical p-values to facilitate comparison with the temporal phase mapping results.

Correlation analysis of phase mapping data

The fMRI responses generated by phase mapping were evaluated using the Hilbert Delay algorithm of the AFNI package (Saad et al., 2001). This algorithm yielded the maximum cross-correlation between each voxel's BOLD response and a reference sinusoidal waveform, and provided the optimal time lag to achieve maximum cross-correlation. The correlation coefficient at the optimal time lag of each voxel was converted to a z-score having the same p-value as the original correlation coefficient.

As mentioned in Section 1.1, the temporal phase shift of a voxel's BOLD response is determined by its preferred location in the visual field plus an additional phase shift caused by the hemodynamics. Phase shifts for all activated voxels were corrected to account for this hemodynamic component prior to determining the preferred stimulus location for each voxel. To make this correction, the preferred stimulus eccentricity and angle for each voxel were first estimated from the uncorrected angle and eccentricity phase data. A symbol (circle) for each voxel was then plotted at the corresponding eccentricity and angle on a diagram of the visual field. This was repeated for all active voxels in medial occipital cortex of each hemisphere. Since medial occipital voxels from one hemisphere primarily respond to the contralateral visual field, the temporal phase measurements were then adjusted by a constant factor so that the symbols from each hemisphere were symmetrically distributed with respect to the vertical meridian. A comparable adjustment was also made to the eccentricity phase data and cross-checked to ensure that V1 voxels activated by foveal stimuli were at or near the occipital pole. The corrected polar angle and eccentricity phase data were then used to identify the final preferred stimulus location for each voxel.

Visual area delineation

To view retinotopic maps, a 3D cortical surface model of each subject's hemisphere was reconstructed from the high resolution anatomical data using Caret software (Van Essen et al., 2001), as shown in Fig. 2A. The surface model was then computationally unfolded to yield a 2D flat map of the cortex (Fig. 2B). ROIs of visual areas were delineated based on consensus of two raters through visual inspection of the phase-encoded retinotopic maps. Specifically, the flat maps were first colored by voxels' phase shifts from polar angle mapping. Boundaries were drawn at the place where the sequential color pattern reversed itself, which indicated a phase reversal. The boundary of each visual area was closed to include a complete hemi or quarter field representation, and cross checked with eccentricity phase maps to ensure that each visual area contained a foveal representation. Depending on the quality of the retinotopy for each subject, up to nine complete visual areas, V1, V2, V3, V3AB, hV4, hMT +, V7, VO1, VO2, were identified using published criteria (Arcaro et al., 2009; Hansen et al., 2007; Harvey and Dumoulin, 2011; Swisher et al., 2007; Wandell et al., 2007; Winawer et al., 2010). A set of these 9 areas for one subject is shown as colored patches on the cortical surface models in Fig. 2. In the worst case, three complete visual areas and five half visual areas (only one hemisphere) were identified.

Fig. 2

Visual area ROIs. A. 3D cortical surface models with legend indicating colors assigned to different visual areas. B. 2D flat maps of the two cerebral hemispheres of one subject. D: Dorsal, M: Medial, A: Anterior, P: Posterior. Visual area names according to Swisher et al. (2007) and Arcaro et al. (2009). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Effects of noise

To test the effects of noise on the ability to correctly classify each voxel's preferred stimulus location, we simulated two new experiments with 20 stimulating wedges throughout the visual field under two assumptions: (1) the BOLD signal is linear (Boynton et al., 1996; Hansen et al., 2004) and (2) the noise is white. Although fMRI responses can be non-linear under certain conditions, it is reasonably linear as long as the stimulation duration is longer than about 4 s (Gu et al., 2005), which is true for most retinotopic mapping experiments. It is known that non-white noise is present in the BOLD signal (Lund et al., 2006). However, these differences will likely be quantitative rather than qualitative. The goal of our analysis was to demonstrate that there was a qualitative difference in how noise affects the retinotopic classification of phase and random mapping.

To demonstrate the differential effects of noise on the two techniques, we wanted all conditions to be constant other than the addition of noise. This can be accomplished most precisely in the context of a pure simulation. For the simulated phase mapping, the 20 wedges were turned ON/OFF using periodic sequences cyclically shifted in time, simulating the temporal stimulating sequences of phase mapping. For the simulated random mapping, the 20 wedges were turned ON/OFF using random sequences, replicating the temporal stimulating sequences of random multifocal mapping.

The ON/OFF sequence from one wedge was used to simulate 100 fMRI time series at SNRs of 5 and 0.5. To do this the selected ON/OFF sequence was convolved with an estimate of the HRF. Next, 100 different series of normally distributed pseudorandom numbers generated by the “randn” function of MATLAB were scaled to a power determined by the desired SNR and added to the sequence generated in the previous step. Finally, the same regression analysis used to analyze the empirical data was used to analyze the simulated time series. The 20 random ON/OFF sequences associated with the stimulus wedges were regressors. From the regression analysis, the wedge whose ON/OFF sequence best matched each simulated time series was identified. A histogram of the number of simulated time series that were assigned to each wedge location was then constructed.

For phase mapping, 20 phase-shifted periodic square waves having the same length as the sequences for random wedges and a period of 40 s were used to simulate the ON/OFF sequences associated with the 20 wedge locations. One hundred fMRI time series were then simulated in the same way as for the random mapping. The same correlation analysis used to analyze the empirical phase-based fMRI data was used to analyze the simulated time series. The resulting delay values were grouped into 2 s bins to generate an equivalent histogram to that prepared for random mapping. For both paradigms at both SNR levels, the response measures (F-statistic, correlation coefficient) were thresholded at an equivalent p-value (7.1e-05) to identify significant responses included in the histogram.

Effects of duty cycle

Duty cycle is an important design parameter for random multifocal stimulation because it controls the number of segments that, on average, are ON or OFF at the same time, which subsequently determines the amount of modulation for a voxel that is stimulated by more than one segment. Therefore we also performed an experiment comparing random 25-segment multifocal mapping with 50% duty cycle, 25% duty cycle and with the phase-mapped polar angle mapping, which produces a 25% duty cycle as described above in Section 2.2. Data acquisition and analysis followed the procedures described in Sections 2.2–2.5. This experiment was performed with two subjects.

Results

We compared random multifocal mapping and temporal phase mapping with respect to several factors that are of particular relevance for clinical use: effectiveness in rendering retinotopy, time efficiency and reaction to increased noise level. We also explored the effect of duty cycle on random multifocal mapping.

Effectiveness

Rendition of retinotopic features

In a clinical context, retinotopic maps assist the physician in understanding the relationship between the functional organization of visual cortex and the integrity of visual perception across the patient's visual field. As illustrated in the top half of Fig. 3 for one subject, both phase mapping and random mapping generated clear retinotopic maps within V1, V2 and V3 (outer boundary of V3 is marked by white curve). Retinotopic maps for extrastriate areas beyond V3 were more homogeneous for phase mapping than for random mapping. Phase mapping revealed at least 6 more visual areas than multifocal mapping including V3AB, hV4, V7, VO1, VO2, and hMT +. This trend was evident in all 5 subjects. (The boundaries of VO1/VO2 in the left hemispheres of Fig. 3 were drawn with some uncertainty based on the limited information provided by the phase-encoded retinotopic maps in that region, so are denoted by dashed outlines.) The global features of all retinotopic maps were accentuated by spatial filtering of the volumetric time course data with a spherical averaging kernel of 3.5 mm radius, as shown in the bottom half of Fig. 3, though this is at the expense of some loss of local detail, and the potential introduction of retinotopic distortions (e.g. compare the unfiltered vs filtered polar angle maps for area MT in the left hemisphere).1 Even with smoothing, the extrastriate retinotopic map obtained with random stimulation still appears less homogeneous than for phase-mapping.

Consistency of voxel-based retinotopy

Overall, multifocal- and phase-mapping tend to yield comparable retinotopic assignments for each voxel in V1–V3. This can be appreciated qualitatively by noting the general similarity of the retinotopic maps within V1–V3 in Fig. 3. To assess this quantitatively, we first identified groups of voxels that responded most strongly to each of the different random stimulus segments. We pooled voxel groups preferring segments at the same polar angle. For each pooled group, we then constructed a histogram of the preferred phase mapped angles. Fig. 4A shows the histograms for each pooled group with the angular range of the corresponding random mapping segments marked by the shading in the background. All histograms shown in Fig. 4 were averaged across five subjects. (This comparison could only be conducted for V1–V3 since random multifocal mapping did not evoke meaningful responses in visual areas beyond V3.) Note that for V1 each phase mapped histogram is roughly symmetrical and overlaps the range of angles spanned by the corresponding group of random stimulus segments, thereby indicating good consistency between the random and phase-mapping retinotopic assignments. The uneven heights of the histograms reflect differences in the number of voxels representing each angular location consistent with other reports from this lab (Janik et al., 2003). For V2 and V3, the correspondence between each histogram and each background shading was roughly conserved but with random rightward or leftward shifts of some histogram peaks, indicating some degradation in the consistency of retinotopic assignments of random multifocal mapping and phase mapping at higher levels of the cortical visual system.

Fig. 4

Correspondence of retinotopic assignments in V1, V2 and V3 between phase and random mapping. A. Correspondence between preferred polar angle for phase mapping (line graphs) and preferred random wedge (color shading). Averaged for five subjects. B. Average phase mapped, polar angle distribution (line graph) vs width of a single random wedge (gray shading). C. Correspondence between preferred eccentricity for phase-mapping (line graphs) and preferred random ring. Note: Ring widths were scaled in proportion to eccentricity. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4B shows the average of the individual histograms in 4A after magnitude normalization and shifting to a common center. The shading in Fig. 4B shows the 45 degree range (± 22.5°) that is spanned by each group of random mapping segments at the same polar angle. It is evident that a single random mapping wedge preferentially activates voxels whose preferred phase-mapped angle extends beyond the range spanned by the segment itself. The extents of this “over-activation” are similar in V1, V2 and V3, indicated by the similar widths of the three normalized and averaged histograms.

The same approach was used to create the eccentricity-related histograms of Fig. 4C. Here the widths of the three histograms in each graph appear to vary because the eccentricity ranges of the corresponding random rings vary with eccentricity in order to roughly compensate for changes in cortical magnification with eccentricity. It is also obvious that the phase-mapped histogram peaks are skewed toward greater eccentricities relative to the corresponding random stimulus rings. This effect is partly accounted for by the non-linear scaling of the stimulus rings but also may reflect some bias in the calibration of the temporal phase data.

Non-uniform representation of visual field locations

Previously, we have shown that temporal phase-based retinotopic maps of V1 often show a significant under-representation of voxels responding preferentially to visual field locations near the vertical meridian as compared to the horizontal meridian (Janik et al., 2003). In this study we revisited this issue using both phase and random mapping. For both methods, the numbers of voxels responding best to different polar angle locations were non-uniform for V1 but less so for V2 and V3 (Fig. 5, averaged for five subjects). More voxels represented the horizontal meridians than the vertical meridian (note color code in Fig. 5).

Fig. 5

Number of voxels having preferred polar angles for (A) random mapping and (B) phase mapping. Averaged for five subjects. Inset: Color code for preferred angle. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Efficiency

For clinical patients who may have attention, alertness or cognitive difficulties, time efficiency can be a critical factor in determining the success or failure of an fMRI exam. Consequently, we inspected the efficiency of the two mapping paradigms with respect to the duration of data acquisition needed to: 1) determine the responsiveness/non-responsiveness of each voxel and 2) stabilize the number of voxels responding to each stimulus location.

Responsiveness of visual areas

For the random mapping paradigm, the F-statistic for the full regression model provides an overall indicator of a voxel's ability to respond significantly to one or more of the stimulus segments. The temporal correlation coefficient provides an analogous statistic for phase mapping. To permit a quantitative comparison between the F- and correlation statistics, we converted both to equivalent z-scores. Fig. 6A shows the mean and the upper and lower quartiles of the z-scores for all voxels within each cortical visual area averaged across five subjects as a function of data acquisition time for phase-mapped polar angle mapping and random 25-segment mapping. Curves for polar angle mapping are colored black while those for random 25-segment mapping are colored red. In all visual areas, the z-scores for polar angle mapping increased faster with acquisition time than those for random mapping. In other words, phase mapping was more efficient with respect to the data acquisition time needed to achieve a specified level of statistical significance. This was confirmed by a 3-way ANOVA that showed significant main effects (p < 0.05) for mapping method (random vs polar angle), acquisition time and visual area as well as significant interaction effects (p < 0.05) for method vs acquisition time. (Interaction tests were only computed using data at time points common to both polar angle and random mapping, i.e. 200, 300 and 400 s.) Other interaction effects were not significant. The graphs of Fig. 6A also show that random mapping was not capable of strongly activating visual areas other than V1, V2, and V3 even at the longest scan duration of 1200 s. This is consistent with the qualitative impression obtained from the cortical maps in Fig. 3. One thing to note is that the phase-mapped eccentricity mapping was always more efficient than polar angle mapping (data not shown).

Fig. 6

Time efficiency represented by z-scores. A. Time efficiency for random vs phase mapping in different visual areas. B. Time efficiency for random stimuli with 13, 25 or 49 segments. Solid lines: mean z-score. Error bar = standard error of the mean. Dashed lines: Lower and upper quartiles. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 6B shows the averaged mean and upper and lower quartiles as a function of acquisition time for random mapping with different numbers of stimulus segments (13, 25, and 49). ANOVA showed significant main effects (p > 0.05) for number of stimulus segments, acquisition time and visual area as well as significant interaction effects between all three pairs of factors. It is also noteworthy that in V1–V3, z-scores improved slightly faster (steeper slope) for smaller numbers of segments. In visual areas beyond V3, z-scores didn't increase much for any number of segments, but they tended to be higher for larger rather than smaller numbers of segments (Note the relative vertical order of the different colored curves for V1–V3 versus the other areas.)

We repeated this same analysis using the largest partial F-statistic for each voxel rather than the full F but obtained virtually identical results suggesting that for the vast majority of voxels, the response may be primarily driven by a single random segment. Since both our qualitative and quantitative assessments indicated that random mapping as performed here was not effective beyond V3 for all five subjects, all subsequent comparisons focus on data just from V1–V3.

Stabilization of retinotopic assignment

Another measure of fMRI mapping efficiency is the time needed to reach a stable assignment of retinotopic preference for each voxel. This measure goes beyond simply identifying a voxel as responsive. The focus here is on the time needed to stabilize the retinotopic assignment of each voxel. To assess this factor, we tracked the number of voxels2 responding to each stimulus location (random segment or phase-mapped polar angle location) as a function of data acquisition time. In Fig. 7A and B, the top row graphs show the increase in the total number of voxels identified as responsive in visual areas V1–V3 versus data acquisition duration. Graphs in the second rows track the number of voxels responding preferentially to each of the 25 stimulus segments or the 20 phase-mapped polar angle positions (each segment or position is represented by a different colored line). Finally, the third rows of graphs depict rate-of-change (1st derivative) of the voxel counts. We defined the time-to-stabilize as the time point when the mean unsigned rate-of-change first fell below 20% of the peak rate-of-change (Fig. 7, red, blue dashed lines). In panel A, the graphs show that for random mapping the numbers of voxels responding to each stimulus segment tend to stabilize within approximately 300 s for V1, 360 s for V2, and 510 s for V3. However, the ‘settling’ behavior for phase-mapped polar angle mapping is noticeably different, as shown in panel B. The total number of voxels identified as active is almost constant after the first 1.5 cycles of stimulus presentation (60 s). But, the numbers of voxels preferring each polar angle position vary widely during the first 120 s for all three visual areas. On the whole, voxel counts for phase-mapped polar angle mapping were less stable than for random mapping as best appreciated in the rate of change graphs. After stabilization, the mean unsigned rate-of-change was 4 times larger for polar angle mapping compared to random mapping.

Fig. 7

Time evolution of voxel activation for (A) random vs (B) phase mapping in visual areas V1, V2, V3, A1, and B1: Total number of activated voxels for all visual field locations. A2 and A3: Number of activated voxels for each individual visual field location (separate curve for each). A3 and B3: Rate of change of voxel counts for each stimulus location (separate curve for each). Dotted lines: Stabilization times. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

As described in detail in the Material and methods section, we used a simulation to test the effects of noise on the retinotopic assignments obtained with random multifocal versus phase mapping. The retinotopic assignments of the two paradigms for 100 simulated time series under two different SNR levels are shown as histograms in Fig. 8. When SNR = 5, all simulated time series passed the statistical threshold of p < 7.108e-05 and were correctly assigned to the stimulating wedge by both paradigms. When SNR = 0.5, phase mapping mis-assigned 23% of time series to the neighboring wedges even though all 100 time series passed the same threshold. Importantly, for random multifocal mapping, the addition of random noise does not cause the time series to be assigned to an incorrect location. Rather, 5% of the time series fail to reach statistical significance thus dropping out of the analysis, in effect, reflecting reduced sensitivity.

Fig. 8

Effect of noise on preferred wedge location for simulated random and phase mapped responses to a single wedge location (#10). Histograms include only responses that are statistically significant (p < 7.108e-05) A. SNR = 5. B. SNR = 0.5.

Effect of duty cycle on the multifocal paradigm

Fig. 9 illustrates the effect of duty cycle on the z-scores of random multifocal mapping averaged across two subjects. In V1, z-scores of random multifocal mapping with a 25% duty cycle (blue curve) increase slightly faster than for a 50% duty cycle (red curve). This difference was increasingly apparent in V2, V3, V3AB, and hV4. However, phase-mapped polar angle mapping was still more efficient than random mapping at either duty cycle. A 3-way ANOVA showed significant effects (p < 0.05) for method (25% duty cycle, 50% duty cycle, and phase), visual area and acquisition time. Interaction effects between all three pairs of factors were also significant (p < 0.05).

Fig. 9

Effect of duty cycle on time efficiency of random stimulation across visual areas. Solid lines: Mean z-score. Error bar = standard error of the mean. Dashed lines: Lower and upper quartiles. Polar angle phase mapped curves shown for reference. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Discussion

We compared a random, multifocal visual stimulation paradigm and associated GLM analysis with conventional temporal phase mapping and correlation analysis, both common paradigms. Random multifocal mapping generated good quality retinotopic maps only in V1–V3, while phase mapping rendered retinotopic features in three to six more extrastriate visual areas for different subjects. The quantitative analysis of averaged z-scores in each visual area confirmed this observation. Within V1–V3, retinotopic assignments obtained with both methods were consistent for the majority of voxels. Phase mapping was more efficient with respect to both the acquisition time needed to achieve a statistically significant response across voxels, and to achieve stable retinotopic assignments of those voxels, though the asymptotic stability of those assignments was poorer than for random multifocal mapping. With increased noise levels, phase mapping detected the same number of significantly responding voxels but allowed small errors in retinotopic assignment. In contrast, random multifocal mapping kept all correct retinotopic assignments but with a slightly reduced number of significant voxels. Decreasing the duty cycle of the random stimuli marginally improved their efficiency, but the multifocal paradigm still was not as efficient as phase mapping.

Randomized stimulation sequences

Previously m-sequences and quadratic residue sequences have been proposed for fMRI-based visual cortex mapping (Buracas and Boynton, 2002; Vanni et al., 2005). Mathematically, m-sequences are considered to be maximally efficient (Buracas and Boynton, 2002). However, m-sequences are based on rigid sequence lengths (powers of 2 minus 1) and other constraints that are not necessarily optimal for clinical applications. Therefore, our approach was to use a more flexible method and focus on selecting a unique set of random sequences that maximally reduces the variance of the fit coefficients in our linear regression model, thereby reducing the error in the retinotopic maps. This approach allows the scan length and number of stimulus segments to be flexible, which can then be optimized as most appropriate for a particular clinical application where minimum MR scan time and assurance of scan success may be more important than highly detailed rendition of retinotopic features other than the identification of foveal vision (DeYoe et al., 2011).

Insufficient activation by random mapping

The generally higher z-scores obtained with temporal phase mapping compared to random mapping may reflect an advantage due to the retinotopic organization of visual cortex which may favor focal, spatially sequential stimulation. Perhaps response “recruitment” within relatively large zones of cortex activated by phase mapping stimuli yields enhanced fMRI responses compared to smaller distributed foci of activation evoked by the random multifocal stimuli. Such “recruitment” could be either neural or hemodynamic in origin. In addition, visual attention focused on a single phase mapping wedge or ring may enhance the cortical response compared to random mapping in which multiple segments spread throughout the visual field provide no consistent target for focal attention.

One consistent observation was that random mapping failed to activate extrastriate visual areas beyond V3 that were well activated by temporal phase mapping. This effect may be related to the size of the random stimulus segments relative to the sizes of voxel pRFs in different cortical visual areas. Small pRFs such as those in V1 will often be contained entirely or largely within the retinotopic domain of a single random stimulus segment as depicted in Fig. 10B (top). As a result, the voxel will be fully modulated as the segment turns ON and OFF (Fig. 10A middle waveform). In contrast, a voxel in, say, V4 or VO can have a pRF that covers portions of several adjacent stimulus segments as depicted in Fig. 10B (bottom). Since each segment is turning ON and OFF randomly, the pRF may rarely be fully stimulated or fully unstimulated thereby leading to poorer modulation as illustrated in Fig. 10A (bottom waveform). This effect will be exacerbated as the size of the pRF increases and the size of the stimulus segments becomes smaller. We observed both such effects in our empirical data.

Fig. 10

Effect of pRF size on response modulation. A. Simulated pRF responses to segments 17, 18 or both. B. Simulated pRF's (white outlines) relative to random segment stimuli (black outlines) color coded by response magnitude (regression coefficient). Small pRF responds only to segment 17. Large pRF responds to both 17 and 18. Note effect on response modulation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We further explored our conjecture about the effect of pRF size by computationally predicting the fMRI response to our random stimulus using a modified version of the pRF modeling approach described by Dumoulin and Wandell (2008). Briefly, voxel pRFs were modeled by a 2D Gaussian sensitivity profile but with a range of sizes. The pRFs were then spatially convolved with 2D “binary” images of our random stimulus. The resulting time series was then convolved with an estimate of the HRF (Harvey and Dumoulin, 2011). Finally, “noise” was added to the predicted response waveform in order to achieve a signal-to-noise ratio of 5.0. These simulated fMRI responses were then subjected to the identical analysis procedure used previously for the empirical data. Fig. 10B shows a diagram of the random stimulus segments colored to represent the magnitude of response (regression coefficients) evoked in a voxel having the correspondingly sized pRF (represented by the white circles). As one might expect, the more segments the pRF encompassed, the lower was the response magnitude (compare colors in Fig. 10B). The full F-statistics also dropped from 28.3966 to 16.4522 as the pRF model increased in size from 2° to 6° vis. ang. indicative of a poorer fit. This simulation supports our interpretation that voxels in higher visual areas with large pRFs are likely to be poorly modulated by the random multifocal stimuli used in this study. We note that the results from the duty cycle experiment are also consistent with this interpretation in that sequences with a 25% duty cycle could allow better modulation of voxels having pRFs that encompass more than one stimulus segment.

Time efficiency

Temporal phase mapping proved to be more efficient than the random multifocal mapping paradigm used here on the basis of signal quality (z-score) versus time. This advantage was reduced, but not eliminated, by the need to run temporal phase mapping twice to obtain both angular and eccentricity data, both of which are obtained in a single run with the multifocal paradigm. Temporal phase mapping also has the disadvantage that to achieve accurate retinotopic assignments, the inherent BOLD hemodynamic delay that can vary widely across individuals (Handwerker et al., 2004) must be “corrected”. Currently accurate calibration requires either additional data acquisition with stimulus patterns moving in opposite directions, or a time consuming manual calibration process. Although calibration of the hemodynamic delays is not necessary for delineation of visual areas, other potential uses of retinotopic mapping do require calibration (DeYoe et al., 2011). Moreover, calibrated data should always be provided if possible, since use of such data in ongoing assessments of a patient clinically might inadvertently assume so. In sum, temporal phase mapping may be more efficient than random mapping in terms of the total time required in most, but not all applications.

Non-uniform representation of visual field locations

For both clinical and academic uses of fMRI, an often overlooked issue is whether fMRI maps are “functionally complete”. For visual cortex, this means that the maps should accurately reflect the underlying neural representation of the visual field within each cortical visual area. Fig. 5 reveals that both random-mapped and phase-mapped data for V1 show 2 or even 3 fold differences in the number of voxels responding preferentially to visual field locations near the horizontal meridians compared to the vertical meridians. In an ongoing study from our lab, we and others (Winawer et al., 2010) have argued that this irregularity likely reflects the effects of midline vasculature structures that alter BOLD signals by introducing excessive noise or atypical temporal phase delays. Such distortions may be inevitable as long as the signals arise from vascular-based mechanisms. The uneven representation of the visual field can obscure the detection of a true scotoma especially for random multifocal mapping which divides the visual field into a coarse array of stimulus segments. For instance, if a segment stimulates a region encompassing a scotoma plus a portion of responsive cortex, then the net number of voxels responding to this segment might appear similar to the number of voxels responding to some other segment that normally activates fewer voxels. On the other hand, phase mapping provides a continuous encoding of the visual field locations. The higher spatial resolution could facilitate the detection of a scotoma.

Effect of noise

The effect of noise on the retinotopic assignment of voxels for the random multifocal mapping and phase mapping are different based on our simulation. The effects can be summarized as either incorrect voxel assignment or as a reduction in the overall number of significantly activated voxels. At low SNR (0.5) the random multifocal mapping produced less significant voxels but made correct retinotopic assignments to all significant voxels. The reduction in the number of significant voxels could cause false detection of scotoma or incomplete mapping of the visually responsive cortex. At the same SNR, phase mapping activated more significant voxels but made some mis-assignments to the neighboring retinotopic locations. Currently a small number of incorrect, but neighboring assignments do not cause any practical problem, although conceptually they might shrink or enlarge a representation of a scotoma.

Identification of spatially complex pRFs

A region in the visual field that stimulates a voxel in the visual cortex is termed as its population receptive field (pRF) (Dumoulin and Wandell, 2008). Pathologies such as albinism and achiasma can cause abnormal routing of optic fibers, which can induce changes in visual cortex maps. For instance, voxels have been shown to have two mirror-imaged pRFs in the two hemifields (Hoffmann et al., 2003, 2012; Sinha and Meng, 2012). Such spatially complex receptive fields can also be associated with voxels that straddle the two banks of a cortical sulcus. Indeed, in a previous study using voxel sizes of 3.75 × 3.75 × 4.0 mm, we were able to identify voxels responsive to multiple separate visual field locations (Ward et al., 2012). In contrast, we were not able to find such voxels in this study using voxel sizes of 2.5 mm3. This finding might be explained by the results from another study on layer-specific BOLD activations in human V1 that indicated BOLD signal is mostly derived from cortical layer IV, which is about 1 mm under the cortical surface (Koopmans et al., 2010).

If such complex receptive fields are not detected in separate scans, temporal phase mapping will typically assign the pRF to a single erroneous location. This is due to the fact that the sum of multiple sinusoidal waveforms with the same period but different phases is, typically, a sinusoidal waveform with a phase different from either of the original phases. Random multifocal stimulation can identify the different components of a complex pRF without the need for a priori assumptions concerning the pRF structure. Other non-phase-dependent stimulation paradigms can also resolve such spatially complex pRFs. For example, multiple drifting-bar stimuli sweeping through the visual field in different orientations and directions, successfully identified two mirror-imaged responding locations for single voxels when used in conjunction with pRF modeling (Hoffmann et al., 2012). However, a priori knowledge of the spatial structure (e.g. two mirror imaged locations) of the pRF is required for this paradigm, while multifocal mapping is capable of resolving pRFs with any type of spatial structure without making any assumption of the number or distribution of the responding locations.

Conclusion

The use of fMRI for clinical applications such as pre-surgical mapping impose constraints such as time efficiency and mapping completeness which may be of higher priority than for academic research. MRI scan time is a major factor affecting scan success as well as the number of patients that can be scanned per day. Completeness is important in order to ensure that a functionally critical structure near a site of pathology is mapped as fully as possible so that the neurosurgeon can avoid damaging eloquent cortex and causing post-operative deficits. Compared to random multifocal mapping, temporal phase mapping was fundamentally more efficient and able to map extrastriate visual areas more completely. Its shortcomings such as retinotopic mis-assignments caused by noise and hemodynamic variations are modest and likely not to be critical for most clinical applications unless retinotopy with very high precision is needed, as might be the case when trying to assess rehabilitation-induced recovery of cortical function near the margin of a scotoma. Random multifocal mapping provides simultaneous eccentricity and angle mapping, and does not need hemodynamic calibration. But spatial resolution is limited by the number of stimulus segments and its inability to fully map extrastriate visual cortex may limit its utility when a pathology involves lateral or dorsal occipital cortex. With increased noise, random multifocal mapping will virtually always provide a correct retinotopic assignment but the number of significant voxels may be reduced. In short, though phase mapping is probably the better choice for clinical applications such as pre-surgical brain mapping, random multifocal mapping also offers some unique advantages that might be optimal for some specific applications.