Detecting abnormal connectivity in schizophrenia via a joint directed acyclic graph estimation model

A Bayesian network is a probabilistic graphical model that encodes the random variables as its nodes and their conditional dependencies as its directed edges. The node set $V$ and the edge set $E=V\times V$ make up with a directed acyclic graph (DAG) $G=(E,V)$ in which each edge $(i,j)$ represents a directed edge from variable $i$ to $j$. $i$ is called the parent node of $j$ and $j$ is the child node of $i$ accordingly. Given one variable, say $j$, and the set of all of its parent nodes $Pa_{j}$, we can always write the conditional probability distribution of $j$ as $P(j|Pa_{j})$.Without loss of generality, we can always factorize the joint distribution of random variables $Y=\{y_{1},y_{2},\dots,y_{d}\}$ as $P(y)=\prod_{i=1}^{d}P(y_{i}|Pa_{i})$.

Given the data matrix $X\in\mathbb{R}^{n\times d}$ as the $n$ observations of random variables $y\in\mathbb{R}^{d}$, learning the structure of the Bayesian network or directed acyclic graph (DAG) means finding a proper distribution $P(X)$ defined on the graph $G=(V,E)$. More specifically, we want to determine all the conditional dependencies via the observation $X$. We model these conditional dependencies via a structural equation model(SEM) defined by $X=XW+z$, where $W\in\mathbf{R}^{d\times d}$ is the weighted adjacency matrix and $z$ is the random noise vector. Without further assumption on the noise distribution, we follow the approach in [16] and use the least square loss function for the linear SEM. It should be noted that the loss can be changed to any other smooth functions over $\mathbf{R}^{d\times d}$. Moreover, several previous studies have shown that minimizing the least square loss can guarantee revealing the true DAG with high probability in finite-samples and even in high dimensions setting($d\gg n$) and this result is consistent for both Gaussian SEM[17], [18] and non-Gaussian SEM[19].

In the literature[16], the authors firstly characterized the adjacency matrix of a DAG algebraically and turned the fussy searching problem in traditional GES methods into a continuous optimization issue. Their characterization is formulated with the theorem

where $\circ$ is the Hadamard product and $e^{\cdot}$ is the matrix exponential operator. Fortunately, this function is differentiable with a simple gradient

To overcome the sample size limitation and take advantage of data from multiple groups, we propose a joint DAG estimation model to incorporate the similarity prior in the data. The primal problem can be formulated as

where $l(\mathbf{W};\mathbf{X})$ is the loss function of SEM, i.e. $l(\mathbf{W};\mathbf{X})=\sum_{k=1}^{K}\frac{1}{2n_{k}}||X^{(k)}-X^{(k)}W^{(k)}||_{F}^{2}$. $P(\mathbf{W})$ is the penalty term of the weighted adjacency matrices and usually chosen to encourage $W^{(1)},W^{(2)},\dots,W^{(K)}$ to share certain characteristics such as the similar pattern of nonzero elements. In addition, the sparsity of those matrices are usually encoded as a prior knowledge which benefits both model training and interpretation of the results. Considering the similar causal structure underlining different observations, we borrow the idea of group regularization from the undirected graphical model to encourage the shared DAG structures. This technique has been proven useful in many gene expression network studies. The group regularization term is often formulated as:

where $\lambda_{1}$ and $\lambda_{2}$ are nonnegative parameters. $\lambda_{1}$ controls the sparsity of $W^{(k)}$s and the larger $\lambda_{1}$ is, the sparser the solution of (1) will be. On the other side, $\lambda_{2}$ restrains the patterns of the nonzeros in $W^{(k)}$s and the larger $\lambda_{2}$ is, the more identical $W^{(k)}$s will be. If $\lambda_{1}$ and $\lambda_{2}$ are set to $0$, the model is degraded to the original notears model.

By applying the augmented Lagrange, the primal problem is written as

We use the dual ascent to solve problem (3) and the lagrange multiplier is derived as

Hence, the corresponding dual function is

The dual problem is then

