# Spatio-temporal model for dynamic functional connectivity in resting state fMRI analysis

Massyl Moudoud

ICube, UMR 7357, CNRS,
University of Strasbourg
Strasbourg, France
mmoudoud@unistra.fr

Céline Meillier

ICube, UMR 7357, CNRS,
University of Strasbourg
Strasbourg, France
meillier@unistra.fr

Marion Sourty

ICube, UMR 7357, CNRS,
University of Strasbourg
Strasbourg, France
sourty@unistra.fr

Vincent Mazet

ICube, UMR 7357, CNRS,
University of Strasbourg
Strasbourg, France
vincent.mazet@unistra.fr

**Abstract**—Functional connectivity (FC) in fMRI exhibits spatio-temporal dynamics. A new fine-grained, single subject modeling of this dynamic FC is introduced based on functional connectivity units (FCUs), defined as a small set of anatomical brain regions. The estimation of FCUs is considered as a dictionary learning problem which is solved with an ADMM algorithm. A complex synthetic data generation procedure is developed to mimic brain dynamics and proved the performance of the proposed algorithm. The analyses of resting state fMRI signals of mice revealed that both the spatial structure and the temporal activity of the obtained FCUs are biologically valid and show high inter-subject coherence.

**Index Terms**—Resting state fMRI, dictionary learning, dynamic functional connectivity

## I. INTRODUCTION

Resting-state functional magnetic resonance imaging (rs-fMRI) allows studying brain functional connectivity (FC) at rest. The objective is to estimate brain networks, defined as 3D spatial structures in brain volume, along with their temporal signature but also the relations between these networks. The interactions are both temporal with synchronous activity of several networks over a given time interval and spatial with networks merging (totally or partially) or dividing into sub-networks.

Mathematically, this reduces to expressing the fMRI data as a matrix $Y \in \mathbb{R}^{T \times V}$ where $T$ is the number of time samples and $V$ is the number of voxels:

$$Y = US \quad (1)$$

where $S \in \mathbb{R}^{R \times V}$ is the matrix containing the spatial maps of the $R$ functional networks in its rows and $U \in \mathbb{R}^{T \times R}$ contains in each column the temporal signature of a network. In the literature, the joint estimation of $U$ and $S$ is performed through spatial independent component analysis [1], dictionary learning [2], and more recently anatomical segmentation of the brain (atlas) to define functional networks or regions of interest (ROIs). The correlation matrix of the temporal signatures of these ROIs gives a stationary representation of the FC.

In the last decade, the advent of methods to analyze the dynamic FC (dFC) challenges this stationary definition (see Fig. 1 in [3]). In those methods, the temporal

dynamics of the inter-network connectivity modeled by the matrix $U$ is studied. The timecourses of the $R$ networks are analyzed through a sliding window of size $\Delta_w$ on which the correlation matrices between the networks is computed. A K-means algorithm is then applied on the matrices to identify different configurations of the connectivity at the whole-brain level for a group of subjects [1]. The centroid correlation matrix of each identified class constitutes what is called a brain state or connectivity state, it is then possible to study the transitions between states in time for each subject. This is a first step towards a modeling of the temporal dynamics of brain connectivity. Nevertheless, the spatial definition of brain networks stays static since the model (1) used in the different studies is based on the stationarity assumption of the brain activity at rest. This assumption is questioned in several recent studies (see [3] for a review). Moreover, by supposing that the functional networks do not change in time, the estimated connectivity states are defined at the scale of the whole brain, while it can be reasonably assumed that all networks do not necessarily exhibit the same temporal dynamics.

The goal of this paper is to introduce a spatio-temporal model of the dFC and validate its interest. An algorithm to express fMRI time series in the proposed model is derived based on dictionary learning. The proposed model is evaluated on a synthetic dataset created for this purpose and is challenged on real rs-fMRI mice data.

## II. REPRESENTING SPATIO-TEMPORAL DYNAMICS

To overcome the limitations of defining brain states, we define functional connectivity units (FCUs) at the scale of (sub)-networks based on our previous exploratory works [4], [5]. An FCU is composed of few anatomical ROIs that are strongly connected together at given time intervals. An FCU can represent a network known to be related to some brain activity or can be defined at a smaller scale with at least 3 ROIs. This definition of the FCU allows for great flexibility in the brain dFC. A pair of ROIs can be part of multiple FCUs that can be active at the same time, in which case the correlation of the pair of ROIs is the sum of its activations within the concurrently active FCUs weighted by the activation of the FCUs. The model is defined as:

