BOLD movement scrubbing#

Why scrub BOLD data based on movement parameters?#

When a subject moves in the scanner during BOLD acquisition, a number of things happen that result in possible movement artifacts-induced changes in signal intensity that are not related to actual changes in BOLD signal due to underlying neural activity, but rather due to movement itself. Some consequences of movement can be corrected using movement correction—individual volumes within the time series are realigned so that each voxel/grayordinate represents the same part of the brain. Regressing movement parameters out of the signal can also reduce the extent of the artifactual, movement-driven change in the signal. Large movements, however, still lead to relatively large artifactual changes in the signal.

When running activation analyses, a large part of the artifactual signal will most probably be averaged out and the final estimate of signal magnitude due to task or event of interest won't be too far off. Movement artifacts will however lead to increased error estimate and reduced power to detect activation within a single session, and might also increase variability of response magnitudes and reduce power in second level analyses.

When running functional connectivity, any noise signal that is shared between voxels or regions of interest will induce positive or negative correlation. The larger the ratio between correlated signal changes due to movement and due to signal of interest, the larger will be the artifactual correlations observed in the final results.

In order to increase the power of activation analyses and reduce the artifactual correlations in functional connectivity analyses, one has to minimize the effect of subject movement on BOLD signal. The best way of achieving that is to have subjects be as still as possible during scanning. Not all is lost, however, if there are just a few large movements. In this case, it is possible to identify the volumes (frames) in the BOLD image that might be affected by movement artifacts and ignore them either in activation analyses or when computing functional connectivity.

What follows is a description of how to go about identifying, visualizing and ignoring problematic volumes using QuNex tools.

The strategy#

The strategy in movement scrubbing is the following:

  1. Compute statistics that reflect movement and its artifactual properties

  2. Identify problematic frames

  3. Ignore or interpolate problematic flagged frames in analyses

The tools described use two sources of information that can be used to identify bad frames. The first one are movement correction parameters. Movement correction is a key part of basic preprocessing of BOLD data. Any package used for movement correction provides an option to report movement correction parameters that describe the amount of shift and rotation that was needed to align each volume of the BOLD data to the reference volume. Tracking the changes in these (most often 6) parameters provides an estimate of frame-to-frame movement of the subject. Having that information, it is possible to combine it into one parameter reflecting the overall amount of movement and use it to identify the frames where movement exceeded an arbitrary threshold. The parameter computed and used by the QuNex tools is frame displacement as described in Power et al. (2012). The default threshold value is 0.5.

The second source of information is the acquired BOLD timeseries itself. By comparing intensity values from frame to frame, one can track the signal changes across the BOLD timeseries and identify frames in which signal change exceeds a chosen threshold. The specific parameter used to measure signal change by QuNex tools is dvars, described by Power et al. (2012). In short, dvars is computed by calculating difference between two timeframes for each voxel within the brain separately. The differences are then squared and a square root of the mean is computed. To normalize it, the computed RMS (root mean square) is divided by mean intensity and multiplied by 100 to reflect mean percent change of the signal. Power et al. (2012) suggest to use 3 as a sensible threshold value for the resulting parameter.

Dvars reflects not only movement effects on the signal but also the effects of signal of interest and overall noise, and quite possibly other properties such as voxel size, TR, brain size, tissue composition and others. For that reason some BOLD timeseries have higher mean dvars then others, reflecting higher baseline frame-to-frame differences. When baseline dvars is higher, relatively smaller changes in the signal can lead to crossing the chosen threshold. To account for the differences in baseline frame-to-frame variability QuNex tools also express dvars at each frame as a ratio to baseline dvars, computed as median across the BOLD timeseries. The resulting dvarsme statistic for each frame is interpreted as the factor of change compared to baseline (e.g. dvarsme of 2 means that the change in that frame was 2 times higher than baseline change between successive frames across the run). The proposed default dvarsme threshold is 1.5.

Using both sources frames with excessive movement and/or signal change are easily identified. As the changes can be smeared in time due to filtering and other manipulations of the BOLD timeseries Power et al. (2012) suggest to reject not only the identified frame, but also one preceding and two following frames. QuNex tools allow specification of how many frames before and after identified artifactual frames to be also marked as bad.

