The connection between Linus Pauling and fMRI

When I think of Linus Pauling, two things come to mind: his work on  the nature of the chemical bond (for which he was awarded the 1954 Nobel Prize in Chemistry) and, in his later years, his advocacy of Vitamin C. But Pauling is one of the great names of  20th century science. He contributed to a variety or domains and other connections can be made. His campaigning for a nuclear test ban won him a Nobel Peace prize. And his description of the molecular structure of blood -- of sickle cell anemia in particular -- is highly regarded.

 Pauling, in fact, had had a long interest in the properties of blood. In 1936 he published an important paper on the magnetic properties of hemoglobin [1]. Hemoglobin, which is a major component of blood, can be oxygenated or deoxygenated and it is due to this easy coupling and decoupling that it plays its key role as a carrier of oxygen. Pauling's experiments had shown that deoxyhemoglobin is paramagnetic (i.e.,  in layman's terms, has magnetic properties) and oxyhemoglobin is diamagnetic. The original paper is worth reading; I was struck by the clear, concise and objective manner in which these results are  presented.

Pauling did not envision an application at the time but from the mid-1960s on much work was done on engineering MR field gradients and studying the contrast they provided for different kinds of tissues. In 1990, Seiji Ogawa published a set of papers [2] showing that local inhomogeneities in the magnetic field were formed in response to  changing levels of deoxygenated blood and could be measured. This property became the basis of BOLD or blood oxygenation level dependent fMRI.

The connection  between neural activity and the MR signal  is roughly as follows: 

        increased local brain activity
               → increased blood flow 
                   → change (dilation) in size of blood vessel
                          → change (increase) in the deoxyhemoglobin concentration
                                → changes in the local magnetic field as measured by the MRI (T2*) signal. 

There is an excellent lecture on YouTube where Geoffrey Aguirre explains this. Another talk, this time by Christof Koch, also briefly mentions Pauling's paper.


References  
(1) The Magnetic Properties and Structure of Hemoglobin, Oxyhemoglobin and Carbonmonoxyhemoglobin.
Linus Pauling, Charles D. Coryell. Proc Natl Academy of Sciences. 1936 April; 22(4): 210–216.

(2) Brain magnetic resonance imaging with contrast dependent on blood oxygenation.
Seiji Ogawa et al.  Proceedings of the National Academy of Sciences. 1990; 87(24): 9868-9872.

Processing a diffusion phantom with FSL

This is a how-to style post intended for FSL users. It may also be of interest to a broader diffusion imaging audience.

In diffusion imaging, a carefully constructed diffusion phantom may be used as a ground truth reference to evaluate competing diffusion models and tractography algorithms. The phantom constructed for the 2009 MICCAI Fiber Cup challenge [1] is one such reference; it is a coronal slice model of fibers that cross, touch and split into different configurations and with different curvatures. The phantom dataset is now publicly available.

Fiber pathways in the diffusion phantom
Image credit: Fibercup paper [1]

  I recently processed the raw phantom dataset to evaluate some algorithms. These are some of my notes on how to preprocess the data on FSL through bedpostx. bedpostx is an FSL program that computes the distribution of diffusion parameters at each voxel. It generates all the files necessary to run the probabilistic tractography tool probtrackx.

• The raw data for the phantom can be obtained from the Fiber Cup download page. There are 6 datasets (2 image resolutions each acquired at 3 different b-values). The data is stored in 4D NIfTI format and has already been preprocessed to correct for eddy currents and movements so one can run bedpostx directly.
• In order to run bedpostx, you need the bvals and bvecs gradient files, as well as the input data file and the nodif_brain_mask.
  • data.nii is the 4D phantom image file.
  • bvals is an ASCII text file of b-values corresponding to the gradient directions; the b-values are listed as a single row.
  • bvecs are the gradient vectors; the file format has 3 rows which correspond to the x-, y- and z-values of the gradient directions.
  • nodif_brain_mask.nii.gz is a 3D binary mask file; the file is generated using the FSL bet (or bet2) command.

Sample commands for the dwi-1500.nii dataset (3x3x3 mm resolution, b=1500):

1. Download the 3x3x3 dwi-b1500.nii file from the Fibercup download page. Rename the file data.nii.
2. Download the diffusion_directions_corrected.txt file. This is a file of coordinates with 130 entries. The first 65 lines (baseline + 64 directions) are repeated (65 x 2 = 130) since two passes were taken. This was most probably done to average over noise.
   This file, which I have renamed bvecs.tmp, needs to be formatted into the bvecs file so that the x-, y- and z-columns are transposed as rows (see 5 below for sample code).
3. To generate the bvals file note that there should be 130 values since the number of entries should correspond to the number of gradients. The first value is 0 (for b0); the next 64 values should be 1500 (for b=1500); the final 65 values are a repeat of the first 65. This is the unix command-line perl code I used to create this file:
   perl -e 'print "0\t", "1500\t" x 64, "0\t", "1500\t" x 64' > bvals
4. The nodif_brain_mask file is generated using the following BET (brain extraction tool) command:
   bet2 data.nii nodif_brain_mask.nii -m -n -f 0.15
   The -m flag specifies a binary mask; the -f flag is the mask threshold.
5. Run dtifit, an FSL command that fits a diffusion tensor model at each voxel.
   dtifit -k data.nii -m nodif_brain_mask.nii.gz -o out -r bvecs -b bvals
   This extra step, suggested by my colleague Julio Duarte, can be used to check the directions in the bvecs file. This is done by viewing the output file out_V1.nii.gz on FSLview and verifying that the primary fiber direction at each voxel aligns with the phantom's.
   
   
   
   
   
   
   Here we see that the horizontal x-direction is aligned with the phantom but the vertical y-direction needs to be replaced by -y. 
   
   
   
   
   
   
   
   
   To correct this I used the following command-line perl code to regenerate the bvecs file:
   perl -ane 'print "$F[0] "' bvecs.tmp > bvecs
   perl -e 'print "\n"' >> bvecs
   perl -ane 'print $F[1]*-1," " ' bvecs.tmp >> bvecs
   perl -e 'print "\n"' >> bvecs
   perl -ane 'print "$F[2] "' bvecs.tmp >> bvecs
   perl -e 'print "\n"' >> bvecs
   
   If we re-run dtifit, we will see that the primary vectors align with the phantom's fiber pathways.
   
   
   
6. Run bedpostx using the directory where the data, bvals, bvecs and nodif_brain_mask files are located:
   bedpostx ./ 
7. As an optional final check compare the bedpostx vector output (dyads*.nii) to the dtifit output (out_V1.nii). 

The bedpostx dyads1 output (blue arrows) is overlaid on the dtifit out_V1 (red lines) vectors. We can see that they align with each other.









Reference
[1] Quantitative evaluation of 10 tractography algorithms on a realistic diffusion MR phantom. Fillard P, Descoteaux M, Goh A, Gouttard S, Jeurissen B, Malcolm J, Ramirez-Manzanares A, Reisert M, Sakaie K, Tensaouti F, Yo T, Mangin JF, Poupon C. Neuroimage 2011; 56(1):220-34.

Some resources
Takuya Hayashi has a program to swap the bvec values.  (Note: I have not verified that it works.)

The see-through brain

A group at Stanford led by Karl Deisseroth have pioneered a chemical process that leaves brain (and other) tissue intact but optically transparent. The procedure, aptly named CLARITY, exposes internal structure by replacing the opaque lipid layers that surround cells with a transparent shape-preserving hydro-gel. Neurons, glia, genes and other molecules of interest can then be selectively highlighted in this transparent brain by molecular markers. The results are spectacular and reminiscent of swimming among the jellyfish.

This is a link to Nature which published the article. The accompanying set of video clips is compelling; in one frame, a single nerve projection is traced out among a tangle of other neurons. The technique will provide a boost for the field of connectomics; researchers who would otherwise have to reconstruct neural wiring from small piece-meal sections of brain will be able to view whole-brain connections in one glance.

Update August 10 2013:   A short clip of Karl Deisseroth on CLARITY.

The use of a spatial distribution model in labeling sulci

While accurate sulcal identification can be a challenge even for expert neuroanatomists, there are sulci that are to some degree more consistent, and for which anatomical correspondence can be established across subjects. These are the larger primary sulci. The localization of these sulci allows us to generate a spatial distribution or probabilistic map which can be used to label candidate sulci. A graph that maps the structural relationships between sulci can also be constructed and unlabeled sulci (or the more variable secondary and tertiary sulci) can be identified against this reference.

These two ideas, the use of the probabilistic atlas and the graph, have been incorporated into automated and semi-automated labeling methods in various ways. In this post I will present the basic idea behind the use of the spatial distribution model.

The use of a probabilistic atlas
Probabilistic maps compute the probability for each tissue class at every voxel location using a large database of segmented and labeled anatomical structures. Evans et al. [3] coined the name Statistical Probabilistic Anatomical Maps or SPAM for these models. Paul Thompson has a nice description of these SPAM models and the Brainvisa website has a nice visualization of a sulcal atlas which is reproduced below:

A straightforward implementation of the probabilistic atlas paradigm can be seen in Le Goualher et al. [1] [2]. SPAM models give the probability for each sulcal class so that at any given location, unlabeled sulci are assigned the most probable label for that location. In other words, to label a new sulcus :

Let:           be a sulcal label

Compute:     where p is the probability from a SPAM atlas
Assign:      


The use of a point distribution model
A different spatial distribution model is used by Lohmann et al. [4]. A point distribution model introduced by Cootes et al. [5] computes the shape of sulcal basins across a training set. Any unlabeled sulcus can be expressed as a linear combination of the eigenvalues generated from the PCA of this shape covariance matrix; an optimization over the linear function gives the best label.

Spatial distribution models give spatial bounds but this is not adequate to discriminate between the sulci in a local region. They are usually combined with graphs which model connections between sulci thus giving local structural context. In the combined strategy, the spatial information is used to supply spatial priors [6], localization constraints or to narrow the search space in an optimization or graph matching process [7].

I will write about the use of graph models in my next post.

References
1) Georges Le Goualher, D. Louis Collins and Christian Barillot, Alan C. Evans, "Automatic Identification of Cortical Sulci Using a 3D Probabilistic Atlas," In MICCAI, 1998, pp. 509-518.
2) Georges Le Goualher, E. Procyk, D.L. Collins, R. Venugopal, Christian Barillot, "Automated Extraction and Variability Analysis of Sulcal Neuroanatomy," IEEE Trans. Med. Imag., 18(3), 1999, pp. 206-217.
3) A. C. Evans, D. L. Collins, P. Neelin, M. Kamber, S. Marrett, "Three-dimensional correlative imaging: Applications in human brain mapping," Advances in Functional NeuroImaging: Technical Foundations,(ed. R. Thatcher and M. Hallett and T. Zeffiro and E. John and M. Huerta) Academic Press, 1994, pp. 145-162.
4) Gabrielle Lohmann and Y. von Cramon, "Automatic labeling of the human cortical surface using sulcal basins," IEEE Trans. Med. Imag., 4, 2000, pp. 179-188.

5) Timothy F. Cootes, Christopher J. Taylor, David H. Cooper, Jim Graham, "Active Shape Models-Their Training and Application," Computer Vision and Image Understanding, 61(1), 1995, pp. 38-59.
6) M. Perrot, D. Rivière, J.-F. Mangin, "Identifying cortical sulci from localizations, shape and local organization," ISBI, 2008, pp. 420-423.
7) Yang, F & Kruggel, F., "A graph matching approach for labeling brain sulci using location, orientation, and shape," Neurocomputing, 2009, pp. 179-190.

Posts on Sulcal Labeling
1) Why we label sulci
2) Why is sulcal labeling difficult ?
3) The use of a spatial distribution model in labeling sulci

Why is sulcal labeling difficult?

This is a follow-up to an earlier post Why we label sulci. There will be two or three more posts; taken altogether, they will describe the sulcal labeling problem.

The labeling of sulci is a challenging problem. This is because, cortical sulci are highly variable. Sulci vary not just across individuals but even between the hemispheres of a single brain [1]. It might be useful when looking for ways to address this variability to classify this variation as follows:

Variation in physical features
Sulci vary in shape, in scale and in their placement (i.e. position and orientation) The figure below illustrates how the variability can make feature selection difficult.

The boxplot shows the length distribution for 18 subjects. The 10 types or classes of sulci shown cannot be identified solely on a length measurement. This poses a problem for feature selection and classification.
Figure credit: Meena Mani




Variation in branching
19th century illustrations such as those from Horsley [2], trace the wide variations along a sulcal fold. A whole nomenclature has developed since then to account for the branch variations possible along a single sulcus. (An example from the Ono atlas is illustrative--see figure below). For this reason, there is no gold standard in sulcal labeling; one neuroanatomist may disagree with another.

The figure to the left shows the pattern variations for a single sulcus (the posterior end of the superior frontal sulcus). Types A, B, C, D, are possible variations for this sulcus (for the 25 postmortem brains examined, 4 variations were found). The pattern in the two hemispheres of a single subject may differ; the left may be Type B and the right may be Type C. The lengths of the small segments and the connections they make to other sulci may also vary. Reproduced from Ono et al. [1].

Variation in number
Sulci may be continuous (present as one uninterrupted segment) in some individuals, fragmented (exist as multiple segments) in others and altogether absent in yet others. The larger primary sulci which start forming early in fetal development are the most consistent; the secondary and tertiary sulci are not always expressed.

References
1) Ono, M., Kubic, S. & Abernathy, C. (1892) "Atlas of the Cerebral Sulci", (Thieme, New York).
2) Horsley, V. (1892) "On the topographical relations of the cranium and surface of the cerebrum", In "Contribution to the surface anatomy of the cerebral hemispheres", pp.306-355, (Royal Irish Academy).


Posts on Sulcal Labeling
1) Why we label sulci
2) Why is sulcal labeling difficult ?
3) The use of a spatial distribution model in labeling sulci

Medical visualization frontier application

The goal of medical imaging is to present the data in a useful format and in a very interesting TEDx talk Anders Ynnerman demonstrates cool new applications in medical visualization that will be possible in the near future.

Graphics processors have become substantially faster in the last ten years. We can now put together** the gigabytes (terabytes if extended to the time domain) of MRI and CT data generated when scanning a single subject and create 3D (or 4D) images from which relevant information can be selectively extracted. This opens up new and very interesting possibilities. One such application is the virtual autopsy where with ipad-style interactions one can look at cadavers in hard-to-maneuver angles or selectively view metal to, for instance, identify the extent of knife stab injuries or locate bullet shards. Ynnerman also suggests touch-sensitive haptic applications: a surgeon can literally touch the data--a beating heart for example--pre-surgery.

It's really a great 17 minute talk--here's the link.

**(Aside from fast GPUs, there are other ways in which people are hoping to handle the explosion of data from these medical scans. The use of oompressive sensing algorithms is one but the more general idea is to reduce the data before, during or after the scan.)

Shape analysis: on selecting a set of points or landmarks

One approach to shape analysis (after Kendall) uses a fixed number of points to define a shape. The points may describe an object boundary or an interior morphology such as the veins on a leaf. They may be selected randomly; alternatively they may be landmarks--i.e. points of significance. A set of points, so selected, constitutes a shape summary and the original shape data, that had been extracted from the raw image, is discarded.

One problem with representing a shape in this way is that it introduces a source of variability. Different shapes can be reconstructed with the same set of points. This is especially true if the number of points selected for that particular shape are few.

There are, of course, other approaches to shape analysis that do not involve selecting points (at least not at this stage of shape representation). Deformable templates is one such methodology and its use is quite common in medical image analysis. The shape could also be represented by a continuous function (see for example Younes et al. [1]). But these methods are also computationally expensive.

Reference
1) Younes, Laurent; Michor, Peter W.; Shah, Jayant; Mumford, David. A metric on shape space with explicit geodesics. Atti Accad. Naz. Lincei Cl. Sci. Fis. Mat. Natur. Rend. Lincei (9) Mat. Appl. 19 (2008), no. 1, 25--57.