Grouping fibers using cortical surface statistics


(Difference between revisions)
Jump to: navigation, search
Rjpatruno (Talk | contribs)
(Created page with "= Grouping fibers using cortical surface statistics = In this project we use visual cortical activations measured by fMRI to define visual areas and then map these areas throu...")

Latest revision as of 14:08, 18 July 2015


[edit] Grouping fibers using cortical surface statistics

In this project we use visual cortical activations measured by fMRI to define visual areas and then map these areas through the white matter of the brain using white matter tracts suggested by a tractography algorithm.

[edit] Background

It has been known for a long time now that most of the posterior of human's brain responds to visual stimuli. The visual cortex includes the whole occipital lobe and also extends into parietal and temporal lobes. More than a century ago, through lesion studies, neurologists realized a correlation between the location of the lesion in the visual cortex and certain deficits and blind spots. By studying this correlation, they defined the primary visual area in the visual cortex (V1), which covers one hemifield in each hemisphere, and responds to visual stimuli in both foveal and peripheral places. Subsequent studies uncovered multiple other maps, for example another area called V2 surrounding V1, etc. Some of these areas are very well-known and everyone agrees on their existance, but the exact size and location of some them are still in discussion. The act of using fMRI measurements of visual cortical activations and our knowledge of the visual maps to define ROIs on the brain is called "Retinotopy". Also, it is known that there is an architecture connecting these areas: we know that signals travel back and forth between this areas. The most famous map of this hierarchy is the one provided by Felleman and Van Essen,1991; Van Essen and Maunsell, 1983 through a monkey study. The only non-invasive way to see if the fibers in the white matter of the brain support this architecture in human brain is through the existing tractography algorithms that use DWI (Diffusion Weighted Imaging) and DTI (Diffusion Tensor Imaging) to tract fibers and fascicles in human's brain. DWI and DTI scans measure how easily water molecules can diffuse in one direction versus another in the white matter of the brain and use the data to estimate the pathways between different points of the brain. The main goal of this project is to first use fMRI data to perform a retinotopy on the brain and define some of the well-known areas on the brain, and then use a group of fibers estimated by a tractography algorithm and associate statistics to them to create a visualization method that would later enable us to study the architecture behind the visual areas of the brain.

[edit] Methods

[edit] Subjects

The subject for the retinotopy part of the project is different from the one used in the second part. I had to change subject for the second part because the DTI data for the first subject was not available. But, the retinotpy on the second subject has been done similarly to the first subject. Both subjects are healthy adult individuals.

[edit] Data acquisition

For performing a retinotopy one needs two kinds of experiments, here are the details of these experiments:

1- Polar Angle Experiment: Blocks of traveling checkerboard wedges Each Block= 8s, Total time of run: 104 s , TR=2s Each Block has 4 stimulus type :up,right,down,left 2 runs / 10 Total Blank blocks

2-Eccentricity Experiment: Blocks of expanding checkerboard rings Each Block= 8s, Total time of run: 104 s , TR=2s Each Block has 4 stimulus type :up,right,down,left 2 runs / 10 Total Blank blocks

The fibers used in the second part of the projects are imported from a tractography algorithm. The fiber group includes 15,203 fibers and they are all either entering or exiting the occipital lobe of the right hemisphere. Here is what the fibers looked like before any processing were done on them:


[edit] Data Analysis

The MR data was motion corrected and the time series of the two runs of each experiment were averaged. The travelling wave analysis was done in MrVista and the ROIs were defined on the inflated surface mesh, using the eccentricity and the polar angle maps and the information and criteria provided in B.Wandell et al. Neuron review paper,2007. For defining LO areas, the MT area was used as a landmark, it was obtained by computing a contrast map on a standard MT localizer scan. The coherence threshold was set to around 0.17. The fibers were imported into QUENCH (a fiber visualizing and editing software written at Stanford by Anthony Sherbondy and Shireesh Agrawal), the eccentricity and polar angle maps along with the ROIs were imported into QUENCH using codes written in MATLAB. Statistics were assigned to the fibers using MATLAB and the fibers were grouped and coloured based on this statistics in QUENCH. The way to do this is as follows:

[edit] Grouping Deterministic Fibers

Step 1 (Exporting maps from MrVista and importing them into QUENCH)-' Export the maps (eccentricity or polar-angle) from MrVista using Volume view Window-> menu-> Xform-> volume to ITKgray-> write itkGray functional data overlay. If you want to export ROIs go to Volume view Window-> menu-> Xform-> volume to ITKgray-> export ROI as NIFTI. This will give you a nifti image of the map or the ROI. In order to be sure that all the exported nifti images are registered to the same anatomical image use the matlab function dtiFixITKGrayHeader.m. This function will ask you for a reference file name to fix the header of your Nifti file based on this reference file. You can use for example the t1.nii.gz file as the reference file. Finally import the anatomical background and the maps or ROIs in quench using menu-> File-> load Volume. You can then change the color maps and the visibility of these overlays using: menu-> panels-> Background and overlays. Be sure to choose the color range of these maps correctly in order to be able to see them.

