Elastic Registration of Single Subject Task Based fMRI Signals

Task based functional magnetic resonance imaging (fMRI) is a powerful approach for examining time varying changes in brain metabolism as a response to a predetermined stimulus [3]. Two main challenges for fMRI analysis are (i) estimation of the shape or pattern of the response (ii) detection of the blood oxygenation level dependent (BOLD) activity, as it only represents a small percentage of the variance of the signal. Non-neuronal contributions such as physiological noise, thermal noise and hardware instabilities, hamper small changes in brain activity from being detected. Low statistical power is therefore mostly caused due to a low sample size or due to a low effect size of the conditions being studied. In order to increase statistical power, research studies typically increase the sample size and perform a group-level analysis to improve both detection and estimation. In this work we focus on the problem of improving statistical power of activations for single subject fMRI analysis by performing functional alignment.

Recently, there have been several interesting approaches [4, 5, 9, 10] that have proposed the synchronization or alignment of fMRI signals primarily for resting state data. Although the end goal in all these methodologies is aligning fMRI time courses to each other, in our work we adopt a different approach. Particularly in the task-based fMRI setting, we would expect the alignment to achieve increased consensus between parts of the time course that correspond to the same stimulus condition. There are several alternatives for solving this problem. (i) Given a set of time series to be aligned, one could choose any arbitrary time course randomly from this set as a template and align all functions to it. This can introduce a potential bias towards the arbitrary chosen template. (ii) On the other hand, one can achieve all possible registrations of functions and choose a particular time series template that satisfies a certain optimality condition. This issue is similar to the one faced in image/shape registration, where the choice of the template or atlas is critical. Further, the process of alignment itself can be achieved either using a least squares criterion [4, 5] or using the dynamic time warping (DTW) approach [10]. The process of least squares alignment is usually non-elastic and thus doesn’t account for temporal reparameterizations. The DTW approach does provide an elegant way of recovering temporal misalignment, but does not directly allow for the computation of an intrinsic average, which is important for template estimation on the space of time series functions. Thus one still needs to resort to extrinsic (generally Euclidean) minimization criteria to determine the optimal template.

In our work, we provide a single solution to the above two problems. We use an elastic time series matching approach to construct an unbiased template of this shape pattern and align the condition-specific time course activity to this template. We treat the problem of aligning fMRI data as a shape matching problem [9] and to achieve elastic alignment between time courses to bring the BOLD responses into correspondence. On a single subject level, our method significantly detects more clusters and more activated voxels in relevant subcortical regions in healthy controls. To our knowledge, this is the first work that performs nonlinear elastic temporal registration of task-based fMRI signals.

Schematic of the elastic hyperalignment approach. Top row: the original time course color-coded by the fMRI conditions that included happy, fearful, and neutral. Row 2, 3: the unaligned and the aligned time courses for those conditions. Row 4: the elastically aligned and the unaligned mean respectively. Row 5: within-condition aligned time course.

Full size image

This section outlines the main approach. Figure 1 shows a schematic of the procedure. The notation and discussion below is described for a fMRI time series acquisition using a blocked design at a single voxel in the brain. Let the conditions for the task and the number of blocks per condition be denoted by \(C_i, i = 1, \ldots , k\) and \(B_j, j = 1, \ldots , n\) respectively. Generally we assume that the total number of blocks stay fixed across all conditions. Each block is assumed to have a boxcar shape that stays zero throughout the resting period and has a unit magnitude for the duration of the task. Thus for a given task condition \(C_i\), the function for block \(B_j\) is denoted by \(g^{C_i}_{B_j} (t)= U(t- t^i_{j1}) - U(t- t^i_{j2}), \forall t \in [0, T]\), where T is the total acquisition time. Here \(t^i_{j1}\) and \(t^i_{j2}\) are the stimulus on and off times for condition \(C_i\) and block \(B_j\). Thus the blocked design stimulus function for the entire duration of the task is given by \(g^C_B(t)= \sum _{i=1}^k \sum _{j=1}^n g^{C_i}_{B_j} (t), \forall t \in [0, T]\). Let the corresponding acquired fMRI time course due to the BOLD hemodynamic response be given by \(h (t), \forall t \in [0, T]\). Our goal is to achieve within-condition temporal reparameterization or alignment across the entire fMRI time series given by h. Therefore for each condition \(C_i\), we extract the BOLD response functions as the set \(\{ f^i_j (t) \}\), where