$$C = DA \quad (2)$$

This work has been supported by the ANR project DynaSTI: ANR-22-CE45-0008
where $\mathbf{C} = [\mathbf{c}_1, \dots, \mathbf{c}_T] \in \mathbb{R}^{E \times T}$ and $\mathbf{c}_t \in \mathbb{R}^{E \times 1}$ is the vectorized lower triangular part of the correlation matrices of the $R$ timecourses of $\mathbf{U}$ at each time window ($E = R(R-1)/2$). Matrix $\mathbf{D} \in \mathbb{R}^{E \times P}$ is the dictionary of the spatial representation of the FCUs. Each column of $\mathbf{D}$ contains the pairs of ROIs that constitute an FCU with their relative connectivity within the FCU. The activation matrix $\mathbf{A} \in \mathbb{R}^{P \times T}$ holds the temporal pattern of activation of each FCU.

### III. DICTIONARY LEARNING

#### A. Optimization problem

The space of possible candidates for $\mathbf{D}$ and $\mathbf{A}$ is very large. In order to restrict the space to biologically valid solutions, a set of constraints are imposed on the matrices $\mathbf{D}$ and $\mathbf{A}$ resulting in the following optimization problem:

$$ \min_{\mathbf{D},\mathbf{A}} \frac{1}{2} \|\mathbf{C} - \mathbf{D}\mathbf{A}\|_F^2 + \lambda_{\ell_1} \|\mathbf{A}\|_1 + \lambda_{TV} TV(\mathbf{A}) \\ + \mathcal{I}_{\mathbb{R}+}(\mathbf{A}) + \mathcal{I}_{\tilde{\mathbb{D}}}(\mathbf{D}) + \mathcal{I}_{\tilde{\mathbb{R}}+}(\mathbf{D}), \quad (3) $$

where $\mathcal{I}_C(x)$ is the indicator function on the set $C$: it takes the value $+\infty$ if $x \notin C$ and zero otherwise. In continuation with our works [4], [5], we focus on significant positive correlations, thus the positivity of all entries of both $\mathbf{D}$ and $\mathbf{A}$ is enforced. The term $TV(\mathbf{A})$ is the total variation regularization defined by

$$ TV(\mathbf{A}) = \sum_{t=2}^{T} \|\mathbf{a}_t - \mathbf{a}_{t-1}\|_1. \quad (4) $$

It promotes smoother solutions and reduces spurious fluctuation in $\mathbf{A}$. The matrix $\tilde{\mathbf{D}}$ is a binary matrix defining the spatial support of the FCUs, i.e. the pairs of ROIs that compose an FCU. The term $\mathcal{I}_{\tilde{\mathbb{D}}}(\mathbf{D})$ imposes matrix $\mathbf{D}$ to have the same structure as $\tilde{\mathbf{D}}$: if $\tilde{\mathbf{D}}_{ij} = 0$ then $\mathbf{D}_{ij} = 0$. The matrix $\tilde{\mathbf{D}}$ is either provided by experts and contains networks and sub-networks relevant to the objective of the study (that are believed to take part in or be influenced by the pathology for example), or it can be estimated by a thresholding of the correlation matrices to detect the significant correlations that occur in group, involving a hypothesis testing scheme to control the threshold level [4].

#### B. Joint optimization

The problem (3) is jointly non-convex over $\mathbf{D}$ and $\mathbf{A}$ and can be seen as a non-negative matrix factorization (NMF) problem with TV regularization for $\mathbf{A}$ and constrained spatial support for $\mathbf{D}$. Albeit the attractiveness of NMF solutions, an alternating minimization of optimization sub-problems according to $\mathbf{A}$ and $\mathbf{D}$ is used to solve it as described in Algorithm 1. While it is probably less efficient than NMF, we propose a more versatile solution allowing to relax the positivity constraint on $\mathbf{D}$ in future work.

#### C. Estimation of the activation matrix $\mathbf{A}$

If $\mathbf{D}$ is fixed, then the estimation of $\mathbf{A}$ comes down to:

$$ \min_{\mathbf{A}} \frac{1}{2} \|\mathbf{C} - \mathbf{D}\mathbf{A}\|_F^2 + \lambda_{\ell_1} \|\mathbf{A}\|_1 + \lambda_{TV} TV(\mathbf{A}) + \mathcal{I}_{\mathbb{R}+}(\mathbf{A}). \quad (5) $$

### Algorithm 1 Alternating minimization of $\mathbf{A}$ and $\mathbf{D}$

Initialize $\mathbf{D}^0 \leftarrow \tilde{\mathbf{D}}$ $i \leftarrow 0$

repeat

Find $\mathbf{A}^{i+1}$ with Sunsal-TV (Section III-C)

Find $\mathbf{D}^{i+1}$ with ADMM (Algorithm 2)

$\mathbf{C}^{i+1} \leftarrow \mathbf{D}^{i+1} \mathbf{A}^{i+1}$

$i \leftarrow i + 1$

until $\|\mathbf{C}^i - \mathbf{C}^{i-1}\|_F < \epsilon$

return $\mathbf{D}^i$ and $\mathbf{A}^i$

The problem (5) can be solved with Sunsal-TV algorithm [6] by adapting the total variation regularization to operate only on the rows of $\mathbf{A}$ as defined in (4) while the original Sunsal-TV algorithm applied TV regularization on both rows and columns of $\mathbf{A}$. The algorithm is initialized using the regularized pseudoinverse of $\tilde{\mathbf{D}}$. In the subsequent calls, the optimization starts with the result of the last call (warm start).

#### D. Estimation of the dictionary $\mathbf{D}$

If $\mathbf{A}$ is fixed, then the estimation of $\mathbf{D}$ comes down to:

$$ \min_{\mathbf{D}} \frac{1}{2} \|\mathbf{C} - \mathbf{D}\mathbf{A}\|_F^2 + \mathcal{I}_{\tilde{\mathbf{D}}}(\mathbf{D}) + \mathcal{I}_{\tilde{\mathbb{R}}+}(\mathbf{D}). \quad (6) $$

The data fidelity term can be recast in terms of the columns of $\mathbf{A}$ and $\mathbf{C}$ as $\frac{1}{2} \sum_{t=1}^{T} \|\mathbf{c}_t - \mathbf{D}\mathbf{a}_t\|_2^2$. Then (6) is reformulated as a global variable consensus problem [7, section 7]:

$$ \begin{array}{ll} \min & \displaystyle \sum_{t=1}^{T} f_t(\mathbf{D}_t) + g(\mathbf{\Delta}) \\ \text{subject to} & \mathbf{D}_t = \mathbf{\Delta} \quad \forall t = 1 \dots T \\ \text{with} & f_t(\mathbf{D}_t) = \frac{1}{2} \|\mathbf{c}_t - \mathbf{D}_t \mathbf{a}_t\|_2^2 \\ \text{and} & g(\mathbf{\Delta}) = \mathcal{I}_{\tilde{\mathbf{D}}}(\mathbf{\Delta}) + \mathcal{I}_{\tilde{\mathbb{R}}+}(\mathbf{\Delta}). \end{array} \quad (7) $$

The alternating direction method of multipliers (ADMM) is well suited for this problem, it results in the formulation

$$ \mathbf{D}_t^{k+1} = \underset{\mathbf{D}_t}{\operatorname{argmin}} \left( f_t(\mathbf{D}_t) + \frac{\rho}{2} \|\mathbf{D}_t - \mathbf{\Delta}^k + \mathbf{U}_t^k\|_F^2 \right) \quad (8a) $$

$$ \mathbf{\Delta}^{k+1} = \underset{\mathbf{\Delta}}{\operatorname{argmin}} \left( \mathcal{I}_{\tilde{\mathbf{D}}}(\mathbf{\Delta}) + \mathcal{I}_{\tilde{\mathbb{R}}+}(\mathbf{\Delta}) + \frac{T\rho}{2} \|\mathbf{\Delta} - \tilde{\mathbf{D}}^{k+1} - \tilde{\mathbf{U}}^k\|_F^2 \right) \quad (8b) $$

$$ \mathbf{U}_t^{k+1} = \mathbf{U}_t^k + \mathbf{D}_t^{k+1} - \mathbf{\Delta}^{k+1}. \quad (8c) $$

