Understanding the Relationship Between Human Brain Structure and Function by Predicting the Structural Connectivity From Functional Connectivity

Over the past decade, a growing number of studies have investigated the relationship between the structure and function of human brain by predicting the resting-state functional connectivity (rsFC) from structural connectivity (SC). Yet how the whole-brain patterns of FC emerge from SC still remains incompletely understood. Unlike previous studies, here we propose an alternative approach for addressing this issue by predicting SC from rsFC. We first hypothesize that the functional couplings among brain areas at rest are shaped at least in three phases temporally: the initial direct interplay between brain areas, the communications within and between network modules, and followed by the indirect interactions ascribed to indirect structural pathways. We then introduce a network deconvolution (ND) algorithm inspired from the mechanism of cell differentiation, named CDA, to distinguish the direct dependencies from the functional network followed by a weight trimming algorithm based on Euclidean distance kernel function for shrinking the modular effects. Finally, we keep those region pairs with shorter shortest path length (SPL) together with shorter Euclidean distance as the structural connections. We apply the model and the algorithms to three intensively studied group averaged empirical connectome datasets with different parcellation resolutions and the results demonstrate that the predicted intrahemispheric structural connections and the weights distribution are highly consistent with the empirical SC derived from diffusion magnetic resonance imaging (dMRI) and probabilistic tractography, thus strongly supporting the model and algorithms proposed.


I. INTRODUCTION
Unravelling the relationship between the relatively static anatomical topology and the diverse functional repertoire of the human brain is highly crucial for understanding the mechanisms of how brain cognition and diseases are developed [1]- [3]. Since the beginning of 1990s, diffusion MRI The associate editor coordinating the review of this manuscript and approving it for publication was Haipeng Yao . and functional MRI (fMRI) have become two emerging types of in vivo noninvasive neuroimaging means for mapping the connectivity of human brain [4], [5]. Diffusion MRI especially diffusion tensor imaging (DTI) and diffusion spectrum imaging (DSI) are able to measure the number and density of the anatomical white matter fiber tracts between brain areas while the functional MRI is capable of recording the blood oxygenation level-dependent (BOLD) signals that indirectly reflect the true neural activities occurring in human brain.
Both imaging modalities can be used to estimate one kind of brain connectome rendered as a square symmetric connectivity matrix [6]- [8]. The matrix inferred by dMRI is called structural connectivit (SC) in which each entry refers to the anatomical connection strength or density measured between two brain regions, while the matrix estimated by fMRI is termed functional connectivity (FC) in which each element is usually obtained by computing the statistical correlation between BOLD time series.
In the past decade, the relationship between SC and FC has received increasing attention in the field of neuroscience. Numerous statistical and computational modeling studies have been conducted to uncover how SC gives rise to FC and how FC emerges from SC by predicting FC from SC. Some statistical studies suggested that FC can be predicted but to a limited extent by some anatomical connectivity features such as fiber length, fiber counts or physical distance between brain regions [9]- [13]. For example, a study using SC to predict FC [9] observed that the relationship between SC and FC is robust when structural connections are present. However, when direct SC is absent, the FC varies over a wide range and could be accounted by indirect linkages which were shown a significant predictor of FC, especially for the visual cortices of each hemisphere. Another study on macaque also supported this finding that the generation of FC in the absence of direct SC can be determined by signal flows via connections with common afferents mediated by a third area [10]. The configuration of motifs, such as length2-SC indirect patterns, revealed an important role of SC-FC relation. Note that the brain functional and structural networks also have topological properties, such as modularity. Intramodular connections, which predominate in highly modular brain networks, are generally short distance. However, between-modular connections are generally long-distance, even though they can help achieve high global efficiency of brain networks. Penalizing SC connections by a function of distance was demonstrated efficient to simulate FC [12]. Others explored the relations using network topological features such as shortest path, search information, path transitivity [14], [15], as well as unweighted degree product [16], [17]. While computational modeling approaches such as linear and nonlinear neural mass models (NMM) tend to simulate the neural dynamics among neuronal populations in human brain [3], [9], [18]- [21], in which hundreds of differential equations and tens of physiological parameters are involved, making it a daunting task to find an optimal solution. More recent studies relate SC and FC through machine learning approaches such as matrix or spectral mapping [22], [23], connectome embedding [24], temporal multiple kernel learning [25], and routing with linear programming [26], etc. The key difference between machine learning based approaches and traditional methods is that the machine learning based approaches are able to highly capture the relation between SC and FC from the connectome data directly, yet the model parameters need to be retrained for different connectome data. In fact, the fundamental role of how the whole-brain pattern of FC particularly FC between region pairs without direct anatomical link is shaped still remains fully unknown.
Different from previous studies, here we present an opposite approach to address this issue by predicting the SC from rsFC. We first establish a FC model that allows for three different types of effects including the direct dependencies between neuronal populations, information exchange within and between network modules, and the indirect interactions via indirect pathways. We then build network inversion algorithms to remove the indirect and modular effects from FC in a reverse order and finally predict the structural connections from the remaining connections with shorter shortest path length and Euclidean distance. Finally, we test the model and algorithms on three group averaged empirical connectome datasets with different resolutions and the performances are evaluated by comparing the predicted structural connections with the empirical structural connections in terms of the number of the correctly or falsely predicted links as well as the weight distributions of the connections.

