Keywords

1 Introduction

Task-based functional magnetic resonance imaging (tfMRI) has been a powerful and popular noninvasive neuroimaging methodology for the study of brain activity patterns and cognitive behaviors of the human brain [1]. To model the very informative but complex tfMRI time series data, a variety of brain network reconstruction and hemodynamic response modeling techniques have been developed in the literature. These methods include model-driven methods such as the general linear model (GLM) [2], and data-driven methods like principal component analysis (PCA) [3], independent component analysis (ICA) [4] and sparse representation/dictionary learning methods [5,6,7]. In general, these methods reconstructed hundreds of meaningful functional brain networks and corresponding hemodynamic response patterns from tfMRI datasets and greatly advanced our understanding of the regularity of brain activities [2, 5].

However, the ability to represent and describe the tfMRI data still limited the performance of hemodynamic response patterns modeling. Therefore, developing a descriptive model that can sufficiently deal with diverse and complex tfMRI data, as well as large noises, is the key towards automatic, effective and accurate modeling of those hemodynamic response patterns in tfMRI data. Under tfMRI condition, participants need to participate in predefined sequential tasks during the whole scan session and the functional brain activity is modulated from the interactions of brain networks in different time periods, which quite coincides with the characteristics of recurrent neural network (RNN) models [8]. Therefore, it is straightforward and justified to adopt RNNs to explore and represent hemodynamic response patterns. RNNs are feed forward neural networks, which can use their internal memory units to process arbitrary sequences of inputs and model the sequential and time dependencies. They are connectionist models with the ability to selectively pass information across sequence steps. In order to characterize the tfMRI brain activities, a deep recurrent neural network (DRNN) model [9] was proposed to reconstruct the whole brain tfMRI signals from stimulus task design patterns. This framework not only identified typical brain networks by traditional methods (e.g., GLM), but also simultaneously obtain a variety of temporal brain activity patterns at multiple time scales. These results proved the great advantage of RNN model in charactering the temporal dependency signals in tfMRI data. However, the DRNN model still highly relies on the prior knowledge of task stimulus patterns which greatly limited the analysis power of the model.

In order to overcome current limitations in DRNN model, in this study, we proposed a novel unsupervised framework of deep recurrent autoencoder (DRAE) for modeling diverse and complex hemodynamic response patterns in tfMRI data. The basic idea is combing the DRNN model and autoencoder to automatically estimate the optimal task stimulus patterns of the tfMRI data and reconstruct the meaningful functional brain networks simultaneously. Autoencoder [10] is an unsupervised model that automatically learns a latent or compressed representation of the input data by minimizing the error between the input and its reconstruction. In this study, we take advantage of the autoencoder to automatically estimate the task stimulus patterns from the original tfMRI data. When the model is converged, the learned weight matrix between the FC layers and reconstructed signals represents the spatial distributions of functional brain networks underlying the tfMRI data and the output of the top RNN layer in decoding part represents the diverse and complex hemodynamic response patterns under the task condition. We adopted the motor tfMRI dataset of HCP 900 subjects release and the whole HCP Q1 release tfMRI datasets as test beds. Extensive experimental results demonstrated that the proposed DRAE model can not only automatically estimate the task stimulus patterns, but also reconstruct the meaningful functional brain networks and corresponding complex and concurrent hemodynamic response patterns with different time delays.

2 Materials and Methods

2.1 Data Acquisition and Pre-processing

The Human Connectome Project (HCP) dataset has been considered as one of the most systematic and comprehensive neuroimaging datasets. Importantly, this dataset is publicly available which makes it a good test bed for different research studies. The design paradigms are available in [11]. There are 68 subjects in HCP Q1 release dataset and over 800 subjects in HCP 900 subjects release dataset. The detailed acquisition parameters of these HCP tfMRI data are as follows: 220 mm FOV, in-plane FOV: 208 × 180 mm, flip angle = 52, BW = 2290 Hz/Px, 2 × 2×2 mm spatial resolution, 90 × 104 matrix, 72 slices, TR = 0.72 s, TE = 33.1 ms. The preprocessing of the task fMRI data sets includes skull removal, motion correction, slice time correction, spatial smoothing, and global drift removal (high-pass filtering). All these preprocessing steps were implemented in FSL FEAT. All of these individual fMRI datasets are first registered to the MNI common space for further study.

