I Costantini

P Filipiak

K Maksymenko

R Deriche

Deconvolution S Deslauriers-Gauthier

Keywords: BOLD, Deconvolution, Edge detection, Functional MRI, Hemodynamic Response Function, Image Regularization, Paradigm Free, Partial Diferential Equations, Resting-state v

The advent of new brain imaging techniques such as resting-state functional MRI (fMRI), has led to the need for new approaches to recover brain functional activations without a prior knowledge on the experimental paradigm, as it was the case for task-fMRI. Conventional methods, i.e. the general linear model, requires the knowledge of the task paradigm to estimate the contribution of each voxel's time course to the given task. To overcome this limitation, approaches to deconvolve the blood-oxygen-level-dependent (BOLD) response and recover the underlying neural activations without necessity of prior information has been proposed. Supposing the brain activates in constant blocks, first we propose a temporal regularized deconvolution technique which uses an exponential operator, whose shape and performance can be adjusted, into a least absolute shrinkage and selection operator (LASSO) model solved via the Least-Angle Regression (LARS) algorithm. We reduced the number of parameters to be set by the user, when compared with the state of the art. Second, we introduce a paradigm-free regularization algorithm that applies on the 4-D fMRI image, acting simultaneously in the 3-D space and the 1-D time dimensions. The approach is based on the idea that large image variations should be preserved as they occur during an activation, whereas small variations should be smoothed to remove noise. It allows to smooth the whole fMRI image with an anisotropic regularization, thus blindly recovering the location of the brain activations in space and their timing and duration. Both approaches were tested on phantom and real data and were demonstrated to improve the results obtained in the state of the art.

Resumé

L'avènement de nouvelles techniques d'imagerie du cerveau comme l'IRM fonctionnelle (IRMf) au repos a conduit à la nécessité de nouvelles méthodes pour récupérer les activations fonctionnelles du cerveau sans connaissance du paradigme expérimental, comme dans l'IRMf basée sur tâche. Les méthodes conventionnelles, par example le modèle linéaire général, nécessitent la connaissance de la tâche pour pouvoir estimer la contribution du signal de chaque voxel à la tâche donnée. Pour surmonter ces limitations, des méthodes de déconvolution de la réponse dépendant du niveau d'oxygène dans le sang et de récupération des activations neurales sans avoir besoin d'informations préalables ont été proposées. Dans cette thèse, nous proposons d'abord une technique de déconvolution avec une regularisation temporelle qui utilise un opérateur exponentiel, dont la forme et la performance peuvent être ajustées. Avec cette méthode, nous avons réduit le nombre de paramètres à régler par l'utilisateur, par rapport à l'état de l'art. Ensuite, nous avons introduit un algorithme de régularisation qui s'applique à l'image IRMf 4-D, agissant simultanément dans les dimensions spatiale et temporelle. La méthode est basée sur l'idée que les grandes variations de l'image doivent être préservées car elles se produisent lors d'une activation et les petites variations doivent être lissées pour éliminer le bruit. Elle permet de lisser l'image IRMf avec une régularisation anisotrope, récupérant ainsi aveuglément la localisation des activations cérébrales et leur durée. Les deux méthodes ont été testées sur des données synthetiques et réelles et ont démontré une amélioration des résultats de l'état de l'art. which induces the activation for each voxel i. u s (i, t) is the derivative of u(i, t) which represents the activation and it is assumed to be a block-type signal. h(t) is the hemodynamic response function, that convolved with the activityinducing signal u(i, t) gives the activity-related signal x(i, t). y(i, t) represents the measuref fMRI signals for each voxel i, which is obtained by adding noise to the activity-related signal. (Figure adapted from Karahano glu et al. [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF]). . time series and all the other brain voxels is calculated. The outcome is a map that shows which areas show a high level of functional connectivity with the selected seed [START_REF] Van Den | Exploring the brain network: a review on resting-state fmri functional connectivity[END_REF]. I α is the exponential accumulation function that leads from s(t) to the activity-inducing signal u(t), which is piece-wise constant. H represents the hemodynamic response function, x(t) is the activity-related signal, a is the additive noise and y(t) is the simulated acquired fMRI time course. . . . . . . [START_REF] Buxton | Modeling the hemodynamic response to brain activation[END_REF] 3.3 Impulse response of the accumulation function I α . The plot shows how the accumulator changes increasing α, e.g. from 0.2 (green) to 1 (blue). . . . . . . [START_REF] Buxton | Dynamics of blood flow and oxygenation changes during brain activation: the balloon model[END_REF] 3.4 L-curve. The graph represents for one time series of one voxel the plot of the two quantities, one related to the norm of the solution (abscissa) and the other related to the residuals (ordinate), i.e. (||ŝ|| 1 , ||y -Aŝ|| 2 2 ), paramtrized by the regularization parameter λ. Each black dot in the curve is related to each λ outputted by the Least Angle Regression (LARS) algorithm [START_REF] Efron | Least angle regression[END_REF]. The red dot represents the optimal solution as the nearest to (0, 0) ∈ R 2 . . . . . . . . . . . . 56 3.5 fMRI signal simulation. u(t) is the activity-inducing signal, a block of constant activation represented with a piece-wise constant signal. b is the block type noise and m the model noise which were consecutively added to the activityinducing signal u(t) thus leading to u n (t) before the convolution with the hemodynamic response function H. x(t) is the activity-related signal, a is the additive gaussian noise and y(t) is the simulated acquired fMRI time course. [START_REF] Lund | Non-white noise in fmri: does modelling have an impact?[END_REF] 3.6 Each plot shows the root mean square errors (RMSE) and the roots of standard deviations computed between the ground truth activation and the estimated one using the Least Angle Regression (LARS) combined with the L-curve (in red) and with the leave-one-out cross-validation (LOO-CV, in light blue).

Mot-clés: IRMf -Activation Cérébrale -Reconstruction Régularisée -Regularisation Anisotropique

List of Figures

Each plot refer to a different experiment, as reported in Table 3.3. For each experiment, the RMSE were averaged across 100 repetitions of the simulated time series. The errors (in y-axes) are plotted with respect to the α value (in x-axes) used for the α-filter. In the 5 experiments, we added block-type noise ( b ) for false activations, and both m and a to the data, increasing it from experiment 1 to 5. The plots show that the LOO-CV-bases approach is more sensitive to noise while increasing it and results are more disperse if compared with the L-curve-based ones. In all experiments the results obtained with Lcurve show errors that are significantly lower than those given by the LOO- Each plot refer to a different experiment, as reported in Table 3.3. For each experiment, the RMSE were averaged across 100 repetitions of the simulated time series. The RMSE (in y-axes) are plotted with respect to the α value (in x-axes) used for the α-filter. In experiments from 6 to 10, we increasingly added a to the synthetic times series. The plots show that the LOO-CV-based approach is less sensitive to noise if compared to experiments 1 to 5 in Figure 3.6, but still more sensitive to noise while if compared with the L-curve-based ones. In all experiments the results obtained with L-curve show errors that are significantly lower than the ones obtained with the LOO-CV for α values lower than 0.5. For α > 0.75 results are similar. . . . . . . . . . . . . . . . . . . 3.8 The graphs in the plots show examples of the recovered activation using the mixed Least Angle Regression (LARS)-L-curve approach (in red) and the mixed LARS and leave-one-out cross-validation (LOO-CV) (in light blue) superimposed on the ground truth activation in black and the noisy synthetic fMRI time series (in green). Each row shows a different experiment, meaning different noises applied to the ground truth activation and each columns refer to a specific α used in the α-filter. In experiments 1, 3 and 5 the added noise was coming from different sources: block-type noise ( b ) to simulate falseactivations, model noise m and additive noise a . The standard deviations of m and a increases from experiment 1 to 5. The plots show that the recovered activations using the L-curve-based approach are much closer to the ground truth compared with LOO-CV. The L-curve-based approach shows results that are closer to the ground truth in terms of amplitude for bigger alphas. For α = 0.75, 3, results are similar, while in contrast with those given by α = 0.3.

3.9

The graphs in the plots report examples of the recovered activation using the mixed Least Angle Regression (LARS)-L-curve approach (in red) and the mixed LARS and leave-one-out cross-validation (LOO-CV) (in light blue) superimposed on the ground truth activation in black and the noisy synthetic fMRI time series (in green). Each row shows a different experiment (see Table 3.3), meaning that different amount of noises were applied to the ground truth activation Each columns refer to a specific α used in the α-filter. In experiments 6 and 10 time series are corrupted with a with increasing standard deviation.

The different amount of noise are described in Table 3.3. The recovered activations using the L-curve-based approach are closer to the ground truth compared with LOO-CV. The L-curve-based approach shows results that are closer to the ground truth in terms of amplitude for greater alphas α = 0.75, 3, in contrast with α = 0. different durations, as proposed by Farouj et al. [START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF]. The activation on the right, B, is made of a single long block of simulated activity. . . . . . . . . . . 3.12 Reconstructed activity-inducing signal û(t) obtained with our approach (red) superimposed on the ground truth (u(t), black) and simulated fMRI signal (y(t), green). The plots in the blue square are taken from exemplificative voxels and are related to activation A with pSNR = 5.17 dB; the ones in the yellow square are related to activation B with SNR = 5.12 dB. Vertically, for both A and B, the plots show that our approach is able to discern between an activation and a non-activation. . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.13

Reconstructed activity-inducing signal û(t) obtained with our approach (red) and the total activation (TA, blue) superimposed on the ground truth activation (black) and simulated fMRI signal (green). The plot on the top is related to activation A with pSNR = 5.17 dB. The plot on the bottom is related to activation B with pSNR = 5.12 dB. . . . . . . . . . . . . . . . . . . . . . . . .

3.14

The plot shows the û(t) obtained with the Human Connectome Project (HCP) data. Results obtained using our approach (red) and the TA (blue) superimposed on the real fMRI signals (green) were all averaged within the region of interest located in the Broadmann Area 4p. In the x-axes the time

(a)

From left to right: spatial maps of the simulated functional MRI (fMRI) image y, ground truth activation u, recovered activation using the Total Activation (TA) approach ( ûTA ) and our approach ( ûPFf MRI ). Each row corresponds to a a different peak-SNR (pSNR): 6.54 dB, 5.99 dB, 5.9 dB, 3.93 dB from the top to the bottom. (b) Reconstructed time series û(t) obtained with our approach ûPF-f MRI (t) (red) and the TA approach ûTA (t) (blue) superimposed on the ground truth activation u(t) (black) and fMRI signal y(t) (green). A zoom of the plot is shown in Figure 4.4. . . . . . . . . . . . . .

4.4

The plot is a zoom of Figure 4.3 representative of a trend that we found in all the four experiment showed in Figure 4.3. From top left: spatial maps of a slice of the simulated functional MRI (fMRI) image y, ground truth activation u, recovered activation maps obtained using the Total Activation (TA) approach ( ûTA ) and the PF-fMRI ( ûPFf MRI ). The map obtained using TA ( ûTA ) had lower amplitude compared with the ground truth, nonetheless the scale between the values were kept. The map obtained using PF-fMRI ( ûPFf MRI ) showed amplitude comparable to the ground truth, the scale between the values were kept. The time series on the bottom shows that the recovered signal ûPF-f MRI (t) (red) was closer to the ground truth u(t) (black) compared to ûTA (t) (blue). . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

