Structural pipeline¶

Diffusion weighted MRI processing of a subject starts with the main function structural_pipeline that executes processing steps with reconstruction parameters specified on the command line or in a configuration file in JSON-format. Parameters in the user provided configuration file overwrite the default parameters and parameters specified on the command line overwrite both. The Structural configuration page gives an overview of all parameters in the structural configuration file.

Structural preprocessing¶

The preprocessing step takes FreeSurfer and DWI data and generates intermediate files for the consecutive processing steps. The structural preprocessing step executes preprocessing specified in a bash script (specified by the parameter preprocessingScript).

Three example preprocessing scripts (preprocess_topup_eddy.sh, preprocess_eddy.sh and preprocess_minimal.sh) are provided with the structural pipeline.

Script preprocess_topup_eddy.sh uses

FSL Topup and FSL eddy to correct for susceptibility induced distortions, eddy current distortions and motion artifacts in the DWI data [1],
(using indexFile and acqpFile, and generating dwiProcessedFile)

updates the b-vectors to adjust for the DWI corrections [2],
(updated b-vectors and b-values are saved in processedBvalsFile and processedBvecsFile)

computes a DWI reference image based on the corrected diffusion-unweighted (b0) volumes,
(generating dwiReferenceFile)

computes the registration matrix between DWI reference image and the T1 image using bbregister [3] and
(generating registrationMatrixFile)

registers the Freesurfer segmentation to the DWI reference image.
(generating segmentationFile)

The script preprocess_eddy.sh performs the same steps except correction for susceptibility induced distortions using FSL topup. The minimal-preprocessing script preprocess_minimal.sh does not perform any DWI data corrections (see Table below).

Show advanced notes

The example scripts use by default eddy_openmp executable. A different version of eddy can be specified in the eddyVersion parameter.

Sometimes the b-values file (rawBvalsFile or processedBvalsFile) contains nonzero b-values for scans that are intended as diffusion unweighted b0-scans. However, diffusion reconstruction methods (often) assume b-values equal to 0 for diffusion unweighted scans. The pipeline accomodates this by putting all b-values that are smaller or equal to bValueZeroThreshold to 0.

The pipeline checks that the norm of all b-vectors is very close to one (i.e. all b-vectors are unit vectors). If the norm of any b-vector deviates from 1 more than bValueScalingTol, then a warning is given on the command line. No b-value scaling or b-vector scaling is performed.

Custom preprocessing script¶

The user can also specify a custom preprocessing script using the preprocessingScript configuration parameter. The custom preprocessing script will get all general and structural_preprocessing configuration parameters as input (in the form: --parameterName=parameterValue). For example, when using the configuration parameters

{
    "structural_preprocessing":{
        "dwiFile": "DTI/DTI.nii.gz",
        "customParameter": "test",
        "preprocessingScript": "CONFIGDIR/my_custom_script.sh"
    }
}

then the structural_preprocessing script executes the script my_custom_script.sh that is located in the same directory as the configuration file with customParameter as one of the input arguments:

/some/directory/my_custom_scipt.sh \
        --dwiFile="DTI/DTI.nii.gz" \
        --customParameter="test" \
        ...other parameters...

Tip: Make sure that the custom preprocessing script is executable.

Overview of preprocessing steps in provided preprocessing scripts.¶
	Minimal	Eddy	Topup Eddy	Custom
Correct susceptibility induced distortions (using FSL topup)	no	no	yes	?
Correct eddy current distortions and motion artifacts in the DWI data (using FSL eddy) [1]	no	yes	yes	?
Update vectors to adjust for the DWI corrections [2]	no	yes	yes	?
Compute a DWI reference image based on the diffusion unweighted (b0) volumes (using FSL)	yes	yes	yes	?
Compute registration matrix between DWI reference image and the T1 image (using Freesurfer) [3]	yes	yes	yes	?
Register the Freesurfer segmentation to the DWI reference image (using Freesurfer).	yes	yes	yes	?

Anatomical steps¶

Parcellation is an optional step that calls the FreeSurfer software suite to create additional cortical parcellations of the surface with respect to reference atlases. Standard atlases in CATO include:

‘Desikan-Killiany’ atlas present in FreeSurfer [4],
120, 250 and 500 regions Cammoun sub-parcellations of the Desikan-Killiany atlas [5] and
‘Von Economo-Koskinas’ cortical region and cortical-type atlas [6, 7]

CATO has a template directory (TOOLBOXDIR/templates). In this template directory, each atlas has a directory (named TEMPLATE) with a parcellation script (parcellate_TEMPLATE.sh), a ROIs file (ROIs_TEMPLATE.txt) and template specific files. All filenames are defined as variables (templatesDir, templateScript, ROIsFile) in the configuration file and can be adjusted.