Denote $\mathbf{W}^{*}_{\mathbf{\alpha}}$ as the local minimizer of problem (5), i.e. $g(\alpha_{1},\alpha_{2},\dots,\alpha_{K})=L(\mathbf{W}^{*}_{\mathbf{\alpha}},\alpha_{1},\alpha_{2},\dots,\alpha_{K})$. By noting that $g(\alpha_{1},\alpha_{2},\dots,\alpha_{K})$ is a linear function of $\alpha_{k}$s, the partial derivatives are given by $\partial g(\alpha_{1},\alpha_{2},\dots,\alpha_{K})/\partial\alpha_{k}=h(W^{(k)^{*}})$. Hence we can perform the dual ascent by updating $\alpha_{k}$s with

The rest of the concern will be the unconstrained optimization of subproblem (5). Due to its high dimensionality and non-convexity, we follow the similar idea in [16] and solve it with L-BFGS algorithm[20] when $\lambda_{1}=\lambda_{2}=0$. When $\lambda_{1},\lambda_{2}>0$, the problem can be solved with proximal quasi-Newton (PQN) method[21]. As an iteration algorithm, at the $k$-th step, the descent direction is searched through a quadratic approximation of the smooth term:

where $\mathbf{g}_{k}$ is the gradient of $f(w)$ and $B_{k}$ is the L-BFGS approximation of the Hessian matrix of $f(w)$. In addition, for each coordinate $j$, the solution of problem (8) has a closed form update $\mathbf{d}\leftarrow\mathbf{d}+z^{*}e_{j}$ in which $e_{j}$ is the unit vector in standard basis and

It should be noted that the low-rank structure of $B_{k}$ makes it sparse and fast to compute during the coordinate update. The sparse regularization in the model can further enable us to speed up the computation. Instead of updating for all coordinates, an active set of coordinates can be chosen based on the subgradients and we just need to update regarding the active set. More details of the subproblem optimization can be find in Appendix A.

As it shows in the joint estimation model, there are two parameters, i.e., $\lambda_{1}$ and $\lambda_{2}$ that need to be specified before we use it for causality inference. Different setting of those parameters will lead to different causal structures. Training the model with a large $\lambda_{1}$ will lead to a sparse $\mathbf{W}$ but may underfit the data. On the other hand, the model with a smaller $\lambda_{1}$ gives rise to a denser $\mathbf{W}$ but may overfit the data. Consequently, improper choice of the parameter will lead to inaccurate inference result: a sparser causal graph tends to have more false negatives while a denser graph usually has severe false positives problem. Hence we would like to balance both the simplicity of the model and the goodness of fit. Similarly, the parameter $\lambda_{2}$ reflects how similar the DAG structures will be between groups. A larger $\lambda_{2}$ implies more similar structures between groups and vice versa. When $\lambda_{2}$ is trivially set to $0$, the joint estimation model is then degraded to $K$ seperate NOTEARS models for $K$ groups; As $\lambda_{2}$ goes to $\infty$, group variations will be wiped out gradually and the resulted $K$ DAGs will be identical to each other which is exactly the case that we regard the data as one group and apply the NOTEARS model directly. To find a proper combination of $(\lambda_{1},\lambda_{2})$, we did a 5-fold cross-validation using the value of $l(\mathbf{W})$ as the performance measure.

