{'fraction_non_alphanumeric': 0.047070258261594, 'fraction_numerical': 0.015551235767842266, 'mean_word_length': 4.471511775132945, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The role of noise modeling in the estimation of resting-state brain effective connectivity

Causal relations among neuronal populations of the brain are studied through the so-called effective connectivity (EC) network. The latter is estimated from EEG or fMRI measurements, by inverting a generative model of the corresponding data. It is clear that the goodness of the estimated network heavily depends on the underlying modeling assumptions. In this present paper we consider the EC estimation problem using fMRI data in resting-state condition. Specifically, we investigate on how to model endogenous fluctuations driving the neuronal activity.


INTRODUCTION
A deeper understanding of the complex interactions which take place inside our brain, which currently represents a major open research challenge, could possibly open new avenues to novel terapeutic treatments as well as early diagnosis tools. It is thus not surprising that this area is attracting the interest of scientists from many different fields, as it calls for diverse expertise ranging from medical science to physics, from biology to engineering and statistics. In particular, the estimation of effective connectivity, accounting for the causal relationships among different neuronal populations, has been the subject of several research contributions in the last decade (Valdes-Sosa et al., 2011;Razi and Friston, 2016).
Effective connectivity is typically inferred from EEG or fMRI measurements, by inverting a generative model of the corresponding data. With the language of System Identification, estimating an effective connectivity model can be framed as the problem of inferring a dynamic network model from indirect measurements. In fact, such models should not only describe the causal interactions taking place at neuronal level (described by the so-called effective connectivity), but also capture the mapping between such activity and the available measurements. Specifically, when the fMRI technique is adopted, this mapping is called hemodynamic response and links the underlying synaptic dynamics to the so-called BOLD (Blood Oxygenation Level Dependent) signal, which is measured through fMRI. This signal reflects an increase of the neuronal activity on a change of the relative levels of oxyhemoglobin and deoxyhemoglobin.
The neuroscience literature has proposed different characterizations of the aforementioned generative models, ranging from linear to non-linear models, from discrete-time (DT) to continuous-time (CT) frameworks (see Valdes-This work has been partially supported by the project BIRD162411/16 and the FIRB project RBFR12M3AC funded by MIUR. Sosa et al. (2011) for a comprehensive review). In this work we consider the Dynamical Causal Model (DCM) proposed by Friston et al. (2003), which is nowadays a widespread approach in the neuroscience community . A DCM is a CT nonlinear statespace model, where the hemodynamic response is described by the so-called Balloon-Windkessel model (Buxton et al., 1998). Whereas the original version of DCM was suited for fMRI data recorded while the subject is performing a task (Friston et al., 2003), the DCM has also been recently extended to account for resting-state data (Friston et al., 2014). While in task-dependent DCM, neuronal dynamics is driven by deterministic inputs (external stimuli), in resting-state DCM it is elicited by stochastic sources, representing brain endogenous fluctuations. In the present paper, we will focus on DCM for resting-state data.
In Friston et al. (2014) and Razi et al. (2015), brain endogenous fluctuations are assumed to have a power-law spectrum. Such choice not only leads to a low frequencyconcentrated spectrum, but also accounts for the scalefree properties of neuronal fluctuations which have been observed in previous works (Stam and De Bruin, 2004;Linkenkaer-Hansen et al., 2001;Shin and Kim, 2006). Other contributions by the same authors (Friston, 2008;Friston et al., 2008) suggest the use of non-Markovian noise processes when modeling biophysical fluctuations, such as those taking place in the brain.
In our previous work (Prando et al., 2017) we tackled the estimation of resting-state effective connectivity by adopting a linearized DCM and postulating that neuronal activity is driven by white Gaussian noise; the effective connectivity matrix was inferred by imposing a sparsity inducing prior (Wipf and Nagarajan, 2010). The remarkable advantage of such a choice is that we are able to avoid the combinatorial bottleneck which typically characterizes the DCM inversion algorithms. The main contribution of this work is twofold: first, following the modeling assumptions for the endogenous fluctuations suggested in the aforementioned literature, we extend the model used in Prando et al. (2017) allowing for a CT autoregressive (AR) model of the random fluctuations driving the neuronal activity. Second, using both real and synthetic data, we compare white vs. AR(1) noise models as sources for the resting-state neuronal activity. This comparison is carried out both studying the effect of noise models in capturing functional correlations as well as evaluating the difference between effective connectivity estimated under different modeling assumptions. Our preliminary results suggest that, indeed, a low-pass AR(1) noise is more suited in our setup.
The manuscript is organized as follows. Modeling aspects are treated in Sec. 2, with a specific focus on the modeling of brain endogenous fluctuations. An extensive simulation study is reported in Sec. 3. Concluding remarks and future investigations are outlined in Sec. 4.
Notation: I n denotes the identity matrix of size n and 0 m×n the m × n zero matrix; (blk)diag(·) denotes the (block-)diagonal operator. A T indicates the tranpose of matrix A.

EFFECTIVE CONNECTIVITY MODELING
Inference of brain effective connectivity requires the definition of a generative model of the observed BOLD timeseries measured by fMRI scanners. Here we use the DCMs introduced in (Friston et al., 2003(Friston et al., , 2014; these are nonlinear state-space models which jointly describe the neuronal activity and the map from the neural activity to the BOLD signal. Denoting with x(t) = [x 1 (t), ..., x n (t)] T ∈ R n the activity of n neuronal populations at time t ∈ R, its evolution is described aṡ (1) where A ∈ R n×n is the effective connectivity matrix. The process w(t) represents random fluctuations, whose modeling is the main focus of this work. Previous contributions (Friston et al., 2014;Razi et al., 2015) have postulated that w(t) has a power-law power spectrum: (2) with ν and β being parameters to be inferred and ω denoting the angular frequency. The DCM specification is completed by the so-called Balloon-Windkessel model (Buxton et al., 1998), a 4thorder non-linear SISO system which takes as input each x i (t) and outputs the measured BOLD signal y i (t): (4) The model depends on the parameters θ b := {κ i , γ i , τ i , ρ i , ξ} that are characterized in the literature by empirical priors and have to be estimated from data (Friston et al., 2003). The states {r i , f i , v i , q i } are biophysical quantities that are affected by the neuronal activity: r i denotes the vasodilatatory signal, f i is the blood inflow, v i and q i are respectively the blood volume and the deoxyhemoglobin content. The output equation (4) includes the measurement noise e i , while the signal component depends on the resting blood volume fraction V 0 (typically V 0 = 0.4) and on the constants k 1 , k 2 and k 3 . We refer the reader to Stephan et al. (2007) for more details.
In our previous work (Prando et al., 2017) we proposed an alternative formulation of the DCM in Eqs. (1)-(4). First, since fMRI data are acquired with a sampling time T R = 2 s., we adopt a sampled version of Eq.
(1) and model w(t) as a CT white noise with intensity σI n . Then, defining x(k) := x(kT R ), we discretize Eq. (1) as (6) Second, we have proposed a statistical linearization of the Balloon-Windkessel model, exploiting the empirical priors for the hemodynamic parameters θ b . Specifically, we have reformulated model (3) with a prior on h i ∈ R s derived by the prior on θ b . The length s of the impulse response h is chosen large enough to capture all the relevant dynamics (see Sec. III in Prando et al. (2017)). Using the linearisation (7) and defining Q w :=Var{w d (k)} = blkdiag(Q w , 0 n(s−1)×n(s−1) ), the resting state DCM becomes a linear stochastic statespace model; namely: Exploiting the linearity of model (11) as well as a suitable sparsity inducing prior on A, in Prando et al. (2017) we have introduced an EM-type iterative procedure which is much less computationally demanding than variational methods typically employed to invert classical DCMs. For reasons of space we refer the interested reader to Prando et al. (2017) for details on the estimation algorithm.
Choosing the most suitable statistical description of the endogenous fluctuations w(t) in Eqs. (1) and (5)  and Razi et al. (2015) postulate that the power spectrum of w(t) is concentrated at low frequencies. Following this line of work, in this paper we reformulate system (11) assuming w(t) to be a 1st-order CT autoregressive process.

Vector Autoregressive (VAR) Process Noise
According to the discussion in the previous section, we could either model w(t) in (1) as a CT AR process, or we could directly consider the discrete time counterpart Eq. (5) and model w d (k) as DT coloured noise. Following the first option we model w(t) using the CT model is a CT white noise with intensity δI n and Λ ∈ R n×n is an Hurwitz-stable matrix. Alternatively, w d (k) can be described by the discrete-time model These two options lead in general to different models. However, in the special case Λ = λI n , the models in (12)-(13) are equivalent under a proper choice of ∆ andQ v . Indeed, sampling Eq. (12), we get dsdτ. Now, assuming Λ = λI n and recalling Eq. (6), we can rewrite Eq. (15) as the DT AR model e Aτ e 2λs e A T τ dsdτ (17) which directly depends on the dynamics of model (1). Hence, the two DT AR models in Eqs. (13) and (16) coincide only if ∆ = e λT R andQ v is set equal to Q v .
Given previous neuroscience studies (Stam and De Bruin, 2004) and the fact that w(t) directly represent the brain endogenous fluctuations, in the following we shall work with the AR(1) CT noise model (12) where Λ = diag(λ 1 , ..., λ n ), λ i ∈ R − . Under this hypothesis Eq. (1) can be rewritten as ẇ(t) (18) where the last equation defines M ∈ R 2n×2n and η(t) ∈ R 2n . The CT model (18) can be discretised obtaining: where, again, the notation x(k) := x(kT R ) has been used (analogously for w(k) and η d (k)) and η d (k) is a stationary white noise with covariance Q η = T R 0 e M τ Σe M T τ dτ, Σ = blkdiag(0 n×n , δI n ) (20) Using the statistical linearization of the Balloon-Windkessel model as above, the linear stochastic state space model The process noise η d (k) := [η T d (k) 0 1×n(s−1) ] T is white with covariance Q η := blkdiag(Q η , 0 n(s−1)×n(s−1) ). It should be noticed that the order of model (21) has increased only by n, w.r.t. that of (11). The algorithm in Prando et al. (2017) can also be applied to this scenario with obvious modifications. We thus refer the reader to Prando et al. (2017) for details. Remark 1. We should stress that modeling w(t) as an AR(1) process coincides with a specific instance of the hierarchical model proposed by Friston (2008) (see Eqs. (2)-(3)), which arises when the dynamics is expressed only in terms of the first two generalized coordinates of motion.

EXPERIMENTS
This section aims at validating (or invalidating) the assumption that the endogenous fluctuations should be modeled as an autoregressive process. To this purpose we shall perform experiments on both synthetic as well as real data. Concerning the synthetic experiment, we shall simulate fMRI signals exciting the simple brain model (1)-(4) with endogenous brain fluctuations w(t) which are either white or AR(1) noise. Both for simulated and real data we estimate the effective connectivity (EC), using either model (11) and (21) (assuming Λ diagonal). We shall then compare the results with the goal of understanding which modelling assumption is most suited to describe the real fMRI data.

Synthetic Data.
Synhetic neural responses are simulated using the discretizations (5) and (19) with T R = 0.05s, with the "true" EC matrix given by The BOLD signal is then obtained giving these simulated neural responses as inputs to the Balloon-Windkessel model (3)-(4). The output is then down-sampled to T R = 2s (the sampling rate commonly used in fMRI scanners), obtaining estimation (D i ) and test (D te i ) datasets , i = 1, .., 50, each containing a BOLD time-series of length N = 300. For each dataset the Balloon model parameters θ b are drawn from their empirical priors reported in Friston et al. (2003) and the measurement noise variance (see Eq. (4)) is chosen so as to guarantee SNR=10. When using AR input noise (see (19)) the coefficients Λ are drawn for each Monte-Carlo run from a uniform distribution U(−1, 0).

Real Data.
Real BOLD time-series from 333 brain regions (ROIs) are measured with T R = 2s on two subjects at rest on 10 minutes scanning sessions, corresponding to BOLD timeseries of length N = 300. The 333 ROIs are then clustered in so-called functional networks. Here we focus on the Visual Network (VIS), responsible for the human vision and on the Default Mode Network (DMN), which collects the most active areas during resting-state condition. In Sec. 3.5 we will use the BOLD time-series in VIS and DMN (denoted with D VIS and D DMN respectively) to infer the EC restricted to these networks. Finally, the 333 BOLD time series are reduced to 53 components by averaging on functionally homogenous areas, obtaining the dataset D wb which is used to infer the whole-brain EC.

Modeling Assumptions
We estimate the EC matrix A under three different modeling assumptions (denoted below with I := {W, AR, VAR}); in brackets the corresponding symbol used for the estimated EC:

Evaluation Metrics
The performance are evaluated according to two metrics: one directly considers the estimated EC, while the other involves an indirect measure based on the so-called functional connectivity (FC). As we shall see below, FC is nothing but the correlation between BOLD signals measured in different regions, and thus is a measure of statistical dependencies among brain regions. The interested reader is referred to Razi and Friston (2016) for further details on the distinction and relation between EC and FC.

Effective Connectivity Metrics.
In the synthetic case, where ground-truth EC is available, we can evaluate how close the estimated EC is to the "true" one using: • The Root-Mean-Squared-Error (RMSE) RMSE(Â I ) := A −Â I F / n(n − 1), I ∈ I (24) where the notation A denotes the matrix A with its diagonal set to 0 1 . • The number of errors on the sparsity pattern of the EC matrix, that is On real data, no ground-truth is available and thus the two metrics above cannot be computed. Instead, we can quantify the dissimilarity arising from the different modeling assumptions listed in Sec. 3.3 by means of the Pearson's correlation coefficient among the estimated ECs. For two generic vectors u ∈ R m and v ∈ R m , this is defined as whereâ I denotes the vectorization ofÂ I . For synthetic data the same dissimilarity measure can also be computed w.r.t. the ground-truth, i.e. ρ EC (I) := ρ(a,â I ) between the true A and the estimatedÂ I .
Functional Connectivity Metrics. FC between two brain regions is the Pearson's correlation coefficient between the corresponding BOLD time-series. Hence, the FC matrix computed from a dataset D containing n time-series Y i = [y i (1), ..., y i (N )] T , i = 1, ..., n, is given by with ρ(·, ·) as defined in Eq. (26). We also consider the FC matrix F I obtained from the output correlations of the estimated models, where the EC has been estimated according to modeling assumption I ∈ I, namely where Σ y = C Σ x C T +ˆ I n and Σ x is the solution of the Lyapunov equation Σ To assess how well the estimated EC (and thus the model), under the modeling assumption I ∈ I, is able to capture the empirical FC, we evaluate the following Pearson's correlation coefficient wheref I and f (D) respectively denote the vectorization of the upper diagonal parts of F I and of the FC matrix computed with the estimation data D.
When dealing with the synthetic experiment we can also use the test dataset D te to compute ρ FC (D te , I) := ρ(f (D te ),f I ).

Results
Synthetic Data. Figs. 1 and 2 show the results obtained in the synthetic datasets in terms of metrics (24) and (25), respectively. The three modeling assumptions I = {W, AR, VAR} give rise to similar performance, independently of how the data have been generated. The unique exception is observed on the data generated from VAR process noise when a white noise is postulated for the brain endogenous fluctuations w(t): the RMSE is significantly larger than that obtained assuming a colored process noise (that is, under assumptions AR and VAR). However, such poor performance does not negatively impact the metrics based on FC, that is ρ FC (D, I) and ρ FC (D te , I), as can be noticed in Figs. 3-4. The similarity of the performance produced by the three modeling assumptions is particularly evident in Fig. 4, which reports the correlation between the output FC of the estimated model and the empirical FC computed on test dataset. These outcomes show that all the estimated models are equally able to reproduce both estimation and test data, even if the EC matrix is not correctly inferred, thus suggesting a possible identifiability problem. Furthermore, these results seem to discredit the reliability of metrics ρ FC (D, I) in assessing the goodness of the estimated EC. However, such performance index is commonly used for this purpose when a ground truth is not available, that is when real data are treated. If we consider the correlation among the estimated ECs the situation is quite different: when data are generated by a white process noise (Fig. 5(a)), ρ EC (I, J) takes values which are very close to 1, regardless of the modeling assumption. On the other hand, when a VAR process is responsible for the excitation of model (21) (Fig. 5(b)), the EC estimated under hypothesis W is much less similar to those obtained from hypothesis AR and VAR. In this  case ρ EC (AR, VAR) is still large, but smaller than what is reported in Fig. 5(a). The correlations ρ EC (I) between the true EC matrix A and the estimated ones (first row of plots in Fig. 5) reflect what we observe in Figs. 1 and 2: while ρ EC (I), I ∈ I, are very similar in Fig. 5(a), in Fig. 5(b) the ECs estimated postulating a white process noise appear quite dissimilar from the true one, w.r.t. to those returned by AR and VAR modeling assumptions.
Real Data. Considering real data, we have to rely only on the metrics ρ FC (D, I) and ρ EC (I, J) since a ground truth is not available. Table 1 reports the correlations between the FCs (F (D)) computed from the BOLD time-series and the model FCs ( F I ) produced by the systems estimated under the modeling assumptions {W, AR, VAR}. Not surprisingly, in all the three datasets, these hypothesis give rise to equal ρ F C (D, I) values, in agreement with the results obtained on synthetic data. However, if we inspect the correlation ρ EC (I, J) between the estimated ECs (Fig. 6), the situation more closely resembles what reported in Fig. 5

CONCLUSIONS
In this work we have tackled the estimation of brain effective connectivity from both synthetic and real restingstate fMRI data. Following our previous work (Prando et al., 2017), we describe the dynamics of the neuronal activity by means of a linear state-space model, which also includes a linearized hemodynamic response, in charge of mapping the neuronal activity to the so-called BOLD signal, measured by fMRI. At rest, brain activity is elicited by endogenous fluctuations, which play the role of process noise in the aforementioned state-space system. While in Prando et al. (2017), this was considered white, here we reformulate our generative model in order to account for the presence of colored process noise. Specifically, we assume it to be a 1st-order autoregressive process. Such hypothesis is motivated by previous works in the neuroscience community, where brain fluctuations are characterized as lowfrequency signals (Friston et al., 2014;Linkenkaer-Hansen et al., 2001). Not surprisingly, using an AR process noise in the generative model seems a more robust choice, independently of the noise type which has actually generated the data. Furthermore, despite its widespread use in neuroscience, we argue that functional connectivity does not allow to discriminate between the modeling assumptions of white or AR process noise. This is to be expected since FC is just a zero lag output correlation, which is largely insufficient to identify the model parameters. Finally, evaluating the correlation among the ECs estimated under different modeling assumptions, we find a partial agreement between the results on the real data and those obtained on synthetic data generated by an AR process noise. However, we believe that deeper investigations are required to validate such statement. To this purpose we are conducting bootstrap simulations to compare synthetic and real experiments on solid statistical grounds.