A. CONNECTOME DATASETS
Three connectome datasets with different resolutions were used for testing and assessing the proposed models and algorithms.
The 90 ROIs (regions of interest) dataset was obtained at Weill Cornell Medicine and employed in [18], [19], in which the structural and diffusion MRI together with the restingstate fMRI data were collected on 8 healthy adults and parcellated into 90 cerebral regions after diffusion tractography processing. The structural connection weight between any two regions was estimated by a weighted sum of fiber tracts going between them and corrected by topological distance. The functional connectome and the corresponding matrices were obtained by computing the Pearson correlation between BOLD time series derived from resting-state fMRI. Both resting-state FC and SC matrices were averaged across the 8 individual participants. To be in common with the other two datasets, we rearranged the lobe order as: frontal, parietal, occipital, temporal, and subcortical.
The 246 ROIs dataset is a subset of the connectivitybased brain imaging research database (C-BIRD) at Beijing Normal University [27], which contains multi-modal MRI data from 49 young healthy subjects. Informed consent to all participants in the Institutional Review Board (IRB) of the State Key Laboratory of Cognitive Neuroscience and Learning of Beijing Normal University and approved study was obtained. All participants agreed to share data freely on the Internet in an anonymous form. In our study, the human Brainnetome atlases of 246 brain ROIs was used (http://atlas.brainnetome.org) [28]. This atlas composed of 210 cortical and 36 sub-cortical subregions, which contains the information on both anatomy and functional connections. Probabilistic tracking was performed and the fiber direction was determined according to the probability to obtain the structural connection of 246 ROIs. The structural connectivity matrix of each participant was obtained through PANDA software (https://www.nitrc.org/projects/panda/), in which rows and columns represent brain nodes, and element values represent connection probabilities between nodes. The resting-state FC matrix was calculated using pair wise Pearson's correlation coefficients of BOLD time series obtained from each brain area. Both rsFC and SC matrices were averaged across the 49 individual participants.
The high resolution dataset with 998 ROIs has been extensively studied in previous research [9], [14], [29], [30], in which the structural and diffusion MRI together with the resting-state fMRI data were collected on 5 healthy right-handed male subjects. DSI was performed using a diffusion-weighted single-shot echo planar imaging (EPI) sequence. Following diffusion spectrum and T1-weighted MRI acquisitions, the segmented gray matter was partitioned into 998 ROIs. White-matter tractography was performed with a custom streamline algorithm and fiber connectivity was aggregated across all voxels within each of the 998 predefined ROIs. The resting-state FC matrix was calculated using pair wise Pearson's correlation coefficients of BOLD time series obtained from each brain area. Both rsFC and SC were averaged across the 5 individual participants. Refer to [14] and [29] for more details.
All the aforementioned connectivity matrices are symmetric and arranged that the upper left and lower right quadrants map the right hemisphere and left hemisphere of the brain respectively, while the off-diagonal quadrants map the interhemisphere between the two.
Because weak and non-significant links may represent spurious connections in FC [31], therefore, we removed all the connections whose strength were smaller than 0.1 and the negative connections (the functional anticorrelations still remain elusive) in FC. Furthermore, all self-connections (diagonal elements in the FC matrix) were also excluded.
The three empirical group-averaged connectomes are demonstrated in Fig.1.

B. MODEL OF FC
We start with the assumption that the observed resting-brain functional network, G obs , can be formulated as the sum of at least three parts in temporal order of formation: the direct network G dir , the modular network G mdl , along with the indirect network G ind : where G obs , G dir , G mdl , and G ind are four M × M symmetric matrices, containing non-negative elements only. M indicates the number of brain areas. The direct network, G dir , stemming from the relatively invariant anatomical structure, which reflects the initial interplay among brain areas or neuronal populations. The modular network, G mdl , arising from the communication within and between network modules, which allows for local segregation and global integration. While the indirect network, G ind , induced by the indirect paths between brain areas, which can facilitate the signal transmission along structural pathways. Define G dirm = G dir +G mdl as the whole direct effects after the modular effects set up. According to previous studies [9], [14], [22], [23], the indirect effects on FC cannot be exclusively attributed to the indirect paths of length 2, but there is no denying that longer paths contribute less to the whole FC. Here, for computational convenience, we suppose that the indirect effects induced by the longer indirect paths are approximately proportional to the direct effects, thus we can model the indirect effects as where 0 < α < 1, indicating different contributions to the whole indirect effects by indirect paths of length 2 and longer paths.

C. NETWORK DECONVOLUTION WITH CDA
Distinguishing the direct network G dirm from the observed network G obs in (3) can be formulated as a network deconvolution problem [32]. Obviously, there are no real and analytical solutions to (3) due to its nonlinearity and high dimensionality. Here, we formulate the network deconvolution problem as finding an optimal direct connection matrix by optimization. Specifically, for a given direct connection matrix G dirm , we can obtain a prediction of the observed correlation matrix from (3), denoted by G p obs . Then network deconvolution process can be converted into searching for an optimal sparse direct connection matrix G dirm , subject to arg min The Frobenius 2norm term in (4) measures how well the direct connection matrix describes the observed connection matrix, and here we define this as the sum of the square of all the entries in the error matrix G obs − G p obs : whereas the 1norm term in (4) is used to expedite the iteration process during optimization and defined by the sum of the absolute value of all the entries in the error matrix G obs − G p obs : It is not easy to find an global optimal direct connection matrix G dirm from (3) and (4) using the traditional gradient descent algorithms due to the high dimensionality of the problem and will get even harder when the size of the network increases. While bio-inspired optimization algorithms, such as genetic algorithms (GAs) [33], particle swarm optimization (PSO) [34], and ant colony optimization (ACO) [35], etc., can only handle problems with dozens of dimensions. It is worth noting that Zhong et al. [36] incorporated the multi-agent concept into GAs and proposed a multiagent genetic algorithm (MAGA) which can handle high dimensional function optimization problems with 20-1000 dimensions. However, the performance of MAGA worsens rapidly when the dimensions of the solution space approaches 10000. Normally, MAGA needs tens of thousands of iterations to get the optimal solution even with 1000 dimensional functions.
Studies in developmental biology show that cell differentiation is a fundamental and central process for shaping tissues and organs with specific functions in our body. The resulting morphogenesis involves several cellular behaviors such as division, differential growth, migration, fusion, differential adhesion, contraction, as well as apoptosis (programming death) [37]. These behaviors allow for local cell-cell interactions and gene regulations occurring under precise spatiotemporal coordination [38], [39]. Inspired by the developmental mechanisms of cell differentiation at the microscale, in our previous work, we developed a new biologically inspired optimization algorithm for handling super-high dimensional optimization problems [40], in which only four behaviors were simply modeled and the algorithm can only be applied to numerical function optimization. In this paper, we redefine all the cellular behaviors mathematically and extend the algorithm to matrix optimization cases because a matrix can be simply symbolized as a cell. Additionally, the elements in the matrix can be modeled and expressed by genes and the matrix operations can be viewed as the morphological changes of a cell, thus the matrix optimization process can be realized by performing the behaviors of cell differentiation. A number of child cells descended by their parent interact with each other by regulating their internal gene expressions to evolve generation by generation by way of exhibiting different differentiation behaviors aforementioned. On the basis of the principle of survival of the fittest, cells with higher activity values will be more likely to continue to be alive after each differentiation behavior, while those with smaller activity values will die and be replaced by other more robust cells. Therefore, the quality of the cells will be better and better. After a certain number of evolving generations, the cell swarm will be smoothly differentiated towards the best shape representing a global optimal solution to the problem.
Consider C is a set of N cells, C = {cell 1 , cell 2 , · · · , cell N }, and define C s as the stem cell of C. Where cell i (i = 1, 2, . . . , N ) is the ith cell denoted by a n × n matrix and cell i−1 , cell i+1 are its two neighbors. Let P mi , P ad , P fu and P ap , denote the probability of cell migration, differential adhesion or contraction, fusion, and apoptosis during cell development and differentiation, respectively; act i and age i refer to the activity value and age of cell i, respectively.

1) CELL DIVISION AND GROWTH
Cell division is the first developmental behavior for nearly all types of cells, particularly asymmetric mitosis, which occurs autonomously and can be modeled by dividing each cell asymmetrically into descendent cells distinct from their parents [37].
Specifically, assume k is a random number within (1, n), the division of a cell can be defined as dividing the mother cell into two daughter cells: where the first k rows of cell i are duplicated into the corresponding rows of cell i1 , the remaining n − k rows of cell i are copied into the corresponding rows of cell i2 , while the rest of the rows in cell i1 and cell i2 are all set to be 0. Then make the two daughter cells grow by replacing the last n − k rows of cell i1 and the first k rows of cell i2 with random values generated from the stem cell, C s , respectively, i.e. u(0, C s ), note that u(·, ·) is a uniform random number generator and u(0, C s ) represents a random matrix.
For computational consideration, if the activity value of either of the daughter cells is higher than their parent, the parent cell will be replaced by the daughter cells with higher activity value. Otherwise, the parent cell will remain to live. VOLUME 8, 2020 2) CELL MIGRATION Cell migration means the cells inside the tissue will move from one location to another. Migration can be random or towards a preferred direction. Here, we use two strategies for modeling cell migration in CDA: directed migration and random migration. Directed migration means the cell will move towards where the best cell (with the highest activity value currently) in the current generation is located, while random migration means the cell moves to some location between its current location and the best cell randomly.
Specifically, for each cell i , if u(0, 1) > P mi , the cell will perform directed migration strategy: Otherwise, the cell will perform random migration strategy: where cell best is the best cell in the current generation, Mstep i stands for the migration step length which is defined by, where act i stands for the activity value of cell i , which can be defined by the objective function in (4).

