Data-driven modelling of brain activity using neural networks, Diffusion
  Maps, and the Koopman operator
We propose a machine-learning approach to model long-term out-of-sample
dynamics of brain activity from task-dependent fMRI data. Our approach is a
three stage one. First, we exploit Diffusion maps (DMs) to discover a set of
variables that parametrize the low-dimensional manifold on which the emergent
high-dimensional fMRI time series evolve. Then, we construct
reduced-order-models (ROMs) on the embedded manifold via two techniques:
Feedforward Neural Networks (FNNs) and the Koopman operator. Finally, for
predicting the out-of-sample long-term dynamics of brain activity in the
ambient fMRI space, we solve the pre-image problem coupling DMs with Geometric
Harmonics (GH) when using FNNs and the Koopman modes per se. For our
illustrations, we have assessed the performance of the two proposed schemes
using a benchmark fMRI dataset with recordings during a visuo-motor task. The
results suggest that just a few (for the particular task, five) non-linear
coordinates of the high-dimensional fMRI time series provide a good basis for
modelling and out-of-sample prediction of the brain activity. Furthermore, we
show that the proposed approaches outperform the one-step ahead predictions of
the naive random walk model, which, in contrast to our scheme, relies on the
knowledge of the signals in the previous time step. Importantly, we show that
the proposed Koopman operator approach provides, for any practical purposes,
equivalent results to the FNN-GH approach, thus bypassing the need to train a
non-linear map and to use GH to extrapolate predictions in the ambient fMRI
space; one can use instead the low-frequency truncation of the DMs function
space of L^2-integrable functions, to predict the entire list of coordinate
functions in the fMRI space and to solve the pre-image problem.
Ioannis K. Gallos, Daniel Lehmberg, Felix Dietrich, Constantinos Siettos
Data-driven modelling of brain activity using neural networks, Diffusion Maps, and the Koopman operator

###### Abstract

We propose a machine-learning approach to model long-term out-of-sample dynamics of brain activity from task-dependent fMRI data. Our approach is a three stage one. First, we exploit Diffusion maps (DMs) to discover a set of variables that parametrize the low-dimensional manifold on which the emergent high-dimensional fMRI time series evolve. Then, we construct reduced-order-models (ROMs) on the embedded manifold via two techniques: Feedforward Neural Networks (FNNs) and the Koopman operator. Finally, for predicting the out-of-sample long-term dynamics of brain activity in the ambient fMRI space, we solve the pre-image problem coupling DMs with Geometric Harmonics (GH) when using FNNs and the Koopman modes per se. For our illustrations, we have assessed the performance of the two proposed schemes using a benchmark fMRI dataset with recordings during a visuo-motor task. The results suggest that just a few (for the particular task, five) non-linear coordinates of the high-dimensional fMRI time series provide a good basis for modelling and out-of-sample prediction of the brain activity. Furthermore, we show that the proposed approaches outperform the one-step ahead predictions of the naive random walk model, which, in contrast to our scheme, relies on the knowledge of the signals in the previous time step. Importantly, we show that the proposed Koopman operator approach provides, for any practical purposes, equivalent results to the FNN-GH approach, thus bypassing the need to train a non-linear map and to use GH to extrapolate predictions in the ambient fMRI space; one can use instead the low-frequency truncation of the DMs function space of \(L^{2}\)-integrable functions, to predict the entire list of coordinate functions in the fMRI space and to solve the pre-image problem.

Brain dynamics task-fMRI Machine learning Diffusion maps Geometric harmonics Koopman operator Reduced order models Numerical Analysis

## 1 Introduction

Understanding and modelling the emergent brain dynamics has been a primary challenge in contemporary neuroscience (Bullmore and Sporns, 2009; Ermentrout and Terman, 2010; Jorgenson et al., 2015; Siettos and Starke, 2016;
Breakspear, 2017; Sip et al., 2023). Towards this goal, different approaches at a variety of scales have been introduced, ranging from individual neuron dynamics to the behaviour of different regions of neurons across the brain (Izhikevich, 2007; Deco et al., 2008; Bullmore and Sporns, 2012; Papo et al., 2014; Sietto and Starke, 2016; Breakspear, 2017; Sip et al., 2023). Some of them are physics-biologically informed models, thus derived directly from first principles (such as the Hodgkin-Huxley (Hodgkin and Huxley, 1952; Nelson and Rinzel, 1998; McCormick et al., 2007; Spiliotis et al., 2022), the FitzHugh-Nagumo (FitzHugh, 1955; Nagumo et al., 1962), phase-oscillators (Kopell and Ermentrout, 1986; Laing, 2017; Skardal and Arenas, 2020)) and neural mass models (Segneri et al., 2020; Taher et al., 2020), the Baloon model (Buxton et al., 1998) and Dynamic Causal Modelling (DCM) (Penny et al., 2004a)) (for a review of the various dynamic models see also (Izhikevich, 2007)). On the other hand, there is the data-driven approach, exploiting a wide range of methods (see also (Heitmann et al., 2018)) extending from independent component analysis (Beckmann and Smith, 2004), to Granger causality-based models (Seth, 2010; Seth et al., 2015; Friston et al., 2013; Protopapa et al., 2014; Kugiumtzis and Kimiskidis, 2015; Protopapa et al., 2016; Kugiumtzis et al., 2017; Almpanis and Siettos, 2020), phase-synchronization (Mormann et al., 2000; Rudrauf et al., 2006; Jirsa and Muller, 2013; Zakharova et al., 2014; Mylonas et al., 2016; Scholl, 2022) and network-based ones (Kopell and Ermentrout, 2004; Spiliotis et al., 2022; Petkoski and Jirsa, 2019). For an extended review of the various multiscale approaches, see (Siettos and Starke, 2016).

Most of the multiscale data-driven modelling approaches are "blended", i.e., they rely on EEG, MEG and fMRI data to calibrate biological-inspired models or to build surrogate models for the emergent brain activity. Within, this framework, machine learning (ML) has also been exploited to build surrogate dynamical models from neuroimaging data (Grossberg and Merrill, 1992; Niv, 2009; Richiardi et al., 2013; Suk et al., 2016; Gholami Doborjeh et al., 2018; Sun et al., 2022). However, when it comes to fMRI data, most-often consisting of hundreds of millions of measurements/features along space and time, one has to confront the so-called "curse of dimensionality". Problems arising from the "curse of dimensionality" in machine learning can manifest in a variety of ways, such as the excessive sparsity of data, the multiple comparison problem, the multi-collinearity of data, or even the poor generalization of the constructed surrogate models (Phinyomark et al., 2017; Altman and Krzywinski, 2018). To deal with these problems, methodological advances have been made involving random field theory and strict (e.g. corrected for multiple comparison problem) hypothesis testing to make reliable inferences in a unified framework (e.g. the Statistical Parametrical Mapping (SPM) (Penny et al., 2011)). To this day, this approach serves as the standard approach of processing fMRI time-series, relying mainly on a Generalized Linear Model approach to infer about the activity of certain voxels or regions of them. Despite the fact that SPM has been instrumental in analyzing task-related fMRI data, at the same time imposes critical assumptions about the structure and properties of data and, thus, setting limits to the problems that can be solved (Madhyastha et al., 2018).

Another effective way of dealing with problems related to the curse of dimensionality of brain signaling activity, is to represent the dataset on a low-dimensional embedding space (Belkin and Niyogi, 2003; Coifman and Lafon, 2006; Ansuini et al., 2019) using non-linear manifold learning (Gallos et al., 2021). For example, Qiu et al. (2015) used manifold learning for the prediction of dynamics of brain networks with aging by imposing a Log-Euclidean Riemannian manifold structure on the brain and then using a framework based on Locally Linear Embedding (Roweis and Saul, 2000) to uncover the manifold. Pospelcov et al. (2021) employed both linear and non-linear dimensionality reduction techniques with the aim of discriminating fMRI resting state recordings of subjects between acquisition and extinction of an experimental fear condition. The differences in terms of interclass-intraclass distance and the classification process took part on the low dimensional space, retaining up to 10 dimensions. Non-linear manifold learning algorithms outperformed the linear ones with Laplacian eigenmaps (Belkin and Niyogi, 2003) producing the best results. Gao et al. (2021) used Diffusion Maps to demonstrate that task-related fMRI data of different nature span to a similar low-dimensional embedding. In a series of works, (Gallos et al., 2021, 2021, 2022), we have applied linear and non-linear manifold learning techniques (such as Multi-Dimensional Scaling (Cox and Cox, 2008), ISOMAP (Tenenbaum et al., 2000), Diffusion Maps (Coifman and Lafon, 2006; Nadler et al., 2006, 2008)) to construct embedded functional connectivity networks from fMRI data. All the above studies provided evidence that dimensionality reduction may lead to a better classification performance, i.e., a finer detection of biomarkers. However, the above studies focused more on the clustering/classification problem than building dynamical models to predict the brain dynamics. Thiem et al. (2020) presented a data-driven methodology based on DMs to discover coarse variables, ANNs and Geometric Harmonics to learn the dynamic evolution laws of different versions of the Kuramoto model.

