COMMIT2

This tutorial illustrates the basics to begin using COMMIT2.

Important

Commands in this tutorial refer to the latest version of COMMIT. Older versions may not be compatible. The list of changes and updates can be found here.

Download data

This tutorial uses the same dataset as the Getting Started tutorial:

• MRI data
• tractogram

Then, download into the same folder the following file:

• GM_1x1x1.nii.gz: parcellation of the gray matter for defining the bundles of interest.

After downloading all the data, the content of the folder should look like this:

Initial setup

Open the Python interpreter, go to the folder containing all the files and import the following required libraries:

import numpy as np
import os
import commit
from commit import trk2dictionary
# commit.setup() # NB: this is required only the first time you run COMMIT

Group the streamlines in bundles

First, we calculate the connectome corresponding to the input tractogram, i.e. demo01_fibers.tck:

os.system( 'dice_connectome_build fibers_assignment.txt connectome.csv -t demo01_fibers.tck -a GM_1x1x1.nii.gz -f')

Using this information, we can extract from the tractogram only the connecting streamlines, i.e. those that do connect two valid gray-matter ROIs as defined in the file GM_1x1x1.nii.gz and re-organize them as contiguous groups in demo01_fibers_connecting.tck:

os.system( 'dice_tractogram_sort demo01_fibers.tck GM_1x1x1.nii.gz demo01_fibers_connecting.tck -f')

Import the tractogram

Import the tractogram just created into COMMIT:

trk2dictionary.run(
    filename_tractogram = 'demo01_fibers_connecting.tck',
    filename_mask       = 'WM.nii.gz',
    fiber_shift         = 0.5
)

Initial evaluation with standard COMMIT

COMMIT2 requires an initial guess for the relative importance of each bundle, to consider each bundle equally regardless of its cardinality (see Eq. 6 in the corresponding publication here). This is required to not penalize small bundles just because they are small or favor big ones because they are big. We can get such an initial guess by running the standard COMMIT:

# convert the bvals/bvecs pair to a single scheme file
import amico
amico.util.fsl2scheme( 'bvals.txt', 'bvecs.txt', 'DWI.scheme' )

# load the data
mit = commit.Evaluation( '.', '.' )
mit.load_data( 'DWI.nii.gz', 'DWI.scheme' )

# use a forward-model with 1 Stick for the streamlines and 2 Balls for all the rest
mit.set_model( 'StickZeppelinBall' )
d_par       = 1.7E-3             # Parallel diffusivity [mm^2/s]
d_perps_zep = []                 # Perpendicular diffusivity(s) [mm^2/s]
d_isos      = [ 1.7E-3, 3.0E-3 ] # Isotropic diffusivity(s) [mm^2/s]
mit.model.set( d_par, d_perps_zep, d_isos )

mit.generate_kernels( regenerate=True )
mit.load_kernels()

# create the sparse data structures to handle the matrix A
mit.load_dictionary( 'COMMIT' )
mit.set_threads()
mit.build_operator()

# perform the fit
mit.fit( tol_fun=1e-3, max_iter=1000 )
mit.save_results( path_suffix="_COMMIT1" )

Preparing the anatomical prior on bundles

First, we retrieve the number of streamlines present in each bundle from the connectome previously calculated, i.e. connectome.csv:

C = np.loadtxt( 'connectome.csv', delimiter=',' )
C = np.triu( C ) # be sure to get only the upper-triangular part of the matrix
group_size = C[C>0].astype(np.int32)

With this information, we create an array with the indices of the streamlines composing each bundle:

tmp = np.insert( np.cumsum(group_size), 0, 0 )
group_idx = np.fromiter( [np.arange(tmp[i],tmp[i+1]) for i in range(len(tmp)-1)], dtype=np.object_ )

Now, we initialize a dictionary of additional parameters for the group lasso regularisation. In this way the individual weight for each bundle is computed internally to avoid the introduction of a penalization bias due to different bundles' size (see Eq. 6).

params_IC = {}
params_IC['group_idx'] = group_idx
params_IC['group_weights_cardinality'] = True
params_IC['group_weights_adaptive'] = True

Finally, we create the regularization term implementing the anatomical prior that axons are organized in bundles. Please recall that this requires a percentage of the maximum value of the regularization parameter lambda to be defined, i.e. perc_lambda:

perc_lambda = 0.00025 # change to suit your needs

# set the regularisation
mit.set_regularisation(
    regularisers   = ('group_lasso', None, None),
    is_nonnegative = (True, True, True), 
    lambdas        = (perc_lambda, None, None),
    params         = (params_IC, None, None)
)

Important

The value for perc_lambda used in this example is only an empirical value for the purpose of this demo and should not be taken as a reference or ground truth value. This value should be adapted to your data!

Evaluation with the new COMMIT2

We can now fit the streamlines to the diffusion MRI signal using COMMIT2 (i.e. COMMIT with bundle regularization):

mit.fit( tol_fun=1e-3, max_iter=1000 )
mit.save_results( path_suffix="_COMMIT2" )

Comparison between COMMIT and COMMIT2

The following figure shows the (binary) connectomes corresponding to the raw input tractogram (left) as well as after the filtering with COMMIT (center) and COMMIT2 (right):

As we can see, the connectomes are sparser with COMMIT and COMMIT2, as indeed the network density (computed with the brain connectivity toolbox) is 0.66, 0.54, and 0.35, respectively.

However, despite sparser, the streamlines that were kept by COMMIT/COMMIT2 exhibit a similar distribution (i.e., streamline density) and explain equally well the reference diffusion MR signal (e.g., signal RMSE), as shown in the following figure: