intro to working with volume and surface brain data

When preparing to update my surface (gifti) tutorial knitr, I realized that some of the confusion I was trying to address wasn't due to the code so much as the general differences between volume and surface datasets. This post tries to fill in some of those gaps for people using my R scripts. Much of the information in this post is general, but there are many topics I am not going to cover at all - what I describe here is not the only way to work with fMRI brain data; there are many other file formats and software packages.

first, volumes.

Volumetric brain data files are first, both historically (all fMRI analyses were originally volumetric) and fundamentally: we (including our brains) are physically in three dimensions, as is the data that comes off the scanner (often as DICOM files, which neuroimagers usually promptly convert to NIfTI, such as via dcm2niix when making a BIDS dataset).

Consider a standard anatomical image, such as the HCP 1200 Group Average, which makes a good default underlay for adult MNI datasets (S1200_AverageT1w_81x96x81.nii.gz is a downsampled version used in my tutorials). Opening S1200_AverageT1w_81x96x81.nii.gz with MRIcron shows us slices  for each of the three planes (below left), and the schematic at right is how the data comes out of a NIfTI image:

This image's data is a big cube subdivided into a grid of little cubes (voxels; volumetric pixels), each 2.4 mm on a side. It has 81 voxels along the i axis, 96 along the j, and 81 along the k, for 81*96*81 = 629,856 voxels total. There is a number in every one of these 629,856 voxels: the voxels fill the entire cube, not just the brain part in the middle. Thus, when you plot slices through a NIfTI it hopefully makes a brain-shaped picture, but that's not guaranteed - you can put any shape you want into a NIfTI file. (And hence why you should always include a "does it look like a brain" QC step.)

The code below opens the NIfTI in R with the splendid RNifti package and shows us the same information as the MRIcron screen capture above: the voxel under the cursor (at 25, 56, 37; red scribble) has value 605.887 (green scribble). (Warning: these match since both R and MRIcron are 1-based; if using something 0-based like afni or python you will need to shift indices.)

 library(RNifti); # package for NIfTI image reading https://github.com/jonclayden/RNifti  
   
 fname <- paste0(in.path, "S1200_AverageT1w_81x96x81.nii.gz");  # https://osf.io/6hdxv/  
 if (file.exists(fname)) { under.img <- readNifti(fname); } else { stop(paste("missing:", fname)); }  
   
 str(under.img)  
 # 'niftiImage' num [1:81, 1:96, 1:81] 0 0 0 0 0 0 0 0 0 0 ...  
 # - attr(*, "pixdim")= num [1:3] 2.4 2.4 2.4  
 # - attr(*, "pixunits")= chr [1:2] "mm" "Unknown"  
 # - attr(*, ".nifti_image_ptr")=<externalptr>   
 # - attr(*, ".nifti_image_ver")= int 2  
   
 under.img[25,56,37]  
 # [1] 605.887  

Note that under.img is a 3d array of numbers with various attributes; its pixdim and pixunits match those in the MRIcron "NIfTI Header Information" dialog in the screen capture (open with Window -> Information). There are many more NIfTI header fields; the header information is what makes it possible for software like MRIcron to match the 3d array of numbers in the NIfTI with actual space. 

We don't need to implement every NIfTI header field or work out warping and realigning on our own, thankfully! The key is to think of the volumetric images as a 3d (or 4d, if functional timeseries - lots of 3d images stuck together along the 4th dimension) array of numbers with important properties attached. If  you have multiple NIfTI files you can treat them as normal arrays of numbers, so long as all the files have matching properties (same voxel sizes, orientation, etc. - compatible headers). Thankfully, all the main fMRI software programs have functions for making image properties match; here's my afni example.

why surfaces?

I wrote above that neuroimaging data (and our brains!) start in 3d volume space; why would we want to think about it any other way? Well, when doing fMRI the goal is often to analyze the cortical grey matter BOLD signal, and we know the human cortex is physically a folded sheet (cortical "ribbon"); a convoluted surface.

Figure 2 from doi:10.1007/s00429-012-0411-8, PMID: 22488096 (below) shows a cortical surface (lots more info in the paper and elsewhere); the blue scribble I added to pane c follows the grey matter ribbon.