The steps (8a) and (8c) are performed $\forall t = 1 \dots T$. The step size $\rho$ was empirically set to 1 for all the experiments. Matrix $\tilde{\mathbf{D}}$ denotes the average of the matrices $\mathbf{D}_t$, $\forall t = 1 \dots T$. The estimation of $\mathbf{D}$ is given in Algorithm 2.

### IV. SYNTHETIC DATA GENERATION

The proposed method is validated on synthetic data. Most of the synthetic datasets proposed in the literature only include time dynamics [3]. They are generated by creating a set of discrete “brain connectivity states” (correlation matrices) and create $\mathbf{C}$ by alternating the activation of the states. We propose to go further and include a spatial dynamic in the generated
**Algorithm 2 ADMM for D**

Initialize $D_t^0 \leftarrow 0$ and $U_t^0 \leftarrow 0 \quad \forall t = 1...T$

Initialize $\Delta_t^0 \leftarrow D_t^i$ from Algorithm 1

$k \leftarrow 0$

repeat

for $t = 1$ to $T$ do

$$D_t^{k+1} \leftarrow [c_t a_t^T + \rho (\Delta_t^k - U_t^k)] (a_t a_t^T + \rho I_P)^{-1}$$

end for

$$\Delta_t^{k+1} \leftarrow \mathcal{P}_{\bar{D}} \circ \mathcal{P}_{\mathbb{R}+}(\bar{D}_t^{k+1} + \bar{U}_t^k)$$

for $t = 1$ to $T$ do

$$U_t^{k+1} \leftarrow U_t^k + D_t^{k+1} - \Delta_t^{k+1}$$

end for

$k \leftarrow k + 1$

$$\mathbf{until} \quad ||\Delta_t^{k+1} - \Delta_t^k||_F < \delta$$

$\mathbf{return} \quad \Delta^k$

<figure><img src="image_1.png" /><figcaption>Fig. 1: Synthetic data generated from $C = DA$ with $R = 10$ ROIs ($E = 45$ pairs) and $P = 11$ FCUs. Time varying correlation $C$ (left) is the product of the dictionary of FCUs $D$ (center) and the time activation of the FCUs $A$ (right).</figcaption></figure>

data to better match the spatio-temporal dynamic assumption of the dFC. Both a dictionary of FCUs $D$ and an activation matrix $A$ are generated, the ground truth correlation matrix is the product $C = DA$.

### A. The dictionary D

We set $R=10$ ROIs, i.e. $E=45$ pairs of ROIs. Three main networks and some of their possible sub-networks are created to represent $P=11$ FCUs. Each column of the dictionary $D$ contains the weights of the pairs of ROIs that compose an FCU. The center panel of Fig. 1 shows the generated $D$.

### B. The activation matrix A

We generate the support $\tilde{A}$ of $A$ and the amplitudes $\tilde{A}$ separately. The support $\tilde{A}$, with $T = 1000$ time points, is a binary matrix indicating whether a given FCU is active at a given time. The initial state of each FCU is given by a binomial distribution. The state then alternates from active to inactive; its duration follows a uniform distribution $U(40, 150)$. A post-processing guarantees that if and only if a network is active then all of its sub-networks are inactive. The amplitudes $\tilde{A}$ are sampled from a normal distribution $N(1, 0.5)$ which are low-pass filtered along the rows (time), and finally masked by the elementwise product with the support $\tilde{A}$.

The resulting matrix $\tilde{D}\tilde{A}$ may contain correlations greater than 1. Moreover, due to the acquisition noise, pairs of ROIs can not be perfectly correlated in practice. For these reasons, in the simulations, the correlations are limited to be less than 0.9. This is done by optimizing the following problem using sequential least squares programming:

$$A = \underset{A}{\operatorname{argmin}} \|\tilde{A} - A\|_F^2 \quad \text{s.t.} \quad DA \le 0.9 \text{ and } A \ge 0. \quad (9)$$

An example of generated data is given in Fig. 1.

## V. RESULTS AND DISCUSSION

### A. Simulation

The proposed model is evaluated on the synthetic dataset presented in Fig. 1. Previous work [5] showed that the choice