Here, we propose a three-step machine learning approach for the construction of surrogate dynamical models from real task fMRI data in the physical/ambient high dimensional space, thus dealing with the curse of dimensionality by learning the dynamics of brain activity on a low-dimensional space. In particular, we first use Diffusion Maps to learn an effective non-linear low-dimensional manifold. Then we learn and predict the embedded brain dynamics based on the set of variables that span the low-dimensional manifold. For this task, we use both FNNs and the Koopman Operator framework (Mezic, 2013; Williams et al., 2015; Brunton et al., 2016; Li et al., 2017; Bollt et al.,
2018; Dietrich et al., 2020; Lehmberg et al., 2021). Finally, we solve the pre-image problem, i.e. the reconstruction of the predictions in the original ambient fMRI space, using Geometric Harmonics (GA) (Coifman and Lafon, 2006; Dsilva et al., 2013; Papaioannou et al., 2021; Evangelou et al., 2023) for the FNN predictions, and the Koopman modes per se (Li et al., 2017; Lehmberg et al., 2021) for the Koopman operator predictions.

For our illustrations, we used a benchmark fMRI dataset, recorded during an attention to a visual motion task. The proposed approach provides out-of-sample long-term predictions, thus reconstructing accurately in the ambient space the brain activity in response to the visual/attention stimuli in the five brain regions that are known to be the key ones for the particular task from previous studies. Furthermore, we show that the Koopman operator numerical-analysis/based scheme provides numerical approximations with similar accuracy to the FNN-GH scheme, while bypassing both the construction of a non-linear surrogate model on the DMs, and the use of GH to extrapolate predictions in the ambient fMRI space.

## 2 Benchmark dataset and Signal extraction

For our illustrations, we focus on a benchmark fMRI dataset of a single subject during a task that targets attention to visual motion. This particular dataset has served as a benchmark in many studies (Buchel and Friston, 1997; Friston et al., 2003; Penny et al., 2004b; Almpanis and Siettos, 2020) and is publicly available for download from the official SPM website [https://www.fil.ion.ucl.ac.uk/spm/data/attention/](https://www.fil.ion.ucl.ac.uk/spm/data/attention/). The images are already smoothed, spatially normalised, realigned and slice-time corrected.

In the experiment (Buchel and Friston, 1997), subjects observed a black screen displaying a few white dots. The experimental design consisted of specific epochs where the dots behave as either static or moving. In between these epochs, there was also an intermediate epoch where only a static picture without dots was present (i.e. a fixation phase, which is treated as baseline). The subjects were guided to keep focused for any change concerning the moving dots, even, when no changes in the intermediate epochs really existed. Thus, there were four distinct experimental conditions: "fixation", "non attention" where the dots were moving, but subjects did not need to pay attention on the screen, "attention" where subjects needed to pay attention on existing (if any) changes regarding the movement of dots (i.e., possible changes in velocity or acceleration) and "static" where the dots were stationary. The dataset consists of 4 runs stuck together, where the 10 first scans of each run were discarded for the elimination of non-desirable magnetic effects. Thus, the length of the fMRI dataset is 360 scans.

Here, we have first parcellated the brain into 116 brain regions of interest (ROIs) as derived by the use of Automated Anatomical Labeling (AAL) (Tzourio-Mazoyer et al., 2002) from the fMRI data. After the formulation of the ROIs by AAL, the average of the Blood Oxygen Level Dependent (BOLD) signal from all voxels located in each region was calculated at each time point. Out of all 116 ROIs and due to the limited field of view of the fMRI acquisition, some parts of the cerebellum were excluded as there was no existing signal (i.e. there was barely any change in the BOLD signal throughout the experiment). Specifically, these regions were the right and left part of Cerebellum_10, Cerebellum_7b and left part of Cerebellum_9. The remaining 111 time series were linearly detrended and standardized for further analysis.

## 3 Methodology

Our data-driven machine-learning methodology for modelling grain activity deploys in three main steps. First, as a zero-step, we follow the standard procedure regarding fMRI data preprocessing with the general linear model to identify the voxels of statistically significant BOLD activity. Then, at the first step, we use parsimonious Diffusion Maps to identify a set of variables that parametrize the low-dimensional manifold where the embedded fMRI BOLD time-series evolve. Based on them, we construct reduced-order models (ROMs), based either on FNNs or Koopman operator. Finally, at the third step, we use either Geometric Harmonics (GH) (Coifman and Lafon, 2006), when we use FNNs, or Extended Dynamic Model Decomposition (EDMD) (Williams et al., 2015) when we use Koopman operator, to solve the pre-image problem, i.e. to provide out-of-sample predictions in the ambient phase space of BOLD signals. In the following sections, we describe the above steps in detail. A schematic of the three-step methodology is also shown in Fig. 1.

### Step 0. Data Preprocessing using the General Linear Model

For the past 20 years, the General Linear Modelling (GLM) approach has been the cornerstone of standard fMRI analysis (Friston et al., 1994, 1995; Worsley and Friston, 1995). Within the GLM framework, time series associated with BOLD signal at the voxel level are modeled as a weighted sum of one or more predictor variables (i.e. in
relationship to the experimental conditions) including an error term for the unobserved/unexplained information in the model. The ultimate goal of the aforementioned approach is to determine whether a predictor variable does contribute (and to what extent) to the variability observed in the data (e.g., a known pattern of stimulation induced by specific experimental settings).

Specifically, we can model the amplitude of the BOLD signal with fMRI in the \(i\)-th brain region, as a weighted sum of some known predictor variables \(\mathbf{u}_{1},...,\mathbf{u}_{p}\) each scaled by a set of parameters \(\mathbf{\beta}\):

\[\begin{array}{l}\mathbf{x_{i,1}}=\mathbf{u_{1,1}}\mathbf{\beta_{1}}+\mathbf{u_{1,2}}\mathbf{ \beta_{2}}+...\mathbf{u_{1,p}}\mathbf{\beta_{p}}+\mathbf{\epsilon_{i,1}}\\ \mathbf{x_{i,2}}=\mathbf{u_{2,1}}\mathbf{\beta_{1}}+\mathbf{u_{2,2}}\mathbf{\beta_{2}}+...\mathbf{u_{2,p}}\mathbf{\beta_{p}}+\mathbf{\epsilon_{i,2}}\\ \mathbf{...}\\ \mathbf{x_{i,N}}=\mathbf{u_{N,1}}\mathbf{\beta_{1}}+\mathbf{u_{N,2}}\mathbf{\beta_{2}}+...\mathbf{u_{N,p}}\mathbf{\beta_{p}}+\mathbf{\epsilon_{i,N}}\end{array} \tag{1}\]

A more compact notation leads to the following formula:

\[\mathbf{x_{i}}=\mathbf{U}\mathbf{\beta}+\mathbf{\epsilon}, \tag{2}\]

where \(\mathbf{x_{i}}=[\mathbf{x_{i,1}},\mathbf{x_{i,2}},\ldots,\mathbf{x_{i,N}}]^{T}\in\mathbb{ R}^{N}\) is the column vector that stores the BOLD responses in all time instances, \(\mathbf{U}=\begin{bmatrix}\mathbf{u_{1,1}},&\ldots&\mathbf{u_{1,p}}\\ \cdots&\ldots&\cdots\\ \mathbf{u_{N,1}}&\ldots&\mathbf{u_{N,p}}\end{bmatrix}\in\mathbb{R}^{N\times p}\) is the design matrix with the column \(\mathbf{j}\)-th representing the time series of the predictor variable \(\mathbf{j}\) and \(\mathbf{\beta}=[\mathbf{\beta_{1}},\mathbf{\beta_{2}},\ldots,\mathbf{\beta_{p}}]^{T}\in\mathbb{ R}^{p}\) is a vector of the unknown parameters that set the contribution of each of the predictors to the BOLD signal \(\mathbf{x_{i}}\). Finally, \(\mathbf{\epsilon_{i}}=[\mathbf{\epsilon_{i,1}},\mathbf{\epsilon_{i,2}},\ldots,\mathbf{ \epsilon_{i,N}}]^{T}\in\mathbb{R}^{N}\) is a vector containing the corresponding (modelling) errors.

In this study, the standard GLM approach was used inside the SPM framework (Penny et al., 2011) and its toolbox extension AAL (Tzourio-Mazoyer et al., 2002). Out of the 3 available methods for automatic anatomical labelling of the global functional activation map, we used the cluster labeling. This is a rational choice since we wanted to localize the activated clusters in the level of anatomical regions. The final results contain which percentage of the activated cluster belongs to a specific anatomical region, and also which percentage of the region is part of the activated cluster (Tzourio-Mazoyer et al., 2002).

Here, the experimental variables \(\mathbf{x_{i}}\) are the following: "fixation", "static", "attention", "non-attention". We considered as activated clusters, 10 or more nearby voxels that passed the contrast of "attention", "non attention" and "static" vs

Figure 1: Schematic of the proposed machine learning based methodology for modelling and predicting the dynamics of brain activity from task-depended fMRI data.
"fixation" which is here treated as a low level baseline (Buchel et al., 1998). We considered two different thresholds, one of \(\mathbf{p<0.05}\) corrected for a family wise error (FWE) and another more liberal threshold of \(\mathbf{p<0.01}\) (uncorrected). The analysis was restricted only on voxels included in the predefined regions of AAL (eg. by using the atlas as an inclusive mask). In this way, we wanted to determine the anatomical regions that are actually "active" throughout the experiment in respect with any of the experimental variables other than fixation (e.g., when either when subject is paying or not paying attention to moving or static dots).

Step 1: Identification of an Embedded Manifold in the fMRI BOLD signals via Parsimonious Diffusion Maps

Diffusion maps is a non-linear dimensionality reduction/ manifold learning algorithm introduced by Coifman and Lafon (2006a) that exploits the inherent relationship among heat diffusion and random walk Markov chain. The main idea is that jumping from one data point to another in a random walk regime, it is more likely to jump to nearby points than points that are further away. The algorithm produces a mapping of the dataset in the Euclidean space, whose coordinates are computed utilizing the eigenvalues and the eigenvectors of a diffusion operator on the data. An embedding in the low-dimensional diffusion maps space is obtained by the projections on the eigenvectors of a normalized graph Laplacian (Belkin and Niyogi, 2003). Specifically, given a dataset of \(\mathbf{N}\) points, \(\mathbf{X}=[\mathbf{x_{1}}\quad\mathbf{x_{2}}\quad\mathbf{\dots}\quad\mathbf{x_{N}}]\), \(\mathbf{x_{i}\in\mathbb{R}^{M}}\) is the column vector containing the BOLD signal in \(\mathbf{N}\) instances of the \(\mathbf{i}\)-th brain region, \(\mathbf{M}\) corresponds to the number of brain regions, we construct an affinity matrix \(\mathbf{W}\in\mathbb{R}^{\mathbf{N}\times\mathbf{N}}\), defined by a kernel, any square integrable function, \(\mathbf{w}\in\mathbf{L^{2}:\mathbf{X}\times\mathbf{X}\rightarrow\mathbb{R}+^{0}}\), between all pairs of points \(\mathbf{x_{i}}\) and \(\mathbf{x_{j}}\). The affinity matrix is most often computed though the use of the so-called Gaussian heat kernel with elements:

\[\mathbf{w_{i,j}=\exp\left(-\frac{||\mathbf{x_{i}}-\mathbf{x_{j}}||_{2}^{2}}{2\mathbf{\sigma}} \right), \tag{3}\]

where \(||.||_{2}\) computes the distance between two points using the Euclidean \(\mathbf{L_{2}}\) norm (any metric could be used) between points, and \(\mathbf{\sigma}\) serves as a scale parameter of the kernel. A Gaussian heat kernel in Eq. 3 produces both symmetric and positive semi-definite affinity square matrix \(\mathbf{W}\). The property is fundamental and the one that allows for the interpretation of weights as scaled probabilities. We thus formulate the diagonal normalization matrix \(\mathbf{K}\) with elements

\[\mathbf{k_{ii}}=\sum_{\mathbf{j=1}}^{\mathbf{N}}\mathbf{w_{ij}}. \tag{4}\]

Here, a family of anisotropic diffusions can be parameterized by adding a parameter \(\mathbf{\alpha}\) which controls the amount of influence of the density in the data. This constant can be used to normalize the affinity matrix \(\mathbf{W}\) as:

\[\tilde{\mathbf{W}}=\mathbf{K}^{-\mathbf{\alpha}}\mathbf{W}\mathbf{K}^{-\mathbf{\alpha}} \tag{5}\]

For \(\mathbf{\alpha=0}\), we get the normalized graph Laplacian, \(\mathbf{\alpha=0.5}\), we get the Fokker-Planck diffusion, and for \(\mathbf{\alpha=1}\), we get the Laplace-Beltrami operator (Thiem et al., 2020).

Next, we renormalize the matrix \(\tilde{\mathbf{W}}\) with the diagonal matrix \(\mathbf{\tilde{K}}\), with elements \(\mathbf{\tilde{k}_{ii}=\sum_{\mathbf{j=1}}^{\mathbf{N}}\tilde{\mathbf{w}_{ij}}.\ \mathbf{P \in\mathbb{R}^{\mathbf{N}\times\mathbf{N}}}}\) is a Markovian row stochastic matrix which is also called diffusion matrix given by:

\[\mathbf{P=\tilde{K}^{-1}\tilde{\mathbf{W}}.} \tag{6}\]

In this stage, the elements of the diffusion matrix can be seen as the scaled probabilities of one point jumping to another in a random walk sense. Consequently, taking the power \(\mathbf{t}\) of the diffusion matrix \(\mathbf{P}\) is essentially identical of observing \(\mathbf{t}\) steps forward of a Markov chain process on the data points. The element \(\mathbf{P^{t}(\mathbf{x_{i}},\mathbf{x_{j}})}\) denotes the transition probability of jumping from point \(\mathbf{x_{i}}\) to point \(\mathbf{x_{j}}\) in \(\mathbf{t}\) steps.

Applying singular value decomposition (SVD) on \(\mathbf{P}\), we get:

\[\mathbf{P=\Psi\ \Lambda\Psi^{T},} \tag{7}\]

with \(\mathbf{\Lambda}\) being the diagonal matrix that stores the \(\mathbf{M}\) eigenvalues and \(\mathbf{\Psi}\in\mathbb{R}^{\mathbf{N}\times\mathbf{N}}\) the eigenvectors of \(\mathbf{P}\).

These eigenvectors are the discrete approximations of the eigenfunctions of the corresponding to the parameter \(\mathbf{\alpha}\) (see above) kernel operator on the manifold. Thus, such eigenfunctions provide an orthonormal basis of the \(\mathbf{L^{2}}\) space, i.e. of the space of square-integrable functions (Coifman and Lafon, 2006a).

The eigenvalues that correspond to the eigenvectors in descending order are \(\mathbf{\lambda_{0}=1\geq\lambda_{1}\geq\lambda_{2}\...\geq\lambda_{k}}\) with the first being the trivial eigenvalue that is equal to 1 (\(\mathbf{P}\) is a Markovian matrix) and \(\mathbf{k}\) being the target embedding
dimension. The embedding is obtained by mapping the points from the original/ ambient space to the diffusion maps space by preserving the diffusion distance among all data points in the dataset (Nadler et al., 2006). Since the diffusion distance takes into account all possible paths between data points, it is thus more robust against noise perturbations (Coifman and Lafon, 2006a) (comparing for example with the geodesic distance).

The coordinates of the nonlinear projection of a point \(\mathbf{x}_{i}\) in the \(\mathbf{k}\) dimensional space spanned by the first \(\mathbf{k}\) DMs \(\{\mathbf{\psi_{1}},\mathbf{\psi_{2}},\dots,\mathbf{\psi_{k}}\}\) read:

\[\mathbf{\mathcal{R}}(\mathbf{x}_{i})\equiv\mathbf{y}_{i}=[\mathbf{y}_{i,1},\mathbf{y}_{i, 2},\dots,\mathbf{y}_{i,k}]^{T}=[\mathbf{\lambda}_{i}^{t}\mathbf{\psi}_{i,1},\mathbf{\lambda}_{ 2}^{t}\mathbf{\psi}_{i,2},\dots,\mathbf{\lambda}_{k}^{t}\mathbf{\psi}_{i,k}]^{T},\mathbf{i}= \mathbf{1},\mathbf{2},\dots,\mathbf{N}. \tag{8}\]

where, \(\mathbf{\psi}_{i,\mathbf{l}}\) is the \(\mathbf{i}\)-th element of \(\mathbf{\psi_{1}},\mathbf{l}=\mathbf{1},\mathbf{2},\dots,\mathbf{k}\). For all practical purposes, the embedding dimension is determined by the spectral gap in the eigenvalues of the final decomposition. A numerical gap (or a sharp decrease) between the first few eigenvalues and the rest of the eigenspectrum would indicate that a few eigenvectors could be adequate for the approximation of the diffusion distance between all pairs of points (Coifman and Lafon, 2006a). A compact version of the algorithm is the following:

1. **Input:** linearly detrended data set \(\mathbf{X}=[\mathbf{x}_{1}\quad\mathbf{x}_{2}\quad\dots\quad\mathbf{x}_{N}]\), \(\mathbf{x}_{i}\in\mathbb{R}^{M}\).
2. Compute the Gaussian heat kernel with elements \(\mathbf{w}_{i,\mathbf{j}}=\mathbf{exp}\left(-\frac{||\mathbf{x}_{i}-\mathbf{x}_{j}||_{ 2}^{2}}{2\mathbf{\sigma}}\right),\mathbf{i},\mathbf{j}=\mathbf{1},\mathbf{2},\dots,\mathbf{N}\).
3. Apply the normalization on \(\mathbf{W}\) (Eq. 5) and renormalization (Eq. 6) to compute the diffusion matrix \(\mathbf{P}\);
4. Apply SVD on matrix \(\mathbf{P}\) (Eq. 7) and keep the first \(\mathbf{k}\) eigenvectors corresponding to the \(\mathbf{k}\) largest eigenvalues.
5. **Output:**\(\mathbf{\Psi}_{\mathbf{k}}=[\mathbf{\psi_{1}}\quad\mathbf{\psi_{2}}\quad\dots\mathbf{\psi}_{k}] \in\mathbb{R}^{N\times\mathbf{k}}\), \(\mathbf{y}_{i}=[\mathbf{y}_{i,1},\mathbf{y}_{i,2},\dots,\mathbf{y}_{i,k}]^{T}\), \(\mathbf{i}=\mathbf{1},\mathbf{2},\dots,\mathbf{N}\) (see Eq.(8).

Next, from the initial set of eigenvectors, we retain only the \(\mathbf{d}\) parsimonious eigendimensions as proposed in (Dsilva et al., 2018). Specifically, using a local linear function,

\[\mathbf{\psi}_{i,\mathbf{l}}\approx\mathbf{c}_{i,\mathbf{l}}+\mathbf{\beta}_{i,\mathbf{l}}^{T}\mathbf{ \Psi}_{i,\mathbf{l}-1},\quad\mathbf{\Psi}_{i,\mathbf{l}-1}=[\mathbf{\psi}_{i,1},...,\mathbf{\psi} _{i,\mathbf{l}-1}]^{T}\in\mathbb{R}^{\mathbf{l}-1}, \tag{9}\]

\(\mathbf{\alpha}_{i,\mathbf{l}}\in\mathbb{R},\mathbf{\beta}_{i,\mathbf{l}}\in\mathbb{R}^{\mathbf{l }-1}\) we approximate each point of \(\mathbf{\psi}_{i,\mathbf{l}}\) based on the remaining data points \(\mathbf{\psi}_{i,\mathbf{l}-1}\). This leads to the following optimization problem (Thiem et al., 2020):

\[\underset{\mathbf{c},\mathbf{\beta}}{\text{argmin}}\,\sum_{i\neq j}\mathbf{K}\Big{(}\mathbf{ \Psi}_{i,\mathbf{l}-1},\mathbf{\Psi}_{j,\mathbf{l}-1}\Big{)}\Big{(}\mathbf{\psi}_{j,\mathbf{l}}- \left(\mathbf{c}+\mathbf{\beta}^{T}\mathbf{\Psi}_{j,\mathbf{l}-1}\right)\Big{)}^{2}, \tag{10}\]

where \(\mathbf{K}\) is the Gaussian kernel. Finally, the computation of the normalized error \(\mathbf{c}\mathbf{r}_{\mathbf{k}}\) is done through the use of leave-one out cross validation for each of the local linear fit. The normalized error is defined as:

\[\mathbf{c}\mathbf{r}_{\mathbf{l}}=\sqrt{\frac{\sum_{i=1}^{N}\left(\mathbf{\psi}_{i,\mathbf{l}}- \left(\mathbf{\alpha}_{i,\mathbf{l}}+\mathbf{\beta}_{i,\mathbf{l}}^{T}\mathbf{\Psi}_{i,\mathbf{l}-1} \right)\right)}{\sum_{i=1}^{N}(\mathbf{\psi}_{i,\mathbf{l}})^{2}}}. \tag{11}\]

A small or negligible error \(\mathbf{c}\mathbf{r}_{\mathbf{l}}\) denotes that \(\mathbf{\psi}_{\mathbf{l}}\) can be actually predicted from the remaining eigenvectors \(\mathbf{\psi_{1}},\mathbf{\psi_{2}},...,\mathbf{\psi}_{\mathbf{l}-1}\) and thus is a repeated eigendirection (i.e., \(\mathbf{\psi}_{\mathbf{l}}\) is considered a harmonic of the previous eigenmodes). Therefore, only the eigenvectors that exhibit a large \(\mathbf{c}\mathbf{r}_{\mathbf{l}}\) are selected in a way of seeking the most parsimonious representation. More information regarding the choice of the most parsimonious eigenvectors can be found in (Dsilva et al., 2018), (Lee et al., 2020) and (Galaris et al., 2022).

A compact pseudocode for identifying parsimonious eigenvectors is the following:

1. **Input:** Set of \(\mathbf{\Psi}_{\mathbf{k}}\in\mathbb{R}^{N\times\mathbf{k}}\), \(\mathbf{k}<<\mathbf{M}\) corresponding to the largest eigenvalues from the application of Diffusion Maps and \(\mathbf{d}\) the number of parsimonious eigenvectors to retain.
2. Solve the optimization problem given by Eq.(10) and calculate the normalized error \(\mathbf{c}\mathbf{r}_{\mathbf{l}}\) given by Eq.(11).
3. **Output:**\(\mathbf{\Psi^{\prime}}_{\mathbf{d}}\in\mathbb{R}^{N\times\mathbf{d}}\), with columns the parsimonious DMs \(\mathbf{\psi}_{\mathbf{j}}^{\prime}\in\mathbb{R}^{N},\mathbf{j}=\mathbf{1},\mathbf{2},\dots,\mathbf{d}\) corresponding to the \(\mathbf{d}\) largest \(\mathbf{c}\mathbf{r}_{\mathbf{l}}\)s.

### Step 2. Construction of Diffusion Maps-based Reduced order models

After the reduction of space dimension to a set of \(\mathbf{d}\) of parsimonious DMs components, we proceed to the construction of reduced order (regression) models (ROMs).
Here, we have used two approaches for the construction of ROMs and lifting back to the original ambient space: (a) via FNNs and Geometric Harmonics, and, (b) via the Koopman modes/EDMD. In both cases, the general form of a ROM reads:

\[\boldsymbol{\psi^{\prime}_{i+1,t}}=\boldsymbol{f}(\boldsymbol{\Psi^{\prime}_{i,t }},\boldsymbol{I_{i,s}},\mathbf{b})+\boldsymbol{\epsilon_{i,t}},\quad\boldsymbol {l=1,2,\ldots,d,i=1,2,\ldots,N-1}, \tag{12}\]

where \(\boldsymbol{\Psi^{\prime}_{i,d}}=[\boldsymbol{\psi^{\prime}_{i,1}},\ldots, \boldsymbol{\psi^{\prime}_{i,t}},\ldots,\boldsymbol{\psi^{\prime}_{i,d}}]^{T} \in\mathbb{R}^{d}\), is the vector containing the \(\boldsymbol{i}\)-th elements of the \(\boldsymbol{d}\) parsimonious DMS (\(\boldsymbol{\psi^{\prime}}_{i,t}\) denotes the \(\boldsymbol{i}\)-th element of the \(\boldsymbol{l}\)-th parsimonious DM), \(\boldsymbol{b}\) denotes a set of parameters of the model, \(\boldsymbol{\mathbf{I_{i,s}}}\) the external stimuli \(\boldsymbol{s}\) indicating specific experimental conditions at the time instant \(\boldsymbol{i}\), and \(\boldsymbol{\epsilon_{i,t}}\) is the modelling error. We note that after training, predictions were done iteratively: we set the initial conditions at the last point of the training set, and then iterate the ROMs to produce the long-term predictions.

In order to compare the efficiency of the two approaches, we applied also a simple naive random walk (NRW) model in a direct way for one-step ahead predictions: predictions at time \(\boldsymbol{i}+\boldsymbol{1}\) are equal to the last observed values, i.e., \(\boldsymbol{\psi_{i+1,t}}=\boldsymbol{\psi_{i,t}}\).

Thus, in contrast to the predictions via the embedded ROMs that are produced iteratively, i.e. without any knowledge of the previous values, but the initial conditions, the predictions made via the NRW model are based on the knowledge of the values of the previous points of the test set.

#### 3.3.1 Reduced order models with FNNs

As a first reduced order modelling approach on the discovered manifold \(\boldsymbol{\mathcal{M}}\), we used \(\boldsymbol{d}\) FNNs, for each one of the embedded coordinates \(\boldsymbol{\psi^{\prime}_{i}}\), \(\boldsymbol{l=1,2,\ldots,d}\), with one hidden layer, \(\boldsymbol{H}\) units and a linear output layer with \(\boldsymbol{n}\) units that can be written compactly as

\[\boldsymbol{\psi^{\prime}}_{i+1,t}=\boldsymbol{\mathbf{w^{T}_{1,o}}}\mathbf{S }(\boldsymbol{\mathbf{W^{T}_{1,t}}}\boldsymbol{\mathbf{\Psi^{\prime}_{i,d}}} +\mathbf{b_{1,t}})+\boldsymbol{b_{o,t}},\quad\boldsymbol{l=1,2,\ldots,d}. \tag{13}\]

\(\mathbf{S}(\boldsymbol{\cdot})\) is the activation function (based on the above formulation it is assumed to be the same for all networks and nodes in the hidden layer), \(\mathbf{b_{1,t}}\in\mathbb{R}^{\boldsymbol{H}}\) is the vector containing the biases of the nodes of the first hidden layer, \(\mathbf{W_{1,t}}\in\mathbb{R}^{\boldsymbol{d}\times\boldsymbol{H}}\) is the matrix containing the weights between the input and the first hidden layer, \(\mathbf{w_{o,t}}\in\mathbb{R}^{\boldsymbol{H}}\) denotes the vector containing the weights between the hidden layer and the output layer, and finally \(\boldsymbol{b_{o,t}}\) is the bias of the output.

The total number of input units were matched to the embedding coordinates of the manifold plus the external stimuli \(\boldsymbol{\mathbf{I}}\) (i.e., \(\boldsymbol{d}\) embedding coordinates plus 4 external stimuli: attention, non attention, static and fixation) while there was one output unit. We used the logistic transfer function \(\boldsymbol{S(x)}=\frac{\boldsymbol{1}}{\boldsymbol{1+e}^{-x}}\) as the activation function for all neurons in the hidden layer. A learning rate decay parameter was also used as a regularization parameter to prevent over-fitting and improve generalization (Krogh and Hertz, 1992) of the final model. The neural networks were trained on the first 280 time points (which accounts for roughly the 77% of the data points) using repeated (10 times) 10-fold cross validation approach. All parameters like different number of neurons \(\boldsymbol{a}\) and learning rate decay parameter \(\boldsymbol{c}\) were optimized via grid search inside the cross validation procedure. Finally, we optimized five final models (i.e., one model for each one of the parsimonious eigenvectors selected) to make one-step predictions in the following way: we give \(\boldsymbol{\psi_{i,t}}\) as input to the FNN in order to make prediction on \(\boldsymbol{\psi_{i+1,t}}\). Using the best candidate model for each one of the parsimonious components \(\boldsymbol{\psi_{t}}\), we predicted iteratively the left out/ unseen test points (i.e., the next 80 unseen time points).

#### 3.3.2 Reduced order models with the Koopman operator

In the Koopman operator framework, predictions are performed in the function space over the data manifold, and so state prediction turns into coordinate function prediction (Mezic, 2005; Budisic et al., 2012; Dietrich et al., 2020). The defining property in this framework is that the operator always acts linearly on its domain, which makes it amenable to spectral analysis and approximation techniques from numerical analysis linear algebra. Given the flow map \(\boldsymbol{f}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}\) of the embedded dynamics on the manifold, such that \(\boldsymbol{\Psi^{\prime}_{i+1,t}}=\boldsymbol{f}(\boldsymbol{\Psi^{\prime}_{i, t}})\), the Koopman operator \(\boldsymbol{U:\mathcal{F}}\rightarrow\mathcal{F}\) acts on observables \(\boldsymbol{g}\in\mathcal{F}\) by composition with \(\boldsymbol{f}\), such that

\[[\boldsymbol{U}\boldsymbol{g}]\left(\boldsymbol{\Psi^{\prime}_{i,t}}\right)= \boldsymbol{g}(\boldsymbol{f}(\boldsymbol{\Psi^{\prime}_{i,t}}))=\boldsymbol{g }(\boldsymbol{\Psi^{\prime}_{i+1,t}}). \tag{14}\]

If a function \(\boldsymbol{\phi_{j}}\) is an eigenfunction of \(\boldsymbol{U}\) with eigenvalue \(\boldsymbol{\omega_{j}}\), then \([\boldsymbol{U}\boldsymbol{\phi_{j}}]\left(\boldsymbol{\Psi^{\prime}_{i,t}} \right)=\boldsymbol{\phi_{j}}(\boldsymbol{\Psi^{\prime}_{i+1,t}})=\boldsymbol{ \omega_{j}}\boldsymbol{\phi_{j}}(\boldsymbol{\Psi^{\prime}_{i,t}})\). For many systems, eigenfunctions of the Koopman operator span a large subspace in \(\mathcal{F}\), so that many observables \(\boldsymbol{g}\in\mathcal{F}\) can be written as a linear combination of eigenfunctions, \(\boldsymbol{g}=\sum_{\boldsymbol{j}}\boldsymbol{c_{j}}\boldsymbol{\phi_{j}}\), with coefficients \(\boldsymbol{c_{j}}\in\mathbb{C}\). These coefficients are often called "Koopman modes" of the observable \(\boldsymbol{g}\).

The most prevalent numerical algorithm to approximate the Koopman operator is dynamic mode decomposition (DMD) (Schmid, 2010, 2022). It approximates a linear map between the given observables of a system and their
future state, assuming that the observables span a function subspace invariant to the dynamics. In our case, we apply DMD to five of the embedding diffusion map coordinates \(\mathbf{\Psi^{\prime}_{d}}\) (this can be interpreted as a vector-valued observable, i.e., the Koopman operator is acting on five coordinates simultaneously (Budisic et al., 2012)).

### Step 3. Predictions in the Ambient fMRI Space: Solution of the Pre-Image Problem

The final step is to "lift"predictions back to the high dimensional ambient space (here a 111-dimensional space of brain regions), i.e. solve the so-called "pre-image" problem.

For linear manifold learning methods such as the Principal Component Analysis (PCA), this is trivial. However, for non-linear manifold learning algorithms such as DMs, there is no explicit inverse mapping and no unique solution (see Chiavazzo et al. (2014); Papaioannou et al. (2021); Evangelou et al. (2023)) and one has to construct a lifting operator \(\mathbf{\mathcal{L}}=\mathbf{\mathcal{R}}(\mathbf{X})\rightarrow\mathbf{X}\) for new unseen samples on the manifold \(\mathbf{y_{\star}}\notin\mathbf{\mathcal{R}}(\mathbf{X})\).

The problem regarding an unseen data point is often referred to as the "out-of-sample extension" problem. While this problem traditionally addresses the direct problem (i.e., the extension of a function \(\mathbf{g}\) on the ambient space for new unseen data points \(\mathbf{x_{\star}}\notin\mathbf{X}\)), here, we are interested in the inverse problem (i.e., the lifting of predictions made on the manifold back to the ambient space). Here, as discussed, the out-of-sample solution of the pre-image problem is solved using Geometric Harmonics when using FNNs for ROMs, or with Extended Dynamics Mode Decomposition when using the Koopman operator. Below, we present both methods.

#### 3.4.1 Geometric Harmonics

Typically, the out-of-sample extension problem is solved through the Nystrom extension methodology (Nystrom, 1929) derived from the solution of the Fredholm equation of the second kind:

\[\mathbf{f}(\mathbf{t})=\mathbf{g}(\mathbf{t})+\mathbf{\mu}\int_{\mathbf{a}}^{\mathbf{b}}\mathbf{k}(\mathbf{t},\mathbf{s })\mathbf{f}(\mathbf{s})\mathbf{ds}, \tag{15}\]

where \(\mathbf{k}(\mathbf{t},\mathbf{s})\) and \(\mathbf{g}(\mathbf{t})\) are known functions while \(\mathbf{f}(\mathbf{t})\) is the unknown. The Nystrom method approximates the integral

\[\int_{\mathbf{a}}^{\mathbf{b}}\mathbf{y}(\mathbf{s})\mathbf{ds}\approx\sum_{\mathbf{j}-\mathbf{1}}^{\mathbf{N} }\mathbf{w}_{\mathbf{j}}\mathbf{y}(\mathbf{s}_{\mathbf{j}}), \tag{16}\]

where \(\mathbf{N}\) are the collocation points and \(\mathbf{w}\) the weights determined by, for example, the Gauss-Jacobi quadratic rule. Using Eq.16 in 15 and evaluating \(\mathbf{f}\) and \(\mathbf{g}\) at the collocation points \(\mathbf{N}\), we can make the following approximation:

\[(\mathbf{I}-\mathbf{\mu}\mathbf{\tilde{K}})\mathbf{\hat{f}}=\mathbf{g}, \tag{17}\]

where \(\mathbf{\tilde{K}}\)=\(\mathbf{\tilde{k}_{ij}}=\mathbf{k}(\mathbf{s_{i}},\mathbf{s_{j}})\mathbf{w_{j}}\). Solution to the homogenous Fredholm problem (where \(\mathbf{g}=\mathbf{0}\)) is provided through the solution of the eigenvalue problem:

\[\mathbf{\tilde{K}}\mathbf{\hat{f}}=\frac{\mathbf{1}}{\mathbf{\mu}}\mathbf{\hat{f}}, \tag{18}\]

i.e.,

\[\sum_{\mathbf{j}=\mathbf{1}}^{\mathbf{N}}\mathbf{w}_{\mathbf{j}}\mathbf{y}(\mathbf{s_{i}},\mathbf{s_{j}})\mathbf{ \hat{f}}_{\mathbf{j}}=\frac{\mathbf{1}}{\mathbf{\mu}}\mathbf{\hat{f}}_{\mathbf{i}},\qquad \mathbf{i}=\mathbf{1},\mathbf{2},\ldots,\mathbf{N} \tag{19}\]

where \(\mathbf{\hat{f}_{\mathbf{i}}}=\mathbf{\hat{f}}(\mathbf{s_{i}})\) the \(\mathbf{i}\)-th component of \(\mathbf{\hat{f}}\). The Nystrom extension in the full domain for \(\mathbf{N}\) number of collocation points at a point \(\mathbf{x}\) is given by:

\[\mathbf{\mathcal{E}}(\mathbf{f}(\mathbf{x}))=\mathbf{\hat{f}}(\mathbf{x})=\mathbf{\mu}\sum_{\mathbf{j}- \mathbf{1}}^{\mathbf{N}}\mathbf{w}_{\mathbf{j}}\mathbf{k}(\mathbf{x},\mathbf{s_{j}})\mathbf{\hat{f}_{\mathbf{j}}}. \tag{20}\]

The above Eq.(20) provides a map from high dimensional to the reduced order space and back.

Since the eigenvectors of the DMs form a basis for the embedded manifold, the extension is formulated (Coifman and Lafon, 2006) by the expansion of \(\mathbf{f}(\mathbf{x_{i}})\) using the first \(\mathbf{d}\) parsimonious eigenvectors of the diffusion matrix \(\mathbf{P^{t}}\):

\[\mathbf{\hat{f}}(\mathbf{x_{i}})=\sum_{\mathbf{i}=\mathbf{1}}^{\mathbf{d}}\mathbf{a}_{\mathbf{t}}\psi^{ \prime}_{\mathbf{i}}(\mathbf{x_{i}}),\qquad\mathbf{i}=\mathbf{1},\mathbf{2},\ldots,\mathbf{ N}, \tag{21}\]
where \(\mathbf{a_{l}}=\left\langle\psi^{\prime}{}_{l},\mathbf{f}(\mathbf{x_{i}})\right\rangle\), \((\left\langle\cdot,\cdot\right\rangle\) denotes inner product), \(\mathbf{f}\in\mathbb{R}^{N}\) is the vector containing the values of the function \(\mathbf{f}\) at the \(\mathbf{N}\) points. Thus, the Nystrom extension of \(\mathbf{f}\) at a new unseen point \(\mathbf{x_{i}^{\star}}\) using the same coefficients \(\mathbf{a_{k}}\), reads:

\[\mathcal{E}(\hat{\mathbf{f}}(\mathbf{x_{i}^{\star}}))=\sum_{l=1}^{d}\mathbf{a_{l} }\hat{\psi_{l}}(\mathbf{x_{i}^{\star}}), \tag{22}\]

where,

\[\hat{\psi}(\mathbf{x_{i}^{\star}})=\frac{1}{\mathbf{\lambda_{i}^{\star}}}\sum_{j=1 }^{N}\mathbf{k}(\mathbf{x_{i}^{\star}},\mathbf{x_{j}})\psi^{\prime}{}_{l}(\mathbf{ x_{j}}),\qquad l=1,2,\ldots,d, \tag{23}\]

are the geometric harmonics extension of each one of the parsimonious eigenvectors to the unseen data. Thus, for a new set \(\mathbf{X^{\star}}\) of \(\mathbf{n}\) points we get:

\[\mathcal{E}(\hat{\mathbf{f}})=\mathbf{K}_{n\times N}\Psi^{\prime}{}_{N \times d}\mathbf{\Lambda}_{d\times d}^{-1}\Psi^{\prime}{}_{d\times N}^{T}\mathbf{ f}_{N\times 1}, \tag{24}\]

where \(\mathbf{K}\) is the \(\mathbf{n}\times\mathbf{N}\) kernel matrix, \(\mathbf{\Lambda}\) is a \(d\times d\) diagonal matrix that stores the eigenvalues of the eigenvectors.

Here, for the solution of the pre-image problem, we used the approach of "double Diffusion maps" (Evangelou et al., 2023; Papaioannou et al., 2021; Patsatzis et al., 2023)). A pseudocode for the solution of the pre-image problem is given below:

1. **Input:**\(\Psi^{\prime}{}_{d}\in\mathbb{R}^{N\times d}\) matrix whose columns are the \(d\) parsimonious diffusion maps (from Eq.(3.2)).
2. For each new point \(\{\mathbf{y_{i}^{\star}}\}_{i=1}^{n}\) in the DMs embedded space computed using the Nystrom extension (Eq.(22), compute the kernel matrix \(\mathbf{K}\in\mathbb{R}^{n\times N}\) with elements \(\mathbf{k_{i,j}}=\exp\left(-\frac{||\mathbf{y_{i}^{\star}}-\mathbf{y_{j}}||_{2}^{ 2}}{2\sigma}\right)\), \(\mathbf{i}=\mathbf{1},\mathbf{2},\ldots,\mathbf{n}\), \(\mathbf{j}=\mathbf{1},\mathbf{2},\ldots,\mathbf{N}\).
3. **Output:** Reconstructed points in the ambient space using Eq.(24).

#### 3.4.2 Koopman Operator and Extended Dynamic Mode Decomposition (EDMD)

When using the Koopman operator approach, the DMs coordinates considered as a truncated basis of a function space over the original measurements, can be used effectively similarly as when using Dynamic Mode Decomposition (DMD) (see Williams et al., 2015), but on the original measurements in the ambient fMRI space.

Thus, unlike standard DMD, in our case, the Koopman modes are computed for the original coordinates, not the DM eigenfunctions, so that we obtain a linear map from Koopman eigenfunctions to original measurements. A compact version of the EDMD algorithm we use is as follows:

1. **Input:** Diffusion maps coordinates \(\mathbf{\Psi^{\prime}_{i,d}}\in\mathbb{R}^{N\times d}\) at time steps \(\mathbf{i}=\mathbf{1},\ldots,\mathbf{N}\) from Section 3.2.
2. Approximate the Koopman matrix \(\hat{\mathbf{U}}=\mathbf{\Psi^{\prime}_{1}}\Psi^{\prime}{}_{-1}\in\mathbb{R}^{d \times d}\), where \(\dagger\) denotes the pseudo-inverse and \(\mathbf{\Psi^{\prime}_{-}}=\left[\mathbf{\Psi^{\prime}_{1,d}},\mathbf{\Psi^{\prime}_{2,d}},\ldots,\mathbf{\Psi^{\prime}_{N-1,d}}\right]\in\mathbb{R}^{d\times(N-1)}\) and \(\mathbf{\Psi^{\prime}_{+}}=\left[\mathbf{\Psi^{\prime}_{2,d}},\mathbf{\Psi^{\prime}_{3,d}},\ldots,\mathbf{\Psi^{\prime}_{N,d}}\right]\in\mathbb{R}^{d\times(N-1)}\) are time shifted data matrices with \(\mathbf{\Psi^{\prime}_{i,d}}\in\mathbb{R}^{d}\) as columns.
3. Compute all eigenvectors \(\mathbf{\hat{\phi}_{j}}\in\mathbb{R}^{d}\) and eigenvalues \(\mathbf{\hat{\omega}_{j}}\) of \(\hat{\mathbf{U}}\) by solving \(\hat{\mathbf{U}}\mathbf{\hat{\phi}_{j}}=\mathbf{\hat{\omega}_{j}}\mathbf{\hat{\phi}_{j}},\mathbf{j}= \mathbf{1},\mathbf{2},\ldots,\mathbf{d}\).
4. Obtain the Koopman modes \(\mathbf{c_{j}}\in\mathbb{R}^{M}\) for the original measurements \(\mathbf{x_{i}}\in\mathbb{R}^{M}\), by solving the linear systems \(\mathbf{x_{i}}=\sum_{j}\mathbf{c_{j}}\mathbf{\hat{\phi}_{i,j}}\) with data from all available time steps \(\mathbf{i}=\mathbf{1},\ldots,\mathbf{N}\) simultaneously.
5. **Output:** Eigenvectors \(\mathbf{\hat{\phi}_{j}}\), eigenvalues \(\mathbf{\hat{\omega}_{j}}\), and Koopman modes \(\mathbf{c_{j}}\) for the original measurements.
6. **Prediction:** To obtain an approximation of \(\mathbf{x_{i+1}}\), we evaluate \(\hat{\mathbf{x_{i+1}}}=\hat{\mathbf{U}}\sum_{\mathbf{j}}\mathbf{c_{j}}\mathbf{\hat{\phi}_{i,j}} =\sum_{j}\mathbf{\hat{\omega}_{j}}\mathbf{c_{j}}\mathbf{\hat{\phi}_{i,j}}\).

## 4 Results

For our computations, we have used the Python packages Datafold (Lehmberg et al., 2020) with scikit-learn (Pedregosa et al., 2011). For the implementation of the FNNs we utilized the R package "nnet" (Ripley et al., 2016).
### Data Preprocessing

Our analysis starts with the data-preprocessing using GLM. In Figure 2, the thresholded image for the activation during the contrast "attention" + "non-attention" + "stationary" vs "fixation" is presented as an overlay on the brain extracted structural T1 MNI image (2mm space). We considered activations of clusters of 10 or more voxels with two different levels of significance: A) \(\mathbf{p<0.05}\) with Family Wise Error correction (FWE) and B) \(\mathbf{p<0.001}\) (uncorrected). In Table 1, the results considering the more stringent threshold of \(\mathbf{p<0.05}\) (FWE corrected) are presented in detail. We report the local maximum of each cluster, its size in number of voxels, the brain regions that contribute to each cluster and finally the percentage of the contribution. Similarly, in Table 2 the results concerning the activations of clusters using a more liberal threshold of \(\mathbf{p<0.001}\) (uncorrected) are reported. The brain regions that showed activations were mainly parts of the Occipital Lobe such as the Calcarine, Cuneus, Occipital gyrus, the Fusiform gyrus and parts of the Cerebellum. When a more liberal threshold was applied (\(\mathbf{p<0.001}\)), the clusters would naturally get larger. Consequently, two more regions appeared to have clusters of activation, namely, the Postcentral gyrus and the inferior parietal cortex. Focusing on these "active" regions, we then proceeded to the prediction of their behaviour through macroscopic variables extracted by the Diffusion Maps.