Step 2 (Fiber tracking)- This step consists of two separate sub-steps, first we have to generate some masks that we use in our tractography process and then we have to track some fibers using dti data. These steps are described in detail below:

  • Step 2.1 (Generating Masks)- In order to track fibers that are exiting occipital lobe in the right hemisphere, you need a white matter exit plane of the occipital lobe that is used to limit the fibers to the occipital lobe. You can easily create this mask by segmenting the b0 image and saving the mask in nifti file. Let's call it wm.nii.gz.
  • Step 2.2 (Generating Fibers)- In order to create fibers to use in the next steps, we have used Camino tractography structure to generate a set of fibers by following a deterministic tracking algorithms. The motivation to use Camino structure is that it provides many tractography options (deterministic or probabilistic; connectivity segmentation tractography etc.) and also it has a comprehensive set of commands that are easy to follow and can be used with a variety of data types.

The goal here, is to create some fibers in the right hemisphere occipital lobe that either enter or exit the occpitial lobe (in other words, pass the occipital lobe exit plane generated above). The way to do this is by following a set of commands described in the following pipeline:

# make Camino output directory Copy all the masks and your raw data to this directory

$ mkdir camino

# Make scheme file: Before we can process the data, we must create a scheme file from bvec- and bval- files that tells Camino about the acquisition sequence. The scheme file is a text file containing data in this format: The diffusion time in seconds The number of gradient directions

gradient 1 x

gradient 1 y

gradient 1 z

gradient 2 x

gradient 2 y


the bscale option enables you to specify the scaling factor to convert the b-values into different units

$ fsl2scheme -bvecfile raw/dti_g150_b2500_aligned_trilin.bvecs -bvalfile raw/dti_g150_b2500_aligned_trilin.bvals -bscale 1E9 -flipx > camino/subject1.scheme2

% Get raw data into camino format: now we should get the raw data into a voxel ordered, camino format file the shredder command extracts periodic chunks of data from a file. Here we have 2-endian short integer datas, so we set the period to be 2 , and the offset to 0.

$ fslchfiletype ANALYZE camino/raw.nii.gz $ cat camino/raw.img | shredder 0 -2 0 | scanner2voxel -voxels 652536 -components 156 -inputdatatype short > camino/raw.Bfloat

% Create tensors: then we fit the diffusion tensor in each voxel using the data file created above and save it in a Bdouble file.

$ dtfit camino/raw.Bfloat camino/subject1.scheme2 > camino/dt.Bdouble

% Checking tensor orientations: Before you begin the whole process of tractography, it is a good idea to check and see if everything is working correctly so far or if any direction needs to be flipped. Camino has visualization command (pdview) that allows you to directly visulize the principle direction of the fitted tensor in each voxel. You have to provide the dimensions of the image.Because pdview needs the diagonalized tensors rather than the tensors themselves, we must pipe the data through an additional command, dteig that diagonalizes the fitted tensors

$ cat camino/dt.Bdouble | dteig | pdview -datadims 81 106 76

Here is for example an image of the coronal view:

Test dt.jpg

% Tracking: Performing tractography requires that we specify one or more seed points, from which tractography is initiated as well as the diffusion tensor file. Here we use the right hemisphere white matter exit plane mask generated in the previous step as the seed point and save the fibers in a Bfloat file.

$track < camino/dt.Bdouble -inputmodel dt -seedfile camino/wm.nii.gz -anisfile camino/wm.nii.gz -anisthresh 0.1 -outputfile camino/tracts.Bfloat

% From all the fibers keep only the streamlines that travel within the occipital lobe: Finally, we process the fibers using the procstreamlines command and the -exclusionfile option. The way it works is that any fiber that is entering the exclusion file is discarded. We use the complement of the white matter mask as the exclusion file. If you want to save the discarded fibers, use the -truncateinexclusion option and save them in a separate file.

cat camino/tracts.Bfloat | procstreamlines -exclusionfile camino/wmNOT.nii.gz -truncateinexclusion > camino/tracts_trunc.Bfloat

% Convert the .Bfloat file to .pdb file to be able to see it in QUENCH The way we do this is by reading the seed file, then getting the xyz and ijk transform matrices from it. Then using them to bring the tracts.Bfloat file into acpc space using dtiXformFiberCoords.m and then exporting the fibers to a pdb file using mrExportFibers (an option will be added to our add statistics function which allows you to input .Bfloat fibers directly into that):

% Load seed volume

seedFile = 'camino/wm.nii.gz';

seedNI = readFileNifti(seedFile);