In principle, if the preprocessing algorithms find this grey matter ribbon we can concentrate our analyses there and have better precision and stronger results. There are a lot of complications in finding and transforming the grey matter ribbon of any particular person into a surface (a few of my musings; Andrew Jahn's), but we won't explore any of those here. The explanations below use standard images from FreeSurfer; to get a surface version of a new dataset I recommend preprocessing with software (e.g., fmriprep) which will output gifti files (or whatever format you want) directly (don't, for example, preprocess completely as volumes and then convert to surfaces for analysis).

introduction to gifti files

A single NIfTI file is enough for an entire brain, but not a single gifti file: gifti files come in paired sets (see also). The "paired" is because each gifti file has info for one hemisphere only: for the whole cortex you need a pair of files, one left and one right. And the "sets" is because gifti files come in multiple types, labeled by the bit of the filename before the .gii. Here we'll only use two types: .surf.gii surface "underlay" images, and .label.gii parcel labels.

For my gifti knitr tutorial I shared several sets of paired gifti files; here we'll look at the Schaefer2018 parcellation fsaverage5 ones in R, using the gifti library's readGIfTI function. (I converted the original files to gifti using FreeSurfer's mris_convert function.)

First, read four files: a pair each of .surf.gii and .label.gii:

 library(gifti);  # for readGIfTI  https://github.com/muschellij2/gifti  
   
 # surface "underlay anatomy"  
 fname <- paste0(in.path, "fsaverage5_lh.inflated_avg.surf.gii");  # https://osf.io/nxpa7/  
 if (file.exists(fname)) { surf.L <- readGIfTI(fname); } else { stop("can't find fsaverage5_lh.inflated_avg.surf.gii"); }  
 fname <- paste0(in.path, "fsaverage5_rh.inflated_avg.surf.gii");  # https://osf.io/24y9q/  
 if (file.exists(fname)) { surf.R <- readGIfTI(fname); } else { stop("can't find fsaverage5_rh.inflated_avg.surf.gii"); }  
    
 # Schaefer parcellation "overlay" surfaces for each hemisphere   
 fname <- paste0(in.path, "Schaefer2018_400Parcels_7Networks_order_10K_L.label.gii");  # https://osf.io/7wk8v/  
 if (file.exists(fname)) { parcels.L <- readGIfTI(fname); } else { stop(paste("missing:", fname)); }  
 fname <- paste0(in.path, "Schaefer2018_400Parcels_7Networks_order_10K_R.label.gii");  # https://osf.io/7tms8/  
 if (file.exists(fname)) { parcels.R <- readGIfTI(fname); } else { stop(paste("missing:", fname)); }   
   

The next code looks a bit at the list objects made with readGIfTI(). For plotting we only need the $data fields, but I suggest investigating the other entries as well to gain a sense of the type of information in them.

 str(surf.L$data);  
 # List of 2  
 # $ pointset: num [1:10242, 1:3] -0.685 22.359 41.099 9.658 -36.773 ...  
 # $ triangle: int [1:20480, 1:3] 0 0 0 0 0 3 3 3 4 4 ...  
   
 str(surf.R$data)  
 # List of 2  
 # $ pointset: num [1:10242, 1:3] -5.74 8.43 43.5 12.6 -40.33 ...  
 # $ triangle: int [1:20480, 1:3] 0 0 0 0 0 3 3 3 4 4 ..  
   
 str(parcels.L$data)  
 # List of 1  
 # $ labels: int [1:10242, 1] 57 81 49 180 40 132 14 22 0 143 ...  
   
 str(parcels.R$data)  
 # List of 1  
 # $ labels: int [1:10242, 1] 61 81 46 110 111 200 21 95 104 119 ...  

Notice that the two surf objects each have two subpart arrays ($pointset and $triangle), while each parcel object only has a vector of integers. The parcel vectors' lengths are 10,242, which matches the number of rows in the surface pointset arrays; for plotting, the gifti data sizes must match like this. Why 10,242? That's the number of vertices per hemisphere in fsaverage5; other surface spaces have a different number of vertices.

Triangles and vertices? These terms come from how surface meshes can be represented in computer science: points (vertices, orange arrow) connected by edges of specific lengths, making lots of little triangles, as in this part of the FreeSurfer website brain I enlarged:

using and plotting giftis

The key idea when working with gifti images using my scripts is that we manipulate the vectors of vertex-wise data, not the .surf.gii files (nor any other spatial information). In a given surface space (here, freesurfer5) the location of each vertex as ordered in the 10,242-element vectors is fixed, so we can manipulate its value as needed, as shown in the next examples.

This code produces the image below. The left hemisphere is blank: all its vertices (vals.L) are 0. For the right I again started with a vector of 10,242 0s, but then changed its 1,800th entry to 1, producing a brain picture with a red splotch at the location of vertex 1800. 

 vals.L <- rep(0, 10242); # start with all vertices 0  
 map.L <- gifti.map(surf.L.infl$data$pointset, surf.L.infl$data$triangle, vals.L); # map vals.L onto left surface  
   
 vals.R <- rep(0, 10242); # start with all vertices 0  
 vals.R[1800] <- 1; # change vertex 1800 to 1  
 map.R <- gifti.map(surf.R.infl$data$pointset, surf.R.infl$data$triangle, vals.R);  # map vals.R onto right surface  
   
 # plot lateral and medial of both hemispheres  
 plot.surface(map.L, ttl="left hemisphere, blank (all vertices 0)", pos.lims=c(0.5,5.5), which.surface='fsaverage5');   
 plot.surface(map.R, ttl="right hemisphere, vertex 1800 in red", pos.lims=c(0.5,5.5), which.surface='fsaverage5');   

And a final example (more are in the knitr tutorial and other demo code): the next code plots (Schaefer 400x7) parcel 60 on the left hemisphere, and all 200 parcels on the right hemisphere. How to know which of the 10,242 vertices correspond to left hemisphere parcel 60? The key is parcels.L (created above from Schaefer2018_400Parcels_7Networks_order_10K_L.label.gii), which has a 10,242-element vector of parcel labels: we just need to find the 60s.

 vals.L <- rep(0, 10242); # start with all vertices 0  
 inds <- which(parcels.L$data[[1]] == 60); # find indices; where 60s are in the 10,242-element vector parcels5.L$data  
 vals.L[inds] <- 4;  # change indices in vals.L for parcel 60 to 4  
 map.L <- gifti.map(surf.L.infl$data$pointset, surf.L.infl$data$triangle, vals.L);  # map vals.L onto left surface  
 map.R <- gifti.map(surf.R.infl$data$pointset, surf.R.infl$data$triangle, parcels.R$data[[1]]); # map all parcels.R vertices onto right surface  
   
 # plot lateral and medial of both hemispheres  
 plot.surface(map.L, ttl="left hemisphere, parcel 60 in orange", pos.lims=c(0.5,5.5), which.surface='fsaverage5'); # plot lateral and medial, left  
 plot.surface(map.R, ttl="right hemisphere, parcels 1 to 200", pos.lims=c(0.5,200.5), which.surface='fsaverage5'); # plot lateral and medial, right  

"crescent" artifacts revisited

Back in 2018 I posted twice about "crescent" artifacts (first, second), and have kept a casual eye out for them since, wondering how much they might affect fMRI analysis results. 

The artifact isn't totally consistent across people, but when present, is very stable over time (i.e., sessions several years apart), and bright in temporal standard deviation QC images. The artifact's location varies with encoding direction: front for PA encoding, rear with AP encoding; the PA encoding versions tend to be much brighter and obvious than the AP. 

Below is an example of the artifact, pointed out by the pink arrows. This is the temporal standard deviation image for sub-f8570ui from the DMCC55B dataset (doi:10.18112/openneuro.ds003465.v1.0.6), the first (left, AP encoding) and second (right, PA encoding) runs of Sternberg baseline (ses-wave1bas_task-Stern), after fmriprep preprocessing:

These two runs were collected sequentially within the same session, but the artifact is only visible in the PA encoding run (right). (Briefly, DMCC used a 3T Siemens Prisma, 32-channel headcoil, CMRR MB4, no in-plane acceleration, 2.4 mm iso, 1.2 s TR, alternating AP and PA encoding runs; details at openneuro, dataset description paper, and OSF sites, plus DMCC-tagged posts on this blog.)

In the previous "crescent" posts we speculated that these could be N/2 ghosts or related to incomplete fat suppression; I am now leaning away from the Nyquist ghost idea, because the crescents don't appear to line up with the most visible ghosts. (Some ghosts are a bit visible in the above image; playing with the contrast and looking in other slices makes the usual multiband-related sets of ghosts obvious, but none clearly intersect with the artifact.) It also seems odd that ghosts would be so much brighter and change their location with AP vs. PA encoding; I am no physicist, though!

link to cortex volume?

This week I gave three lab members a set of temporal standard deviation images (similar to the pair above) for 115 participants from the DMCC dataset. The participant images were in random order, and I asked the lab members to rate how confident they were that each showed the "crescent" artifact or not. My raters agreed that 34 participants showed the artifact, and 39 did not. (Ratings were mixed or less confident on the others; I asked them to decide quickly from single pairs of QC images, not investigate closely.)

We didn't measure external head size in the participants, but did run freesurfer during preprocessing, so I used its CortexVol and eTIV statistics as proxies (a different stat better?): and the group my colleagues rated as having the artifact tended to have smaller brains than those without:

If the appearance of this artifact is indeed somewhat related to head size, then it's logical that it would (as I've observed) generally be stable over time. DMCC's population was younger adults; it'd be interesting to see if there's a relationship with a wider range of head sizes.

only with DMCC or its MB4 acquisition?

Is the artifact restricted to this particular acquisition or study? No, not somehow related to DMCC; I've checked a few DMCC participants with the artifact who later participated in other studies and they have it (or not) in all of the datasets.

To see if it's restricted to the MB4 acquisition, I looked at a few images from a different study, which also has adult participants, a 3T Prisma, 32-channel headcoil, 2.4 mm iso voxels, and CMRR sequences, but with MB6 TR 0.8 s, and PA encoding for all runs. Below are standard deviation images for three different people from this MB6 study, one run each of the same task, after fmriprep preprocessing. (I chose these three because of the artifacts; not all are so prominent.)

Since this study has all PA runs I can't directly compare artifacts across the encoding directions, but there are clearly some "crescents", and more sets of rings than typical with MB4 (which makes sense for MB6). The rings are especially obvious in person #3; some of these appear to be full-brain ghosts. I suspect the artifacts would be much clearer in subject space; I haven't looked (I'm not directly involved in the study). But a substantial minority of these participants' standard deviation images resemble #1, whose artifacts strike me as quite similar to the "crescents" in some DMCC MB4 PA images.

but does it matter?

Not all artifacts that look strange in QC images actually change task-related BOLD statistics enough to be a serious concern. (Of course, how much is too much totally depends on the particular case!) I suspect that this artifact does matter for our analyses, though, both because of where it falls in the brain and because it affects BOLD enough to be visible by eye in some cases.

The artifact's most prominent frontal location with PA runs puts it uncomfortably close to regions of interest for most of my colleagues, and is one reason I have advised shifting to all AP encoding for new studies I'm involved with. Preprocessing, motion, transformation to standard space, and spatial smoothing blurs the artifact across a wider area, hopefully diluting its effect. But the artifact's location is somewhat consistent across participants, and present in a sizable enough minority (a third, perhaps, in the datasets I've looked at), that it seems possible it could reduce signal quality in our target ROIs.

For showing that it does indeed actually affect task-related BOLD enough to matter, so far I mostly just have qualitative impressions. For example, below left is the standard deviation of one DMCC person's PA Sternberg run, with the cursors on the artifact. The right side is from averaging together (after voxel-wise normalizing and detrending) frames after pressing a button with the right hand. Squinting, the statistical image is brighter in sensible motor-related grey matter areas, marked with green. But the "crescent" may also be faintly visible, as pointed out in pink.

I can imagine quantitative tests, such as comparing the single-run (so separating AP and PA encoding runs) GLM output images from the group of participants with and without the artifact. Differences in estimates in parcels/searchlights/areas overlapping the artifact would be suggestive, particularly as the estimates vary with encoding direction and participant subgroup (with-artifact or without).

thoughts?

I'm curious, have you seen similar? Do you think this artifact is from N/2 ghosting, incomplete fat suppression, or something else? (What should I call it? "Crescent" is visually descriptive, but not standard. 😅) Seem reasonable it could be related to head size? And that it can significantly affect BOLD? Other reactions? Thanks! (And we can chat about this at my OHBM 2025 poster, which will be on this topic.)

tutorial: making new versions of NIfTI parcellation files

This post is a tutorial on making new versions of NIfTI files, particularly parcellation images. For example, suppose you ran an analysis using a set of parcels from the Schaefer 400x7 volumetric parcellation, and you want to make a figure showing just those parcels. One way to make such a figure is to write a new NIfTI image the same size and shape as the full parcellation, but with only the parcels you want to show, and then display it as an overlay in a visualization program or via code.

Or, suppose you want to make a parcel-average timecourse using afni's 3dROIstats function, but of a group of parcels together, not each parcel individually. If you make a new mask image with the same number to all the voxels in all the parcels you want treated as a group, 3dROIstats will average them together into a single timecourse.

The final example assigns a new value to each parcel, such as to show the t-values resulting from a statistical test run on each parcel individually.

The code below uses the same files as my NIfTI knitr tutorial, Schaefer2018_400x7_81x96x81.nii.gz  and  S1200_AverageT1w_81x96x81.nii.gz. To run the code, save the two files somewhere locally, and then change path to their location.

new NIfTI with one parcel only

This first block of code sets the variables for all examples, and then makes a new NIfTI file, p24.nii.gz, which has all the voxels in parcel 24 set to 1, and all the other voxels set to 0. 

 library(RNifti); # package for NIfTI image reading https://github.com/jonclayden/RNifti  
   
 rm(list=ls()); # clear memory  
 path <- "//storage1.ris.wustl.edu/tbraver/Active/MURI/BraverLab/JoWorking/tutorial/"; # path to input images  
 p.img <- readNifti(paste0(path, "Schaefer2018_400x7_81x96x81.nii.gz"));  # read the parcellation image  
   
 # make a new nifti with only parcel 24   
 new.img <- array(0, dim(p.img));  # blank "brain" of zeros, same size as p.img  
 new.img[which(p.img == 24)] <- 1;  # voxels with value 24 in p.img set to 1 
 writeNifti(new.img, paste0(path, "p24.nii.gz"), template=paste0(path, "Schaefer2018_400x7_81x96x81.nii.gz"));  
 # note the template= in the writeNifti function: it specifies that the new file should have the same 
 # header settings as the original parcellation image; important for both to match properly.  

Below are two MRIcron windows, both with S1200_AverageT1w_81x96x81.nii.gz as the anatomical underlay. At left is p24.nii.gz; the title bar text giving the value of the voxel under the crosshairs as 1. The right pane overlay is Schaefer2018_400x7_81x96x81.nii.gz, and the voxel under the crosshairs has the value 24, as expected.

new NIfTI with several parcels given the same value

This code is very similar to the first, but instead of only setting the voxels in parcel 24 to 1, it changes the values of the voxels in all the parcels listed in p.ids to 1 (effectively, making a new parcel out of the previous four):

 # make a new nifti with parcels 24, 26, 370, and 231 all given the value 1.  
 p.ids <- c(24, 26, 231, 360);    # parcel ids we want to combine  
 
 new.img <- array(0, dim(p.img));  # blank "brain" of zeros, same size as p.img  
 for (i in 1:length(p.ids)) { new.img[which(p.img == p.ids[i])] <- 1; }  # set voxels in p.ids to 1   
 writeNifti(new.img, paste0(path, "four.nii.gz"), template=paste0(path, "Schaefer2018_400x7_81x96x81.nii.gz"));  
   

new NIfTI with unique numbers for each parcel

To give different numbers to each parcel, loop over all the parcels and assign the corresponding value:

 # assign a value to each parcel, such as from the results of a statistical test performed on each parcel individually.  
 # note: 400 is the total number of parcels in the example parcellation. It is usually best not to hard-code  
 # the number of parcels, but it is here to keep the code as short as possible.  
 stat.vals <- rnorm(400); # random numbers to plot, one for each parcel
   
 new.img <- array(0, dim(p.img));  # blank "brain" of zeros, same size as p.img  
 for (i in 1:400) { new.img[which(p.img == i)] <- stat.vals[i]; }  # set voxels in parcel i to stat.value i  
 writeNifti(new.img, paste0(path, "stat.nii.gz"), template=paste0(path, "Schaefer2018_400x7_81x96x81.nii.gz"));  

Below are views of four.nii.gz (left) and stat.nii.gz (right). Since the vector of "statistics" for each parcel in the example code is random, it will be different each time, but the extreme positive and negative values were plotted below in hot and cool colors as is usual for fMRI statistical images. 

Corresponding Schaefer2018 400x7 and 400x17 atlas parcels by number

My "default" cortical parcellation is the 400 parcels by 7 networks version of Schaefer2018. I like these parcellations because they're available in all the spaces we use (volume, fsaverage, fsLR (HCP)) and are independent from our analyses; "all parcellations are wrong, but some are useful".

The Schaefer parcellations come in several versions for each resolution, however, and a parcel described by its 7Networks number likely has a different 17Network number. Note that the parcel boundaries are the same for all same-resolution versions: there is only one set of 400 parcels, but which of those 400 parcels is #77 varies between the 7 and 17 network versions. This post describes how to translate parcel numbers between network versions, using a bit of (base) R code.

Logic: since there is only one set of parcels at each resolution, there is only one set of centroid coordinates at each resolution. Thus, we can match parcels across network orderings by centroids.

First, set sch.path to the location of your copy of the Schaefer2018 Parcellations directory and load the 7 and 17 network files:  (apologies for the wonky code formatting)

sch.path <- "/data/nil-bluearc/ccp-hcp/DMCC_ALL_BACKUPS/ATLASES/Schaefer2018_Parcellations/";

cen7.tbl <- read.csv(paste0(sch.path, "MNI/Centroid_coordinates/Schaefer2018_400Parcels_7Networks_order_FSLMNI152_1mm.Centroid_RAS.csv"));
cen17.tbl <- read.csv(paste0(sch.path, "MNI/Centroid_coordinates/Schaefer2018_400Parcels_17Networks_order_FSLMNI152_1mm.Centroid_RAS.csv"));

Next, make vectors for translating the 7Network number to the 17Network number (and the reverse):

x7to17 <- rep(NA, nrow(cen7.tbl));
x17to7 <- rep(NA, nrow(cen7.tbl));
for (i in 1:nrow(cen7.tbl)) { 
  x7to17[i] <- which(cen17.tbl$R == cen7.tbl$R[i] & cen17.tbl$A == cen7.tbl$A[i] & cen17.tbl$S == cen7.tbl$S[i]); 
  x17to7[i] <- which(cen7.tbl$R == cen17.tbl$R[i] & cen7.tbl$A == cen17.tbl$A[i] & cen7.tbl$S == cen17.tbl$S[i]); 
}

Now the vectors can be used to translate parcel numbers: parcel #77 in the 7Networks ordering is parcel #126 in the 17Networks ordering.

x7to17[77]   # [1] 126

cen7.tbl[77,]
#   ROI.Label                     ROI.Name   R   A  S
#77        77 7Networks_LH_DorsAttn_Post_9 -33 -46 41

cen17.tbl[126,]
#    ROI.Label                  ROI.Name   R   A  S
#126       126 17Networks_LH_ContA_IPS_5 -33 -46 41

x17to7[126]    # [1] 77

Note that the parcel's R A S coordinates are the same, but its name and label (number) vary between the two network versions.

OHBM 2023 links

I was a co-organizer (with Paul Taylor, Daniel Glen, and Richard Reynolds, all of afni fame) of the "Making Quality Control Part of Your Analysis: Learning with the FMRI Open QC Project" course at OHBM 2023. FMRI Open QC Project participants described how they approached the Project datasets, and we had several general QC discussions. Thanks to all participants and attendees! I hope we can keep making progress towards a common vocabulary for fMRI QC and encourage researchers to include it if they do not already.

The session was unfortunately not recorded live, though speakers submitted recordings in advance. OHBM is supposed to make all of these recordings available; in the meantime I've linked material already public.

• Brendan Williams, "Reproducible Decision Making for fMRI Quality Control"
• Céline Provins, "Quality control in functional MRI studies with MRIQC and fMRIPrep"
• Daniel Glen, "Quality control practices in FMRI analysis: philosophy, methods and examples using AFNI"
• Dan Handwerker, "The art and science of using quality control to understand and improve fMRI data"  slides
• Chao-Gan Yan, "Quality Control Procedures for fMRI in DPABI"
• Rasmus Birn, "Quality control in resting-state functional connectivity: qualitative and quantitative measures"
• Xin Di, "QC for resting-state and task fMRI in SPM"
• Jo Etzel, "Efficient evaluation of the Open QC task fMRI dataset"  video 
• Rebecca Lepping, "Quality Control in Resting-State fMRI: The Benefits of Visual Inspection"
• Francesca Morfini, "Functional Connectivity MRI Quality Control Procedures in CONN"

I also presented a poster, "Which Acquisition? Choosing protocols for task fMRI studies", #700. Here's some of the introduction and conclusion for an abstract. The test data is already public; the code isn't written up properly, but I could share if anyone is interested.

When planning a task fMRI study, one necessary choice is the acquisition sequence. Many are available, with recommendations varying with hardware, study population, brain areas of interest, task requirements, etc.; it is rare to have only one suitable option. Acquisition protocols for task studies can be difficult to evaluate, since metrics like tSNR are not specific to task-related activity. But task designs can make choosing easier, since there is a known effect to compare the candidate acquisition protocols against. 

The procedure illustrated here will rarely make the choice of acquisition completely unambiguous, but can indicate which to avoid, and give the experimenters confidence that the chosen sequence will produce usable data. After choosing the acquisition, more piloting should be performed with the study tasks to confirm that image quality and response clarity are sufficient and as expected.   

... it's taken me so long to finish this post (started August 2023!) that I'm publishing it without adding a proper explanation of the protocol-choosing logic. Take a look at the poster pdf, and please ask questions or otherwise nudge me for more information if you're interested.

US researchers: use extra caution with participant gender-related info

Last summer I wrote several times about pregnancy-related data suddenly becoming much more sensitive and potentially damaging to the participant if released. Unfortunately, now we must add transgender status, biological sex at birth (required, e.g., for GUID creation), gender identity, relationship status, and more to the list of data requiring extra care. 

The NIH Certificate of Confidentiality thankfully means that researchers can't be required to release data on something like someone's abortion or transgender history for criminal or other proceedings. We are responsible for ensuring that our data is handled and stored securely, so that it isn't accidentally (or purposely) shared and then possibly cause the participant harm. I suggest researchers review their data collection with an eye towards information that may be more sensitive now than it was in the past; is this information required? If so, consider how to store and use it securely and safely.

Also consider what you ask potential participants during screening (and how screening is done in general): people may not wish to answer screening questions about things like pregnancy or hormone treatments. If these possibly-sensitive questions must be asked, consider how to do so while minimizing potential discomfort or risk.

For example, one of our studies can't include pregnant people, so we must ask about pregnancy during the initial phone screenings. We used to ask potential participants about pregnancy separately, but then changed the script so that this sensitive question was in a list, and the participant is asked if any apply. This way, the participant doesn't have to state explicitly that they are pregnant, and the researcher doesn't have any specific notes, or even respond in a specific way (e.g., we don't want them to say something like "sorry Jane Doe, but you can't be in our study now because you're pregnant").

Here's the relevant part of the new screening script:

In order to determine your eligibility for our study, I need to ask you some questions. This will take about 15 minutes. 

Before we collect your demographic information, I will ready you a list of four exclusionary criteria. If any of these describe you, please answer yes after I read them all; if none apply, please answer no. You do not need to say which ones apply. 

• You are a non-native English speaker (you learned to speak English as an adult); 
• You are over the age of 45 or under the age of 18; 
• You are pregnant or breastfeeding; 
• You were born prematurely (before 37 weeks, or if twin, before 34 weeks) 

 Do any of these describe you? Yes (I am sorry, you do not qualify to be in our study) or No (continue with questions)

We switched to this "any of the above" style screening script for pregnancy last summer, and it has been working well. We recently reviewed our procedures again, and confirmed that we do not ask questions about sex or gender status or history. But if we did, we'd be looking closely about how exactly the questions were asked and responses recorded, with the aim of collecting the absolute minimum of information required.

bugfix for niftiPlottingFunctions.R, plot.volume()

If you use or adapted the my plot.volume() function from niftiPlottingFunctions.R, please update your code with the version released 30 March 2023. This post illustrates the bug (reversed color scaling for negative values) below the jump; please read it if you've used the function, and contact me with any questions. 

Huge thanks to Tan Nguyen for spotting and reporting the color discrepancy!