of the length $\Delta_w$ of the sliding window on which correlations are computed has a great impact on the estimation of the true $C$. To avoid this influence and be able to evaluate the proposed method, it is applied directly on the generated matrix $C$. To simulate a more realistic situation, Gaussian noise is added to $C$ with signal to noise ratio (SNR) of 0 dB, $C$ being clipped to the interval $[-1, 1]$. We observe empirically that the optimization converges in tens of iterations.

The effect of the regularization parameters $\lambda_{\ell_1}$ and $\lambda_{TV}$ on the reconstruction error is evaluated by running the model for 2500 combinations of the two parameters (50 values for each) with $\lambda_{\ell_1} \in [0, 0.8]$ and $\lambda_{TV} \in [0, 2]$. Fig. 2 shows the Frobenius norm of the reconstruction error for the matrices $A$, $D$ and $C$. TV regularization has little effect on $C$ and $D$ and moderate effect on $A$, while most of the variation in the error is due to the $\ell_1$ regularization. It should be noted that the best regularization parameters will depend on the structure of the data (sparsity level and smoothness). The best estimations of $A$ and $C$ are obtained for $\lambda_{\ell_1} = 0.05$ and $\lambda_{TV} \in [0.5, 1]$. The error in $D$ is minimized for no or little $\ell_1$ regularization. These results show that the couples $(\lambda_{\ell_1}, \lambda_{TV})$ that estimate correctly $A$ and $D$ are close to the hyperparameters that minimizes the error in $C$. In addition, each error increases slowly in the neighborhood of the minimum, meaning that the estimation is not sensitive to the regularization in this neighborhood. This is of great interest since only the error in $C$ is available in practice. To assess the stability of the optimization scheme, the grid search is evaluated on 200 datasets generated with the same procedure. It is found that the results are highly reproducible, with a standard deviation in error through the datasets less than 2 (to be compared with the error range in Fig. 2). Finally, we notice that for large values of $\lambda_{\ell_1}$, a scale factor is introduced between $A$ and $D$ as the $\ell_1$ penalty drives $A$ toward small values which are compensated by higher values in $D$.

### B. Real dataset

The proposed model is used to decompose the rs-fMRI data of 8 mice. The fMRI signals of each mouse were preprocessed by the classical pipeline (inter-slice drift correction, motion correction, normalisation to Allen brain atlas<fnref n="1" />, smoothing,

<footnote n="1">¹Allen Mouse Brain Atlas available from mouse.brain-map.org</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 2: Grid search on $\lambda_{\ell_1}$ and $\lambda_{TV}$ (displayed in the best performing ranges). Frobenius norm of estimation error of the matrices $\mathbf{C}$ (left), $\mathbf{D}$ (center) and $\mathbf{A}$ (right) of synthetic data with SNR 0 dB. Marker $\mathbf{C}_{\min}$ (resp. $\mathbf{D}_{\min}$ and $\mathbf{A}_{\min}$) indicates the couple $(\lambda_{\ell_1}, \lambda_{TV})$ minimizing the reconstruction error for $\mathbf{C}$ (resp. $\mathbf{D}$ and $\mathbf{A}$).</figcaption></figure>

band pass filtering, regression of head motion and ventricles signal). The voxel signals were grouped into anatomical ROIs based on the atlas. The timecourse of each ROI is the spatial average of the timecourses of the voxels composing the region. An expert selected $P = 11$ networks (7 known to be highly active and 4 known to have moderate activity) composed of $R = 32$ ROIs. The support of the dictionary $\tilde{\mathbf{D}}$ is constructed from the definition of the FCUs as the complete graph of the ROIs of each network. The matrix $\mathbf{C}$ is obtained by applying sliding window correlation on the timecourses of the ROIs with a step of 1 sample ($T = 435$). The window is a door function of $\Delta_w = 55$ samples convolved with a Gaussian $\sigma = 3$ to avoid border effects.

The regularization parameters of the model, $\lambda_{\ell_1}$ and $\lambda_{TV}$, were optimized for each mouse by a grid search onto $\mathbf{C}$ similar to the procedure presented in Section V-A. For all mice, reconstruction error changes little with $\lambda_{\ell_1}$ and $\lambda_{TV}$, and the lowest error was achieved for no regularization. However, the regularization is not useless. Mathematically, the correlations in $\mathbf{C}$ are estimated on small windows (55 samples in this case), then the variance of the sample correlation is very high for small sample sizes, resulting in a very noisy correlation matrix, so the best reconstruction of $\mathbf{C}$ will overfit the noise. From a biological point of view, regularization drives the model to a more relevant solution because the fMRI data have smooth evolutions and not spurious spikes.

