A new causal centrality measure reveals the prominent role of subcortical structures in the causal architecture of the extended default mode network

Abstract

Network representation has been an incredibly useful concept for understanding the behavior of complex systems in social sciences, biology, neuroscience, and beyond. Network science is mathematically founded on graph theory, where nodal importance is gauged using measures of centrality. Notably, recent work suggests that the topological centrality of a node should not be over-interpreted as its dynamical or causal importance in the network. Hence, identifying the influential nodes in dynamic causal models (DCM) remains an open question. This paper introduces causal centrality for DCM, a dynamics-sensitive and causally-founded centrality measure based on the notion of intervention in graphical models. Operationally, this measure simplifies to an identifiable expression using Bayesian model reduction. As a proof of concept, the average DCM of the extended default mode network (eDMN) was computed in 74 healthy subjects. Next, causal centralities of different regions were computed for this causal graph, and compared against several graph-theoretical centralities. The results showed that the subcortical structures of the eDMN were more causally central than the cortical regions, even though the graph-theoretical centralities unanimously favored the latter. Importantly, model comparison revealed that only the pattern of causal centrality was causally relevant. These results are consistent with the crucial role of the subcortical structures in the neuromodulatory systems of the brain, and highlight their contribution to the organization of large-scale networks. Potential applications of causal centrality—to study causal models of other neurotypical and pathological functional networks—are discussed, and some future lines of research are outlined.

Appendices

Appendix 1: Causal centrality derivation

In this appendix, the identifiable expression in Eq. 6 is derived from the definition of causal centrality in Eq. 4, and the normalization constant (Z in Eq. 4) is clarified in this process. To start, we expand the KL-divergence between the pre- and post-intervention distributions:

$$\begin{aligned} D_{{\text{KL}}} \left( {P\left( {\theta ,Y|m} \right)\parallel P^R \left( {\theta ,Y|m^R } \right)} \right) & = \smallint \nolimits_\theta P\left( {Y,\theta |m} \right) \log \frac{{P\left( {Y,\theta |m} \right)}}{{P^R \left( {Y,\theta |m^R } \right)}} \\ & = \smallint \nolimits_\theta P\left( {\theta |Y,m} \right) P\left( {Y|m} \right) \log \frac{{P\left( {\theta |Y,m} \right) P\left( {Y|m} \right) }}{{P^R \left( {\theta |Y,m^R } \right) P^R \left( {Y|m^R } \right)}} d\theta \\ \end{aligned}$$

(11)

The first equality is just the definition of KL-divergence: $D_{{\text{KL}}} \left( {P_1 ||P_2 } \right) \triangleq \smallint \limits_z P_1 \left( z \right)\log \frac{P_1 \left( z \right)}{{P_2 \left( z \right)}}{\text{d}}z$. The second equality uses the product rule of probability: $P\left( {Y,\theta |m} \right) = P\left( {\theta |Y,m} \right) P\left( {Y|m} \right)$. Herein, R as a superscript refers to the reduced/post-intervention model, whereas the full/pre-intervention model bears no superscript. With some re-arrangement, Eq. 11 simplifies to:

$$\begin{aligned} & = P\left( {Y|m} \right) \smallint \limits_\theta P\left( {\theta |Y,m} \right)\left[ {\log \frac{{P\left( {\theta |Y,m} \right) }}{{P^R \left( {\theta |Y,m^R } \right)}} + \log \frac{{P\left( {Y|m} \right)}}{{P^R \left( {Y|m^R } \right)}}} \right] {\text{d}}\theta \\ & = P\left( {Y|m} \right) \left[ { \smallint \limits_\theta P\left( {\theta |Y,m} \right)\log \frac{{P\left( {\theta |Y,m} \right) }}{{P^R \left( {\theta |Y,m^R } \right)}} {\text{d}}\theta + \smallint \limits_\theta P\left( {\theta |Y,m} \right) \log \frac{{P\left( {Y|m} \right)}}{{P^R \left( {Y|m^R } \right)}} {\text{d}}\theta } \right] \\ & = P\left( {Y|m} \right)\left[ { D_{{\text{KL}}} \left( {P\left( {\theta |Y,m} \right) P^R \left( {\theta |Y,m^R } \right)} \right) + \log \frac{{P\left( {Y|m} \right)}}{{P^R \left( {Y|m^R } \right)}} \smallint \limits_\theta P\left( {\theta |Y,m} \right){\text{d}}\theta } \right] \\ & = P\left( {Y|m} \right) \left[ { D_{{\text{KL}}} \left( {P\left( {\theta |Y,m} \right) P^R \left( {\theta |Y,m^R } \right)} \right) + \log \frac{{P\left( {Y|m} \right)}}{{P^R \left( {Y|m^R } \right)}} \times 1} \right] \\ \end{aligned}$$

(12)