We now seek an optimal temporal alignment or registration among the set of functions given by \(\{ f^i_j (t) \}\). This is achieved using the square root velocity parameterization [6, 7, 12] for the representation of functions. We describe this idea in brief below and refer the reader to [6, 7, 12] for more details. For the discussion below, we drop the subscripts i, j for convenience. For a given fMRI time course signal \(f:I \equiv [0, 1] \rightarrow \mathbb {R}\), and its velocity \(\dot{f}(t) = \frac{df}{dt}\) and magnitude \(|\dot{f}(t)|\), we define its functional representation by the square-root velocity field (SRVF) map q given by \(q:[0,1] \rightarrow \mathbb {R}, \ q(t) = \frac{\dot{f}(t)}{\sqrt{|\dot{f}(t)|}}\). For an absolutely continuous f, the SRVF transformation ensures that q is square integrable. The set of SRVFs is then given by \(\mathbb {L}^2([0,1],\mathbb {R})\), which is a Hilbert space. To recover the time domain fMRI signal, we use \(f(t) = f(0) + \int _0^t q( \tau ) |q(\tau )| d \tau \). The SRVF mapping is invertible up to a given f(0). We assume \(f(0) = 0\) as the initial condition of the fMRI signal at time \(t=0\). One can impose an analogous unit length constraint on the q function by obtaining \(\tilde{q} = \frac{q}{||q||}\). This unit length transformation forces q to lie on a Hilbert sphere denoted by \({\mathbb {S}_f}\). The space \({\mathbb {S}_f}\) is defined as \({\mathbb {S}_f} \equiv \left\{ q \in \mathbb {L}^2|\int _{0}^{1}(q(s),q(s))_{\mathbb {R}^2}ds=1,q(s): [0, 1] \rightarrow \mathbb {R}^2 \right\} .\)

Examples of fMRI hyperalignment for time courses in two voxels from the amygdala and the putamen respectively.

Full size image

2.1 Unbiased Within-Condition Time Series Template Estimation

To match the shape of the responses within conditions, we use the idea of time reparameterization to shift the signals to maximize peaks and troughs of functional activity. The time reparameterization function is represented by diffeomorphic function \(\gamma : I \rightarrow I\), where \(\dot{\gamma } > 0, \forall t \in I\). To change the time domain parameterization, one can compose f with \(\gamma \) as \({f \circ \gamma }\). In the SRVF domain, this is given by \(q \cdot \gamma = \frac{(\dot{f}\circ \gamma )\dot{\gamma }}{\sqrt{|(\dot{f}\circ \gamma )\dot{\gamma }|}} = (q\,\circ \,\gamma )\sqrt{\dot{\gamma }}\). To compute an unbiased optimal template that represents an average functional response related to the condition of the task, we exploit the geometry of the space of these time course representations by defining the notion of a metric on the space of q functions. The Riemannian metric on the tangent space of this sphere \(T_q( {\mathbb {S}_f} )\) assumes an elastic form (invariance to temporal reparameterization) and is given by \( \langle q_1, q_2 \rangle = \int _t q_1 q_2 dt\). The spherical geometry of the space of functions provides a convenient closed form solution for geodesics between two functions \(q_{1}\) and \(q_{2}\), which are given by \(\chi _{t}({q_{1}};{v})=\cos \left( t\cos ^{-1}\langle q_{1},q_{2}\rangle \right) {q_{1}}+\ \sin \left( t\cos ^{-1}\langle q_{1},q_{2}\rangle \right) {v}\), where \(t \in [0, 1]\) and the initial tangent vector \({v}\in T_{q_{1}}(\mathcal{Q})\) is given by \(v=q_{2}- \langle q_{1},q_{2} \rangle q_{1}\). Then the geodesic distance between the two shapes \(q_{1}\) and \(q_{2}\) on \({\mathbb {S}_f}\) is given by \(d(q_{1},q_{2})=\int _{0}^{1}\sqrt{\langle \dot{\chi _{t}},\dot{\chi _{t}}\rangle }dt\). Following, the optimal temporal alignment is obtained as \(\hat{\gamma } = \mathop {{\mathrm {argmin}}}\nolimits _{\gamma \in \varGamma } d(q_{1}, q_{2}\cdot \gamma )\) by a combination of dynamic programming and gradient descent. This optimal \(\hat{\gamma }\) yields the temporal reparameterization (alignment) between two time course signal portions.

Now for the time series, for a given task condition k, we find the Karcher mean [8] of the set \(\{ f^k_j \}\) as solving \(\mu _k = \mathop {{\mathrm {argmin}}}\nolimits _{q} d(q, q^k_j )^2\), which gives us the unbiased within-condition time series template. We align all the time courses for the condition to this mean \(\mu \) to achieve temporal alignment of all the partial time courses for that condition. Algorithm 1 outlines the procedure for performing within-subject blocked condition-wise optimal reparameterization for task based fMRI time courses.