The regularization parameters were thus manually fixed by observing the estimated activations $\hat{\mathbf{A}}$. It was found that setting $\lambda_{TV} = 0.1$ is a good trade-off. Concerning $\ell_1$ penalty, $\lambda_{\ell_1}$ was set to zero as only 11 FCUs are used and 7 are believed to be highly active. The $\ell_1$ regularization would be of great importance in case the support of the dictionary $\tilde{\mathbf{D}}$ is not provided by expert and obtained by thresholding the correlation matrices through hypothesis testing (as in our previous work [4]). In that case, a very large number of FCUs would be used and the $\ell_1$ penalty would prevent redundancy in $\hat{\mathbf{A}}$ by cancelling the activity of spatially similar FCUs.

The result of applying the model on the timecourses of mouse 0 is displayed in Fig. 3. The actual correlation matrix computed on the timecourses (Fig. 3 (a)) clearly shows dynamics in the correlations of the pairs of ROIs at different time scales. The correlations are highly heterogeneous between the

<figure><img src="image_3.png" /><figcaption>Fig. 3: (a) Time evolution of correlation matrix $\mathbf{C}$ computed on real data for $\Delta_w = 55$ samples, (b) $\hat{\mathbf{C}} = \tilde{\mathbf{D}}\hat{\mathbf{A}}$. (c) Estimated activation matrix $\hat{\mathbf{A}}$. (d) $\hat{\mathbf{C}}_{\text{K-means}}$ estimated as the centroids of the states assigned to each instant in K-means [1] with $K = 7$ states, dashed lines indicate transitions between states.</figcaption></figure>

different pairs of ROI. The reconstructed matrix $\hat{\mathbf{C}}$ (Fig. 3 (b)) is very close to the original $\mathbf{C}$ with a RMSE of 0.21. The errors mainly reside in the negative correlations present in $\mathbf{C}$ that are not explained by the model. It should be noted that only the pairs of ROIs that are included in at least one FCU are displayed and used to compute the error as the other pairs can not be explained by the model (121 pairs out of 496). The activation matrix $\hat{\mathbf{A}}$ (Fig. 3 (c)) is coherent with expectations. The 7 FCUs believed to be highly active are indeed highly active (7 top rows in $\hat{\mathbf{A}}$) and the 4 FCUs that are supposed to be less active have a moderate activity.

Data also show significant anti-correlations. They could be explained by the left/right division of the brain regions in the atlas but are not considered by our model. This reinforces the idea of removing the positivity constraint on $\mathbf{D}$ thus allowing negative coefficients in $\hat{\mathbf{C}}$.

To check if a group trend can be found in the activations between subjects, the average activity of each FCU over time is computed, the distribution of these means across mice for each FCU are reported in the box plot of Fig. 4. There is little variance in the global activity of the FCUs across subjects. There is, on the contrary, a big variation in the activity between FCUs. Indeed, the model is able to clearly distinguish the FCU with high activity from those with low activity, suggesting the validity of the estimated dynamics.

From Fig. 4, FCU 2 is the most active, corresponding to the default mode network (DMN). The dynamics of this network are displayed in Fig. 5 by concatenating its estimated
<figure><img src="image_4.png" /><figcaption>Fig. 4: Box plot of the distribution of the time average activation of the FCUs.</figcaption></figure>

activation (third row of **A** for each mouse). For all mice, the DMN is highly active, with a slow oscillation between active and inactive states. If the oscillation look similar between mice, the experiment is too short to conclude on its frequency.

Fig. 6 reports for each mouse, the relative connectivity strength of the pairs of ROIs that compose the DMN (column of $\tilde{D}$, only the pairs that compose the FCU are displayed). For all subjects, two pairs are the most active, (Right ACAv, Right ACAd) and (Right RSP, Left RSP). The other pairs are also connected but at a relatively lower strength.

### C. Comparison with K-means