3) CELL FUSION
Cell fusion indicates the cell will merge with the best cell around its neighbors. There are also two strategies of cell fusion, one is dominated by the cell per se, the other by the best cell from its neighbors.
For each cell i , if u(0, 1) < P fu , the cell will perform the following fusion strategy: Otherwise, the cell will perform strategy: where cell best is the cell with maximum activity value among its 4 neighbors:

4) CELL DIFFERENTIAL ADHESION AND CONTRACTION
Cell differential adhesion and contraction are crucial steps in cell differentiation and defined by the gene networks. Cell differential adhesion and contraction change the size and shape of cells and then generate long-range forces between the differentiated cells, and results in cell migration and mechanical stresses on cells, which may trigger cell division, growth, and death [39]. In CDA, cell differential adhesion and contraction are modeled to change the values of each gene by increasing or decreasing a small random value related to the stem cell through heating or cooling the cell by changing the temperature T .
For each cell i , if u(0, 1) < P ad , the cell will perform the following adhesion and contraction strategy: (13) where T denotes the temperature, which will increase by 1 degree after each iteration, if it exceeds a predefined value it will be reset to zero. At present, the predefined temperature is set to be 100 degrees Celsius for CDA.

5) CELL APOPTOSIS
Cell apoptosis is also crucial in pattern formation and morphogenesis, which can transform one pattern into another.
In CDA, we model cell apoptosis by just resetting some genes in the cell to be zero or replacing them by a small value generated randomly if the activity values of these genes are still below some value after a certain number of generations. The life span of a cell, 100, is chosen in the current CDA. Specifically, for each cell i , if its age exceeds the life span and u(0, 1) > P ap , those genes with smaller activity values will be set to zero; or if u(0, 1) < P ap , these genes will be randomly generated using u(0, 0.1).
Based on the above description, the network deconvolution algorithm based on CDA can be implemented by the following steps: Step 1: Initialize the cell swarm. Generally, in an evolutionary algorithm, the initial population is usually created by randomly assigning each gene to a binary or a real number. While for CDA-ND, all the N initial cells are descended from the stem cell, denoted by the observed functional network, G obs , i.e.: It is worth noting that each cell represents a candidate solution and can be specified by a M × M matrix.
Step 2: Estimate the activity value of each cell, act(cell i ), based on (4); Step 3: For each cell, exhibit the division behavior first, then followed by growth, migration, fusion, adhesion/ contraction, and apoptosis differentiation behaviors as well. And then evaluate all the cells after each behavior by comparing the activity values between the mother cells and the child cells and keep the cells with higher activity values; Step 4: Update the cell swarm and find the one with the highest activity value, i.e., the best cell in the current generation t, cell t best ; Step 5: If act(cell t best ) > act(cell t−1 best ), leave the cell t best in the next generation; otherwise cell t best ← cell t−1 best . This step will ensure the CDA to converge to the global optimum because the best solution is always maintained in the swarm; Step 6: Decide whether it satisfies the termination criterion or not. If the stopping criterion holds, the algorithm will end, otherwise repeat from Step3.
The overview of CDA-ND is schematically illustrated in Fig.2.