Templates provided with the pipeline use FreeSurfer’s mris_ca_label to create the annotation files, mri_aparc2aseg to merge the cortical atlases with the sub-cortical ASEG parcellation and mris_anatomical_stats to provide for every template an anatomical stats file describing a.o. for each region:

the number of vertices
surface area (mm²)
gray matter volume (mm³)
average thickness (mm)

The modular character of CATO allows easy implementation of other surface based or volume based atlases by adding a template to the template directory and adding this new template to the general templates parameter that specifies which templates are used in the pipeline.

Show advanced notes

Standard both functional and structural pipelines include the parcellation step. To avoid the re-computation of FreeSurfer atlas files and stat files and (very small) differences in regional properties between both pipelines, CATO checks whether FreeSurfer atlas files already exist. If a file already exists, then CATO does not overwrite this file. To force CATO to overwrite FreeSurfer files, set parameter forceFreesurferOverwrite to true.

The collect_region_properties step collects volumetric and surface data from the regions of interest defined in the ROIsfile and summarizes this data in the regionPropertiesFile. The following statistics are included for every brain region:

center of mass of each region (calculated from the parcellation file parcellationFile)
the number of vertices (from FreeSurfer output created in the parcellation step)
surface area (mm²) (from FreeSurfer)
gray matter volume (mm³) (from FreeSurfer)
average thickness (mm) (from FreeSurfer)

Reconstruction diffusion¶

The reconstruction_diffusion step estimates the white matter fiber organization in each voxel from the measured DWI data. Structural connectivity modeling is based on the principle that white matter fibers restrict the movement of water molecules resulting in peaks (preferred directions) in the diffusion patterns of water molecules. Pipeline step reconstruction_diffusion provides three methods to infer diffusion peaks from DWI data, including:

Diffusion Tensor Imaging (DTI)
Constrained Spherical Deconvolution (CSD)
Generalized Q-sampling Imaging (GQI)

The applied methods are specified by the reconstructionMethods parameter:

"reconstructionMethods": ["DTI", "CSD", "GQI", "GQI_DTI", "CSD_DTI"]

The DTI option models the measured signal of a voxel by a tensor and estimates one preferred diffusion-direction per voxel. CATO uses the informed RESTORE algorithm [8, 9] to reduce the impact of physiological noise artifacts on the DTI modeling, performing tensor tensor estimation while identifying and removing outliers during the fitting procedure. The Levenberg-Marquardt method as implemented by Gavin et al. [10] is used to solve the nonlinear least squares problem.

Show advanced notes

The iRESTORE method identifies measurements that are outliers and that are excluded from the tensor fitting. To ensure that enough information is preserved for reliable tensor estimation the iRESTORE algorithm checks that the B-matrix is well conditioned and directionally balanced. If this is not the case then all measurements are used to perform the nonlinear least-squares fitting. The B-matrix is considered well-conditioned if the condition number is lower than thresCondNum and directionally balanced if the variation in average projection scores is lower than thresVarProjScores. The condition number and variation in average projection scores thresholds are specific for each gradient acquisition scheme and are therefore estimated by the function thresholdAssistant. This function estimates both thresholds using a bootstrapped sample of gradient schemes obtained by randomly removing gradient directions from the original scheme. Following experience, as described in de Reus (2015), thresholds where chosen from the distribution such that:
The removal of 5% of the gradients was allowed for 75% of the samples.

And the removal of 50% of the gradients was allowed in 25% of the samples.
For each threshold these two points were obtained across all permutations and the final thresholds were obtained by averaging the two points. 2. The non-linear least squares fitting procedure used parameters suggested by Gavin et al. [10] (\(\epsilon_1\) = 0.001, \(\epsilon_2\) = 0.001, \(\epsilon_3\) = 0.1, \(\epsilon_4\) = 0.1, \(\lambda_{\uparrow}\) = 11, \(\lambda_{\downarrow}\) = 9, \(\lambda_0\) = 0.01 and a maximum of 100 iterations).

Constrained spherical deconvolution reconstructs the diffusion peaks in a voxel by deconvolution of the measured signal with the diffusion profile associated with a fiber [11]. Signal deconvolution is performed in the super-resolved spherical harmonics framework to allow for a natural description of the surface of a diffusion process. The used implementation of spherical deconvolution is constrained to eliminate spherical harmonics with negative values which reduces high frequency noise resulting in the reconstruction of a few well-defined peaks [11]. The implementation in the pipeline follows the implementation as used in the Dipy software package and as described in [11] and was extended to allow super-resolved reconstructions.

