Fast shared response model for fMRI data
Hugo Richard Lucas Martin Ana Luısa Pinho Jonathan Pillow Bertrand Thirion




1]Hugo Richardcorrespondingauthor
[correspondingauthor]Corresponding author
hugo.richard@inria.fr
1]Lucas Martin
1]Ana Luísa Pinho
2]Jonathan Pillow
1]Bertrand Thirion
[1]Parietal team, Inria, Neurospin, CEA, Université Paris-Saclay, France
[2]Princeton Neuroscience Institute, Princeton University, USA


The shared response model provides a simple but effective framework to analyse fMRI data of subjects 
exposed to naturalistic stimuli.

However when the number of subjects or runs is large, fitting the model requires a large amount of memory and computational power, 
which limits its use in practice. 

In this work, we introduce the FastSRM algorithm that relies on an intermediate atlas-based representation. 

It provides considerable speed-up in time and memory usage, hence it allows easy and fast large-scale
analysis of naturalistic-stimulus fMRI data.

Using four different datasets, we show  that our method matches the performance
of the original SRM algorithm while 
being about 5x faster and 20x to 40x more memory efficient.

Based on this contribution, we use FastSRM to predict age from movie watching data on the CamCAN 
sample.

Besides delivering accurate predictions (mean absolute error of 7.5 years), FastSRM extracts topographic
patterns that are predictive of age, demonstrating that brain activity during free perception reflects age.



	Functional magnetic resonance imaging Naturalistic imaging Shared response model Large scale analysis





§ INTRODUCTION

When exposed to naturalistic stimuli (e.g. movie watching), subjects'
experience is closer to their every-day life than with classical
psychological experiments.

This makes naturalistic paradigms an attractive class of
stimulation protocols for brain imaging.

While there is a broad interest in understanding how the brain reacts
in such ecological conditions, the recorded brain activity is
difficult to analyse.

Standards methods such as the general linear model <cit.> require the experimenter to construct a design matrix that models features of the presented stimuli across time. 

Such design matrices are notoriously difficult to construct for naturalistic stimuli as one has to rely on manual annotations (see <cit.>) or deep learning techniques (see e.g. <cit.>, <cit.>, <cit.> or <cit.>) that are hard to use, and provide high-dimensional, cumbersome models of the stimulus.

<cit.> has shown that brains exposed to the same natural stimuli exhibit synchronous activity.

The shared response model (SRM) <cit.> models this behaviour and extracts a common response from different subjects exposed to the same stimuli and subject-specific spatial components.

The resulting shared response can serve as a design matrix, while the spatial components are naturally seen as weighting factors. 


SRM has initially been designed to work within regions of interest using few subjects.

It has been used in <cit.> to transfer knowledge between subjects allowing precise location of a 15s time segment of fMRI data from a left-out subject using data from other subjects exposed to the same stimuli. 

A first big step towards more scalability has been made in <cit.> reducing fitting time and memory requirements by several orders of magnitude thanks to a smart use of the inversion lemma. 

After this improvement, studies using larger regions of interest have emerged such as <cit.>, which uses SRM to predict text embeddings from fMRI data or <cit.>, which shows that shared memories also come with shared structure in neural activity.

However, when using full brain data and a large number of subjects and runs, computational costs are still very high. Memory requirements are difficult to meet since all data have to be loaded in memory and since all full brain spatial components of all subjects are updated at each iteration, this leads to a heavy computational burden.


Fortunately, these high costs can be reduced. Intuitively, since the shared response lives in a reduced space, a compressed representation of the input is good-enough to find a suitable estimate.

With the use of off-the-shelf atlases and careful memory management we implement this idea and build FastSRM.


FastSRM is scalable with the number of runs and subjects. Fitting time and memory requirements are reduced considerably, making it possible to fit FastSRM using a laptop in a reasonable amount of time.