-1 ∑ n=-∞ g(n) is equal to +∞ ∑ n=0 g(n). This means that, if we compute +∞ ∑ n=-∞ g(n) = -1 ∑ n=-∞ g(n) + +∞ ∑ n=0 g(n) this is equivalent to +∞ ∑ n=-∞ g(n) = 2 +∞ ∑ n=0 g (

Problem

The problem we are addressing in this thesis is the recovery of brain functional activations from fMRI data. Traditional techniques for measuring the brain function are based on experimental paradigms, where the subject is asked to perform a task in order to measure the difference between the rest condition, representing a baseline, and the task condition. fMRI data analysis methods, i.e. the general linear model (GLM), require the knowledge of the task paradigm to be able to estimate the contribution of each brain region to the given task. Nevertheless, an experimental setup is not suitable for those patients whose conditions do not allow them to perform tasks. Interestingly, the development of new techniques in the field of fMRI, i.e. the resting-state fMRI (rs-fMRI), provide signals that may give insights into the brain function in the absence of stimuli, when the subject is at rest, and is not required to perform any task. This has emphasized the need to recover the underlying neural activations from fMRI signals in the absence of an experimental paradigm. Therefore new techniques that are able to uncover the information hidden in the signal, avoiding the need for additional priors on the activations, are necessary.

Contributions and List of Publications

In this thesis we present two main contributions with a major goal: recovering the brain activation starting from the acquired fMRI data without a priori information on the timing and the duration of the underlying activation.

First Contribution. fMRI Deconvolution via Temporal Regularization using a LASSO model and the LARS algorithm

The fMRI inverse problem is in general ill-posed and the acquisitions and the forward operator, in this case the hemodynamic response function (HRF), are not sufficient to recover a unique solution. The solution is indeed highly sensitive to small perturbations in the data. In the first contribution of this thesis, we solve the inverse problem thus recovering voxel-wise brain functional activations from the blood-oxygen-level-dependent (BOLD) signal. To do this, we propose a temporal regularized deconvolution technique which uses an exponential operator, whose shape and performance can be adjusted by tuning a parameter α, into a least absolute shrinkage and selection operator (LASSO) model. We solved the problem via the Least-Angle Regression (LARS) algorithm that computes the entire solution path for all possible lambdas. Therefore we had just to reasonably choose the optimal regularization parameter, hence the optimal solution, among all those outputted by the algorithm. This approach was tested and validated both on phantom and real data. Compared with the temporal regularized deconvolution implemented in the Total Activation (TA) approach proposed by Karahano glu and colleagues [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF] we avoided the need of defining a priori the regularization parameter and we reduced the computation time.

Second Contribution. A Paradigm Free Regularization Approach to Recover Brain Activation from fMRI Data

A natural extension of the first contribution is to include also the spatial prior into the regularization problem to achieve a solution that is more accurate. In fact, fMRI images have a high spatial resolution, and the information recovered between voxels that are neighbours to each other can be exploited.

State-of-the-art techniques, in particular the TA approach [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF], consider the problems of spatial and temporal regularization as decoupled tasks, thus doubling the number of parameters to be set and requiring the solver to alternate between the constraints. In this thesis we propose a paradigm-free regularization algorithm based on partial differential equations named PF-fMRI (Paradigm-Free fMRI) that applies on the 4-D fMRI image, acting simultaneously in the 3-D space and the 1-D time dimensions. The PF-fMRI is based on the idea that large image variations should be preserved as they occur during brain activation, whereas small variations should be smoothed to remove noise. Starting from this principle, using the PF-fMRI allows us to smooth the whole fMRI image with an anisotropic regularization, thus recovering the location of the brain activations in space and their timing and duration. The PF-fMRI is validated both on synthetic and on real task-fMRI data, that provided us with an experimental paradigm as ground truth to be able to assess the quality of the results. We compared the PF-fMRI with state-of-the-art techniques and finally applied it to rs-fMRI data as a proof of concept. Other experiences

List of Pubblications

Software Contributions

The second contribution of this thesis has been entirely implemented in a Python package, called Paradigm-Free Functional MRI (PF-fMRI), that is available in a GitLab repository dedicated to the Computational Brain Connectivity Mapping (CoBCoM ERC AdG project).

In the package we implemented all the steps we describe in Chapter 4. It simply requires the acquired and preprocessed fMRI image and returns as output the regularized 4-D data.

Overview of the Thesis Chapter 1. From Neurons to BOLD-contrast imaging

This chapter introduces the neurophysiological bases of the brain function and structure, starting with a description of neurons, as the fundamental elements of the brain. Starting from a description at a cellular level, we gradually arrived at describing what is the brain anatomy and its function. In this chapter we also introduce the neuroimaging modalities and the MRI technique, starting from its basic principles. In this framework, we provide a focus on fMRI to have a comprehensive knowledge of the signal we are measuring and analyzing in the following chapters of this thesis.

Chapter 2. fMRI Data Analysis

This chapter gives an overview of the state of the art concerning the methods for fMRI data analysis. We starts from a classical approach, the GLM which requires a priori information to recover the brain activations. We then introduce the data-driven approaches as a set of methods that aims at grouping together voxels whose signals have a certain type of similarity into networks or clusters. After that, we introduce a group of state-of-the-art approaches based on deconvolution techniques, that do not require a priori information on the brain activation they uncover. We then discuss the limitations of the current state of the art which delineates the context and the reasons that lead to the new techniques that we propose in this thesis.

Chapter 3. fMRI Deconvolution via Temporal Regularization using a LASSO model and the LARS algorithm

This chapter introduce our first main contribution. We started from giving an overview of regularization approaches, then we introduce our temporal regularized deconvolution approach to recover brain activity. To do this, we used a LASSO model and the LARS algorithm that interestingly gives as output all possible regularization parameters and their associate solutions. To validate our approach we tested it on phantom and real task-fMRI data to have a ground truth to which to refer to assess the results. We compare our approach with a state-of-the-art technique and discuss the results.

As for the phantom data, in this chapter we also describe how we simulate the synthetic fMRI images. The novelty we introduced in the data simulation, that we employed for the validation of the method described in Chapter 3, is that in addition to a random Gaussian noise, we added a block type noise to simulate the fMRI measurements, which takes into accounts head motions and/or false neural activations.

In Appendix A we report the design of the proposed so-called α-filter, that implements the accumulation function we employed in our model.

Chapter 4. A Paradigm Free Regularization Approach to Recover Brain Activation from fMRI Data

In this chapter we describe our second main contribution, where we aimed at exploiting both the temporal and spatial features of the fMRI data structure. To do this, we aimed at uncovering the brain activations treating the entire 4-D fMRI image as a whole. Therefore, starting from the idea that the brain activates in constant blocks we propose an approach that keeps big image variations because they represent neural activations and smooth small variations to reduce the signal degeneration due to noise. This novel approach, the Paradigm Free fMRI (PF-fMRI), acts via the 4-D image structure tensor and uses partial differential equations (PDEs) to iteratively and anisotropically regularize the whole 4-D fMRI image.

To validate our approach, we applied it both on synthetic data and on real task-fMRI data Introduction from 51 subjects. We finally successfully compare our results with several state-of-the-art techniques.

We applied the PF-fMRI on task-fMRI data to have a ground truth to be able to assess the results, even if the final aim is the application to rs-fMRI data, to which we also applied the PF-fMRI as a proof of concept. The proposed approach was also applied in a work that aimed at inferring the information flow in the WM of the brain and at recovering cortical activity using a multi-modal approach based on fMRI, dMRI, and MEG, without a manual selection of the WM connections of interest as reported in Appendix B.

Chapter 5. Concluding Remarks and Open Problems

This last chapter contains a discussion and a conclusion summarizing the main contributions of this dissertation and the improvements achieved compared to the state of the art. We also discuss the current limitations of the proposed methods and the main perspectives for future works and clinical applications.

Part II

Background

Chapter 1

From Neurons to BOLD-contrast imaging

The brain is an astonishing less than 1.5 kg organ that regulates all body functions, defines the essence of the human mind and soul and interprets information from the outside world [START_REF] Kuffler | From Neuron to Brain[END_REF]. Intelligence, creativity, emotion, and memory are a few of the many things governed by the brain. The brain uses electrical signals, to process and examine all the information it receives. Even though these signals are virtually identical in all nerve cells, they do not at all look like what they appear to produce in the practical life of the human being. In fact, by means of these signals, the brain gets information by the human five senses -sight, smell, touch, taste, and hearing -and assembles the messages in a way that has meaning for us, and can store that information in our memory. The brain controls our thoughts, memory and speech, movements, and the vital function of many organs within our body. To understand how this is done, it is therefore an essential task to decode these signals and the way they are organized and processed by the brain.

In this chapter we provide an introduction to the brain, giving a comprehensive overview of its biological composition, its anatomy and, because in this thesis we are interested in the recovery of the brain functional activations, a zoom on the brain function and open challenges. After, we describe the neuroimaging modality we are interested in, the fMRI technique, starting from the basic principles of the MRI until we reach the BOLD contrast imaging. We conclude with a a description of the two main fMRI imaging techniques:

task-and rs-fMRI. The aim of this chapter is to give an understanding of the signal we are measuring and the type of information this signal provides.

Brain Biology, Anatomy and Physiology

In this section we introduce the brain and the terminology used in the discipline of neurology. There are three sections which cover: the cellular structure of the brain, its anatomy, and its function. The contents reported in this section are mainly based on the books referenced here: [START_REF] Kuffler | From Neuron to Brain[END_REF][START_REF] Catani | Atlas of human brain connections[END_REF][START_REF] Ashburner | Human brain function[END_REF].

Brain Cellular Structure

The human brain is one of the largest and most complex organs in our body; it is composed by approximately 86 billion neurons, which constitute the brain's fundamental units, or building blocks. Neurons are nerve cells within the nervous system that transmit information to other nerve cells, thus allowing the exchange and "communication" within each others by means of a combination of electrical and chemical signals [START_REF] Caire | Physiology, synapse[END_REF]. Most neurons have a cell body, an axon and dendrites (Figure 1.1). The cell body is the place where the metabolic activity of a neuron take place and the neurotransmitters are synthesized.

The dendrites are the cell body projections, covered by dendritic spines, with a receptive role. The axon is a nerve fiber which constitutes a projection of the neuron, where the conduction of electrical impulses, i.e. the action potentials, from the body cell to the presynaptic terminal, take place. These electrical signals are propagated along the axons, by means of its coating sheath, known as myelin.

Key concepts for the comprehension of the brain function are the so-called synapses. A synapse is a tiny gap across which neurons transmit their energy from the pre-synaptic neuron to the post-synaptic neuron thus "talking" to each other. The process through which a cell collects and aggregates all the incoming input signals that then transmit to other neurons is called integration. Different types of neurons, which are tailored to the job they perform, constitute the central nervous system, that controls most functions of the body and the mind. Signals from sensory receptors over the body feed along the spinal cord to the brain, and signals are sent from the brain to execute a task, for example, muscles contraction. A failure of these complex system may cause a malfunctioning of the entire human body mechanism, leading to medical diseases. Neuro-degenerative diseases have been impacting the health of human beings because of the increased life expectancy. For example in the Parkinson's disease where there is a deficiency of the neurotransmitter dopamine. Therefore it is crucial to study the brain's function and have a better understanding of it to be able to reduce the impact of neurodegenerative diseases.

Human Brain Anatomy

The brain is composed by three main parts: the cerebrum, the cerebellum and the brainstem.

The cerebrum is the largest and uppermost portion of the brain, accounting for the twothirds of the total weight of the brain and is divided into a right and a left hemisphere. It performs higher functions such as thinking, speech, reasoning, emotions, learning, and fine control of voluntary movements as well as interpreting touch, vision and hearing. In most cases, the left hemisphere is functionally dominant and accounts for language and speech control. The right hemisphere instead is involved in the interpretation of visual and spatial information. The cerebellum is positioned underneath the cerebrum and is responsible of coordination of movements and balance. The brainstem is located between the spinal cord and the rest of the brain, thus playing as a transfer center. It performs basic automatic functions such as controlling breathing, heart rate, temperature of the body, sleep, digestion, sneezing, coughing and swallowing. Going more into depth on the cerebrum: it appears as an ensemble of folds, aslo known as gyri, and grooves between folds, also called sulci. The cerebral cortex consists of a thin strip of brain cells, or gray matter (GM), representing the outermost surface of the cerebral hemispheres, and an inner core of myelinated nerve fibres (axons) constituting the WM. Nerve fibres in the WM connect functional areas of the cerebral cortex. The GM of the cerebral cortex contains sensory, motor and important association areas and its usually is divided into four lobes, delineated by major surface folds [START_REF] Javed | Neuroanatomy, cerebral cortex[END_REF]. The lobes are: the frontal lobe, the parietal lobe, the temporal lobe and the occipital lobe; the limbic lobe can be considered as a fifth lobe. The frontal lobe is the most anterior and the largest lobe of the human brain; the parietal lobe is located between the frontal and the occipital lobes, above the temporal one. The limbic lobe is constituted by structures on the medial surface which surround the corpus callosum. In order to describe the brain lobes, it is necessary to define conventions and to point out that the brain surface can be observed from different points of view as illustrated in that extends from the central sulcus and reaches the cingulate sulcus. The parieto-occipital sulcus delimits the separation of the parietal from the occipital lobe, from point 5 to point 12. Finally, from point 6 to point 8, an arbitrary line running from the preoccipital notch to the lower tip of the parieto-occipital sulcus, separates the occipital and the temporal lobes.

As for the two hemispheres, they are connected via a thick band of WM, i.e. the corpus callosum, which allows integration of sensory input and functional responses from both sides of the body. Other cerebral structures include the hypothalamus, that regulates the metabolism and preserves the homeostasis, and the thalamus, that constitutes the principal sensory relay center. The brain structures are surrounded by ventricles, which are spaces filled with the cerebrospinal fluid (CSF). The role of the CSF is to supply the brain cells with nutrients and to provide mechanical support and absorb eventual shocks. Finally, the brain is surrounded by a layer of tissue called the meninges and above these the skull, that encloses and protects the brain from injuries.

Human Brain Function

The concept of brain function generally refers to the brain's ability to perform a cognitive or physiological task [START_REF] Papo | Gauging functional brain activity: from distinguishability to accessibility[END_REF]. A brain functionality is achieved by the cooperation of multiple neurons adjacent to each other as well as by the association of neurons which are instead segregated in space and located in different brain regions. The brain function is increasingly regarded as the result of widely interconnected neurons arranged both laterally and hierarchically within the cerebral cortex and deep brain nuclei. This architecture essentially allows continuous input, integration and output of several multi-modal sensory and physiological flows simultaneously.

As it emerged already from the previous lines, there exists in fact two main concepts of brain function: the specialization and the integration. From a macro-scale point of view, this is what is defined as functional specialization. In fact, the frontal lobe contains control centers for motor functions, problem solving, speech and judgment as well as personality, behaviour and emotions, concentration, intelligence and self-awareness [START_REF] Pirau | Frontal lobe syndrome[END_REF][START_REF] Javed | Neuroanatomy, cerebral cortex[END_REF]. The parietal lobe is involved in somatic sensations, such as touch, pain and temperature, as well as spatial and visual perceptions, interpretations and processing of language, visual, hearing, motor, sensory and memory stimuli [START_REF] Ghoneim | Balint syndrome (chronic visualspatial disorder) presenting without known cause[END_REF]. The temporal lobe is responsible for auditory reception and processing, and memory [START_REF] Javed | Neuroanatomy, wernicke area[END_REF]. The occipital lobe controls the visual acquisitions and processing and the limbic lobe governs the smell, taste, and emotions.

The brain specialization has been found already two centuries ago, when the physician Gall claimed that the brain is the seat of the mind and that the mind is made of different mental faculties residing in different specific brain regions [START_REF] Flourens | Recherches experimentales sur les propretes et les fonctions fu systeme nerveux[END_REF][START_REF] Lashley | Brain mechanisms and intelligence: A quantitative study of injuries to the brain[END_REF]. After that, in the second half of the 19th century, several scientists, e.g. Brodmann and Broca, worked in the localization of brain regions related to precise functions and by the early 20th century, it was commonly proved and accepted by the scientific community that sensory and motor functions are located in specialized cerebral areas [START_REF] Kanwisher | Functional specificity in the human brain: a window into the functional architecture of the mind[END_REF][START_REF] Mahon | The specialization of function: Cognitive and neural perspectives[END_REF]. Nonetheless, evidences showed that the brain is much more complex than this and indeed there was already a debate on the fact that the brain is not only specifically organized. Among those who were sustaining this idea there was for example Lashley, who was not able to find located cortical areas responsible for memory and cognition [START_REF] Lashley | Functional determinants of cerebral localization[END_REF][START_REF] Lashley | In search of the engram[END_REF][START_REF] Kanwisher | Functional specificity in the human brain: a window into the functional architecture of the mind[END_REF]. Besides the functionally specialized regions responsible for sensory and motor functions, there are evidences that the cortical regions which activates for some functions appear also when performing other high-level cognitive functions integrating with other brain regions, for example for language [START_REF] Gazzaniga | The cognitive neurosciences[END_REF].

Modern neuroscientists challenge the exclusive idea of the brain as being just functionally specialized and recognize that a single region participates in multiple and diverse functions [START_REF] Huettel | Functional magnetic resonance imaging[END_REF]. When we talk about functional integration we refer to the brain as a whole, meaning that different brain regions cooperate together and process information to achieve a certain function. Brain areas, i.e. neurons, are indeed connected at different degrees via the WM fiber bundles, that allows the propagation of the information from a region to another.

The comprehension of the brain function gives rise to two complementary objectives: on one hand, understanding how the brain anatomical structures and dynamics control its functions; on the other hand, comprehend how actions and behaviours produce functional brain subdivisions [START_REF] Papo | Gauging functional brain activity: from distinguishability to accessibility[END_REF]. We can then assert that the brain is difficult to study not only because of its inherent complexity -the billions of neurons, the hundreds or thousands of types of neurons, the trillions of connections -but also because it works at a number of different scales, both in the physical sense and in the time domain. To capture and understand the brain's electrical activity at these scales, it is crucial to collect information about the brain as a whole and nowadays no single technology is extensively enough. In fact researchers are limited in the sort of approaches they can use to study human brain activity, because they suffer from a lack of detail. Nevertheless, even if we are still far away from the full understanding of how the brain implements given functions, the study of brain dynamics using different neuroimaging data analysis techniques give hints on how these functions may be explained.

Introduction to NeuroImaging

Comprehending the nature of the human brain, the biological basis of learning, memory, behavior, perception and consciousness has been described by Eric Kandel as the "ultimate challenge" of biological sciences [START_REF] Kandel | Principles of neural science[END_REF]. In recent decades there has been a continuous development of neuroimaging techniques, increasingly used in scientific research and in clinical practice. Standard neuroimaging techniques provide non-invasive access not only to the anatomy of the human brain but also to its physiology, its functional architecture and its dynamics. Functional neuroimaging techniques provide an excellent opportunity for investigating the human brain in vivo [START_REF] Grady | Cognitive neuroscience of aging[END_REF]. In particular, modalities such as positron emission tomography (PET), single photon emission computed tomography (SPECT), functional magnetic resonance imaging (fMRI), and magneto-encephalography (M/EEG) led to a new era in the study of the brain functioning.

EEG measures brain electrical activity via electrodes set on the scalp; MEG records the magnetic field by means of sensors placed above the head. Both EEG and MEG have high temporal resolution, on the order of milliseconds; nonetheless their spatial resolution, on the order of centimeters, is low. Over EEG, MEG has the advantage of showing a better signal localisation, however it is expensive and it is limited in the detection of deeper brain structures events. PET gives measures of metabolic processes, while fMRI measures increased neural activities reflected in the changes of blood oxygenation. fMRI has the advantage of having a high spatial resolution, on the order of a few millimeters, even if it has low temporal resolution, between hundreds of milliseconds and seconds. Compared with older techniques, such as the single-or multi-units recordings used to investigate the neurons' physiology, the activity recorded with M/EEG, PET and fMRI, also known as functional neuroimaging modalities, have the advantage of being non-invasive, hence exploitable for the investigation of brain functions in the living human beings. Furthermore, they allow a whole-brain reconstruction, rather than the image of a small brain portion and this is crucial when the aim is investigating the functional integration of different brain areas. However, because these techniques outputs just indirect measures of brain activity -e.g. the ratio of the blood oxygenation, firing rates or membrane potentials -there is a gap between what is actually measured and how these measures reflect the brain functioning architecture and mechanism. In other words what is observed is an example of dynamics, from which functional brain activity should be inferred. This is the context where neuroscientists, computer scientists, mathematicians, psychologists as well as neurologists are needed, by developing mathematical models, experimental designs and interpreting data to make evidences. Besides, distinguishing between bare dynamic and authentic function is an extremely non-trivial task, but of fundamental importance when comparing experimental observations and interpreting their meaning [START_REF] Papo | Gauging functional brain activity: from distinguishability to accessibility[END_REF]. For this reason, it is important to provide more and more precise techniques for the recovery of the brain functional activations.

Fundamentals of MRI

To Regarding brain imaging, MRI was then mainly used to diagnose neurological disorders such as atrophies related to epilepsy, cancerous tumors, or Alzheimer disease [START_REF] Brant-Zawadzki | Mr imaging of the aging brain: patchy white-matter lesions and dementia[END_REF]. Only in the early 90's, MRI started to be used to study brain function [START_REF] Kwong | Dynamic magnetic resonance imaging of human brain activity during primary sensory stimulation[END_REF].

Magnetization

Most of the atoms constituting the human body, such as the human brain, are characterized by having nuclear magnetic resonance properties since they behave like very small magnets and rotate -spin -around their axes. If a constant magnetic field is applied to these atoms, they start precessing around axes parallel to the magnetic field thus maintaining a gyroscopic motion with different angles and frequencies. When they are magnetized, after a few seconds, protons reach equilibrium and spin around their own axis, all aligned along the same direction and around the principal axis of the magnetic field, also known as B 0 .

They can assume a state parallel or anti-parallel to the magnetic field. These states are also referred to as a low energy state, in the case of the parallel one and a high energy state for the anti-parallel. The difference of the numbers of atoms in parallel and anti-parallel state is usually summarized by a vector that represents the sum of the magnetic momenta of all atoms, i.e. the net magnetization. This is a steady state, with static condition where there are no changes. The number of protons aligned to B 0 depends on the strength of the magnetic field, which is measured in Tesla (T); this translate into the fact that MRI scanners with higher strength provides images with better signal-to-noise ratios (SNRs).

Excitation

When a small magnetic field, perpendicular to B 0 is introduced in the system, this causes spins to move away from their principal direction. This process is known as excitation, and it consists in using a radio-frequency (RF) pulse at the resonance frequency of the magnetized atoms that induces the magnetization to oscillate between the value they had at equilibrium and its opposite one. In practice the spinning protons are pushed out of the alignment of the static magnetic field B 0 . For example, when a 90 degree RF pulse is emitted, the spins are pushed away from B 0 of 90 degrees.

Relaxation

When this RF pulse is interrupted, spins will follow the already described precession to reach again the equilibrium and move from a high energy state to a low energy one and spin along the principal axis of B 0 .

There are two aspects of this precession that we can measure: one is the longitudinal relaxation and the other is the transverse relaxation (Figure 1.4). The longitudinal relaxation, also known as T 1 relaxation, is the process by which the net magnetization returns to its original alignment parallel to B 0 . So, from Figure 1.4, if we consider B 0 aligned to the z-axis, starting at 0 from the xy plane, T 1 is the time for the z-component of the magnetization M to reach the 63% of its maximum value. The transverse relaxation is the process by which the transverse component of the magnetization, M xy , decays. The transverse relaxation is also referred to as T 2 , that is the time required for the transverse magnetization to fall to approximately 37% of its initial value. The T 2 physically describes the decay of the atoms' phase coherence due to magnetic interactions between the nuclei of the system. If local inhomogeneities in the magnetic field are considered in the phase decay, this leads to the T * 2 relaxation. T * 2 relaxation refers to an exponential decrease in M xy following the initial excitation pulse as a function of time T * 2 . So, at the beginning this is a large vector, while it decreases as the spin relaxes and realigns to the static field B 0 , essentially inversely to T 1 .

These relaxations, that comprehensively can be defined as the time that the spins take to recover after being disturbed from equilibrium by a RF pulse, generate a current that can be measured by means of a receiver coil positioned in the MRI scanner. These relaxations are the most commonly used contrasts and represent the MR signal. According to their compositions, different types of biological tissues, and in particular in the case of our interest WM, GM and CSF, have very different but consistent T 1 and T 2 relaxation times, also when varying the strength of the magnetic field, from 1.5 T to 3 T for example. In Table 1.1 the typical values of T 1 and T 2 relaxation times for brain tissues are reported.

The frequency that is needed to excite the protons is known as Larmor frequency .

ω 0 = γB 0 (1.1)
where ω 0 is the Larmor frequency in MHz, γ is the gyromagnetic ratio in MHz/T and B 0 is the strength of the static magnetic field in T. ω 0 is hence defined by the gyromagnetic ratio which is inherent to the considered tissue, since different particles have specific and well known gyromagnetic ratios, and the strength of the magnetic field. In the case of MRI, because the human body is mostly composed by water, the targets are the hydrogen nuclei, which have γ = 42.58.

In order to introduce spatial specificity and therefore locate and image and the different parts of a tissue, each portion of tissue has to experience a unique magnetic field. To do this, linear gradients in the magnetic field are applied in the three orthogonal directions: x, y and z. By choosing the strength of the gradients properly, a one-to-one correspondence can be defined between a spatial volume and the MR signal it generates. An example of T1-weighted and T2-weighted MRI images are reported in Figure 1.5.

The Echo Planar Imaging (EPI) Sequence

The most common imaging sequence for MRI acquisitions, and in particular fMRI, is the echo planar imaging (EPI), a fast method introduced by Mansfield in 1977 [START_REF] Mansfield | Multi-planar image formation using nmr spin echoes[END_REF] that allows to acquire plane-by-plane and in few seconds the whole raw data representing a brain volume, namely the Fourier transform in the k-space [START_REF] Liang | Principles of magnetic resonance imaging: a signal processing perspective[END_REF]. The time required for recording all the lines of the k-space corresponds to the repetition time (TR). In single-shot EPI, all the spatial-encoding data of an image can be obtained after a single RF excitation. As for the conventional spin-echo sequence (SE), a SE EPI sequence starts with 90°and 180°RF pulses.

The initial k-space position is set by means of a G x and a G y gradient pulse and by the phase inversion of the 180°RF pulse. Then, the high values of the frequency encoding gradient T1-weighted T2-weighted The time required for the acquisition of the whole plane (k x , k y ) after each RF pulse is about 100 ms. This technique allows to freeze the dynamics of physiological processes, which contribute to motion artifacts in slower MRI methods, leading to a significantly higher temporal stability.

Basic Principles of Functional MRI

fMRI is a widespread and powerful imaging technique for investigating the brain's functionality in vivo. It provides metabolic information to be integrated to the anatomical information given by the conventional MRI technique. In the XIX century, it was noticed that there was a mechanism through which the blood supply in any cerebral tissue area varies in accordance with chemical changes related to a functional activation [START_REF] Roy | On the regulation of the blood-supply of the brain[END_REF]. The principle behind this mechanism is nowadays successfully used in several neuroimaging modalities, including fMRI, that are based on hemodynamic responses to neuronal activity [START_REF] Logothetis | Neurophysiological investigation of the basis of the fmri signal[END_REF]. fMRI is an imaging technique able to record the effects of altered cerebral electric activity and it can be made sensitive to regional blood perfusion changes, to blood volume and to blood oxygenation as consequence of the neuronal activity [START_REF] Matthews | Functional magnetic resonance imaging[END_REF].

The Blood Oxygenation Level Dependent (BOLD) Contrast Imaging

BOLD contrast imaging is one of the most commonly used methods to obtain information about brain functionality. As mentioned before, the neural activity is not directly measured by the BOLD signal, which instead relies on a cascade of physiological events linking synaptic activity to the generation of the MRI signal [START_REF] Iannetti | Bold functional mri in disease and pharmacological studies: room for improvement?[END_REF]. The BOLD contrast is based on MR images sensitive to oxygenation changes in the hemoglobin [START_REF] Ogawa | Oxygenation-sensitive contrast in magnetic resonance image of rodent brain at high magnetic fields[END_REF]. If the hemoglobin is bound with the oxygen (oxyhemoglobin), it is diamagnetic with negligible effects on the magnetic field; conversely, the deoxygenated hemoglobin (deoxyhemoglobin) is paramagnetic [START_REF] Amaro | Study design in fmri: basic principles[END_REF], introducing a magnetic distortion on the scale of capillaries and venules, and a inhomogeneity of the magnetic field at the voxel scale, leading to a decreased BOLD signal [START_REF] Ogawa | Oxygenation-sensitive contrast in magnetic resonance image of rodent brain at high magnetic fields[END_REF][START_REF] Turner | Echo-planar time course mri of cat brain oxygenation changes[END_REF] (Figure 1.7). One expects that the deoxyhemoglobin ratio would increase and that the BOLD signal would decrease due to oxygen consumption following neuronal activity; however, the vascular effect causes overcompensation of oxygen, which leads to an increase in the BOLD signal [START_REF] Malonek | Interactions between electrical activity and cortical microcirculation revealed by imaging spectroscopy: implications for functional brain mapping[END_REF]. The above effect depends on the degree of regional activation in the GM: due to auto regulation mechanisms, the regional blood flow increases according to the level of neural activity; the augmented oxygenated blood largely exceeds the local metabolic uptake thus augmenting the saturation level of capillaries and venous compartment, as well; the decreased inhomogeneity ultimately results in increased BOLD signal by T2* weighted sequences. It is interesting to notice that the brain weights 2% of total body weight but it receives about 15% of the cardiac output [START_REF] Javed | Neuroanatomy, cerebral cortex[END_REF].

In fMRI, the oxy-and the deoxyhemoglobin are used as endogenous contrast agents.

According to the percentages of each within each voxel, representing a small part of the brain, it is implied how the MR signal will behave in BOLD image: areas where the oxyhemoglobin is highly concentrated give a higher signal -a brighter image -than areas where the oxyhemoglobin is poorly concentrated [START_REF] Amaro | Study design in fmri: basic principles[END_REF].

The time course of the BOLD signal after a short neural activation is known as hemodynamic response function (HRF, 1.8) and is approximately modeled as the impulse response of a linear dynamic system [START_REF] Glover | Deconvolution of impulse response in event-related bold fmri1[END_REF]. The initial part of the HRF reveals a transient increase in deoxyhemoglobin concentration, highlighted by an initial dip [START_REF] Yacoub | Investigation of the initial dip in fmri at 7 tesla[END_REF]. Then, an increase in the oxy/deoxyhemoglobin ratio occurs, leading to a consequent increase in the MR signal, proportionally to the underlying neural activity [START_REF] Logothetis | On the nature of the bold fmri contrast mechanism[END_REF]. This signal increase can reach a plateau if the stimulus is maintained for a sufficient time [START_REF] Buxton | Modeling the hemodynamic response to brain activation[END_REF] and it is of a few percentage parts, typically in the order of 0.5 -3% at 1.5 T. After the cessation of the stimulus, the MR signal returns to the baseline, and eventually underpasses it. This last event is named the "undershoot effect" [START_REF] Buxton | Dynamics of blood flow and oxygenation changes during brain activation: the balloon model[END_REF]. An example of a fMRI acquisition is showed in Figure 1.9.

Sources of Noise in fMRI Images

The BOLD variations observed in fMRI images are not always due to pure neural activation, because fMRI images are often very noisy images. There exist multiple sources of noise.

fMRI images may be corrupted by fluctuations registered as a consequence of physiological noise due to cardiac and respiratory cycles. Head motions during the scans may lead to changes in voxel intensity much greater than the BOLD activation response [START_REF] Clare | Functional mri: methods and applications[END_REF]. Thermal noise due to eddy currents and heating provoked by the motion of electrons in the subject and in the scanner [START_REF] Friston | Event-related fmri: characterizing differential responses[END_REF][START_REF] Friston | Movementrelated effects in fmri time-series[END_REF] is another source of noise. Furthermore there is a variability of the evoked hemodynamic response across subjects and within the same subject acquired in different scanning sessions [START_REF] Aguirre | The variability of human, bold hemodynamic responses[END_REF]. If the motion artifacts and low frequency drifts are removed via an appropriate data preprocessing, the residual errors can be assumed as white noise [START_REF] Lund | Non-white noise in fmri: does modelling have an impact?[END_REF][START_REF] Karahanoglu | Sparsity-promoting spatiotemporal regularization for data mining in functional magnetic resonance imaging[END_REF].

Standard Minimal Preprocessing of fMRI Data

A broad and variate set of preprocessing pipelines are provided in the literature and implemented in multiple available toolboxes. Nonetheless, the standard preprocessing steps of fMRI data included in a standard pipeline are enumerated below.

• Realignment of the functional volumes of a subject for head motion correction.

• Removal of non-brain tissues from images of the whole head. • Spatial smoothing via a Gaussian kernel to reduce the spatial variability and increase the SNR.

• High-pass temporal filtering to remove low frequency drifts.

• Registration of the functional images to the anatomical image.

• Coregistration of functional images to a standard space, e.g. the Montreal Neurological Institute (MNI) to allow comparison between subjects and group studies.

fMRI Techniques: Task and Resting-State fMRI

fMRI studies can be divided according to two acquisition approaches: task-dependent and task-independent. In standard task-dependent fMRI acquisition, fMRI aims at examining the brain function by detecting metabolic/hemodynamic responses consequent to brain activity and, consequently, at identifying areas of increased or decreased neuronal activity [START_REF] Logothetis | The underpinnings of the bold functional magnetic resonance imaging signal[END_REF][START_REF] Raichle | Brain work and brain imaging[END_REF] while subjects perform motor, sensory, cognitive or emotion-provoking tasks [START_REF] Huettel | Functional magnetic resonance imaging[END_REF].

An experimental task of interest is presented alternately with a control task and the BOLD signal during the experimental task is compared to the one measured during the control task.

However, the difference between baseline and task-related activation, in terms of changes observed in the BOLD contrast, accounts for about 1-5% of the total signal [START_REF] Damoiseaux | Greater than the sum of its parts: a review of studies combining structural connectivity and resting-state functional connectivity[END_REF] and statistics over repeated activations (either block-design or event-related) are necessary to provide response images as statistical parametric maps. On the other hand, the task-independent approach allows to study the brain function by analyzing the signal fluctuations depending on the level of blood oxygenation (i.e., BOLD again) when the subject does not perform any particular task. This practice is also known as Resting State fMRI (rs-fMRI) and poses challenging data analysis problems since no external reference is given and features attain to the internal functional connectivity are blindly extracted.

This practice is based on the idea that the regions typically modulated by stimuli or tasks (e.g., verbal, physical, visual, auditory, etc.) display vigorous and persistent functional activity, even when patient is asked to keep their eyes closed and to relax [START_REF] Buckner | The brain's default network: anatomy, function, and relevance to disease[END_REF]. This activity is detected as spontaneous low-frequency (< 0.1 Hz) BOLD signal fluctuations and inter-regional correlations of these fluctuations can be estimated as measures of functional connectivity [START_REF] Friston | Functional connectivity: the principal-component analysis of large (pet) data sets[END_REF]. Functional connectivity is defined as the temporal dependency between spatially remote neurophysiological events [START_REF] Friston | Functional and effective connectivity: a review[END_REF] and is used to describe the relationship between neural activation patterns related to anatomically separate brain regions, reflecting their level of communication [START_REF] Van Den | Exploring the brain network: a review on resting-state fmri functional connectivity[END_REF]. The coherent activity of functionally related brain areas can be captured as temporally correlated fluctuations in BOLD signal during rs-fMRI acquisitions.

The first rs-fMRI study, conducted by Biswal and colleagues in 1995 [START_REF] Biswal | Functional connectivity in the motor cortex of resting human brain using echo-planar mri[END_REF], examined the correlation between the time course of the signal in a seed region of interest (ROI), placed in the motor area, with the signal of all other brain voxels. The result showed that, at rest, the left and right hemispheres of the primary motor network are not silent, but show a high temporal correlation between their rs-fMRI BOLD time courses. This result suggested that at rest there is an ongoing information exchange and an ongoing functional connectivity between these regions [START_REF] Biswal | Simultaneous assessment of flow and bold signals in resting-state functional connectivity maps[END_REF][START_REF] Cordes | Mapping functionally related regions of brain with functional connectivity mr imaging[END_REF][START_REF] Greicius | Functional connectivity in the resting brain: a network analysis of the default mode hypothesis[END_REF][START_REF] Lowe | Correlations in low-frequency bold fluctuations reflect cortico-cortical connections[END_REF]. When fMRI studies started to examine the possibility of measuring functional connectivity between brain regions as the level of co-activation of spontaneous fMRI time series recorded during rest [START_REF] Biswal | Simultaneous assessment of flow and bold signals in resting-state functional connectivity maps[END_REF][START_REF] Greicius | Functional connectivity in the resting brain: a network analysis of the default mode hypothesis[END_REF][START_REF] Lowe | Correlations in low-frequency bold fluctuations reflect cortico-cortical connections[END_REF], it was found that, at rest, the brain is organized into networks related to the functions of vision, hearing, motor planning and attention [START_REF] Calhoun | A method for making group inferences from functional mri data using independent component analysis[END_REF][START_REF] Beckmann | Investigations into restingstate connectivity using independent component analysis[END_REF][START_REF] Damoiseaux | Consistent resting-state networks across healthy subjects[END_REF][START_REF] Smith | Correspondence of the brain's functional architecture during activation and rest[END_REF]. These networks showed surprisingly consistent patterns of activation, though not identical, very similar to those of task-induced activation and deactivation [START_REF] Calhoun | Modulation of temporally coherent brain networks estimated using ica at rest and during cognitive tasks[END_REF][START_REF] Harrison | Consistency and functional specialization in the default mode brain network[END_REF][START_REF] Smith | Correspondence of the brain's functional architecture during activation and rest[END_REF] and are called Resting State Networks (RSNs). The RSNs showed high reproducibility [START_REF] Damoiseaux | Consistent resting-state networks across healthy subjects[END_REF] and moderate/high reliability [START_REF] Franco | Interrater and intermethod reliability of default mode network selection[END_REF][START_REF] Shehzad | The resting brain: unconstrained yet reliable[END_REF][START_REF] Zuo | The oscillating brain: complex and reliable[END_REF]. This suggests that rs-fMRI allows an exhaustive study of the intrinsic functional architecture of the human brain connectome [START_REF] Biswal | Toward discovery science of human brain function[END_REF]. The analysis of the rs-fMRI data can be used to study a wide range of neurological and psychiatric disorders, for example Alzheimer's disease [START_REF] Damoiseaux | Functional connectivity tracks clinical deterioration in alzheimer's disease[END_REF], dementia with Lewy bodies [START_REF] Galvin | Resting bold fmri differentiates dementia with lewy bodies vs alzheimer disease[END_REF], frontotemporal dementia [START_REF] Whitwell | Altered functional connectivity in asymptomatic mapt subjects: a comparison to bvftd[END_REF], epilepsy [START_REF] Luo | Disrupted functional brain connectivity in partial epilepsy: a resting-state fmri study[END_REF], Parkinson's disease [START_REF] Kwak | Altered resting state cortico-striatal connectivity in mild to moderate stage parkinson's disease[END_REF], stroke [START_REF] Park | Longitudinal changes of resting-state functional connectivity during motor recovery after stroke[END_REF], depression [START_REF] Sheline | Resting-state functional mri in depression unmasks increased connectivity between networks via the dorsal nexus[END_REF], schizophrenia [START_REF] Cole | Variable global dysconnectivity and individual differences in schizophrenia[END_REF], the obsessive-compulsive disorder [START_REF] Fitzgerald | Developmental alterations of frontal-striatal-thalamic connectivity in obsessive-compulsive disorder[END_REF], attention-deficit/hyperactivity disorder [START_REF] Fair | Atypical default network connectivity in youth with attention-deficit/hyperactivity disorder[END_REF], Tourette's syndrome [START_REF] Church | Control networks in paediatric tourette syndrome show immature and anomalous patterns of functional connectivity[END_REF], autism [START_REF] Kennedy | The intrinsic functional organization of the brain is altered in autism[END_REF], and others.

More recently, dynamics of rs-fMRI provided a new insights onto the organizational principles of brain function: from studies of dynamic functional connectivity [START_REF] Preti | The dynamic functional connectome: State-of-the-art and perspectives[END_REF] to the study of different spatio-temporal structures of functional components that dynamically assemble the RSNs [START_REF] Karahano | Transient brain activity disentangles fmri resting-state dynamics in terms of spatially and temporally overlapping networks[END_REF].

The task-independent approach is very useful especially when subjects are children or patients with disabling pathologies. Compared to task-fMRI, rs-fMRI also allows to study the activity of more brain networks and more functions in a single acquisition. The disadvantages of the rs-fMRI, compared to the standard task-activation fMRI, are the absence of a priori hypothesis on the underlying activation and of externally triggered temporal references.

Exploratory and methodological research is still active and no standard method for the study of the brain activity has been established, so far. There is a broad range of methods with their respective pros and cons, fostering efforts to optimize the currently available methods.

Increasing evidences suggest that coherent intrinsic brain activity is an important feature of healthy brain functioning and a translation to clinical neurology would have a deep impact [START_REF] Fox | Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging[END_REF].

Conclusion

In this chapter we presented the brain giving its cellular composition, its anatomy, and its function. We pointed out the complexity of the brain that poses open challenges for the understanding of the brain functionality because the brain regions are specialized as well as integrated at different scales. In this chapter we also introduced the neurogimaging and the principles of MRI as well as fMRI, presenting the BOLD signal and fMRI data. We finally introduced the two typologies of fMRI experiments: task and rs-fMRI and we highlight the necessity of new methodologies to investigate the whole-brain fMRI data, especially in absence of a priori knowledge on the experimental paradigm.

Part III

State of the Art Chapter 2 fMRI Data Analysis

There exist multiple goals in the analysis of fMRI data; among them, there are the localization of cortical areas activated by a certain task, the determination of brain functional networks -i.e. group of brain regions which activates or deactivates together -and the identification of biomarkers for a better understanding of neurological and psychological diseases. The choice of the method for the analysis of fMRI data depends strictly on the goal of the study. There are two main categories of approaches: (i) confirmatory approaches, i.e. the GLM, that aims at testing a hypothesis on the data from the experimental paradigm, and

(ii) exploratory or data-driven approaches that are intended to uncover the underlying and unknown data structures. In the next section we discuss the state of the art concerning both the confirmatory and exploratory approaches currently adopted for fMRI studies. We will start by explaining the signal recovery methods, such as the general linear model that requires a priori knowledge on the experimental paradigm and the deconvolution approaches that instead does not. After that we will illustrate and discuss the connectivityoriented approaches, such as the voxel/seed-based mothod, the ICA and the clustering methods.

The General Linear Model

The most widespread approach, for the characterization and localization of brain areas which are activated during a certain task, relies on approaches based on the GLM, adapted by Friston and colleagues in 1998 in the context of fMRI data analyses [START_REF] Friston | Event-related fmri: characterizing differential responses[END_REF].

The GLM approach treats the data as a linear combination of several model predictors, or regressors, defined as explanatory variables corresponding to some experimental effects [START_REF] Flandin | Bayesian fmri data analysis with sparse spatial basis function priors[END_REF], and the error, which is given by the acquired noise and artifacts. The shape of the model functions are known and given by the experimental paradigm while their amplitudes are unknown, therefore they have to be estimated. A GLM needs a priori assumption on the task parameters and timing of events, as well as assumptions about neural and hemodynamic responses.

In practice, the GLM aims at modeling at each voxel the observed time series Y as a linear combination of explanatory variables plus an error term. Supposing that we have a fMRI time series observed in a voxel v at time t

y v t = x v t1 β 1 + ... + x v tl β l + ... + x v tL β L + v t = ∑ i x v ti β i + v t (2.1)
where x tl and β l are the explanatory variables at time t and their corresponding weights, respectively, for l = 1, . . . , L. For t = 1, . . . , N time points, the above equation can be rewritten as follows:

          y v 1 . . . y v t . . . y v N           =           x v 11 . . . x v 1l . . . x v 1L . . . . . . . . . . . . . . . x v t1 . . . x v tl . . . x v tL . . . . . . . . . . . . . . . x v N1 . . . x v Nl . . . x v NL                     β 1 . . . β l . . . β L           +           v 1 . . . v t . . . v N          
In matrix notation this corresponds to the following linear model (Figure 2.1)

y v = X v β v + v (2.2)
where for the voxel v, y v is a the column of observations of lenght N, β is a column vector of lenght L representing the weights and the vector , of lenght N, the noise. X v is a N × L matrix of regressors: the so-called design matrix. The design matrix can contain covariates, such as age or subject motion [START_REF] Friston | Movementrelated effects in fmri time-series[END_REF] or systemic drifts [START_REF] Lund | Non-white noise in fmri: does modelling have an impact?[END_REF] or cardiac and respiratory effects and temporal derivatives of the HRF, as well as variables that reflect the experimental design.

The GLM approach also include the HRF in the method, meaning that the task paradigm is convolved with the HRF and then included in the design matrix. The vector β contains the unknown weights, or parameters. Only some of the parameters in the vector β will be of interest, meaning that they are the effect of some condition or task. Typically, in fMRI these conditions consists of repeated blocks of different tasks over time. The others, which are referred to as confounding factors will be of no interest for the data because due to nonexperimental sources of variability.

Solving the problem corresponds to estimate the parameter weights β that minimize the error . This is done via generalized least squares as follows

β = argmin β (y -Xβ) T Σ -1 (y -Xβ) (2.3) where β = (X T Σ -1 X) -1 X T Σ -1 y (2.4)
that necessitates an estimation of the unknown covariance matrix of the errors Σ. Assuming an identity matrix for the covariance, Σ = I, the solution of (2.

3) is β = (X T X) -1 X T y (2.5)
Once the weights are estimated, t-tests may be used to verify the hypothesis and in fact, the GLM has been proposed in association to the generation of Statistical Parametric Maps (SPMs) [START_REF] Friston | Event-related fmri: characterizing differential responses[END_REF]. SPMs are defined as statistical processes which generates spatial maps about regional specific effects of a task.

It is worth again to emphasize that this approach can be used for task-fMRI data, where the stimulus paradigm is given. On the other hand, the GLM approach is not applicable in the analysis of rs-fMRI, when no a priori information on the timing of the activation is given, or in the study of neural activations which cannot be modeled with predefined stimulus functions, as for example in the case of epileptic crisis [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF].

Deconvolution of the fMRI BOLD Signal

Deconvolution approaches have been developed to address the problem of studying and uncovering brain activations hidden below fMRI time series at the voxel level.

The base problem of deconvolution, that is widely used in signal and image processing, is to reverse the effect of a convolution, that explain the acquired data. For example, given the convolution

y = h * x (2.6)
where y is the acquired data, x is a signal that we aim at recovering and h is the transfer function of an operator, the deconvolution is usually performed via the Fourier Transform on the acquired signal y and the transfer function h. The problem, in absence of noise, becomes

X = Y/H (2.7)
where Y, X and H are the Fourier Transforms of y, x and h respectively. Finally, by applying the inverse Fourier Trasnform on X the deconvolved signal x will be estimated.

fMRI deconvolution was firstly introduced by Glover in 1999 [START_REF] Glover | Deconvolution of impulse response in event-related bold fmri1[END_REF], who investigated the performance of Wiener deconvolution for deblurring the fMRI response to consecutives finger tapping episodes to reduce image distortions. The above-mentioned approach resulted in smooth recovered activation [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF] and was limited by the fact that it required the independent measurement of the noise spectral density [START_REF] Gitelman | Modeling regional and psychophysiologic interactions in fmri: the importance of hemodynamic deconvolution[END_REF]. Gitelman and colleagues [START_REF] Gitelman | Modeling regional and psychophysiologic interactions in fmri: the importance of hemodynamic deconvolution[END_REF] developed an approach based on linear deconvolution and attempted the explanation of brain responses located in one area, as consequence of the interplay with another area as effects of psycho-physiological interactions.

Furthermore, dynamical filter methods, such as Kalman and Bayesian filtering, and Local Linearization filters have been developed and applied to fMRI [START_REF] Friston | Dem: a variational treatment of dynamic systems[END_REF][START_REF] Havlicek | Dynamic modeling of neuronal responses in fmri using cubature kalman filtering[END_REF][START_REF] Riera | A state-space model of the hemodynamic approach: nonlinear filtering of bold signals[END_REF]. However, because these approaches are based on non-linear models in continuous time, they are limited by the high computational cost compared with linear models, therefore they are convenient just for the analysis of predefined ROIs [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF].

Few years later, researchers started to make assumptions on the underlying signals, from a spatial and a temporal point of view, thus adding priors in convex optimization problems.

L1-norm regularization approaches have been developed to exploit sparse temporal features of the hidden neural activation. This was done by means of the majorization-minimization of a cost function to find an optimal solution to the inverse problem [START_REF] Hernandez-Garcia | Neuronal event detection in fmri time series using iterative deconvolution techniques[END_REF]. In addition, Caballero Gaudes et al. developed firstly a ridge-regression regularization by minimizing both the variance and the power of the residuals [START_REF] Gaudes | Detection and characterization of single-trial fmri bold responses: Paradigm free mapping[END_REF] and secondly a sparse regression [START_REF] Gaudes | Paradigm free mapping with sparse regression automatically detects singletrial functional magnetic resonance imaging blood oxygenation level dependent responses[END_REF] by imposing just few coefficients having significantly non-zero amplitudes, assuming short neuronal activation. Also temporal regularized optimization problems based on wavelets were explored [START_REF] Khalidov | Activelets: Wavelets for sparse representation of hemodynamic responses[END_REF].

Recently, a spatio-temporal regularized deconvolution of fMRI BOLD signal was proposed by Karahano glu and colleagues [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF]. In this work, a variational approach, called Total Activation (TA), has been proposed to reconstruct the signals related to the whole brain activation in a voxel-wise framework by imposing informative priors on the signal of interest. They imposed the following two priors, that led to a spatio-temporal regularized deconvolution problem. Supposing the activation is block-type and is called "activityinducing" signal, in one regularization term they exploited the sparsity of the innovation signal, i.e. the derivative of the activity-inducing signal, in the temporal dimension (Figure 2.2). As a second regularization term, because fMRI data presents different kinds of spatial correlations, they used a mixed norm on the spatial dimension. They indeed supposed spatial coherence inside a brain region, favoring smooth patterns of activity inside the region, though brusque changes were assumed to happen across different brain regions.

Later, Farouj et al. [START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF] proposed another spatio-temporal regularization approach for fMRI data and overcame the necessity of using a predefined brain atlas. The regularization proposed in this work [START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF] is expressed as

x = argmin x 1 2 y -x 2 2 + R T (x) + R S (x) (2.8)
where y is the fMRI noisy data, R T and R S are the temporal and spatial regularization, respectively.

More into details, in the problem expressed in (2.8) Farouj and colleagues [START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF] implemented a generalized total variation (TV) for the temporal regularization and a GM-driven TV-based spatial regularization. The temporal regularization is a total variation norm of a 1-D signal in the temporal dimension, which is defined as the l 1 -norm of the derivative of the abovementioned signal. 2011 [START_REF] Karahanoglu | A signal processing approach to generalized 1-d total variation[END_REF] introduced an additional linear differential operator in the l 1 -norm that accounts for the presence of the linear system modelling the hemodynamic response. The spatial regularization instead is implemented as a TV limited to the voxels that can be activated, i.e. those belonging to the GM, and it does not require a predefined division of brain regions, e.g. via atlases. To solve the optimization problem, they used a generalized forwardbackward proximal (GFB) splitting algorithm [START_REF] Raguet | A generalized forward-backward splitting[END_REF] that provides a weighted average of the two solutions obtained by solving one regularization at the time. In this partiular case, each regularization problem was choosen to be solved using the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. The solution estimated using a GFB algorithm for iteration k would than be expressed as follows:

Algorithm 1 GFB splitting algorithm for solving the problem in 2.8

Input: corrupted fMRI data y, (ω

T , ω S ) ∈ [0, 1] 2 with ω T + ω S = 1 Output: Estimate x 1: for k = 1 : k MAX do 2: x k T = argmin x 1 2 y -x 2 2 + R T (x)
3:

x k S = argmin x 1 2 y -x 2 2 + R S (x)
4:

x k = ω T x k T + ω S x k S 5: end for 6: x = x k MAX
where k is the number of the current iteration and k MAX is the maximum number of iterations. ω T is the weight given to the temporal regularization and x k T is the solution of the temporal regularization for iteration k, whereas ω S is the weight given to the spatial regularization and x k S is the solution of the spatial regularization for iteration k. Solving one regularization at the time also means that the regularization parameters lambdas has to be differently set for both problems. In this case the parameter lambda for the temporal regularization was tuned as the median absolute deviation of fine-scale 3 th order Daubechies wavelet coefficient for each voxel while the lambda for the spatial regularization was empirically set. As for the weights, ω T was set to 0.75 and ω S to 0.25, thus giving a higher weight to the temporal regularization with respect to the spatial one.

This approach was exploited to recover transients in rs-fMRI BOLD activitions (positive spikes) and de-activations (negative spikes). To do this, in the work proposed by Karahano glu and Van De Ville in 2015 [START_REF] Karahano | Transient brain activity disentangles fmri resting-state dynamics in terms of spatially and temporally overlapping networks[END_REF], the recovered innovation signals (Figure 2.2) of the whole brain undergo a temporal clustering that allows to find the so-called innovativedriven co-activation patterns (iCAPs), which are consistent and well localized spatial maps of transients during rest (Figure 2.3). This new approach provides new insights into the study of the brain dynamic organization, revealing the dynamics of the well known RSNs, in terms of networks which segregate and aggregate over time.

Voxel/Seed-based Methods

The primary voxel-based (or seed-based) method was introduced by Biswal et al. (1995) [START_REF] Biswal | Functional connectivity in the motor cortex of resting human brain using echo-planar mri[END_REF]. This method evaluates the similarity (e.g., the temporal cross-correlation) between a reference signal taken as seed and the time series of each brain voxel, or the average of the signals contained in specific ROIs [START_REF] Biswal | Functional connectivity in the motor cortex of resting human brain using echo-planar mri[END_REF][START_REF] Martino | Functional connectivity of human striatum: a resting state fmri study[END_REF][START_REF] Taylor | Two systems of resting state connectivity between the insula and cingulate cortex[END_REF]. The outcome is a map of the brain regions significantly correlated with the selected seed, or a quantitative assessment of the correlation strength within a specific ROI [START_REF] Golestani | Regions of interest for resting-state fmri analysis determined by inter-voxel cross-correlation[END_REF] (Figure 2.4). Typically, the seeds are selected using anatomical images as reference [START_REF] Martino | Functional connectivity of human striatum: a resting state fmri study[END_REF][START_REF] Stein | Functional connectivity in the thalamus and hippocampus studied with functional mr imaging[END_REF][START_REF] Taylor | Two systems of resting state connectivity between the insula and cingulate cortex[END_REF][START_REF] Zhou | Altered resting-state functional connectivity and anatomical connectivity of hippocampus in schizophrenia[END_REF], or they are based on the location of activity during a task [START_REF] Biswal | Functional connectivity in the motor cortex of resting human brain using echo-planar mri[END_REF][START_REF] Jiang | Modulation of functional connectivity during the resting state and the motor task[END_REF][START_REF] Lowe | Resting state sensorimotor functional connectivity in multiple sclerosis inversely correlates with transcallosal motor pathway transverse diffusivity[END_REF][START_REF] Wu | Changes of functional connectivity of the motor network in the resting state in parkinson's disease[END_REF][START_REF] Xiong | Interregional connectivity to primary motor cortex revealed using mri resting state images[END_REF], or on standardized coordinates space [START_REF] Maldjian | An automated method for neuroanatomic and cytoarchitectonic atlas-based interrogation of fmri data sets[END_REF]. However, the volume of anatomical regions can vary between individuals [START_REF] Van Essen | Surface-based and probabilistic atlases of primate cerebral cortex[END_REF],

with aging [START_REF] Terribilli | Age-related gray matter volume changes in the brain during non-elderly adulthood[END_REF] and in presence of neurological disease [START_REF] Aydin | Smaller gray matter volumes in frontal and parietal cortices of solvent abusers correlate with cognitive deficits[END_REF], and functional boundaries of brain regions may not be well defined. Hence, the seed time series can be affected by undesired signal, or functionally relevant voxels may be excluded. Even if this method is a powerful and effective tool for the identification and characterization of the resting state, the obtained networks depend on the arbitrary choice of the seed [START_REF] Cole | Identifying the brain's most globally connected regions[END_REF]. Furthermore, this technique allows only the evaluation of the connectivity between the brain and the seed and considers one seed at a time. However, in absence of an a priori hypothesis, as in rs-fMRI case, it is more advisable to simultaneously detect and characterize brain activations in a comprehensive manner, rather than one at a time.

Data-driven Approaches

Data-driven methods have been proposed to analyse images obtained in resting-state conditions, when no information about the occurrence of the activation is available.

After the first acquisitions of rs-fMRI time courses, which before were considered as background noise [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF], it was shown that there was an ongoing information exchange between cortical regions [START_REF] Biswal | Functional connectivity in the motor cortex of resting human brain using echo-planar mri[END_REF][START_REF] Biswal | Simultaneous assessment of flow and bold signals in resting-state functional connectivity maps[END_REF][START_REF] Greicius | Functional connectivity in the resting brain: a network analysis of the default mode hypothesis[END_REF][START_REF] Cordes | Mapping functionally related regions of brain with functional connectivity mr imaging[END_REF][START_REF] Lowe | Resting state sensorimotor functional connectivity in multiple sclerosis inversely correlates with transcallosal motor pathway transverse diffusivity[END_REF]. Indeed it was found that, at rest, the brain is organized into networks related to the different functions, such as: visual and auditory function, motor planning and attention [START_REF] Beckmann | Investigations into restingstate connectivity using independent component analysis[END_REF][START_REF] Damoiseaux | Consistent resting-state networks across healthy subjects[END_REF][START_REF] Smith | Correspondence of the brain's functional architecture during activation and rest[END_REF]. These networks showed consistent patterns of activation very similar, though not identical, to those of task-induced activation and deactivation [START_REF] Smith | Correspondence of the brain's functional architecture during activation and rest[END_REF][START_REF] Harrison | Consistency and functional specialization in the default mode brain network[END_REF]; they are called Resting State Networks (RSNs). The RSNs showed high reproducibility [START_REF] Damoiseaux | Consistent resting-state networks across healthy subjects[END_REF] and from moderate to high reliability [START_REF] Franco | Interrater and intermethod reliability of default mode network selection[END_REF][START_REF] Shehzad | The resting brain: unconstrained yet reliable[END_REF][START_REF] Zuo | The oscillating brain: complex and reliable[END_REF]. This suggests that rs-fMRI allows an exhaustive study of the intrinsic functional architecture of the human brain connectome [START_REF] Biswal | Toward discovery science of human brain function[END_REF].

Data-driven approaches include blind source separation approaches such as the independent component analysis (ICA) [START_REF] Beckmann | Probabilistic independent component analysis for functional magnetic resonance imaging[END_REF][START_REF] Mckeown | Analysis of fmri data by blind separation into independent spatial components[END_REF][START_REF] Calhoun | Unmixing fmri with independent component analysis[END_REF], the principal component analysis (PCA) [START_REF] Andersen | Principal component analysis of the dynamic response measured by fmri: a generalized linear systems framework[END_REF][START_REF] Baumgartner | Comparison of two exploratory data analysis methods for fmri: fuzzy clustering vs. principal component analysis[END_REF], the Temporal Clustering Analysis [START_REF] Liu | The temporal response of the brain after eating revealed by functional mri[END_REF][START_REF] Morgan | Development of 2dtca for the detection of irregular, transient bold activity[END_REF], and clustering methods [START_REF] Salvador | Neurophysiological architecture of functional magnetic resonance images of human brain[END_REF][START_REF] Cordes | Hierarchical clustering to measure connectivity in fmri resting-state data[END_REF][START_REF] Golland | Data-driven clustering reveals a fundamental subdivision of the human cortex into two global systems[END_REF][START_REF] Lee | Clustering of resting state networks[END_REF].

These methods are useful if the aim is grouping together voxels showing the same spatial or temporal features; in particular the spatial ICA [START_REF] Calhoun | Spatial and temporal independent component analysis of functional mri data containing a pair of task-related waveforms[END_REF] is the approach the most commonly used in rs-fMRI data analysis.

The Principal Component Analysis

The PCA, also known as Karhunen-Loeve Transform [START_REF] Jolliffe | Principal component analysis: Wiley online library[END_REF], is a statistical technique that linearly transforms a dataset into a smaller and uncorrelated set of variables that captures most of the variance in the original data. Usually, only few directions are considered to be of interest and are taken as a new coordinate system to reveal the underlying structure of the data. The resulting set of uncorrelated variables is more informative than the larger set of correlated variables. The PCA is often used as a pre-processing step to prepare data for subsequent analyses, such as clustering analysis and the ICA [START_REF] Kim | Comparison of data-driven analysis methods for identification of functional connectivity in fMRI[END_REF].

The Independent Component Analysis

The ICA was developed by Bell and Sejnowski in 1995 [START_REF] Bell | An information-maximization approach to blind separation and blind deconvolution[END_REF]. It is a powerful explorative analysis technique that allows to express a set of random variables as linear combination of statistically independent latent components. In the context of blind source separation problems [START_REF] Hyvärinen | What is independent component analysis? independent component analysis[END_REF], the ICA tries to discover hidden and statistically independent source signals only from the measured observations [START_REF] Comon | Independent component analysis, a new concept?[END_REF]. In the context of fMRI, the ICA consists of two different approaches: the spatial-ICA and the temporal-ICA. Spatial-ICA was the first ICA approach to be applied to fMRI data [START_REF] Mckeown | Independent component analysis of fmri data: examining the assumptions[END_REF]; the temporal-ICA was used for the first time in 1999, in a study conducted by Biswal and Ulmer [START_REF] Biswal | Blind source separation of multiple signal sources of fmri data sets using independent component analysis[END_REF]. The description reported below has as main reference the one proposed by Calhoun and colleagues in 2001 [START_REF] Calhoun | Spatial and temporal independent component analysis of functional mri data containing a pair of task-related waveforms[END_REF]. where Ŵs is the N × N estimated linear mixing matrix found using ICA, and C s is an N × M matrix that contains the N independent components. We can write

Spatial

X = Ŵ-1 s C s (2.10)
where Ŵs the spatially independent components (spatial maps) are in the rows of C s and the associated spatially independent time series are located in the columns of Ŵ-1 s (Figure 2.5). The sources C s are estimated by iteratively optimizing the unmixing matrix Ŵ-1 s , so that C s contains mutually independent rows, using the information-maximization (Infomax) algorithm [START_REF] Bell | An information-maximization approach to blind separation and blind deconvolution[END_REF]. Infomax algorithm [START_REF] Bell | An information-maximization approach to blind separation and blind deconvolution[END_REF] is one of the most popular neural network-based approaches. It uses a global measure of minimizing the mutual information of the components C si by maximizing the entropy of the output of a network [START_REF] Correa | Performance of blind source separation algorithms for fmri analysis using a group ica method[END_REF].

Temporal-ICA -X is a M × N matrix, corresponding to the transpose of the X matrix used for spatial-ICA. The temporal-ICA decomposition is

C t = Ŵt X (2.11)
where Ŵt is the M × M estimated linear mixing matrix found using ICA, and C t is an M × N matrix with the M independent components. We can then write

X = Ŵ-1 t C t (2.12)
where the temporally independent time courses are found in the rows of C t and the associated temporally independent maps (images) are found in the columns of Ŵ-1 t .

Temporal-ICA is typically much more computationally demanding than spatial-ICA for fMRI applications because of a higher spatial than temporal dimension and can grow quickly beyond practical feasibility. Recent studies on task-related events showed that spatial-ICA and temporal-ICA provide similar results [START_REF] Petersen | On the independent components of functional neuroimages[END_REF]. However, spatial-ICA has dominated the functional imaging literature so far, as mentioned above.

Probabilistic-ICA -An extension of spatial-ICA is the probabilistic-ICA model, which assumes that the p-dimensional vectors of observations (time series) are generated from a set of q < p (i.e., there are fewer source processes than observations in time) statistically independent non-Gaussian sources (spatial maps) via a linear mixing process corrupted by additive Gaussian noise, η(t)

x i = As i + µ + η i (2.13)
where x i expresses the individual measurements at voxel location i, s i denotes the non-Gaussian source signal contained in the data and η i denotes the Gaussian noise. The covariance of the noise is allowed to be voxel dependent to comply with the various noise covariances observed in different tissue types [START_REF] Woolrich | Mixture models with adaptive spatial regularization for segmentation with an application to fmri data[END_REF]. The vector µ defines the mean of the observations x i and the matrix A, with dimensions p × q, is assumed to be non-degenerate, i.e. of rank q. Solving the blind separation problem requires finding a linear transformation matrix W such that

ŝ = Wx (2.14)
provides a good approximation of the true source signal s [START_REF] Beckmann | Probabilistic independent component analysis for functional magnetic resonance imaging[END_REF][START_REF] Beckmann | Modelling with independent components[END_REF]. A problem of this approach is the correct choice of the number of components to extract. Dimensional underestimates may discard valuable information and result in sub-optimal signal extraction. Overestimates, however, results in a large number of spurious components and will over fit the data, harming later inference and increasing computational costs.

Nevertheless, a study by Abou-Elseoud et al. ( 2010) [START_REF] Abou-Elseoud | The effect of model order selection in group pica[END_REF] demonstrated that an optimal detection of components needs a model order within the range between 20 and 80. The final outcome is a set of independent components, some of which are clearly related to the signals of interest for rs-fMRI study, the so-called RSNs (Figure 2.6), while some are related to other physiological processes (e.g. respiration, etc) or to imaging artifacts (e.g. motion, ghosting, slice dropout, noise, etc).

Clustering approaches

Clustering methods attempt to group voxels that show a high level of similarity in their BOLD time courses into sub-groups. For example, hierarchical clustering builds a treelike structure of all the data-points [START_REF] Cordes | Hierarchical clustering to measure connectivity in fmri resting-state data[END_REF][START_REF] Salvador | Neurophysiological architecture of functional magnetic resonance images of human brain[END_REF], while in k-means [START_REF] Golland | Data-driven clustering reveals a fundamental subdivision of the human cortex into two global systems[END_REF] and c-means [START_REF] Lee | Clustering of resting state networks[END_REF] clustering algorithms, all voxels are assigned to one of several clusters on the basis of their distances from the cluster centers, which, in turn, are calculated from an average of their members.

The iterations to update memberships and cluster centers continue until convergence is achieved.

Medial visual cortical areas Lateral visual cortical areas

Auditory system Sensory-motor system Visuo-spatial system Executive control Dorsal visual stream Dorsal visual stream

Conclusion

In this chapter we proposed a state-of-the-art review of the main approaches used to study the fMRI data. The confirmatory approaches, i.e. the GLM, need a priori assumptions on the activations they are estimating. The exploratory approaches (i.e. data-driven methods) aims at blindly extract some features hidden in the data or, as in the case of deconvolution approaches, at recover the underlying neural activation without knowledge on the experimental paradigm.

Limitations of the Current State of the Art

The approaches described in this chapter present several limitations. The GLM is conceived for the recovery of voxel-wise time series of fMRI data when the experimental paradigm is given. It is then limited to the application to task-fMRI data, where the subject is asked to perform a motor or cognitive task at known time and for given duration.

As for the data-driven approaches, i.e. the ICA and clustering approaches, they do not require any prior knowledge on the experimental paradigm, therefore they are largely used for the analysis of rs-fMRI images. Regardless, these methods are useful if the aim is grouping together voxels showing the same spatial or temporal features those voxels are grouped together according to some sort of similarity, in the time or in the space domains.

They allow to recover networks, i.e. the RSNs or clusters, but they cannot be used if the goal is identifying activations at the voxel level. Furthermore, they are also limited by the necessity of choosing a priori the number of components or clusters and by the interpretation of the results [START_REF] Gaudes | Detection and characterization of single-trial fmri bold responses: Paradigm free mapping[END_REF] and they do not consider including any hemodynamic effect. In fact, the aforementioned methods cannot be used if the goal is identifying how and when the voxel's signals are activated or deactivated.

Deconvolution approaches are nowadays powerful methods to blindly recover voxel-wise fMRI activations, and they have been recently applied to rs-fMRI images. Despite the fact that they do not require any prior on the underlying activations, they have to handle the problem of space and time in the same optimization problem. In fact, these approaches split the optimization problem into two decoupled spatial and temporal regularization. In practice, this results in doubling the number of parameters to be set and requires the solver to alternate between the constraints. The two problems, i.e. one with a spatial constrain and another with a temporal constrain, are indeed solved separately. Meanwhile they are inserted in a general forward-backward splitting algorithms, and the results associated to the two regularizations are then iteratively summed-up according to user-defined weights.

In the next chapters we propose the two main contributions of this thesis, where we aim at the blind estimation of the neural activation without necessary assumptions on it.

Part IV

Contributions

First Contribution

The Forward Model

The observation model that we propose in this thesis, relating the parameters of the model to the observations is the following linear mapping

y = Ax + (3.1)
where y is the discrete acquired fMRI signal of length N, obtained starting from an unknown input signal x, of length N, that goes through some process A, of size N × N, also known as the observation model and makes the link between the observed data and the predicted data. is the noise that corrupts the observed data during the considered process and it is typically modeled as an additive noise.

The goal of this work is solving the inverse problem by recovering the original signal x from the noisy measurements y.

The Regularization Approach

The inverse problem we are addressing in this work is to find the original signal x from Eq.

(3.1) starting from the acquired noisy measurements y and the operator A.

A problem is considered well-posed if a solution exists, if this solution is unique and if the solution is stable, which means that small changes in the measurements cause small changes in the solution. Because of the presence of the operator A, if A has small eigenvalues, small perturbations on the data y may give big changes in the estimated solution x and the problems in Eq. (3.1) is ill-posed.

The fundamental solution of the inverse problem we are dealing with is given by minimizing the data fitting error, also known as least squares solution, given by

x = argmin x y -Ax 2 2 . (3.2)
Because the above minimization formulation has a solution that strongly depends on the operator A, additional constraints are introduced into the formulation to reduce the set of solutions and find a unique solution. The idea consists on finding two functionals: a data fidelity term, F (x), that measures the distance between y and Ax, and a regularization term, R(x), that favours appropriate minimizers penalizing potential solutions with undesired structures [START_REF] Benning | Modern regularization methods for inverse problems[END_REF]. The regularized problem is then defined as

x = argmin x F (x) + R(x). (3.3)
The problem in Eq. (3.3) can then be written as

x = argmin x y -Ax 2 2 + λR(x) (3.4)
where: y -Ax 2 2 is the fidelity term, R(x) is a penalization (or regularization) term that encodes some constraints on x, which lead to desirable properties of the solutions, and λ is the regularization parameter that controls the impact of each term into the optimization process: i.e. the higher the λ the more the solution is regularized, the lower the λ the more the solution is given only by the contribution of the fidelity term.

There exists a wide range of methods which propose different choices for the regularization term R(x) to encode different a priori hypotheses on the signal x. Many studies on signal and image processing have been focused on optimization problems based on regularization approaches to exploit some features of the desired solution and to keep the computational costs acceptable. The most common approach is based on the regularization of the Euclidean norm, also known as Tikhonov regularization or ridge regression [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF] that is formulated as follows:

x = argmin x y -Ax 2 2 + λ Fx 2 2 (3.5)
where λ is again the regularization parameter and F is a linear mapping. This problem has a unique solution that is obtained by deriving the expression in Eq. (3.5), since both the fidelity term and the regularization term are quadratic and differentiable. The minimization then is: x = (A T A + λF T F) -1 A T y. This approach aims at minimizing the energy using a l 2 -norm thus yielding to smooth solutions.

Because in this work we hypothesized that the brain activates in blocks of constant activity, we will now focus on sparsity-pursuing methods, that indeed allow to represent the signal with least amount of elements.

The highest sparsity can be achieved via a l 0 -(quasi)norm, that promotes solutions with only few large coefficients as in

x = argmin x y -Ax 2 2 + λ x 0 . (3.6)
Nonetheless, the regularization problem is not convex and hence costly in terms of computation.

The l 1 -norm in Eq. (3.7), which has been introduced as a convex approximation of the l 0 -(quasi)norm, is instead computationally feasible.

x = argmin x y -Ax 2 2 + λ x 1 (3.7)
where

x 1 = ∑ i |x i |.
This approach is also known as Least Absolute Shrinkage and Selection Operator or LASSO regression [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF].

The introduction of sparsity-pursuing regularization methods introduces non-quadratic functions, therefore differentiability is not always assured. Several methods exist to solve the problem in Eq. (3.7). The Least Angle Regression (LARS) [START_REF] Efron | Least angle regression[END_REF] is an algorithm that outputs the solution for all lambdas. The whole path of LASSO solutions is generated for "practically the cost of one least square calculation on the data" [START_REF] Efron | Least angle regression[END_REF]. Other approaches such as proximal gradient descent (Iterative Shrinkage Thresholding Algorithm (ISTA) [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] and the Fast Iterative Shrinkage Thresholding Algorithm (FISTA) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]), unlike LARS, can be used for LASSO generalizations as for example the case of TV-norm regularized problem