Figure 2: Thresholded image of the contrast “attention” + “non-attention” + “stationary” vs “fixation” is presented as an overlay on the reference structural MRI (T1) image.A) Level of significance p<0.05 (FWE corrected), B) \(\mathbf{p<0.001}\) (uncorrected).
\begin{table}
\begin{tabular}{l l l l} \multicolumn{3}{c}{Level of significance p \(<\) 0.05 (FWE corrected)} \\ \hline Local maximum x,y,z mm & Number of voxels in the cluster & Region Label & \% in the cluster \\ \hline
3, -93, 0 & 191 & Calarine L & 11.96 \\  & 191 & Cuneus L & 10.17 \\  & 191 & Occipital Sup L & 8.89 \\  & 191 & Calcarine R & 4.72 \\  & 191 & Occipital Mid L & 0.21 \\  & 191 & Lingual R & 0.15 \\ \hline -33, -84, -21 & 28 & Cerebellum Crus1 L & 39.29 \\  & 28 & Lingual L & 39.29 \\  & 28 & Fusiform L & 17.86 \\  & 28 & Occipital Inf. L & 3.57 \\ \hline
21, -66, -3 & 20 & Lingual R & 100 \\ \hline
27, -81, -18 & 11 & Cerebellum 6 R & 45.45 \\  & 11 & Cerebellum Crus1 R & 27.27 \\  & 11 & Fusiform R & 18.18 \\  & 11 & Lingual R & 9.09 \\ \end{tabular}
\end{table}
Table 1: Clusters of 10 or more activated voxels for the contrast “attention” + “non-attention” + “stationary” vs “fixation” using SPM. The level of significance was set to \(\boldsymbol{p<0.05}\) (FWE corrected). For each cluster, the coordinates of the local maximum (in standard MNI coordinates) is shown along with the total number of voxels in the cluster and the brain regions that these voxels are found. Finally, the percentage of a region’s voxels in each cluster is also presented.
### Manifold Learning with Diffusion Maps