D. KERNEL-BASED WEIGHT TRIMMING (KWT)
The deconvolved direct correlation network G dirm after network deconvolution with CDA consists of two parts: the initial direct dependencies between brain areas and the modular effects owing to the modular network structure. Studies on network topology and properties of the anatomical connectivity of human brain indicated that the brain structure shows small-world characteristics, i.e., with high local clustering coefficient and small average shortest path length, allowing for local functional segregation and global functional integration [2], [12], [29]. This can be evidenced by the rich-club organization of the structural connectivity of the human brain [41], in which many small group of modules clustered by short-range connections between cortical region pairs are highly connected through hubs. The smallworld and modular structure of human brain network suggest that there more likely exists an anatomical link between region pairs with shorter physical distance while the probability of having a link between two distant regions is much smaller. Thus, we can construct a Gaussian kernel matrix based on Euclidean distance between brain areas and map it onto the deconvolved direct network G dirm after CDA to hold the short-range and weaken the long-range connection strength.
Specifically, for each region pair (i, j), we choose the Gaussian kernel function as, where D(i, j) denotes the Euclidean distance between region i and j, which can be estimated using the mean Talairach coordinates of voxels comprising an ROI [9]; σ is a volume parameter, controlling the size of the modules and varying with the network size. Finally, the resulting direct network can be obtained by kernel-based weight trimming (KWT), as follows,