x = argmin x y -x 2 2 + λ ∆ D x 1 (3.8)
where ∆ D x = x(t n )x(t n-1 ) is the finite difference operator. The TV-norm regularization promotes piece-wise constant signals [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF].

In the context of fMRI, this idea has been generalized, thus leading to the generalized TV (gTV) for fMRI data in the work proposed by Karahano glu and colleagues in 2015 [START_REF] Karahanoglu | A signal processing approach to generalized 1-d total variation[END_REF]. The problem then can be reformulated as

x = argmin x y -x 2 2 + λ ∆ L x 1 (3.9)
where ∆ L x 1 = ∆ D ∆ L h accounts for a differential operator ∆ D and an additional linear operator ∆L h representing the hemodynamic inverse filter [START_REF] Karahanoglu | A signal processing approach to generalized 1-d total variation[END_REF]. The fundamental idea is to explore sparseness in a transformed domain, which in the case of fMRI is given by the presence of the BOLD signal that is accounted for via the convolution between the HRF and the signal representing the activation.

fMRI Image Structure

In the previous section we provided a state-of-the-art review of regularization approaches and optimization algorithms to have an overview of the tools we used in our work. In the next paragraphs we will describe our first main contribution. To do this, we initially give a representation of the fMRI data that we will be employing later on. fMRI scans produce a set of 3-D images recorded over time, therefore a fMRI image is represented as a 4-D data with a spatial resolution, given by the size of the voxels in the acquired brain volume and a temporal resolution, that is directly related to the TR. The TR is the time between two given successive RF pulses and it corresponds in fMRI to the time necessary to collect one entire brain volume. Thus fMRI images have two main complementary information: one given by the 3-D spatial dimension, and the other given by the 1-D temporal dimension, that is given because each fMRI volume is acquired several times to have access to the dynamic of the brain activity over time. An illustration of the 4-D fMRI data is showed in Figure 3.1.