Manifold learning with Diffusion maps was applied on the first 280 (out of the whole 360 time points) time points of the 111 linearly detrended time series. The last 80 points were left out to evaluate later the prediction error based on ROMs and the "lifting" process. In Figure 3, we present the eigenspectrum of the DMs decomposition. The red vertical line shows the number of the extracted eigenvectors (more eigenvectors would not affect the final outcome). The five most parsimonious DMs components that we used for further analysis are marked with red arrows, namely \(\mathbf{\psi_{1}}\), \(\mathbf{\psi_{5}}\), \(\mathbf{\psi_{8}}\), \(\mathbf{\psi_{13}}\) and \(\mathbf{\psi_{15}}\). As it can be seen in Figure 3, after the first 30 eigenvalues, most of the variance has already been captured. In other words, consequent eigenvectors have a very small contribution to the final embedding. The parameters we actually set in our computations are \(\mathbf{t}=\mathbf{0}\), \(\mathbf{\sigma}=\mathbf{65}\), \(\mathbf{\alpha}=\mathbf{1}\), \(\mathbf{\kappa}=\mathbf{30}\) and \(\mathbf{d}=\mathbf{5}\).

### Construction of Reduced Order Models and Predictions on the Manifold

Out-of-sample predictions on the embedded manifold were based on the 5 parsimonious DMs, namely \(\mathbf{\psi_{1}}\), \(\mathbf{\psi_{5}}\), \(\mathbf{\psi_{8}}\), \(\mathbf{\psi_{13}}\) and, \(\mathbf{\psi_{15}}\) are shown in Figure 4. Based on the these, we constructed ROMs using FNNs and the Koopman operator as described in subsection 3.3.