III. RESULTS
Our results are verified by following: first, evaluating the optimization algorithm-CDA by minimizing the error between observed matrix (i.e., functional networks) and simulated observed matrix to find the optimal deconvolved direct matrix on example datasets; second, validating the role of the shortest path length (SPL) to reconstruct empirical SC by showing the receiver-operating characteristic curves (ROC) and the area under curves (AUC); third, verifying the predicted SC by checking both the rate of correctly predicted connections and the similarity of connection strength with the empirical SC; fourth, comparing proposed algorithm with two existing methods by showing the prediction rates. For convenience, we list the important acronyms used in Table 1.

A. EVALUATION OF CDA
Before applying the proposed CDA to the empirical connectome datasets. We first apply it to search for an optimal symmetric matrix with 100×100 dimensions, which satisfies (3), where α = 0.6. The simulated direct matrix and the observed matrix are shown in Fig. 3 (A) and (B), respectively. Some model parameters need to be predetermined before the CDA works. The population size of cells, Nc, can usually be chosen as 20-50. The migration probability P mi determines whether CDA exploits new locations between the current cell and the current best cell or explores around the current global best cell, when P mi < 0.5, CDA mainly acts on searching new solutions (exploitation). The fusion probability P fu decides how often the cells exchange their information with their local best (within their neighbors). While the adhesion or contraction probability P ad ensures the diversity of cells, obviating the premature of the algorithm and the apoptosis probability After a large number of tests, in our study, the parameters of CDA are chosen as: N = 20, P mi = 0.5, P fu = 0.5, P ad = 0.2, P ap = 0.1, the life span of a cell is 100, the temperature T is 100. Two termination criteria are exploited: (1) to ensure each entry in the error matrix,, G obs − G p obs , is less than 10 −5 , then the total error calculated by (4) satisfies ε ≤ (M * (M − 1) /2) × 10 −5 , where M is the order of the matrix to be optimized; (2) the maximal number of generations exceeds a predefined number.
The optimal direct matrix after CDA is shown in Fig.3 (C) and the convergence performance is shown in Fig.3 (D). It can be seen that the deconvolution error drops down to the predefined precision (0.0495) from 14.517 after 55,597 numbers of generations ( Fig.3 (D), about 1.634 hours on a Dell laptop computer with Intel (R) Core (TM) i7-7560 CPU @ 2.40GHz, 16G memory, and Windows10 operating system). The Pearson correlation between the simulated direct matrix and the deconvolved direct matrix is 0.9987, indicating that the simulated direct network with size 100 × 100 can be fully distinguished from the observed functional network with the proposed CDA.
We then apply the CDA to the empirical 90-ROI and 998-ROI datasets, respectively. The recordings of the convergence curves for the 45 × 45 matrix optimization (RH, 90-ROI dataset) and 500 × 500 matrix optimization (RH, 998-ROI dataset) are shown in Fig. 4 (A) and (B), respectively. It can be seen that the deconvolution error drops down to 0.0017 from 5.77 after 10,000 numbers of iterations for the 45 × 45 matrix optimization (Fig. 4 (A), about 3 minutes on the same Dell laptop computer), while from 258.38 to 5.76 after 20,000 numbers of iterations for the 500 × 500 matrix optimization (Fig. 4 (B), about 10 hours on the same  Dell laptop computer ), showing the proposed CDA is able to find an optimal matrix with size up to 500 × 500, i.e. 250,000 dimensions, after a certain number of iterations.