In variational Bayesian inference, the optimized variational density and free energy replace the posterior density and log model evidence: $Q\left( {\theta |Y,m} \right) \approx P\left( {\theta |Y,m} \right)$ and $F \approx \log P\left( {Y|m} \right)$. Hence, Eq. 12 simplifies to:

$$\approx \exp \left( F \right)\left[ {D_{{\text{KL}}} \left( {Q\left( {\theta |Y,m} \right)\parallel Q^R \left( {\theta |Y,m^R } \right)} \right) + \left( {F - F^R } \right)} \right]$$

(13)

Causal centrality is estimable as the normalized version of Eq. 13, where the normalization constant $Z = P\left( {Y|m} \right) \approx \exp \left( F \right)$. As such, we arrive at Eq. 6:

$$\begin{aligned} {\text{Causal centrality}} & \triangleq \frac{1}{{P\left( {Y|m} \right)}}D_{{\text{KL}}} \left[ {P\left( {\theta ,Y|m} \right)\parallel P^R \left( {\theta ,Y|m^R } \right)} \right] \\ & \approx D_{{\text{KL}}} \left[ {Q\left( {\theta |Y,m} \right)\parallel Q^R \left( {\theta |Y,m^R } \right)} \right] + \left( {F - F^R } \right) \\ \end{aligned}$$

(14)

where the posterior and free energy of the reduced model can be computed analytically using Bayesian model reduction, as elaborated in “Appendix 2”. Moreover, under Gaussian assumptions on the posterior distributions in variational Laplace (Friston et al. 2007; Zeidman et al. 2023), the KL-divergence between the two (full and reduced) Gaussian posteriors can be readily computed as:

$$\begin{aligned}D_{{\text{KL}}} \left( {{\rm{\mathcal{N}}}_{\text{F}}||{\rm{\mathcal{N}}}_{\text{R}} } \right) &= \frac{1}{2}\left\{{\text{tr}}\left( {\Sigma_{\text{R}}^{ - 1} \Sigma_{\text{F}} }\right) - d \right.\\ &\quad\qquad\left.+ \left( {\mu_{\text{R}} -\mu_{\text{F}} } \right)^{\text{T}} \Sigma_R^{ - 1} \left({\mu_{\text{R}} - \mu_{\text{F}} } \right) + \ln \left( {\frac{{\det \Sigma_{\text{R}} }}{{\det \Sigma_{\text{F}} }}} \right) \right\}\end{aligned}$$

(15)

where $\mu$ and ${\Sigma }$ denote the mean vector and the covariance matrix, respectively; $d$ is the dimensionality of the multivariate Gaussians; $F$ and $R$ subscripts refer to the full and reduced models, respectively; tr stands for the trace of a matrix, and $\ln = \log_{\text{e}}$.

Appendix 2: Bayesian model reduction derivations

Bayesian model reduction refers to the analytic inversion of reduced models using the priors and posteriors of a full model. Reduced models are nested within a full model; that is, they include only a subset of the parameters of the full model after “switching off” the other parameters (by imposing very precise null priors, which shrinks them to zero). Consider Bayes rule replicated for the full and reduced models (denoted by $m$ and $m^R$, respectively) (Friston et al. 2016, 2018; Friston and Penny 2011):

$$P\left( {\theta |y,m} \right) = \frac{{P\left( {y|\theta ,m} \right)P\left( {\theta |m} \right)}}{{P\left( {y|m} \right)}}$$

(16)

$$P^R \left( {\theta |y,m^R } \right) = \frac{{P\left( {y|\theta ,m^R } \right)P\left( {\theta |m^R } \right)}}{{P\left( {y|m^R } \right)}}$$

(17)

Since the models differ only in terms of their priors, the likelihood terms are identical:

$$P^R \left( {y|\theta ,m^R } \right) = P\left( {y|\theta ,m} \right)$$

(18)

Hence, equating the expressions in Eqs. 16 and 17 over the likelihood gives:

$$\frac{{P^R \left( {\theta |y,m^R } \right)P^R \left( {y|m^R } \right)}}{{P^R \left( {\theta |m^R } \right)}} = P\left( {y|\theta ,m} \right) = \frac{{ P\left( {\theta |y,m} \right)P\left( {y|m} \right)}}{{P\left( {\theta |m} \right)}}$$

(19)

By re-arranging Eq. 19, we get the posterior distribution over the parameters of the reduced model:

$$P^R \left( {\theta |y,m^R } \right) = P\left( {\theta |y,m} \right)\frac{{P^R \left( {\theta |m^R } \right)}}{{P\left( {\theta |m} \right)}}\frac{{P\left( {y|m} \right)}}{{P^R \left( {y|m^R } \right)}}$$

(20)

On the R.H.S. of Eq. 20 only $P^R \left( {y|m^R } \right)$ is unknown, which is the model evidence for the reduced model. This term can be obtained by integrating both sides of Eq. 20 over $\theta$:

$$\int {P^R } \left( {\theta |y,m^R } \right) {\text{d}}\theta = 1 = \int P \left( {\theta |y,m} \right)\frac{{P^R \left( {\theta |m^R } \right)}}{{P\left( {\theta |m} \right)}}\frac{{P\left( {y|m} \right)}}{{P^R \left( {y|m^R } \right)}} {\text{d}}\theta$$

(21)

$$\to \frac{{P^R \left( {y{|}m^R } \right)}}{{P\left( {y{|}m} \right)}} = \smallint^P\left( {\theta {|}y,m} \right)\frac{{P^R \left( {\theta {|}m^R } \right)}}{{P\left( {\theta {|}m} \right)}} {\text{d}}\theta$$

(22)

$$\to \log P^R \left( {y|m^R } \right) = \log P\left( {y|m} \right) + \log \smallint^P\left( {\theta |y,m} \right)\frac{{P^R \left( {\theta |m^R } \right)}}{{P\left( {\theta |m} \right)}} {\text{d}}\theta$$

(23)

$$\to F^R = F + \log \smallint^Q\left( {\theta |y,m} \right)\frac{{P^R \left( {\theta |m^R } \right)}}{{P\left( {\theta |m} \right)}} {\text{d}}\theta$$

(24)

In the last line, the optimized variational density and free energy have replaced the posterior density and log model evidence: $Q\left( {\theta |Y,m} \right) \approx P\left( {\theta |Y,m} \right)$ and $F \approx \log P\left( {Y|m} \right)$. Under Gaussian assumptions for the distributions (in variational Laplace), the reduced posterior and free energy take simple forms, as elaborated in Friston et al. (2016, 2018) and Friston and Penny (2011).

Appendix 3: Graph theoretical centralities

This appendix includes the mathematical expressions for several graph-theoretical centrality measures. These measures are defined for weighted networks, based on Opsahl et al. (2010), which generalizes the definitions originally proposed for binary networks (Freeman 1979).

Node strength is the generalization of node degree (i.e., the number of connections of a node) for weighted networks, which is defined for node $x_i$ as:

$$C_s \left( {x_i } \right) = \mathop \sum \limits_{j = 1}^n \left| {a_{ij} } \right|$$

(25)

where $n$ is the number of nodes, |.| denotes the absolute value function, and $a_{ij}$ is the $\left( {i,j} \right){\text{th}}$ entry of the weighted adjacency matrix $A$.

Closeness centrality is defined as:

$$C_c \left( {x_i } \right) = \left( {\mathop \sum \limits_{j = 1}^n d\left( {i,j} \right)} \right)^{ - 1}$$

(26)

where $d(i,j)$ denotes the shortest path between nodes $x_i$ and $x_j$, which is the path that minimizes the cost of travelling from $x_i$ to $x_j$; the cost of travelling is the inverse of the weight of the (directed) edges connecting the two nodes—as encoded in the adjacency matrix.

Betweenness centrality is defined as:

$$C_B \left( {x_i } \right) = \mathop \sum \limits_{j = 1}^n \mathop \sum \limits_{k = 1}^n \frac{{g_{jk} \left( {x_i } \right)}}{{g_{jk} }}$$

(27)

where $g_{jk} (x_i )$ is the number of shortest paths between nodes $x_j$ and $x_k$ that go through node $x_i$; and $g_{jk}$ is the total number of shortest paths between nodes $x_j$ and $x_k$.

Note that in the formulation of DCM (Eq. 7), $A$ encodes the directed causal influences; hence, it corresponds to the adjacency matrix of a (directed weighted) signed graph, in which excitatory (positive) and inhibitory (negative) effects are equally important. In this case, the absolute values of the causal connections are typically computed (Cai et al. 2021; Dablander and Hinne 2019). So, to compute the closeness and betweenness centralities (based on the notion of shortest paths), $A$ was converted to an unsigned matrix first, using the absolute function: $A \to |A|$.

Finally, for a node to have high eigenvector centrality, it should be connected to many other nodes that also have high eigenvector centrality. In other words, connection to an important node counts more than connection to a less important node. This notion is formalized using the first eigenvector of the adjacency matrix (Bonacich 2007; Bonacich and Lloyd 2001):

$$C_E \left( {x_i } \right) = \frac{1}{\lambda } \sum \limits_{j \in N\left( {x_i } \right)} \nu_j$$

(28)

where $\lambda$ is the largest eigenvalue of an undirected adjacency matrix, and $\nu_j$ is the jth element of the corresponding eigenvector; $N(x_i )$ denotes the set of nodes that are adjacent/neighbor to $x_i$. Since eigenvector centrality is defined for undirected graphs, the (unsigned weighted) directed adjacency matrix $\left| A \right|$ was first symmetrized as follows: $\left| A \right| \to \frac{{\left| A \right| + \left| {A^{\text{T}} } \right| }}{2} = \left| A \right|_{{\text{sym}}}$; where $T$ denotes the transpose operator.

The graph-theoretical centrality measures were computed using the centrality function in MATLAB R2021b.