The ROMs were trained via one time step ahead predictions using the first 280 time points, which account to the 77% of the total data points. Validation was done using 10-fold cross validation, repeated 10 times. Thus, the performance of the ROMs was assessed by simulating the trained ROMs iteratively: setting the initial conditions at the time point 280, we iterated the ROMs to produce the next 80 points (i.e. up to the 360th time point).

\begin{table}
\begin{tabular}{l l l l}  & \multicolumn{2}{c}{Level of significance p \textless{} 0.001 (uncorrected)} \\ \hline Local maximum (x,y,z mm) & Number of voxels in the cluster & Region Label & \% in the cluster \\ \hline
3, -93, 0 & 917 & Calcarine L & 20.17 \\  & 917 & Occipital Mid L & 13.52 \\  & 917 & Lingual R & 12.65 \\  & 917 & Occipital Sup L & 11.01 \\  & 917 & Calcarine R & 10.58 \\  & 917 & Fusiform R & 10.25 \\  & 917 & Cuneus L & 7.96 \\  & 917 & Lingual L & 6 \\  & 917 & Cerebellum Crus1 R & 2.94 \\  & 917 & Cerebellum 6 R & 1.74 \\  & 917 & Cuneus R & 1.64 \\  & 917 & Occipital Inf. L & 1.2 \\  & 917 & Occipital Sup R & 0.22 \\  & 917 & Occipital Inf. R & 0.11 \\ \hline -33, -84, -21 & 68 & Cerebellum Crus1 L & 44.12 \\  & 68 & Lingual L & 25 \\  & 68 & Fusiform L & 19.12 \\  & 68 & Occipital Inf. L & 11.76 \\ \hline -39, -90, -3 & 10 & Occipital Mid. L & 90 \\  & 10 & Occipital Inf. L & 10 \\ \hline
39, -39, 54 & 15 & Parietal Inf. R & 80 \\  & 15 & Postcentral R & 20 \\ \hline -33, -69, -18 & 35 & Fusiform L & 94.29 \\  & 35 & Cerebellum 6 L & 2.86 \\  & 35 & Lingual L & 2.86 \\ \hline
27, -90, 20 & 11 & Occipital Sup. R & 72.73 \\  & 11 & Occipital Mid.R & 27.27 \\ \end{tabular}
\end{table}
Table 2: Clusters of 10 or more activated voxels for the contrast “attention” + “non-attention” + “stationary” vs “fixation” using SPM. The level of significance was set to p < 0.001 (uncorrected). For each cluster, the coordinates of the local maximum (in standard MNI coordinates) is shown, along with the total number of voxels in the cluster and the brain regions that these voxels are found. Finally, the percentage of a region’s voxels in each cluster is also presented.
For the FNNs, we used the logistic function as the activation function for all neurons in the hidden layer and a learning rate decay parameter or regularization to prevent over-fitting and improve generalization of the final model. Parameter tuning was done via grid search during the repeated cross validation process.