% Calc xform from camino space to surface data

xformMMToSeed = seedNI.qto_ijk;

xformMMToSeed(1:3,4) = 0;

xformSeedToAcPc = seedNI.qto_xyz;

xformAcPcToSurface = surfaceNI.qto_ijk;

xformMMToSurface = xformAcPcToSurface * xformSeedToAcPc * xformMMToSeed;

% Load camino path file

tractsFile = 'camino/tracts.Bfloat';

fg = mtrImportFibers(tractsFile,eye(4));

% Get fibers into AcPc

fg_acpc = dtiXformFiberCoords(fg, xformSeedToAcPc * xformMMToSeed);

% write the fibers into a pdb file


Step 3 (Adding statistics to the fibers)- Use the matlab function dtiAddCorticalSurfaceStatsToFibers.m to add statistics to your fibers (the .pdb file). This function uses the NIFTI files exported from MrVista and adds statistics to your fibers and overwrites the .pdb file with the fibers containing statistics. Read the help document of the matlab function dtiAddCorticalSurfaceStatsToFibers.m for more information on how to do this. (When you have multiple ROIs and want to have a single NIFTI file to group the fibers based on it, the easy way to do this is to export them one by one from MrVista, then create a new NIFTI file with the same fields as one of the ROIs, and as its data you can add the data of all the ROIs with different weights and then group the fibers based on this new NIFTI file containing all the ROIs. For example: =*2+*3+...)

Step 4 (Fiber grouping)- Import the fibers containing the statistics in Quench using menu->load pathways. Then, select all the fibers, go to menu-> panels-> refine selection and choose the first group based on the statistics. Then you can choose another color, select all the fibers, go to the same menu option and choose the second group, etc. You can also choose a colormap for your fiber groups so that they would match your overlay using menu-> view-> choose colormap for fiber group.

Step 5 (Saving the state)- Be sure to save the state in QUENCH after grouping the fibers, so that you can access them at a later time. menu-> File-> save state.

[edit] Generating Probabilistic Fibers

Now if you follow these steps and look at the results in QUENCH (shown below) you may notice that the fibers don't look that good and there are a lot of empty spots in the white matter. So, in order to fix this, we tried creating better looking fibers using a probabilistic tractography algorithm again in Camino. Also this time, we tried using better masks in our tractography algorithm. By better masks, we mean masks generated from T1 anatomical image(which has a much better resolution than the dti image) instead of a b0 image. Se here are the steps we followed:

Step 1 (Segment and label the high resolution T1 image)- This step is performed in Itkgray. You can segment and label the white matter and the gray matter in the right hemisphere and save them in a nifti file, let's call it t1_class.nii.gz.

Step 2 (take the segmentation nifti image to the same xyz space as the b0 image)- If you look at the qto_xyz and qto_ijk transform matrices of t1 and b0 images, you would see that they are different, in other words, each of them are transferred to the xyz space using a different transform matrix. Since we use the dti data in our tractography process, we want our masks to be in the identical xyz space as dti data, otherwise we would have misaligned data. So, what we did here was writing a matlab function called dtiCreateHighResImageForTractography.m, which takes the two images (for example t1 and b0) and re-slices the t1 image to have the same xyz transform matrix as our dti image using the mrAnatResliceSpm.m function. It is crucial to run this function on the masks you create from t1 segmentation image, before you can use them in your tractography process. After running this function, if you want to check that your new file is in the same xyz space as b0, you can use mricron which is a image visualization platform. If you can see the two images on top of each other, then they are in the same xyz space. Note that for the reslice function, we use a spline kernel that differs from its default value ([0 0 0 7 7 7]). We use [0 0 0 0 0 0] so that it does a 0th order or nearest neighbor interpolation. Also you HAVE to specify your bounding box as we have done it in our function.

Step 3 (Generate masks)- Now we are ready to generate some masks. We have written a Matlab script called s_generate_masks.m which has the pipeline for doing this. What this script does is that it creates a gray matter mask in the occipital lobe and a white matter exit plane mask from the segmentation file. Let's call them gm.nii.gz and wm_exit.nii.gz.