The Hemodynamic Response Function

The BOLD model employed in this work comes from the one proposed by Friston and colleagues [START_REF] Friston | Nonlinear responses in fmri: the balloon model, volterra kernels, and other hemodynamics[END_REF]. It combines the Balloon/Windkessel model [START_REF] Buxton | A model for the coupling between cerebral blood flow and oxygen metabolism during neural stimulation[END_REF][START_REF] Buxton | Dynamics of blood flow and oxygenation changes during brain activation: the balloon model[END_REF][START_REF] Mandeville | Evidence of a cerebrovascular postarteriole windkessel with delayed compliance[END_REF] with a model of how the neural activations trigger changes in blood flow.

The model is made of a non-linear state-space system composed by 4 state variables (Eq.

3.10) and the observed BOLD signal y(t) = λ(v, q, E 0 ) (Eq. 3.11).

The system mentioned above is the following: The observed BOLD signal is

                   τ 0 v = f in -v 1/α τ 0 q = f in 1 -(1 -E 0 ) 1/ f in E 0 v 1-α α q ḟin = s ṡ = u - s τ s - f in -1 τ f (3.10)
y(t) = λ(v, q, E 0 ) = V 0 (k 1 (1 -q) + k 2 (1 - q v ) + k 3 (1 -v)). (3.11)
The constants and the variabels appearing in the hemodynamic model are summurized in The first equation in (3.10) represents the rate of the volume change; the second equation represents the change in the deoxyhemoglobin which reflects the amount of deoxyhemoglobin arriving at a venous compartement and the amount exiting from it; the third equation links the synaptic activity and the regional cerebral blood flow; the fourth one describes the signal as generated by the neural activity and other neurogenic and diffusive signal subcomponents [START_REF] Friston | Nonlinear responses in fmri: the balloon model, volterra kernels, and other hemodynamics[END_REF]. As proposed by Khalidov and colleagues [START_REF] Khalidov | Activelets: Wavelets for sparse representation of hemodynamic responses[END_REF], to be able to use the hemodynamic model for linear convolution, it is linearized and considered as time-invariant. The system in (3.10) was linearized around the resting point: {x 1 , x 2 , x 3 , x 4 = 0, 0, 0, 0} [START_REF] Khalidov | Activelets: Wavelets for sparse representation of hemodynamic responses[END_REF] by defining the following variables

Symbol

k 1 BOLD constant 1 7E 0 k 2 BOLD constant 2 2 k 3 BOLD constant 3 2E 0 -0.2
             x 1 = s x 2 = 1 -f in x 3 = 1 -v x 4 = 1 -q (3.12)
Therefore the system in (3.10) becomes

                         ẋ1 = u - 1 τ s x 1 + 1 τ f x 2 ẋ2 = -x 1 ẋ3 = 1 τ 0 x 2 - x 3 α ẋ4 = cx 2 - 1 -α ατ 0 x 3 - 1 τ 0 x 4 (3.13)
where

c = 1 + (1 -E 0 ) ln(1 -E 0 )/E 0 τ 0

This linearization led to the definition of the hemodynamic response function (HFR) whose

Fourier transform is

H(ω) = K (jω -z) (jω -p 1 )(jω -p 2 )(jω -p 3 )(jω -p 4 ) (3.14)
where z = -

(k 1 + k 2 ) 1 -α ατ 0 - c α - k 3 -k 2 τ 0 -(k 1 + k 2 )cτ 0 -k 3 + k 2 p 1 = - 1 τ 0 p 2 = - 1 ατ 0 p 3,4 = - 1 2τ s 1 ± j 4τ 2 s τ f -1 K = V 0 τ 0 -(k 1 + k 2 )cτ 0 -k 3 + k 2
From now on, when we refer to the operator H we refer to the circulant matrix for the hemodynamic response (h(t)) of size N × N where N is the length of a fMRI time series.

The Forward Model of fMRI Data

As already proposed in the state of the art [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF][START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF], in this work we modeled the fMRI time series supposing blocks of neural activations. We name the neural activation u(t) and we refer to it as activity-inducing signal, out of previous works [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF][START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF]. The activity-inducing signal u(t) is modeled as a piece-wise constant function. The activity-related signal x(t) is the convolution between the HRF and the activity-inducing signal

x(t) = u(t) * h(t) (3.15)
where * is the operator of the convolution, h(t) is the HRF, and u(t) is the activity-inducing signal. The acquired fMRI time course is given by the activity-related signal corrupted with noise due to artifactual sources such as cardiac rate and respiratory fluctuations, movements in the scanner, and thermal noise [START_REF] Karahanoglu | Sparsity-promoting spatiotemporal regularization for data mining in functional magnetic resonance imaging[END_REF][START_REF] Lund | Non-white noise in fmri: does modelling have an impact?[END_REF]. The noisy fMRI time courses were therefore modelled as

y(t) = x(t) + a = u(t) * h(t) + a (3.16)
where a is the additive random Gaussian noise with zero mean and standard deviation σ a .

L1-norm Temporal Regularized Deconvolution of the fMRI BOLD Signal

In this section we will introduce the first contribution of this thesis, where we develop a temporal-regularized deconvolution of fMRI time series.

In this contribution, the goal we are addressing is the paradigm-free recovery of the fMRI brain activation supposing the brain activates in constant blocks. To do this we propose a deconvolution technique that has the advantage of not requiring necessary and a priori information on the timing and duration of the underling activation, in contrast with the confirmatory approach discussed in the background.

Definition of LASSO Optimization Problem

In our first contribution, we introduce a regularized formulation in a minimization of a cost function; the problem is therefore composed by a data fitting term and a regularization term. We solved the inverse problem thus reconstructing the activity-inducing signal u(t), that refers to the neural activation, using a deconvolution technique in absence of external references represented by the stimulation paradigm.

To do this, we started from the idea proposed by Karahano glu and colleagues (2013) [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF] and lately revisited by Farouj et al. (2017) [START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF] who recover the brain activation exploiting the fact that the derivative of a piece-wise constant signal, i.e. the brain activation, is a sparse signal (see Section 2.2 for more details). This signal is named innovation signal s(t).

We implemented the deconvolution exploiting the sparsity of the innovation signal s(t) [START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF],

and we built the following LASSO model with a l 1 -norm regularization on the temporal dimension of the fMRI data

ŝ = argmin s 1 2N ||y -As|| 2 2 + λ||s|| 1 (3.17)
where

• y is the acquired and noisy fMRI signal;

• ŝ is the recovered innovation signal of length N, that is expected to be sparse because it is the derivative of a piece-wise constant activity;

• N is the length of the acquired and noisy fMRI signal;

• s is the innovation signal of size N;

• A in As = HI α s is a linear operator where H is the N × N circulant matrix for the hemodynamic response and I α is the N × N circulant matrix for the α-filter that is explained in the next paragraph. A is the resulting operator of size N × N;

• λ is the regularization parameter;

• || • || 2 indicates a l 2 -norm; • || • || 1 indicates a l 1 -norm.
All the signals and operators we are naming in this section are illustrated in Figure 3.2.

Starting from the estimated optimal solution ŝ(t), the activity-inducing signal, û(t), in other words the recovered neural activation, was hence obtained by the convolution of ŝ(t) with the Heaviside function H(t)

û(t) = H(t) * ŝ(t) (3.18)
where û(t) and ŝ(t) represent the recovered activity-inducing signals and innovation signals respectively.

The α-Filter Design

To exploit the sparsity of the innovation signal s(t), and to obtain the activity-inducing signal u(t) starting from s(t), we implemented the operator I α that refers to an exponential accumulation function, normalized by the factor S, whose z-transform is The α-filter design is extensively represented in Appendix A. The interesting point of this filter is that the shape of I α (Figure 3.3) depends on the parameter α. Hence, for large α, typically greater than 1.5, according to the level of noise, the operator works as a finite difference ∆ t .

I α (z) = S e -α z -1 (1 -e -α z -1 ) 2 1 (1 -z -1 ) 2 - e -α z (1 -e -α z) 2 1 (1 -z) 2 (3.19)
We choose to design the α-filter, and not to use an operator implementing a finite difference, ∆ t , as in [START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF][START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF] because we expect it to further smooth the noise still present on the data, for example, in the presence of blocks of spurious and false activation, given for example by artifacts that has not been removed during the preprocessing [START_REF] Power | Spurious but systematic correlations in functional connectivity mri networks arise from subject motion[END_REF]. It must be stressed that the activity-inducing signal u(t), showed in Figure 3.2, is represented as a piece-wise constant function as it was previously proposed in literature by Karahanoglu et al. ( 2013) [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF]. The effect of I α could smooth the activity-inducing signal u(t), but this happens only for small values of the parameter α (typically less than 0.75). Whereas, for big values of α (typically greater than 1.5), the filter described by the function I α (t) corresponds to a sharp accumulation function, therefore u(t) corresponds exactly to a piece-wise constant signal.

For great values of α, the function I α (t) is a sharp accumulation function, that, if reverted, corresponds to a finite difference operator. Therefore, the innovation signal s(t) results to be a sparse signal. The parameter α plays a role such that the appropriate filter smooths the noise while it keeps the solution as close as possible to a piece-wise constant signal. (Please pay attention to Appendix A in the thesis, where we clearly show that the derivative filter operator converges thanks to its normalization factor to a simple difference operator for α

2 ), paramtrized by the regularization parameter λ. Each black dot in the curve is related to each λ outputted by the Least Angle Regression (LARS) algorithm [START_REF] Efron | Least angle regression[END_REF]. The red dot represents the optimal solution as the nearest to (0, 0) ∈ R 2 .

large (typically greater than 1.5), i.e [ f (n + 1)f (n -1)]/2).

Solution of the LASSO problem using the LARS Algorithm and L-Curve

To solve the optimization objective introduced in Eq. (3.17), we minimized the objective function using the LARS algorithm [START_REF] Efron | Least angle regression[END_REF], that outputs at once all lambdas of interest and their associated solutions. Practically all regularization methods for computing stable solutions to inverse problems involve a trade-off between the how much the solution is regularized and the how much it fits the original data. The various regularization approaches aims at finding a optimal measure of these quantities, and the optimal trade-off between them. If too much weight is given to the regularization, the solution will not fit the data properly and the residuals ||y -As|| 2 2 (3.17) will be too large. On the contrary, if the regularization is insufficient then the fit will be good but the solution will be dominated by data errors and ||s|| 1 (3.17) will be too big [START_REF] Hansen | The l-curve and its use in the numerical treatment of inverse problems[END_REF]. Therefore, given the important roles of the weights assigned to the two norms presented in Eq. (3.17), one related to the solution and one to the residuals, and given that the LARS algorithm outputted all the solutions associated to all the important lambdas, we plot for each lambda the two quantities (||y -Aŝ|| 2 2 , ||ŝ|| 1 ) in the so-called L-curve [START_REF] Hansen | The l-curve and its use in the numerical treatment of inverse problems[END_REF]. We choose as optimal λ the one associated with the closest point to the origin of the Cartesian axes in the L-curve. Finally, the solution ŝ corresponds the one associated to the optimal λ. An exemplification of L-curve for one voxel is shown in Figure

Leave-one-out Cross-Validation

We compared results using the L-curve with the ones obtained using a cross validation approach. Given a time series of N time points, we implemented a Leave-one-out cross validation (LOO-CV) [153] as shown in Algorithm 2. The LOO-CV was implemented such that for a time series y(t), composed by N time points, the i-th sample y[i] was removed, thus resulting in a new time course y (t) composed by N -1 time points. The LARS algorithm was applied to solve the new objective function

ŝ (t) = argmin s 1 2(N -1) ||y (t) -A s (t)|| 2 2 + λ ||s (t)|| 1 (3.21)
where A represents the operator A, after removing the i-th row from it. Because LARS gives a set of solutions associated to a set of lambdas, for each of these solutions, the removed sample of the original time series y(t) was estimated starting from the new solution ŝ (t)

y r [i] = A[i, :] • ŝ . (3.22)
After that for each lambda, meaning for each solution, the squared errors between the original time sample y[i] and the estimated one y r [i] were computed. This process was repeated N times, ones for each time sample removed. Because LARS gives a new set of lambdas (λ) for each i-th removed sample, the errors were then interpolated to have error values associated to an array of fixed lambdas in the range [0, 0.2]. All the procedure indicated above was repeated V times, for each simulated fMRI time series, and the errors related to each fixed lambdas were averaged across the V repetitions of the simulated corrupted fMRI time courses. The optimal λ was then chosen as the one corresponding to the least mean squared error. Finally, the LARS algorithm was run on the original optimization problem and the solution

ŝ = argmin s 1 2N ||y -As|| 2 2 + λ OPT ||s|| 1 (3.23)
was the one corresponding to the lambda the closest to the optimal one previously computed.

To test the difference between the lambda selection methods that we considered, i.e. the Lcurve and the LOO-CV, we simulated fMRI acquired time series corrupted with different amount of noise and we tested different alphas for the α-filter. In the next section we will describe how we simulated the fMRI data to test and validate the above-mentioned approaches.

As for the accumulation function I α (t) (see Section 3.4.2), as it was previously explained, it has a shape and therefore a performance that can be adjusted by tuning the parameter α. This parameter applies such that: the lower the α the smoother the function I α (t); the higher the α the sharper the function I α (t). The choice of the parameter α would play a role on the signal sparsity. The accumulation function I α (t) has been proposed here because it allows smoothing the noise, nevertheless a too small α would impact on the sparsity of the resulting activation. For this reason, a trade-off between the value of α and the sparsity of the resulting signal is necessary. The parameter α was indeed chosen such that it does not impact considerably on the sparsity of the resulting innovation signal s(t).

Algorithm 2 Leave-one-out cross-validation (LOO-CV) Input: V fMRI corrupted time series y of length N.

Output: Regularization parameter λ OPT that provides the optimal solution ŝ for the problem in Eq. (3.17).

1: for v = 1 : V do 2:

for i = 1 : N do 3:

Computation of the vector y (i) by removing the i-th element (y[i]) from y.

4:

Computation of the operator A(i) by removing the i-th row from A.

5:

Application of LARS to solve the new objective function

ŝ λ (i) = argmin s 1 2(N -1) ||y (i) -A (i)s (i)|| 2 2 + λ (i)||s (i)|| 1 gives L solutions ŝ λ (i) for each λ ∈ Λ L (i))
where Λ L (i) is the set of L λ .

6:

for l = 1 : L do 7:
Computation of the removed sample y r (i, l)

y r (i, l) λ = A[i, :] ŝ λ (i, l) 8:
Computation of the squared errors for each

λ ∈ Λ L λ (i, l) = (y r (i, l) λ -y[i]) 2 9:
end for 10:

Given the set of λ (i, l) obtained for for each λ ∈ Λ L (i), interpolation of λ (i, l) and sampling of the obtained interpolant at the values contained in Λ F v

I (i) = interp( (i, l), λ L (i, l))(Λ F )
where Λ F is an array of fixed lambdas in the range [0, 0.2]. Given the fact that LARS outputs different lambdas for each i, the interpolation allows to have the errors as a function of Λ F .

µ v ( v I (i)) : Λ F ∈ R +
13: end for 14:

λ * = argmin λ∈Λ F µ v ( v I (i))(λ) 15: Given Λ N from LARS applied to solve ŝ = argmin s 1 2N ||y -As|| 2 2 + λ OPT ||s|| 1 16:
The optimal solution ŝ is the one associated with

λ OPT = argmin λ∈Λ N (|λ -λ * |)

Simulation of fMRI Time-Courses

To reproduce the acquired fMRI signals we simulated the activity-inducing signal as a boxcar function

u(t) = H(t -a) -H(t -b) (3.24)
where H(t) is the Heaviside step function.

We corrupted them with block-type noise, b , to simulate false activations related to motion artifacts [START_REF] Power | Spurious but systematic correlations in functional connectivity mri networks arise from subject motion[END_REF] and added noise to u(t) representing the random intrinsic electrical fluctuations within neuronal networks. These fluctuations are not associated with a response to internal or external stimuli. This noise aims at modelling the synaptic noise that refers to the constant bombardment of synaptic activity in neurons. To do this, we corrupted the activity-inducing signal u(t) with an additive random Gaussian noise with zero mean and standard deviation σ m that we called "model noise" m , yielding

u n (t) = x(t) + b + m . (3.25)
We simulated the activity-related signal x(t), consequent to the neural activation as the convolution of the activity-inducing signal with the HRF h(t) [START_REF] Khalidov | Activelets: Wavelets for sparse representation of hemodynamic responses[END_REF]:

x(t) = u n (t) * h(t). (3.26)
Real time series acquired using the fMRI technique are corrupted by different kinds of noise and artifacts given by mechanisms which are not reflecting any neurophisiological function, such as the heart rate and the respiratory fluctuations, motion artifacts, thermal noise and scanner drifts [START_REF] Lund | Non-white noise in fmri: does modelling have an impact?[END_REF]. For this reason we added noise to x(t) thus simulating the acquired fMRI signals

y(t) = x(t) + a = u n (t) * h(t) + a (3.27)
where a is the additive random Gaussian noise with zero mean and standard deviation σ a .

A scheme representing the generation of the phantom fMRI time series is shown in Figure b is the block type noise and m the model noise which were consecutively added to the activity-inducing signal u(t) thus leading to u n (t) before the convolution with the hemodynamic response function H. x(t) is the activityrelated signal, a is the additive gaussian noise and y(t) is the simulated acquired fMRI time course.

L-curve vs Cross-Validation

In this section we will test our approach in order to assess the performance of L-curve and compare the results with the ones obtained using the LOO-CV approach, both adopted for the selection of the most appropriate regularization parameter λ for solving the problem in (3.17).

Methods

We tested the difference on the solution to the optimization problem in (3.17) with the two different approaches for the selection of the regularization parameter λ: L-curve and LOO-CV.