The widespread K-means method [1] was implemented by clustering the columns of the matrix **C**. The estimated matrix $\tilde{\mathbf{C}}_{\text{K-means}}$ is constructed by assigning a brain state (K-means cluster) to each time point and concatenating the centroid vectors of the corresponding states. We faced the problem of setting the number *K* of states, using *K* between 5 and 10 do not explain the variability in the data and using more states leads to incoherent results. When applied on the simulation data of Fig. 1, the states switches identify well most of the transitions at the cost of a large number *K* of states. For the mice fMRI data, the clustering was performed at subject level to estimate the states. Fig. 3 (b) and (d) display $\tilde{\mathbf{C}}$ and $\tilde{\mathbf{C}}_{\text{K-means}}$ with *K* = 7 states for mice 0. Our approach has the advantage that each FCU activity can exhibit its own time scale, unlike brain states, which synchronize the representation of the whole brain activity. Indeed, $\tilde{\mathbf{C}}_{\text{K-means}}$ is piecewise constant. In addition, the states are not repeated even for the relatively small *K*. A closer look at the K-means centroids revealed that they are very similar leading to switches between nearly identical states. Moreover, it is not possible to estimate the evolution of the activity of a network like the DMN using K-means, since the pairs of ROIs composing the DMN were found to be active in most of the states.

<figure><img src="image_5.png" /><figcaption>Fig. 5: Concatenation of the activity of the default mode network (DMN) FCU for each mouse.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 6: Concatenation of the estimated relative connectivity between pairs of ROI (column of $\tilde{D}$) of the default mode network (DMN) FCU for each mouse.</figcaption></figure>

## VI. CONCLUSION

This paper introduced a new representation of dynamic functional connectivity in the form of FCUs for single subject analysis. The expression of fMRI data into the proposed model is formulated as an optimization problem incorporating biologically inspired constrains. The problem was solved using joint optimization with sunsal-TV and ADMM algorithms. An original method for generating synthetic data allows validating the convergence of the algorithm, the stability of the solution and the positive effect of regularization. When applied on mice rs-fMRI data, our model was able to correctly identify the relative activity of the brain networks and achieved high coherence between individuals. The code for reproducing the experiments and generating data is available at the link https://github.com/massylmoudoud/DynaSTI_Eusipco2024

Future works will consider the definition of a heuristic to choose the regularization parameters, based on the analysis window length (for the total variation) and biological priors on the activity and sparsity of the FCUs. To overcome the scale indeterminacy factor, additional regularization term could be added on the matrix **D** to have all elements in [0, 1] or in [-1, 1] if negative correlations are considered, thus each column of **D** will represent a correlation pattern. The proposed method can be applied on the analysis of neurological pathologies both on preclinical and human data, to identify intergroup differences in the dFC, to study outlier-subjects within a group, or for longitudinal studies at a single-subject level.

## REFERENCES

[1] E. A. Allen, E. Damaraju, S. M. Plis, E. B. Erhardt, T. Eichele, and V. D. Calhoun, "Tracking whole-brain connectivity dynamics in the resting state," *Cerebral Cortex*, vol. 24, no. 3, pp. 663-676, 2014.

[2] A. Mensch, J. Mairal, B. Thirion, and G. Varoquaux, "Dictionary learning for massive matrix factorization," in *International Conference on Machine Learning*. PMLR, 2016, pp. 1737-1746.

[3] D. J. Lurie, D. Kessler, D. S. Bassett, R. F. Betzel, M. Breakspear, S. Kheilholz et al., "Questions and controversies in the study of time-varying functional connectivity in resting fMRI," *Network Neuroscience*, vol. 4, no. 1, pp. 30-69, 2020.

[4] A. Adam, C. Meillier, S. Achard, G. Becq, A. Bhanot, and A. Leborgne, "Modélisation parcimonieuse de la dynamique spatio-temporelle de la connectivité fonctionnelle en IRMF cérébrale." GRETSI, 2022.

[5] V. Portmann, C. Meillier, and V. Mazet, "Analyse de la dynamique spatio-temporelle de la connectivité fonctionnelle cérébrale : données synthétiques et modélisation." GRETSI, 2023.

[6] M.-D. Iordache, J. M. Bioucas-Dias, and A. Plaza, "Total variation spatial regularization for sparse hyperspectral unmixing," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 50, no. 11, pp. 4484-4502, 2012.

[7] S. Boyd, "Distributed optimization and statistical learning via the alternating direction method of multipliers," *Foundations and Trends in Machine Learning*, vol. 3, no. 1, pp. 1-122, 2010.