B. STRUCTURAL CONNECTION IDENTIFICATION
Note that G dir represents the positive couplings between not only the structurally connected brain regions but also some regions that are indirectly structurally coupled, which could be simulated by a neural mass model.Nevertheless, it is an intractable issue to infer the true SC by inversion from a neural mass model. Here we adopt an alternative way to build the bridge between G dir and the actual SC.
We find that the SPL together with the number of steps along the shortest path between brain areas can highly capture the relation between both connected and unconnected pairs of regions, including structural links showing negative correlations. We validate this by extracting the SPL matrices from the empirical SC matrices of the three datasets and keeping those region pairs with shorter weighted path length together with fewer steps along each shortest path. Log-weighted length was chosen due to the log-normal distribution of the connection strengths [9], [29], [42]. We evaluate the reconstructed results using the ROC curves (Fig.5), in which all the AUC value of the three datasets are close to 1, meaning nearly all the intrahemispheric structural connections, including those connections with negative functional correlations, can be correctly reconstructed when detecting the same number of connections as the actual number of connections within each hemisphere of the empirical SC. The reconstructed results are assessed by true positive rate (TPR) and false positive rate (FPR), as listed in Table 2, showing that 93.26%, 84.45%, and 91.7% of the mean intrahemispheric links are correctly reconstructed for the 90-ROI, 246-ROI, and 998-ROI datasets respectively.

C. SC PREDICTION
In what follows, we apply the model and algorithms to the three empirical datasets to predict the structural links of the 209932 VOLUME 8, 2020 right hemisphere (RH) and the left hemisphere (LH), respectively. We do not evaluate the interhemisphere (IH) because studies on dMRI suggested that the structural connections of the IH are noisy and some structural connections might be missed or falsely measured [43]- [47]. As the predicted SC is weighted, we apply two methods to verify the quality of the predicted SC. Firstly, we check the rate of correctly predicted links (prediction rate) by detecting the same number of connections as the actual number of connections within each hemisphere of the empirical SC and plotting ROC curves under a range of continuous predicted SC thresholds. The empirical structural network is binarized where connections are either absent (strength is 0) or present (strength is positive) before ROC analysis. A high prediction rate indicates that the structural links are well distinguished from a dense connected functional network regardless of the connection strength. Secondly, to assess the connection strengths of the predicted SC, we also graphically demonstrate the predicted connectomes as compared with the empirical weighted SC for each dataset, and the Pearson correlations which linearly depict the strength similarity between the two are given as well. It is worth noting that, in the current model, the optimal FC model parameter α was determined as 0.6 after extensive tests across all the three datasets, while the optimal kernel parameter σ was chosen as 30 for the 90-ROI dataset and 15 for the other two datasets with higher resolutions, respectively.
The prediction results are averaged with 5 runs and compared among CDA alone, CDA followed by SPL (CDA -SPL), as well as CDA followed by KWT and then SPL (CDA-KWT-SPL) ( Table 3). It can be seen that, compared with CDA alone, the prediction rate of the RH and LH are increased by 4.34% and 3.98% for the 90-ROI dataset, 2.88% and 3.35% for the 246-ROI dataset, and 6.44% and 7.22% for the 998-ROI dataset with CDA-SPL, respectively. While with CDA-KWT -SPL, the prediction rate of the RH and LH are raised by 18.79% and 13.58% for the 90-ROI dataset, 7.09% and 6.7% for the 246-ROI dataset, and 14.83% and 16.29% for the 998-ROI dataset respectively as compared with CDA alone. The mean SC prediction rates are more than 91% for the 90-ROI dataset, 83% for the 246-ROI dataset, and 74% for the 998-ROI dataset, respectively. These results suggest that the indirect paths and the Euclidean distance between brain areas are two contributing factors and play different roles in yielding the whole brain FC at rest, while the shortest paths with shorter lengths reflect the key features of SC.
Furthermore, Fig.6 demonstrates the predicted intrahemispheric SC connectomes in contrast to the corresponding empirical connectomes for all the three datasets.

D. COMPARED WITH EIGEN-DECOMPOSITION BASED (ED-ND) ALGORITHM
We further evaluate the performance of the proposed CDA-ND in contrast to the ED-ND algorithm [30], [32] as well as the FC thresholding (FC-TH) method [9], [11] in terms of ROC curves, as illustrated in Fig.7. The results show that the ED-ND approach outperforms the FC-TH method for the 90-ROI dataset, but both ED-ND and FC-TH methods have the similar performance for the 246-ROI and 998-ROI datasets. Whereas the performance of the proposed CDA-ND approach far exceeds the other two methods across all the three datasets.