To do this, as explained in section 3.5, we simulated fMRI time series of 100 s with a TR of 1 s. We set the activation with an onset at 20 s and the offset at 40 s, for a duration of 20 s. The amplitude of the activation was set to 1. We randomly added up to 10 blocks of false activations ( b ), with a maximum amplitude of 0.7, a minimum duration of 3 s and a maximum of 7 s. We set 10 experiments in which the data were corrupted by different amounts of noise ( m and a ) as shown in Table 3.3. For each experiment, the time courses corrupted with random noise were simulated 100 times to have a sample for statistical analysis. In Table 3. where u, u n , y and x are the ground truth and the noisy activity inducing signal, the synthetic acquired fMRI time course and the activity-related signal, respectively. As shown in Table 3.3, we tested the data either with and without adding model noise.

In order to assess the effect of the α-filter, we tested different values for α = (0.3, 0.4, 0.5, 0.75, 1, 2, 3).

To compare the results, for each experiment we computed the roots of the mean squared errors (RMSE) and its related standard deviations (STD) between the ground truth, simulated as a piece-wise constant activation u(t), and the recovered activity-inducing signal û(t) averaged across the 100 repetitions.

Results

Results are shown in Figure 3.6 and 3.7. They show that the RMSEs obtained using the Lcurve (in red) are lower than the ones obtained using the LOO-CV-based approach. This highlights the fact that the solution to (3.17) using the L-curve-based approach is closer to the ground truth if compared with the ones obtained using the LOO-CV-based approach. Results also highlight that the solution obtained using the LOO-CV (in light blue) is much more sensitive to model noise, i.e. m from experiment 1 to 5 (Figure 3.6), with respect the solutions obtained with the L-curve. However, Figure 3.7 shows that the LOO-CV-based approach gives lower errors if only a is added to the data, compared with the cases where m and a (Figure 3.6) were both added to the simulated data. These trends are also shown in the exemplifying time series plotted in Figure 3.8 and 3.9. In particular, the plots in Figure 3.8,

6 0 → ∞ 0.1 28
7 0 → ∞ 0.3 18
8 0 → ∞ 0.5 14
9 0 → ∞ 0.7 11
10 0 → ∞ 0.9 8
where both m and a were added to simulate the acquired time series (in green), show that the recovered activations using L-curve (in dotted red) are less sensitive to noise, compared to the ones recovered using LOO-CV (in light blue). The plots in Figure 3.9, show instead good results for the LOO-CV-based approach, nevertheless the L-curve assures a solution with a amplitude closer to the ground truth than those obtained with the LOO-CV. As for the effects of the α-parameter, globally it emerges, from Figure 3.6 and 3.7, and Figure 3.8 and 3.9, that the error is higher and more disperse for little α and it decreases while increasing α for both approaches. It is also shown that globally at a α of 0.75 the error stabilizes to a certain value, meaning that a sharper α-filter would not significantly change the results.

Both Figure 3.6 and 3.7 refer to synthetic data. For the approach involving the combination of the LARS algorithm and L-curve (red curves) and the one involving LARS and the LOO-CV (blue curves), we observe as expected that the errors and its standard deviation decreases while decreasing the noise. We also observe that if we increase the value of α involved in the α-filter the errors decrease. This happens because, for very small α the α-filter gives a smooth effect that is then encountered in the error, therefore the errors for low α are greater with respect to smaller α, that instead favors sharper solutions. Given that the ground truth is supposed to be a piece-wise constant function, the α-filter with higher values for α gives a solution that is closer to the ground truth. Nevertheless, as mentioned before, at a value of α = 0.75 the error stabilizes and the compromise between the noise reduction and the sparsity of the solution is found. If we observe Figure 3.8 and 3.9, the resulting û(t) obtained using the LARS algorithm combined with the L-curve (in red) presents very sharp transitions and follows the ground truth with a very small error, this suggests that the recovered innovation signal ŝ(t) is sparse and corresponds to the onset and offset of the underlying simulated activation. 3.3. For each experiment, the RMSE were averaged across 100 repetitions of the simulated time series. The errors (in y-axes) are plotted with respect to the α value (in x-axes) used for the α-filter. In the 5 experiments, we added block-type noise ( b ) for false activations, and both m and a to the data, increasing it from experiment 1 to 5. The plots show that the LOO-CVbases approach is more sensitive to noise while increasing it and results are more disperse if compared with the L-curve-based ones. In all experiments the results obtained with L-curve show errors that are significantly lower than those given by the LOO-CV. With respect to the parameter α, results obtained via the LOO-CV are improved by increasing it, instead the ones obtained with L-curve are more stable with respect to it. and the roots of standard deviations computed between the ground truth activation and the estimated one using the Least Angle Regression combined with the L-curve (in red) and with the leave-one-out cross-validation (LOO-CV, in light blue). Each plot refer to a different experiment, as reported in Table 3.3. For each experiment, the RMSE were averaged across 100 repetitions of the simulated time series. The RMSE (in y-axes) are plotted with respect to the α value (in x-axes) used for the α-filter. In experiments from 6 to 10, we increasingly added a to the synthetic times series. The plots show that the LOO-CV-based approach is less sensitive to noise if compared to experiments 1 to 5 in Figure 3.6, but still more sensitive to noise while if compared with the L-curve-based ones. In all experiments the results obtained with L-curve show errors that are significantly lower than the ones obtained with the LOO-CV for α values lower than 0.5. For α > 0.75 results are similar. The standard deviations of m and a increases from experiment 1 to 5. The plots show that the recovered activations using the L-curve-based approach are much closer to the ground truth compared with LOO-CV. The L-curvebased approach shows results that are closer to the ground truth in terms of amplitude for bigger alphas. For α = 0.75, 3, results are similar, while in contrast with those given by α = 0.3. Experiment 10

u(t) y(t) u*_CV(t) u*_L-curve(t) _LOO-CV (t) _L-curve (t)
FIGURE 3.9: The graphs in the plots report examples of the recovered activation using the mixed Least Angle Regression (LARS)-L-curve approach (in red) and the mixed LARS and leave-one-out cross-validation (LOO-CV) (in light blue) superimposed on the ground truth activation in black and the noisy synthetic fMRI time series (in green). Each row shows a different experiment (see Table 3.3), meaning that different amount of noises were applied to the ground truth activation Each columns refer to a specific α used in the α-filter. In experiments 6 and 10 time series are corrupted with a with increasing standard deviation. The different amount of noise are described in Table 3.3. The recovered activations using the L-curve-based approach are closer to the ground truth compared with LOO-CV. The Lcurve-based approach shows results that are closer to the ground truth in terms of amplitude for greater alphas α = 0.75, 3, in contrast with α = 0.3.

Validation on Phantom fMRI Images

In this paragrah we validate our approach on the phantom fMRI data and we compare it with the state-of-the-art approach proposed by Karahanoglu and colleagues [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF], namely the TA. Given the results illustrated in the previous section, from now on when we refer to our approach we intend the one given by the combination of the LARS algorithm and the L-curve.

Methods

To simulate a 4-D fMRI image, similarly to the work proposed by Farouj and colleagues [START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF],

we scaled a 3-D activation map computed with the Functional Magnetic Resonance Imaging of the Brain (FMRIB) Software Library (FSL) simulation tool Possum (https://fsl.fmrib.

ox.ac.uk/fsl/fslwiki/POSSUM/UserGuide) in the range [0, 3], with a resolution of 2 × 2 × 2 mm 3 (Figure 3.10).

We multiplied this activation map by two different block-type signals of 200 s, u(t): the first one, A, with 4 onsets; the second one, B, with one long onset (Figure 3.11). The resulting phantom image is then made of "activated" voxels and background voxels where the activation was set to 0.

Starting from this two activation paradigms, we simulated the time series for each brain voxel as described already in Section 3. Experiments set for testing our temporal regularized deconvolution on phantom fMRI image. The first column indicates the experiment number. The second and third columns report the standard deviation of the model noise (σ m ) and the additive gaussian noise (σ a ), respectively. The fourth and fifth columns show the pSNR computed as in (3.30), expressed in dB, for activations A and B.

background ones and the "active" ones. The goal is to see if the approach that we propose is able to distinguish between an activated and not activated voxel and to assess how much the recovered activations are close to the ground truth.

We set up five experiments with different standard deviations for the added noises m and a as shown in Table 3.4. We compute the pSNR as the maximum pSNR computed across the voxels belonging to the activation map showed in Figure 3.10, that is

pSNR = max 10 log 10 max(u) 2 Var[y -u] . (3.30)
To evaluate the results obtained with the procedure described above, we computed the roots of the mean square errors (MSEs) and standard deviation (STDs) between u(t) and û(t) averaged among the voxels belonging to the GM masked activation. On the light of the results discussed in the previous chapter, the results we show here are related to α = 0.75 employed in the α-filter. We compared our results with those obtained with the temporal regularization implemented in the Total Activation (TA) toolbox [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF]. For the moment, because we proposed a temporal regularized deconvolution approach, we disregarded the spatial regularization and considered only the temporal regularization implemented in the TA approach in order to allow a comparison between the two approaches.

Results

In Figure 3.12 we show that using our approach we are able to clearly discern between a voxel that is activated (light blue box) and a voxel that instead is not (yellow box).

Table 3.5 shows that the roots of MSEs ± STDs change for different pSNRs and that they are significantly lower than the ones obtained using TA. Figure 3. [START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF] shows examples of the reconstructed activity inducing signal using our approach ( ûL-curve (t)) and the TA approach( ûTA (t)).

The value of lambda used in the TA for the temporal regularization, tuned for each voxel as the median absolute deviation of fine-scale 3rd order Daubechies wavelet coefficients, has been found to be greater than the one that we propose here using our approach, that combines the LARS algorithm with the L-curve. If too much weight is given to the

A B time [s] time [s] time [s] time [s] u(t) y(t)
u*_L-curve(t) _L-curve (t) FIGURE 3.12: Reconstructed activity-inducing signal û(t) obtained with our approach (red) superimposed on the ground truth (u(t), black) and simulated fMRI signal (y(t), green). The plots in the blue square are taken from exemplificative voxels and are related to activation A with pSNR = 5.17 dB; the ones in the yellow square are related to activation B with SNR = 5.12 dB. Vertically, for both A and B, the plots show that our approach is able to discern between an activation and a non-activation.

.13: Reconstructed activity-inducing signal û(t) obtained with our approach (red) and the total activation (TA, blue) superimposed on the ground truth activation (black) and simulated fMRI signal (green). The plot on the top is related to activation A with pSNR = 5.17 dB. The plot on the bottom is related to activation B with pSNR = 5.12 dB.

regularization, the solution does not fit the data properly, therefore the residuals will be too large. The solution will be indeed far from the original data, and as in this case the amplitude of the recovered activity-inducing signal obtained using the TA is lower than the one of the original data. It has been shown in literature that the total variation denoising, formulated as l1-norm penalty term, tends to underestimate the amplitude of the signal. This is indeed a limitation of the l1-norm regularization. A large value of the regularization parameter lambda, as it is the case for the TA, makes the l1-norm solution sparser, but it reduces the coefficient amplitudes [START_REF] Selesnick | Sparse regularization via convex analysis[END_REF]. As for the TA approach, that in this thesis was used as provided in the TA toolbox, it implements a temporal regularization that aims at minimizing the total variation of the signal. It must be stressed that all the parameters were not set by the user, but for the thesis we run the TA approach implemented in the TA toolbox and all the parameters were set by the tool itself and not by the user.

Regarding Table 3.5 and Figure 3.13, if we would compute the Pearson correlation coefficients between the ground truth u(t) and the recovered activation û(t) using the TA and using our approach results would be similar. Nevertheless, if we are looking for a voxel that is active/not-active at a certain time and for example we would apply a threshold to define when a voxel is activated/deactivated, using our approach would allow to use a higher threshold and this would allow to not to encounter spurious and false activations with lower amplitudes.

Application on Real task-fMRI Data

In this section we validate our approach on real task-fMRI data and we compare our results with the TA. The reason why we choose to use task-fMRI data, even if the final goal is the application to rs-fMRI data is because in this way we have a ground truth to which to refer to in order to assess the performance of our approach.

Methods

We used the preprocessed task-fMRI image of subject 100307 taken from the Human Connectome Project (HCP) database [START_REF] Van Essen | The wu-minn human connectome project: an overview[END_REF]. The data were acquired using a 3T SIEMENS MAGNETOM Connectome Syngo MR D11 scanner, using a multi-band accelerated EPI sequence ( Acquisition Time = 204.48 s; TR = 720 ms; Echo Time (TE) = 33,10 ms; flip angle (FA) = 52°; 72 slices; Field of View (FOV) = 208 mm × 180 mm; 2 mm isotropic resolution; 284 volumes)1 . The data underwent a minimal pre-processing pipeline [START_REF] Glasser | The minimal preprocessing pipelines for the human connectome project[END_REF] which includes: correction of gradient-non linearity-induced distortions; registration of each image frame to the signal-band reference image to achieve motion correction; phaseencoding distortion correction; EPI image distortion correction; registration of the fMRI volumes to the structural data; coregistration of the fMRI data to the MNI space; masking and fMRI image intensity normalization to the 4-D whole global mean of 10 3 . Furthermore, as additional pre-processing steps, each voxels' time series were then detrended to remove linear trends and finally normalized to 0 mean and unit standard deviation.

The task involves, for known timing and duration, the movement of the tongue, as shown in Table 3.6.

Task-fMRI data provided us with a ground truth that allowed us to assess the quality of results. The ground truth is given by the onset, offset and duration of the task. The reconstructed û(t) were averaged in a ROI of 6 × 6 × 6 mm 3 centered in the Brodmann Area 4p (MNI coordinates: 62, -14, 30) [START_REF] Kiviniemi | Independent component analysis of nondeterministic fmri signal sources[END_REF].

Results

Results on real task-fMRI data are showed in Figure 3.14. They show that, given the real task-fMRI data (in green), with our approach (in red) we were able to clearly recover the brain activations, supposingly happening during the gray areas, compared with the ones recovered using TA (in blue). As we showed for the simulated data, in Figure 3.8, 3.9, 3.12 and 3.13 using our approach allows a deconvolution of the BOLD signal from fMRI data, thus our approach allows to recover the simulated block-type activation. As for real data, in Figure 3.14, given that the task performance is subject-dependent, this could explain why the results are subject to a certain variability, meaning that the recovered activation û(t) is longer than the task duration. Moreover, we can observe in the figure that the recovered solution using our approach (in red) follows the variance of the real task-fMRI time course in green and retrace an activation that is underlined below it.

Discussion and Contributions of this Chapter

In this chapter we introduced our first main contribution. Initially inspired by the work proposed by Karahano glu et al. in 2013 [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF], we proposed an approach to solve an optimization problem with the purpose of recovering the brain functional activation in time, supposing the brain activates for a period of time and finally deactivates.

To solve the inverse problem, we proposed a combined approach that involves the use of the LARS algorithm and the L-curve. As a first point, the LARS algorithm outputs at once all important lambdas and their associated solution. Given the set of lambdas, we proposed two approaches for choosing the optimal one: the L-curve and the LOO-CV. Results showed that the L-curve allows a faster and better estimation of the appropriate regularization parameter lambda for the given data, compared with the LOO-CV. In fact the LOO-CV gives still a piece-wise constant solution but it is more sensitive to noise thus leading to the recovery of false activations. The L-curve instead, in almost all the test we performed in this thesis, is much less sensitive to noise and is able to select a solution that is very close to the ground truth, both in terms of amplitude and activation's onset and offset. From the computational side, using the LOO-CV requires to run several times the LARS algorithm to solve the optimization problem and to select the optimal lambda, whereas the combination between the LARS algorithm and the L-curve has to be run just once, therefore it is less costly.

In the experiments we showed in this chapter, the use of different α parameters affected mostly the LOO-CV-based approach than the L-curve-based one.

Our approach applies both on synthetic data, giving as output a curve that is very close

y(t) u*_L-curve(t) _L-curve (t) u*_L-curve(t) _TA (t)
[TR] FIGURE 3.14: The plot shows the û(t) obtained with the Human Connectome Project (HCP) data. Results obtained using our approach (red) and the TA (blue) superimposed on the real fMRI signals (green) were all averaged within the region of interest located in the Broadmann Area 4p. In the x-axes the time is expressed in TRs and the gray areas represent the duration of the tongue movements.

to the ground truth, even for very noisy data, corrupted with several kind of added noise, i.e. the block-type noise, the model and the additive noise. Furthermore, the comparison with the TA approach showed that we can obtain results that are closer to the ground truth, with a decreased computation time and without the need of defining lambda a priori. The lambda chosen in the TA approach has revealed to be greater than those obtained using our approach, this means the recovered solution obtained with the TA is too far from the data.

Results on real fMRI data highlighted a variability if we compare the task duration with the recovered activation. Therefore, it would be interesting to see results on an extended set of data, to observe if this pattern is confirmed or if it is due to the subject performance, and to better elucidate and confirm the validity of our approach.

In conclusion, our findings show that the joint use of the LARS algorithm and the L-curve for solving our optimization problem allowed us to choose the optimal lambda and its associate solution among all those outputted at once by the algorithm. In this way, we avoided a need of defining lambdas a priori, allowing to improve the state of the art and therefore brain dynamics recovery for future clinical application.

Publications Arising from this Contribution

Spatio-Temporal Deconvolution of the fMRI BOLD

Response

By supposing the brain activates in constant blocks, Karahano glu and colleagues [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF],

later revisited by Farouj et al. [START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF], developed a deconvolution approach which involves both spatial and temporal regularization called TA. However, these approaches split the optimization problem into two decoupled spatial and temporal regularizations. This increases the number of parameters to be set to 4, i.e. two regularization parameters and two weights in the context of the forward-backward splitting algorithm used to solve the optimization problem. This requires the solver to alternate between the constraints as explained in Section 2.2.

To overcome these limitations, starting from the idea that large image variations should be preserved as they occur during brain activation, whereas small variations should be smoothed to remove noise, in this chapter we propose a novel approach, based on partial differential equations (PDEs), named Paradigm-Free fMRI (PF-fMRI). The PF-fMRI applies a diffusion process whose diffusivity is steered by derivatives of the evolving image, in order to smooth the fMRI image and to simultaneously enhance important features such as spatial edges and temporal functional activations.

Let us start before with a physical background of diffusion processes.

Introduction to the Diffusion Process

In the book Anisotropic diffusion in image processing [START_REF] Weickert | Anisotropic diffusion in image processing[END_REF], Weickert described the diffusion in physics as a process that aims at maintaining the equilibrium between concentration variations without creating new mass, or destroying it as expressed in Fick's law

j = -D • ∇c (4.1)
where ∇c is the concentration gradient, j is the flux induced to compensate for this gradient and D is the diffusion tensor that describes the relation between ∇c and j. The case where j and ∇c are parallel is called isotropic, whereas if they are not parallel is called anisotropic.

The continuity equation allows us to write the observation that diffusion does not create or destroy mass

∂c ∂t = -divj (4.2)
where t denotes the time. So, if we combine the two equations (4.1) and (4.2) we finally obtain the diffusion equation

∂c ∂t = div(D∇c) (4.3)
that appears in many physical transport processes, for example, in the context of heat transfer where it is called heat equation. Eq. ( 4.3) can be linked and exploited in image processing if we consider the gray values of an image at a certain location as a concentration.

If the diffusion tensor is a function of the differential structure of the image that evolves with time, this leads to non-linear diffusion filters.

Image Regularization with Partial Differential Equations

Regularization methods have been enriched by the use of non-linear PDEs in several context for the last 30 years. Firstly applied to physics and fluid mechanics, it has been showed that non-linear PDEs allow smoothing the data while preserving large global features, such as discontinuities of the signal [START_REF] Tschumperle | Vector-valued image regularization with pdes: A common framework for different applications[END_REF], which can be found, for example, in image contours and corners [START_REF] Tschumperlé | Anisotropic diffusion partial differential equations for multichannel image regularization: Framework and applications[END_REF].

In their pioneering work presented in 1990, Perona and Malik [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] were the first to exploit the link between image regularization and anisotropic diffusion. They employed anisotropic diffusion PDEs for the restoration of noisy and blurred digital data to overcome the limitations associated to linear filtering approaches [START_REF] Tschumperle | Diffusion pdes on vector-valued images[END_REF]. This approach is rooted on the isotropic diffusion equation, i.e. heat flow, which has subsequently been extended to other theoretical contributions. Among them there are the anisotropic smoothing [START_REF] Weickert | Anisotropic diffusion in image processing[END_REF][START_REF] Sapiro | Geometric partial differential equations and image analysis[END_REF] and the PDEs-based gradient descent used to solve energy functionals minimizations [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Charbonnier | Deterministic edgepreserving regularization in computed imaging[END_REF][START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF][START_REF] Kornprobst | Mathematical problems in image processing[END_REF][START_REF] Kimmel | Images as embedded maps and minimal surfaces: movies, color, texture, and volumetric medical images[END_REF]. To date, PDEs-based regularization algorithm has been applied to 2-D scalar images [START_REF] Aubert | Mathematical problems in image processing: partial differential equations and the calculus of variations[END_REF][START_REF] Nielsen | Regularization, scale-space, and edge detection filters[END_REF][START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF][START_REF] Weickert | Anisotropic diffusion in image processing[END_REF] and vector-valued images [START_REF] Tschumperle | Diffusion pdes on vector-valued images[END_REF]. Interestingly, inspired by the physics of fluids, many authors assimilated the process of image regularization with the diffusion of chemical concentrations and propose to apply the diffusion PDE [START_REF] Tschumperle | Vector-valued image regularization with pdes: A common framework for different applications[END_REF][START_REF] Tschumperlé | Anisotropic diffusion partial differential equations for multichannel image regularization: Framework and applications[END_REF][START_REF] Weickert | Anisotropic diffusion in image processing[END_REF][START_REF] Tschumperle | Diffusion pdes on vector-valued images[END_REF][START_REF] Weickert | Coherence-enhancing diffusion of colour images[END_REF]]

∂I ∂t = div(D∇I) (4.4)
where I is the input image, ∇ is the gradient operator, t is the time, div(•) is the divergence and

D = λ 1 uu T + λ 2 vv T (4.5)
is the diffusion tensor of the image I., also called structure tensor [START_REF] Tschumperle | Vector-valued image regularization with pdes: A common framework for different applications[END_REF][START_REF] Förstner | A feature based correspondence algorithm for image matching[END_REF][START_REF] Förstner | A fast operator for detection and precise location of distinct points, corners and centres of circular features[END_REF].

The diffusion tensor D, which is a symmetric and positive definite matrix, has λ 1 and λ 2 as positive eigenvalues and u, v as corresponding orthogonal eigenvectors that drive the regularization process; the amount of diffusion in the directions u and v will be weighted by λ 1 and λ 2 . Furthermore, PDEs smooth the image at each step with a notion of scale-space [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF][START_REF] Nielsen | Regularization, scale-space, and edge detection filters[END_REF][START_REF] Lindeberg | Scale-space theory in computer vision[END_REF];

this means that, at each iteration, the image is smoothed and fine-scale properties, such as noise in the case of our interest, are gradually suppressed. At the same time, besides the regularization term, a fidelity term prevents the solution to be too far from the input data.

Using this approach actually means performing minimization of image variations as well as an image deconvolution.

The PF-fMRI: Theory

For the sake of simplicity the approach will be firstly explained in 3-D and then it will be where ∇I is the gradient of the image I, while I x , I y and I z are the partial derivatives of I with respect to x, y and z, respectively. By definition, D has three positive eigenvalues (λ 1 ≥ λ 2 ≥ λ 3 ≥ 0) and their associate three orthogonal eigenvectors (θ 1 , θ 2 and θ 3 ) explain the distribution and orientation of the gradient ∇I = (I x , I y , I z ) of I in a given neighborhood:

D = λ 1 θ 1 θ T 1 + λ 2 θ 2 θ T 2 + λ 3 θ 3 θ T 3 . (4.7)
A structure tensor allows to distinguish between an anisotropic and a isotropic diffusion. If λ 1 >> λ 2 , λ 3 it means that the gradient has a principal orientation (in this case θ 1 ) and the diffusion is anisotropic and can be represented with an ellipsoid (Figure 4.1). Whereas, if

λ 1 ≈ λ 2 ≈ λ 3 it
means that the gradient is not oriented in a main direction and θ 1 , θ 2 and θ 3 are eigenvectors of D with equal weight, i.e. isotropic diffusion that can be represented with a sphere.

A isotropic structure tensor could be then written as a weighted identity matrix

I d λ 1 θ 1 λ 2 θ 2 λ 3 θ 3
∂I ∂t = (1 -λ R ) H T (I 0 -HI) I 0 2 + λ R div( D∇I) div( D0 ∇I 0 ) 2 (4.9)
where the term on the left, ∂I/∂t, is the regularization flow, the first term on the right is the data fitting term, and the second term on the right minimizes image variations. λ R ∈ [0, 1] is the user-defined regularization parameter that represents a trade-off between the data fidelity term and the regularization term. Starting from the initial image I 0 , while increasing the scale variable t, the restored image I becomes more and more simplified, with respect to the image at small t. At the same time, no new structures are introduced in the image [START_REF] Kornprobst | Mathematical problems in image processing[END_REF].

Going more into the details of the data fidelity term F (I) = H T (I 0 -HI) I 0 2 (4.10) I 0 and I are the original and the regularized image respectively, I 0 2 in the denominator is the normalization factor, H is the HRF [START_REF] Khalidov | Activelets: Wavelets for sparse representation of hemodynamic responses[END_REF] operator (see Section 3.2.1) and H T is its transpose. The multiplication of H T with (I 0 -HI) corresponds to correlation and can be implemented via convolution with the time-reversed HRF. Note that the products with the HRF, which corresponds to a convolution with the time-reversed HRF, were computed only along the time dimension. The regularization term is defined as where Q contains the orthogonal eigenvectors (θ 1 , θ 2 , θ 3 , θ 4 ) of D and Λ contains their associated eigenvalues (λ 1 ≥ λ 2 ≥ λ 3 ≥ λ 4 ). We then recomputed the matrix