To choose the most appropriate threshold and have a feel for your data, it is prudent to visually check the movement parameters and the extent of the suggested frames to be removed. QuNex tools provide means for both visual inspection as well as tabular review of the data. Do make use of them instead of just accepting the defaults.

Once the frames are identified the software for running the relevant analyses needs to be informed about which frames to ignore. The tools provide two venues for specifying which frames to remove or ignore in analyses. The first option is to record the relevant frames in standard .fidl file format. This option creates a .fidl file in which each line records the beginning of contiguous frames to be ignored and the number of frames to be ignored as a negative value. This is the format that fidl as well as QuNex scripts that work with event-related data will accept. To use the information specified in such generated .fidl files, they first need to be combined appropriately with existing .fidl files. This is also accomplished with the use of QuNex tools and existing .conc and .fidl files.

The second option in specifying the frames to be ignored is the .scrub file that compute_bold_stats writes beside each BOLD run, as bold[N].scrub. This is a regular plaintext table in which each frame is marked either with 0 (do not ignore) or 1 (ignore). There is a separate column for each of the discussed criteria (mov for frame displacement, dvars and dvarsme) as well as for their intersection (idvars, idvarsme — a frame is considered bad if both criteria are met) and their union (udvars, udvarsme — a frame is considered bad if either of the criteria is reached). A final use column marks the frames to keep under the criterion the run was scrubbed with. These column names are exactly the values that --mov_bad, --mov_post and --mov_fidl take.

In the following sections tools for performing each of the listed steps will be presented in more detail. First, the two methods used to compute the statistics and scrubbing will be presented, followed by the functions and commands that use them to process image files in bulk.

Supported methods implemented in MATLAB / Octave#

For computation of BOLD statistics, identification of bad frame img_stats_time and img_compute_scrub methods are used, respectively. Note that these are 'core' functions implemented in MATLAB (or Octave). These functions are called automatically by 'higher-order' QuNex BOLD de-noising tools. However, a user can envoke them specifically within MATLAB.

img_stats_time method#

img_stats_time is a method that computes a desired list of statistics across all specified voxels for each frame of the BOLD image. To use the method run:

[stats, done] = img.img_stats_time(do, mask);

Both parameters are optional. Do should be a cell array of strings that specify what statistics to compute and mask the object that defines a mask for voxel selection. Statistics are computed only across voxels that have a non-zero value in the mask image. If do is absent or empty, all the statistics will be computed by default. If no mask is provided, statistics are computed across all voxels in a volume. The statistics to be computed are specified with the following strings:

  • n — number of non-NaN values in the mask for each frame

  • m — mean across all voxels in the mask for each frame

  • me — median across all voxels in the mask for each frame

  • max — maximum value within all voxels in the mask for each frame

  • min — minimum value within all voxels in the mask for each frame

  • sum — a sum of values of voxels in the mask for each frame

  • sd — standard deviation computed across all voxels in the mask for each frame

  • var — variance computed across all voxels in the mask for each frame

  • dvars — signal change across all voxels in the mask between each frame and the preceding frame

The results are stored in a structure stats that consists of one-dimensional vectors in named fields. The names of the fields are the same as the commands to compute them. The only exception is computation of dvars, which computes two additional statistics: dvarsm is computed by dividing the value of dvars for each frame by the mean for that frame and multiplying by 100; dvarsme is computed by dividing the value of dvarsm for each frame by the median dvarsm computed across all frames. For explanation see the text above. A call:

stats = img.img_stats_time({'m','var','dvars'}, mask);

will return the following fields:

   stats.m
   stats.var
   stats.dvars
   stats.dvarsm
   stats.dvarsme

Besides the structure with the statistics, the method also returns a cell array done that lists the executed commands.

img_compute_scrub method#

img_compute_scrub is a method that evaluates BOLD motion and image statistics to identify bad frames. It is used by calling:

img = img_compute_scrub(img, comm);

img is the image to be evaluated, and comm is the command description of how to compute scrubbing. The command is provided as a string of 'key:value' pairs separated by a pipe (|) character. The following keys can be specified, with the defaults in []:

  • radius ... head radius in mm [50]

  • fdt ... frame displacement threshold [0.5]

  • dvarsmt ... image intensity normalized root mean squared error (RMSE) (dvarsm) threshold [3.0]

  • dvarsmet ... median normalized RMSE (dvarsme) threshold [1.6]

  • after ... how many frames after the identified bad ones are to be marked as well [0]

  • before ... how many frames before the identified bad ones are to be marked as well [0]

  • reject ... which criteria is to be used for identification of bad frames [udvarsme]