FastSRM makes large-scale analysis of movie-watching fMRI fast and easy. We demonstrate its usefulness on CamCAN data where we predict age from movie watching data with a mean absolute error (MAE) of 7.5 years.

We also show in an encoding experiment that FastSRM's ability to transfer data between subjects is superior to current implementations (<cit.>).  

Our code is freely available at <https://github.com/hugorichard/brainiak/tree/fastsrm>.


There is a long history of using latent factor models for fMRI data analysis.

ICA was first applied on fMRI data in <cit.> as an alternative to the generalized linear (GLM) model.

A few years later, <cit.> proposed Probabilistic ICA reducing overfitting problems of original ICA by introducing a Gaussian noise model and low rank structure.

In order to be able to compare patterns across subjects, ICA can be applied on time-wise concatenated data <cit.>. 

In 2010, <cit.> introduces CanICA that yields more stable group decompositions than previous group ICA methods (<cit.>, <cit.> and <cit.>).

Since then, a number of different factor models have been applied to fMRI data such as non-negative matrix factorization <cit.> or dictionary learning <cit.>.

All models enforce constrains on the data that are more or less realistic. <cit.> shows that total variation constrains on spatial components yield accurate and stable decompositions across subjects. However, with high quantities of data, the impact of such regularizations on the result vanishes (see <cit.>) and in this case one should therefore favor efficient algorithms (such as the online dictionary learning implementation of <cit.>).

In practice, most factor algorithms have been used to derive atlases from rest data. But other applications exist. In <cit.>, dictionary learning is used to derive an atlas from task data (contrast maps) and ICA is used in <cit.> and <cit.> to study respectively the effect of alcohol and speed on driving behavior thanks to a simulated driving protocol.

SRM has first been introduced in <cit.>. Like <cit.>, it learns subject-specific atlases and extracts a shared functional representation from the fMRI data of multiple subjects exposed to the same stimuli.

It relies on the hypothesis that orthonormality constrains put on spatial components allow for the extraction of more interpretable spatial components from naturalistic-stimulus fMRI data. 

SRM is not sensitive to the spatial variability of activated area between subjects since the only common features across participants are time courses making up the shared representation.

It is thus related to template estimation methods such as hyperalignment <cit.> or optimal transport <cit.>.

Many variants of SRM exists. <cit.> uses a convolutional autoencoder to preserve spatial locality, <cit.> uses matrix normal prior on spatial components and shared response,  <cit.> adds a subject-specific component and <cit.> gives a semi-supervised version of SRM.

All of these methods suffer from the efficiency issues already mentioned for standard SRM.

Note that, unless specified otherwise, whenever we refer to SRM, we refer to the current most efficient implementation of ProbSRM as provided by <cit.> which is freely available in the brainiak library <https://github.com/brainiak/brainiak>.

FastSRM needs atlases with a sufficient number of regions to work. Good atlases of various size are available in the literature such as Basc (up to 444 parcels) <cit.>, Shaeffer (up to 800 parcels) <cit.> or MODL (up to 1024 parcels) <cit.>.

We show that taking any such atlas yields quantitatively similar results.



§ MATERIAL AND METHODS



 §.§ The shared response model (SRM)

The shared response model is a latent factor model. The brain images of subject i during run s are stored in a matrix 𝐗^(s)_i ∈ℝ^t × v where v is the number of voxels and t the number of acquired brain images.

At time τ, the brain volume 𝐗^(s)_i[τ]  is modeled as a weighted sum of k orthonormal spatial components stored in 𝐖_i ∈ℝ^k × v.


Keeping things simple, we assume that all n subjects have the same
number v of voxels and all m runs have the same number t of
timeframes.

The extension to the more general case where each run has its own
number of timeframes and each subject its own number of voxels is
straightforward. In our implementation runs can have different number
of timeframes but each subject has the same number of voxels.

Typical values for the number of voxels v, the number of timeframes
t, the number of runs m, the number of subjects n and the
number of components k are given in Table <ref>.