R(I) = div( D∇I) div( D0 ∇I 0 ) 2 . ( 4
D = Q ΛQ T (4.14)
such that for each voxel the highest eigenvalue

λ1 = exp - λ 2 1 max(Λ) 2 1 2σ 2 D (4.15)
was set according to a gaussian function with standard deviation σ D , such that if λ 1 is big, the current voxel may be located on a edge or activation and the diffusion tensor D is steered to be anisotropic, by setting λ1 λ2 , λ3 , λ4 . Since we aim at performing a smoothing only along the other three directions to smooth preferably along the coherence directions, the three eigenvalues λ 2 ≈ λ 3 ≈ λ 4 are indeed set to 1, therefore λ2 = λ3 = λ4 = 1. Whereas, if λ 1 is small, the diffusion will be isotropic in the four directions because λ1 ≈ 1 and λ2 = λ3 = λ4 = 1. Using the function in Eq. (4.15) corresponds to reassigning to each voxel different eigenvalues constituting the matrix Λ, before recomputing the operator D as in Eq. (4.14). In fact, if λ 1 /max(Λ) is the highest eigenvalue λ1 of the considered voxel will tend to zero. This steers the geometrical regularization to be anisotropic, because the smoothing will apply equally in the remaining three directions but it will be negligible in the perpendicular to the detected contour. Otherwise, if λ 1 /max(Λ) is small, the greatest eigenvalue λ1 will tend to 1 and this leads to an isotropic regularization almost in all the four directions (x, y, z, t). In both cases indeed λ2 , λ3 , λ4 are set to 1.

This procedure is applied for each voxel of the entire 4-D image such that at each iteration the image I computed in Eq. (4.9) is gradually removed from the image obtained at the previous iteration. In this way, supposing the brain activates in constant blocks, we were able to smooth the image simultaneously in space and time. Furthermore we were able to keep large image variations occurring during a brain activation or a spatial edge, whereas small variations, corresponding to noise, were smoothed and gradually removed, while conserving and enhancing coherent structures of the fMRI image.

In the following sessions we will explain how we tested and validated the PF-fMRI on phantom and real fMRI data.

Simulation of a Whole Brain fMRI Image

Similarly to what we did in Section 3.5, to simulate the whole brain fMRI image, we modeled the activity-inducing signal as a boxcar function u(t) for each voxel v.

We added noise to u(t) representing the random intrinsic electrical fluctuations within neuronal networks. To do this, we corrupted the activity-inducing signal u(t) with an additive random Gaussian noise with zero mean and standard deviation σ m that we called

"model noise" m u n (v, t) = u(v, t) + m .
We modeled the activity-related signal x(t), consequent to the brain activation as the convolution of the activity-inducing signal with the HRF, h(t) [START_REF] Khalidov | Activelets: Wavelets for sparse representation of hemodynamic responses[END_REF]:

x(v, t) = u n (v, t) * h(t).
We finally added noise to x(t) to model heart rate, respiratory fluctuations, motion artifacts, thermal noise and scanner drifts [START_REF] Lund | Non-white noise in fmri: does modelling have an impact?[END_REF], thus obtaining the acquired fMRI signals

y(v, t) = x(v, t) + a = u n (v, t) * h(t) + a
where a is the additive random Gaussian noise with zero mean and and standard deviation σ a .

Validation on phantom fMRI data

In this section we will test and validate the PF-fMRI approach on simulated fMRI images corrupted with different amount of noise.

Methods

As a proof of concept, to test and validate the PF-fMRI, we scaled a 3-D activation map computed with FSL Possum in the range [0,3], with a 2-mm isotropic resolution (Figure

4.2.a).

We multiplied it by a piece-wise constant signal u(t) of 100 s, with one onset of 40 s, from 20 s to 60 s (Figure 4.2.b). After that, the image was corrupted with model noise and subsequently each time course was convolved with the HRF [START_REF] Khalidov | Activelets: Wavelets for sparse representation of hemodynamic responses[END_REF]. Finally we added gaussian noise thus simulating the fMRI time-courses y(t). The TR was set to 1 s.

We tested the PF-fMRI approach on several simulated images obtained by adding different amount of noise for each experiment.

We regularized the whole image using our approach as showed in Section 4.4, and we recovered the voxel-wise activity-inducing signals û(t). Finally, to evaluate the results, we computed the root of the MSE and STD as well as the Pearson correlation (r) and its STD between u(t) and û(t) averaged among the voxels belonging to the GM. We compared our results with those obtained using the spatio-temporal deconvolution proposed by Farouj et al. [START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF], implemented in the TA toolbox1 . To better investigate the performance of the PF-fMRI, additional tests were executed given the shape of the red curve in Figure 4.6, that are almost constant and close to 1 with a variance that is very low and constant. We would expect a trend where the correlation increases with the pSNR and the variance instead decreases with the pSNR, as it is the case for TA. In Figure 4.7 we report experiments for more corrupted images where we applied the PF-fMRI. In the top row results, we applied the PF-fMRI on the simulated fMRI image y, after corrupting the ground truth image u both with model noise ( m ) and additive noise ( a ) with a standard deviation of 2, the resulting pSNR was 0 dB; in the bottom row, we applied the PF-fMRI after corrupting the ground truth image u both with model noise ( m ) and additive noise ( a )

Results

with a standard deviation of 5, leading to a negative pSNR (pSNR = -6.51 dB). Note that, because images were too noisy, the colorbar between the spatial maps in panel (a) and the time series in panel (b) are mismatching; otherwise we would not have been able to visualize the spatial maps. Figure 4.8 shows the roots of MSEs, and relative standard deviations, between the ground truth u(t) and the recovered signals ûPF-f MRI (t) averaged across the voxels belonging to the GM, for all the experiments that we run in the thesis, represented in On the light of results emerged with these experiments, we show that for very noisy data our results become more sensitive to noise. On one hand, we observe that the values of the MSEs are higher for very noisy data, and stabilize at a pSNR 0 dB. On the other hand, as expected, the Pearson correlation coefficients become lower for very noisy data, and increase with the pSNR whereas the standard deviations decrease with the pSNRs. In conclusion, we show the expected two trends: the errors decreases while the pSNR increase and the Pearson correlations increase with the pSNR. Nonetheless, if we compare these results to those of Figure 4.5 and 4.6, and in particular observing the values of the errors computed on the results obtained using the TA approach, that we run using the provided toolbox without setting any parameter, we observe that the errors obtained using the PF-fMRI on very noisy data are still lower than the ones obtained using the TA. As for the Pearson correlation coefficients, we observe that the values obtained using the PF-fMRI are still greater than the ones observed with the TA for very corrupted acquired images. In this regard, the parameter λ R is adjusted with respect to noise, and this compensate for the presence of it. From top left: spatial maps of a slice of the simulated functional MRI (fMRI) image y, ground truth activation u, recovered activation maps obtained using the Total Activation (TA) approach ( ûTA ) and the PF-fMRI ( ûPFf MRI ). The map obtained using TA ( ûTA ) had lower amplitude compared with the ground truth, nonetheless the scale between the values were kept. The map obtained using PF-fMRI ( ûPFf MRI ) showed amplitude comparable to the ground truth, the scale between the values were kept. The time series on the bottom shows that the recovered signal ûPF-f MRI (t) (red) was closer to the ground truth u(t) (black) compared to ûTA (t) (blue).

PF-fMRI TA

Validation on real fMRI data

In this section we will test and validate the PF-fMRI approach on real data. Even if the PF-fMRI is conceived to be applied to rs-fMRI data, where no experimental paradigm is given, we firstly used task-fMRI images to be able to have a ground truth, i.e. the timing of the tasks, to which to refer to test the validity of our approach. We also applied the PF-fMRI on real rs-fMRI data as a proof of concept.

Methods

Validation on task-fMRI Data

The study was conducted on the motor task-fMRI data from 51 subjects from the HCP dataset [START_REF] Van Essen | The wu-minn human connectome project: an overview[END_REF]. The data acquisition parameters and preprocessing steps are reported in Section 3.8.1.

The task involves, for known timing and duration, the movement of the left foot and left hand, the right foot and right hand and the tongue; all the tasks were following a visual cue of 3 s, and their timing and durations are shown in After applying the PF-fMRI on the whole brain images of each subject, we obtained the reconstructed activity-inducing signals û(t), without prior knowledge on the onset times of the evoked stimuli. The regularization parameter λ R was set experimentally to 0.9997, σ G was set to 1, sigma D was set to 0.2 and we computed up to 40 iterations.

To highlight the ability of the PF-fMRI to recover brain activations without knowledge of the experimental paradigm, we qualitatively compared brain regions recovered using the PF-fMRI to those recovered using the GLM as implemented in the FSL library. To estimate the results given by the PF-fMRI, we computed the voxel-wise correlation maps, by estimating the Pearson correlation coefficient (r) between the recovered activations and the 5 tasks (Table 4.1). The tasks, i.e. the tongue, the right and left hand, the right and left foot, were simulated as a piece-wise constant signal with unit amplitude when the subject is performing the task and zeros elsewhere. As for the GLM, we included all the 5 tasks (Table 4.1) in a design matrix and we estimated the regressor weights using FSL. Results showing differences and similarities of both approaches were qualitatively assessed.

Subsequently, we quantitatively compared the results obtained using the PF-fMRI with the ones given by the TA approach. To do this, we firstly defined 4 ROIs located in brain regions which are involved in the 5 considered tasks. The ROI related to the tongue was bilater whereas for the hands and the feet we defined two separate ROIs: one for the left side and one the right side of the brain. To do this, we started by defining the ROIs from the work proposed by Roux and colleagues in 2018 [START_REF] Roux | Functional architecture of the somatosensory homunculus detected by electrostimulation[END_REF], who mapped the somatosensory homunculus MNI coordinates using the electrostimulation. For each coordinate center, we built a spherical 3 mm-radius ROI and we grouped the multiple ROIs related to each task into a unique ROI. Coordinates' centers are shown in After we had defined the ROIs, similarly to the previous comparison between the PF-fMRI method and the GLM, we computed the whole-brain voxel-wise correlation maps between the time course related to each task and the recovered activity-inducing signals û(t) obtained using the PF-fMRI and the TA approaches. After that, for each subject we firstly computed the average of the Pearson correlation coefficients (r) inside GM-masked ROIs, and then we calculated the mean and the standard deviation of these averaged correlation values among the 51 subjects belonging to the sample data.

Furthermore, in order to prove that the PF-fMRI is able to differentiate between a region that is activated and one that is not, the time courses û(t) of one representative subject (100307), were averaged in two ROIs of 6 × 6 × 6 mm 3 : one which is expected to be active during the task, and one located in a brain area that is not involved in the task. We selected the task related to the tongue, and we chose one ROI centered in the Brodmann Area (BA) 4p (rBA4p; MNI coordinates: 62, -14, 30) which is activated during a tongue motor task, and another centered in the primary auditory cortex (TE1.2; MNI coordinates: 56, 4, 10) [START_REF] Kiviniemi | Independent component analysis of nondeterministic fmri signal sources[END_REF],

that is not involved in the tongue movement. After that, Pearson correlation coefficients (r)

were firstly computed between the tongue activation and the recovered û(t) for each voxel, and secondly averaged among the voxels belonging to the two GM-masked ROIs. The tongue activation was again simulated as a piece-wise constant signal with unit amplitude when the subject is performing the task and zeros elsewhere. We compared results obtained using the PF-fMRI with those obtained using the TA toolbox.

Application on Resting-State fMRI Data

Finally, as a proof of concept, we applied the PF-fMRI on the rs-fMRI image of subject 100307 from the HCP database. Data were acquired with a SIEMENS MAGNETOM Connectome Syngo MR D11 using a gradient-echo EPI sequence (TR = 720 ms; TE = 33.1 ms; FA = 52°; FOV = 208 × 180 mm; slice thickness 2.0 mm; number of slices = 72; 2.0 mm isotropic voxels).

The subject was asked to lay in the scanner without thinking about anything in particular.

The number of acquired frames were 1200 and the duration of the acquisition was 14:33 min. In the case of rs-fMRI data, the task paradigm is unavailable since the subject does not perform any task in the scanner. The data underwent the same minimal preprocessing of the task fMRI data as proposed in the HCP pipeline [START_REF] Glasser | The minimal preprocessing pipelines for the human connectome project[END_REF]. Besides this, time series were detrended to remove linear drifts and normalized to zero mean and unit standard deviation.

After applying the PF-fMRI on the whole rs-fMRI image, we chose a set of a-priori ROIs representing the Default Mode Network (DMN). To do this we used as reference the BAs corresponding to each anatomical area of the DMN, as specified by Buckner

Results

As for the real data analyses, and specifically the comparison between the PF-fMRI approach and the GLM, we showed that the correlation maps related to each task computed with the PF-fMRI approach were well overlapped to the values of the regressors coefficients obtained using the GLM one as shown for one illustrative subject (100307) in Figure 4.10. The GLM shows results which follows the GM, while activations found with the PF-fMRI, which were performed across the whole brain, and not masked with the GM mask, cover also voxels belonging to the WM. Interestingly, the found activations overlap the areas found to be active in the motor Homunculus brain [START_REF] Penfield | Somatic motor and sensory representation in the cerebral cortex of man as studied by electrical stimulation[END_REF].

The comparison with the GLM was proposed in this thesis because one of the aim of the thesis was to find an alternative to the GLM, if the user wants to analyze data where the experimental paradigm is not available, as it is the case, for example, for resting-state fMRI data. The GLM aims at modelling at each voxel the observed fMRI time series as a linear combination of explanatory variables (or regressors) plus an error term. The beta maps obtained with the toolbox FSL tell how much each brain voxel contributed to the given task, when solving a GLM. In other words, the beta maps tell how much the time series of that voxel fits the given task, given all the tasks that the subject performed as regressors. Even if the estimation of the beta maps passes through a convolution with the HRF, that is applied to match the input function, i.e. the stimulus waveform, and the output fMRI data function, the GLM with FSL takes as input the tasks time courses, and the aim is to give a measure of "similarity" between the given task and the voxel's time course.

Quantitative comparison between the activity-inducing signal recovered using the PF-fMRI approach and the TA are shown in the number of iteration, until it converges after 25 iterations for the hands, 5 iterations for the feet and about 35 iterations for the tongue. Moreover, starting from the first iteration, we

show higher correlation values compared to the ones obtained using the TA method.

In particular for the comparison between the PF-fMRI and the TA, Figure 4.12 shows the reconstructed signal û(t) and the correlations values given by the PF-fMRI and the TA for a single subject. We show a clear difference between the correlation values estimated in the area involved in the task and the one which is not involved. In fact, we showed a higher correlation between the tongue activation (Figure 4.12.b) and the recovered activation û(t) in the ROI rBA4p, which we expect to be involved in the motor task, while a low correlation is shown with the recovered signal in the ROI rTE12 which is not involved in the task.

Whereas, the TA approach was not able to clearly distinguish between an active and a not active region since it showed low correlation values for both ROIs (Figure 4.12.a and 4.12.c).

As for the ranges in the colorbars, in The blue curves are related to the TA approach, while the red one to the PF-fMRI approach. TR: repetition time.

Discussion and Contributions of this Chapter

In this chapter we introduced our second main contribution, an innovative approach that aims at the recovery of brain functional activations, in space and in time, without the necessity of a priori knowledge of the paradigm, namely the PF-fMRI. The PF-fMRI is conceived to be exploited in cases where the subject is not able to perform a task, for example because affected by some neurological disease.

The approach proposed in this chapter was demonstrated to be successfully able to recover and locate the fMRI activations via a geometric regularization. The interesting point is that this is done on the whole 4-D fMRI image, meaning that we had to solve a unique problem, coupling the space and the time dimensions. Other approaches have been proposed so far to analyse fMRI data. Among them: (i) the GLM, that fit a linear model to the fMRI time series, but it assumes prior knowledge of the tasks occurrences [START_REF] Friston | Event-related fmri: characterizing differential responses[END_REF]; (ii) deconvolution methods, which are used to uncover the brain activations from the BOLD response without prior information on the underlying activity [START_REF] Gaudes | Detection and characterization of single-trial fmri bold responses: Paradigm free mapping[END_REF][START_REF] Gaudes | Paradigm free mapping with sparse regression automatically detects singletrial functional magnetic resonance imaging blood oxygenation level dependent responses[END_REF][START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF][START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF]. In particular, the deconvolution approach proposed in the TA splits the optimization problem into a spatial and temporal regularization problems, meaning that the user has to specify two regularization parameters (i.e. λ S and λ T ) and the two weights used for the forward-backward splitting algorithm (i.e. ω S and ω T ) in order to have a solution which is given by a weighted sum of the two separates regularization processes.

In contrast to these, the PF-fMRI overcomes several limitations found in the previous literature. When comparing the regions recovered using the GLM and PF-fMRI, we noted overall very good agreement between the two methods. It should again be emphasized that, while the GLM requires knowledge of the experimental paradigm, PF-fMRI does not.

These results highlight that PF-fMRI can be used to recover brain activity in the absence of an experimental paradigm, such as rs-fMRI. When comparing correlation maps obtained for PF-fMRI and TA, correlation values obtained with PF-fMRI were significantly higher than those obtained with the TA suggesting an improved recovery of brain activity. We also demostrated that the PF-fMRI can be used for different purposes, e.g., to recover the brain activations and their location in a task experimental paradigm as well as in a rs-fMRI study, where the subject is asked not to perform any task while lying in the MRI scanner.

In particular, our recovered signals from rs-fMRI data could be further employed in the approach proposed by Karahano glu and Van De Ville [START_REF] Karahano | Transient brain activity disentangles fmri resting-state dynamics in terms of spatially and temporally overlapping networks[END_REF] to investigate the dynamics of resting-state networks and reveal transients in spontaneous brain activity, as described in Section 2.2.

Publications Arising from this Contribution

Application of the PF-fMRI Approach

The PF-fMRI approach described in this chapter has been applied also in a different context as reported in the following paper:

Software Contribution

In the context of our second contribution, we created a Python package, called Paradigm-Free Functional MRI (PF-fMRI), that implements all we had described in this chapter.

The Python package is available in a GitLab repository dedicated to the Computation Brain Connectivity Mapping (CoBCom) Project. The package is easy to use because it simply requires the acquired noisy or preprocessed fMRI image and returns as output the regularized 4-D fMRI data.

Part V Conclusion

Chapter 5

Concluding Remarks and Open

Problems

In this dissertation, we have proposed two main contributions:

1. A novel fMRI deconvolution approach via Temporal Regularization using a LASSO model and the LARS algorithm to recover brain functional activations without necessity of a priori information of the timing and the duration of the underlying activation. In what follows, we discuss the main achievements of both contributions and we explore the potential future research directions, and identify possible extensions of our method.

fMRI Deconvolution via Temporal Regularization using a LASSO model and the LARS algorithm

In our first main contribution, we provided a temporal regularized deconvolution approach to deconvolve the BOLD response given the fMRI noisy time-courses, thus recovering the brain functional activations. This was done without the necessity of a priori information on the experimental paradigm, as it was for the GLM approach.

To solve the inverse problem, we proposed an approach that combines the use of the LARS algorithm and the L-curve. As a first point, the LARS algorithm gives as output all important lambdas and their associated solutions. Secondly, the combined approach of LARS and Lcurve allowed us to choose the optimal regularization parameter lambda and its associated solution among all those outputted by the algorithm. In this way we overcame the need of choosing lambda a priori.

Our approach was tested and validated both on synthetic data and on real task-fMRI data, that provided us of a ground truth to which to refer in order to assess the goodness of results. The comparison with the state-of-the-art approach, namely the TA, showed that we can obtain results which are closer to the ground truth, with a decreased computation time and without the need of defining lambda a priori. Our results suggests that using our approach improves the brain dynamics recovery for future clinical application.

Nevertheless, the LARS algorithm proposed in this contribution, is limited by the fact that it cannot handle a whole fMRI image, made of almost a million of voxels because the problem becomes too heavy in terms of computation.

A Paradigm-Free Regularization Approach to Recover Brain Activations from fMRI Data

In the second contribution, as a natural extension of the previous one and to overcome its limitations, we aimed at regularize the whole fMRI image exploiting both the temporal and the spatial information contained in the image itself. Therefore, we proposed and described an innovative method to analyse fMRI images thus recovering the location and the occurrence in time of the functional brain activations.

The proposed approach, namely the PF-fMRI, aims at geometrically regularize the fMRI image such that: it saves and highlights big variations of the image as they are present at the occurrence of a brain activation or in the presence of a spatial edge with respect to little image variations, which instead are removed to reduce noise. To do this, we used the PDEs in a iterative algorithm and exploited the 4-D image structure tensor, that defines the directions of the gradient in the neighborhood of a voxel and directs towards an anisotropic or isotropic regularization. This gradient contains all the four principal directions of the image, which are: the 3-D space and the time direction, which suggests that the whole 4-D fMRI image was smoothed contemporaneously in space and time at once.

The PF-fMRI was tested and validated both on phantom and real task-fMRI data. In fact, even if the PF-fMRI is conceived to be applied to rs-fMRI data, where no experimental paradigm is given, we firstly used task-fMRI images to be able to have a ground truth to which to refer to assess the performances of our approach. When compared to the state-of-the-art approaches, the PF-fMRI showed a very good agreement with the results obtained using the GLM methods, a widespread and well-known approach that instead requires a priori information on the task paradigm. This confirmed that the PF-fMRI can be successfully applied to uncover brain activations, in the absence of a known and given experimental paradigm, such as in the case of rs-fMRI data. We also showed significantly better results when the PF-fMRI was compared to the TA, meaning that we improved the recovery of brain functional activitions. It must be stressed again, that the PF-fMRI is in contrast with the deconvolution approaches, namely the TA, because it treats the image as a whole. The deconvolution approaches instead uses to split the problem into a spatial and temporal regularization problems, meaning that the user has to specify two regularization parameters (i.e. lambdas) and two weights used to have a unique solution which is given by a weighted sum of the results given by the two separates regularization processes. Finally, as a proof of concept, we applied the PF-fMRI on rs-fMRI data, while the subject is not exposed to stimuli and is asked not to perform any task or not to think about anything in particular during the MRI acquisition. This represents a powerful tool in the case the subject is not able to perform tasks because of a neurological disease.

The PF-fMRI is proposed in a Python package available in a GitLab repository dedicated to the Computation Brain Connectivity Mapping (CoBCom) Project. The package was made easy to use, because it requires the fMRI measurements and returns as output the regularized 4-D fMRI data. Nonetheless, at the same time it was made such that it is possible for the user to easily navigate into it and manipulate the code for further improvements.

The PF-fMRI has also already been applied in a multi-modal approach that aims at inferring the information flow in the white matter without a manual selection of the connections of interest. The fMRI recovered activations using the PF-fMRI were used to locate active brain regions to be used as priors on cortical regions, to filter structural connections (Appendix B).

Future Perspectives and Open Problems

In this dissertation we proposed and validated two approaches to deconvolve the fMRI images and recover the brain activity from fMRI signals without prior knowledge. The first one exploited the information contained in the temporal structure of the fMRI data; the second one, as a logic extension of the first one, exploited all the fMRI image, that contains both a temporal and a spatial information.

Our findings show that both contributions improved the state-of-the-art approaches. In particular, the PF-fMRI have marked a step forward in the recovery of brain activation because it treats the fMRI image as a whole. Indeed, the PF-fMRI enabled us to solve a unique problem, coupling the spatial and the temporal dimension and to recover brain activation overlapping the ones obtained with the GLM and improving those obtained with the TA. Furthermore, in contrast with the ICA approach, both proposed methods allow a voxel-wise recovery of brain activations and incorporate the hemodynamic model.

Nevertheless, there are aspects of our method that may need further research. For example, the PF-fMRI is limited by the fact that it does not give a very sharp piece-wise solution representing the activity-inducing signal, as output, and this could be due to the choice of the regularization parameter.

The PF-fMRI could be exploited in the analysis of fMRI images, since it provides regularized time-courses as well as in more complex investigations. For example, because it provides a voxel-wise map of recovered brain activations, it could be used to explore the dynamics of resting-state networks as well as to investigate possible WM fMRI activations, that is an emerging debated topic on the neuroimaging field.

Future works could include into the solution of the problem the dMRI data. The streamlines would define a more complex neighborhood, not only given by surrounding voxels, but also by voxels which are segregated in space but structurally connected.

In conclusion, in this dissertation we opened a new channel for the analysis of rs-fMRI data, to improve the recovery of brain activations and their dynamics and to be used for the investigation in future clinical application.

In order now to design an exponential accumulation operator once computed the ztransformed H(z), we define the accumulation property of the z-transform n ∑ k=-∞

x(k) = 1

(1z -1 ) X(z).

(x * u)(n) = +∞ ∑ k=-∞ x(k)u(n -k) = n ∑ k=-∞ x(k)
The z-transform of a convolution is

Z {x(n) * u(n)} = Z {x(n)} • Z {u(n)} therefore Z n ∑ k=-∞ x(k) = X(z) • z (z -1) = 1 1 -z -1 • X(z).
We accumulated H(z) twice and we obtained the transfer function

I α+ (z) = S e -α z -1
(1e -α z -1 ) 2 1

(1z -1 )

1 (1 -z -1 )
= S e -α z -1

(1e -α z -1 ) 2 1

(1z -1 ) 2 .

(A.10)

The anticausal part of the filter is then derived from I α+ (z) as

1. 1

 1 Neurons structure and process of synaptic transmission. The figure shows how the neurons are composed and connected to each other. The cell body contains the nucleus and cytoplasm. The axon extends from the cell body and gives rise to small branches before ending at nerve terminals. Dendrites extend from the neuron cell body and carry the information between neurons. Synapses are the points where the "communications" between neurons take place. Figure taken from [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Conventional terminology which indicates the different orientations representing the different brain views: in surface (A), section (B) and connectional anatomy (C). Surface neuroanatomy is the one describing the appearance of grooves (sulci) and folds (convolutions or gyri) of the brain. Sectional neuroanatomy illustrates the cortical and subcortical structures and their relationship, commonly described along the axial, coronal and sagittal plane. The connectional neuroanatomy desribes the connecting fiber tracts' origin, paths and terminations. Figure taken from Catani and Thiebaut de Schotten [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Brain lobes of a dorsolateral (A) and medial (B) surface of the left cerebral hemisphere. Solid lines show the main sulci dividing the lobes, while dashed lines represent arbitrary lines of separations which are following a sulcus. The gyri of each lobe are also illustrated. Figure taken from Catani and Thiebaut de Schotten [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Longitudinal and transverse relaxation. T1 corresponds to the time in which the longitudinal component regain a value which is the 37% lower of its equilibrium value; T2 corresponds to the time in which the transverse component is reduced to the 37% of it. Figure taken from Pizzolato [3]. . . . . 1.5 Examples of MRI brain scans in the axial plane: T1-weighted (left) and T2weighted images (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Echo Planar Imaging (EPI) Sequence -(A) EPI sequence with phase-encoding along y direction using blipped gradient pulses. (B) Trajectory of the EPI sequence in the k-space. The different colors of the trajectories refer to the gradients G x and G y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Schematic illustration of the generation of the blood-oxygen-level-dependent (BOLD) signal. Figure adapted from Iannetti and Wise [4]. . . . . . . . . . . . 1.8 Hemodynamic response function (HRF) to a stimulus with a short duration, illustrated with the red bar. The peak related to the blood-oxygen-leveldependent (BOLD) effect starts to increase after approximately 3 s from the stimulus starting point. Figure taken from Amaro and Burker [5]. . . . . . . . xiv 1.9 Example of functional MRI (fMRI) brain scans. Each image represent a scan of fMRI data, consecutively acquired in time (t). The images are represented, from left to right, in their sagittal, coronal and axial planes. . . . . . . . . . . . 2.1 The approach of the general linear model (GLM). Each time series of fMRI data (Y) is a linear combination of regressors, composed both of task-related regressors (green) and nuisance regressors (light green) weighted by β, and the error . The goal of the GLM is to minimize the error thus estimating the parameters weights β. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 fMRI signal model. From left, u s (i, t) is the spike train (innovation signal)

2. 3

 3 [START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF] innovative-driven co-activation patterns (iCAPs) ordered with respect to their occurrence during a resting-state fMRI acquisition. The iCAP 1 contains auditory regions; the iCAP 2 includes regions of the fronto-parietal attention network; the iCAP 3 as well as the iCAP 4 cover the primary and secondary visual areas; the iCAPs 5 reveals the precuneus, the posterior cingulate cortex and the thalamus; the iCAP 6 represents the visuospatial/dorsal attention network; the iCAP 7 covers the motor network and the medial fromtal gyrus; the iCAP 8 corresponds to posterior part of the DMN; the iCAP 9 includes the anterior executive network; the iCAP 10 shows again a posterior part of the DMN; the iCAP 11 reveals the anterior salience network; the iCAP 12 is composed by the combination of different regions which are located in the limbic and the subcortical area, superior and middle temporal and occipital gyrus; the iCAP 13 contains the frontal gyrus, the anterior cingulate cortex and caudate. Figure adapted from Karahano glu and Van De Ville [7]. . . . . . 2.4 Seed-based method -To examine the time series and extract meaningful information from fMRI data, i.e. the functional connectivity level between a seed voxel and another brain voxel i, the correlation between the two time series showed on the bottom right of the figure are computed. Highly correlated time courses reflect a high level of functional connectivity. To map all the functional connections of the seed and obtain the map showed at the bottom left of the figure, a voxel-wise correlation between the seed voxel's

2 . 6 49 3. 2

 26492 Figure taken from Hu and Zeng [9]. . . . . . . . . . . . . . . 2.5 Matrix representation of the spatial indepenedent component analysis (spatial-ICA). In spatial-ICA, the algorithm try to estimate spatially independent components with related time series. . . . . . . . . . . . . . . . . xv Different probabilistic independent component analysis (probabilisti-ICA)estimated resting state outputs. Estimated from a group of 10 subjects, the eight spatial maps, coregistered and superimposed to the Montreal Neurological Institute (MNI) template, shows the sagittal, coronal and axial views of different components associated with low-frequency resting-state patterns. R: right; L: left. Figure adapted from Beckmann and colleagues [10]. 39 3.1 Representation of the 4-D functional MRI (fMRI) data. A fMRI image is composed by a set of 3-D volumes recorded over time, thus leading to a 4-D structure. TR = repetition time. Figure adapted from [11]. . . . . . . . . . . From top left, following the black lines connecting the different plots, s(t) is the spike train (innovation signal) which induces the activations.

64 xvi 3 . 7

 6437 CV. With respect to the parameter α, results obtained via the LOO-CV are improved by increasing it, instead the ones obtained with L-curve are more stable with respect to it. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Each plot in the figure show the root mean square errors (RMSE) and the roots of standard deviations computed between the ground truth activation and the estimated one using the Least Angle Regression combined with the L-curve (in red) and with the leave-one-out cross-validation (LOO-CV, in light blue).

  3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 3-D activation map obtained from an auditory task superimposed on the standard Montreal Neurological Institute (MNI) brain. The three images, from left to right, represent the axial, the sagittal and the coronal view. Voxels' intensity is ranged between 0 and 3. . . . . . . . . . . . . . . . . . . . . . . . . xvii 3.11 Piece-wise constant activity-inducing signals simulated for validation. The activation on the left, A, is composed by four blocks of activations with

  is expressed in TRs and the gray areas represent the duration of the tongue movements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Ellipsoidal representation of the 3-D diffusion tensor. . . . . . . . . . . . . . . 4.2 Ground truth for the functional MRI (fMRI) simulated data: (a) activation map. (b) Simulated activation u(t), with a repetition time (TR) of 1 s. . . . . .

xviii 4 . 5

 45 The graph shows, for different peak-SNRs, the roots of the mean square errors (MSE) and standard deviation (STD) between u(t) and û(t) averaged among the voxels belonging to the gray matter. . . . . . . . . . . . . . . . . . . . . . .4.6The graph shows, for different peak-SNRs (pSNRs), the Pearson correlation coefficient computed between u(t) and û(t) and averaged among the voxels belonging to the gray matter and their standard deviation. (µ r : mean correlation coefficient; σ r : standard deviation of the correlation coefficients.) . 4.7 (a) From left to right: spatial maps of the simulated functional MRI (fMRI) image y, ground truth activation u, recovered activation using our approach ( ûPFf MRI ). Each row corresponds to a different peak-SNR (pSNR): 0 dB, -6.51 dB, from the top to the bottom. (b) Reconstructed time series obtained with our approach ûPF-f MRI (t) (red) superimposed on the ground truth activation u(t) (black) and the original fMRI signal y(t) (green). . . . . . . . . 4.8 The graph shows, for different peak-SNRs (pSNR), the roots of the mean square errors (MSEs) and standard deviation (STD) between the ground truth signal u(t) and the recovered signals ûPF-f MRI (t) averaged across the voxels belonging to the gray matter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 The graph shows, for different peak-SNRs (pSNRs), the Pearson correlation coefficient computed between the ground truth signal u(t) and the recovered activation ûPFf MRI (t) averaged over the voxels belonging to the gray matter and their standard deviation. (µ r : mean correlation coefficient; σ r : standard deviation) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Qualitative comparison between the General Linear Model (GLM) and our approach (PF-fMRI). On the left column, in a blue-lightblue color-map, superimposed to the standard Montreal Neurological Institute (MNI) brain, the β-regressors map obtained using the GLM implemented in the FMRIB Software Library (FSL) tool. On the right column, in a red-yellow colormap, the whole-brain voxel-wise correlation maps obtained using the PF-fMRI superimposed to the standard MNI brain. The Pearson correlation was computed voxel-wise among the whole brain, between the reconstructed activity inducing signals û(t) and the five motor tasks simulated as piece-wise constant signals with ones in the time points where the subject is executing a task and zeros elsewhere. The values r of the correlations are indicated by the color-bars. Each row corresponds to a specific motor task, from top to bottom: the tongue, the right and left hand, and the right and left foot. A: anterior; P: posterior; S:superior; I: inferior; R: right; L: left. . . . . . . . . . . . . . . . . . . 4.11 Barplots of the mean (µ r ) ± standard deviations (σ r ) of the Pearson correlation coefficients (r) computed on the sample data of 51 subjects in 5 ROIs related to the tasks of the left and right hand, the tongue end the left and right foot. For each task, the bars in red represents the results using the PF-fMRI for an increasing number of iterations (from 1 to 40). The blue bar represents the results obtained with the Total Activation (TA) approach. The black lines are the standard deviations (σ r ). (lHAND: left hand; rHAND: right hand; lFOOT: left foot; rFOOT: right foot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix 4.12 a) Reconstructed signals û(t) obtained with our approach (PF-fMRI, red) and the total activation tool (TA, blue) superimposed on the real acquired fMRI signals (green). The ground truth in this case corresponds to the task (black), simulated as piece-wise constant signals with ones where the subject was asked to perform the task and zeros elsewhere. The plot on the top is related to the region of interest (ROI) located on the Brodmann Area 4p (rBA4p), the plot on the bottom is associated to the ROI positioned on the primary auditory cortex (rTE1.2). All the signals were averaged across the voxels belonging to the gray-matter (GM)-masked ROIs. The grey areas represent the occurrence and the duration of the tongue movements. b) Simulated tongue activation. c) Mean Pearson correlation coefficients (µ) and their associated standard deviations (σ) computed between the tongue activation and the recovered signals û(t) averaged across the voxels belonging to the GM-masked ROIs (rBA4p on the left, rTE1.2 on the right). The blue curves are related to the TA approach, while the red one to the PF-fMRI approach. TR: repetition time. . . 4.13 Each row in the figure shows for different iterations of the PF-fMRI approach, specifically the 5 th , the 15 th and the 45 th iterations, the reconstructed signals (in red) obtained with our approach superimposed on the real resting state fMRI signals (in green). The activity-inducing signals û(t) recovered via the PF-fMRI were averaged in a ROI located on the default mode network and masked using the gray matter mask. . . . . . . . . . . . . . . . . . . . . . . . . A.1 Function h(n) defined in Eq. (A.1), for S = -1 and α = 1. . . . . . . . . . . . . A.2 Function g(n) defined in Eq. (A.15), for S = -1 and α = 1. . . . . . . . . . . . A.3 Function I α (n) defined in Eq. (A.13), for S = -1 and α = 1. . . . . . . . . . . A.4 Function g(n) defined in Eq. (A.15), for S = -1 and α = 1. The plot shows that

  n). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.5 The plot shows for different values of α the function I α (n) defined in Eq. (A.13), normalized by the factor S defined in Eq. (A.35). The graphs show that for greater values of α, the filter is sharper, whereas if α is smaller, the filter has a smoother shape. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.1 Illustration of the Desikan-Killiany atlas (A), functional MRI priors (B), and thresholded functional MRI priors with activation above 0.83 (C). The remaining regions are those used to select connections introduced into the proposed Bayesian network which are illustrated in D. . . . . . . . . . . . . . xxiList of Tables1.1 T 1 and T 2 relaxation times for brain tissues, gray matter and white matter, at 3 Tesla [14]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Variables involved in the hemodynamic model [15, 16]. . . . . . . . . . . . . . 3.2 Constants involved in the hemodynamic model [15, 16]. . . . . . . . . . . . . .3.3 Experiments set for testing the L-curve and the LOO-CV to compare the two approaches for the selection of optimal lambda. The first column indicates the experiment number, the second column the standard deviation (σ m ) of the gaussian random noise added as model noise, the third column represents the peak signal to noise ration (pSNR m ) computed as in (3.28) in dB. The fourth column reports the σ a for the additive gaussian noise and the fifth column the pSNR a estimated as in (3.29) expressed in dB. . . . . . . . . . . . . . . . . 3.4 Experiments set for testing our temporal regularized deconvolution on phantom fMRI image. The first column indicates the experiment number. The second and third columns report the standard deviation of the model noise (σ m ) and the additive gaussian noise (σ a ), respectively. The fourth and fifth columns show the pSNR computed as in (3.30), expressed in dB, for activations A and B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Summary of rooted MSEs and STDs obtained for phantom fMRI data. . . . . 3.6 Tongue task paradigm of HCP data. Note that the time for activation is expressed in seconds and the repetition time for the HCP data is 0.72s. . . . . 4.1 Motor task paradigm of HCP data. The time is expressed in seconds and the repetition time for the HCP data is 0.72s. . . . . . . . . . . . . . . . . . . . . . . 4.2 MNI coordinates centers of the brain areas found to respond to the somatosensory stimulation. The coordinates, adapted by Roux et al. [17] are expressed in MNI standard space. . . . . . . . . . . . . . . . . . . . . . . . . . introduced with the advent of neuroimaging techniques, new insights into the understanding of the functioning of the human brain have been possible. The human brain is the most interesting, fascinating as well as mysterious organ of the human body. It is where the mind and the intelligence, the soul and the personality, the thoughts and the behaviour as well as the functions that regulates and allows the movement of the entire human body are created and directed. Nonetheless, the brain does this starting from its basic component, the neuron, in a scale of microns and moreover it does it via simple electrical signals. Furthermore, the understanding of the brain becomes crucial to treat the neurodegenerative diseases due to the increased life expectancy of the population in the world. The ultimate goal of the neuroimaging community is to decode these signals and understand how the brain works and is organized while it is working, in a living human being, in other words in vivo. Nowadays, there are several imaging modalities available and they can be classified into two main categories: the structural and functional neuroimaging modalities. (i) Magnetic Resonance Imaging (MRI) is a unique tool that provides anatomical images of the brain tissues. (ii) Functional MRI (fMRI) allows to obtain information about the brain function and how the different regions of the brain interact and integrate by exploiting a natural endogenous contrast of the properties of the blood oxygenation. (iii) Diffusion MRI (dMRI) instead allows the study of the white matter (WM) fiber bundles, hence the tracts connecting one brain region to another. (iv) Near-infrared spectroscopy (NIRS) gives information about regional oxygenation and blood flux via near-infrared region of the electromagnetic spectrum. (v) Magneto-and Electro-encephalography (MEG/EEG) perform measurements of the brain electro-magnetic activity. (vi) Computed Tomography (CT) allows 3-D scannings of the brain using x-rays. (vii) Position Emission Tomography (PET) gives information about the metabolic activity of the brain.Neuroimaging modalities give indirect and degraded measures of the brain activity or structure, meaning that the acquired signals are corrupted by noise, and have to be interpreted using sophisticated analysis methods which allow to extract meaningful information.

  Papers • I. Costantini, S. Deslauriers-Gauthier, and R. Deriche, A Paradigm Free Regularization Approach to Recover Brain Activation from Functional MRI Data. Manuscript in preparation. • S. Deslauriers-Gauthier, I. Costantini and R. Deriche, Non-invasive inference of information flow using diffusion MRI, functional MRI, and MEG, Journal of Neural Engineering. Manuscript submitted for publication. Participation in Conferences • I. Costantini, S. Deslauriers-Gauthier, and R. Deriche, A Paradigm Free Regularization Approach to Recover Brain Activations: Validation on Task fMRI, International Society for Magnetic Resonance in Medicine, ISMRM 2020, Sydney, Australia. (Oral Presentation) • I. Costantini, S. Deslauriers-Gauthier, and R. Deriche, Deconvolution of fMRI Data using a Paradigm Free Iterative Approach based on Partial Differential Equations, Organization for Human Brain Mapping Annual Meeting, OHBM 2019, Rome, Italy. • I. Costantini, S. Deslauriers-Gauthier, and R. Deriche, Novel 4-D Algorithm for Functional MRI Image Regularization using Partial Differential Equations, International Society for Magnetic Resonance in Medicine, ISMRM 2019, Montreal, Canada. (Power Pitch) • I. Costantini, P. Filipiak, K. Maksymenko, S. Deslauriers-Gauthier, and R. Deriche, Temporal Regularized Deconvolution of fMRI Data using a LASSO Model and the LARS Algorithm, C@UCA, June 2018, Frejus, France. • I. Costantini, P. Filipiak, K. Maksymenko, S. Deslauriers-Gauthier, and R. Deriche, fMRI Deconvolution via Temporal Regularization using a LASSO model and the LARS algorithm, 40th International Conference of the IEEE Engineering in Medicine and Collaboration with other authors • M. Frigo, I. Costantini, R. Deriche, and S. Deslauriers-Gauthier, (2018, September). Resolving the crossing/kissing fiber ambiguity using Functionally Informed COMMIT. In International Conference on Medical Image Computing and Computer-Assisted Intervention (pp. 335-343). Springer, Cham. (Appendix C). • M. Frigo, GG. Diez, I. Costantini, A. Daducci, D. Wassermann, R. Deriche, and S. Deslauriers-Gauthier, Reducing false positive connection in tract function filtering, Organization for Human Brain Mapping, OHMB 2018, Singapore.

