Accumulation System: Distributed Neural Substrates of Perceptual Decision Making Revealed by fMRI Deconvolution

Stimuli

All the stimuli were digital photographic grayscale images of 640 (horizontal) × 480 pixels, whose original versions were selected from clip-art collections (Sozaijiten; Image Navi Corporation) and personal photographs. The images were presented using the software Presentation (Neurobehavioral Systems) onto a back-projection screen mounted outside the MRI scanner bore behind the participants' head, with subtending 8.0 × 6.0° of visual angle. In the decision section, we used 90 images consisting of 30 stimulus images of each of the three categories: human, scene, and tool. A “human” image depicted humans with both the head and body parts. A “scene” image included architectural scenery objects, such as houses and bridges. A “tool” image included handheld tools such as sports goods and stationery. We produced 15-level degraded versions of each image by means of spatial low-pass Fourier filtering, whose high cutoff frequency was manipulated from 1.0 to 4.5 cycles per degree (on the screen) with 15 levels of equal increments. Using the gradual unmasking technique (James et al., 2000; Carlson et al., 2006; Eger et al., 2007), the 15 versions were sequentially presented every 1.4 s from the most to the least degraded level (Fig. 1A). Participants were required to categorize the stimulus images as one of the three categories of human, scene, and tool. In the practice before the data collection, we also used 15-level degraded versions of another 20 images of each category.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Examples of stimulus images. A, Images used in the decision section, which were ones of three categories: human, scene, and tool. Although the original versions of the images were easy to categorize, each of them was low-pass filtered to produce 15-level degraded versions, which were sequentially presented by every 1.4 s from the most degraded level (level 1) to the least (level 15). The times below the levels indicate the presentation times from the stimulus onset. B, Images used in the localization section, which were ones of five categories: face, body, human, scene, and tool. The backgrounds of the target objects were filled with uniform gray, except for the scene images. The mean luminance of each image was adjusted to be the same among all the stimulus images to avoid any effects caused by changes of stimulus luminance.

In the localization section, we prepared 200 images consisting of 40 images of each of the five categories: face, body, human, scene, and tool. Each image included an exemplar object in the assigned category (Fig. 1B). A “face” image included a human face except other body parts, while a “body” image included parts or a whole of a human body except a head. A “human” image included a whole human body with head. A “scene” image includes natural scenery objects such as mountains or the same types of architectural scenery objects as in the decision section. A “tool” image included the same types of handheld tools as in the decision section. The backgrounds of the target objects were filled with uniform gray, except for the scene images, whose backgrounds were kept as the important content of the expanse of scenery. Unlike the decision section, all the images were not degraded but were clear enough to recognize the depicted objects easily and unambiguously.

Acquisition and preprocessing of MRI data

MRI images were acquired using a 3T MR scanner (Trio, Siemens). For functional imaging, an ascending T2*-weighted gradient-echo echoplanar imaging (EPI) procedure was used to produce trans-axial slices covering the entire cerebrum and cerebellum, excluding the eyeballs, by oblique scanning. Imaging parameters in the decision section were repetition time (TR) = 2000 ms, echo time (TE) = 30 ms, flip angle (FA) = 79°, field of view (FOV) = 192 mm, 64 × 64 matrix, 34 slices, and voxel size = 3.0 × 3.0 × 3.45 mm. In the localization section, TR = 2500 ms, TE = 30 ms, FA = 80°, FOV = 192 mm, 64 × 64 matrix, 38 slices, and voxel size = 3.0 × 3.0 × 3.0 mm. To ensure stabilization of the functional imaging data, the first five scans (10 s) of the decision section and the first four scans (10 s) of the localization section were discarded (not used in the subsequent analyses).