### Predictions in the fMRI space and comparison with the Naive Random Walk Model

At the final step, the out-of-sample predictions on the manifold made by FNNs and the Koopman operator were "lifted" back to the original space using GH and Koopman modes per se, respectively (see 3.4). Specifically, the predicted values were of size \(\mathbf{80\times d=5}\) (reduced space) and the lifting of those predictions led to a new set of size \(\mathbf{80\times 111}\) (original space). In Figure 5, we depict those predictions overlaid on the unseen (out-of-sample) test data (red color) for the two approaches, namely for the FNN-GH and the Koopman operator. Indicatively, we show the first four regions of interest (ROIs) corresponding to the biggest activated clusters (see Table 2), specifically, the left Calcarine sulcus (A), the left Middle Occipital gyrus (B), the right Lingual gyrus (C) and the left Superior Occipital gyrus (D).

In Table 3, we present the errors on the test set for each ROIs as derived for each of the two schemes (FNN- GH and the Koopman operator).

To assess the prediction efficiency of the proposed methodology, we also provide the errors obtained when applying the naive random walk (NRW) model in a direct way (the predictions at time \(\boldsymbol{i+1}\) are the last observed values at time \(\boldsymbol{i}\)). Hence, to predict the values of the DMs at the time point \(\boldsymbol{i+1}\) of the test set with the NRW model, we assumed that _we know the actual values of the DMs at the time \(\boldsymbol{i}\) of the test set_. The errors are reported in terms of the RMSE and the \(\boldsymbol{L_{2}}\) norm; the smallest errors are marked with bold. In general, the two strategies behave similarly well, outperforming the NRW predictions for all ROIs except one, the Cerebelum Crus 1 on the right hemisphere. This is important, revealing the efficiency of the proposed approach, since as discussed, the predictions of the ROMs were produced iteratively, i.e. _without any knowledge of the values in the test set_, while the predictions made with the NRW model _were based on the knowledge of the values of the previous points of the test set_.