Let us introduce the following notations:



  * For each subject i, the concatenation 𝐗_i ∈ℝ^mt × v of the acquisition data for all runs:

    𝐗_i = 
    [ 𝐗^(1)_i; 𝐗^(2)_i;       ⋮; 𝐗^(m)_i ]



  * The concatenation 𝐗∈ℝ^mt × nv of the brain acquisition of all subjects for all runs:

    𝐗 = 
    [ 𝐗^(1)_1 𝐗^(1)_2       ⋯ 𝐗^(1)_n; 𝐗^(2)_1 𝐗^(2)_2       ⋯ 𝐗^(2)_n;       ⋮       ⋮       ⋱       ⋮; 𝐗^(m)_1 𝐗^(m)_2       ⋯ 𝐗^(m)_n;         ]




  * The concatenation 𝐖∈ℝ^k × nv of the spatial components of all subjects:

    𝐖 = 
    [ 𝐖_1 𝐖_2   ⋯ 𝐖_n ]

where 𝐖_i∈ℝ^k × v contains the k spatial components of subject i.


  * The concatenation of the weights 𝐒∈ℝ^mt × k of all runs:

    𝐒 = 
    [ 𝐒^(1); 𝐒^(2);     ⋮; 𝐒^(m) ]

where 𝐒^(s)∈ℝ^t × k contains the k weights of run s across time.  




Formally SRM is defined by: 

    𝐗 = 𝐒𝐖 + 𝐄,

Where the spatial components 𝐖, the shared response 𝐒 and the noise 𝐄 need to be estimated from 𝐗.

An illustration of this definition is given in Figure <ref>.

Framed this way, group versions of dictionary learning, ICA, blind signal separation or matrix factorization can be seen as particular instances of shared response model.

Most versions of SRM impose orthonormal constrains on spatial components of each subject:

    ∀ i ∈{1, ⋯, n} 𝐖_i 𝐖_i^T = 𝐈_k

where 𝐈_k ∈ℝ^k × k is the identity matrix of size k.






 §.§ Deterministic SRM model (DetSRM)

The deterministic SRM model assumes Gaussian noise with the same variance for all subjects and orthonormal spatial components.

Formally the model reads:

    ∀ i ∈{1 ⋯ n},  𝐗_i[τ] ∼𝒩(𝐒[τ] 𝐖_i, σ^2 𝐈_v)  such that 𝐖_i 𝐖_i^T = 𝐈_k,


Maximizing the log-likelihood we obtain the following optimization problem:

    min_𝐖_1, ⋯, 𝐖_n, 𝐒∑_i=1^n ||𝐗_i - 𝐒𝐖_i||^2 
       such that ∀ i ∈{1,⋯, n} 𝐖_i 𝐖_i^T = 𝐈_k


This can be solved efficiently using alternate minimization on 
(𝐖_1, ⋯ , 𝐖_n ) and 𝐒. At each iteration we have two problems to solve that have a closed form solution:

    ∀ i ∈{1 ⋯ n}, 
         argmin_{𝐖_i  such that 𝐖_i𝐖_i^T = 𝐈_k}∑_j=1^n ||𝐗_j - 𝐒𝐖_j||^2 = 𝐔_i 𝐕_i,  
         where 𝐔_i, 𝐃_i, 𝐕_i = SVD(𝐒^T 𝐗_i)

where SVD stands for singular value decomposition, and

    argmin_𝐒∑_i=1^n ||𝐗_i - 𝐒𝐖_i||^2 = 1/n∑_i=1^n 𝐗_i 𝐖_i^T


Assuming k ≪ v, t, the time-complexity of this approach such as implemented in the brainiak library is in O(n_iter nmtvk) and storage requirements are in O(nmvt). This means that the method becomes expensive whenever the number of subjects n or runs m becomes large.



 §.§ Probabilistic SRM model (ProbSRM)