Step 4 (Track fibers using a probabilistic tractography algorithm- Again you can use Camino structure to create fibers using a probabilistic tractography algorithm. The pipeline to do this is as follows:

% Do all the deterministic tractography steps up to fitting diffusion tensors to data

% Estimate SNR in b=0 image: the SNR is required to evaluate the uncertainty in principal directions.

$ estimatesnr -inputfile camino/raw.Bfloat -schemefile camino/subject1.scheme2 -bgmask camino/gm.nii.gz

% Define uncertainty lookup table: Before defining PDFs first we define the lookup table that maps tensor shape to the estimated uncertainty for that tensor. We do this for all the tensors in the scheme file.

$ dtlutgen -schemefile camino/subject1.scheme -snr 16.0 -inversion 1 > camino/picoTable.dat

% Create the PDFs: This step actually takes a long time.

$ picopdfs -inputmodel dt -luts camino/picoTable.dat > camino/brainDT_pico.Bdouble

% Track fibers: Now we are ready to track fibers.

$ track -inputmodel pico -inputfile camino/brainDT_pico.Bdouble -seedfile camino/wm_exit.nii.gz -iterations 1000 -outputtracts voxels > picoVoxels.Bshort

% Keep only the fibers in the occipital lobe: We do this by using the exclusionfile command and using gm mask as our exclusionfile but this time we want only the discarded fibers.

$ cat camino/picoVoxels.Bshort | procstreamlines -exclusionfile camino/wm_exit.nii.gz -truncateinexclusion > camino/tracts.Bshort

[edit] Connectivity Segmentation

Now if you look at these fibers, they should look much better than the previous set of fibers. When you create these new fibers, another cool thing you can do is to perform a connectivity segmentation on the occipital lobe exit plane using Camino. A connectivity segmentation allows you to generate an image of this plane using an atlas created from a retinotopic map, that for each voxel shows to which part of the map it is connected with the highest probability. For example you could generate an at atlas from your retinotopic ROIs and then do this connectivity segmentation to see which part of the exit plane is mostly connected to V1, which part to V2, etc. The way to do this is as follows:

Step 1 (Generating the atlas)- By atlas, we mean the gray matter mask that has the retinotopic map information in it. For example for the ROI case, you can use the nifti file that contains all your ROIs, and then your occipital lobe gray matter mask (all created above) and the label your gray matter mask and create an atlas. Note than because the ROIs were originally defined in t1 space, to use them in your tractography process, you have to first convert your ROIs file to b0 space. Now save this atlas as a nifti file, say atlas.nii.gz.

Step 2 (Connectivity Segmentation)- We do this by using the streamline processing command of camino. The seed file shows where the fibers start, the target file shows where they are connected to (in our case the atlas) and the waypoint file constrains streamlines to pass through certain volumes within the image.

$ procstreamlines -inputfile camino/tracts.Bshort -inputmodel voxels -outputcp -iterations 1000 -seedfile camino/wm_exit.nii.gz -targetfile camino/atlas.nii.gz -waypointfile camino/gm.nii.gz -outputroot camino/cc_cbs_ -outputcbs

This command generates an image of the exit plane mask with the required data.

[edit] Results

[edit] Retinotopy Results

Here are the ROIs on both hemispheres resulted from performing retinotopy on the subject:

ROI RH.jpg

ROI LH.jpg

Here is the mean time series, the fft of the mean time series, and the coverage map for the V1 area in the right hemisphere that can be used as a verification of the results.


[edit] Fiber grouping Results

You can make different visualizations in QUENCH with the available MATLAB codes, following the steps described above. Here are some examples: Here is a series of screenshots of the ROIs exported from MrVista, that can provide a good 3D visualization of the position of the ROIs in the visual cortex:

ROI movie map.jpg ROIS.jpg

Here is a screenshot of the eccentricity map imported to QUENCH from MrVista:


Here is a screenshot of the fibers grouped based on their eccentricity values:

                                       Group eccs.jpg

Here is a series of screenshots showing the endpoints of the fibers grouped and coloured based on their eccentricity values as we travel from back of the occipital lobe to the front in the right hemisphere:


Here is a screenshot of the fibers passing V1 and V2 areas of the brain:

V1 fiber.jpg

V2 fibre.jpg

[edit] Summary

This project provides a mean for exporting different retinotopic maps as well as ROIs from MrVista as well as a platform for visualizing them on a 3D anatomical image. You can also import fibers resulted from a tractography algorithm and assign statistics to them and group them based on these maps and show them with the maps, all on one anatomical image. This way we can create some interesting visualizations as the ones provided above. You can see some interesting results in some of the visualizations: For example: In the screenshots showing the endpoints of the fibers grouped and coloured based on their eccentricity values you can see how a cluster of eccentricity mapping fibers are forming and heading out of the occipital lobe. Or in the following screenshot showing only the foveal and peripheral mapping fibers, you can see that the down going stream tends to be more foveal, and the up going stream seems to be more peripheral. When one considers the theory of ventral and dorsal streams in the brain, which says that the ventral stream carries information to places responsible for object recognition and the dorsal stream carries information to places responsible for spatial orientation recognition, this can be an interesting result.


Also, we have introduced some of the processing and visualizations that one can do using Camino structure to help studying the mapping of retinotopic maps through the white matter.

[edit] References

[1] B. A. Wandell, S.O. Dumoulin and A. A. Brewer , "Visual Field Maps in Human Cortex", Neuron, V. 56 , p. 366-383, 2007

Personal tools