Here, we should note that predictions made for some brain regions, like the Cuneus of the right hemisphere, might seem to exhibit a poor or spurious predictions. This is due to the fact that in this particular case, the volume of voxels that are "activated" (only for the uncorrected threshold of \(\boldsymbol{p<0.001}\)) is found to contribute only 1.64 % in

Figure 3: The eigenspectrum of DMs on the 111 time series that correspond to different brain regions. The red vertical line indicates how many eigenvectors we computed. The five parsimonious eigenvectors that were finally used as macroscopic variables are also marked with red arrows.
the activated cluster (see Table 2). Since in this study we try to predict the behaviour of a whole brain region, each time series is averaged over all voxels of the predefined region (here, using AAL). Thus, when trying to predict the "average" behaviour of voxels that constitute a region when only a small proportion of them is actually "activated", can't be reliable. Therefore, we would normally expect the average signal of such a region to be noisy and without a pattern which could be attributed to the task imposed to the subject (here, the attention to visual motion).

Figure 4: Predictions based on the 5 parsimonious eigenvectors (namely \(\mathbf{\psi_{1},\psi_{5},\psi_{8},\psi_{13}}\) and \(\mathbf{\psi_{15}}\)) when applying FNNs and the Koopman operator. Actual points for each one of the DMs are presented with black color up to the 280th point, which is the last point of the training set. The predicted values are presented with light green color up to the end of the time series (as derived iteratively by the ROMs from the point 281 to the end.).
## 5 Discussion

Our work addresses a three-step machine-learning based approach for the modelling of brain activity based on task-depended fMRI data in order to deal with the curse of dimensionality. In the first step, we used parsimonious Diffusion Maps to create a basis that spans a low dimensional subspace, a manifold, in which the information of the brain activity is retained. In a second step, we trained ROMs using this low-dimensional basis, thus relaxing the curse of dimensionality, using both FNNs and the Koopman operator. Finally, we predicted brain activity in the fMRI space by solving the pre-image problem using Geometric Harmonics and the Koopman modes. By doing so, we demonstrated that for the particular benchmark task-experiment, the brain response to the attention to visual motion task is contained in a 5-dimensional manifold. By lifting our predictions to the ambient fMRI space, we found that the key brain regions during attention to visual motion are actually four: the Calcarine Sulcus, the Cuneus, the Superior Occipital Gyrus and Lingual Gyrus. These regions are the same regions reported also in (Buchel et al., 1998) "re-discovered" by our methodology.

Our proof-of-concept study is the first to propose such a machine-learning framework to model and importantly accurately predict brain activity from task-based fMRI data, thus facing challenges that are not met when producing controlled experiments from synthetic data produced just by model simulations.

In the numerical experiments detailed above, we show that the performance of the two schemes FNNs-GH and the Koopman operator is comparable, while both outperform the naive random walk whose predictions assume the knowledge of brain activity at the previous time step. Importantly, the comparison between the two schemes reveals that it is not always necessary to construct a non-linear embedding (here, with Diffusion Maps) and then also construct surrogate non-linear models using, e.g., FNNs to predict future states. Effectively, one can utilize the low-frequency truncation of the function space of square-integrable functions (\(\mathbf{L^{2}}\)) over the original data, as obtained with Diffusion Maps, within the Koopman operator framework, to predict the entire list of coordinate functions in a linear fashion. Conversely, the Neural Network approach treats the embedding coordinates as a base space for the dynamics and predicts them point-wise.

Figure 5: Predictions in the amplitude of the BOLD signal in the original fMRI space for four of the most “activated” regions as found by the classical GLM methodology A) left Calcarine Sulcus, B) left Middle Occipital Gyrus, C) right Lingual Gyrus, D) left Superior Occipital Gyrus. The red color marks the actual values of the test set, while the other colors correspond to the predictions based on the proposed methodology.
We believe that the proposed methodology may trigger further developments in the field, thus providing the base for a general framework for modelling the dynamics of brain activity, from high-dimensional time series, including also the solution of the source localization problem combining for example simultaneous EEG-fMRI recordings.

## 6 Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## 7 Acknowledgements

F.D. was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - project no. 468830823.