A. SUMMARY OF THE MAIN FINDINGS
This study aims to explore how the long-time averages of the whole brain functional connectivity emerges from the relatively invariant structural topology in the absence of any external stimuli, which intrigues a growing num-VOLUME 8, 2020 ber of researchers in neuroscience during the past decade. Traditional approaches usually laid emphasize on predicting FC from SC using various of computational models. However, none of these models are able to completely cast light on the nature of the phenomenon. There is still much unknown for how the correlations emerged between regions that are not directly structurally coupled. The key distinction between our research and previous studies is that our study shows how one might go in the other direction as opposed to the prediction of FC, to uncover the structural connections from the resting-state fMRI correlations between brain regions.
We first put forward a FC model in which the observed whole resting-brain FC is formulated as the weighted combination of direct dependencies and indirect dependencies majorly owing to the indirect paths of length 2. Then we introduce a bio-inspired network deconvolution algorithm as well as a kernel-based weight trimming algorithm to search for the best structural connections that produce the observed fMRI correlations by minimizing the wiring strength and the shortest path length.
The major finding of our study is that the predicted structural connection strengths are well fit by lognormal distributions (Fig.8). The predicted connection strength ranges from 3.5476e-9 to 0.4635 for the 90-ROI dataset, 3.4087e-9 to 0.3438 for the 246-ROI dataset, and 3.1997e-5 to 0.3625 for the 998-ROI dataset, respectively, showing that the connection strengths span several orders of magnitude. Furthermore, the mean Pearson correlation between the predicted intrahemispheric SC and the corresponding empirical SC are 0.9048 for the 90-ROI dataset, 0.8735 for the 246-ROI dataset, and 0.7388 for the 998-ROI dataset, respectively. These results are highly consistent with the empirical SC derived with dMRI and the probabilistic tractography [9], [14], [17], [29], [45], [46]. The minimized wiring cost together with the maximized efficiency (shortest path length or shorter Euclidean distance between structurally connected node pairs) provide strong accounts for supporting the view that the human brain network is organized in a tradeoff between wire cost and efficiency [48].
After CDA-ND, we found that more than 62% of the mean intrahemispheric FC was caused by the indirect effects for the 90-ROI dataset, 64% for the 246-ROI dataset, and 73% for the 998-ROI dataset, respectively. While after KWT, we found that the modular effects remain relatively steady, accounting for ∼24% of the whole FC for the 90-ROI dataset, ∼29% for the 246-ROI dataset, and ∼20% for the 998-ROI dataset, respectively. Finally, our findings indicated that only a small fraction of FC directly reflects SC, accounting for ∼14% of the mean intrahemispheric FC for the 90-ROI dataset, ∼7% for the 246-ROI and 998-ROI datasets, respectively, indicating that the higher the resolution of parcellation is, the smaller the direct effects influenced by SC.
The main methodological contributions of our work are the integration of indirect effect and modular effect in the contributing factors of shaping FC to uncover SC; and the utilization of new method CDA to search optimal results, leading to several benefits. First of all, the indirect effect is mainly attributed to the length2-paths, therefore, providing a simple representation in the model of FC which reduces computing and time complexity. Besides, SPL, the minimum path distance, is used to remove the indirect effect further, providing new insight and methodology to uncover SC. The modular effect, described by physical distance, is trimmed via a Gaussian kernel, suggesting a possible relationship between network topology and functional network. The integration of above effects not only helps to shed light on how much they contribute to generating FC but can offer a way of comparison between the two by removing effects one by one. The second benefit of our method is the possibility to apply the CDA to other optimization studies, especially for the high-dimension solution. CDA has been demonstrated efficient with the result of low deconvolution error and good performance when the dimensions of the solution space are 500 × 500 = 250,000, compared with other optimization algorithms [33]- [36].

