Estimating effective connectivity in linear brain network models

TL;DR

This paper introduces a novel EM-based method for estimating effective brain connectivity from resting state fMRI data, accounting for the nonlinear hemodynamic response and promoting sparse network structures.

Contribution

It presents a new linear population model combined with an iterative EM procedure to accurately infer neuronal interactions from fMRI data, improving upon existing methods.

Findings

Effective estimation of brain connectivity demonstrated on rs-fMRI data.

Method outperforms state-of-the-art routines like SPM12.

Sparsity in network connections is effectively enforced.

Abstract

Contemporary neuroscience has embraced network science to study the complex and self-organized structure of the human brain; one of the main outstanding issues is that of inferring from measure data, chiefly functional Magnetic Resonance Imaging (fMRI), the so-called effective connectivity in brain networks, that is the existing interactions among neuronal populations. This inverse problem is complicated by the fact that the BOLD (Blood Oxygenation Level Dependent) signal measured by fMRI represent a dynamic and nonlinear transformation (the hemodynamic response) of neuronal activity. In this paper, we consider resting state (rs) fMRI data; building upon a linear population model of the BOLD signal and a stochastic linear DCM model, the model parameters are estimated through an EM-type iterative procedure, which alternately estimates the neuronal activity by means of the Rauch-Tung-Striebel (RTS) smoother, updates the connections among neuronal states and refines the parameters of the hemodynamic model; sparsity in the interconnection structure is favoured using an iteratively reweighting scheme. Experimental results using rs-fMRI data are shown demonstrating the effectiveness of our approach and comparison with state of the art routines (SPM12 toolbox) is provided.

Tables