FIGURE 1. 1 :

 1 FIGURE 1.1: Neurons structure and process of synaptic transmission. The figure shows how the neurons are composed and connected to each other. The cell body contains the nucleus and cytoplasm. The axon extends from the cell body and gives rise to small branches before ending at nerve terminals. Dendrites extend from the neuron cell body and carry the information between neurons. Synapses are the points where the "communications" between neurons take place. Figure taken from [1]

Figure 1. 2 :

 2 from the front (frontal or anterior view), from a side (lateral view), from the back (occipital or posterior view) and from the middle (medial view). The brain lobes are exhaustively illustrated in Figure 1.3. On the dorsolateral surface, from point 1 to point 2 in Figure 1.3.A, the line of separation between the frontal lobe and the parietal lobe is the central sulcus. The lateral sulcus, from point 2 to point 3, is the division line between the frontal lobe and the temporal lobe; it also partially separates the temporal lobe from the parietal lobe, from point 2 to point 4. An arbitrary line that goes from the dorsal tip of the parieto-occipital sulcus and the preoccipital notch, from point 5 to point 6, separates the occipital lobe from the parietal and temporal lobes. Another arbitrary line, which runs from the anterior edge of the occipital lobe to the posterior tip of the lateral sulcus, from point 7 to point 4, is dividing the ventral surface and the posterior parietal lobe. On the medial surface of the brain, as illustrated in Figure 1.3.B, from point 8 to point 9, the anterior, posterior and dorsal borders of the limbic lobe are demarcated by the cingulate, olfactory, parieto-occipital and subparietal sulci. The ventral borders of the limbic lobe are instead separated by the rhinal and the collateral sulci, from point 8 to point 10. The frontal lobe is segregated from the parietal one via the line extending from point 1 to point 11,

FIGURE 1. 2 :FIGURE 1 . 3 :

 213 FIGURE 1.2: Conventional terminology which indicates the different orientations representing the different brain views: in surface (A), section (B) and connectional anatomy (C). Surface neuroanatomy is the one describing the appearance of grooves (sulci) and folds (convolutions or gyri) of the brain. Sectional neuroanatomy illustrates the cortical and subcortical structures and their relationship, commonly described along the axial, coronal and sagittal plane. The connectional neuroanatomy desribes the connecting fiber tracts' origin, paths and terminations. Figure taken from Catani and Thiebaut de Schotten [2].

  date, Magnetic Resonance Imaging is the most diffused technique based on Nuclear Magnetic Resonance (NMR) for medical imaging. The NMR was discovered by Isidor Isaac Rabi, who received the Nobel Prize in Physics in 1944. In the late 1940s and early 1950s, FelixBloch at Stanford[START_REF] Bloch | Nuclear induction[END_REF][START_REF] Purcell | Resonance absorption by nuclear magnetic moments in a solid[END_REF] and Edward Purcell at the Massachusetts Institute of Technology, independently developed the NMR spectroscopy, an analytical chemistry technique that examines the structure of molecules. The technique has rapidly evolved since then. As for MRI, the first acquisition realized using this technique was done in 1976 by Damadian et al.[START_REF] Damadian | Nmr in cancer: Xvi. fonar image of the uve human body[END_REF], but the technique was approved for clinical use almost ten years later[START_REF] Gangarosa | Operational safety issues in mri[END_REF]. In 2003 the Nobel Prize in Physiology or Medicine was awarded jointly to the two scientists Paul Lauterbur and Peter Mansfield, who played key roles in the development of MRI.

FIGURE 1. 4 :

 4 FIGURE 1.4: Longitudinal and transverse relaxation. T1 corresponds to the time in which the longitudinal component regain a value which is the 37% lower of its equilibrium value; T2 corresponds to the time in which the transverse component is reduced to the 37% of it. Figure taken from Pizzolato [3].

FIGURE 1. 5 :

 5 FIGURE 1.5: Examples of MRI brain scans in the axial plane: T1-weighted (left) and T2-weighted images (right).

G x cause gradient

  Figure 1.6.

FIGURE 1. 6 :

 6 FIGURE 1.6: Echo Planar Imaging (EPI) Sequence -(A) EPI sequence with phase-encoding along y direction using blipped gradient pulses. (B) Trajectory of the EPI sequence in the k-space. The different colors of the trajectories refer to the gradients G x and G y .

FIGURE 1. 7 :

 7 FIGURE 1.7: Schematic illustration of the generation of the blood-oxygenlevel-dependent (BOLD) signal. Figure adapted from Iannetti and Wise [4].

FIGURE 1. 8 :

 8 FIGURE 1.8: Hemodynamic response function (HRF) to a stimulus with a short duration, illustrated with the red bar. The peak related to the bloodoxygen-level-dependent (BOLD) effect starts to increase after approximately 3 s from the stimulus starting point. Figure taken from Amaro and Burker [5].

FIGURE 2

 2 FIGURE 2.1: The approach of the general linear model (GLM). Each time series of fMRI data (Y) is a linear combination of regressors, composed both of task-related regressors (green) and nuisance regressors (light green) weighted by β, and the error . The goal of the GLM is to minimize the error thus estimating the parameters weights β.

FIGURE 2 . 3 :

 23 FIGURE 2.3: 13 innovative-driven co-activation patterns (iCAPs) ordered with respect to their occurrence during a resting-state fMRI acquisition. The iCAP 1 contains auditory regions; the iCAP 2 includes regions of the frontoparietal attention network; the iCAP 3 as well as the iCAP 4 cover the primary and secondary visual areas; the iCAPs 5 reveals the precuneus, the posterior cingulate cortex and the thalamus; the iCAP 6 represents the visuospatial/dorsal attention network; the iCAP 7 covers the motor network and the medial fromtal gyrus; the iCAP 8 corresponds to posterior part of the DMN; the iCAP 9 includes the anterior executive network; the iCAP 10 shows again a posterior part of the DMN; the iCAP 11 reveals the anterior salience network; the iCAP 12 is composed by the combination of different regions which are located in the limbic and the subcortical area, superior and middle temporal and occipital gyrus; the iCAP 13 contains the frontal gyrus, the anterior cingulate cortex and caudate. Figure adapted from Karahano glu and Van De Ville [7].

FIGURE 2 . 4 :

 24 FIGURE 2.4: Seed-based method -To examine the time series and extract meaningful information from fMRI data, i.e. the functional connectivity level between a seed voxel and another brain voxel i, the correlation between the two time series showed on the bottom right of the figure are computed. Highly correlated time courses reflect a high level of functional connectivity.To map all the functional connections of the seed and obtain the map showed at the bottom left of the figure, a voxel-wise correlation between the seed voxel's time series and all the other brain voxels is calculated. The outcome is a map that shows which areas show a high level of functional connectivity with the selected seed[START_REF] Van Den | Exploring the brain network: a review on resting-state fmri functional connectivity[END_REF]. Figure taken from Hu and Zeng [9].

FIGURE 2 . 5 :

 25 FIGURE 2.5: Matrix representation of the spatial indepenedent component analysis (spatial-ICA). In spatial-ICA, the algorithm try to estimate spatially independent components with related time series.

FIGURE 2 . 6 :

 26 FIGURE 2.6: Different probabilistic independent component analysis (probabilisti-ICA)-estimated resting state outputs. Estimated from a group of 10 subjects, the eight spatial maps, coregistered and superimposed to the Montreal Neurological Institute (MNI) template, shows the sagittal, coronal and axial views of different components associated with low-frequency resting-state patterns. R: right; L: left. Figure adapted from Beckmann and colleagues [10].

FIGURE 3 . 1 :

 31 FIGURE 3.1: Representation of the 4-D functional MRI (fMRI) data. A fMRI image is composed by a set of 3-D volumes recorded over time, thus leading to a 4-D structure. TR = repetition time. Figure adapted from [11].

FIGURE 3 . 2 :FIGURE 3 . 3 :

 3233 FIGURE 3.2: From top left, following the black lines connecting the different plots, s(t) is the spike train (innovation signal) which induces the activations.I α is the exponential accumulation function that leads from s(t) to the activity-inducing signal u(t), which is piece-wise constant. H represents the hemodynamic response function, x(t) is the activity-related signal, a is the additive noise and y(t) is the simulated acquired fMRI time course.

FIGURE 3 . 4 :

 34 FIGURE 3.4: L-curve. The graph represents for one time series of one voxel the plot of the two quantities, one related to the norm of the solution (abscissa) and the other related to the residuals (ordinate), i.e. (||ŝ|| 1 , ||y -Aŝ||22 ), paramtrized by the regularization parameter λ. Each black dot in the curve is related to each λ outputted by the Least Angle Regression (LARS) algorithm[START_REF] Efron | Least angle regression[END_REF]. The red dot represents the optimal solution as the nearest to (0, 0) ∈ R 2 .

FIGURE 3 . 5 :

 35 FIGURE 3.5: fMRI signal simulation. u(t) is the activity-inducing signal, a block of constant activation represented with a piece-wise constant signal.

2 Var

 2 [START_REF] Pizzolato | Computational diffusion & perfusion MRI in brain imaging[END_REF] we also report the peak-SNR (pSNR) for each experiment averaged across the 100 repetitions of the considered experiment. The pSNR m and pSNR a for each repetition were computed as follows pSNR m = 10 log 10 max(u)

FIGURE 3 . 6 :

 36 FIGURE 3.6: Each plot shows the root mean square errors (RMSE) and the roots of standard deviations computed between the ground truth activation and the estimated one using the Least Angle Regression (LARS) combined with the L-curve (in red) and with the leave-one-out cross-validation (LOO-CV, in light blue). Each plot refer to a different experiment, as reported in Table3.3. For each experiment, the RMSE were averaged across 100 repetitions of the simulated time series. The errors (in y-axes) are plotted with respect to the α value (in x-axes) used for the α-filter. In the 5 experiments, we added block-type noise ( b ) for false activations, and both m and a to the data, increasing it from experiment 1 to 5. The plots show that the LOO-CVbases approach is more sensitive to noise while increasing it and results are more disperse if compared with the L-curve-based ones. In all experiments the results obtained with L-curve show errors that are significantly lower than those given by the LOO-CV. With respect to the parameter α, results obtained via the LOO-CV are improved by increasing it, instead the ones obtained with L-curve are more stable with respect to it.