In the probabilistic SRM, spatial components are assumed orthonormal, the shared response is modeled by its covariance matrix Σ∈ℝ^k × k and the variance of the Gaussian noise σ_1, ⋯ , σ_n is assumed different for different subjects.

In the original paper <cit.>, an intercept is also learned but since we remove the mean of each time-course as a preprocessing step, it is of no use here.

Formally the model reads:

    𝐒[τ] ∼𝒩 (0_k, Σ) 
       ∀ i ∈{1 ⋯ n}, 
         𝐗_i[τ] ∼𝒩 (𝐒[τ]𝐖_i, σ_i^2 𝐈_v) 
         such that  𝐖_i𝐖_i^T = 𝐈_k


The optimization is done using an expectation maximization algorithm described in <cit.> and <cit.>. Assuming k ≪ v, t, the time-complexity of this approach such as implemented in the brainiak library is O(n_iter k nv tm) and storage requirements are in O(nmvt). So ProbSRM is also expensive whenever the number of subjects n or runs m becomes large.



 §.§ FastSRM model

When we deal with large datasets (n, m, v, t are large), above implementations require huge computational power. We introduce FastSRM, a fast and memory-efficient algorithm.

In a first step, we project the data 𝐗∈ℝ^tm × nv onto a chosen atlas 𝐀 with c parcels.

The atlas 𝐀∈ℝ^c × v can be either probabilistic or a strict partition of the set of voxels but the number of parcels c of the atlas should be larger than the number of components k of the FastSRM model and small compared to the number v of voxels.

In typical settings c reaches hundreds to one thousand.

Projection onto the atlas yields reduced data 𝐗̂∈ℝ^tm × nc.


    ∀ i ∈{1 ⋯ n}, s ∈{1 ⋯ m } 𝐗̂^(s)_i = 𝐗^(s)_i 𝐀^T (𝐀𝐀^T)^-1


In cases where the atlas is a partition, projecting onto the atlas is equivalent to averaging brain activation in each parcel of the atlas. 

Formally, denoting R_j the parcel j of atlas 𝐀, we compute the projection over the atlas by: 

    ∀ i ∈{1 ⋯ n}, s ∈{1 ⋯ m } 𝐗̂^(s)_i[τ, j] = ∑_x ∈ R_j𝐗^(s)_i[τ, x]/|R_j|


In a second step we apply our preferred SRM algorithm on the reduced data to find the shared response in reduced space 𝐒̂ (in our implementation we use a deterministic SRM). Since c is small compared to v, this step is very fast even if the number of iterations is high.
 

    𝐖̂, 𝐒̂ =DetSRM(𝐗̂)

where the reduced maps 𝐖̂ and the shared response in reduced space 𝐒̂ are output of the deterministic SRM algorithm (DetSRM).

The spatial components of each subject are recovered by orthonormal regression using the shared response in reduced space 𝐒̂ and the data 𝐗:

    ∀ i, 𝐖_i= 𝐔_i 𝐕_i  
       where 𝐔_i, 𝐃_i, 𝐕_i = SVD(∑_s=1^m 𝐒̂^(s)^T𝐗^(s)_i )


In FastSRM as well as in ProbSRM and DetSRM, fitting the algorithm only means learning the spatial components. If a temporal model is needed it has to be recomputed a posteriori.

In order to compute the shared response from fMRI data of n subjects 𝐗^(s)_1 ⋯𝐗^(s)_n in a particular run s, one has just to average the projection of the data onto the basis of each subject:

    𝐒^(s) = 1/n∑_i=1^n 𝐗^(s)_i 𝐖_i^T


Note that when the shared response of train data is needed one cannot use directly the FastSRM shared response  𝐒̂ in reduced space, since it does not have the right scale.

However, the spatial components obtained by orthonormal regression from the ill-scaled shared response are valid. 