Table 1. TABLE I: Performance of Algorithm 1 and classical reweighted ℓ 1 subscript ℓ 1 \ell_{1} method over 50 Monte-Carlo runs (using measured neuronal activity x ( k ) , k = 1 , . . , 600 x(k),\ k=1,..,600 . Median and std (between brackets) of metrics ( 53 ) and ( 54 ) are reported.

	$ERR\left(\hat{A}\right)$	$RMSE\left(\hat{A}\right)$
Algorithm 1	3 (2.9)	0.05 (0.11)
Linear reweighted ${\ell }_1$ algorithm	10 (1.4)	0.18 (0.006)

Table 2. TABLE II: Performance of the SPM12 routines spm_dcm_fmri_csd— ( sDCM ) and spm_dcm_estimate— ( DEM ) over 20 Monte-Carlo runs, when they are given the 21 sparsity patterns of Fig. 2 . Median and std (within brackets) of index ( 53 ) are reported, together with the no. of times ( # ​ C ​ h ​ o ​ s ​ e ​ n # 𝐶 ℎ 𝑜 𝑠 𝑒 𝑛 \#Chosen ) each sparsity pattern is chosen according to the free-energy metric. Bold font highlights the most frequently selected models.

		spectral DCM		DEM
	$ERR\left(\hat{A}\right)$	$\#Chosen$	$RMSE\left(\hat{A}\right)$	$\#Chosen$	$RMSE\left(\hat{A}\right)$
(a)	0	2	0.12 (0.05)	0	0.21 (0.004)
(b)	4	0	0.25 (0.02)	0	0.21 (0.005)
(c)	6	0	0.23 (0.03)	0	0.21 (0.003)
(d)	8	0	0.18 (0.02)	0	0.21 (0.004)
(e)	6	1	0.14 (0.05)	0	0.21 (0.005)
(f)	8	0	0.22 (0.04)	0	0.21 (0.005)
(g)	10	0	0.28 (0.09)	0	0.22 (0.007)
(h)	5	0	0.22 (0.02)	0	0.20 (0.003)
(i)	10	0	0.24 (0.07)	0	0.21 (0.004)
(j)	10	0	0.17 (0.004)	0	0.22 (0.002)
(k)	8	0	0.27 (0.03)	0	0.21 (0.005)
(l)	8	0	0.18 (0.01)	0	0.22 (0.002)
(m)	19	1	0.33 (0.04)	3	0.23 (0.003)
(n)	22	14	0.31 (0.06)	2	0.23 (0.005)
(o)	19	0	0.36 (0.04)	14	0.23 (0.003)
(p)	15	0	0.33 (0.04)	0	0.22 (0.002)
(q)	21	2	0.32 (0.04)	0	0.23 (0.003)
(r)	16	0	0.28 (0.02)	0	0.22 (0.002)
(s)	24	0	0.31 (0.03)	0	0.23 (0.004)
(t)	17	0	0.28 (0.03)	1	0.23 (0.002)
(u)	23	0	0.27 (0.02)	0	0.22 (0.002)

Table 3. TABLE III: Performance of Algorithm 3 over 20 Monte-Carlo runs. Median and std (within brackets) of the metrics ( 53 ) and ( 54 ) are reported.

	$ERR\left(\hat{A}\right)$	$RMSE\left(\hat{A}\right)$
Algorithm 3	4 (0.74)	0.13 (0.008)

Equations

$$\dot{x}(t)=\left(A+\sum _{j=1}^m{u}_j(t)B_j\right)x(t)+Cu(t),t\in \mathbb{R}$$

$$\dot{x}(t)=Ax(t)+Cu(t)+w(t),t\in \mathbb{R}$$

$${\dot{s}}_i(t)$$

$${\dot{f}}_i(t)$$

$${\tau }_i{\dot{v}}_i(t)$$

$${\tau }_i{\dot{q}}_i(t)$$

$$y_i(t)$$

$$+e_i(t)$$

$$k_1=4.3{\vartheta }_0{\rho }_i{T}_E,k_2=\epsilon r_0{\rho }_i{T}_E,k_3=1-\epsilon .$$

$$\left\{\begin{array}[]{rcl}\dot{z}(t)&=&f(z(t),u(t),w(t),\theta)\\ y(t)&=&g(z(t))+e(t),\qquad\qquad t\in\mathbb{R}\end{array}\right.$$

$$x(k+1)=e^{A{T}_R}x(k)+w(k)$$

$$Q={\sigma }^2{\int }_0^{T_R}{e}^{A\tau }{e}^{A^T\tau }\mathrm{d}\tau .$$

$$y_i(k)=\sum _{l=0}^{s-1}{h}_l{x}_i(k-l)$$

$${\Sigma }_h=\frac{1}{N_s}\sum _{j=1}^{N_s}(h^{(j)}-\bar{h})(h^{(j)}-\bar{h}{)}^T.$$

$$h:=H\alpha =\left[\begin{array}{ccccc}\bar{h}& u_1& u_2& \cdots & u_p\end{array}\right]\alpha ,\alpha \in {\mathbb{R}}^{p+1}.$$

$$\mathbf{x}(k)$$

$$\mathbf{A}:=\left[\begin{array}{cc}{e}^{A{T}_R}& 0\\ I_{n(s-1)}& 0\end{array}\right],\mathbf{C}:=(H\alpha {)}^T\otimes I_n,$$

$$\left\{\begin{array}[]{rcl}\mathbf{x}(k+1)&=&\mathbf{A}\mathbf{x}(k)+\mathbf{w}(k)\\ y(k)&=&\mathbf{C}\mathbf{x}(k)+e(k).\end{array}\right.$$

$$\mathbf{w}(k)$$

$$e(k)$$

$$X_{+}$$

$$W$$

$$X_{+}=X{e}^{A^T{T}_R}+W.$$

$$\Delta X=X{A}^T{T}_R+W$$

$$\mathrm{x}=\Phi a+\mathrm{w}$$

$$a^{(j+1)}={\left({\Phi }^T{Q}^{-1}\Phi +[{\Gamma }^{(j)}{]}^{-1}\right)}^{-1}{\Phi }^T{Q}^{-1}\mathrm{x}.$$

$$y(k)=\mathbf{Cx}(k)+e(k),k=1,...,N.$$

$$y(k)=\phi (k)H\alpha +e(k),k=1,...,N$$

$$Y:=\left[\begin{array}{c}y(1)\\ ⋮\\ y(N)\end{array}\right],\mathbf{H}:=\left[\begin{array}{c}\phi (1)\\ ⋮\\ \phi (N)\end{array}\right]H,E:=\left[\begin{array}{c}e(1)\\ ⋮\\ e(N)\end{array}\right]$$

$$Y=\mathbf{H}\alpha +E.$$

$${\mu }_{\alpha }:=[10\cdots 0{]}^T,{\Sigma }_{\alpha }:=\mathrm{diag}(s_1,...,s_p)$$

$$\hat{\alpha }={\mu }_{\alpha }+({\mathbf{H}}^T\mathbf{H}+{\hat{\lambda }}^2{\Sigma }_{\alpha }^{-1}{)}^{-1}{\mathbf{H}}^T(Y-\mathbf{H}{\mu }_{\alpha })$$

$${\hat{\lambda }}^2=\frac{(Y-\mathbf{H}{\hat{\alpha }}_{LS}{)}^T(Y-\mathbf{H}{\hat{\alpha }}_{LS})}{Nn-p-1}$$

$$p(\mathbf{x}(k)|Y,\eta )=\mathcal{N}({\hat{\mathbf{x}}}^s(k),{\mathbf{P}}^s(k)).$$

Full text

Experimental results using rs-fMRI data are shown demonstrating the effectiveness of our approach and comparison with state of the art routines (SPM12 toolbox) is provided.

I Introduction

Recent research in the field of computational neuroscience has shown a growing interest in the estimation of the effective connectivity among brain regions, that is, the existing interactions between neuronal populations located in different brain areas [1]. What makes the problem particularly challenging is the lack of a direct measurement of the neuronal activity: standard neuroimaging techniques (such as EEG, MEG and fMRI) only provide an indirect measurement of such activity. In particular we shall be concerned with fMRI data, which provide noisy measurements of the so-called BOLD (Blood Oxygenation Level Dependent) signal. The dependency between this quantity and neuronal activity is described by the so-called hemodynamic response. In brief, an increment of the synaptic activity in a certain brain region causes an increased blood flow to that area, which in turn augments the ratio of oxygenated haemoglobin w.r.t. the deoxygenated one in that region. By exploiting the different magnetic properties of these two types of hemoglobin, fMRI is able to detect local changes in their relative concentration, thus highlighting the eventual presence of synaptic activations.

Existing methods tackle the problem of effective connectivity estimation by postulating a parametric generative model of the observed fMRI data, which not only accounts for the coupling among neuronal populations, but also for the aforementioned hemodynamic response. Nonlinear dynamical models are used and generally referred to as Dynamic Causal Models (DCMs) [2]: both deterministic [2, 3, 4] and stochastic [5, 6] formulations exist, with the latter being more suited to describe resting state fMRI data.

Having postulated a DCM, effective connectivity inference is reduced to model parameters estimation, which is typically accomplished using Variational Bayes techniques [7, 8, 9]. However, existing approaches are computationally intensive and are effective only on small brain networks. A large portion of the required computational effort is spent for the determination of the structure of the network connecting the different brain regions. Indeed, this choice involves the comparison of several candidate structures, making the problem combinatorial, [10]. In addition, the inversion of stochastic DCMs appears more involved since it also requires the estimation of the hidden random neuronal activity [9].

DCM is a Bayesian approach originally proposed to investigate from fMRI signal the influence on the neuronal system exerts by external stimuli, i.e. the so-called task-evoked effective connectivity. DCM can be also be used to study the directed connections among distributed brain areas when only endogenous brain fluctuations are present, i.e. the so-called resting-state effective connectivity. Our study will focus on this latter problem. We develop a linear DCM, and adapt filtering techniques to estimate the hidden neuronal activity as well as the DCM parameters. Differently from existing methodologies, which require the comparison of a set of candidate patterns of effective connectivity, our approach estimates directly the connectivity pattern, by resorting to an iteratively $\ell_{1}$ reweighed algorithm which induces sparsity in the estimated connectivity matrix, thus automatically detecting the structure of the brain network.

The paper is organized as follows. We briefly review classical DCMs and parameters inference techniques in Sec. II, while in Sec. III we propose an approximated linear stochastic DCM, which is suited for resting state fMRI data. In Sec. IV we deals with the estimation of the DCM parameters (unrealistically) assuming to measure the neuronal activity. We deal with the realistic situation in which only the BOLD signal is available in Sec. V, where we illustrate the joint estimation of the DCM parameters and of the hidden neuronal activity. Using simulated fMRI data, we evaluate the effectiveness of the proposed method in Sec. VI, while we provide some concluding remarks in Sec. VII.

II Problem Formulation

II-A Generative model

A DCM is a nonlinear multi input-multi output (MIMO) dynamic system, with one output for each monitored brain region. A DCM basically consists of two components: the first describes the dynamic coupling among neuronal populations, while the second maps the neuronal activity to the BOLD signal observed in each brain area, through the hemodynamic response. Accordingly, the system states are grouped in two classes: one contains the $n$ neuronal states (one for each monitored brain region), while the other refers to the biophysical quantities involved in the hemodynamic model.

In the seminal paper [2], the coupling among the neuronal states $x\in\mathbb{R}^{n}$ is modelled through a deterministic bilinear interaction:

[TABLE]

where $u(t)=[u_{1}(t)\ \cdots\ u_{m}(t)]^{T}$, with $u_{j}(t)$ denoting the $j$-th designed external stimulus. Effective connectivity in the absence of external solicitations is encoded by the matrix $A\in\mathbb{R}^{n\times n}$, while the bilinear term $B_{j}\in\mathbb{R}^{n\times n}$ accounts for the change in the neuronal coupling due to the $j$-th input; finally, $C\in\mathbb{R}^{n\times m}$ modules the direct influence of external inputs on the neuronal activity. The system matrices $\theta_{c}=\{A,B_{1},B_{2},...,B_{m},C\}$ are unknown parameters which need to be estimated in order to retrieve the effective connectivity among the considered brain areas.

Extensions of the bilinear model (1) can be found in the literature, see e.g. [3] and [4].

In this paper we will consider the resting state activity, which can be modelled as [6]

[TABLE]

where $w(t)$ is a stochastic process which accounts for the unmeasurable random fluctuations that drive the resting state activity. Even if in this situation the exogenous inputs $u(t)$ are typically set to zero, model (2) can still accommodate the presence of external (non-modulatory) signals. In this paper we assume $u(t)=0$.

The hemodynamic response is generally charachterized through an extension of the so-called Balloon-Windkessel model, [11, 12], a SISO dynamic system with the neural activity $x_{i}(t)$, i.e. the component of $x(t)$ in position $i$, as input and the corresponding observed BOLD signal $y_{i}(t)$ as output. The model is described by $4-th$ order nonlinear state space model depending upon the biophysical parameters $\theta_{h}=\{\kappa,\gamma,\tau,\alpha,\rho\}$ more specifically:

[TABLE]

The so called hemodynamic states $\{s_{i},f_{i},v_{i},q_{i}\}$ are biophysical quantities which are affected by the neuronal activity: $s_{i}$ denotes the vasodilatatory signal $s_{i}$, $f_{i}$ is the blood inflow, $v_{i}$ and $q_{i}$ are respectively the blood volume and the deoxyhemoglobin content, $e_{i}(t)$ is the observation error, which may arise from thermal and physiological causes during data acquisition. The output equation (7) depends on the resting blood volume fraction $V_{0}$ (typically set to 0.4) and on the constants $k_{1}$, $k_{2}$ and $k_{3}$. Despite these have found different characterizations in the literature, a general form is derived in [13]:

[TABLE]

Here, $T_{E}$ is the echo time, while $\vartheta_{0}=40.3s^{-1}$ denotes the frequency offset at the outher surface of the magnetized vessel for fully deoxygenated blood at 1.5 T; $\varepsilon\approx 0.4$ is the ratio of intra- and extravascualr signal, while $r_{0}=25s^{-1}$ is the slope of the relation between the intravascular relaxation rate and oxygen saturation.

Collecting the neuronal and the hemodynamic states into the vector $z\in\mathbb{R}^{5n}$ and the two parameters sets in $\theta$, i.e. $\theta=\{\theta_{c},\theta_{h}\}$, a stochastic DCM can be compactly described as the following nonlinear dynamic system:

[TABLE]

where $y(t)=[y_{1}(t)\ \cdots\ y_{n}(t)]^{T}$ and similarly for $e(t)$. The practical implementation of (9) is as follows 111See http://www.fil.ion.ucl.ac.uk/spm/.. The simulated BOLD data, equivalently fMRI data, are given by discretizing the Balloon model (3)-(7) as well as the trajectory of the neural activity in (2) using the sampling time $T_{R}=2$s, which is the time often used for the data acquisition. With some abuse of notation, $x(k):=x(kT_{R})$ and $y(k):=y(kT_{R})$, with $k\in\mathbb{Z}$, be, respectively, the discretized trajectory of $x$ and $y$. By (2), with $u(t)=0$, it is not difficult to see that

[TABLE]

where $w(k)\sim\mathcal{N}(0,Q)$ and

[TABLE]

II-B Effective connectivity estimation

When considering the resting state condition, the effective connectivity estimation problem consists in estimating $\theta_{c}$, i.e. matrix $A$ (and possibly $C$), using the measured BOLD signal $y(k)$, $k=1\ldots N$. It is worth noting that the estimation of $\theta_{c}$ requires to estimate $\theta_{h}$ as well, i.e. the parameters of the Balloon model. This appears as a challenging task, because of the nonlinearity of system (9), but also of the unavailability of the random states $z$, which have to be somehow inferred. Previous works have tackled this problem following two different routes.

A first approach, referred to as spectral DCM, [6], estimates the second-order statistics (i.e. cross-spectra) of the hidden states $z$, which are assumed to be stationary. In this way, the inference of the hidden time-series $z$ is avoided. Furthermore, by parametrizing the spectral density of both process and measurement noise (typically as a power law distribution), the DCM is made deterministic and its parameters are estimated through standard Bayesian approaches [2, 8], postulating Gaussian priors on the parameters $\theta$. Typically, a stability-inducing prior is formulated for the matrix $A$, while empirical priors are assigned to the hemodynamic parameters $\theta_{h}$. However, the nonlinear nature of system (9) complicates the analytical computation of the posterior probability density function, as well as of the marginal likelihood. To overcome this bottleneck, a variational Bayes scheme is used, under the mean-field and the Laplace assumptions: namely, the approximated marginal posterior distributions are considered independent and Gaussian (see [7, 9] for details).

A second family of approaches still parametrizes the distribution of the process and observation noises in (9), but operates in the time-domain, thus estimating not only the model parameters $\theta$ but also the hidden neuronal states $z$. The most common methods are Dynamic Expectation Maximization (DEM) [8] and Generalized Filtering (GF) [14, 5]: even if both adopt the Variational Bayes procedure, DEM uses the mean-field and the Laplace approximations, while GF only exploits the latter.

Despite these techniques represent the state-of-the-art for the estimation of brain effective connectivity, computational feasibility limits their application to small networks (typically, no more than 10 regions). In particular, this drawback severely affects DEM and GF [15]. The reformulation of the nonlinear stochastic DCM (9) as a linear stochastic dynamic model may represent a way to address this issue. In this paper we pursue this idea, by deriving an approximated linear DCM according to the procedure detailed in Sec. III. This reformulation will allow us to jointly estimate the hidden neuronal states $z$ and the DCM parameters $\theta$. It is worth noting that a recent contribution in this direction focuses on deterministic DCMs and transforms the neuronal state equation (2) (with $w(t)=0$) into the frequency domain; then, using a fixed hemodynamic response, the inversion of a deterministic DCM is formulated as a Bayesian linear regression problem, [16].

III Statistical Linearization of DCMs

We model the hemodynamic response as a Finite Impulse Response (FIR) model which takes as input only a neuronal state $x_{i}(k)$ and as output the BOLD signal $y_{i}(k)$,

[TABLE]

where the length $s$ of the impulse response $h:=[h_{1}\ \cdots\ h_{s}]^{T}$ is chosen large enough to retain all the relevant dynamics components. A prior on the impulse response $h$ is derived, using statistical linearisation techniques, from the Baloon model (3)-(7) where the parameters $\theta_{h}=\{\kappa,\gamma,\tau,\alpha,\rho\}$ are described by prior distributions typically used in the literature [2]. This will lead to linear model of the form $h=H\alpha$ (see eq. (14) below) where $\alpha$ will have to be estimated under a Bayesian perspective using the available fMRI data. To to so, we first compute a population of typical responses generated by the Balloon model (3)-(7), then we define $h$ as the linear combination of their empirical mean and of the first $p$ principal components of their sample covariance matrix. Specifically:

Sample $\theta_{h}^{(j)}$, $j=1,...,N_{s}$ from the empirical Gaussian distributions reported in Table 1 of the seminal paper [2]. 2. 2.

For each $\theta_{h}^{(j)}$ compute the output of the Balloon model (3)-(7), where $x_{i}(k)=\delta(k)$ and $\delta(k)$ is the Kronecker delta function. Let $h^{(j)}$ be the corresponding output truncated at length $s$, which represents a statistic of $h$. 3. 3.

Compute the empirical mean $\bar{h}=\frac{1}{N_{s}}\sum_{j=1}^{N_{s}}h^{(j)}$. 4. 4.

Compute the empirical covariance matrix $\Sigma_{h}\in\mathbb{R}^{s\times s}$:

[TABLE] 5. 5.

Compute the eigenvalues decomposition of $\Sigma_{h}$, $\Sigma_{h}=USU^{T}$ , where $S:=\mathrm{diag}(s_{1},...,s_{s})$ and $U:=[u_{1}\ u_{2}\ \cdots\ u_{s}]$. 6. 6.

Define $h$ as

[TABLE]

Figure 1 illustrates the empirical mean $\bar{h}$ of the sampled population of hemodynamic responses and the singular values of their sample covariance matrix $\Sigma_{h}$, showing that there is limited variability in the generated responses.

By introducing the extended state $\mathbf{x}(k)\in\mathbb{R}^{ns}$

[TABLE]

the matrices

[TABLE]

and in view of (10), we can write the linear DCM as a classical stochastic state-space model

[TABLE]

Here, $\mathbf{w}(k):=[w^{T}(k)\ \mathbf{0}]^{T}\in\mathbb{R}^{ns}$ and we assume

[TABLE]

where $\mathrm{blkdiag(\cdot)}$ is the block-diagonal operator and the covariance matrix $Q$ has been defined in (11). Notice that model (17) only depends on the neuronal states $x\in\mathbb{R}^{n}$ and on the parameters $\eta=\{A,\alpha,\sigma,\lambda\}$. Next sections will describe how these are estimated, starting from the simpler (but unrealistic) situation in which the neuronal states $x$ are measurable, Sec. IV, to the situation in which $x$ are not measurable, Sec. V.

IV Parameters Estimation: Measured Neuronal Activity

In this Section we unrealistically assume that observations of the neuronal activity (encoded in the states $x$) are available and we detail how the parameters defining model (17) are estimated.

IV-A Estimation of matrix $A$

The matrix $A$ is expected to be sparse, encoding the fact that each $x_{i}(t)$ is influenced by the activities of few other areas $x_{j}(t)$. To induce sparsity in the estimated $A$, we use a reweighted $\ell_{1}$ method, [17, 18], originally proposed for linear regression models. However, because of the nonlinearity (in the parameters) of model (10), such approach cannot be straightforwardly applied but needs to be properly adapted.

Introducing the matrices

[TABLE]

we rewrite (10) in the non-linear regression form

[TABLE]

which can be linearised exploiting $e^{A^{T}T_{R}}\simeq I+A^{T}T_{R}$:

[TABLE]

where $\Delta X=X_{+}-X$. Using the vectorization operator, we can rewrite (32) as follows

[TABLE]

where $\mathrm{x}:=\mathrm{vec}(\Delta X)$, $\Phi:=(I\otimes X)T_{R}$, $a:=\mathrm{vec}(A^{T})$, $\mathrm{w}:=\mathrm{vec}(W)$. Adopting the Sparse Bayesian Learning perspective [19], we model $a\sim\mathcal{N}(0,\Gamma)$ with $\Gamma:=\mathrm{diag}(\gamma_{1}\ldots\gamma_{n^{2}})$ and independent of $\mathrm{w}$, so that the covariance matrix of $\mathrm{x}$ is $\mathrm{Var}[\mathrm{x}]=\Phi\Gamma\Phi+Q\otimes I_{N}$. Direct application of the reweighted $\ell_{1}$ method [17] to the approximated linear model (33) is not possible here, since $T_{R}=2$s is far to be small, (33) is only a crude approximation of the continuous time model. For this reason, we propose Algorithm 1

where $a^{(j)}=\mathrm{vec}(A^{(j)})$; $a_{i}^{(j)}$ is the entry in position $i$ in vector $a^{(j)}$; $\Gamma^{(j)}=\mathrm{diag}(\gamma_{1}^{(j)}\ldots\gamma_{n^{2}}^{(j)})$; $\phi_{i}$ is the $i$-th column of $\Phi$. Here the update for $\gamma_{i}$’s is the same as [17, formula (29)] while the main difference rests in update for $A^{(j)}$, which is set as the posterior mean of (31), i.e. using the exact model, rather than the posterior mean of the approximated model (33); more precisely:

[TABLE]

Notice that, the search of the optimal $A^{(j)}$ is restricted to the set $\mathcal{A}$. In this paper, $\mathcal{A}$ is chosen as the set of stable matrices, i.e. with eigenvalues having negative real part.

Remark 1

It is worth to stress that the proposed reweighted scheme automatically selects a sparsity pattern for the matrix $A$. This represents a significant computational savings w.r.t. the Variational Bayes methods mentioned in Sec. II-B, which instead have to be provided with several candidate sparsity patterns, among which the best is selected resorting to model selection techniques, e.g. using the free-energy [10]; this gives rise to a combinatorial explosion of candidates making practical inference infeasible.

IV-B Estimation of the coefficients $\alpha$

In Sec. III, a linear model for the hemodynamic response has been derived. As previously remarked, this depends on some parameters $\alpha$ that have to be estimated from data. Here we assume to have measured $N$ data-pairs $(x(k),y(k))$ related by the output equation of model (17), that is

[TABLE]

The previous equation can be rewritten as

[TABLE]

by defining $\mathrm{vec}(\varphi(k)):=\mathbf{x}(k)$ and observing that $\mathbf{C}\mathbf{x}(k)=\mathrm{vec}(\varphi(k)H\alpha)=\varphi(k)H\alpha$. Consequently, introducing the matrices

[TABLE]

we have the following linear regression form

[TABLE]

We consider the Gaussian prior $\alpha\sim\mathcal{N}(\mu_{\alpha},\Sigma_{\alpha})$, where

[TABLE]

and $s_{i}$, $i=1,...,p$, are the first $p$ singular values of the empirical covariance matrix $\Sigma_{h}$ in equation (13). Recalling that $e(k)\sim\mathcal{N}(0,\lambda^{2}I_{n})$ and exploiting standard rules on conditioned Gaussian variables, we easily derive the MAP estimate of $\alpha$ as

[TABLE]

where

[TABLE]

and $\hat{\alpha}_{LS}=(\mathbf{H}^{T}\mathbf{H})^{-1}\mathbf{H}^{T}Y$.

V Parameters Estimation: Unmeasured Neuronal Activity

This section deals with the realistic situation in which only the resting-state fMRI data are observed and the neuronal states $x$ are hidden. Sec. V-A will describe how we estimate $x$ starting from the measured BOLD signal $y$, having fixed the model parameters $\eta=\{A,\alpha,\sigma,\lambda\}$; the subsequent Sec. V-B proposes two algorithms for the combined estimation of $x$ and $\eta$.

V-A Neuronal Activity Estimation

We assume that we are provided with $N$ observations of the BOLD signal $y$, modelled according to the linear dynamic model (17), specified by the parameters $\eta$ (here supposed to be known). We use the Rauch-Tung-Striebel smoother (RTSS) [20] to estimate the hidden neuronal states $x(k),\ k=1,...,N$, and in particular, the closed form smoothing solution to the linear model (17), that is

[TABLE]

The implementation of the RTSS is summarized in Algorithm 2.

V-B Parameters Estimation

We turn to the problem of estimating the parameters $\eta$ specifying the linear model (17). Adopting the Bayesian perspective, we would like to determine the value of $\eta$ which maximizes the marginal posterior

[TABLE]

where $\mathbf{X}:=[\mathbf{x}^{T}(0)\ \cdots\ \mathbf{x}^{T}(N)]^{T}$. However, the computation of such a high-dimensional integral is typically avoided by exploiting the decomposition $p(\eta|Y)\propto p(Y|\eta)p(\eta)$ and resorting to an iterative algorithm to compute the likelihood

[TABLE]

Here we assume $p(\eta)\propto p(A)p(\alpha)p(\sigma)p(\lambda)$, where $p(\sigma)$ and $p(\lambda)$ are uninformative priors, while $p(\mathrm{vec}(A^{T}))\sim\mathcal{N}(\mathbf{0},\Gamma)$ and $p(\alpha)\sim\mathcal{N}(\mu_{\alpha},\Sigma_{\alpha})$.

In this work we use the Expectation-Maximization (EM) algorithm to maximize $\ln p(Y|\eta)$, which is equivalent to maximizing $p(Y|\eta)$. The EM method performs such optimization by iteratively maximizing the lower bound

[TABLE]

of $\ln p(Y|\eta)$, w.r.t. an arbitrary distribution $q(\mathbf{X})$ and $\eta$. In the statistical learning literature $\mathcal{F}(q(\mathbf{X}),\eta)$ is also known as free-energy. At the $l$-th iteration of the algorithm, $\mathcal{F}(q(\mathbf{X}),\eta^{(l)})$ is maximized by $q^{(l+1)}(\mathbf{X})=p(\mathbf{X}|Y,\eta^{(l)})$. Plugging this into (45), one obtains

[TABLE]

which now needs to be maximized w.r.t. $\eta$. Noticing that the second term of $\mathcal{F}(q^{(l+1)}(\mathbf{X}),\eta)$ does not depend on $\eta$, this coincides with optimizing

[TABLE]

w.r.t. $\eta$. Using the Markovian property of system (17), $\mathcal{Q}(\eta,\eta^{(l)})$ can be rewritten as ([21], Ch.12)

[TABLE]

Plugging the smoothing distributions $p(\mathbf{x}(k)|Y,\eta^{(l)})=\mathcal{N}(\hat{\mathbf{x}}^{s}(k),\mathbf{P}^{s}(k))$ (see Eq. (42)) and the pairwise smoothing distributions ([21], Ch.12)

[TABLE]

into (48) we get

[TABLE]

where

[TABLE]

The complete routine we adopt to estimate the parameters $\eta$ through the EM method is summarized in Algorithm 3. The procedure is initialized with $\eta^{(0)}=\{A^{(0)},\sigma^{(0)},\alpha^{(0)},\lambda^{(0)}\}$:

[TABLE]

where $\bar{Y}=\frac{1}{Nn}\sum_{i=1}^{Nn}Y_{i}$.

Notice that, at Step 7 of Algorithm 3, the search of $\eta$ is restricted to the set $\Omega=\{\mathcal{A},\ \mathbb{R}_{+},\ \mathbb{R}^{p+1},\ \mathbb{R}_{+}\}$. Since this optimization may be computationally intensive, we have also tested a simpler method, where the uncertainty of the estimates returned by the smoother (encoded in the matrix $\mathbf{P}^{s}(k)$) is ignored. Briefly, this alternative routine iteratively estimates $A$ and $\sigma$ through Algorithm 1, $\alpha$ and $\lambda$ through Eqs. (40), (41), replacing the hidden neuronal states with the smoother estimates $\hat{\mathbf{x}}^{s}(k)$. The results are inferior to those of Algorithm 3 and thus not reported.

VI Experimental Results

The methods illustrated in Secs. IV and V are experimentally evaluated using simulated rs-fMRI data. In particular: (a) we analyze the effectiveness of the reweighted $\ell_{1}$ technique detailed in Sec. IV-A in detecting the correct brain network structure and in estimating the strength of the neuronal connections; (b) we evaluate the performance of Algorithm 3 which only uses the resting state BOLD signals, when the true network structure has to be inferred from the fMRI data. In the latter setup, the proposed methods are compared with the Variational Bayes routines implemented in the SPM12 toolbox (http://www.fil.ion.ucl.ac.uk/spm/).

VI-A Experimental Setup

The rs fMRI data used in these experiments are generated through the SPM12 routines spm_int_J, spm_fx_fmri and spm_gx_fmri. Specifically, we consider a brain network with $n=7$ neuronal nodes, whose activity is modulated according to Eq. (2), where $w(t)$ is a Gaussian white noise process with variance $\sigma^{2}=0.01$ and

[TABLE]

$N=600$ data are generated with sampling time $T_{R}=2$s. No observation noise is added to the BOLD signal generated by the routine spm_gx_fmri. It is worth stressing that the self-connections are a-priori fixed to $-0.5$ to prevent instability issues.

VI-B Performance Evaluation

Let $\hat{A}$ be an estimate of the effective connectivity matrix: we define the Root-Mean-Squared-Error (RMSE)

[TABLE]

where the notation $\underline{A}$ denotes the matrix $A$ with its diagonal set to 0. Indeed, the diagonal elements of $A$ (i.e. the self-connections) do not give any information on the effective connectivity.

A second performance index measures the number of errors committed by a certain method in the recovery of the true sparsity pattern of matrix (52). Namely,

[TABLE]

where $SP(\cdot):\mathbb{R}^{n\times m}\rightarrow\mathbb{R}^{n\times m}$ and

[TABLE]

VI-C Estimation of effective connectivity using measured neuronal activity

In this section we assess the validity of the reweighted $\ell_{1}$ method illustrated in Sec. IV-A. To this purpose we pass as inputs to Algorithm 1 the neuronal time-series $x(k)$, $k=1,...,600$, generated by the SPM12 routine spm_int_J when provided with spm_fx_fmri and the matrix (52). In addition, we compare Algorithm 1 with the linear reweighted $\ell_{1}$ based on approximation (33).

Table I compares the results achieved by Algorithm 1 and by its linear version on the data generated in 50 Monte-Carlo runs, using different realizations of the noise $w(t)$: the medians (and the standard deviations) of indexes (53) and (54) are reported. Before evaluating these indexes, a tresholding is applied to the estimated matrices $\hat{A}$: namely the entries with absolute value less than 0.1 are set to zero.

The results reported in Table I show that Algorithm 1 benefits from using the true generative model at Step 6, in place of the linear approximation (33).

VI-D Estimation of effective connectivity using rs fMRI data

We now evaluate the performance of the proposed Algorithm 3 in the recovery of matrix (52) when the BOLD data $y(k)$, $k=1,...,600$, generated under the setup of Sec. VI-A, are available. In particular, we conduct a Monte-Carlo study over 20 different realizations of the noise $w(k)$. The performance of Algorithm 3 is compared with the one returned by the SPM12 routines spm_dcm_fmri_csd and spm_dcm_estimate, both provided with null exogenous inputs. These functions implement two of the methods mentioned in Sec. II-B: specifically, spm_dcm_fmri_csd runs the spectral DCM (sDCM) [6], while spm_dcm_estimate is invoked activating the option DCM.options.stochastic in order to implement the DEM procedure [8].

To guarantee a fair comparison between our methods and the SPM12 routines, the latter should perform model selection among all the possible sparsity patterns of matrix $A$. However, this would require the estimation of a combinatorial number of matrices, making this type of approach computationally infeasible. Attempting to still conduct a fair comparison, we pass 21 candidate sparsity patterns $\mathsf{m}_{i}$, $i=1,...,21$, to the SPM12 routines, including the true one (denoted by (a)). These are illustrated in Figure 2. For each of them, Table II reports the value of index (54), as well as the median and the standard deviation of index (53) corresponding to the matrices estimated by the routines spm_dcm_fmri_csd and spm_dcm_estimate. In addition, the columns denoted with $\#Chosen$ report the number of times each model is chosen according to the free-energy metric. The most frequently selected structures are very different from the true one (as indicated by the value of the metrics $ERR(\hat{A})$). According to our experimental experience, the SPM routines tend to select structures with a low degree of sparsity.

Table III contains the median and the standard deviation of indexes (53) and (54) computed for the matrices $\hat{A}$ estimated by the proposed Algorithm 3. The performance of our approach is superior to the one achieved by the SPM12 routines, both in terms of $RMSE(\hat{A})$ and $ERR(\hat{A})$. In addition, the sparsity patterns detected by Algorithm 3 differ at most by 7 entries from the true one. It is worth noting that $RMSE(\hat{A})$ returned by Algorithm 3 is comparable to that achieved by spectral DCM when the latter is provided with the true sparsity pattern.

On the conducted study, Alg. 3 takes an average computational time equal to 558s, while spm_dcm_fmri_csd and spm_dcm_estimate respectively take on average 721s and 4044s to process all the 21 patterns.

VII Conclusions and Future Work

We have proposed a new approach for the estimation of the effective connectivity in brain networks using resting state fMRI data. Existing methods tackle the problem by formulating a parametric generative model (DCM) of the observed data and by inverting it using Bayesian approaches. However, due to the nonlinear and stochastic nature of the standard DCMs, the current techniques are highly computationally intensive. In addition, when no a-priori knowledge is available on the structure of the underlying network connecting the monitored brain regions, these approaches require to compare several candidate interaction patterns in order to choose the most plausible one.

The technique we have developed is based on a linear stochastic generative model (DCM), which we have built by linearly approximating the Balloon model describing the hemodynamic response. Exploiting the linearity of the derived DCM we have applied classical smoothing techniques to estimate both the hidden neuronal activity and the model parameters starting from the observed BOLD signal. A key feature of our algorithms is the ability to automatically detect the structure of the underlying brain network, thanks to the use of a sparsity inducing method, thus avoiding a combinatorial search over candidate structures.

Future work will include an extended experimental study of the proposed method and, in particular, its application to real fMRI data.