B. COMPARISONS WITH RELATED METHODS
So far, systematic comparisons of structural connectivity and resting-state functional connectivity have gained significant success in using the observed dMRI structural strengths to predict the resting-state fMRI correlations, including correlations between regions that are not directly structurally coupled. One previous research evidenced that the fMRI correlations can be better predicted by the topology of the shortest (and presumably most efficient) structural paths [14]. The results showed that both search information (accessibility of a path) and path transitivity (the density of local detours along a path) along the shortest path can predict the strength of FC among both connected and unconnected node pairs well. Here our study showed that the presence of FC between brain regions cannot be exclusively attributed to the shortest paths. A large portion of FC is induced by indirect paths, particularly by two-step paths linking two brain areas rather than the local detours along the path. Furthermore, the shortest paths can highly capture the anatomical structure of the human brain (Fig.5), indicating that there is more likely a direct structural link between two regions with shorter shortest path length even though the functional correlation between the two regions is relatively weak or even negative.
In a study by Feizi and colleagues [32], they also formulated the inference of the direct relationship from the observed correlations as a network deconvolution problem and derived a closed-form solution based on eigen-decomposition. Meanwhile, Robinson et al [30] presented the similar solution based on linear neural field theory and applied it to the determination of effective brain connectivity from FC. The eigen-decomposition based network deconvolution (ED-ND) models the observed correlation as the combination of direct and indirect dependencies owing to transitive effects of correlations of all indirect paths with arbitrary length, and therefore the solution can be solved analytically with infinite-series sums. However, their results can be achieved only on condition that the largest absolute eigenvalue of the direct connection matrix is strictly less than one, otherwise the observed correlation matrix needs rescaling. Unfortunately, it does not apply to the human brain networks. In fact, the human brain functional network is much more intricately organized, in which the observed functional correlations are related not VOLUME 8, 2020 FIGURE 8. Distributions of the negative logarithm of the predicted structural connection strengths between region pairs for the three datasets. The predicted connection strengths range from 3.5476e-9 (-log 10 (3.5476e-9) = 8.4501) to 0.4635 (-log 10 (0.4635) = 0.3340) for the 90-ROI dataset, 3.4087e-9 (-log 10 (3.4087e-9) = 8.4674) to 0.3438 (-log 10 (0.3438) = 0.4637) for the 246-ROI dataset, and 3.1997e-5 (-log 10 (3.1997e-5) only to the indirect effects by transitivity but also to the modular effects owing to the modularized network topological structure. To address such a daunting issue, we proposed a bio-inspired network deconvolution algorithm based on the mechanism of cell differentiation for distinguishing the direct dependencies from the observed functional correlations. Results on all the three connectome datasets demonstrate that the CDA-ND far outperforms the ED-ND in predicting the human brain SC from FC (Fig.7).
In addition, it should be noted that our findings revealed different statistical results on the ratio of each component to the whole FC in contrast to previous studies [21], [49], in which FC was simulated from increasingly computational models. Our approach obtained the statistical results through removing each type of effect from the whole FC. Moreover, our study also showed that the establishment of the indirect effects on FC rests on the modular effects rather than the anatomical connectivity directly.

C. LIMITATIONS AND FUTURE WORK
Although the modeling approach and the network deconvolution algorithms proposed in this paper have shown promising results, they also have some limitations. First, the relationship of FC and SC is so complex that there are still many unresolved challenges and open debates. Some studies even challenged the point that the structural and functional connectivity are related in a straightforward manner. For example, O'Reilly and colleagues [50] demonstrated relatively intact interhemispheric functional connectivity in a macaque brain in the absence of major commissural fibers; Uddin et al. [51], [52] partly characterize residual functional connectivity between two hemispheres in a complete commissurotomy patient. All these findings indicate that the mechanisms shaping the relationship between structure and function especially the FC between distant cortical regions still remain open issues. Therefore, the proposed FC model may not fully capture the relationship between the direct and indirect effects of the functional network. However, it can be improved with more insights gleaned from both empirical and theoretical findings.
Secondly, the performance of our approach is prone to the accuracy of the fMRI data acquisition and dMRI tractography. More noise may be produced by fMRI and dMRI for datasets with higher resolutions.
In addition, the negative functional correlations were not considered in the current study since the mechanism of anticorrelations still remains elusive.
However, it is worth highlighting that the proposed CDA-ND algorithm is promising although it takes longer time when dealing with large-scale network. In fact, the CDA really works when tackling super high dimensional optimization problems. To our best knowledge, an optimization algorithm that can find an optimal solution globally in a nonlinear space with at least 250,000 dimensions has never been reported before and it will find more unexpected applications in revealing direct dependencies for some other types of complex networks.
In the future, we will apply the proposed model and algorithms to optimize some noisy or falsely measured as well as to infer some missed interhemispheric structural connections measured with dMRI, which will produce a better human brain structural map in combination with the current diffusion MRI tractography.

D. CONCLUSIONS
We present an opposite approach to investigating the structure-function relationship by predicting SC from rsFC. Our method relied on one presupposition-the direct anatomical links, indirect pathways, and module topology interact with one another forming temporally dependent FC. Based on CDA optimization, we demonstrated that empirical SC can be reliably predicted, provided that both SPL and distance kernel function (i.e, KWT) constraints are set in the deconvolution process. The indirect length2-paths and module interplay are found to account for a large proportion effect with above 60% and 20% of intrahemispheric FC across three datasets, respectively. The direct SC, on the contrary, is suggested to play a small role (less than 10%) in shaping intrahemispheric FC. By analyzing the different contributing factors, our study leads to a better understanding of how the underlying anatomy configures functional networks and points out the importance of different effects.