Indeed if we multiply the shared response by a scaling factor f, for any subject i the singular value decomposition of ∑_s=1^m f Ŝ^(s)^T𝐗̂^(s)_i is given by 𝐔_i f 𝐃_i 𝐕_i where 𝐔_i 𝐃_i 𝐕_i is the singular decomposition we would obtain with f=1 and therefore the spatial components of subject i are given by 𝐖_i = 𝐔_i 𝐕_i which is independent from the scale factor.

The orthonormal projection can easily be replaced by any kind of regression. One can easily impose sparsity, non-negativity, smoothness or similarity constrains on spatial components instead of orthonormality to obtain more interpretable patterns.

This however leads to more heavier computations.

In our implementation, we use orthonormal regression.

Assuming k ≪ c ≪ v, t, the bottleneck is the projection onto the atlas which yields a time complexity in O(n m t v c).

We need to keep spatial components, the atlas as well as one run of one subject in memory which yields a memory complexity of O(nvk + vc + vt)).

If parallelization is used memory requirements become O(nvk + vc + n_jobs vt)) and fitting time becomes O(nmtvc/n_jobs) where n_jobs is the number of jobs. 

It is also possible to write spatial components on disk instead of keeping them in memory which adds about O(vk) read/write operations but reduces the memory complexity to O(v(c + n_jobst)). In our experiments we write spatial components on the disk.


 




 §.§ Experiments



  §.§.§ Datasets


We use five fMRI datasets of subjects exposed to naturalistic stimuli. When needed, datasets are preprocessed with FSL <http://fsl.fmrib.ox.ac.uk/fsl> using slice time correction, spatial realignment, coregistration to the T1 image and affine transformation of the functional volumes to a template brain (MNI). Using nilearn <cit.>, preprocessed data are resampled, masked (using a full brain mask available at <http://cogspaces.github.io/assets/data/hcp_mask.nii.gz>), detrended and standardized after a 5 mm smoothing is applied.  



  
SHERLOCK
In SHERLOCK dataset, 17 participants are watching "Sherlock" BBC TV show (episode 1). 

These data are downloaded from <http://arks.princeton.edu/ark:/88435/dsp01nz8062179>. 

Data were acquired using a 3T scanner with an isotropic spatial resolution of 3 mm. 

More information including the preprocessing pipeline is available in <cit.>.

Subject 5 is removed because of missing data leaving us with 16 participants.

Although SHERLOCK data contains originally only 1 run, we split it into 4 runs of 395 timeframes and one run of 396 timeframes for the needs of our experiments. 




  
FORREST
In FORREST dataset 20 participants are listening to an audio version of the movie Forrest Gump.

FORREST data are downloaded from OpenfMRI <cit.>. 