The possible criteria for identification of bad frames are:

  • mov ... frame displacement threshold (fdt) is exceeded

  • dvars ... image intensity normalized root mean squared error (RMSE) threshold (dvarsmt) is exceeded

  • dvarsme ... median normalized RMSE (dvarsmet) threshold is exceeded

  • idvars ... both fdt and dvarsmt are exceeded (i for intersection)

  • udvars ... either fdt or dvarsmt are exceeded (u for union)

  • idvarsme ... both fdt and dvarsmet are exceeded

  • udvarsme ... either fdt or udvarsmet are exceeded

The method is used on an mrimage img object that has both movement parameters data (img.mov field structure) and per frame statistics data (img.fstats structure) already computed and present. Movement data is read automatically when the image is loaded if the data is present in the relevant folder (images/functional/movement). Frame statistics data is also read if present. It is computed using img_stats_time method and is added and saved when using general_compute_bold_stats command (see below).

Note, the image is not scrubbed, only the scrubbing parameters are computed and the chosen criterium is used to create the img.use mask that in other functions and methods defines what frames to use or not.

Processing BOLD images#

QuNex tools provide multiple ways of processing BOLD for image statistics and scrubbing. One possibility is to use the methods described above directly. The other method is to use QuNex MATLAB and QuNex python functions to process single BOLD image files (general_compute_bold_stats QuNex MATLAB command), a list of BOLD image files (general_compute_bold_list_stats QuNex MATLAB command) or BOLD image files across sessions (qunex compute_bold_stats terminal command).

general_compute_bold_stats function#

When you want to perform complete processing of a BOLD image file, i.e. reading in the relevant data, computing statistics and scrubbing, and saving the results, general_compute_bold_stats is the command to use. The function is run as:

general_compute_bold_stats(img, mask, target, store, scrub, verbose);

with the following parameters:

  • img ... A path to a BOLD file to process.

  • mask ... A path to a mask file to use.

  • target ... A folder to save results into.

    If target is set to '' (the default), the results are saved where the BOLD image is, if target is specified as 'none', results are not saved to an external file.

  • store ... Whether to store the resulting data in the image file.

    If store is set to 'same', the computed statistics and scrubbing information are embedded and saved in the original BOLD image itself. If store is set to <ext>, the information is stored embedded in a copy of the original BOLD image, with extension <ext>.

    If store is left blank (i.e. '', the default), the results are not embedded and saved in an image file.

    NOTE: Due to changes of the default format, this option is currently disabled and might be removed.

  • scrub ... A string describing whether and how to compute scrubbing information.

    See img_compute_scrub above. If 'none' is specified, no scrubbing computation is to be performed.

  • verbose ... Whether to report the progress or not (the default).

The command initially identifies relevant brain voxels in two manners. First, it identifies voxels with intensity higher than 300 on the first BOLD frame. If there are more than 20000 valid voxels, it then selects those for which the intensity is above the specified threshold throughout the timeseries and selects those for computation of image statistics. Second, if the first method fails (e.g. in the case when images were demeaned), it identifies all the voxels for which the variance across the frames is more than 0. After the voxels are identified, the image is additionally masked if a mask was specified.

To compute the per frame statistics on the identified and masked brain voxels, the command runs img_stats_time method and saves the computed statistics in a simple text file named the same as the BOLD file it relates to with .bstats extension. The structure of the results is rather straightforward. The file starts with a header line that lists all the variables, whereas all the following lines list statistics, one frame per line. All the values are separated by tabs. The columns are provided in the following order: frame (frame number), n, m, min, max, var, sd, dvars, dvarsm, dvarsme.

If scrubbing was specified it also runs img_compute_scrub method and creates a .scrub file with detailed scrubbing information, as well as .use file with information on which frames are to be used, and which not. Both are simple text files. .scrub file starts with the header line that specifies the columns: frame, mov, dvars, dvarsme, idvars, idvarsme, udvars, udvarsme, use. The following lines provide the information on frame number, whether the frame is marked as bad (1) or not (0) by each specified criterium, and whether the frame is to be used (1) or not (0) based on the criterium specified when running the function.