Data: Imaging data was acquired from ten healthy subjects on a Siemens 3T PRISMA scanner. Task fMRI scans were acquired using a multiband (MB) EPI sequence (TE = 37 ms, TR = 800 ms, voxel size = 2 \(\times \) 2 \(\times \) 2 mm\(^3\), FA = 52\(^\circ \), MB accl. factor = 8) and a T1w structural scan (TE = 4.6 ms, TR = 9.9 ms, voxel size = 0.8 \(\times \) 0.8 \(\times \) 0.8 mm\(^3\), FA = 2\(^\circ \)) was also acquired. The functional task consisted of a previously validated face-matching paradigm [1], where four different types of images are presented to the subject: neutral, happy, and fearful faces, and objects. In this work, we exclusively considered all the 60 subcortical regions of interest (ROIs) extracted from a widely used automated segmentation protocol for our analysis [2]. After performing standard preprocessing and within-subject function to structure registration steps [11], Algorithm 1 was applied to each individual subject for four contrasts; happy > neutral, fearful > neutral, happy > fearful, and fearful > happy face stimuli. Figure 2 shows examples of fMRI hyperalignment for time courses from the amygdala and the putamen. Single subject GLM analyses were conducted using FEAT [11] and statistical maps were generated for both the original and the warped volumes.

Activated clusters in two different subjects for two different conditions (happy > neutral faces and fearful > neutral faces) for the original data (blue clusters) and for the warped data (orange clusters). The overlapping clusters between the two datasets are shown in green (e.g cluster b). Cluster locations: (a) bilateral putamen; (b) right putamen; (c) left caudate; (d) left caudate; (e) bilateral amygdala.

Full size image

Single Subject Analysis: Results show robust differences in the activation maps of all 10 individual subjects before and after alignment. As an example, Fig. 3 shows the activation maps for two different subjects for the original data (in blue) and for the warped data (in orange) for two different contrasts (happy > neutral faces and fearful > neutral faces).

Here, for the happy> neutral contrast, Subject 1 shows activated clusters in the right and left putamen (cluster label a in Fig. 3) for the warped data, but not for the original data. For the same subject, both the original and the warped data show a cluster in the right putamen (b in Fig. 3), but only the warped data shows an activated cluster in the left caudate (c in Fig. 3) for the fearful > neutral contrast. For subject 2, only the warped data revealed clusters in the left putamen for happy > neutral and in the bilateral amygdala for fearful > neutral. The original data did not contain any activated cluster for both contrasts. These results agree with previously published literature implicating the putamen in positive emotional processing, and the role of amygdala in negative facial emotional processing [1] and suggest greater signal to noise for detecting regional differences in BOLD response after applying the warping procedure. Table 1 reports the significant clusters represented in Fig. 3.

Within Group Mean Activations: To evaluate the group differences between the original data and the warped data we computed cluster thresholded averaged GLM maps across subjects for two contrasts (happy > neutral and fearful > objects) (Fig. 4). From these maps both the original data and the warped data show activated clusters in the cerebellum (cluster a in Fig. 4), in the right amygdala (cluster b in Fig. 4), in the bilateral putamen (cluster c in Fig. 4), and in the bilateral hippocampus (cluster d in Fig. 4). Interestingly only the warped data showed activated clusters in the bilateral caudate (cluster e in Fig. 4) as well as an extended cluster in the right putamen that followed its anatomical boundary.

Mean activation maps cluster thresholded for the original data (top row) and for the warped data (bottom row) for two different contrasts (happy faces > objects and fearful faces > objects). Cluster locations: (a) cerebellar cortex; (b) right amygdala; (c) bilateral putamen; (d) hippocampus; (e) right caudate.

Full size image

Paired Statistics Between Warped and Unwarped Images: Finally, paired t-tests were computed to evaluate the change in the number of clusters and the number of voxels activated before and after performing elastic alignment for the four different conditions (happy > neutral; fearful > neutral; happy > objects and fearful > objects). There was a significant increase (Bonferroni corrected) in the number of clusters and number of voxels activated after applying warping (Table 2).

We proposed an elastic time series matching approach to construct an unbiased condition-specific template of the BOLD shape patterns. While the warped data shows increased number of clusters and voxels, we want to caution against the possibility of registration and enhancement of noisy activations. This can be potentially validated with repeated randomization tests and test-retest reliability statistics in the future. Further, the Riemannian metric used here promotes both shrinking, stretching and bending. Experiments that control the parameters of this metric will also yield useful data addressing this question. Finally, additional testing of this method on event based fMRI paradigms or simple motor tasks that generate reliable activations will be necessary to understand the behavior and potential utility of this method.