%\documentclass[final,authoryear ,5p,times,9pt, twocolumn]{elsarticle} % for final
\documentclass[authoryear,times,9pt,onecolumn]{elsarticle}% for working copy
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphics}
\usepackage{amsthm}
\usepackage{amsmath}
\usepackage{longtable}
\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{pifont}
\usepackage{tikz}
\usepackage[]{algorithm2e}
\usepackage{todonotes}
\usepackage{pdflscape}
\usepackage{floatpag,mwe}
\usepackage{changepage}
\usepackage{adjustbox}
\usepackage{geometry}
\usepackage{parskip}
\usepackage{hyperref}
\usepackage{todonotes}
\graphicspath{{figures/}}
\newcommand{\cmark}{\ding{51}}
\newcommand{\xmark}{\ding{55}}
\def\x{{\mathbf x}}
\def\L{{\cal L}}
\def\str{1}
\DeclareMathOperator*{\argmin}{\arg\!\min}
\DeclareMathOperator*{\argminB}{argmin}
\DeclareMathOperator*{\maxi}{max}
\DeclareMathOperator*{\corr}{corr}
\DeclareMathOperator*{\dist}{dist}
\newcommand{\X}{{\mathbf{X}}}
\newcommand{\xu}{\mathbf{x^{u}}}
\newcommand{\xl}{\mathbf{x^{l}}}
\newcommand{\bxu}{\mathbf{tilde{x}^{u}}}
\newcommand{\bxl}{\mathbf{tilde{x}^{l}}}
\usepackage{lineno}
\journal{NeuroImage}
\begin{document}
\begin{frontmatter}
\title{TbCAPs: A ToolBox for Co-Activation Patterns Analysis}
\author[label1,label2]{Thomas AW Bolton \corref{cor1}}
\author[label3]{Constantin Tuleasca}
\author[label4]{Gwladys Rey}
\author[label5]{Diana Wotruba}
\author[label4]{Julian Gavirina}
\author[label6]{Herberto Danis}
\author[label6]{Eva Blondiaux}
\author[label6]{Baptiste Gauthier}
\author[label6]{Lukasz Smigielski}
\author[label1,label2]{Dimitri Van De Ville}
\address[label1]{Institute of Bioengineering, \'{E}cole Polytechnique F\'{e}d\'{e}rale de Lausanne (EPFL), Lausanne, Switzerland}
\address[label2]{Department of Radiology and Medical Informatics, University of Geneva (UNIGE), Geneva, Switzerland}
\address[label3]{Centre Hospitalier Universitaire Bic\^{e}tre, Service de Neurochirurgie, Paris, France}
\address[label4]{Department of Neuroscience, University of Geneva (UNIGE), Geneva, Switzerland}
\address[label5]{Program for Sustainable Development of Mental Health, University of Zurich, Zurich, Switzerland}
\address[label6]{Laboratory of Cognitive Neuroscience, \'{E}cole Polytechnique F\'{e}d\'{e}rale de Lausanne (EPFL), Lausanne, Switzerland}
\cortext[cor1]{thomas.bolton@epfl.ch}
\begin{abstract}
Functional magnetic resonance imaging provides rich spatiotemporal data of human brain activity during task and rest. Many recent efforts have focussed on characterizing dynamics of brain activity. One notable instance is co-activation patterns (CAPs) analysis, a frame-wise analytical approach that disentangles the different functional brain networks interacting with a user-defined seed region. While promising applications in various clinical settings have been demonstrated, there is not yet any centralised, publicly accessible resource to facilitate the deployment of the technique.
Here, we release a working version of TbCAPs, a new toolbox for CAPs analysis, which includes all steps of the analytical pipeline, introduces new methodological developments that build on already existing concepts, and enables a facilitated inspection of CAPs and resulting metrics of brain dynamics. The toolbox is available on a public academic repository \url{https://c4science.ch/source/CAP_Toolbox.git}.
In addition, to illustrate the feasibility and usefulness of our pipeline, we describe an application to the study of human cognition. CAPs are constructed from resting-state fMRI using as seed the right dorsolateral prefrontal cortex, and, in a separate sample, we successfully predict a behavioural measure of continuous attentional performance from the metrics of CAPs dynamics (R=0.59) .
Functional magnetic resonance imaging (fMRI) has enabled to track temporal changes in activity levels at the whole-brain scale by means of the blood oxygenation level-dependent (BOLD) contrast, a proxy for neural activation~\citep{Logothetis2001}. In addition to more traditional task-based studies in which BOLD changes are mapped to a paradigm of interest~\citep{Friston1994}, the characterisation of statistical interdependence between remote brain locations (termed \textit{functional connectivity}~\citep{Friston1994b}) in the resting-state, and the concomitant definition of large-scale \textit{resting-state brain networks} (RSNs), has been a popular endeavour~\citep{Biswal1995,Fox2005,Damoiseaux2006,Power2011}, with great benefits for the understanding of cognition and disease~\citep{VanDenHeuvel2010,Greicius2008,Fox2010}.
Over the past years, it has become increasingly appreciated that cross-regional relationships do not remain static over the course of a full scanning session~\citep{Chang2010}: instead, a given region rearranges its interactions along time, in ways that have been addressed with very diverse analytical tools (see~\citet{Preti2017} for an exhaustive review of the \textit{dynamic functional connectivity} field).\todo{Also add Hutchinson dynFC review here}
In one family of approaches that has been developed, it is assumed that only few salient time points contain the information of interest that shapes whole-brain correlational relationships; selecting only these frames, by means of a seed-based thresholding process, already enables to derive accurate RSN maps, even if as little as 10\% of data points is retained~\citep{Tagliazucchi2012b}. The analysis then moves from a second-order correlation-based characterisation to a first-order activation viewpoint, and reduces computational load, a desirable feat in light of the numerous large-scale acquisition initiatives embraced by the fMRI community~\citep{VanEssen2013,Nooner2012,Holmes2015}.
Building on this point-process analysis concept, and inspired by the dynamic viewpoint on resting-state brain function,~\citet{Liu2013} hypothesised that at different moments in time, the seed region of interest would display distinct interactions with the rest of the brain. A k-means clustering step was thus appended to frame selection, so that fMRI volumes with a large enough seed activity would be partitioned into a limited set of \textit{co-activation patterns} (CAPs).
Since then, co-activation pattern analysis has started to gain momentum as a potent tool to reveal functional brain dynamics subtleties: analyses taking the posterior cingulate cortex (PCC) as a seed revealed alterations of spatial intensity level and occurrence in specific CAPs~\cite{Amico2014,DiPerri2018}, while in adolescent depression,~\citet{Kaiser2019} showed that the time spent in a specific frontoinsular-default network CAP positively correlated with symptoms severity. CAPs analysis also enabled to track the renormalisation of CAP occurrences in patients with essential tremor following surgical intervention~\citep{Tuleasca2019}.
In parallel to clinical applications, the technical details of the approach have also been addressed, in terms of retaining activation versus deactivation time points~\citep{Di2013}, extending the approach to the whole brain~\citep{Liu2013b}, designing novel metrics of interest~\citep{Chen2015}, or constraining the extent of spatial overlap across CAPs~\citep{Zhuang2018}. For more details, the reader is pointed at the recent review of~\citet{Liu2018}.
Here, we wish to further foster the development of CAPs analysis by releasing a dedicated toolbox, which enables to easily navigate through the steps of the analytical pipeline through a graphical user interface, and also offers additional technical developments regarding frame selection and metrics computation. While the mathematical underpinnings of CAPs analysis are relatively straightforward, we hope that providing such a resource will encourage practitioners to embrace the method, and it will become easier to compare CAPs analyses based on subtle, but sometimes important, differences in the processing pipeline.
Further, to exemplify the use of our toolbox, we describe an application of CAPs analysis in the yet unaddressed setting of predicting cognitive skills: in a battery of healthy individuals, we show that continuous performance in a visual attention and vigilance task correlates with the expression profile of task-positive network (TPN) CAPs.
Let us consider the data matrix $\mathbf{X_s}\in\mathbb{R}^{V \times T}$ for subject $s$, where $V$ is the number of voxels to consider in the analysis and $T$ the number of time points. Each voxel-wise time course is temporally z-scored, so that $\mu_i=\frac{\sum_{t=1}^{T}X_{s}(i,t)}{T}=0$ and $\sigma_i=\sqrt{\frac{\sum_{t=1}^{T}(X_{s}(i,t)-\mu_i)^2}{T-1}}=1$, for all $i=1,2,\cdots,V$.
Co-activation patterns analysis requires the definition of a seed region, whose interactions with the rest of the brain will be probed. Formally, a set of voxels $\mathcal{S}$ that one wishes to consider is specified, and a time point $t$ of the seed activation time course is then given by:
Only time points when the seed time course takes sufficiently extreme values (denoting significant seed (de)activation) are considered. Let the activation threshold be $T$, we then construct the set $\mathcal{T}_s$ of time points that satisfy $S_{s}(t)>T$ (if we wish to consider solely activation moments) or $S_{s}(t)<T$ (if we are interested in seed deactivation time points).\todo{Rather $<-T$ for deactivation?}
In this work, in addition to the above, we propose an extension in which more than one seed region can be considered: for each seed $j$, a set of time points $\mathcal{T}_{s,j}$ is derived. Assuming $J$ separate seeds, one can then consider the time points when all seed time courses jointly take extreme values:\todo{Should be the other way around: here intersection, and for the next one, union?}
Finally, other additional criteria can be incorporated at the time point selection step: for instance, given the deleterious impact that head motion exerts on BOLD signals even following standard preprocessing~\citep{Power2012,VanDijk2012,Satterthwaite2012}, it may be desirable to only retain the frames for which framewise displacement does not exceed a threshold $T_{m}$.
After having selected the frames to keep for each subject, the next step is the population-level clustering of data points into CAPs. K-means clustering is used for this purpose, to optimise:
where $K$ is the number of co-activation patterns to derive, $\mathbf{\mathcal{C}}=\{\mathcal{C}_1,\cdots,\mathcal{C}_K\}$ summarises the hard assignment of the frames to each CAP, and $\mathbf{c}_k$ is the spatial map for co-activation pattern $k$. The $\dist$ function depends on the type of distance to use in the algorithm. In addition, since k-means clustering is an iterative process with no guaranteed convergence towards the global optimum, the algorithm is run $n_{rep}$ times.
In several previous works using CAPs, it was also suggested to solve Equation~\ref{CAPs_EQ1} after setting to 0 the voxel intensity values that, for each frame of interest, would not be part of the largest $P_P$ or $P_N$ percents (for positive-valued and negative-valued voxels, respectively)~\citep{Liu2013,Liu2013b}.
Table~\ref{CHAPTER5_TABLE0} summarises the different parameters that are defined for CAPs analysis, and also highlights the default values that we used in this work.
$\mathcal{S}$& Voxel set to use as seed & Right dorsolateral prefrontal cortex \\
Polarity & Sign of the seed excursions to consider & Activation \\
Seed combination & Whether all or at least one seed should be (de)active to retain a time point & n.a. \\
$T$& Threshold for frame selection &$1.5$\\
$T_{mot}$& Threshold of framewise displacement above which to scrub &$0.3$mm \\
$K$& Number of clusters to use & 16 \\
$n_{rep}$& Number of replicates of k-means &$50$\\
$P_P$& Percentage of positive-valued voxels to keep in each frame for clustering &$100$\\
$P_N$& Percentage of negative-valued voxels to keep in each frame for clustering &$100$\\
$\dist$& Distance measure used for clustering &$\corr$\\\hline\hline
\end{tabular}
\normalsize
\caption{\textbf{Parameters to define for co-activation pattern analysis.} The seed region was extracted from a TPN Independent Component map derived in~\citet{Shirer2012}. Investigated cluster number values $\{K\}$ were determined through consensus clustering, a subsampling-based assessment of clustering robustness~\citep{Monti2003}, and the final choice for the analyses was defined based on an exploratory assessment of the brain/behaviour correlation significance.}
Once all retained frames have been assigned to CAPs representatives, it becomes possible to construct, for each subject, an empirical transition probability matrix $A_s$ that summarises the likelihood to transit from a given CAP at time $t$ to another at time $t+1$. Another available piece of information regards the likelihood to transit from and back to the baseline state (when the seed was not significantly (de)active). Further, if separate subject populations are used in computing CAPs and deriving associated metrics (as in our example applications below), there are also occurrences to enter an extra state associated to frames that could not be matched to any CAP with sufficient certainty.
An indicative example of averaged transition probability matrix across subjects is displayed in Figure~\ref{CHAPTER5_P1F1}A (left column). Individual elements of the transition probability matrix may be considered as such~\citep{Chen2015}, which would amount to a total of $K^2$ values per subject.
To meaningfully lower the amount of features of interest, we propose to rather view the available information as a directional graph representation, from which a series of summarising metrics can be derived~\citep{Rubinov2010}. First, by sampling the diagonal elements of the matrix, we obtain a measure of \textit{resilience} for each CAP: that is, the likelihood to remain in the same configuration from time $t$ to $t+1$. Second, after having set the diagonal elements of the matrix to 0, we can define the \textit{in-degree}$k_{in}$ (how likely is a CAP visited from any other), the \textit{out-degree}$k_{out}$ (how likely is a CAP exited towards any other), and the \textit{betweenness centrality} (how important is a CAP regarding the shortest paths between other pairs)~\citep{Freeman1979}. In total, the feature space has thus been reduced from $K^2$ to $4K$. This alternative viewpoint is exemplified in Figure~\ref{CHAPTER5_P1F1}A (right column).\todo{About Figure 1 -- For the barplot: could you add a std across subjects? For the graph plot, add colorbar with encoding of green?}
\caption{\textbf{Generation of CAP metrics.}\textbf{(A)} Transitions across CAPs can be viewed in terms of individual transition probabilities for a total of $K^2$ features (left), or alternatively, a directional graph representation can be constructed (right) to extract in-degree (red bars), out-degree (blue bars), betweenness centrality (green color coding of the nodes in the bottom plot) and resilience (size of the nodes) information, for a reduced total of $4K$ features.\textbf{(B)} The extracted metrics contain equivalent information compared to the more traditional occurrences, and each such metric characterises partly different aspects of CAP dynamics, as seen by moderate correlations. Each box plot/violin plot representation depicts correlation values across $K$ CAPs. Illustrations from this figure are generated from the data presented in our example application, with $K=32$ and $T_P=5$\%.}
\label{CHAPTER5_P1F1}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In several works, \textit{counts} or \textit{occurrences} (that is, how many times a given CAP is expressed) were used as metrics of interest~\citep{DiPerri2018,Kaiser2019,Tuleasca2019}. We verified that our suggested metrics also include the information rendered by the counts: as seen in Figure~\ref{CHAPTER5_P1F1}B, the average correlation across CAPs between counts and in-degree, out-degree or resilience exceeded 0.8 (respectively $\rho=0.83\pm0.11$, $\rho=0.85\pm0.08$ and $\rho=0.81\pm0.07$). From pair-wise comparisons between our four metrics, it can also be seen that in-degree and out-degree are strongly correlated ($\rho=0.87\pm0.1$), while resilience and betweenness centrality capture separate information given their more moderate correlations (for resilience: $\rho=0.45\pm0.11$, $\rho=0.5\pm0.11$ and $\rho=0.3\pm0.19$ with in-degree, out-degree and betweenness centrality, respectively; for betweenness centrality: $\rho=0.59\pm0.13$ and $\rho=0.59\pm0.13$ with in-degree and out-degree, respectively). Despite their overall similarity, we decided to keep both in-degree and out-degree as they still yielded different values in specific CAP cases.
In addition to the above metrics that summarise the transitory behaviour across different CAPs, an interesting complement is the assessment of which CAPs are entered from the baseline state of seed activity, as well as of which CAPs are the ones expressed prior a return to baseline activity. With this additional information, a total of $6K$ features of interest per subject is available. These are the summarizing measures that we use in our example application.
\subsection{TbCAPs: implementation}
We implemented the CAPs processing pipeline as a toolbox in Matlab version 2017a (The MathWorks, Natick, USA). This software is freely accessible at \url{https://c4science.ch/source/CAP_Toolbox.git}. It contains a graphical user interface to facilitate the different steps of the pipeline. In addition, we also provide a scripted version of a typical analysis pipeline for power-users. An illustrative display of the graphical user interface at the end of a typical analysis is provided in Figure~\ref{CHAPTER5_P1F2}. Next, we describe the steps to be performed by the user, and the available options at each stage of the analysis:
\caption{\textbf{Illustration of TbCAPs graphical user interface.} The typical output of a CAPs analysis is displayed for a set of healthy volunteers, with three seeds chosen within a functional circuitry associated to essential tremor (this data was explicitly analysed in~\citet{Tuleasca2019}).}
\label{CHAPTER5_P1F2}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Data preparation}
The toolbox requires four different input files, to be selected at once when clicking on the \textbf{A. Load data} button: first, the preprocessed fMRI data should be provided as a Matlab structure of size $1\times S$. Each cell should contain a matrix of size $T \times V_{all}$, with $V_{all}$ all the acquired voxels. Second, framewise displacement over time for all subjects should be provided as a matrix of size $T \times S$. Third, a NIFTI file header containing the dimensionality information about the data at hand should be included. Fourth, a mask of size $V_{all}$ should be passed as a logical vector, setting to $1$ the $V$ voxels to include for CAPs analysis.
\subsubsection{Spatiotemporal selection}
Following data loading, the user is prompted to select one or more seeds to use in the analysis: all seed files should be entered at once, each as logical vectors of size $V \times1$. The brain area covered by the seeds can be inspected in brain slice representations, which can be navigated through by means of dedicated sliders. For now, we allow up to three separate seeds to be entered for the analyses. In addition, the interested user can also plot the seed-based correlation map associated to the first selected seed.
The next step is to select which types of events should be retained (activation versus deactivation), and if more than one seed was selected, whether all seeds should show an extreme event at once to select a time point (\textbf{Union} option), or if a frame should be kept as long as at least one does so (\textbf{Intersection} option).\todo{Same remark as above} The user is also prompted to determine the threshold $T$ to use in selecting frames (or alternatively, a percentage $P$ of most (de)active frames to retain), and the threshold $T_m$\todo{I would rather call this $\text{FD}_\text{limit}$ or something like that} above which frames will be deemed excessively corrupted by head motion, and scrubbed out.
\subsubsection{Generation of co-activation patterns}
Regarding the subsequent generation of CAPs, if the optimal number of clusters $K$ to select is not known a priori, we offer the possibility to run consensus clustering~\citep{Monti2003}, where clustering is run many times from $K=2$ to a user-specified $K_{max}$ using a subsample of the data (the percentage of data points to use is specified by $P$, and the number of iterations by $N$). A good clustering solution is one for which across folds, two frames are either always clustered together, or never clustered together (but not an intermediate case). We quantify this by the Percentage of Ambiguously Clustered pairs (PAC)~\citep{Senbabaoglu2014}, and display the stability measure $1-PAC$.
When k-means clustering is run by clicking the \textbf{Cluster} button, the first $5$ CAPs with most occurrences across the subject population are displayed, and can be visually inspected. If using the \textit{Intersection} option in a multi-seed analysis, the user is also shown pie charts summarising, for each CAP, what fraction of frames was selected in a given seed combination configuration.
\subsubsection{Metrics generation and frame assignment}
Finally, upon clicking the \textbf{Compute metrics} button, displays of CAP expression time courses, cumulative CAP expression along time, and violin plots summarising the distribution of computed metrics across subjects for each CAP, are provided. Transition probability matrices can also be inspected for each subject.
In some settings, the user may wish to compare different subject populations: this can be done by sequentially loading up to $4$ different populations at the start of the analysis. CAPs will be derived from the first population, and there is then the option to assign the frames from the other populations to the CAPs by a matching process. In doing so, the spatial correlation between a frame and the CAP to which it is most similar is compared to the distribution of spatial correlations of the frames from population 1 that belong to the CAP in question: if the $T_P$\textsuperscript{th} percentile of this distribution is exceeded, assignment is performed; else, the frame is left unassigned and belongs to an extra $(K+1)$\textsuperscript{th} cluster.
\subsection{Application to experimental fMRI data}
\subsubsection{Functional data preprocessing}
As a proof of feasibility and application of TbCAPs, we considered a sample of $181$ subjects from the \textit{Human Connectome Project}~\citep{VanEssen2013}, selected because they had at least one fully exploitable resting-state scanning session on which to apply the method, and less than 5\% of recorded behavioural entries that were missing. Details regarding acquisition parameters can be found elsewhere~\citep{Smith2013}, but briefly, the data was acquired at a TR of 0.72s over 15 minutes (for a total of 1200 fMRI volumes), with a spatial resolution following initial preprocessing steps of 2 x 2 x 2mm\textsuperscript{3}.
We started from the \textit{minimally preprocessed} resting-state data (first session, LR acquisition direction). The first 10 samples of the data were discarded. We then performed linear detrending, and regressed out low-frequency components of the discrete cosine transform basis with a cutoff frequency at 0.01Hz. Due to collinearity with this basis, we did not regress out average white matter or cerebrospinal fluid time courses. We also chose not to regress motion parameter time courses, as motion is handled within the co-activation pattern pipeline by scrubbing, and because recent evidence points at the fact that motion regression schemes may not always be beneficial in the context of brain/behaviour analyses~\citep{Bolton2019c}. As for global signal regression, given the lack of a clear consensus~\citep{Murphy2017}, we preferred to leave the data as untouched as possible and did not include it.
Following the regression step, the data was scrubbed at a framewise displacement threshold of $0.3$mm, and excised volumes were estimated with cubic spline interpolation. Although scrubbing is performed within TbCAPs, we reasoned that if we wished to try and assign scrubbed frames to CAPs in our additional analyses regarding head motion, it would make more sense to have previously corrected these volumes to the best of our abilities.
Then, individual fMRI volumes were smoothed at a full-width at half maximum value of 5mm, and in order to make the analyses computationally more affordable, spatial resolution was downsampled at 3 $\times$ 3 $\times$ 3mm\textsuperscript{3}. Finally, each voxel-wise time course was z-scored to obtain the resting-state functional imaging input data to our pipeline.
\subsubsection{Selection of seed and behaviour of interest}
As a behaviour of interest to study, we selected the Short Penn Continuous Performance Test (SCPT), which quantifies continuous sustained attention~\citep{Gur2010}. In more details, participants see vertical and horizontal red lines flash on screen, and from block to block, must respond either when the lines form a number, or a letter. The lines are displayed for 300ms, followed by a 700ms inter-trial interval.
We started from raw behavioural entries provided by the HCP, for $951$ different subjects. There are $8$ available SCPT measures: amount of true positives, false positives, true negatives or false negatives, median response time for true positive responses, sensitivity, specificity and longest run of non-responses. In order to reduce this information into one summary measure while filling in missing behavioural entries, we performed probabilistic PCA~\citep{Bishop1999}. The output composite score positively correlated with true positives, true negatives, sensitivity and specificity ($\rho=$0.24, 1.00, 0.25 and 1.00, respectively), thus summarising overall task performance. We z-scored this output measure across subjects, in order to quantify performance with respect to the overall population. We then extracted the behavioural data related to the $181$ subjects considered in this work.
To study sustained attention, we focussed on a right dorsolateral prefrontal cortex seed from the task-positive network, which we extracted from the associated Independent Component map provided by~\citet{Shirer2012}. Our hypothesis was that the expression of different TPN configurations would relate to sustained attention performance.
In both analyses, we resorted to a threshold $T=1.5$ to select active frames. In the multi-seed case, we retained the time points for which at least one of the seeds was active. Further, we performed scrubbing with a framewise displacement threshold of $0.3$mm.
To avoid double dipping~\citep{Kriegeskorte2009}, CAPs were extracted from a randomly selected subset of $100$ subjects, while we performed correlations with behaviour for the remaining $81$ only. To determine the optimal number of clusters, we used consensus clustering~\citep{Monti2003}. We then ran k-means $n_{rep}=50$ separate times, keeping the best solution. We included all voxels in the analyses ($P_{P}=P_{N}=100$\%), and used spatial correlation as our distance measure; given two similarly-sized vectors $\mathbf{x}$ and $\mathbf{y}$, this thus yields $\dist(\mathbf{x},\mathbf{y})=1-\corr(\mathbf{x},\mathbf{y})$.
Following the extraction of the CAPs on our $100$\textit{training subjects}, we determined which CAP was expressed at each retained fMRI volumes of the other $81$ subjects. To do so, we used the aforementioned assignment process with $T_{P}$ ranging from 0 to 100\%.
\subsubsection{Assessment of brain/behaviour relationships}
As imaging metrics of interest, we considered in-degree, betweenness centrality and resilience for each CAP, and also included the amount of excursions from the baseline state, and the amount of excursions back to the baseline state. Thus, we generated a total of $5K$ imaging features per subject.
After having obtained the behavioral scores $\mathbf{b}\in\mathbb{R}^{81\times1}$ and metrics ${M}\in\mathbb{R}^{81\times6K}$ for our population of subjects, we used Partial Least Squares (PLS) analysis~\citep{McIntosh2004,Krishnan2011} to probe the existence of a brain/behaviour relationship.
Briefly, consider a matrix of behavioural features $B\in\mathbb{R}^{S\times n_B}$ and a matrix of imaging metrics $M\in\mathbb{R}^{S\times n_M}$. Assuming that $n_B < n_M$, and using the singular value decomposition, the covariance between these two sets is given by:
\begin{equation*}
R = M^{\top} B = U \Sigma V^{\top} = \sum_{i=1}^{n_{B}}\sigma_i \mathbf{u_{i}}\mathbf{v_{i}}^{\top},
\end{equation*}
where each column in $U$ and $V$ contains the weights (so called \textit{saliences}) that respectively multiply imaging and behavioural markers to yield a maximised covariance between both sets. The associated singular value $\sigma_i$ is proportional to the fraction of covariance explained by the component at hand.
In our case, since $n_{B}=1$ (we only consider one behavioural measure), only one covariance component is retrieved, which implies $v_{1}=1$. The interesting information lies in $\mathbf{u_{1}}$: positive-valued saliences highlight metrics that are larger in subjects who show a greater cognitive ability, and negative-valued saliences are associated to metrics that, when larger, impede attentional performance.
Prior to running the algorithm, each of the $6K$ features was z-scored across subjects. In order to assess significance, we reran PLS $1000$ times after having randomly shuffled the subject entries in one of the two matrices; to non-parametrically derive a p-value, the singular value of the actual covariance component was compared to the null distribution constructed from this permutation process.
Further, to establish the stability of the salience weights, we reran PLS $1000$ times using a subsample of 80\% of the data each time, and computed a bootstrap score for each salience weight as its mean across folds divided by its standard deviation.
\subsubsection{Influence of head motion on quantified metrics}
While scrubbing enables to minimise the deleterious impacts of motion on the analysis and compute clean CAPs, discarding frames also has the potential to distort transition probability estimates. For example, a succession of three frames in the same state (which would amount to a higher resilience for the CAP in question) would not be captured if the middle frame is scrubbed out.
To verify that our findings were minimally sensitive to this effect, we ran another series of analyses in which we also performed the aforementioned assignment process on scrubbed frames, with a similar $T_{P}$ range as for assigning test subject frames. This way, frames strongly distorted by head motion still do not enter the analysis, but more mildly affected fMRI volumes can still be matched to their CAP. We verified the reproducibility of our findings upon this additional analytical step.
%%%% Application of CAPs analysis on the HCP data
\section{Results}
Consensus clustering results are displayed in Figure~\ref{CHAPTER5_P1F3}A for $K$ values ranging from 10 to 40. Positive peaks highlight good candidate values (see figure legend for details); such values are present for diverse numbers of clusters (more notably at $K=16,22,32$). While a lower number of clusters yields a reduced feature space and more interpretable outcomes, CAPs may not be segmented finely enough to resolve insightful dynamic properties regarding cognition. Our strategy was thus to first perform an exploratory assessment, in which we evaluated the significance of the brain/behaviour correlation across a set of candidate $K$ values ($K_{opt}=\{14,16,22,32\}$) and assignment thresholds (forcing the assignment of all frames, or using $T_P=[0:5:100]$, to select the relevant parameters to proceed forward with.
The results of this exploratory assessment are displayed in Figure~\ref{CHAPTER5_P1F3}B, when scrubbed frames are discarded (left panel) or also assigned to the CAPs at threshold $T_P$. Both settings yield very similar significance values, which is good evidence that remaining head motion effects only have a minor influence on the analyses. As the assignment threshold increases (that is, less and less frames are assigned because the criterion becomes more and more stringent), significance generally decreases. A smooth spot can be observed for $K=\{16,22\}$ and $T_P < 15$\%, which indicates that this granularity is optimal in the context of behavioural prediction. We selected $K=16$ and $T_P=0$\% as values for more detailed subsequent analyses.
\caption{\textbf{Parameter selection.}\textbf{(A)} Consensus clustering results for a range of K values from 10 to 40. Outputs from the algorithm are percentages of ambiguously clustered pairs (PAC); since this measure requires the definition of an interval of consensus values within which clustering is deemed \textit{ambiguous}, we investigated a range of values as colour coded from black to yellow. For each value, we fitted a decreasing exponential to capture the overall tendency of PAC values going down with larger cluster number. The y-axis of the plot depicts the difference between this fitted value and the actual value: thus, a positive difference means that the considered cluster number is more satisfying than what would be predicted in terms of the overall behaviour. \textbf{(B)} Across candidate cluster number values selected from consensus clustering, and assignment threshold values, significance of the relationship between CAP metrics and attentional abilities, as quantified with the p-value obtained by PLS analysis. The left panel depicts the results for which scrubbed frames are not considered at all, while in the right panel, scrubbed frames were also assigned to the CAPs. The infinity symbol is used to depict a case in which assignment is done for all frames, even if a frame is not sufficiently close to any CAP when comparing its spatial correlation to the associated correlation distribution.}
\label{CHAPTER5_P1F3}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CAPs are displayed, for this chosen parameter set, in Figure~\ref{CHAPTER5_P1F4}A, while their involvement in driving the brain/behaviour relationship, as quantified by salience weights across our range of investigated metrics, is depicted in Figure~\ref{CHAPTER5_P1F4}B. The correlation between actual and predicted attentional performance was strongly significant (R=0.587, p<0.001; Figure~\ref{CHAPTER5_P1F4}C). The associated covariance component found by PLS analysis was significant at p=0.003. Note that this relationship is derived from only subjects that were not used to construct the CAPs.
CAP\textsubscript{1} depicts co-activation of a range of resting-state networks, including the auditory, somatomotor, visual and salience ones. Attentional performance was better in the subjects that transited more frequently from the baseline state of seed activity to this CAP. CAP\textsubscript{1} was also more often the entry point towards other CAPs in high performance subjects, as seen from a strongly positive out-degree salience weight.
Good subjects in terms of continuous performance also more often entered CAP\textsubscript{2} and CAP\textsubscript{7} from other states (as seen from positive in-degree salience weights), and these same 2 CAPs were also more influential in the transitory behaviour of CAP dynamics (since betweenness centrality salience weights also showed large positive values). In both CAPs, the seed region co-activates with a restricted set of areas including the right inferior parietal cortex (for both), the posterior cingulate cortex and medial prefrontal cortex (for CAP\textsubscript{2}), and the right anterior prefrontal cortex (for CAP\textsubscript{7}).
CAP\textsubscript{3} and CAP\textsubscript{5} were associated to good attentional abilities from the viewpoint of several metrics, which emphasises the importance of their expression: for CAP\textsubscript{3}, it involved resilience, in-degree and out-degree, while for CAP\textsubscript{5}, it was return to baseline, resilience, in-degree and betweenness centrality. Both CAPs include strong co-activation with the right inferior parietal cortex, and for CAP\textsubscript{3}, also with the left cerebellum lobule VI and a subpart of the occipital cortex.
CAP\textsubscript{4} and CAP\textsubscript{14} were the only states whose expression was detrimental for attentional performance, in terms of betweenness centrality for the former, and of return to baseline for the latter. CAP\textsubscript{4} displayed bilateral right superior cortex and anterior prefrontal cortex co-activation with the seed, while for CAP\textsubscript{14}, involved areas were the fusiform gyrus, parahippocampal cortex, and a diffuse right lateralised spot covering parts of the auditory, secondary somatosensory and posterior insular cortices.
The majority of the other CAPs that did not show any link to attentional performance involved co-activations with regions that were not part of the attentional networks: for instance, CAP\textsubscript{6} includes the anterior cingulate cortex and anterior insula; CAP\textsubscript{8} contains the anterior cingulate, visual and right somatosensory cortices; CAP\textsubscript{9} showcases primary visual, and auditory cortices; CAP\textsubscript{10} shows the angular gyrus and part of the precuneus; CAP\textsubscript{11} and CAP\textsubscript{15} mostly highlight ventral medial prefrontal cortex signal.
\caption{\textbf{Co-activation patterns, and relationship to attentional abilities.}\textbf{(A)} The 16 obtained CAPs are plotted, with coloured rectangles symbolising the CAPs that are particularly important (bootstrap score of the associated salience weight larger than 2) for the brain/behaviour relationship as quantified from a given metric (grey: entries from baseline state, black: exits to baseline state, yellow: resilience, red: in-degree, blue: out-degree, green: betweenness centrality). Filled rectangles highlight beneficial CAPs, while hollow rectangles depict detrimental CAPs (negative bootstrap score). \textbf{(B)} Salience weights across all 16 CAPs and the 6 investigated CAP metrics. \textbf{(C)} Actual attentional performance score (x-axis) versus predicted values with PLS (y-axis).}
\label{CHAPTER5_P1F4}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Discussion
\section{Discussion}
\label{sec:discussion}
In this work, we have introduced TbCAPs, a toolbox for co-activation patterns analysis, which provides practitioners with an intuitive dedicated graphical user interface as well as a powerful scripting equivalent. It provides an easy control over all key analytical parameters of the technique, novel methodological additions for augmented analyses, and facilitated visualisation of the resulting CAPs and associated metrics. Although we have focussed on the usefulness of CAPs analysis in the resting-state setting, we also remark that nothing precludes the use of the technique in task-based investigations.
As most CAP studies to date have revealed the potential of the approach in clinical settings~\citep{Amico2014,DiPerri2018,Kaiser2019,Tuleasca2019}, we sought to demonstrate the relevance of the technique in another context; i.e., rather than considering a \textit{classification} problem in which two or more distinct subject populations are separated, we considered a \textit{regression} task in which we attempted to explain attentional abilities within a more homogeneous population in a continuous vigilance task by means of CAPs dynamics.
We observed that of all the extracted CAPs showing co-activation\todo{This is confusing due to double-use of term. I would name the ``second-order'' effect COUPLING, similar to what we use in the iCAPs context.} with the right dorsolateral prefrontal seed, the large majority either did not appear to be involved in attentional abilities, or showed positive salience weights indicating a positive impact of their expression. This is not so surprising, given that our analysis was focused on a region of the attention network in the first place. The common feature of beneficial CAPs appeared to be the co-activation of an array of other regions previously pinpointed in continuous performance tasks, including the inferior parietal cortex, cerebellum lobule VI or occipital cortex~\citep{Hager1998,Ogg2008,Tana2010}. At the same time, these beneficial CAPs also barely involved co-activation of other functionally distinct networks.
The one CAP for which the above reasoning does not hold is CAP\textsubscript{1}: despite the involvement of a very diverse set of regions, it was also retrieved as beneficial for attentional performance. More precisely, contrarily to most of the others, the expression of this CAP appears to be essential at the start of a seed activation sequence: indeed, salience weights were large specifically for the entries from baseline, and out-degree metrics. In other words, there is first a transition from baseline to this CAP, followed by the exit of that configuration to reach more spatially well-defined states. This involvement of short-lived periods of extensive cross-network interactions in mediating some aspects of human cognition has recently started to be appreciated as an insightful functional brain mechanism~\citep{Betzel2016,Fukushima2018b}.
The fact that all the probed metrics significantly contributed to explain attention is good evidence in favour of the temporal complexity of functional brain dynamics: instead of an instantaneous characterisation or a one-frame expression of a telling functional state, what truly matters is a complex mix between how activation starts (captured by the from-baseline and to-baseline metrics), how transitions occur across distinct functional states (as seen from in-degree, out-degree and betweenness centrality), and how lasting a given state is (as quantified by resilience). Our characterisation can be placed in the broad family of temporal modelling approaches, of which other notable examples include the use of hidden Markov models~\citep{Vidaurre2017,Bolton2017b} or graph-theoretical analysis for energy landscape~\citep{Kang2019}.
Prediction of continuous performance abilities from resting-state fMRI recordings has been shown possible in previous functional connectivity work relying on second-order correlational measures across brain regions~\citep{Rosenberg2016}. More recently, this characterisation has been pushed to the dynamic level by~\citet{Fong2019}, who showed that prediction can also be successfully achieved when \textit{temporal variability}, which quantifies fluctuations in functional connectivity over the course of a scanning session, is used as a metric of interest. Our prediction accuracy is on par with the one achieved in this whole-brain analysis, despite focussing on one seed region. Interestingly, the authors described the fact that lowered temporal variability was beneficial for better attentional performance, and that many of the most important features for prediction involved the executive-control brain networks: this is fully consistent with our results, in which increased resilience (which can be expected to yield lowered temporal variability) of CAPs featuring executive-control areas is beneficial.
In comparison to clinical applications of CAPs analysis, in which between $3$ to $8$ CAPs are typically considered, a finer granularity was required in the present work (see Figure~\ref{CHAPTER5_P1F3}). This is not surprising given that the regression problem at hand here is more challenging than a classification task, as we need to predict a value within a continuum. Furthermore, the functional underpinnings of inter-individual differences in cognitive abilities are likely more subtle than when comparing subjects across consciousness levels or disease severity levels. In fact, if a too low number of CAPs is extracted, patterns with a different cognitive relevance are averaged together as a single cluster, which impedes prediction; conversely, if a too large number of CAPs is extracted, meaningful configurations become further segmented, and statistical power is lost due to the smaller amount of frames constituting each CAP.
In future work, it will be interesting to examine clinical or cognitive research hypotheses at the broader focus level of more than one seed region. As alluded to above, this is already feasible with our current toolbox version, and may enable to better bridge the gap between seed-based and whole-brain analyses. We also foresee additional technical developments in the near future, such as the possibility to extract \textit{co-activation sequences} (that is, series of successive fMRI volumes) rather than CAPs.
\todo[inline]{About Figure 4. For the B panel. I would add dashed horizontal lines at the positive/negative threshold of the bootstrap score when you consider the variable as robust. Maybe, you could also indicate those that survive this threshold by making the box filled with color, vs fill in another way when they don't. }