The acquired MRI data of raw DICOM format were converted to NIFTI format and preprocessed using the software Statistical Parametric Mapping 12 (SPM12; Friston et al., 1995; Wellcome Trust Center for Neuroimaging, London, United Kingdom; https://www.fil.ion.ucl.ac.uk/spm), implemented in MATLAB 2016b (MathWorks). For the data of each session, the displacement of scans caused by the participants' head motion was realigned relative to the mean of the images using the affine transformation, modifying the headers of the images to reflect their relative displacement. In the determination of regions of interest (ROIs), to improve the accuracy of blocked-design SPM analysis, we applied the slice timing correction, which concerns differences in signal acquisition times because of the positions of slices and modifies each time series in the slice to have values that would have been obtained if the slice had been acquired at the same time as the reference slice. The reference slices were the middles of the scans, which were 17 and 19 in the decision and localization sections, respectively. In the temporal-profile analysis of data of the decision section, however, we did not apply the slice timing correction; instead, we conducted the slice-timing modeling by which we obtained not-modified values of fMRI data with their accurate acquisition times (Kiebel et al., 2007) as described below in the subsection, Profile classification of deconvolved signals. Spatial normalization was applied to the data with nonlinear three-dimensional transformations into the ICBM template of East Asian brains, which we used because all the participants were native Asian (Japanese). The normalized images were spatially smoothed in a three-dimensional Gaussian kernel with a half-maximum width of 6 mm.

Experimental sections

Each participant had the decision section first and the localization section on a different day to prevent fatigue. The decision section consisted of 90 trials, which were given in six sessions of 15 trials each. One stimulus image was used in each trial, and 90 trials covered 30 images for each of the three categories (human, scene, and tool). As shown in Figure 2, each trial lasted for 36 s, which began with a 2-s rest epoch in which a white fixation point in the center of a uniform gray background was presented, followed by a 1-s cue epoch in which a black word of the cue for the oncoming stimulus was presented. Intervened by a 2-s rest epoch, a stimulus epoch lasted from 5 s to 26 s, in which a set of 15-level degraded versions of a stimulus image were sequentially presented every 1.4 s according to the gradual unmasking technique. The stimulus epoch was followed by a 12-s rest epoch, which was composed of the remaining 10 s of the present trial and 2 s at the beginning of the next trial. This length of the rest epoch was expected to be long enough for the BOLD signal elevated in the preceding stimulus epoch to approach the baseline levels (Tremel and Wheeler, 2015). The beginning of a trial was precisely synchronized with the fMRI scan with a TR of 2 s, so that a trial comprised 18 scans and the stimulus onset (at 5 s) was at the midpoint of the third scan of the trial. The least common multiple of the TR (2 s) and the unmasking step interval (1.4 s) was large enough (14 s) so that the image-change moments of the unmasking steps were little synchronized with the scan timing. A cue word was given ahead of each stimulus epoch to investigate the influence of prior knowledge on categorical decisions. We expected that a cue congruent with its following stimulus would facilitate decision while an incongruent cue would obstruct the decision to the contrary, so that cue congruency was expected to influence the accumulation of the decision process. Either “Human,” “Scene,” “Tool,” or “****” (indicating no cue) was displayed in characters at the center of the screen with a size around 5° of visual angle. (Quotation marks, “ ”, are not displayed.) To evaluate the effects of cues, we gave stimulus-congruent cues, incongruent cues, and no cues, each of which were in one third of all the trials; for example, 30 epochs of human stimuli were preceded by “Human” in 10 trials, “Scene” and “Tool” in five trials each, and “****” in 10 trials. The presentation order of the cue-stimulus combinations was randomized. Participants were required to discriminate the stimulus category among human, scene, and tool and respond by clicking the assigned button as soon as they could discriminate it. The buttons for the index, middle, and ring fingers were assigned to human, scene, and tool, respectively. We encouraged participants not to hesitate to respond by allowing the second clicks to revise their answers. Participants were told that only clicks during a stimulus epoch were valid. As a result, 68% of all trials across participants responded correctly in the first response, and in 96% of all the trials correct answers were obtained in the final (including the second and later) responses, showing that most of the stimuli with the least degradation could be correctly categorized. In the following fMRI time series analysis, we used only the trials with the first correct responses not to intermingle the activity of the second-decision process. To monitor the subjective verification of the participants' recognition, the instruction asked participants to double-click the same button as their previous final response at any time within the stimulus epoch when they became confident of the answer. The results of the present study showed that participants verified their responses in 86% of all the correct-first-response trials, indicating that participants consistently recognized the images through the stimulus epochs of most of the correct trials. Relative temporal positions of the verification-RTs were shown to be almost uniformly distributed in the intervals between the correct first response and the end of stimulus epoch of trials with a mean of 0.54 and SD of 0.25 (across trials and participants) by comparing with a theoretical uniform distribution that has the mean of 0.50 and SD of 0.29. Thus, influences of the verification responses on activity profiles seem to be minimal. One session comprised 278 scans (556 s), the initial five scans of which were discarded. The remaining 273 scans, which consisted of 270 scans of 15 trials by 18 scans and an additional three scans after the last trial, were used in the following analysis. A decision section typically took ∼100 min for the data collection, including breaks inside the MRI scanner and ∼50 min for instructions, practice using another 20 images of each category, and postexperiment debriefing outside the scanner.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Design of a trial of the decision section. Each trial lasted for 36 s, which was synchronized with the 18 scans of fMRI with a TR of 2 s. R and C indicate the epochs of rest and cue, respectively. In the stimulus epoch indicated by S, whose onset (at 5 s) was at the midpoint of the third scan, a set of 15-level degraded versions of a stimulus image were sequentially presented with a changing interval of 1.4 s causing little synchrony with the scan timing.

The localization section consisted of four sessions, each of which included 10 trials. Each trial was dedicated to one of the five categories of face, body, human, scene, and tool, and eight trials for each category were conducted in a randomized order. Each trial began with a 2-s rest epoch (only a fixation point), followed by an 18-s stimulus epoch and a subsequent 10-s rest epoch, constituting a 30-s duration. During a stimulus epoch, the first image was presented for 2 s followed by a 1-s fixation point, and 10 images of the same category were presented in successive turns every 1.5 s. The first image appeared once or twice again in the subsequent 10 images of randomly selected turns. All the other nine or eight images were different. To enhance attention to the stimuli, participants were required to remember the first image and to press an assigned button when the same image appeared again within the stimulus epoch. Each of the 200 images was used twice in this section, covering the 400 images used in the 40 trials in total. Each trial began in synchrony with the fMRI scan, with a TR of 2.5 s. One session comprised 126 scans (315 s), the initial four scans of which were discarded. The remaining 122 scans, which consisted of 120 scans of 10 trials by 12 scans (30 s) and an additional two scans after the last trial, were used in the following analysis. A localization section for a participant typically took ∼70 min including instructions, practice, short breaks, and debriefing.

Determination of ROIs

As for regions associated with perceptual decision making, we determined ROIs that showed the decision task-related positive activation. To determine coordinates of each ROI, we adopted “individual peaks within group activation” method (O'Reilly et al., 2012), in which we selected the voxel that showed the strongest task effect in the individual-level analysis of each participant with applying an inclusive mask that was the activation ROI detected by the group-level analysis. This approach provides the most sensitive results of the individual-level effect, allowing for interindividual differences within the functionally identical ROI. Using the preprocessed (slice-timing corrected) fMRI data, the following three steps were conducted using SPM12. In the first step, using the decision section data, we determined positively activated regions during the decision process regardless of the stimulus category, to use it as an inclusive mask of positive activation. The reason to restrict ROIs to positive activation was that the accumulator regions, which were the primary targets of this study, should show positively-growing activity during the decision process. We defined a model of task-related positive activation as being positive relative to the baseline throughout all the stimulus epochs of the decision section; thus, the model of positive activation was common to all the sessions, whose stimulus epochs had the same temporal structures. This model was convolved with the canonical hemodynamic response function (HRF; Friston et al., 1998) for use as the regressor of positive activation in a blocked design SPM analysis. In the individual-level (first-level) analysis, the fMRI data of each participant were analyzed by a general linear model (GLM) using the regressors, applying high-pass filtering with a cutoff period of 128 s, and serial autocorrelation estimation with a first-order autoregressive model to eliminate artifactual low-frequency trends and aliased biorhythms. A single-condition contrast of the positive activation was applied to the parameter weights (β values) of the regressor, and a t-contrast image of an individual participant was constructed. The t-contrast images of all participants were incorporated into a random-effects model to make inferences of a group-level (second-level) analysis (Holmes and Friston, 1998). Activation was detected with a voxel-wise threshold of p < 0.001 (uncorrected for multiple comparisons) and a cluster-wise threshold of the familywise error rate (FWER) p < 0.05, corrected for multiple comparisons based on random field theory (Friston et al., 1994b), which corresponded to a cluster-size threshold of 1257. Then, a spatial mask of positive activation was made of all the voxels detected with this cluster-wise threshold of the group-level analysis (Fig. 3A), and used as an inclusive mask common to all the individual participants in the following steps.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Procedure to determine coordinates of ROIs. A, The mask of positive activation, which consists of all the voxels detected in the group-level analysis of positive activation throughout stimulus epochs of the decision section (cluster-wise threshold of the FWER p < 0.05). B, An example of coordinates of a CSR: left body-selective region. The light green region indicates a body-selective region, which was obtained in a group-level analysis (FWER p < 0.05) with inclusive application of the mask of positive activation. Following the “individual peaks within group activation” method (O'Reilly et al., 2012), within this region, the peak voxels of body-selective activation of all the participants were extracted (small dark red spheres) from the individual-level t-contrast images (voxel-wise threshold p < 0.001, uncorrected). Their mean coordinates (large red sphere) represent the location of this body-selective region. The SD of coordinates of participants was 7.4 mm in this example of our data. C, An example of coordinates of an NCSR: the left anterior insula. The light blue region indicates the left anterior insula of the TPM atlas masked inclusively by both the masks of positive activation and noncategory selectivity. Within this region, the peak voxels of decision task-related activation of participants were extracted (small dark red spheres) from the individual-level t-contrast images (voxel-wise threshold p < 0.001, uncorrected) obtained from the decision-section data. Their mean coordinates (large red sphere) represent the location of the left anterior insula. The SD of coordinates of participants was 4.8 mm in this example of our data.

In the second and third steps, choice-selective and nonchoice-selective regions were determined. The choice-selective regions of this study were visual CSRs because the task of the decision section was to discriminate the object categories. The nonchoice-selective regions, to which we refer as non-CSRs (NCSRs), were defined as all regions that did not show any selectivity for the concerned categories. The second step was to determine the coordinates of the CSRs of individual participants using data from the localization section. In the SPM analysis, we modeled five regressors at the individual level representing the stimulus epochs of the five categories: face (F), body (B), human (H), scene (S), and tool (T). The GLM of these regressors was applied to the preprocessed data of the localization section, including the high-pass filtering and serial autocorrelation estimation in the same way as in the first step. The t-contrast images of category selectivity were obtained using four contrasts of the weighted regressors, which were designed according to the previous studies of CSRs as follows: F > (S + T)/2 for face selectivity (Kanwisher et al., 1997), B > (S + T)/2 for body selectivity (Downing et al., 2001), S > (H + T)/2 for scene selectivity (Epstein and Kanwisher, 1998), and T > (H + S)/2 for tool selectivity (Grill-Spector et al., 1999; Chao et al., 1999). Here, we used the regressor of H to define scene and tool selectivity because it is not clear whether the F-regressors and B-regressors, objects of whose categories are evidently related, can be used equivalently with the T-regressor or S-regressor in making the contrasts. With the H-regressor, all the kinds of selectivity could be defined by contrasting with two other categories in a consistent manner as above. For each category selectivity, the individual-level t-contrast images of all participants were incorporated into a group-level analysis with inclusive application of the positive-activation mask. The selective regions were detected with a voxel-wise threshold of p < 0.001 (uncorrected) and a cluster-wise threshold of the FWER (p < 0.05, corrected for multiple comparisons), which corresponded to cluster-size thresholds of 90, 80, 161, and 921 for selectivity of face, body, scene, and tool, respectively. Since separate (discontiguous) clusters of voxels detected at the group level with even the same contrast could be distinct functional regions of that category selectivity, each contiguous cluster was regarded as a distinct CSR, resulting in 3, 5, 5, and 1 regions selective for face, body, scene, and tool, respectively. Each of these CSRs was applied as an inclusive mask to the individual-level t-contrast image of each participant regarding the corresponding category selectivity, and we extracted all the local peak voxels whose t-values were above the voxel-wise threshold p < 0.001 (uncorrected). In six of the 14 CSRs above, two or more of the 14 participants showed no significant voxels; thus, we excluded these six CSRs from the following analysis. For the remaining eight CSRs, which were two face-selective, two body-selective, three scene-selective, and one tool-selective region, the Montreal Neurologic Institute (MNI) coordinates of each region were determined for each participant by selecting the peak voxel that showed the maximum t-value for its selectivity condition. Since in the later stages we averaged activity signals from the voxels of the participants in the same region, we expected that reducing the dispersion of the locations among participants of an identical region would contribute to lowering functional variability and clarifying common activity profiles among participants. Thus, to diminish the spatial dispersion of the coordinates among participants, we conducted an additional operation in which the location of a participant with the largest deviation from the mean location of all the participants was replaced with the peak of the next greatest t-value of that participant (only when the participant had two or more peaks) if this replacement reduced the variability of the locations. The replacement was allowed to repeat for a participant who showed the largest deviation among participants until no reduction in the variability was obtained. As a result, 11 peaks (10%) of 112 regions (8 CSRs by 14 participants) were replaced by the second greatest peaks, and this operation reduced the SD of the coordinates among participants by 0.7 mm, resulting in 7.8 mm on average of the eight CSRs. The location of each region was represented by the mean coordinates of all the participants of that region (Fig. 3B).

In the third step, by inversely using the category-selective maps, we defined NCSRs that showed positive activation without any significant category selectivity. The second step produced 56 images of category-selective activation in the individual-level analyses for four categories by 14 participants with a voxel-wise threshold of p < 0.001 (uncorrected). Using the Image Calculator of SPM12, these 56 images were summed up and inverted (by taking 1 for 0 and 0 for non-zero values of the summed image) to obtain a binary image, all of whose voxels did not show any category-selective activation of any participant with such a liberal threshold. This image was used as a spatial mask for noncategory selectivity. To parcel the noncategory-selective range, which extends over the whole brain, into the anatomic structures, we employed the anatomic atlas of the tissue probability map (TPM) available with SPM12. The TPM provides the MNI-space volume data of 136 tissue structure-based regions of the brain, including both neural and non-neural types. Excluding 15 structures that lack neuronal cell bodies, such as white matter, optic chiasm, vessels, and ventricles, we used the volume data of 121 atlas regions of cortical areas and nuclei. The intersection of each of these TPM-atlas regions and both masks of noncategory selectivity and positive activation provided a TPM-based NCSR that showed positive activation during the decision task. We refer to this volume of intersection simply as an NCSR. Each of these 121 NCSRs was applied as an inclusive mask to the individual-level t-contrast images of the task-related positive activation that were obtained from the decision-section data in the first step, and we extracted all the local peak voxels whose t-values were above the voxel-wise threshold p < 0.001 (uncorrected). In 66 of the 121 NCSRs, two or more of the 14 participants showed no significant voxels; thus, we excluded these 66 regions from the following analysis. We determined the MNI coordinates of each of the remaining 55 NCSRs for each participant in the same way as the CSRs; that is, the coordinates of the voxel with the maximum t-value were adopted, followed by an additional operation in which the voxel of a participant deviated most from the mean location of all participants was replaced with a voxel of the second greatest t-value of that participant if this replacement reduced the location variability. As a result, 137 peaks (18%) of 770 regions (55 NCSRs by 14 participants) were replaced by the second greatest peaks, and 10 peaks (1%) were replaced by the third or fourth greatest peaks. This operation improved cohesiveness of the peak locations among participants, reducing the coordinates SD by 2.5 mm resulting in 5.0 mm on average of the 55 NCSRs. The location of each region was represented by the mean coordinates of all the participants of that region (Fig. 3C).

Deconvolution of BOLD signals

Using the preprocessed, normalized, but not slice-timing corrected fMRI data of the decision section, BOLD signals were extracted from voxels at the coordinates of the eight CSRs and 55 NCSRs determined above for each participant. The BOLD time series of each session (273 data points) was linearly detrended, keeping its mean value, and we estimated its baseline level with the mean of 45 signal values of the first, second, and 18th scans of the 15 trials, which were temporally the most remote rest scans from the previous stimulus onsets. We used the 271st and 272nd scans just after the last (15th) trial instead of the first two scans of the first trial, expecting the baseline activity to be more stabilized. Using this baseline, the time series was scaled to percent signal change (PSC) from the baseline. Here, we introduce deconvolution of the PSC of BOLD signals aiming to improve availability of temporal information of the signals on the following grounds. Because of the linearity of the hemodynamic system, a BOLD signal can be well approximated by convolving the neural activity with a HRF that is an impulse response (Friston et al., 1994a; Boynton et al., 1996). Since the convolution makes the BOLD signal delayed and blurred from the original neural activity, it is difficult to determine the timing of the original activity using the BOLD signal in a straight manner. On the other hand, the accumulators are characterized by their temporal profiles of activity regarding the onsets and peaks relative to the stimulus onsets and subject's responses. Therefore, recovery of the neural activity from BOLD signals can be expected to be effective in searching for the accumulators, and we employed Wiener deconvolution as one of the powerful methods for such recovery. The PSC time series was used to estimate the original neural activity using the Wiener deconvolution (Papoulis, 1977; Glover, 1999; David et al., 2008; Surhone et al., 2011; Wu et al., 2013) with the following details.

Let s(t) and h(t) be the time series of neural activity and HRF, respectively. The measured time series of the BOLD signal in PSC, m(t), originating from the neural activity, is expressed as following: m(t)=h(t)∗s(t) + n(t), where ∗ denotes the convolution operation of the neural activity and HRF, and n(t) denotes the noise in the measurement. Using the Wiener deconvolution filter w(t), the neural activity is estimated as ŝ(t) from the BOLD signal m(t) as follows: ŝ=w(t)∗m(t).

According to the Wiener deconvolution theory, W(f)=1H(f)[|H(f)|2|S(f)|2|H(f)|2|S(f)|2+|N(f)|2], and ŝ(t)=IFT[W(f)M(f)],(1) where W(f), M(f), H(f), S(f), and N(f) are the Fourier transforms of w(t), m(t), h(t), s(t), and n(t), respectively, and the squared forms such as |N(f)|2 are the energy spectral density (ESD) of the corresponding quantity. IFT indicates the inverse Fourier transformation. For simplicity, we also refer to W(f), the Fourier transform of the Wiener filter w(t), as the Wiener filter. Postulating that the neural activity (signal) and noise are uncorrelated, then |M(f)|2=|H(f)|2|S(f)|2 + |N(f)|2, and using this relation we obtain W(f)=1H(f)[1−|N(f)|2|M(f)|2].(2)

This indicates that the Wiener filter attenuates the Fourier components with larger ratios of the noise relative to the measured signal, consequently diminishing noise in the deconvolved signal. Thus, using the appropriate functions of H(f) and N(f), we can obtain an approximation of s(t) from the measured data m(t).

The measured data M(f) at high frequencies seem to be dominated by the noise component N(f) because the HRF works as a temporal low-pass filter that substantially cuts off high-frequency components of neural activity. As shown in Figure 4A, this property was verified by the ESD |M(f)|2 of the PSC time series of this study, averaged across all the sessions over regions and participants, which displays the flat energy content over high frequencies above ∼0.2 Hz suggestive of the constant noise components. This grand-averaged ESD (Fig. 4A) was used to estimate the mathematical form of |M(f)|2, which was applied to data of individual regions below. Following Glover (1999), we assumed a constant spectrum of noise (white noise) as |N(f)|2=N02, whose magnitude was estimated from the high-frequency components of |M(f)|2. We estimated the noise energy logN₀² by averaging the values at the 20 highest frequencies above 0.214 Hz, as indicated by the thick line in Figure 4A. As a typical function describing this type of data distribution, we fitted an exponential function with the asymptote logN02: log|M(f)|2=Cexp(−Df)+log N02, to the ESD data by means of the least squares method, excluding the DC component (f = 0 Hz) and the harmonics of trial frequency (1/36, 2/36, …, 6/36 Hz because of a trial period of 36 s). The results showed that the exponential function fitted well with the grand-averaged ESD of the measured data (goodness of fit, r² = 0.99). We postulated that this mathematical form could be applied to individual data, and the same estimation method was used to obtain the parameters of the exponential function best fitted to the ESD data of each region of each participant (averaging the six sessions). In this ESD estimation, we excluded two regions of a single participant (#11) whose ESDs had smaller values at low frequencies than the noises so that the exponential functions could not fit. In the other cases, the exponential functions were successfully fitted to the ESD data of individual regions of participants with r² = 0.64 on average, as shown by the examples in Figure 4B. The values of parameters C and D obtained for each region of each participant provided the function: log(N02|M(f)|2)=−Cexp(−Df), which was used to construct the Wiener filter of Equation 2.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Construction of the Wiener filter. A, The red symbols with lines indicate the one-sided ESD |M(f)|2 of the signals in PSC on a logarithmic scale, grand-averaged across all the sessions over regions and participants. The frequency ranged from 0 Hz (DC component) to 0.2491 Hz with the interval of 1/546 Hz, given by the Fourier transformation of the time series of 273 sample points with TR of 2 s. The values at the 20 highest frequencies (from 0.2143 to 0.2491 Hz), whose range is indicated by the black thick line, were averaged to estimate the noise energy logN02, whose magnitude is indicated by the black dotted line. The period of a trial (36 s) caused the outstanding increases at the harmonics of 1/36 Hz. Excluding these harmonics of the trial frequency (3 points centering at each peak of the first to sixth harmonics) and the DC component (at 0 Hz), we fitted log|M(f)|2−logN02 with an exponential function Cexp(−Df) by means of the least square's method. The result showed a successful fitting (r² = 0.986) with C = 1.32 and D = 14.3 s. B, Examples of the ESDs of individual regions of individual participants. The noise energy estimation and fitting of the exponential functions were conducted in the same way as A. First example: the left body-selective region of participant #02. C = 1.23, D = 9.8 s, r² = 0.63. Second: the left anterior insula of participant #10. C = 1.97, D = 16.1 s, r² = 0.86. C, Broken line, The average time course of the model functions fitted to the 62 samples of the Handwerker dataset of HRFs. Red line, The model function fitted to the average time course, which was used as the dataset-based empirical HRF in the Wiener deconvolution (see the main text for the parameters). Blue line, For reference, the canonical HRF (Friston et al., 1998; provided by SPM12) with the default shape parameters (α1= 6, α2 = 16) and normalized amplitude (A₁ = 5.70, A₂ = −0.95). The dotted line indicates the zero-baseline. All the functions were calculated with the time step of 0.1 s and normalized to start at 0 and peak at 1. D, The two-sided ESD of the empirical HRF of (C) |H(f)|2, which was used in the Wiener filter. (The function is shown on the lower half of the frequency axis.) E, The magnitude of the Wiener filter |W(f)| of Equation 3, where δ = 24. This panel shows the Wiener filter obtained by using the function N02/|M(f)|2 of the grand-averaged ESD (panel A), whereas in the deconvolution of data the filter was obtained for each region of each participant.

In contrast to M(f) and N(f), which can be estimated using the measured data, we do not know the particular H(f) of the concerned regions in most cases, so that we need to presume a certain form of HRF. Whereas the magnitude of convolution h(t)∗s(t) is constrained by the measured BOLD data, the absolute magnitude of the HRF is arbitrary because the magnitude of the neural activity is not defined in general. Thus, we cannot compare the magnitudes of the deconvolved signals among different regions with a presumed HRF whose magnitude is arbitrarily given. It should be noted here that effective information provided by the deconvolution method in this study was not magnitude (except comparing conditions within the same region) but temporal profile (shape) of activity. It is also critical to consider the fact that the shapes of observed HRFs, which are obtained as BOLD responses to an impulse neural activity, vary across regions and individuals (Aguirre et al., 1998; Handwerker et al., 2004). Our strategy to overcome this problem of the HRF variability was to estimate the sizes of the variability effects on deconvolution results and to show that the effects can be smaller than the fMRI sampling timescale. As the HRF h(t) used in the Wiener deconvolution, we adopted a typical form of HRF, which is characterized by a positive peak and the postpeak undershoot that returns to the baseline within a proper timescale according to the properties described by many studies. The typical HRF was determined by averaging substantial amount of empirical data of human HRFs with excluding apparently atypical cases such as ones whose undershoots do not return to the baseline. Then the Wiener deconvolution using this typical HRF was applied to empirical HRF data without excluding any atypical or outlier cases to estimate the effects of the HRF variability on the deconvolution results. By courtesy of Daniel Handwerker (National Institute of Mental Health) and his collaborators, we were permitted to use the dataset of BOLD impulse responses presented by Handwerker et al. (2004), which includes 80 samples of the BOLD responses empirically obtained from four different regions [the primary motor cortex (M1), primary visual cortex (V1), FEFs, and supplementary eye fields (SEFs)] of 20 participants. Each sample of the BOLD RT series was scaled to 0 at the initial point and 1 at the maximum. It is shown that a reasonably good fit to the BOLD impulse response is provided by a model function comprising a sum of two γ distribution functions (GDFs; Boynton et al., 1996; Friston et al., 1998), and thus we fitted each sample with this type of model. A GDF was given by G(t,α,θ)=tα−1e−tθθαΓ(α), where α and θ are shape-determining parameters and Γ(α) is a γ function. We fixed θ to 1 for simplicity and omitted it from the fitting scans. The model function was given by: h(t)=A1G(t,α1)−A2G(t,α2), where the first and the second GDFs model the positive peak and the undershoot of the HRF, respectively. In our fitting, the four parameters determining amplitudes (A1, A2) and shapes (α1, α2) were scanned with a precision of 0.01, to minimize the sum of squared residuals. Because the peak of a GDF (with θ = 1) is obtained at α−1, we set the initial values of the scan parameters for the first and second GDF using the maximum and minimum peaks of the sample data, respectively. Ten of the 80 samples, whose undershoots did not return to the baselines, had no optimal fits with α2 under the maximum scanning limit of 23; therefore, these 10 samples were excluded from the following averaging. Another 8 samples had outliers (apart from the quartiles by 1.5 times the interquartile range) of at least one of the four parameters, so that they were also excluded. The fittings were as good as r² = 0.940 on average for the remaining 62 samples, with an improvement from 0.934 for including the 18 excluded samples above. The model functions of the 62 samples were calculated from 0 to 32 s with a precision of 0.1 s, normalized to start at 0 and peak at 1, and averaged at each time point. The average model function was again normalized and fitted by a model function, which provided us the dataset-based empirical HRF with parameters α1 = 5.10, A1 = 5.21, α2 = 11.55, and A2 = −1.89 (Fig. 4C). We used this empirical HRF to construct the Wiener filter. Relating to H(f), whose ESD |H(f)|2 is shown in Figure 4D, Equation 2 has the problem of diverging (very large) magnitudes of the filter W(f) at high frequencies, where the magnitude of the denominator H(f) is close to zero. With diverging W(f), even small fluctuations in the measured data M(f) can cause considerable disturbances in the deconvolved signals in Equation 1. To prevent this problem, we introduced a positive term δ to increase the minimum of the denominator from zero, as follows: W(f)=H*(f)|H(f)|2+δ[1−N02|M(f)|2],(3) where * denotes a complex conjugate. It should also be noted that δ emphasizes the particular frequency components of W(f), where |H(f)|2 = δ because the function x/(x²+δ) has a maximum peak at x² = δ. To prevent this distortion, we used δ = 24, which is larger than the maximum |H(f)|2, so that the influence of δ becomes monotonic over frequency. In this study, we used Equation 3 as the Wiener filter, whose magnitude |W(f)| depends on frequency as shown in Figure 4E.

We estimated the influences of the HRF variability over regions and individuals on the activity signals obtained by the deconvolution filter of Equation 3. As shown in Figure 5, we used a model neural activity of a Gaussian function (mean of 0 s and SD of 2 s) and convolved it with the 80 different HRFs of the dataset of Handwerker et al. (2004) to generate BOLD signals. Subsequently, the Wiener deconvolution (Eq. 3) was applied to these BOLD signals to recover the original (model) activity. A two-way ANOVA showed that the deconvolved signals had peak times dependent on both region (F_(3,57) = 6.60, p < 0.001) and individual (F_(19,57) = 8.03, p < 0.001). The peak times of regions averaged over individuals were 0.52 ± 0.86 s (SD) for M1, and −0.16 ± 0.98 s for V1, 0.30 s ± 0.92 s for FEF, and −0.08 ± 0.91 s for SEF. Using these values, we estimated 95% confidence intervals of the mean peak times of the regions averaged over 14 individuals, which were expressed as mean±1.96×SD/n, where n = 14: [0.07, 0.97] (s) for M1, [−0.67, 0.35] for V1, [−0.18, 0.78] for FEF, and [−0.56, 0.40] for SEF. Thus, we could expect to obtain the mean peak times in the bin of the true (model) activity peak time (0 s), which was [−1, 1] with the same width of the sampling time of fMRI (2 s). While the precision of the deconvolved peak time, which is given by the standard error of the mean, can be improved by increasing the number of participants (n), the accuracy of the deconvolved peak time of a region, which is the closeness of the mean to the true value, is limited by the properties of HRF intrinsic to the region owing to multiple factors such as vasculature and baseline cerebral blood flow (Handwerker et al., 2004). The estimation above, however, showed that the intrinsic peak-time differences among regions are not significantly large relative to a typical fMRI sampling time (2 s), so that detected differences can be significant in most cases. We also estimated the variability of the rise time of the deconvolved signal, which was defined as the time to reach at 0.3 (30% of the amplitude). The mean difference of rise time from the true (model) value averaged over individuals and regions was −0.76 ± 0.71 s (±SD), indicating that the deconvolved signal rose earlier than the model activity. It seems that this earlier rise was caused by the low-frequency transparency of the Wiener filter (Fig. 4E), which is an unavoidable property for the filter to eliminate the high-frequency noise. The above estimations aimed to theoretically show that the differences in the deconvolved signals could reflect significant differences in the original activity despite the HRF variability over individuals and regions. The experimental data, however, include variability caused by many other factors, so that taking a peak of averaged time course could be more effective than averaging peaks taken from time courses, the former of which was the way this study adopted in the following.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Procedure to estimate precision of the present deconvolution method. A model neural activity was given by a Gaussian function with mean 0 s and SD 2 s, whose peak time and rise time (defined as a time for activity to reach 30% of the amplitude) were 0 and −3.1 s, respectively. The model activity was convolved with 80 HRFs obtained from the dataset of Handwerker et al. (2004) to generate BOLD signals, which reflected the HRF variability over participants and regions. Here, we used all the 80 HRFs of the dataset, including HRFs of the outliers, for the purpose of estimating the effects of HRF variability (the figure displays 20 representative HRFs from the dataset). We did not add noise to the model activity to emphasize influences of HRF variability. These BOLD signals were deconvolved by the Wiener filter Equation 3, which was constructed with the empirical HRF. All the time series of this figure were calculated with the time step of 0.1 s and normalized between 0 at the baseline and 1 at the maximum.

Profile classification of deconvolved signals

As illustrated in Figure 6A, a time series of the PSC of a session of the decision section was extracted from each ROI voxel of each participant, and was deconvolved by the Wiener filter. From the deconvolved time series, we extracted segments of trials in which the participant's first responses were correct (68% of all the trials as described above) to avoid intermingling activity of the second-decision process. To evaluate the temporal profiles of activity, the trial segments of each region were averaged separately for participant's RTs of the trials because the activity dependence on decision moment is a critical property of the accumulator regions, and the RT approximates the decision moment. The trial segments were grouped by RT of the trial, with intervals of [5, 7), [7, 9), [9, 11), [11, 13), and [13, 15) (s). (The notation […) indicates an interval with a closed low end and open high end.) We did not use trials with RTs shorter than 5 s or longer than 15 s to avoid high error rates and low observed numbers (see Results). The RT-grouped trial segments were averaged across sessions and participants in both stimulus-locked and response-locked alignments, which were defined by the signal acquisition times relative to the stimulus onset and the participant's response, respectively (Fig. 6B). Signal acquisition times were obtained using the slice-timing modeling. In this method, the original slice position of the ROI voxel before the spatial normalization of the image was recovered by reading the normalization parameter file (y_*.nii) produced by SPM12, and the accurate (modeled) slice acquisition time of the voxel was calculated based on the mean relative position of the original slice in scans. For averaging the trial segments, we used time coordinates with 18 bins with a 2-s width (corresponding to the length of a trial) whose centers were multiples of the TR (2 s). The 18 data points of each trial segment were assigned to the bins according to their acquisition times relative to either the stimulus onset or the response. In the stimulus-locked alignment, the bins for all the ROI voxels were made completely matched the scans because the stimulus onset was at the midpoint of the third scan of the trial (as described above); thus, each of the 18 time bins whose centers were −4, −2, 0, …, and 30 s covered acquisition times of all the ROI voxels within a scan relative to the stimulus onset (Fig. 6B, left panel). In the response-locked alignment, the 18 time bins were made to have their centers at −4−RT₀, −2−RT₀, −RT₀, …, 0, …, 28−RT₀, and 30−RT₀ (s) for trials of an RT group with center RT₀ (s), in a manner of stimulus-locked coordinates shifted by RT₀ (Fig. 6B, right panel). However, the assignments of data points to the bins were not the same as the simply-RT₀-shifted stimulus-locked assignments because the response-locked time bins depended on both the acquisition times of data points and RT of the trial. For example, consider a trial with an RT of 7.8 s, which is grouped into RT₀ of 8 s. The bin of 0 s of the response-locked alignment can include data points with acquisition times ranging from RT − 1 (= 6.8 s) to RT + 1 (= 8.8 s), which can be assigned to bins of 6 and 8 s of the stimulus-locked alignment. By a similar trace, in the case of RT 8.2 s of the same group of RT₀ 8 s, data points of the same bin of 0 s of the response-locked alignment can be assigned to bins of 8 s and 10 s of the stimulus-locked alignment. Generally, a bin of the response-locked alignment can include data points from three consecutive stimulus-locked bins, reflecting temporal relations between the signal acquisitions of the voxels and the response moments more precisely than the simply-RT₀-shifted stimulus-locked alignment. Therefore, employing both the two kinds of alignments seems to be effective to obtain accurate properties of activity profiles. The averaged signal segments were normalized to have values between 0 and 1, as shown in Figure 6B. We refer to this “normalized average deconvolved PSC” as “normalized activity” for short, because the deconvolved PSC is the estimated neural activity.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

From PSC time series to profile parameters. A, An example of time series of PSC (blue line) and its deconvolved signal (red line) obtained from the left anterior insula of the fourth session of participant #05. The deconvolved signal was segmented in a trial-based manner, and each segment was labeled by stimulus, participant's first answer and RT of the trial. Segments of trials with correct answers and the assigned-range RTs (indicated in red letters) were averaged in the RT groups (indicated by small disks with colors in common with panel B). B, The stimulus-locked averages (left) and the response-locked averages (right) of the deconvolved PSC segments of each RT group across sessions and participants. (For CSRs, the averages were also separated by the region's selectivity to stimuli.) The averaged signals were normalized to have values between 0 and 1 to highlight temporal profiles rather than magnitudes of activity. For classification of activity patterns, we obtained three kinds of parameters to characterize the profiles: Peak_stim and Peak_rsp, which were the peak times of time courses of the two alignments, and Slope_rsp, which was the slope of the line fitted to four points (at signal values of 0.5, 0.6, 0.7, and 0.8) on the rising side of response-locked time course. The black thick line segment indicates the fitted line as an example. We also obtained Rise_stm, which was the time of the rising side of a stimulus-locked time course reached at a signal value of 0.3, to use not for classification but for rationalizing the classification in Results.

To classify the activity profiles of regions, we basically followed the approach of Ploran et al. (2007), whereas we also took a heuristic approach when it was necessary. We initially postulated peak and rise times of activity signals as the parameters to characterize the activity profiles because the region classes of Ploran et al. (2007) were defined by dependence of these parameters on the decision moment: both of the two parameters show the dependence in the moment-of-recognition regions, only the peak time in the accumulators, and neither in the sensory processors. We also evaluated validity of the present method by reproducibility of the classification results of Ploran et al. (2007, 2011), the latter of which was an improved version of the former study showing substantially similar classification results with modified methods of stimulation and clustering. [In this study, our approach basically followed Ploran et al. (2007), while we referred to the results of Ploran et al. (2011) as well.] In a heuristic manner, we found that the rising slope was more effective for reproducing the classification results of Ploran et al. (2007, 2011) than the rise time as described in Results, so that we used the rising slope of the response-locked time course (referred to as Slope_rsp) instead of the rise time. The three profile parameters, which were the peak times of the stimulus-locked and response-locked time courses (referred to as Peak_stm and Peak_rsp, respectively) and Slope_rsp, were extracted from each of the RT-group-averaged signal segments, as shown in Figure 6B. For Slope_rsp, four points on the rising side of each response-locked time course were obtained at heights of 0.5, 0.6, 0.7, and 0.8 by interpolating between the consecutive sampling points that sandwiched these heights, and Slope_rsp was defined as the slope of the line fitted to these four points. As shown in Figure 6B, the rise time of the stimulus-locked time course (Rise_stm) was defined as the time of the rising side of the time course at a height of 0.3, which was obtained by interpolating between the consecutive sampling points that sandwiched 0.3 in height. Rise_stm was not used for the classification but to rationalize the present method. To reproduce similar results as the classification of Ploran et al. (2007, 2011), we used four statistical quantities of the profile parameters over RTs as follows: SD of Peak_stm (Peak_stm_SD), SD of Peak_rsp (Peak_rsp_SD), mean of Slope_rsp (Slope_rsp_MN), and mean of Peak_rsp (Peak_rsp_MN). The first three quantities were used to reflect the different dependence on decision moment (i.e., RT), while the fourth was employed in a heuristic manner as described in Results.

We classified the regions by means of a hierarchical cluster analysis in the four-dimensional space of the above four quantities using Ward's method (Ward, 1963) following Ploran et al. (2007). In this method, at the initial step, each data point is regarded as a cluster comprising one member; in the next step, a new cluster is created by merging a pair of clusters that minimizes the increase in the sum of the squared Euclidean distances of all the member data points of the new cluster from their centroid (average coordinates). This cluster-merging step is repeated until all data points are included in one cluster to form a clustering tree (dendrogram) that represents the proximity among data points in a hierarchical manner. Because Ward's method uses the Euclidean distances in the multidimensional space of the parameters, the clustering results can be affected by the absolute scales of the parameters. To evaluate the parameters with equal weights, we scaled the SD of the values of each parameter over the regions to 1 by dividing by their original SD before applying Ward's method. In this study, we used Ward's algorithm implemented in MATLAB Statistics and Machine Learning Toolbox (version 11.0; function linkage).

In the data analyses of this study, we used applications including SPM12 and homemade programs working on MATLAB 2016b. Statistical tests were performed using SPSS version 21 (IBM Corp.) and G*Power (Faul et al., 2007). Visualization of brain mapping in Figures 3, 6, 12 was conducted with the BrainNet Viewer (version 1.53; Xia et al., 2013; http://www.nitrc.org/projects/bnv/).