Data were acquired using a 7T scanner with an isotropic spatial resolution of 1 mm (see more details in <cit.>.

More information about the forrest project can be found at <http://studyforrest.org>.

Subject 10 and run 8 are discarded because of missing data.

We therefore use full brain data of 19 subjects split in 7 runs of respectively 451, 441, 438, 488, 462, 439 and 542 timeframes.
 


  
RAIDERS
In RAIDERS dataset, 10 participants are watching the movie "Raiders of the lost ark".

The RAIDERS dataset pertains to the Individual Brain Charting dataset (<cit.>).

They acquired at NeuroSpin using a 3T scanner with an isotropic spatial resolution of 3 mm.

The RAIDERS dataset reproduces the protocol described in <cit.>.

We use full brain data of 10 subjects split in 9 runs of respectively 374, 297, 314, 379, 347, 346, 350, 353 and 211 timeframes.



  
CLIPS
In CLIPS dataset, 10 participants are exposed to short clips. 

The CLIPS dataset also pertains to the Individual Brain Charting dataset (<cit.>).

It reproduces the protocol of original studies described in <cit.> and <cit.>.

In our experiments we use the data of 10 participants acquired in 17 runs of 325 timeframes.

At the time of writing, the CLIPS and RAIDERS dataset from the individual brain charting dataset <https://project.inria.fr/IBC/> are not yet public, but they will be in the future. Protocols on the visual stimuli presented are available in a dedicated repository on Github: <https://github.com/hbp-brain-charting/public_protocols>. The informed consent of all subjects was obtained before scanning.



  
CamCAN
In CamCAN dataset, 647 participants aged from 18 to 88 years are watching Alfred Hitchcock's "Bang! You're Dead" (edited so that it lasts only 8 minutes).

CamCAN consists of data obtained from the CamCAN repository (available at <http://www.mrc-cbu.cam.ac.uk/datasets/camcan/>) (see <cit.> and <cit.>).

We use all available subjects and runs yielding 647 participants and 1 run of 193 timeframes.

A summary about the size of each dataset is available in Table <ref>.




  §.§.§ fMRI reconstruction: Evaluate the ability to recover BOLD signal on left-out runs



When subjects are exposed to the same stimuli, SRM algorithms posit that the recorded fMRI data can be modeled as a product of two matrices, one of which is fixed across time but subject-specific (the spatial components) while the other varies across time but is common to all subjects (the shared response).

Under this framework, once spatial components are known, we can generate an accurate estimation of the data of one subject given the data of all others.

In this experiment we test whether we can recover data of a left-out subject using previous data of the same subject as well as data from other subjects. 
We denote 𝐗^(-s) brain recordings of all runs but run s:
 
    𝐗^(-s) = 
    [   𝐗^(1);       ⋮; 𝐗^(s-1); 𝐗^(s+1);       ⋮;   𝐗^(m) ]


Similarly, 𝐗_i^(-s) refers to the brain recordings of subject i using all runs but run s. The data from all subjects but i acquired during run s are denoted:

    𝐗^(s)_-i = 
    [   𝐗_1^(s)         ⋯ 𝐗^(s)_i-1 𝐗^(s)_i+1         ⋯   𝐗^(s)_n ]


We evaluate our model using a cross validation scheme known as co-smoothing (see <cit.>). First, all brain recordings for all runs but one 𝐗^(-s) are used for learning subjects spatial components 𝐖^(-s). The exponent in 𝐖^(-s) indicates that these spatial components were learned using all runs but run s. 


    𝐖^(-s), 𝐒^(-s) = SRM(𝐗^(-s))


Then we focus on the left-out run s and use all subjects but one 𝐗^(s)_-i to compute a shared response for the left-out run 𝐒^(s).


    𝐒^(s) = ∑_z=1, z≠ i^n 𝐗^(s)_z𝐖^(-s)^T_z/n - 1


From 𝐒^(s) and 𝐖_i^-s we compute 𝐗^(s)_i which stands as an estimate of brain activity of the left-out subject i during the left-out run s.


    𝐗_i^(s) = 𝐖^(-s)_i 𝐒^(s)






An illustration of our reconstruction experiment is available in Figure <ref>.

The performance is measured voxel-wise using the R^2 score between 𝐗_i^(s) and 𝐗_i^(s) as a similarity measure.
For any two time-courses x ∈ℝ^t and y ∈ℝ^t we define the R^2 score by:


    R^2(x, y) = 1 - ∑_z=1^t (x[z] - y[z])^2/∑_z=1^t (y[z] - y[z])^2


Where y = 1/t∑_z=1^t y[z].

Following the leave-one-out cross validation scheme, all our experiments are done several times with a different left-out subject to reconstruct. 

We measure the average R^2 score across  all left-out subjects. Note that we obtain one such value per voxel.



  §.§.§ Predict age from spatial components

Since spatial components are subject-specific they should be predictive of subject-specific features such as age.

In this experiment we try to predict subject's age from movie-watching data using SRM algorithms.

Functionally matched spatial components are obtained using an SRM algorithm.

They are divided into two groups (train and test data) where the train set contains 80% of the data and the test set 20%.

Within the train set we split again our data into two groups: the first group is used to train one Ridge model per spatial components, the second group is used to train a Random Forest to predict age from Ridge predictions. This way of stacking models is similar to the pipeline used in <cit.>.

We use 5 fold cross validation to split the train set (so that the number of samples used to train the Random Forest is the number of elements in the train set).

Then the train set is used to train one Ridge model per spatial component.

On the test set each Ridge model makes a prediction and the predictions are aggregated using the Random Forest model.

An illustration of the process is available in Figure <ref>. 

In each Ridge model, the coefficient that determines the level of l2 penalization is set by generalized cross validation, an efficient form of leave-one-out cross validation (see the RidgeCV implementation of Scikit Learn <cit.>).

The train and test sets are chosen randomly. 5 different choices for the train and test set are made. We report the average mean absolute error (MAE) on the test set averaged over the 5 splits.






§ RESULTS AND DISCUSSION




 §.§ fMRI Reconstruction

We perform the reconstruction experiment on the FORREST, CLIPS, RAIDERS and SHERLOCK datasets. We compare brainiak's implementation of ProbSRM / DetSRM to our implementation of FastSRM in terms of fitting time, memory usage and performance. In order to be fair, we do not use parallelization (n_jobs = 1) and we set the number of iterations to 10 (n_iter=10) which is ProbSRM's default. Note that the fitting time of ProbSRM is roughly proportional to the number of iterations while it has a limited impact on the fitting time of FastSRM. 


We run our experiments on the full brain and report an R^2 score per voxel.

However we measure the performance in terms of mean R^2 score inside a region of interest (in order to leave out regions where there is no useful information).

In order to determine the region of interest, we focus on the results of ProbSRM with 10, 20, 50 and 100 components and keep only the intersection of regions where the R^2 score is above 0.05.

This means of selecting regions favors ProbSRM. For completeness, full brain R^2 images obtained on the four datasets with ProbSRM and FastSRM using 20 components averaged across subjects are available in Figure <ref>.

In Figure <ref>, we plotted the mean R^2 score against the number of components (k) for ProbSRM, DetSRM and FastSRM algorithm with different atlases.

The R^2 score tends to increase with the number of components (which is what is expected as more information can be retrieved when the number of components is high).

FastSRM matches the performance of ProbSRM and DetSRM. This holds for any atlas we chose (BASC (444 parcels), SHAEFFER (800 parcels), MODL (512 and 1024 parcels)) and for all datasets we tried (SHERLOCK, RAIDERS, CLIPS and FORREST). 

In Figure <ref>, we compare the running time of FastSRM, DetSRM and ProbSRM on the four different datasets.

FastSRM is on average (across datasets) about 5 times faster.

On the FORREST dataset, we compute a shared response in about 3 minutes when it takes about 20 minutes with ProbSRM or DetSRM.

In Figure <ref>, we compare the memory (RAM) consumption of FastSRM, DetSRM and ProbSRM on the four different datasets.

FastSRM is 20 to 40 times more memory friendly than ProbSRM and 10 to 20 times more memory friendly than DetSRM. On the FORREST dataset the memory usage of FastSRM is between 1 and 3 Go depending on the number of components and the atlas used. Most modern laptops meet these requirements. On the same dataset memory consumption is about 80 Go for ProbSRM and 40 Go for DetSRM which is manageable but costly for small labs.

Overall FastSRM yields same performance as ProbSRM and DetSRM while being much faster and using far less memory. We also show that the atlas used to reduce the data only has a minor impact on performance.







 §.§ Predicting age from spatial components

Because FastSRM is fast and memory efficient, it enables large-scale analysis of fMRI recordings of subjects exposed to the same naturalistic stimuli.

We use all 647 subjects of the CamCAN dataset and demonstrate the usefulness of FastSRM by showing that the spatial components it extracts from movie watching data are predictive of age.

A key asset of FastSRM is that these spatial components can be visualized and therefore provide meaningful insights. 

Figure <ref> shows that FastSRM predicts age with a good accuracy (better than ProbSRM and a lot better than chance) resulting in a mean absolute error (MAE) of 7.5 years.

It also shows that on CamCAN data, FastSRM is 4x faster and more than 150x more memory efficient than ProbSRM. As before and in order to ensure fair comparison the number of iterations is set to 10 and we do not make use of parallelization. Note that the memory requirements of ProbSRM on the CamCAN dataset (186Go) make it difficult to use.

FastSRM does not suffer from memory issues, making it suitable to analyse big datasets.

A key asset of our pipeline is that we can see which spatial components are most predictive of age by using feature importance.

Feature importance is assessed by the Gini importance defined in <cit.> or <cit.>.

It measures for each feature the relative reduction in Gini impurity brought by this feature.

Feature importance varies with different splits. We use the averaged feature importance over the 5 splits of our pipeline.

In Figure <ref> are shown the 3 most important spatial components representing respectively 16%, 12% and 8% of total feature importance.

These spatial components in decreasing order of importance represent the visual dorsal pathway, the precuneus and the visual ventral pathway. 

The fact that averaged spatial components are interpretable and meaningful allows us to study the influence of age on brain networks involved in movie-watching.

In Figure <ref>, we plot the most important spatial component averaged within groups of ages.

We see that these spatial components evolve with age allowing us to visually identify which regions are meaningful. 

It turns out that aging is mostly reflected in brain activity as a
fading of activity in the spatial correlates of movie watching,
particularly in the dorsal visual cortex.












 §.§ Conclusion

As studies using naturalistic stimuli will tend to become more common and large within and across subjects, we need scalable models especially in terms of memory usage.

This is what FastSRM provides.

We show that while FastSRM matches the performance of ProbSRM or DetSRM, it is
significantly faster and requires a lot less memory. While FastSRM's scalability relies on the use of atlases to compress the BOLD signal, we show that the precise choice of the atlas has only marginal effects on performance.


FastSRM allows large scale analysis of fMRI data of subjects exposed to naturalistic stimuli. As one example of such analysis, we show that it can be used to predict age from movie-watching data. Interestingly, although FastSRM is an unsupervised model, it extracts meaningful networks and as such constitutes a practical way of studying subjects exposed to naturalistic stimuli.

We also show that individual information can be extracted from the fMRI activity when subjects are exposed to naturalistic stimuli. Our predictive model is reminiscent of that of <cit.>, that have shown that ICA components obtained from the decomposition of resting state data carry important information on individual characteristics. 


As a side note, we chose to keep the orthonormality assumptions of the original SRM model but slight modifications of our implementation of FastSRM would allow one to build more refined model promoting sparsity, non-negativity or smoothness of spatial components for example.

The remaining difficulty with SRM is to interpret the spatio-temporal decomposition. Reverse correlation <cit.> can be used to clarify the cognitive information captured in the shared response.

Our code is freely available at <https://github.com/hugorichard/brainiak/tree/fastsrm>.






§ APPENDICES

In Table <ref>, we show typical values for the main dataset parameters. In Table <ref>, we describe the main dataset parameters of the real datasets we used. In Figure <ref> we compare the memory usage of DetSRM, ProbSRM and FastSRM with different atlases as a function of the number of components used while performing the encoding experiment described in Section <ref>. In Figure <ref>, we show for the same experiment, the R^2 score per voxels averaged across cross validation folds. 





















§ ACKNOWLEDGMENT

CamCAN Data collection and sharing was provided by the Cambridge Centre for Ageing and Neuroscience (CamCAN). CamCAN funding was provided by the UK Biotechnology and Biological Sciences Research Council (grant number BB/H008217/1), together with support from the UK Medical Research Council and University of Cambridge, UK.

This project has received funding from the European Union's Horizon 2020 Research and Innovation Programme under Grant Agreement No. 785907 (HBP SGA2).




§ REFERENCES