general_compute_bold_list_stats function#

When you want to compute BOLD statistics for a larger number of images using MATLAB, general_compute_bold_list_stats is the right command to use. To use the command execute:

general_compute_bold_list_stats(flist, target, store, scrub, verbose);

with the following parameters:

  • flist ... .list text file providing a list of session's BOLD image or conc files to process.

    See list files for more information.

  • target ... A folder (relative to where BOLD image is stored) to save results into.

    If target is set to '' (the default), the results are saved where the BOLD image is,

    if target is specified as 'none', results are not saved to an external file.

  • store ... Whether to store the resulting data in the image file.

    If store is set to 'same', the computed statistics and scrubbing information are embedded and saved in the original BOLD image itself. If store is set to <ext>, the information is stored embedded in a copy of the original BOLD image, with extension <ext>.

    If store is left blank (i.e. '', the default), the results are not embedded and saved in an image file.

    NOTE: Due to changes of the default format, this option is currently disabled and might be removed.

  • scrub ... A string describing whether and how to compute scrubbing information.

    See img_compute_scrub above. If 'none' is specified, no scrubbing computation is to be performed.

  • verbose ... Whether to report the progress or not (the default).

The command runs general_compute_bold_stats on each BOLD file specified in the provided .list file. The only parameter that has to be provided is filelist, whereas target and verbose are optional. If they are not provided, the general_compute_bold_stats defaults will be used. If filelist specifies an ROI file for a session (see list file format, that ROI will be used as a mask for computing statistics. If an ROI file is not specified, only the general_compute_bold_stats internal method for identifying brain signal will be used.

QuNex compute_bold_stats command#

Note that the compute_bold_stats is a 'higher-level QuNex command that wraps the 'core' functions described above. Thus, the easiest way to compute BOLD image statistics for all the sessions in a study is to use the compute_bold_stats QuNex command. This command will process the specified BOLD image files for all the sessions in the specified batch file, taking into account parameters for filtering or explicitly specifying which sessions and BOLD files to work with. Further, all the specific scrubbing parameters described above can be specified. For detailed information about the use of the command, and the results, please consult inline help by running qunex compute_bold_stats --help. In brief, this is a description of key parameters:

general parameters

When running the command, the following general processing parameters are taken into account:

--batchfile      ... The batch.txt file with all the sessions information.
--sessionsfolder ... The path to the study/sessions folder, where the imaging data is supposed to
                     go [.].
--parsessions    ... How many sessions to run in parallel [1].
--parelements    ... How many elements (e.g. bolds) to run in parallel [1].
--overwrite      ... Whether to overwrite existing data (yes) or not (no). Note that previous data
                     is deleted before the run, so in the case of a failed command run, previous
                     results are lost [no].
--bolds          ... Which BOLD images (as they are specified in the batch.txt file) to work on.
                     It can be a single type (e.g. 'task'), a pipe separated list (e.g.
                     'WM|Control|rest') or 'all' to work on all [rest].
--boldname       ... The default name of the bold files in the images folder [bold].

specific parameters

In addition the following specific parameters define the actual results:

--mov_radius  ... Estimated head radius (in mm) for computing frame displacement statistics [50].
--mov_fd      ... Frame displacement threshold (in mm) to use for identifying bad frames [0.5]
--mov_dvars   ... The (mean normalized) dvars threshold to use for identifying bad frames [3.0].
--mov_dvarsme ... The (median normalized) dvarsm threshold to use for identifying bad frames [1.5].
--mov_after   ... How many frames after each frame identified as bad to also exclude from further
                  processing and analysis [0].
--mov_before  ... How many frames before each frame identified as bad to also exclude from further
                  processing and analysis [0].
--mov_bad     ... Which criteria to use for identification of bad frames [udvarsme].

This command is best used within the general preprocessing workflow of BOLD data for purposes of functional connectivity preprocessing or preparation for task activation analysis.

Examples#

Using the defaults:

qunex compute_bold_stats \
    --batchfile="<path to batch file>" \
    --sessionsfolder="<path to study sessions folder>" \
    --overwrite="<whether to overwrite existing files (yes) or not (no)>" \
    --bolds="<which BOLD files to process, e.g. 'all' or 'rest'>"

Specifying additional parameters for identification of bad frames:

qunex compute_bold_stats \
    --batchfile="<path to batch file>" \
    --sessionsfolder="<path to study sessions folder>" \
    --overwrite="<whether to overwrite existing files (yes) or not (no)>" \
    --bolds="<which BOLD files to process, e.g. 'all' or 'rest'>" \
    --mov_fd="0.9" \
    --mov_dvarsme="1.6" \
    --mov_before="1" \
    --mov_after="2"

The breakdown of the parameters is as follows:

  • mov_fd="0.9" ... sets frame displacement threshold to 0.9

  • mov_dvarsme="1.6" ... sets dvarsme threshold to 1.6

  • mov_before="1" ... exclude also one frame before identified bad ones

  • mov_after="2" ... exclude also two frames after identified bad ones

QuNex create_stats_report#

Once compute_bold_stats has written the statistics and the scrubbing flags, create_stats_report turns them into something you can read: per-session plots and fidl snippets, and group level reports that let you review movement and scrubbing across the whole study at once.

It computes no scrubbing of its own. Which frames are bad was decided by img_compute_scrub and written into the .scrub file, and this command reads that file back — so --mov_post and --mov_fidl name a column of it, and nothing else. The work is done in Python, in process, and the figures are drawn with matplotlib.

For each session, into images/functional/movement:

  • bold<nifti_tail>_<mov_plot>_cor.pdf — the six movement correction parameters, one panel per BOLD run.

  • bold<nifti_tail>_<mov_plot>_dvars.pdf — frame displacement and dvars, with the frames identified as bad marked.

  • bold<nifti_tail>_<mov_plot>_dvarsme.pdf — the same for frame displacement and dvarsme.

  • bold[N]<nifti_tail>_scrub.fidl — a fidl snippet listing the frames to exclude from analysis.

The plots are also copied into sessions/QC/movement/<mov_pdf>/, in a cor, dvars and dvarsme subfolder each, so that the whole study's figures can be flipped through in one place.

For the study as a whole, into sessions/QC/movement:

  • <mov_mreport> — for every session and BOLD run, eight rows describing the six movement parameters. Three describe the movement itself — mean, sd and span (maximum minus minimum) — and four describe the frame-to-frame change: max, md (mean absolute), med (median absolute) and md2_max. The eighth row, frame_dspl, describes framewise displacement and carries four values rather than six: its mean, median, maximum and standard deviation.

  • <mov_preport> — the same statistics with the bad frames excluded, so you can see what scrubbing bought you.

  • <mov_sreport> — for every session and BOLD run, how many frames and what percentage each of the seven criteria rejects. This is the table to consult when choosing which criterion to scrub by, because it shows all of them side by side.

Session statistics are appended to the group files as they are computed, so run the command with --parsessions=1 to keep those files in a sensible order.

For detailed information, run qunex create_stats_report --help. In brief, this is a list of key parameters.

general parameters

When running the command, the following general processing parameters are taken into account:

--batchfile      ... The batch.txt file with all the session information.
--sessionsfolder ... The path to the study/sessions folder, where the imaging data is supposed to
                     go [.].
--parsessions    ... How many sessions to run in parallel. Keep at 1 so that the group reports are
                     appended in order [1].
--overwrite      ... Whether to overwrite existing data (yes) or not (no) [no].
--bolds          ... Which BOLD images (as they are specified in the batch.txt file) to report on.
                     It can be a single type (e.g. 'task'), a pipe separated list (e.g.
                     'WM|Control|rest') or 'all' to work on all [rest].
--boldname       ... The default name of the bold files in the images folder [bold].
--nifti_tail     ... The tail of NIfTI volume images to use, which selects the .bstats and .scrub
                     files the report is built from [].
--bold_variant   ... Optional variant of bold preprocessing, if the BOLD images live in
                     images/functional<bold_variant> [].

scrubbing specific options

In addition the following scrubbing specific parameters define computation of results:

--mov_radius  ... Estimated head radius (in mm) for computing frame displacement statistics [50].
--mov_fd      ... Frame displacement threshold (in mm) to use for identifying bad frames [0.5]
--mov_dvars   ... The (mean normalized) dvars threshold to use for identifying bad frames [3.0].
--mov_dvarsme ... The (median normalized) dvarsm threshold to use for identifying bad frames [1.5].
--mov_after   ... How many frames after each frame identified as bad to also exclude from further
                  processing and analysis [0].
--mov_before  ... How many frames before each frame identified as bad to also exclude from further
                  processing and analysis [0].

reporting specific options

Finally, the following options determine how the results will be reported and saved:

--tr          ... TR of the BOLD files, used when writing the fidl snippets [2.5].
--mov_pref    ... The prefix to be used for the figure plot files [].
--mov_plot    ... The base name of the plot files. If set to empty no plots are generated [mov_report].
--mov_mreport ... The name of the group movement report file. If set to an empty string, no file is
                  generated [movement_report.txt].
--mov_sreport ... The name of the group scrubbing report file. If set to an empty string, no file is
                  generated [movement_scrubbing_report.txt].
--mov_preport ... The name of group report file with stats computed with frames identified as bad
                  exluded from analysis. If set to an empty string, no file is generated
                  [movement_report_post.txt].
--mov_post    ... The criterium for identification of bad frames that is used when generating a post
                  scrubbing statistics group report (mov/dvars/dvarsme/idvars/idvarsme/udvars/
                  udvarsme/none) [udvarsme].
--mov_fidl    ... Whether to create fidl file snippets with listed bad frames, and what criterium to
                  use for the definition of bad frames (mov/dvars/dvarsme/idvars/idvarsme/udvars/
                  udvarsme). Set to none to not generate them [udvarsme].
--mov_pdf     ... The name of the folder in sessions/QC/movement in which to copy the individuals'
                  movement plots [movement_plots].

Again this command is best used within the general preprocessing workflow of BOLD data for purposes of functional connectivity preprocessing or preparation for task activation analysis.

Example#

qunex create_stats_report \
    --batchfile="<path to batch file>" \
    --sessionsfolder="<path to study sessions folder>" \
    --overwrite="<whether to overwrite existing reports and files (yes) or no (no)>" \
    --bolds="<which BOLD files to process, e.g. 'all' or 'rest'>" \
    --parsessions="1"

Reviewing scrubbing across the whole study#

The commands above report on one session at a time. To see the scrubbing of every session and run in the study together — how many frames each run lost, how that compares with the rest of your sample, and which runs are worth a closer look — use run_qc_summary. It reads the same .bstats and .scrub files described here, so the frames it reports are the frames your preprocessing excluded, and it plots the per-frame traces in an interactive report.

Additional supporting functions#

QuNex also provides additional supporting functions.

meltmovfidl#

To be able to use the generated _scrub.fidl files they possibly need to be combined and melted with existing .fidl files. Combination ensures that when working with multiple runs, as specified by .conc file, _scrub.fidl that are generated for each run are combined in the correct order and their timing information is adjusted to correspond to the right frame within the concatenated run. Melting enables scrubbing information to be combined with existing event information in a fidl file. Melting is accomplished using:

qunex meltmovfidl <concfile> <ignore_fidl_pattern> <input_fidl_file> <output_fidl_file>

All the listed parameters have to be provided. --concfile provides a path to the .conc file that describes the sequence of BOLD runs to be concatenated in analyses. --ignore_fidl_pattern specifies the ending for the fidl files to be melted with the original fidl file. Typically the ending would be scrub.fidl. --input_fidl_file lists the path to the original fidl file that the scrubbing information is to be added to. --output_fidl_file gives the path to the file in which the melted data is to be saved.

In the following example:

qunex /data/study1/s01/wm.conc _scrub.fidl /data/study1/s01/wm.fidl /data/study1/s01/wm_scrubbed.fidl

The command will generate a new fidl file, based on /data/study1/s01/wm.fidl and /data/study1/s01/wm.conc, that will include scrubbing information in *_scrub.fidl files and save it to /data/study1/s01/wm_scrubbed.fidl file.

Warning on file names and paths

It needs to be taken into account that to correctly match BOLD images listed in the conc file with the corresponding per-run fidl files, both have to include in their name a sequence b*N where N corresponds to a BOLD number. Examples would be s01_bold2.4dfp.img, ad7325_b04_x3rd_atl.nii, id285_brun3_scrub.fidl. In addition, per-run fidl files to be used are expected to reside in the movement folder that corresponds to the session from which the BOLD run is taken from.