Generalized Q-sampling Imaging extrapolates the diffusion signal after which diffusion directions with the highest signal are selected as diffusion peaks. The implementation follows the method described in Yeh et al. [12] and the implementation as presented on http://dsi-studio.labsolver.org.

Comparing the four reconstruction methods, DTI is a robust and simple method, whereas CSD and GQI are more advanced diffusion models that allow for the reconstruction of fibers in multiple directions, which is beneficial for voxels with a more complex white matter organization. To combine the strengths of simple and advanced reconstruction methods, CATO allows the user to combine DTI with CSD or GQI to perform DTI modeling in voxels with one peak and CSD or GQI in voxels with multiple detected peaks.

The Reconstruction diffusion step provides the user with the additional option to export diffusion measures to a NIFTI volume file diffusionMeasuresFileNifti. To export fractional anisotropy, axial diffusivity, radial diffusivity and mean diffusivity measures to a file, add the following parameters to the reconstruction_diffusion section in the configuration file:

{
        "exportNifti":{
                "exportNifti": true,
                "measures": ["fractional anisotropy", "axial diffusivity", "radial diffusivity", "mean diffusivity"],
                "diffusionMeasuresFileNifti": "OUTPUTDIR/SUBJECT_MEASURE.nii.gz"
        }
}

Reconstruction fibers¶

The reconstruction_fibers function performs fiber tracking based on the diffusion peaks of each voxel. The fiber tracking step utilizes an extended version of the “Fiber Assignment by Continuous Tracking” (FACT) algorithm [13], which is a deterministic tracking algorithm that starts fiber reconstruction from seeds in the white matter and propagates streamlines in the main diffusion axis of the voxel while updating the propagation direction each time the tip of the streamline enters a new voxel. Fiber reconstruction starts from one or multiple seeds (number of seeds set by NumberOfSeedsPerVoxel) in all voxels with a matching segmentation label as provided by startRegions.

Fiber reconstruction stops if a tracker is

in a region with fractional anisotropy lower than minFA.
in a stopping region with segmentation code included in the stopRegions variable
about to revisit the current voxel or previously visited voxel,
about to enter a forbidden region with segmentation code included in the forbiddenRegions variable
or about to make a sharp turn, i.e. an angle > maxAngleDeg.

Show advanced notes

Both the Reconstruction fibers and the Network reconstruction step include parameters maxAngleDeg and minFA, but both parameters are specific for their reconstruction step.

Reconstruction fiber properties¶

The reconstruction_fiber_properties step identifies fiber segments that connect brain regions and calculates fiber measures as preparation for the network reconstruction step. For each atlas and reconstruction method, the reconstruction_fiber_properties step iterates through all fibers and if a fiber crosses two or more brain regions of interest, as defined by the regions of interest file (ROIsFile), then the shortest fiber segment between each region pair is included in the fiber properties file (fiberPropertiesFile). In addition to the start and end point of each fiber segment and the associated region pair, additional measures for fiber segments are stored, including:

maximum turn angle (in radials),

minimum fractional anisotropy

fiber length (physical length in mm),

average fractional anisotropy (FA),

axial diffusivity (AD),

radial diffusivity (RD),

mean diffusivity (MD) and

generalized fractional anisotropy (GFA)

Diffusion measures (FA, AD, RD, MD and GFA) are averaged over voxels weighted by the length of the traversed path through each voxel [14].

Network reconstruction¶

The reconstruction_network step builds the connectivity matrices for each of the selected (sub)cortical atlases and reconstruction methods. Brain regions included in the connectivity matrix, and their order, are defined by the regions of interest file (ROIsFile). The network matrix is constructed by iterating through all fiber segments in the fiberPropertiesFile and fiber segments connecting regions of interest are added to the connectivity matrix. Fibers can additionally be filtered by their projection length (with only fibers longer than minLengthMM being included in the network reconstruction), minimum fractional anisotropy (including only fibers that touch voxels with fractional anisotropy higher than minFA) and maximum angle (including only fibers that make turns in their trajectories smaller than maxAngleDeg).

Connections are assigned weights by

the number of streamlines (i.e. the number of fibers) that connect two regions (NOS)

fiber length (physical length in mm)

fractional anisotropy (FA)

axial diffusivity (AD)

radial diffusivity (RD)

mean diffusivity (MD)

generalized fractional anisotropy (GFA)

streamline volume density (SVD, number of streamlines divided by the average volume of both regions) and

streamline surface density (SSD, number of streamlines divided by the average surface area of both regions)

Diffusion measures (FA, AD, RD, MD and GFA) are averaged over voxels weighted by the length of the traversed path through each voxel [14].