2.2 Deep Recurrent Autoencoder

The proposed deep recurrent autoencoder (DRAE) model is summarized as Fig. 1. It consists of two components, the encoder (Fig. 1(b)) and the decoder (Fig. 1(d)). First, for each subject, the extracted and normalized whole brain tfMRI signals are aggregated into a big signal matrix (m voxels’ signals with t time series, Fig. 1(a)). During the encoding stage, the signal matrix is compressed and mapped into a lower dimensional subspace representing a latent structure through a fully connected layer ([m, k], k < m), and then propagated through stacked RNN layers (k input units and n output units) to extract a feature map (n features with t time series, Fig. 1(c)). Next, the decoder passes the extracted feature map through another group of stacked RNN layers (n input units and k output units) to simulate diverse and complex brain activities and then maps the output of the top RNN layer into higher dimensional space (same as the original signals) by a fully connected layer ([k, m]) to reconstruct the whole brain signals (Fig. 1(e)). Specifically, the sequential output of each unit in the top RNN layer represents a temporal brain activity pattern and the corresponding weight vector in the fully connected layer which connects this unit to the reconstructed signals represents the spatial distribution of a functional brain network. The objective of the DRAE model is to minimize the reconstruction errors over all subjects of the training dataset, and the entire training progress is completely data-driven and unsupervised.

Fig. 1.
figure 1

Pipeline of the DRAE model. (a) Signal matrix of whole brain tfMRI signals. (b) The encoder which consists of a fully connected layer and deep RNN layers. (c) Extracted features maps. (d) The decoder which consists of deep RNN layers and a fully connected layer. (e) Matrix of reconstructed whole brain signals.

Full size image

2.3 Estimation of Hemodynamic Responses

In order to further explore the hemodynamic brain response patterns, for each feature, we replaced the feature map with a testing one (Fig. 2(a)) by keeping one selected feature a single impulse and setting the others to zeros, and propagated it through the trained decoder. The decoder was stimulated by the impulse and simulates the brain activities, and the output of each unit in the top RNN layer (Fig. 2(b)) represents the hemodynamic response of the corresponding functional brain network to the certain feature.

Fig. 2.
figure 2

A sketch map of exploring derive hemodynamic response patterns. (a) Testing feature map which keeps one feature a single impulse and set the others to zeros to stimulate the trained decoder. (b) Output of each unit in the top RNN layer representing the hemodynamic response patterns to the certain feature.

Full size image

3 Experimental Results

In this work, the training was applied on the DRAE model with 2 RNN layers of 32 LSTM units and a feature size of 16. To be specific, we extracted 244,341 voxels’ signals for the motor tfMRI dataset of HCP 900 subjects release and 223,945 voxels’ signals for the HCP Q1 release tfMRI dataset. During the training stage, all subjects’ signals were used during the training stage, since training on grouped subjects’ data will help avoid overfitting and either L1- or L2-norm regularization will increase the training loss rapidly, only MSE was taken as the loss function.

3.1 Interpretation of Feature Maps and Spatial Patterns

After training of the DRAE model, a group of feature maps can be obtained for each subject. That is, the whole brain activities can be divided and represented by several feature activities. Since individual feature maps are unique, we work out group feature maps by calculating the group-average values for further interpretation. Among these group-averaged feature maps, a few feature maps which are quite correlated with the task design patterns were identified, as shown in Fig. 3, which suggests that the DRAE model has the ability to extract the whole brain activities to a lower dimensional representation.

Fig. 3.
figure 3

Visualization of individual feature maps and group-averaged feature maps of motor task. Task design stimuli and their corresponding HRF responses are also shown for comparison.

Full size image

