Multi-way Regression Reveals Backbone of Macaque Structural Brain Connectivity in Longitudinal Datasets

Abstract

Brain development has been an intriguing window to study the dynamic formation of brain features at a variety of scales, such as cortical convolution and white matter wiring diagram. However, recent studies only focused on a few specific fasciculus or several region. Very few studies focused on the development of macro-scale wiring diagrams due to the lack of longitudinal datasets and associated methods. In this work, we take the advantage of recently released longitudinal macaque MRI and DTI datasets and develop a novel multi-way regression method to model such datasets. By this method, we extract the backbone of structural connectome of macaque brains and study the trajectory of its development over time. Graphic statistics of these backbone connectomes demonstrate their soundness. Our findings are consistent with the reports in the literatures, suggesting the effectiveness and promise of this framework.

1 Introduction

Nowadays, connectome has been widely studied in basic brain research and clinical studies due to its importance in understanding structural and functional brain architectures from a global perspective [1, 2]. Studying the development of such connectomes could play a critical role in understanding a variety of time-dependent issues, such as the mechanisms of cortical convolution, the emergence of cognition, the possible onset of brain abnormality, and etc. [3]. However, most previous connectome studies are limited to adult brains because of the unavailability of longitudinal datasets [1, 2].

Recently, a longitudinal UNC-Wisconsin neurodevelopment rhesus MRI database is released (http://www.nitrc.org/projects/uncuw_macdevmri/), which includes T1-weighted MRI and diffusion tensor imaging (DTI) data, providing a unique opportunity for studying how brain structures of newborns develop into their infancy, juvenile and adulthood. Studies based on such longitudinal datasets could significantly complement the previous structural brain development studies in three ways: (1) many previous imaging data studies were voxel based. Voxel-wise attributes such as white matter volume growth rates [4, 5] were well studied but they provided limited clues to the development of axonal pathways; (2) axonal pathway related studies were only limited to a certain fasciculus of interest, such as the developments in Area 17&18 connections [6]. It is difficult to make up a big picture of the evolvement of an axonal network architecture, or a connectome, from these scattered profiles; (3) although dissection studies based on data modalities such as microscopic histology data and meso-scale tract-tracing data were widely available [7], it is not possible at this moment to collect longitudinal data or reconstruct a global scale brain connectome via these techniques.

Therefore, we take the advantages of the UNC-Wisconsin longitudinal macaque dataset and attempt to study the postnatal development of global wiring diagram. We construct longitudinal connective matrices for a subject at four different time points. A novel multi-way regression method is proposed to extract the ‘backbone’ connective network by using the matrix of the first time point as the baseline. This process is repeated for the other time points, such that a trajectory of ‘backbone’ connective networks are obtained. We use a variety of graphic statistics, such as node-wise strength and betweenness and efficiency of networks, to demonstrate the soundness of these ‘backbone’ networks. The effectiveness of the method also suggests its potential in future studies of brain functions and abnormality.

2 Materials and Methods

2.1 Overview

Generally, we use T1-weighted MRI, DTI data and a brain map to construct a global cortical connectome, a connective matrix, for a subject at a time point (the first three columns in Fig. 1). All subjects are divided to groups according to the time points of scans. These matrices are stretched to connective feature vectors, and these vectors within one time point compose a feature matrix (the 4^th column). The feature matrix of the t₁ is taken as the baseline, to which the other feature matrices are group-wisely regressed to produce weight vectors vs via the proposed multi-way regression method. The weight vectors can be transformed back to connective matrices (the right-most row). Because of the sparse term (see Eq. (6) in Sect. 2.3), the nonzero values in these weight matrices could imply that they are ‘preserved’ connections over time, and thus taken as ‘backbone’ connections associated to t₁. Notably, we will not differentiate the weight vectors and their corresponding matrices format in the following sections. It is also worth noting that Fig. 1 only illustrates the pipeline for t₁. The method is progressively applied to the other time points in a similar way, such that a trajectory of such ‘backbone’ connective maps are obtained.

Fig. 1.

Flow chart of the analysis pipeline in this paper. Data of the 1^st time point (t₁ for short), for example, in the green frame is used as baseline, to which data of other time points in the blue frame is regressed simultaneously via the multi-way regression method detailed in Sect. 2.3. The yellow frame highlights the results. Please refer to main text for detailed explanations.

2.2 Data Description and Preprocessing

T1-weighted MRI and DTI data in the UNC-Wisconsin neurodevelopment rhesus MRI database are used. Four time intervals are defined as four time points: 1–4 months, 5–8 months, 9–12 months and 13–16 months. Ten subjects that have data of all the four time points are selected in this work. The basic parameters for diffusion data acquisition are: resolution of 0.65 × 0.65 × 1.3 mm³, a matrix of 256 × 256 × 58, diffusion-weighting gradients applied in 120 directions and b value of 1000 s/mm², ten images without diffusion weighting. For each DTI scan, the corresponding T1 weighted MRI data in the dataset has been registered to the UNC Primate Brain Atlas space [8]. The resolution of this atlas space is 0.27 × 0.27 × 0.27 mm³ and a matrix of 300 × 350 × 250.

DTI data is defined as intra-subject standard space. T1 weighted MRI data is nonlinearly warped to the FA map via FSL-fnirt [9]. Two axonal pathway orientations for each voxel in DTI data is estimated via BedpostX in FSL 5, because it was suggested that b-values of 4000 would be required to resolve a 3 fiber orthogonal system [10]. Deterministic white matter fiber streamlines (1 × 10⁴ fiber tracts for each DTI data) are reconstructed via DSIstudio [11] thereafter. Because the selection of fiber tracking parameters is not the focus of this paper, FA and angular threshold are empirically determined as 0.1 and 60^° in that small FA value for primate brain was suggested [12].

To construct structural connective matrices from DTI data, we adopt the paxinos rhesus monkey atlas [13] as the brain map, the brain sites of which are used as nodes. Because this brain map is present in T1-weighted MRI atlas space, we firstly use FSL-fnirt [9] to register the atlas T1-weighted MRI data to the one of an individual, and further register it to the DTI space (FA map) of the individual. The brain map is transposed to the individual DTI space, accordingly. Cortico-cortical connections are of our major interest here. Thus, 152 cortical brain sites on each hemisphere are used and both ipsilateral and contralateral connections are considered. The number of fiber tracts that connect two brain sites are defined as the connective strength. Examples of such connective matrices, X and Ys, are shown in Fig. 1. The connective strength within each individual is normalized by being divided by the total fiber number. We stretch the upper triangle of a connective matrix into a vector, and define it as a connective feature. A pile of such features of ten subjects within one time point gives a feature matrix (x and ys in Fig. 1).

2.3 Multi-way Regression

We have a collection of feature matrices $(\varvec{x},\varvec{y}_{1} , \cdots ,\varvec{y}_{p} )$ where p is the number of time points. $\varvec{x}$ and $\varvec{y}_{k}$ are feature matrices, consisting of ($x^{1} , \cdots ,x^{n}$) and ($y_{k}^{1} , \cdots ,y_{k}^{n}$), where n is the number of subjects and $x^{k}$ is the feature vector converted from connective matrix (see Fig. (1)). Our objective is to use $\varvec{x}$ as a baseline and find a group of weights ($\varvec{v}_{1} , \cdots ,\varvec{v}_{p}$), to which ($\varvec{y}_{1} , \cdots ,\varvec{y}_{p} )$ are projected (Eq. (1)), respectively, such that new $\varvec{y}_{k}$ s shown below are group-wisely similar to $\varvec{x}$. We will use $\varvec{y}$ in the following equations to denote those before the projection.

$$ \varvec{y}_{k}^{new} \to \left\langle {\varvec{v}_{1} ,\varvec{y}_{k}^{old} } \right\rangle $$

(1)

The similarity is defined as correlation. Therefore, the function to be maximized is

$$ \mathop {max}\limits_{{\varvec{v}_{1} , \cdots ,\varvec{v}_{p} }} \sum\nolimits_{k = 1}^{p} {corr(\varvec{x},\varvec{y}_{k} ) + } \sum\nolimits_{k = 1}^{p} {\sum\nolimits_{l = k + 1,l \ne k }^{p} {corr(\varvec{v}_{k} ,\varvec{v}_{l} )} } $$

(2)

where the second term is to maximize the similarity among $\varvec{v}_{k}$. With this term, $\varvec{y}$ s will be regressed to $\varvec{x}$ in a group-wise manner. We define the correlation as

$$ corr(\varvec{x},\varvec{y}) = \langle \varvec{x},\varvec{y}\rangle /\left\| \varvec{x} \right\|\left\| \varvec{y} \right\| $$

(3)

and covariance matrix of (x, y) as

$$ C\left( {\varvec{x},\varvec{y}} \right) = \left[ {\begin{array}{*{20}c} {C_{xx} } & {C_{xy} } \\ {C_{yx} } & {C_{yy} } \\ \end{array} } \right] $$

(4)

where $C_{xx}$ and $C_{yy}$ are within-set covariance matrices and $C_{xy} = C_{yx}^{'}$ are between-sets covariance matrices, we can rewrite Eq. (2) as

$$ \mathop {max}\limits_{{\varvec{v}_{1} , \cdots ,\varvec{v}_{p} }} \sum\nolimits_{k = 1}^{p} {\frac{{\varvec{u}^{'} C_{{xy_{k} }} \varvec{v}_{k} }}{{\sqrt {\varvec{u}^{'} C_{xx} \varvec{u} \varvec{v}_{k}^{'} C_{{y_{k} y_{k} }} \varvec{v}_{k} } }} + } \sum\nolimits_{k = 1}^{p} {\sum\nolimits_{l = k + 1,l \ne k }^{p} {\frac{{\varvec{v}_{l}^{'} \varvec{v}_{k} }}{{\sqrt {\varvec{v}_{l}^{'} \varvec{v}_{l} \varvec{v}_{k}^{'} \varvec{v}_{k} } }}} } $$

(5)

where $\varvec{u} = 1$, because $\varvec{x}$ is used as a baseline. The corresponding Lagrangian is

$$ L\left( {\varvec{v}_{1} , \ldots ,\varvec{v}_{p} } \right) = \sum\nolimits_{k = 1}^{p} {\varvec{u}^{'} C_{{xy_{k} }} \varvec{v}_{k} - \frac{p}{2}\varvec{u}^{'} C_{xx} \varvec{u} - \frac{1}{2}} \sum\nolimits_{k = 1}^{p} {\varvec{v}_{k}^{'} C_{{y_{k} y_{k} }} \varvec{v}_{k} + } \sum\nolimits_{k = 1}^{p} {\sum\nolimits_{l = k + 1,l \ne k}^{p} {\varvec{v}_{l}^{'} } } \varvec{v}_{k} - \frac{{\left( {p - 1} \right)}}{2}\sum\nolimits_{k = 1}^{p} {\varvec{v}_{k}^{'} \varvec{v}_{k} } - \lambda \sum\nolimits_{k = 1}^{p} {\left\| {\varvec{v}_{k} } \right\|}_{1} $$

(6)

Equation (6) is our objective function. Because we have long and sparse feature vectors $\varvec{x}$ and $\varvec{y}$, $\ell_{1}$ norm, $\lambda \sum\nolimits_{k = 1}^{p} {\left\| {\varvec{v}_{k} } \right\|}_{1}$, is added to Eq. (6) to give sparse $\varvec{V}_{{k^{\text{S}} }}$

By computing derivatives in respect to $\varvec{v}_{k}$, we have

$$ \frac{\partial L}{{\partial \varvec{v}_{k} }} = C_{{y^{k} x}} \varvec{u} - pC_{{y_{k} y_{k} }} \varvec{v}_{k} + \sum\nolimits_{k = 1}^{p} {\sum\nolimits_{l = k + 1,l \ne k }^{p} {\varvec{Ev}_{l} - \left( {p - 1} \right)\varvec{Ev}_{k} = 0} } $$

(7)

$$ (p(C_{{y_{k} y_{k} }} + \varvec{E}) - \varvec{E})\varvec{v}_{k} = C_{{y_{k} x}} \varvec{u} + \sum\nolimits_{k = 1}^{p} {\sum\nolimits_{l = k + 1,l \ne k }^{p} {\varvec{Ev}_{l} } } $$

(8)

Equation (8) has the linear form of $\varvec{Ax} = \varvec{b}$, and is used to iteratively update $\varvec{v}_{k}$s. The algorithm is summarized and executed as follows:

We convert weight vectors $\varvec{v}_{k}$ to matrix format $\varvec{V}_{{k^{\text{S}} }}$ (see Fig. 1). The values of these matrices element represent their contribution to the similarity between $\varvec{Y}_{k}$ and $\varvec{X}$. Because of the sparse term introduced, nonzero elements in $\varvec{V}_{k}$ indicate the corresponding connections in $\varvec{Y}_{k}$ are preserved from $\varvec{X}$. It is worth noting that $\varvec{V}_{{k^{\text{S}} }}$ could be similar to each other due to the constraint ${\text{c}}orr(\varvec{v}_{k} ,\varvec{v}_{l} )$ in the objective function Eq. (2). Therefore, only the $\varvec{V}_{1}$ corresponding to $\varvec{y}_{1}$ is used as a representative (highlighted by an orange frame in Fig. 1) to interpret the results, defined as $\varvec{V}_{t1}$. Nonzero elements in $\varvec{V}_{t1}$ could imply that the corresponding connections are preserved over all time points.

We execute this method in a progressive manner. For the first run, $\varvec{x}$ is the feature matrix of t₁. $(\varvec{y}_{1} ,\varvec{y}_{2} ,\varvec{y}_{3} )$ correspond to those from t₂ to t₃. Unit vectors are used to initialize $\varvec{V}_{{k^{\text{S}} }}$ for this run. For the second run, $\varvec{x}$ is the feature matrix of t₂. $(\varvec{y}_{1} ,\varvec{y}_{2} )$ correspond to those from t₃ to t₄. To initialize $\varvec{V}_{{k^{\text{S}} }}$, we replace an element of unit vectors with 0 if its counterpart in $\varvec{V}_{t1}$ from the previous run is nonzero, such that connections not preserved (or new to t₁) are considered with higher priority. Following this scheme, we obtain a series of representative $(\varvec{V}_{t1} ,\varvec{V}_{t2} ,\varvec{V}_{t3} )$ for t₁ ~ t₃.

2.4 Graphic Statistics

In order to further explore the property and evaluate the validity of the extracted connective maps $\varvec{V}_{\text{S}}$, we adopt a variety of graphic statistics at different scales, such as node-wise measurements (strength and betweenness) and connectome-wise one (network efficacy). For node-wise measurements, we define an overlap ratio to compare two connective matrices. By taking node strength for example, we firstly compute the strength for each node of matrices X and Y, and obtain two strength vectors $s_{X}$ and $s_{Y}$. They are then sorted in a descending order, denoted by $s_{X}^{'}$ and $s_{Y}^{'}$. By taking X as a baseline matrix, we generate an overlap curve, the k^th value of which is defined as the length of intersection of the first k elements of $s_{X}^{'}$ and $s_{Y}^{'}$ divided by the length of $s_{X}^{'}$. The overlap ratio is defined as the area under this overlap curve.

3 Results

3.1 Effectiveness of the Multi-way Regression Method

The $\lambda$ in the objective function (Eq. (6)) regulates the sparsity of weight vectors $\varvec{V}_{{k^{\text{S}} }}$ and is optimized via a five-fold cross-validation scheme such that the objective function is maximized. To demonstrate the stability of the method on the used dataset, we divide the ten subjects into two groups (five subjects in each one with no overlap) and separately execute the algorithms on them. 250 permutations of such divisions are used. On average, the similarity (Pearson correlation) of the resulting weight vectors between two groups is 0.70 ± 0.34, demonstrating the robustness of the methods on this dataset. It is worth noting that the resulting weight vectors on five subjects are definitely different from those on all ten subjects. This robustness experiment implies that our method can produce reproducible results on ten subjects. Therefore, the following analyses are all based on ten subjects.

3.2 Development of Connective Connectome of Infant Macaque Brains

Figure 2(a) shows the weight matrices on the left-most column. They are obtained by using connective matrices of t₁, t₂ and t₃ as baselines, respectively. Connections exclusive to a certain time point are highlighted with different colors by comparisons among them. Backbone connections are represented by nonzero elements in them. They are thus used as masks to screen the backbone connections from the original connective matrices other than the baseline ones. For example, $\varvec{V}_{t1}$ associated to t₁ is used to mask the original connective matrices from t₂ to t₄ (the second row). We use the Pearson correlation coefficients to measure the similarity among the original connective matrices. They are compared to those obtained from the masked counterparts. Taking the second row for example, the Pearson correlation coefficients are computed between the masked matrices and the original t₁ matrix in the first row (Table 1). These results suggest that the similarity between the masked connective matrices and the baseline ones are significantly raised regarding to small p-values via t-test, suggesting the soundness of the extracted backbone connectomes. The soundness of these backbone connectomes are also demonstrated by reports in the literatures. For example, cyan connections which are associated to and exclusive to t₃ strongly originate or terminate in frontal/temporal/parietal cortices (Fig. 2(b)) but not in occipital cortex. These results are consistent with previous reports that occipital cortex and white matters associated with it stop growing at this developmental stage (9 months–12 months) while the other cortices and the white matters associated with them have not reached their growth rate peaks yet [4, 5].

Fig. 2.

(a) The three weight matrices $\varvec{V}_{\text{S}}$ associated to t₁ ~ t₃ are on the left column. Elements exclusive to a certain time point are highlighted by different colors. Original connective matrices of four time points from one example subject are compared to those masked by the weight matrices in the left-most column. Interpretations are referred to the text. (b) Columns and rows of $\varvec{V}_{{\varvec{t}3}}$ in (a) are sorted in lobe order. White solid box separates ipsilateral and contralateral connections. White dashed boxes highlight groups of connections exclusive to t₃.

Table 1. Evaluate the soundness of backbone connectomes by means of comparison of Pearson correlation coefficients between connective matrices. P-values of t-test are reported which are conducted on the corresponding cells on the two panels. O: original; M: masked.

To further demonstrate the validity of those backbone connectomes, a variety of graphic statistics in Sect. 2.4 are measured. We compute node strength/betweenness and network efficiency for all connective matrices in Fig. 2. The node-wise measurement overlap ratios defined in Sect. 2.4 between the masked connective matrices and their original counterparts as well as the global efficiency ratios are reported in Table 2. The network efficiency values obtained from the masked matrices are normalized by being divided by those obtained from their original counterparts. Taking the 2^nd row in Table 2 (the 2^nd row in Fig. 2) for example, we find that all measurements of the masked connective matrices from t₂ to t₄ remain at the same level, suggesting that the backbone connectomes are stable over time. This observation applies to the other rows as well, though these values increase after more connections are incorporated. Altogether, these results suggest that the development of connectome might prefer progressively making changes on the originally blank locations to strengthening the developed connections, which is one of the major neuroscientific insights achieved by this study.

Table 2. Strength & betweenness overlap ratios and network efficiency ratios of the masked connective matrices.

4 Discussions and Conclusions

In this work, we proposed a novel multi-way regression method to extract developmental backbone connectomes on the longitudinal MRI and DTI macaque datasets. The effectiveness of those backbone connectomes is demonstrated by reproducibility studies, graphic statistics and previous reports. Studies on the trajectory of backbone connectomes could help reveal a possible structural wiring development mechanism, which would be further validated by other fine-grained data modalities and applied to abnormal brain datasets in the future. The objective of this study is in line with other developmental –omics studies, such as longitudinal functional connectome and genomics datasets recently released by the Allen Institute for Brain Science. An integrative analysis on those longitudinal –omics studies is also of our interest in the future.