Materials and Methods
Experimental setup and procedure
Participants learned the monetary value of 12 novel visual stimuli in a reinforcement learning paradigm. Learning occurred over the course of four MRI scan sessions conducted on four consecutive days. The novel stimuli were 3-dimensional shapes generated with a custom built MATLAB toolbox (code available at
http://github.com/saarela/ShapeToolbox). ShapeToolbox allows the generation of three-dimensional radial frequency patterns by modulating basis shapes, such as spheres, with an arbitrary combination of sinusoidal modulations in different frequencies, phases, amplitudes, and orientations. A large number of shapes were generated by selecting combinations of parameters at random. From this set, we selected twelve that were considered to be sufficiently distinct from one another. A different monetary value, varying from $1.00 to $12.00 in integer steps, was assigned to each shape. These values were not correlated with any parameter of the sinusoidal modulations, so that visual features were not informative of value.
Participants completed 20 minutes of the main task protocol on each scan session, learning the values of the 12 shapes through feedback. The sessions were comprised of three scans of 6.6 minutes each, starting with 16.5 seconds of a blank gray screen, followed by 132 experimental trials (2.75 seconds each), and ending with another period of 16.5 seconds of a blank gray screen. Stimuli were back-projected onto a screen viewed by the participant through a mirror mounted on the head coil and subtended 4 degrees of visual angle, with 10 degrees separating the center of the two shapes. Each presentation lasted 2.5 seconds and, at any point within a trial, participants entered their responses on a 4-button response pad indicating their shape selection with a leftmost or rightmost button press. Stimuli were presented in a pseudorandom sequence with every pair of shapes presented once per scan.
Feedback was provided as soon as a response was entered and lasted until the end of the stimulus presentation period. Participants were randomly assigned to two groups depending on the type of feedback received. In the RELATIVE feedback case, the selected shape was highlighted with a green or red square, indicating whether the selected shape was the most valuable of the pair or not, respectively. In the ABSOLUTE feedback case, the actual value of the selected shape (with variation) was displayed in white font. Between each run, both groups received feedback about the total amount of money accumulated up to that point. The experimental protocol is further described in detail here \cite{Mattar_2016}.
MRI Data collection and preprocessing
Magnetic resonance images were obtained at the Hospital of the University of Pennsylvania using a 3.0 T Siemens Trio MRI scanner equipped with a 32-channel head coil. T1-weighted structural images of the whole brain were acquired on the first scan session using a three-dimensional magnetization-prepared rapid acquisition gradient echo pulse sequence with the following parameters: repetition time (TR) 1620 ms, echo time (TE) 3.09 ms, inversion time 950 ms, voxel size 1 mm by 1 mm by 1 mm, and matrix size 190 by 263 by 165. To correct geometric distortion caused by magnetic field inhomogeneity, we also acquired a field map at each scan session with the following parameters: TR 1200 ms, TE1 4.06 ms, TE2 6.52 ms, flip angle 60\({}^{\circ}\), voxel size 3.4 mm by 3.4 mm by 4.0 mm, field of view 220 mm, and matrix size 64 by 64 by 52. In all experimental runs with a behavioral task, T2*-weighted images sensitive to blood oxygenation level-dependent contrasts were acquired using a slice accelerated multi-band echo planar pulse sequence with the following parameters: TR 2000 ms, TE 25 ms, flip angle 60\({}^{\circ}\), voxel size 1.5 mm by 1.5 mm by 1.5 mm, field of view 192 mm, and matrix size 128 by 128 by 80. In all resting state runs, T2*-weighted images sensitive to blood oxygenation level-dependent contrasts were acquired using a slice accelerated multi-band echo planar pulse sequence with the following parameters: TR 500 ms, TE 30 ms, flip angle 30\({}^{\circ}\), voxel size 3.0 mm by 3.0 mm by 3.0 mm, field of view 192 mm, matrix size 64 by 64 by 48.
Cortical reconstruction and volumetric segmentation of the structural data was performed with the Freesurfer image analysis suite
\cite{dale1999cortical}. Boundary-Based Registration between structural and mean functional image was performed with Freesurfer bbregister
\cite{greve2009accurate}. Preprocessing of the fMRI data was carried out using FEAT (FMRI Expert Analysis Tool) Version 6.00, part of FSL (FMRIB’s Software Library,
www.fmrib.ox.ac.uk/fsl). The following pre-statistics processing was applied: EPI distortion correction using FUGUE
\cite{jenkinson2004improving}, motion correction using MCFLIRT
\cite{jenkinson2002improved}, slice-timing correction using Fourier-space time series phase-shifting, non-brain removal using BET
\cite{smith2002fast}, grand-mean intensity normalization of the entire 4D dataset by a single multiplicative factor, and highpass temporal filtering via Gaussian-weighted least-squares straight line fitting with
\(\sigma=50.0s\). Nuisance time series were voxelwise regressed from the preprocessed data. Nuisance regressors included (i) three translation (X, Y, Z) and three rotation (pitch, yaw, roll) time series derived by retrospective head motion correction (R = [X, Y, Z, pitch, yaw, roll]), together with the first derivative and square expansion terms, for a total of 24 motion regressors
\cite{friston1996movement}); (ii) the first five principal components of non-neural sources of noise, estimated by averaging signals within white matter and cerebrospinal fluid masks, obtained with Freesurfer segmentation tools, and removed using the anatomical CompCor method (aCompCor)
\cite{behzadi2007component}; and (iii) a measure of a local source of noise, estimated by averaging signals derived from the white matter region located within a 15 mm radius from each voxel, using the ANATICOR method
\cite{jo2010mapping}. Global signal was not regressed out of voxel time series
\cite{chai2012anticorrelations,murphy2009impact}.
We parcellated the brain into 112 cortical and subcortical regions, separated by hemisphere using the structural Harvard-Oxford atlas of the FMRIB (Oxford Centre for Functional Magnetic Resonance Imaging of the Brain) Software Library
(FSL; Version 5.0.4) \cite{smith2004advances,woolrich2009bayesian}. We warped the MNI152 regions into subject-specific native space using FSL FNIRT and nearest-neighbor interpolation and calculated the average BOLD signal across all gray matter voxels within each region. The participant’s gray matter voxels were defined using the anatomical segmentation provided by Freesurfer, projected into subject’s EPI space with bbregister. For each individual scan, we extracted regional mean BOLD time series by averaging voxel time series in each of the 112 regions of interest.
Statistical Methods: Inter-subject correlation (ISC) and inter-subject functional connectivity (ISFC)
We used the Pearson correlation coefficient to compute the inter-subject correlation (ISC) of each brain region for each subject. First, we calculated the region-wise temporal correlation between every pair of subjects as:
\begin{equation}
r_{ij}=\frac{\sum_{t=1}^{N}[(x_{i}(t)-\bar{x}_{i})(x_{j}(t)-\bar{x}_{j})]}{\sqrt{\sum_{t=1}^{N}(x_{i}(t)-\bar{x}_{i})^{2}\sum_{t=1}^{N}(x_{j}(t)-\bar{x}_{j})^{2}}},\par
\\
\end{equation}
where \(N\) is the number of time points of time series data, and where \(r_{ij}\) is the correlation coefficient of a region between the times series \(x_{i}\) and \(x_{j}\) of the \(i\)th and \(j\)th subjects, respectively.
To test the statistical significance of the correlation between the fMRI BOLD signals of a single region from two subjects, we performed a fully non-parametric permutation test with 10,000 randomizations \cite{manly2006randomization}. This test accounts for slow-scale autocorrelation structure in the BOLD time series \cite{handwerker2012periodic} by removing phase information from each BOLD signal, through Fourier phase randomization. We repeated this procedure 10,000 times to obtain a null distribution of the maximum noise correlation values and defined the threshold for a pair correlation as the q*100th percentile of the null distribution of maximum values.
To obtain the ISC, we averaged only significant correlation values out of 190 correlation values, \(r_{ij}\) from all subject pairs, to obtain one ISC for each region\(:\)
\begin{equation}
ISC=\frac{2}{M(M-1)}\sum_{i=1}^{M}\sum_{j=2,j>1}^{M}r_{ij},\par
\\
\end{equation}
where M is the number of subjects, which in this study was 20. The procedure (see SI for more information) was followed for each scan, for each day, and for both rest and task conditions.
We constructed functional brain networks for each experimental session of each subject, by computing a wavelet transform coherence (WTC) measure between every possible pair of regional BOLD time series \cite{torrence1998practical}. To calculate the WTC, we first applied a wavelet decomposition to the BOLD time series of each brain region by successively convolving the time series with the Morlet wavelet in a frequency range of 0.05 – 0.11Hz (see the SI for further information). The wavelet transform of each time series was smoothed in both time and frequency \cite{cazelles2007time,grinsted2004application}. Next, we computed the correlation between the smoothed wavelet transformed signals associated with each possible pair of brain regions. We repeated this procedure for all pairs of regions and tabulated the resulting WTC values in a weighted adjacency matrix A reflecting the functional connectivity for each experimental session and for each subject.
Inter-subject functional connectivity (ISFC) was obtained from the functional connectivity (FC) matrices of all the subjects, and can be thought of as an estimate of the correlation in FC between one subject and all other subjects. Specifically, we computed the ISFC of each subject as the correlation between the single subject FC matrix and the average of all other subjects FC matrices as
\begin{equation}
ISFC_{i}=\frac{1}{N}A_{i}[\frac{1}{n-1}\sum_{j\neq i}A_{j}]\par
\\
\end{equation} ,
where A is the functional connectivity matrix of a subject, and where n and N are the total number of subjects and total number of regions, respectively.
To objectively identify a common set of functional brain systems across all subjects, we partitioned the adjacency matrix into network modules by optimizing a multilayer modularity quality index using a Louvain-like \cite{blondel2008fast} locally greedy algorithm \cite{fortunato2010community, mucha2010community}. We constructed a multi-layer network from our experimental data, by concatenating adjacency matrices from all sessions and subjects for each day. Using this approach, we obtained a partition of brain regions into modules for each scan and for each subject (See Supplementary Information for details). To obtain a single representative partition of brain regions into distinct communities, we computed a module allegiance matrix \cite{bassett2011dynamic,Bassett2013,Bassett2015}, whose ijth entry represents the probability that region i and region j belong to the same community across scans and participants. We then applied single-layer modularity maximization to this module allegiance matrix to obtain a single partition of brain regions into consensus modules \cite{Bassett2013}.