When the model is converged, we can also obtain a trained weight matrix of the fully connected layer from the decoder. Specifically, each vector of the weight matrix represents the spatial distribution of a typical functional brain network. Figure 4 illustrates a few identified functional brain networks on the motor task of HCP tfMRI datasets using the DRAE model. As shown, the DRAE model can identify similar functional brain networks for almost all task designs and high spatial overlap rates suggest the reliability of the proposed method. Specifically, the spatial overlap rate is defined as the intersection of the identified brain networks and corresponding GLM activation results. However, the DRAE model was trained in a completely unsupervised process and the functional brain networks can be obtained without prior knowledge of the task designs. These results demonstrates the superiority of the proposed method in modeling functional brain networks.

Fig. 4.
figure 4

Visualization of identified functional brain networks by the DRAE model. The GLM-derived group-wise activation maps are shown for comparison. M1–M6, S1 and S2 represent different task designs in motor and social tfMRI datasets, respectively. Last column O represents spatial overlap rate.

Full size image

3.2 Hemodynamic Response Patterns

After model training, the feature maps represent the task stimulus patterns of the whole brain signals and these signals can be reconstructed from this feature maps. The decoder simulates the complex brain activities with its hierarchical inner cells of stacked RNN units. Diverse hemodynamic response activities are invoked by these features, thus the output of each unit in the top RNN layer just represents a temporal response pattern.

When the formatted feature maps were passed through the decoder as described in Sect. 2.3, both positive and negative response patterns were obtained. These response patterns look similar to the theoretical HRF responses but have different shapes and time delays. Figure 5 detailed shows a few typical and HRF-correlated hemodynamic response patterns. There are several minor differences among these patterns, for example, they have different raising speed when meeting the impulse; some patterns start falling down after a period of time while some don’t fall down until the falling edge of the impulse; some patterns have significant undershoots while some don’t have. These hemodynamic response patterns are all possible and meaningful brain activity patterns which are more specific but still interpretable, which suggests that the proposed DRAE model can obtain many more meaningful brain activities through an unsupervised method on unlabeled datasets. In order to further analyze the hemodynamic response patterns, we drew the delay estimation maps trained on the DRAE model with different depths of RNN layers (Fig. 6). Table 1 represents the average values and standard deviations of the response delays. Since positive and negative responses have similar response time delays, the negative response patterns were inverted when calculating the response delays. For the DRAE model with just one RNN layer, peak delays are almost before the falling edge (12.29 s delay) of the impulse, since the network is very simple. As RNN depth goes deeper, response peaks have larger time delays around the falling edge and larger standard deviations. These response patterns also have various undershoots delays. In general, with deep recurrent layers, the DRAE model is able to estimate diverse and complex hemodynamic response patterns.

Fig. 5.
figure 5

A few typical and HRF-correlated hemodynamic response patterns.

Full size image
Fig. 6.
figure 6

Estimation of the hemodynamic response delays. X-axis represents delay of peaks, Y-axis represents delay of undershoots; red dots indicate positive response patterns, blue dots indicate negative response patterns; solid lines represent raising edge of the testing impulse and dashed lines represent falling edge of the testing impulse. Subfigures (a) to (f) show the results of 1, 2, 3, 4, 6 and 8 RNN layer(s), respectively. (Color figure online)

Full size image
Table 1. Statistic of the estimated hemodynamic response delays.
Full size table

4 Discussion and Conclusion

In this work, we proposed a novel framework of deep recurrent autoencoder (DRAE) for modeling diverse and complex hemodynamic response patterns and functional brain networks. The proposed DRAE model combines the deep recurrent neural network (DRNN) model and autoencoder to automatically estimate the optimal task stimulus patterns of the tfMRI data, reconstruct the meaningful functional brain networks and characterize the corresponding hemodynamic response patterns underlying tfMRI data simultaneously. Diverse and complex hemodynamic response patterns can be obtained, which brings a new way to reverse engineering of the brain’s response function patterns. Furthermore, with deeper stacked RNN layers, the DRAE model is able to simulate more complex hemodynamic response patterns with different time delay estimations. In general, extensive experiment results demonstrated the superiority and effectiveness of our proposed method.