COMMIT2

This tutorial illustrates the basics to begin using COMMIT2.

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:

NB: in this tutorial we will make use of the MRtrix toolbox for several key operations to prepare the data for COMMIT. However, one could also use his/her preferred software to achieve the same functionalities.

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( 'tck2connectome -force -nthreads 0 -assignment_radial_search 2 -out_assignments fibers_assignment.txt demo01_fibers.tck GM_1x1x1.nii.gz connectome.csv' )

NB: depending on the specific version of MRtrix installed on your computer, the connectome in the file connectome.csv can be stored using a comma or a space as the delimiter; thus, in the following code please adapt the delimiter parameter accordingly.

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:

if not os.path.isdir( 'bundles' ) :
    os.mkdir( 'bundles' )
os.system( 'connectome2tck -force -nthreads 0 -exclusive -files per_edge -keep_self demo01_fibers.tck fibers_assignment.txt bundles/bundle_' )

C = np.loadtxt( 'connectome.csv', delimiter=',' ) # NB: change 'delimiter' to suits your needs
CMD = 'tckedit -force -nthreads 0'
for i in range(C.shape[0]):
    CMD_i = 'tckedit -force -nthreads 0'
    for j in range(i,C.shape[0]):
        if C[i,j] > 0 :
            CMD_i += ' bundles/bundle_%d-%d.tck' %(i+1,j+1)
    os.system( CMD_i + ' bundles/demo01_fibers_connecting_%d.tck' % (i+1) )
    CMD += ' bundles/demo01_fibers_connecting_%d.tck' %(i+1)
os.system( CMD + ' demo01_fibers_connecting.tck' )

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, verbose=False )
mit.save_results( path_suffix="_COMMIT1" )

Retrieve the streamline contributions estimated by COMMIT for later use in COMMIT2:

x_nnls, _, _ = mit.get_coeffs( get_normalized=False )

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.array( [np.arange(tmp[i],tmp[i+1]) for i in range(len(tmp)-1)] )

Now, we compute the individual weight for each bundle in order to penalize all groups of streamlines in the same manner regardless of their cardinality (see Eq. 6):

group_w = np.empty_like( group_size, dtype=np.float64 )
for k in range(group_size.size) :
    group_w[k] = np.sqrt(group_size[k]) / ( np.linalg.norm(x_nnls[group_idx[k]]) + 1e-12 )

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

reg_lambda = 5e-4 # change to suit your needs

prior_on_bundles = commit.solvers.init_regularisation(
    mit,
    regnorms    = [commit.solvers.group_sparsity, commit.solvers.non_negative, commit.solvers.non_negative],
    structureIC = group_idx,
    weightsIC   = group_w,
    lambdas     = [reg_lambda, 0.0, 0.0]
)

NB: the value for reg_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, regularisation=prior_on_bundles, verbose=False )
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: