This Christmas break I decided to try to get the cover of my new article (Update: looks like they chose another very cool cover). In essence, that work shows that the thalamic reticular nucleus on each side of the brain shows functional connectivity with the visual nuclei on both side of the brain. My goal was to try to express this through a simple and beautiful graph (final image above).
To do this, I needed three types of files:
- Functional MRI data (in my case, the mean run of a participant looking at my stimulus).
- Anatomical masks (these were drawn on the thalamic nuclei using proton density weighted images, but could be any atlas).
- Statistical masks (regions of the brain that survive some statistical threshold, i.e., q[FDR] < 0.05).
The goal is to find, within each anatomical region, x arbitrary time series, and correlate them. To get started, I define some imports and useful functions.
This min-max normalizer is useful when visualizing time series with largely different magnitudes. Since correlations are insensitive to these sorts of scaling operations, it won’t affect our image.
This loads our NIFTI data using Nibabel and reshapes it to a voxels x timepoint 2D matrix, which I find much easier to work with. In this case, our masks will become vectors of integers and we will use them to filter out the columns of the fMRI data that we are interested in looking at.
Here I’m just defining some options for later on: let’s look at the right & left LGN as well as the left ventral TRN. These are stored as integer values
[1, 2, 4] in the anatomical mask. From each of these regions of interest (ROIs), let’s grab 20 arbitrary voxels. Since we are going to save the traces, let’s also make a list of output names here.
Next, we’re going to load in the data, find statistically relevant regions within each target ROI, and plot them as a stack of subplots. This is so we can extract the traces in Inkscape / Illustrator later. Finally, for each ROI, take the 20 traces and stack them onto master
out_ts matrix that we will use to draw our correlation matrix later.
Let’s look at the data we’ve extracted so far:
Next we correlate these time series to create a weighted graph of the data.
I turned off most of the matplotlib label features so I could add in a few of my own in Inkscape. This is what the data looks like:
We are removing negative correlations because there is some debate in the literature about how to best interpret these relationships, and they aren’t particularly relevant to the questions we’re asking.
Finally, we export the matrix using networkx to a format Gephi understands. Gefx is a XLM-based format that is relatively nice to read. We’re also going to number the nodes by ROI (1 = right LGN, 2 = left LGN, 3 = left TRN). This will help us visualize the relationships between these regions in topographic space.
We can import this file into Gephi (version 0.8.2). At first things won’t look very informative at all:
We can use the ‘partition’ feature on the left edge of the overview screen to allow us to color-code our graph according to the labels we assigned earlier:
And finally, we can select from one of the many available layouts to arrange the nodes. These algorithms aren’t paying attention to the labels of the nodes, but rather the distribution of correlations contained within the edges of the graph. This is a ‘topographic’ space, turning anatomical (or ‘real world’) distance into statistical distance can give you amazing insight into your data. I chose the Fruchterman Reingold algorithm to group correlated nodes together in a way that I thought would be pleasing for a magazine cover (note that I also scaled up the edge width to 3 and turned off the ‘curved edges’ option):
The force atlas 1 layout tells a similar story in a clear, more ‘scientific’ visualization:
With a little bit of tweaking, these graphs can be beautiful communication tools.
For a script with the above code as well as the output images along the way, check out the cover section of my repo dedicated to this publication.