FIGURE 3 . 7 :

 37 FIGURE 3.7: Each plot in the figure show the root mean square errors (RMSE) and the roots of standard deviations computed between the ground truth activation and the estimated one using the Least Angle Regression combined with the L-curve (in red) and with the leave-one-out cross-validation (LOO-CV, in light blue). Each plot refer to a different experiment, as reported in Table3.3. For each experiment, the RMSE were averaged across 100 repetitions of the simulated time series. The RMSE (in y-axes) are plotted with respect to the α value (in x-axes) used for the α-filter. In experiments from 6 to 10, we increasingly added a to the synthetic times series. The plots show that the LOO-CV-based approach is less sensitive to noise if compared to experiments 1 to 5 in Figure3.6, but still more sensitive to noise while if compared with the L-curve-based ones. In all experiments the results obtained with L-curve show errors that are significantly lower than the ones obtained with the LOO-CV for α values lower than 0.5. For α > 0.75 results are similar.

FIGURE 3 . 8 :

 38 FIGURE 3.8: The graphs in the plots show examples of the recovered activation using the mixed Least Angle Regression (LARS)-L-curve approach (in red) and the mixed LARS and leave-one-out cross-validation (LOO-CV) (in light blue) superimposed on the ground truth activation in black and the noisy synthetic fMRI time series (in green). Each row shows a different experiment, meaning different noises applied to the ground truth activation and each columns refer to a specific α used in the α-filter. In experiments 1, 3 and 5 the added noise was coming from different sources: block-type noise ( b ) to simulate false-activations, model noise m and additive noise a . The standard deviations of m and a increases from experiment 1 to 5. The plots show that the recovered activations using the L-curve-based approach are much closer to the ground truth compared with LOO-CV. The L-curvebased approach shows results that are closer to the ground truth in terms of amplitude for bigger alphas. For α = 0.75, 3, results are similar, while in contrast with those given by α = 0.3.

FIGURE 3 .FIGURE 3 . 11 :

 3311 FIGURE 3.10: 3-D activation map obtained from an auditory task superimposed on the standard Montreal Neurological Institute (MNI) brain. The three images, from left to right, represent the axial, the sagittal and the coronal view. Voxels' intensity is ranged between 0 and 3.

•

  I. Costantini, P. Filipiak, K. Maksymenko, S. Deslauriers-Gauthier, and R. Deriche, Temporal Regularized Deconvolution of fMRI Data using a LASSO Model and the LARS Algorithm, C@UCA, June 2018, Frejus, France. • I. Costantini, P. Filipiak, K. Maksymenko, S. Deslauriers-Gauthier, and R. Deriche, fMRI Deconvolution via Temporal Regularization using a LASSO model and the LARS algorithm, 40th International Conference of the IEEE Engineering in Medicine and Biology Society, EMBC 2018, Honolulu, Hawaii. • I. Costantini, P. Filipiak, K. Maksymenko, R. Deriche, and S. Deslauriers-Gauthier, Deconvolution of fMRI BOLD signal in time-domain using an exponential operator and Lasso optimization, Workshop Computational Brain Connectivity Mapping, November 2017, Juan Les Pins, France.A natural step forward, subsequent to the approach proposed in Chapter 3, was to introduce into the optimization problem, already presenting a temporal regularization, also the spatial information. fMRI data has indeed an excellent spatial resolution if compared to other brain imaging techniques. For this reason here we propose an innovative approach to regularize the entire brain fMRI image, that takes into account both the spatial and temporal dimensions contained in the 4-D fMRI data structure.

generalized to a 4 -

 4 D problem. Let us define a scalar-valued image as a function I : Ω ⊂ R 3 , where Ω is the domain of the image. Let us now define a structure tensor, also called diffusion tensor, D as a 3 × 3 symmetric and positive-definite matrix. D = ∇I∇I T =

FIGURE 4 . 1 :

 41 FIGURE 4.1: Ellipsoidal representation of the 3-D diffusion tensor.

. 11 )

 11 where div( D0 ∇I 0 ) is the normalization term and D is the regularization tensor, distinct from the structure tensor D. In order to elucidate the regularization term, let us start by the definition of the operatorD = ∇I∇I T ∇I 2 * G (4.12)that is the 4-D structure tensor of I smoothed by the gaussian kernel G with standard deviation σ G , via the convolution operator * . The matrix D being the diffusion tensor of the image I, its eigendecomposition gives a set of eigenvalues and eigenvectors such that, if the gradient in one direction is large, the eigenvalue associated to that direction is large, whereas the eigenvalues associated to the other three directions are relatively small. Since we are processing fMRI images with the aim of saving activations and contours that occur concomitant to a large gradient in a certain direction, we aim at reversing the diffusion process, therefore at reversing the effect of D into D to enhance and at the same time simplify coherent structures of the fMRI image. We computed the operator D starting from the operator D in Eq. (4.12). After computing D, we defined the directions of the image variations by an eigendecomposition of D such that D = QΛQ T (4.13)

FIGURE 4 . 2 :

 42 FIGURE 4.2: Ground truth for the functional MRI (fMRI) simulated data: (a) activation map. (b) Simulated activation u(t), with a repetition time (TR) of 1 s.

Figure 4 .Figure 4 . 3

 443 Figure 4.3 shows examples of recovered spatial maps (Figure 4.3.a) and time series (Figure 4.3.b) using the PF-fMRI ( ûPFf MRI ) and the TA ( ûTA ). Both approaches do not require any prior knowledge of the paradigm timing or duration. The regularized spatial maps inFigure 4.3.a in the axial plane show how the denoised fMRI image recovered using the PF-fMRI ( ûPFf MRI ) is closer to the ground truth (u) in terms of signal amplitude with respect

Figure 4 .

 4 Figure 4.5 shows how the roots of MSEs ± STDs computed between the simulated activation u(t) and the recovered activity-inducing signals û(t) change for different pSNRs. We show lower errors with lower standard deviations than the ones obtained using TA.

Figure 4 .

 4 Figure 4.6 shows that the activation recovered with the PF-fMRI is more correlated with the ground truth (r ≈ 1), for different pSNRs. Whereas the results obtained with TA are more sensitive to noise and show better performances for less noisy data.

Figure 4 . 6 ,

 46 Figure 4.6, and the two additional ones that we show in Figure 4.7. Figure 4.9 reports the Pearson correlation coefficients (r) between the ground truth u(t) and the recovered signals ûPFf MRI (t) averaged across the voxels belonging to the GM.

FIGURE 4 . 3 :

 43 FIGURE 4.3: (a)From left to right: spatial maps of the simulated functional MRI (fMRI) image y, ground truth activation u, recovered activation using the Total Activation (TA) approach ( ûTA ) and our approach ( ûPFf MRI ). Each row corresponds to a a different peak-SNR (pSNR): 6.54 dB, 5.99 dB, 5.9 dB, 3.93 dB from the top to the bottom. (b) Reconstructed time series û(t) obtained with our approach ûPF-f MRI (t) (red) and the TA approach ûTA (t) (blue) superimposed on the ground truth activation u(t) (black) and fMRI signal y(t) (green). A zoom of the plot is shown in Figure4.4.

FIGURE 4 . 4 :

 44 FIGURE 4.4: The plot is a zoom of Figure4.3 representative of a trend that we found in all the four experiment showed in Figure4.3. From top left: spatial maps of a slice of the simulated functional MRI (fMRI) image y, ground truth activation u, recovered activation maps obtained using the Total Activation (TA) approach ( ûTA ) and the PF-fMRI ( ûPFf MRI ). The map obtained using TA ( ûTA ) had lower amplitude compared with the ground truth, nonetheless the scale between the values were kept. The map obtained using PF-fMRI ( ûPFf MRI ) showed amplitude comparable to the ground truth, the scale between the values were kept. The time series on the bottom shows that the recovered signal ûPF-f MRI (t) (red) was closer to the ground truth u(t) (black) compared to ûTA (t) (blue).

FIGURE 4 . 5 :

 45 FIGURE 4.5: The graph shows, for different peak-SNRs, the roots of the mean square errors (MSE) and standard deviation (STD) between u(t) and û(t) averaged among the voxels belonging to the gray matter.

FIGURE 4 . 6 :

 46 FIGURE 4.6: The graph shows, for different peak-SNRs (pSNRs), the Pearson correlation coefficient computed between u(t) and û(t) and averaged among the voxels belonging to the gray matter and their standard deviation. (µ r : mean correlation coefficient; σ r : standard deviation of the correlation coefficients.)

FIGURE 4 . 7 :

 47 FIGURE 4.7: (a) From left to right: spatial maps of the simulated functional MRI (fMRI) image y, ground truth activation u, recovered activation using our approach ( ûPFf MRI ). Each row corresponds to a different peak-SNR (pSNR): 0 dB, -6.51 dB, from the top to the bottom. (b) Reconstructed time series obtained with our approach ûPF-f MRI (t) (red) superimposed on the ground truth activation u(t) (black) and the original fMRI signal y(t) (green).

FIGURE 4 . 8 :

 48 FIGURE 4.8: The graph shows, for different peak-SNRs (pSNR), the roots of the mean square errors (MSEs) and standard deviation (STD) between the ground truth signal u(t) and the recovered signals ûPF-f MRI (t) averaged across the voxels belonging to the gray matter.

FIGURE 4 . 9 :

 49 FIGURE 4.9: The graph shows, for different peak-SNRs (pSNRs), the Pearson correlation coefficient computed between the ground truth signal u(t) and the recovered activation ûPF-f MRI (t) averaged over the voxels belonging to the gray matter and their standard deviation. (µ r : mean correlation coefficient; σ r : standard deviation)

Figure 4 . 11 .FIGURE 4 . 10 :FIGURE 4 . 11 :

 411410411 FIGURE 4.10: Qualitative comparison between the General Linear Model (GLM) and our approach (PF-fMRI). On the left column, in a blue-lightblue color-map, superimposed to the standard Montreal Neurological Institute (MNI) brain, the β-regressors map obtained using the GLM implemented in the FMRIB Software Library (FSL) tool. On the right column, in a red-yellow color-map, the whole-brain voxel-wise correlation maps obtained using the PF-fMRI superimposed to the standard MNI brain. The Pearson correlation was computed voxel-wise among the whole brain, between the reconstructed activity inducing signals û(t) and the five motor tasks simulated as piecewise constant signals with ones in the time points where the subject is executing a task and zeros elsewhere. The values r of the correlations are indicated by the color-bars. Each row corresponds to a specific motor task, from top to bottom: the tongue, the right and left hand, and the right and left foot. A: anterior; P: posterior; S:superior; I: inferior; R: right; L: left.

Figure 4 .

 4 10 and 4.12 they are referred to a single subject of the HCP database, while Figure4.11 is related to the sample of 51 subjects. This is the reason why the standard deviation in Figure4.10 is greater than those observed in Figure4.12. Moreover, the ROIs used for Figure4.10 are not the same ones used for Figure4.12, therefore results are not directly comparable.Finally, as for the resting-state data analysis, Figure4.13 shows, for iterations 5, 15 and 45 the recovered activations obtained using the PF-fMRI algorithm. Results shows how using the PF-fMRI we obtained a time courses which is related to the activity of the DMN and how these activation become clearer with increasing the number of iterations.

FIGURE 4 . 12 :

 412 FIGURE 4.12: a) Reconstructed signals û(t) obtained with our approach (PF-fMRI, red) and the total activation tool (TA, blue) superimposed on the real acquired fMRI signals (green). The ground truth in this case corresponds to the task (black), simulated as piece-wise constant signals with ones where the subject was asked to perform the task and zeros elsewhere. The plot on the top is related to the region of interest (ROI) located on the Brodmann Area 4p (rBA4p), the plot on the bottom is associated to the ROI positioned on the primary auditory cortex (rTE1.2). All the signals were averaged across the voxels belonging to the gray-matter (GM)-masked ROIs. The grey areas represent the occurrence and the duration of the tongue movements. b) Simulated tongue activation. c) Mean Pearson correlation coefficients (µ) and their associated standard deviations (σ) computed between the tongue activation and the recovered signals û(t) averaged across the voxels belonging to the GM-masked ROIs (rBA4p on the left, rTE1.2 on the right). The blue curves are related to the TA approach, while the red one to the PF-fMRI approach. TR: repetition time.

•

  I. Costantini, S. Deslauriers-Gauthier, and R. Deriche, A Paradigm Free Regularization Approach to Recover Brain Activation from Functional MRI Data. Manuscript in preparation. • I. Costantini, S. Deslauriers-Gauthier, and R. Deriche, A Paradigm Free Regularization Approach to Recover Brain Activations: Validation on Task fMRI, International Society for Magnetic Resonance in Medicine, ISMRM 2020, Sydney, Australia. (Oral Presentation). • S. Deslauriers-Gauthier, I. Costantini and R. Deriche, Non-invasive inference of information flow using diffusion MRI, functional MRI, and MEG, Journal of Neural Engineering. Manuscript submitted for publication. • I. Costantini, S. Deslauriers-Gauthier, and R. Deriche, Deconvolution of fMRI Data using a Paradigm Free Iterative Approach based on Partial Differential Equations, Organization for Human Brain Mapping Annual Meeting, OHBM 2019, Rome, Italy. • I. Costantini, S. Deslauriers-Gauthier, and R. Deriche, Novel 4-D Algorithm for Functional MRI Image Regularization using Partial Differential Equations, International Society for Magnetic Resonance in Medicine, ISMRM 2019, Montreal, Canada. (Power Pitch)

S

  . Deslauriers-Gauthier, I. Costantini and R. Deriche, Non-invasive inference of information flow using diffusion MRI, functional MRI, and MEG, Journal of Neural Engineering. Manuscript submitted for publication.This work has as objective to infer information flow in the WM of the brain and recover cortical activity using fMRI, dMRI, and MEG without a manual selection of the WM connections of interest. A description of the paper is reported in Appendix C.

2 .

 2 A paradigm-free regularization algorithm based on PDEs named PF-fMRI that applies on the 4-D fMRI image, acting simoultaneously in the 3-D space and the 1-D time dimension.

(A. 9 )

 9 Proof of Eq. (A.9). The accumulation of a discrete time signal x(n) as the convolution of a unit step sequence u(n)u(n) =    1 for n ≥ 0 0 for n < 0 with x(n). This convolution is

FIGURE A. 2 :FIGURE A. 3 :

 23 FIGURE A.2: Function g(n) defined in Eq. (A.15), for S = -1 and α = 1.

Table 3

 3

	.1 and 3.2.
	Variable	Definition
	u	Stimulus
	s	Flow inducing signal
	f in	Blood flow
	v	Normalized venous volume
	q	Normalized deoxyhemoglobin voxel content

TABLE 3 .

 3

1: Variables involved in the hemodynamic model

[START_REF] Friston | Nonlinear responses in fmri: the balloon model, volterra kernels, and other hemodynamics[END_REF][START_REF] Khalidov | Activelets: Wavelets for sparse representation of hemodynamic responses[END_REF]

.

TABLE 3 .

 3

2: Constants involved in the hemodynamic model

[START_REF] Friston | Nonlinear responses in fmri: the balloon model, volterra kernels, and other hemodynamics[END_REF][START_REF] Khalidov | Activelets: Wavelets for sparse representation of hemodynamic responses[END_REF]

.

  Experiment σ m pSNR m [dB] σ a pSNR a[dB]

	1	0.1	20	0.1	28
	2	0.3	11	0.3	19
	3	0.5	6	0.5	15
	4	0.7	3	0.7	14
	5	0.9	1	0.9	13

TABLE 3 .

 3 3: Experiments set for testing the L-curve and the LOO-CV to compare the two approaches for the selection of optimal lambda. The first column indicates the experiment number, the second column the standard deviation (σ m ) of the gaussian random noise added as model noise, the third column represents the peak signal to noise ration (pSNR m ) computed as in (3.28) in dB. The fourth column reports the σ a for the additive gaussian noise and the fifth column the pSNR a estimated as in (3.29) expressed in dB.

  [START_REF] Amaro | Study design in fmri: basic principles[END_REF]. All the voxels were corrupted by noise: the

	Experiment σ m σ a	pSNR[dB] A pSNR[dB] B
	1	0.1 0.1	5.17	7.64
	2	0.3 0.3	4.96	6.84
	3	0.6 0.6	4	5.94
	4	1	1	3.14	3.29
	5	0	1	3.98	5.12

TABLE 3 . 4 :

 34

TABLE 3 .

 3 5: Summary of rooted MSEs and STDs obtained for phantom fMRI data.

	Activation			A					B
	Experiment	1	2	3	4	5	1	2	3	4	5
	pSNR [dB]	5.17 4.96	4	3.14 3.98 7.64	6.84 5.94 3.29 5.12
	OUR (α = 0.75)	rMSE 0.11* 0.15 0.22 0.18 0.34 0.05* 0.1* 0.19 0.31 0.14* rSTD 0.13* 0.14 0.16 0.2 0.15 0.05* 0.07* 0.11 0.17 0.08*
	TA	rMSE 0.25 0.26 0.29 0.36 0.3 rSTD 0.31 0.31 0.32 0.33 0.33 0.23 0.19	0.2 0.24 0.24 0.27 0.26 0.24 0.32 0.25

TABLE 3 .

 3 6: Tongue task paradigm of HCP data. Note that the time for activation is expressed in seconds and the repetition time for the HCP data is 0.72s.

	Task	Starting time End time
	Tongue	41.264 s 101.639 s	53.264 s 113.639 s

Table 4

 4 .1. The tasks starting points were considered equal for each subject and inter-subjects differences of the order of milliseconds were neglected.

	Task	Starting time End time
	Right hand	11.009 s 131.894 s	23.009 s 143.894 s
	Left foot	26.136 s 116.766 s	38.136 s 128.766 s
	Tongue	41.264 s 101.639 s	53.264 s 113.639 s
	Right foot	56.391 s 177.142 s	68.391 s 189.142 s
	Left hand	71.518 s 162.518 s	83.518 s 174.518 s

TABLE 4 .

 4 1: Motor task paradigm of HCP data. The time is expressed in seconds and the repetition time for the HCP data is 0.72s.

Table 4

 4

	.2.

TABLE 4 .

 4 2: MNI coordinates centers of the brain areas found to respond to the somatosensory stimulation. The coordinates, adapted by Roux et al.[START_REF] Roux | Functional architecture of the somatosensory homunculus detected by electrostimulation[END_REF] are expressed in MNI standard space.

  MNI coordinates: 58, -18, -14) were located in BA 21. We merged together the six ROIs into a unique network, we masked it with the GM mask and we computed the average of the recovered time series obtained using the PF-fMRI to obtain a recovered signal for the DMN.

	and colleagues
	(2008) [65]. The medial prefrontal cortex (mPFC; MNI coordinates: 0, 54, -9) was positioned
	in BA 10. The posterior cingulate cortex (PCC; MNI coordinates: 0, -53, 26) was placed in in
	BA 29/30. The left inferior parietal lobule (lIPL; MNI coordinates: -51, -57, 27) and the right
	inferior parietal lobule (rIPL; MNI coordinates: 51, -63, 30) were placed in BA 39/40. The left
	lateral temporal cortex (lLTC; MNI coordinates: -58, -18, -14) and the right lateral temporal
	cortex (rLTC;

https://www.humanconnectome.org/storage/app/media/documentation/s1200/HCP_S1200_Release_ Appendix_I.pdf

https://miplab.epfl.ch/index.php/software/total-activation

Acknowledgements

This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (ERC Advanced Grant agreement No 694665: CoBCoM -Computational Brain Connectivity Mapping). Data were provided [in part] by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.

List of Abbreviations

The analysis of fMRI images represents a challenge for the researchers in the field because they are asked to deal with fMRI time courses that actually constitute the neural signal convolved with a HRF [START_REF] Gitelman | Modeling regional and psychophysiologic interactions in fmri: the importance of hemodynamic deconvolution[END_REF]. Nonetheless, the interactions we want to investigate in the brain are expressed at a neural level, rather than at the level of hemodynamic responses. This leads to the need of recovering the neural signal that allows to model the interplay between brain regions. This is the framework where deconvolution approaches are required. In fact, given the BOLD signals, deconvolution-based methods allow an approximation of the underlying activations given a model of the HRF.

Neverthless, fMRI measurements are corrupted by noise and therefore the addressed inverse problem, commonly present in signal and image processing, is typically ill-posed, because small perturbations in the data may give big changes in the solution. For this reason, additional assumptions are necessary in order to solve the inverse problem.

In this chapter we will describe our first main contribution. We developed a temporal regularized deconvolution of the fMRI BOLD signals that aims at the paradigm-free recovery of brain activations at the neural level, namely without the necessity of a priori informations on the experimental paradigm, as it was for the GLM approach. For completeness, we will start by describing the regularization approach, in particular the LASSO model that we used.

We will explain how we solved the inverse problem and how we validated our approach on phantom and on real task-fMRI data. We will finally compare it with the TA approach [START_REF] Karahano Glu | Total activation: fmri deconvolution through spatio-temporal regularization[END_REF][START_REF] Farouj | Regularized spatiotemporal deconvolution of fmri data using gray-matter constrained total variation[END_REF] already presented in Section 2.2.

Introduction to Inverse Problems and Regularization

Approaches

In this section we will give an overview of the tools we used to define the forward model we are facing in this thesis and the approach we implemented to subsequently solve the inverse problem.

Second Contribution

Appendices

Appendix A

α-Filter Design

In this appendix we report the design of the α-filter employed in the first contribution described in Section 3.4.2 of this dissertation.

To do this we initially defined an antisymmetric function with exponential decay

where S is a normalization factor and α is a parameter that defines the shape, i. Note that h(n) represents the derivative filter for optimal edge detection [START_REF] Deriche | Using canny's criteria to derive a recursively implemented optimal edge detector[END_REF] that can be adapted by tuning only one parameter, i.e. α, such that if α 0 the filter corresponds to the classical definition of a discrete derivative

Otherwise, for smaller α the filter can be exploited in the way that it detects edges while smoothing noise.

In order to define the transfer function of the filter we want to design, we start by defining the convolution between two functions as

Let us focus only on the causal part of the filter because the anti-causal part follows a similar procedure. If n > 0 Let us now define the property of geometric series

for r < 1.

The z-transform of Eq. (A.4) is

that, for the property in Eq. (A.5), it corresponds to

1e -α z -1 2 .

(A.8)

The transfer function in z-domain of the exponential accumulation function results as

(A.12)

Let us now define the normalization factor S such that if

we want to normalize such that

where, given h(n) in Eq. (A.1),

In order to compute the normalization factor S that satisfies Eq. (A.14), let us first define the properties of geometric series Then, for n > 0

If we apply the change of variable defined in Eq. (A.18)

Now, for the property in Eq. (A.5) we can write

For the property in Eq. (A. [START_REF] Khalidov | Activelets: Wavelets for sparse representation of hemodynamic responses[END_REF])

Let us now define

and

Given that we want to satisfy Eq. (A.14) to find the factor S that normalizes the function I α (n), let us now focus on the function g(n) of which a zoom is shown in Figure

This allows us to compute the normalization factor S using only the positive part of g(n).

Note that the value at g(0) is accounted for just once and this is justified in Eq. (A.27).

We can therefore normalize the filter, such that Eq. (A.14) is satisfied, by imposing

that, as stated in Eq. (A.25), we can rewrite as

For Eq. (A.23) and (A.24) the previous equation is equal to

For the property in Eq. (A.17) the first term on the left of Eq. (A.31)

For properties Eq. (A.5) the second term on the left of Eq. (A.31)

The solution to Eq. (A.30) is therefore

and, for the change of variable defined in Eq. (A.18), the normalization factor is

Finally, the normalized function

In conclusion, to estimate in a robust way the first derivative of a function f (n), one have to convolve the function f (n) with the following derivative filter

with S given by Eq. (A.35).

Note that for large values for α, all the coefficients of this filter become equal to zero unless

which leads to Eq. (A.2) when convolving this filter with the function f (n), i.e the filter acts as the finite difference operator usually applied for a discrete derivative. Using this normalized and optimal derivative filter, with small values of α, allows to better estimate in a robust way the first derivative of a noisy signal. The objective of this work is to infer information flow in the WM of the brain and recover cortical activity using a multimodal approach based on fMRI, dMRI, and MEG without a manual selection of the WM connections of interest. The fMRI data was exploited because of its good spatial resolution. In this work, the PF-fMRI approach allowed us to uncover brain activations without using a predefined knowledge on the experimental paradigm, thus finding and locating the active regions. The found active regions showed in Figure B.1 where then used as priors on cortical regions, to filter the structural connections. This was done such that only those connections reaching regions with a fMRI activation higher than a certain threshold were kept. Results show that the proposed method is able to identify connections associated with the a sensory-motor task, without requiring prior knowledge of the experimental paradigm. This would allow, for example, to use the same processing This work is related to tractography, an ill-posed problem that aims at reconstructing the WM fiber tracts from dMRI data. Tractography have been shown to generate a significant number of false positive connections between brain regions [START_REF] Maier-Hein | The challenge of mapping the human connectome based on diffusion tractography[END_REF]. This is due to the fact that whenever two bundles of axons are organized in a crossing or kissing configuration, in-vivo tractography is not able to distinguish their actual trajectories from diffusion data.

This implies the presence of many spurious entries in the structural connectome obtained from a tractogram. In the context of tractography optimization or tractogram filtering, what we propose in this work is to overtake the purely structural-based approach and consider both structural and functional information in the formulation of the tractography optimization problem. We do this by showing that the injection of functional priors coming from rs-fMRI gives an effective answer to the false positives issue. The functional information is exploited in the form of static functional connectivity and it plays the role of promoting the involvement of bundles that connect highly correlated cortical regions in the fitting of the dMRI signal. We encapsulate all of this within the Convex Optimization Modeling for Microstructure Informed Tractography (COMMIT) framework [START_REF] Daducci | Commit: convex optimization modeling for microstructure informed tractography[END_REF], which provides a flexible tool for defining a convex optimization problem that, for a given tractogram, simultaneously (i) promotes sparsity among the bundles, (ii) takes into account the considered functional information and (iii) selects the streamlines that are sufficient to explain the dMRI signal. We call this novel framework Functionally Informed COMMIT (FIC).