Integrating Multimodal and Longitudinal Neuroimaging Data with Multi-Source Network Representation Learning

Abstract

Uncovering the complex network of the brain is of great interest to the field of neuroimaging. Mining from these rich datasets, scientists try to unveil the fundamental biological mechanisms in the human brain. However, neuroimaging data collected for constructing brain networks is generally costly, and thus extracting useful information from a limited sample size of brain networks is demanding. Currently, there are two common trends in neuroimaging data collection that could be exploited to gain more information: 1) multimodal data, and 2) longitudinal data. It has been shown that these two types of data provide complementary information. Nonetheless, it is challenging to learn brain network representations that can simultaneously capture network properties from multimodal as well as longitudinal datasets. Here we propose a general fusion framework for multi-source learning of brain networks – multimodal brain network fusion with longitudinal coupling (MMLC). In our framework, three layers of information are considered, including cross-sectional similarity, multimodal coupling, and longitudinal consistency. Specifically, we jointly factorize multimodal networks and construct a rotation-based constraint to couple network variance across time. We also adopt the consensus factorization as the group consistent pattern. Using two publicly available brain imaging datasets, we demonstrate that MMLC may better predict psychometric scores than some other state-of-the-art brain network representation learning algorithms. Additionally, the discovered significant brain regions are synergistic with previous literature. Our new approach may boost statistical power and sheds new light on neuroimaging network biomarkers for future psychometric prediction research by integrating longitudinal and multimodal neuroimaging data.

Appendix

We present the formulation of iterative optimization to obtain the local optimal solution. Basically, for the five learning parameters, i.e. $V^{f*}_{j}$, $V^{d*}_{j}$, $U_{i,j}, V^{f}_{i,j}$ and $V^{d}_{i,j}$, each update step learns one of them by fixing the rest. The algorithm details are described in Algorithm 1.

Fixing $V^{f*}_{j}$ and $V^{d*}_{j}$, minimize L over $U_{i,j}, V^{f}_{i,j}$ and $V^{d}_{i,j}$ Under the defined condition, objective function L only depends on $U_{i,j}, V^{f}_{i,j}, V^{d}_{i,j}$. For brevity in this subsection, we use U, V^f, V^d, V^f∗ and V^d∗ to represent $U_{i,j}, V^{f}_{i,j}, V^{d}_{i,j}, V^{f*}_{j}$ and $V^{d*}_{j}$. Therefore, the new objective function can be simplified as:

$$ \begin{array}{ll} L_{1}= & \parallel{X^{f}-U(V^{f})^{T}}{\parallel^{2}_{F}}+\alpha\parallel{X^{d}-U(V^{d})^{T}}{\parallel^{2}_{F}} \\ &+\lambda_{1}\parallel V^{f}-V^{f*}{\parallel^{2}_{F}} +\lambda_{1}\parallel V^{d}-V^{d*}{\parallel^{2}_{F}}+\lambda_{2}G. \end{array} $$

(12)

First, we further fix V^f and V^d to update U. For a given subject i and time point j, we could take the derivative of L₁ with respect to U.

$$ \begin{array}{@{}rcl@{}} \frac{\partial L_{1}}{\partial U}&=&2(U(V^{f})^{T}V^{f}-X^{f}V^{f})+2\alpha\left( U(V^{d})^{T}V^{d}\right. \\ &&\left.\left.-X^{d}V^{d}\right)\right)+\lambda_{2}G^{\prime}(U). \end{array} $$

(13)

Here, $G^{\prime }(U)$ is the derivative of U with respect to U. Given a step size l, we update U as $U_{new}=U_{pre}-l*\frac {\partial L_{1}}{\partial U_{pre}}$. Then, we fix V^d and U to update V^f. The objective function in functional network part is related to V^f, thus the gradient of L₁ with respect to V^f is:

$$ \frac{\partial L_{1}}{\partial V^{f}} = 2(V^{f}(U)^{T}U-X^{f}U) +2\lambda_{1}(V^{f}-V^{f*})+\lambda_{2}G^{\prime}(V^{f}). $$

(14)

Similarly, we update V^d with the same procedure as V^f,

$$ \frac{\partial L_{1}}{\partial V^{d}} = 2\alpha(V^{d}(U)^{T}U-X^{d}U) +2\lambda_{1}(V^{d}-V^{d*})+\lambda_{2}G^{\prime}(V^{d}). $$

(15)

Fixing $U_{i,j}, V^{f}_{i,j}$ and $V^{d}_{i,j}$, minimize L over $V^{f*}_{j}$ and $V^{d*}_{j}$ For brevity in this subsection, we use ${V^{f}_{i}}, {V^{d}_{i}}, V^{f*}$ and V^d∗ to represent $V^{f}_{i,j}, V^{d}_{i,j}, V^{f*}_{j}$ and $V^{d*}_{j}$. We observe that for each time j, the framework will generate a group-wise $V^{f*}_{j}$ and $V^{d*}_{j}$. Therefore we can reorganize the objective function L to make it only relate to those two parameters, as below:

$$ L_{2}=\lambda_{1}\sum\limits_{i=1}^{N}\parallel {V^{f}_{i}}-V^{f*}{\parallel^{2}_{F}}+\parallel {V^{d}_{i}}-V^{d*}{\parallel^{2}_{F}} +\lambda_{2}G(V^{f*},V^{d*}). $$

(16)

After updating all individual U_i, ${V^{f}_{i}}$ and ${V^{d}_{i}}$, we could take the derivative of L₂ with respect to V^f∗.

$$ \begin{array}{ll} \frac{\partial L_{2}}{\partial V^{f*}}=2\lambda_{1}{\sum}_{i=1}^{N}(V^{f*}-{V^{f}_{i}})+\lambda_{2}G^{\prime}(V^{f*}). \end{array} $$

(17)

For V^d∗, an equality constraint (V^d∗)^TMV^d∗ = I will regulate the gradient direction of L₂ with respect to V^d∗, which makes the solution difficult. Instead of directly finding an optimal direction with gradient descent on the surface described by original objective function, we construct the descent curves on the constraint-based Stiefel manifold (Hwang et al. 2016). Specifically, V^d∗ will be divided into two submatrixes $V^{d*}=[V^{d*}_{1};V^{d*}_{2}]$, where $V^{d*}_{1}\in \mathbb {R}^{s\times p}$ is the free variable to be solved and $V^{d*}_{2}\in \mathbb {R}^{(n-s)\times p}$ is the fixed variable treated as constants. Then we rearrange the constraint as:

$$ \begin{bmatrix} V^{d*}_{1}\\ V^{d*}_{2} \end{bmatrix}^{T}\begin{bmatrix} M_{11} &M_{12}\\ M_{12}^{T} &M_{22} \end{bmatrix}\begin{bmatrix} V^{d*}_{1}\\ V^{d*}_{2} \end{bmatrix}=I. $$

(18)

It is easy to conclude that M₁₁ is a full rank positive definite matrix. Then a descent curve based on the previous V^d∗ will be constructed and it starts at the point $P_{s}=V^{d*}_{1}+M_{11}^{-\frac {1}{2}}M_{12}V^{d*}_{2}$ which is the initial point for the line search on the generalized Stiefel manifold. Given the descending gradient $-L_{2}^{\prime }(P)=-\frac {\partial L_{2}}{\partial V^{d*}}\circ \frac {\partial V^{d*}}{\partial P}$ at point P, we further project $-L_{2}^{\prime }(P)$ onto the tangent space of the Stiefel manifold by constructing a skew-symmetric matrix:

$$ A=L_{2}^{\prime}(P){P_{s}^{T}}-P_{s}L_{2}^{\prime}(P)^{T}. $$

(19)

This will lead to a curve function Y (τ) by the Crank-Nicolson-like design as in the paper (Wen and Yin 2013).

$$ Y(\tau)=(I+\frac{\tau}{2}AM_{11})^{-1}(I-\frac{\tau}{2}AM_{11})P_{s}. $$

(20)

The above function gives a linear search procedure of updating point P by P_new = Y (τ) for small τ which results sufficient decrease in L₂. Finally, the next feasible $V^{d*}_{new}$ will be given as:

$$ V^{d*}_{new}(P)= \begin{bmatrix} P-M_{11}^{-\frac{1}{2}}M_{12}V^{d*}_{2}\\ V^{d*}_{2} \end{bmatrix}. $$

(21)

Optimization of the variant models. In this study, different variants of the model are studied to test how each coupling term affects the outcome of the learning tasks. Their MP parameters, e.g., dimension P, are deliberately set to be the same as the proposed model MMLC. We find relatively similar patterns of MFCSMC and MFCSLC that both of them have a close performance to MMLC with larger P values, i.e., in high dimensional feature space. As we discussed in Section 4, the proposed MMLC enjoys better computational efficiency for high dimensional feature space. Meanwhile, sMFLC yields a relatively stable pattern as EigLC, which has no significant improvements by increasing P beyond 10. It is partially because the knowledge solely comes from the structural modality containing sparse connections. For the purpose of replicating our investigation, we provide the details on how to optimize these comparison algorithm as follows. Specifically, to optimize MFCSMC, in A, we skip steps in line 11-19 and update ${V_{j}^{d}}*$ as ${V_{j}^{f}}*$ in line 9. To optimize sMFLC, we skip steps in line 5,6,9 to avoid updates of $U_{i,j}^{f}$, $V_{i,j}^{f}$ and ${V_{j}^{f}}*$. As for MFCSLC, we update $U_{i,j}^{f}$ and $U_{i,j}^{d}$ independently but keep the rest of A. For all algorithms, Their learning rates are all set as 1e − 5.