The model was applied to fMRI data in the MIND Clinical Imaging Consortium (MCIC) from Mind Research Network(MRN,www.mrn.org). The fMRI data were collected from 208 subjects, among them 92 SZ patients (age: 34 $\pm$ 11, 22 females) and 116 healthy controls (age: 32 $\pm$ 11, 44 females)during a sensory motor task, a block design motor response to auditory stimulation. Specifically, the images were acquired on a Siemens3T Trio Scanner and 1.5 T Sonata with echo-planar imaging (EPI) sequences taking parameters (TR $=2000ms$, TE $=30ms(3.0T)/40ms(1.5T)$, field of view = $22cm$, slice thickness = $4mm$, $1mm$ skip, 27 slices, acquisition matrix = 64 $\times$ 64, flip angle = $90^{\circ}$). Then the data were pre-processed with SPM5 (http://www.fil.ion.ucl.ac.uk/spm) and were realigned, spatially normalized and resliced to $3\times 3\times 3mm$, smoothed with a $10\times 10\times 10mm^{3}$ Gaussian kernel, and analyzed by multiple regression considering the stimulus and their temporal derivatives plus an intercept term as regressors. Finally the stimulus-on versus stimulus-off contrast images were extracted with $53\times 63\times 46$ voxels and all the voxels with missing measurements were excluded resulting in 41236 voxels. In order to filter irrelevant information, we further implemented a multiple t-test between case and control groups at the voxel level. Finally, $p=9816$ voxels were left for analysis and 116 ROIs were extracted based on the Automated Anatomical Labeling brain atlas. For each ROI, we then averaged the beta values of the voxels within that ROI and finally got an $116\times 208$ data matrix.

The application of functional connectivity to study schizophrenia appears to be a very powerful approach since it provides a way to study aberrant connectivity between sets of regions (networks), which is thought to be a core feature of schizophrenia. However, the directed functional connectivity is rarely used as biomarkers in the schizophrenia study. From the comparison between two groups, we found that the directed brain network is significantly denser in schizophrenia patients. Moreover, the mean clustering coefficient of schizophrenia patients is smaller than that of the healthy controls. With a denser network, we observed that both the global and local efficiency of schizophrenia patients are significantly smaller than those of healthy controls. The functional integration represents the capability to combine information from distributed brain regions. This demonstrates the directed brain network of schizophrenia patients is less efficient in terms of brain region communication either globally or locally in the auditory task.

Among the hub nodes identified in two groups, we also observed disruptions of the hub structures in schizophrenia patients. The gyrus retus (REC.L), as an in-hub on the left hemisphere, is missing in SZ patients. Previous study[29] has demonstrated that pronounced gray matter volume decline at gyrus rectus has been observed in SZ patients, which supports our findings with a structural foundation. Another study[30] also showed that the reduced functional connectivity has been observed between gyrus rectus and other regions within self-referential processing network. Inferior frontal gyrus, triangular at left hemisphere (IFGtriang.L), has been shown to have increased activity during auditory hallucinations[31], and the Theory of Mind(ToM) deficits [32] in schizophrenia patients. Our result further implies that the hyperconnectivity related to this region is potentially an important biomarker of schizophrenia. The behavior of temporal pole on the right hemisphere (TPOmid.R) is also distinct between HC and SZP, which agrees with the finding in [33], [34].

In addition to the abnormal hub structures in SZ patients, we also have several findings in coordination with previous studies by comparing the CtEs extracted from both groups. For example, previous research[35] revealed that olfactory bulbs disconnectivity of olfactory regions in schizophrenia may account for olfactory dysfunction and disrupted integration with other sensory modalities in SZ patients. The missing connections to the olfactory cortex in CtEs of SZ patients also confirm this findings. In the subcortical regions, pallidum, putamen and caudate work together to communicate with the subthalamic nucleus[36]. Our findings illustrated the missing communications within and between putamen and pallidum, which is in line with the findings in [37],[38].

From Table IV, parahippocampus on both sides are identified as hub nodes in undirected functional connectivity for SZ patients. The parahippocampus plays an important role in the encoding and recognition of environmental scenes. Previous research[39] have discovered the abnormal connections between parahippocampus and temporal pole in early stage SZ patients. Several other researches also revealed the aberrant increased connectivity between parahippocampus and other limbic areas[40], including hippocampus[41].

When we compared the finding of directed graphical model and undirected graphical model, we found that the overlapped brain networks indicated by both methods cannot represent the feature of whole brain networks, especially the connection variations between SZ patients and healthy controls. This implies that the features from undirected networks and directed networks may not be regarded as the similar features in depicting the interactions between brain regions. In fact, the underlying distributions of two graphical modelling framework are intrinsically different. Figure 7 demonstrates the relationship between two models. Only the distributions falling in the intersection of directed graphs (D) and undirected graphs (U) are perfectly mapped through both graphical models. It is clearly that we do not know the underlying distribution of brain activities. We speculate that either method can only provide a sketch of the ground truth from one particular perspective. Hence, we reported the pairwise comparisons results of both methods.

Currently, both the NOTEARS and the proposed method were optimized with second order approximation with a global searching over the feasible space, which is time consuming if we want to extend those methods to application with hundreds of variables, a typical case in brain imaging study. For example, it is impossible to directly apply the method for voxel-wise network analysis. It will be desirable if we can find a way to reduce the search space in the model. In fact, the model itself has a flexible setting to incorporate the prior information of the network in the optimization, which will significantly alleviate the computation burden. The prior information includes some known connections/disconnections and the directions of connections. Another limitation lies in the variances of the results that we found between the directed network method and the undirected network method. At this stage, most of whole brain functional connectivity studies are based on either directed or undirected network analysis. It will be meaningful to cross-validate the findings with more data sets.