# Brain Cortical Functional Gradients Predict Cortical Folding Patterns via
Attention Mesh Convolution

Li Yang
School of Automation
Northwestern Polytechnical University
Xi’an, China 710072
<EMAIL_ADDRESS>
Zhibin He
School of Automation
Northwestern Polytechnical University
<EMAIL_ADDRESS>
Changhe Li
School of Automation
Northwestern Polytechnical University
<EMAIL_ADDRESS>
Junwei Han
School of Automation
Northwestern Polytechnical University
<EMAIL_ADDRESS>
Dajiang Zhu
University of Texas at Arlington
<EMAIL_ADDRESS>
Tianming Liu
University of Georgia
<EMAIL_ADDRESS>
Tuo Zhang
School of Automation
Northwestern Polytechnical University
<EMAIL_ADDRESS>

###### Abstract

A close relation between brain function and cortical folding has been
demonstrated by macro-/micro- imaging, computational modeling, and genetics.
Since gyri and sulci, two basic anatomical building blocks of cortical folding
patterns, were suggested to bear different functional roles, a precise mapping
from brain function to gyro-sulcal patterns can provide profound insights into
both biological and artificial neural networks. However, there lacks a generic
theory and effective computational model so far, due to the highly nonlinear
relation between them, huge inter-individual variabilities and a sophisticated
description of brain function regions/networks distribution as mosaics, such
that spatial patterning of them has not been considered. To this end, as a
preliminary effort, we adopted brain functional gradients derived from
resting-state fMRI to embed the "gradual" change of functional connectivity
patterns, and developed a novel attention mesh convolution model to predict
cortical gyro-sulcal segmentation maps on individual brains. The convolution
on mesh considers the spatial organization of functional gradients and folding
patterns on a cortical sheet and the newly designed channel attention block
enhances the interpretability of the contribution of different functional
gradients to cortical folding prediction. Experiments show that the prediction
performance via our model outperforms other state-of-the-art models. In
addition, we found that the dominant functional gradients contribute less to
folding prediction. On the activation maps of the last layer, some well-
studied cortical landmarks are found on the borders of, rather than within,
the highly activated regions. These results and findings suggest that a
specifically designed artificial neural network can improve the precision of
the mapping between brain functions and cortical folding patterns, and can
provide valuable insight of brain anatomy-function relation for neuroscience.

##  1 Introduction

Relation between brain function and cortical folding has long been intriguing
due to its importance in offering profound insight into brain functional
operation mechanisms, brain development and evolution, as well as being of
clinical applications [1]. Gyri and sulci, the convex and concave cortical
folding patterns, are two basic anatomical units of cortical morphology. It
has been found that gyri/sulci align well with cortical functional boundaries
[2], at least in primary functional cortex. Specific gyro-sulcal patterns,
such as the $plis$ $de$ $passage$ gyri [3], topology of local gyri and sulci
(such as H-shape sulci in frontal lobe [1], H-shape gyri in fusiform gyri [4]
and the ventromedial prefrontal cortex (vmPFC) in DMN region [5]), were
suggested to correspond to subtypes of a specific function (such as cognitive
and perception ability associated with frontal lobe and reading ability
associated with fusiform gyri) or relate to the placement of a specific
functional region (such as functional activation location in vmPFC). In [6],
functional MRI signals were suggested to have stronger low frequency
components on gyri, whereas sulci have signals of stronger high frequency
components. Neuronal difference in their distributions [7], types [8], wiring
patterns [9] and even morphology [7] were reported to be present between gyri
and sulci, which was further suggested to lead to differentiated cell-level
functioning patterns. In addition, expression of a huge number of genes was
reported to be different between gyri and sulci [10]. Among them,
interestingly, the same group of genes can lead to the formation of gyri/sulci
on gyrencephalic species while mediating brain function distribution patterns
on lissencephaly species [10], suggesting a possible deeply rooted bond
between cortical folding pattern and brain functions. Importantly, the
functional differentiation between gyri and sulci could also provide clues to
artificial neural network design. For example, gyri were suggested to serve as
functional hubs and sulci as local processors, whereas sulci were isolated by
gyri on the cortical sheet [11]. Recalling that graph architecture of an
artificial neural network in the sweet spot resembles the real neural networks
[12], such a gyro-sulcal (hub-nonhub) spatial distribution might suggest an
alternative approach for neural network design.

In spite of the importance in anatomy-function relation, no explicit
explanation or generic theory of such a relation can be satisfactorily applied
across brain regions and species [2][13]. Also, very few computational models
are found [6] to provide a precise prediction from brain function to cortical
folding. The reasons could be attributed to many folds: 1) the mapping between
brain function and folding pattern is far from linear; 2) inter-individual
variabilities are huge in the layout of both gyro-sulcal patterns and
functional regions [2][13]; 3) brain function was usually delineated by
atlases where “abrupt” boundaries were used to segregate functional “units” or
“networks”, the functional attributes within which were considered to be of
the same. However, the transition of a brain function from one to another was
suggested to be “smooth” and “gradual”, which also gains supports from the
related genetic expression [14]; 4) the spatial patterning of brain function
distribution was not fully used to predict cortical folding.

To answer the aforementioned questions, as a preliminary effort, we proposed
to used artificial neuronal networks to investigate the predictive ability of
brain functional gradients to cortical folding patterns at macro scale. By
means of resting-state fMRI (rsfMRI) and structural MRI (sMRI) in Human
Connectome Project (HCP) dataset [15], we estimated the global functional
gradients, embedding axes that encode the “gradual” differences in mesh
vertices’ connectivity patterns and can be present on a cortical sheet (top-
left in Figure 1), where the prediction of gyro-sulcal segmentation (top-right
in Figure 1) was performed. Technically, we proposed to use a mesh convolution
model[16] based on the classic U-net architecture, which has been widely and
successfully applied to medical image segmentation and is a close analogy with
our task. In addition, to enhance the interpretability of the mapping from
functional gradients to cortical folding maps, we added a channel attention
block [17] to the head of the U-net mesh convolution model, where functional
gradients, as multi-channels, were re-weighted by the attention block before
being fed to mesh convolution. The weights can encode the contribution of
functional gradients to cortical folding prediction. Overall, the main
contributions of this work are:

  * •

We proposed a new concept of adopting functional gradients to encode
functional connectivities and present them on the cortex, such that mesh
convolution can be adopted to predict gyri/sulci from spatial patterns of
functional connectivity embedding.

  * •

Channel attention was designed to enhance the interpretability of the
nonlinear mapping between brain function and cortical folding.

  * •

We explored and defined the activation maps of different layers and found
novel clues to function-anatomy relation, which is valuable for both real and
artificial neural networks.

##  2 Related Works

###  2.1 Mesh Convolution

Meshes are composed of three distinct types of geometric primitives: vertices,
edges, and faces. The classification on meshes is a basic research area which
has made great progress with the development of deep learning. DiffusionNet
[18] introduces a general-purpose approach to deep learning on 3D surfaces.
The networks can be discretized on various geometric representations such as
triangle meshes or point clouds, and can even be trained on one representation
and then applied to another. HodgeNet [19] uses many features that make it an
attractive alternative for learning from meshes. During inference, its
structure resembles that of most spectral geometry processing algorithms:
construct a useful operator, and compute features from its spectrum. However,
the above two approaches, as well as many other methods, including Pointnet
[20], MeshCNN [21], PD-MeshNet [22], MeshWalker [23], PFCNN [24], need
geometric features of meshes, such as coordinates of vertices, which is not
applicable to our task that learns and predicts "texture" features on meshes.
BrainSurfCNN [25] and SubdivNet [16] provides a solution to the problem.
BrainSurfCNN is a method which uses spherical convolutional kernel [26] to
predict the task fMRI contrasts from resting-state functional connectivity
ones under a regression framework. The work in [26] uses only one-ring
neighborhood for convolution, which limits the extraction of long-range
features. SubdivNet [16] uses a regular and uniform downsampling scheme to
establish a fine-to-coarse mesh hierarchy. The convolutions efficiently
support stride and large dilation, allowing the model to better capture long-
range features. At the same time, the SubdivNet is a general approach, in
which one can easily use 3D CNN on the mesh as similar as using the 2D CNN on
the image without geometric features. In our method, we choose the convolution
operation of SubdivNet to learn the "texture" features, rather than geometric
features, for the purpose of predicting cortical folding patterns from brain
cortical functional gradients.

###  2.2 Attention Mechanisms

Attention mechanisms were introduced into computer vision with the aim of
effectively finding the salient regions by imitating the human visual system.
Recently, a large number of attention mechanisms and approaches for deep
learning are proposed to improve the performance of classification. Channel
attention adaptively recalibrates the weight of each channel, and determines
what to pay attention to. Hu and colleagues [17] proposed the concept of
channel attention and presented SENet for this purpose. The core of SENet is a
squeeze-and-excitation (SE) block which is used to collect global information,
capture channel-wise relationships and improve representation ability. An SE
block can be added after each residual unit due to their low computational
resource requirements, which is suitable for our model. Wang and colleagues
[27] proposed the efficient channel attention (ECA) block which uses a 1D
convolution to determine the interaction between channels, instead of
dimensionality reduction, which can directly model correspondence between
weight vectors and inputs. However, it is noticed that the ECA block is not
very suitable for models with a small number of channels in our experiments
(as demonstrated in the Table 1 and Figure 2). [28] proposed the convolutional
block attention module (CBAM) which stacks channel attention and spatial
attention in series. However, the method of spatial attention which only fits
2D image can not be directly used in mesh. As a result, the SEnet was chose in
our work.

###  2.3 Predicting Gyro-sulcus From fMRI

So far, there are very few works in predicting cortical folding pattern from
fMRI via a deep learning model. Liu and colleagues [6] designed a
convolutional neural network (CNN) based classifier, which can differentiate
gyral and sulcal fMRI signals with a reasonable accuracy (0.67 on resting
state fMRI). However, in this model, vertices on surfaces are samples (with a
gyrus/sulcus label) independent from each other, whereas functional regions
are not spatially isolated and independent. Either functional network studies
or functional gradient ones demonstrated that remote or neighboring regions
cooperate with each other to elicit and maintain a brain function [14][29].
Therefore, functional connectivity induced gradients that encode the spatial
embedding of brain neuronal interactivity could be more suitable for
prediction of spatial layout of gyri/sulci. Also, the relatively lower
classification accuracy in [6] could be ascribed to the direct use of noisy
fMRI signals as predictive features, whereas connectivities could mitigate
such side effects and direct the attention to brain functions at a global
scale rather than a vertex-wise one. Zhao and colleagues [30] proposed a novel
Hierarchical Interpretable Autoencoder (HIAE) to explore the hierarchical
functional difference between gyri and sulci. The hierarchical features
learned by autoencoder can be embedded into a one-dimensional vector which
interprets the features as spatial-temporal patterns. However, its major aim
is embedding hierarchical and spatial-temporal functional features into a one-
dimensional vector, whereas the functional difference between gyri and sulci
was only a _post hoc_ analysis.

Figure 1: Attention mesh convolution network architecture

##  3 Method

In general, we designed a novel attention mesh convolution method based on
U-net architecture to predict the gyro-sulcal segmentation map from brain
functional gradient maps. The model’s inputs are multi-channel functional
gradients on icosahedral mesh [31], while its output is a two-channel binary
map on the same icosahedral mesh (Figure 1).

###  3.1 Functional Gradients and Cortical Folding Patterns

In brief, for each subject, a functional connectivity matrix was calculated
using the correlation coefficient for all pairs between vertices (32,492 per
hemisphere) on HCP grayordinate mesh system, introduced in section 4.1. For
each row of this “dense” functional connectome matrix, the values of the top
10% of connections were retained, whereas all others were zeroed. The negative
values in the remaining connections were zeroed as well. The similarity of the
connectivity matrix $\mathbf{A}$ was ensured by calculating the cosine
distance between all pairs of rows, yielding a symmetric affinity matrix
$\mathbf{L}$. On this matrix, diffusion embedding [32] that allows local and
long distance connections to be more effectively projected into a common
space, was employed to yield unitless components, also known as functional
“gradients”. The positions of nodes on these embedding axes encode the
differences in nodes’ connectivity patterns. In this work, the number of
gradients was set to 5. Results of other numbers of gradients were compared.
It is noted that each functional gradient is associated with an eigen value
that suggests its dominance. It has been demonstrated that gradients
decomposed from different subjects are consistent, whereas their dominance
could vary across subjects [33]. Hence, the algorithm in [33] used the group-
mean matrix to yield the group-wise gradients, that were further used as
templates to identify the cross-subject correspondence of gradients from
individuals. The top-left panel of Figure 1 illustrates the aligned 5
functional gradient maps from one subject. Surface mean curvature _k_ was used
to segment surface to gyri (_k_ >0) and sulci (_k_ <0) (the right-top panel of
Figure 1).

###  3.2 Channel Attention Block

In order to boost the representation power of a network, the channel attention
block SEnet [17] was applied in front of mesh convolution. In this attention
block, a self-gating mechanism based on channel dependence yields the
excitation of each channel (functional gradient in this work), such that
channels of input were then re-weighted and fed to the following layers. This
weight vector $\mathrm{w}$ could encodes the respective contribution of
functional gradients to folding pattern prediction.

As shown in Figure 1 and Equation(1&2), the $\mathbf{F}$ $(\cdot)$ is global
average pooling function which compresses the characteristics of each channel
into a real number, where $\mathrm{c}$ denotes the channel number, $N_{f}$
denotes the number of mesh faces, $G_{c}$ denotes the channel’s gradient
contrasts. Then, a $\mathrm{1*c}$ vector was sent to
$\mathbf{F}_{\mathrm{ex}}(\cdot,\mathrm{w})$ which was implemented by a fully
connected layer. In our method, we encoded 5 channels to 2 channels and then
decoded them back to 5, in which $\mathrm{w}$ is a learnable parameter. After
these operations, $\mathrm{w}$ gained the weight for every input channel, and
then $\mathbf{F}_{\mathrm{scale}}(^{\cdot},^{\cdot})$ scaled the
$\mathrm{1*c}$ vector and make the gradient contrasts on every mesh face
multiplied by the corresponding channel weight. Finally, the re-weighted
channels were fed to the mesh U-Net.

$z=\mathbf{F}\left(\mathbf{G}_{c}\right)=\frac{1}{N_{f}}\sum_{i=1}^{N_{f}}G_{c}(i)$
(1)
$\mathbf{s}=\mathbf{F}_{ex}(\mathbf{z},\mathbf{W})=\sigma(g(\mathbf{z},\mathbf{W}))=\sigma\left(\mathbf{W}_{2}\delta\left(\mathbf{W}_{1}\mathbf{z}\right)\right)$
(2)

###  3.3 Mesh CNN with U-net Architecture and Loss Function

The model is based on the U-Net architecture [34] using the mesh convolutional
kernel proposed in [16]. As shown in Figure 1, after the attention block, a
’meshlinear’ operation was used to generate 64 channels. This operation is a
linear operation on meshes which resembles ’1D’ convolution on images. Then,
the results above were sent to the classic U-Net framework. The kernel size
was set to 3 on all mesh convolution operations, which was designed in [16],
and the ’mean pooling’ was used for downsampling. A difference between the
classic U-Net framework and ours is that a ’meshlinear’ operation was used
after the third ’pooling’ operation (yellow arrow in Figure 1), and we found
this change improves the accuracy of prediction. Finally, the deep network
yielded two output channels. The loss was calculated by cross-entropy cost
function, and the final result was obtained by an ’argmax’ function, in which
each mesh face was classified as sulcus or gyrus (Equation(3)).

$L=-[y\log\widehat{y}+(1-y)\log(1-\widehat{y})]$ (3)

The segmentation accuracy was defined in Equation(4), in which the ${N_{f}}$
denotes the number of mesh faces, $\mathbf{prd}_{i}$ and $\mathbf{lab}_{i}$
denotes the $i^{th}$ face’s prediction and label, respectively. $sign$
function outputs 0 and 1, according to whether $\mathbf{prd}_{i}$ is equal to
$\mathbf{lab}_{i}$ or not.

$Acc=\frac{1}{N_{f}}\sum_{i=1}^{N_{f}}sign(\mathbf{prd}_{i},\mathbf{lab}_{i})$
(4)

##  4 Experiments

###  4.1 Dataset

We used the 3-Tesla T1-weighted MRI and resting-state fMRI (rsfMRI) data from
Human Connectome Project (HCP) dataset [15]. These data have been minimally
pre-processed upon release [35]. Important imaging parameters are: T1-weighted
MRI: TR=2400 $ms$, TE=2.14 $ms$, flip angle=8 $deg$, image matrix=260×311×260
and resolution=0.7×0.7×0.7 $mm^{3}$. Diffusion MRI: TR=5520 $ms$, TE=89.5
$ms$, flip angle=78 $deg$, FOV=210×180 $mm^{2}$, matrix=168×144,
resolution=1.25×1.25×1.25 $mm^{3}$, echo spacing=0.78 $ms$; rsfMRI: TR=0.72
$s$, TE=33.1 $ms$, flip angle=52 $deg$, in-plane FOV=208×180 $mm^{2}$,
matrix=90×104, 220 $mm$ FOV, 72 slices, 1200 time points, 2.0 $mm$ isotropic
voxels, BW = 2290 $Hz/Px$. T1-weighted MRI was used to reconstruct a cortical
surface, since it provides more precise anatomical information due to its
higher spatial resolution. Preprocessing of T1-weighted MRI includes skull
stripping, tissue-segmentation and white matter surface reconstruction. The
surfaces were transferred to the grayordinate system, a common surface
template, _via_ cross-subject registration to standard volume and surface
spaces. It is noted that the aligned surface of each subject was re-sampled
such that vertices have cross-subject correspondence. Preprocessing of fMRI
signals includes: ICA-denoise, spatial artifact/distortion removal. Within-
subject cross-modal registration was applied such that each grayordinate
vertex was associated with a preprocessed fMRI signal. More details of data
acquisition and preprocessing are referred to [35]. Our experiments included
348 samples for training/testing (0.8/0.2).

Figure 2: Comparison experiments between ground truth and predicated results
_via_ different methods. A ’U’ shaped gyral area is enlarged for inspection in
detail, where yellow circles highlight the unsatisfied predictions.
Abbreviations: G for gradients

###  4.2 Implementation Details

Firstly, to adapt to the requirement in Subdivnet, we re-meshed the original
surface and transferred the gradient features as well as class labels (0 for
gyri and 1 for sulci) from mesh vertices to the mesh faces. The batch size and
learning rate was set to 6 and 0.02, respectively. ’Bilinear’ upsample
operation was used. The accuracy on testing samples converges (75.85%) after
about 20 epochs (from 15 to 21). Then, we chose the checkpoint of the best
accuracy on test data set and continued fine-tuning for several epochs, with
the learning rate set to 0.008.

###  4.3 Comparison with Mesh Convolution Baselines and Other Models

Since our model was modified from Subdivnet [16], it was used as a mesh
convolution baseline. The BrainSurfCNN [26] predicts the task fMRI contrasts
from resting-state functional connectivity ones by using a regression scheme.
For a fair comparison, we modified its’ network architecture by replacing the
100 input channels with 5 channels, the 47 output channels with 2 channels in
our model, as well as the loss function with cross-entropy loss, to implement
the same segmentation task. The features on the mesh vertices were transferred
to faces as well. A classic linear classification method, support-vector
machine (SVM) [36] was also selected for comparison. Figure 2 shows the
segmentation results from one random testing subject and Table 1 reports the
comparison of results _via_ different methods by means of a variety of
metrics. It is found that nonlinear models overwhelmingly outperform linear
SVM in all metrics. Besides the intrinsic nonlinear relation between function
and folding that could be significantly better encoded and decoded by
nonlinear models, the definition of neighboring mask, a key factor for
convolution operations on a mesh, could play an important role in the
improvement. The importance of the convolution mask was also supported by a
better performance of Subdivnet than BrainSurfCNN, since BrainSurfCNN, unlike
Subdivnet, only used a small neighborhood as convolution mask which limits the
extraction of long-range features. It is also found that sulci were often
oversegmented by Subdivnet and BrainSurfCNN (yellow circles in Figure 2) in
highly bent gyral regions, where the cortical morphology is far more complex
than others.

###  4.4 Ablation Study

We considered the effects of different types of attention models and the
number of functional gradient channels. It is found in Figure 2 and Table 2
that adding channel attention block or switching to another one does not
significantly improve the performance in MIoU, Precision, F1 and Recall,
whereas the accuracy could be affected by the selection of a channel attention
block. The choice or absence of a channel attention could yield wrong
segmentation in "straight" gyral/sulcal regions (the left yellow circles in
third row of Figure 2). Neither fewer (4G) nor more (6G) functional gradients
improve the segmentation performance. Mis-classification with 4G and 6G was
usually found on the border between gyri and sulci (yellow circles in the
fourth row of Figure 2).

Table 1: Comparison with other methods. The MIoU, Precision, F1-score, Recall and Accuracy are reported. Method | MIoU | Precision | F1-score | Recall | Accuracy
---|---|---|---|---|---
Our Method | 0.603 | 0.752 | 0.751 | 0.750 | 75.85%
Subdivnet | 0.595 | 0.744 | 0.745 | 0.746 | 75.11%
BrainSurfCNN | 0.464 | 0.657 | 0.622 | 0.626 | 65.03%
Linear SVM | 0.295 | 0.488 | 0.371 | 0.499 | 59.09%
Table 2: Ablation study on attention block and gradients. The MIoU, Precision, F1-score, Recall and Accuracy are reported. Red and blue denote the best and the second-best results, respectively. Abbreviations: G for gradients. Type | MIoU | Precision | F1-score | Recall | Accuracy |
---|---|---|---|---|---|---
SEnet+5Gs | 0.603 | 0.752 | 0.751 | 0.750 | 75.85% |
ECAnet+5Gs | 0.599 | 0.750 | 0.747 | 0.745 | 75.61% |
No Attention+5Gs | 0.598 | 0.748 | 0.747 | 0.746 | 75.45% |
SEnet+4Gs | 0.587 | 0.744 | 0.737 | 0.734 | 74.89% |
SEnet+6Gs | 0.580 | 0.732 | 0.732 | 0.732 | 73.86% |

##  5 Discussion

Figure 3: (a) Activation maps in different deep layers that resemble the
fronto-parietal functional network, not detected by functional gradients. (b)
Similarity between functional gradients (G for short) and curvature map, and
the weights of Gs learnt from attention block. (c) The segmentation activation
maps (red for more activated regions) in the last two-channel layer and
cortical landmarks (cyan bubbles), including gyral peaks and sulcal pits, that
are independent from curvature.

The attention block re-weights the five gradients. The latter gradients (3G
and 5G) gain relatively higher weights (Figure 3(b)), and their similarity to
curvature map (measured by Pearson correlation) is significantly lower than
the dominant gradients (1G and 2G). These results suggest that the implicit
relation between function and folding might not be encoded in dominant
gradients. It is also interesting to find that some functional networks, which
are not explicitly represented by input gradients, can be identified in
activation maps in deep layers. As an example, the fronto-parietal network,
one of the resting-state networks [37] not present in the 5 gradients, were
found in both encoding and decoding layers, as shown in Figure 3(a). Finally,
the activation maps in the last two-channel layer, were related to folding
landmarks, gyral peaks (local convexity maxima) and sulcal pits (local
concavity maxima) [38], which were obtained independent of curvature used to
yield the output gyro-sulcal maps. These landmarks (blue bubbles in Figure
3(c), group-wise mean locations across subjects) are not located on the highly
activated regions (red regions). On average, activation values on these
landmarks are 0.46±0.62 and 0.60±0.55 (_p_ <0.05, _via_ unpaired _t_ -test) on
other regions (expect the gray background), suggesting that these landmarks
might mark the borders of folding patterns that are highly correlated with
brain functions.

##  6 Conclusion

In this work, we developed an attention mesh convolution model based on the
U-net architecture to predict cortical gyro-sulcal segmentation maps from
brain functional gradients derived from resting-state fMRI. The model
architecture, such as input channel, output channel, feature types and
attention block, can be easily extended or customized according to specific
needs. Unlike the previous study, where functional signals from gyral/sulcal
vertices were independently used as samples to a CNN model and their frequency
band was found to be a discriminative feature [6], the mesh convolution in our
model considers the spatial organization of functional gradients and folding
patterns on a cortical sheet. The prediction performance _via_ our model
outperforms the state-of-the-art. In addition, we found that the dominant
functional gradients contribute less to folding prediction. On the activation
maps of the last layer, we found that some well-studied cortical landmarks are
on the borders of, rather than within, the highly activated regions. These
results and findings suggest that a specifically designed artificial neural
network can improve the precision of the mapping between brain functions and
cortical folding patterns, and can provide valuable clues of brain anatomy-
function relation for neuroscience and artificial neural network design.

## References

  * [1] Vanessa Troiani, Marisa A Patti, and Kayleigh Adamson.  The use of the orbitofrontal h-sulcus as a reference frame for value signals.  European Journal of Neuroscience, 51(9):1928–1943, 2020.
  * [2] Bruce Fischl, Niranjini Rajendran, Evelina Busa, Jean Augustinack, Oliver Hinds, BT Thomas Yeo, Hartmut Mohlberg, Katrin Amunts, and Karl Zilles.  Cortical folding patterns and predicting cytoarchitecture.  Cerebral cortex, 18(8):1973–1980, 2008.
  * [3] Jean Régis, Jean-François Mangin, Taku Ochiai, Vincent Frouin, Denis Rivière, Arnaud Cachia, Manabu Tamura, and Yves Samson.  “sulcal root” generic model: a hypothesis to overcome the variability of the human cortex folding patterns.  Neurologia medico-chirurgica, 45(1):1–17, 2005.
  * [4] Arnaud Cachia, Margot Roell, Jean-François Mangin, Zhong Yi Sun, Antoinette Jobert, Lucia Braga, Olivier Houde, Stanislas Dehaene, and Grégoire Borst.  How interindividual differences in brain anatomy shape reading accuracy.  Brain Structure and Function, 223(2):701–712, 2018.
  * [5] Alizée Lopez-Persem, Lennart Verhagen, Céline Amiez, Michael Petrides, and Jérome Sallet.  The human ventromedial prefrontal cortex: sulcal morphology and its influence on functional organization.  Journal of Neuroscience, 39(19):3627–3639, 2019.
  * [6] Huan Liu, Shu Zhang, Xi Jiang, Tuo Zhang, Heng Huang, Fangfei Ge, Lin Zhao, Xiao Li, Xintao Hu, Junwei Han, et al.  The cerebral cortex is bisectionally segregated into two fundamentally different functional units of gyri and sulci.  Cerebral Cortex, 29(10):4238–4252, 2019.
  * [7] Claus C Hilgetag and Helen Barbas.  Developmental mechanics of the primate cerebral cortex.  Anatomy and embryology, 210(5):411–417, 2005.
  * [8] Brian G Rash, Alvaro Duque, Yury M Morozov, Jon I Arellano, Nicola Micali, and Pasko Rakic.  Gliogenesis in the outer subventricular zone promotes enlargement and gyrification of the primate cerebrum.  Proceedings of the National Academy of Sciences, 116(14):7089–7094, 2019.
  * [9] Gang Xu, Andrew K Knutsen, Krikor Dikranian, Christopher D Kroenke, Philip V Bayly, and Larry A Taber.  Axons pull on the brain, but tension does not drive cortical folding.  Journal of biomechanical engineering, 132(7), 2010.
  * [10] Camino de Juan Romero, Carl Bruder, Ugo Tomasello, José Miguel Sanz-Anquela, and Víctor Borrell.  Discrete domains of gene expression in germinal layers distinguish the development of gyrencephaly.  The EMBO journal, 34(14):1859–1874, 2015.
  * [11] Fan Deng, Xi Jiang, Dajiang Zhu, Tuo Zhang, Kaiming Li, Lei Guo, and Tianming Liu.  A functional model of cortical gyri and sulci.  Brain structure and function, 219(4):1473–1491, 2014.
  * [12] Jiaxuan You, Jure Leskovec, Kaiming He, and Saining Xie.  Graph structure of neural networks.  In International Conference on Machine Learning, pages 10881–10891. PMLR, 2020.
  * [13] Martin A Frost and Rainer Goebel.  Measuring structural–functional correspondence: spatial variability of specialised brain regions after macro-anatomical alignment.  Neuroimage, 59(2):1369–1381, 2012.
  * [14] Julia M Huntenburg, Pierre-Louis Bazin, and Daniel S Margulies.  Large-scale gradients in human cortical organization.  Trends in cognitive sciences, 22(1):21–31, 2018.
  * [15] David C Van Essen, Stephen M Smith, Deanna M Barch, Timothy EJ Behrens, Essa Yacoub, Kamil Ugurbil, Wu-Minn HCP Consortium, et al.  The wu-minn human connectome project: an overview.  Neuroimage, 80:62–79, 2013.
  * [16] Shi-Min Hu, Zheng-Ning Liu, Meng-Hao Guo, Jun-Xiong Cai, Jiahui Huang, Tai-Jiang Mu, and Ralph R Martin.  Subdivision-based mesh convolution networks.  ACM Transactions on Graphics (TOG), 41(3):1–16, 2022.
  * [17] Jie Hu, Li Shen, and Gang Sun.  Squeeze-and-excitation networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 7132–7141, 2018.
  * [18] Nicholas Sharp, Souhaib Attaiki, Keenan Crane, and Maks Ovsjanikov.  Diffusionnet: Discretization agnostic learning on surfaces.  ACM Transactions on Graphics (TOG), 41(3):1–16, 2022.
  * [19] Dmitriy Smirnov and Justin Solomon.  Hodgenet: learning spectral geometry on triangle meshes.  ACM Transactions on Graphics (TOG), 40(4):1–11, 2021.
  * [20] Charles R Qi, Hao Su, Kaichun Mo, and Leonidas J Guibas.  Pointnet: Deep learning on point sets for 3d classification and segmentation.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 652–660, 2017.
  * [21] Rana Hanocka, Amir Hertz, Noa Fish, Raja Giryes, Shachar Fleishman, and Daniel Cohen-Or.  Meshcnn: a network with an edge.  ACM Transactions on Graphics (TOG), 38(4):1–12, 2019.
  * [22] Francesco Milano, Antonio Loquercio, Antoni Rosinol, Davide Scaramuzza, and Luca Carlone.  Primal-dual mesh convolutional neural networks.  Advances in Neural Information Processing Systems, 33:952–963, 2020\.
  * [23] Alon Lahav and Ayellet Tal.  Meshwalker: Deep mesh understanding by random walks.  ACM Transactions on Graphics (TOG), 39(6):1–13, 2020.
  * [24] Yuqi Yang, Shilin Liu, Hao Pan, Yang Liu, and Xin Tong.  Pfcnn: Convolutional neural networks on 3d surfaces using parallel frames.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 13578–13587, 2020.
  * [25] Gia H Ngo, Meenakshi Khosla, Keith Jamison, Amy Kuceyeski, and Mert R Sabuncu.  From connectomic to task-evoked fingerprints: Individualized prediction of task contrasts from resting-state functional connectivity.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 62–71. Springer, 2020.
  * [26] Chiyu Jiang, Jingwei Huang, Karthik Kashinath, Philip Marcus, Matthias Niessner, et al.  Spherical cnns on unstructured grids.  arXiv preprint arXiv:1901.02039, 2019.
  * [27] Q. Wang, B. Wu, P. Zhu, P. Li, W. Zuo, and Q. Hu.  Eca-net: Efficient channel attention for deep convolutional neural networks.  2019\.
  * [28] Vittorio Ferrari, Martial Hebert, Cristian Sminchisescu, and Yair Weiss.  Computer Vision–ECCV 2018: 15th European Conference, Munich, Germany, September 8–14, 2018, Proceedings, Part V, volume 11209.  Springer, 2018.
  * [29] Jessica S Damoiseaux, SARB Rombouts, Frederik Barkhof, Philip Scheltens, Cornelis J Stam, Stephen M Smith, and Christian F Beckmann.  Consistent resting-state networks across healthy subjects.  Proceedings of the national academy of sciences, 103(37):13848–13853, 2006.
  * [30] Lin Zhao, Haixing Dai, Xi Jiang, Tuo Zhang, Dajiang Zhu, and Tianming Liu.  Exploring the functional difference of gyri/sulci via hierarchical interpretable autoencoder.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 701–709. Springer, 2021.
  * [31] John R Baumgardner and Paul O Frederickson.  Icosahedral discretization of the two-sphere.  SIAM Journal on Numerical Analysis, 22(6):1107–1115, 1985.
  * [32] Daniel S Margulies, Satrajit S Ghosh, Alexandros Goulas, Marcel Falkiewicz, Julia M Huntenburg, Georg Langs, Gleb Bezgin, Simon B Eickhoff, F Xavier Castellanos, Michael Petrides, et al.  Situating the default-mode network along a principal gradient of macroscale cortical organization.  Proceedings of the National Academy of Sciences, 113(44):12574–12579, 2016.
  * [33] Reinder Vos de Wael, Oualid Benkarim, Casey Paquola, Sara Lariviere, Jessica Royer, Shahin Tavakol, Ting Xu, Seok-Jun Hong, Georg Langs, Sofie Valk, et al.  Brainspace: a toolbox for the analysis of macroscale gradients in neuroimaging and connectomics datasets.  Communications biology, 3(1):1–10, 2020.
  * [34] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In International Conference on Medical image computing and computer-assisted intervention, pages 234–241. Springer, 2015.
  * [35] Matthew F Glasser, Stamatios N Sotiropoulos, J Anthony Wilson, Timothy S Coalson, Bruce Fischl, Jesper L Andersson, Junqian Xu, Saad Jbabdi, Matthew Webster, Jonathan R Polimeni, et al.  The minimal preprocessing pipelines for the human connectome project.  Neuroimage, 80:105–124, 2013.
  * [36] Marti A. Hearst, Susan T Dumais, Edgar Osuna, John Platt, and Bernhard Scholkopf.  Support vector machines.  IEEE Intelligent Systems and their applications, 13(4):18–28, 1998\.
  * [37] Stephen M Smith, Diego Vidaurre, Christian F Beckmann, Matthew F Glasser, Mark Jenkinson, Karla L Miller, Thomas E Nichols, Emma C Robinson, Gholamreza Salimi-Khorshidi, Mark W Woolrich, et al.  Functional connectomics from resting-state fmri.  Trends in cognitive sciences, 17(12):666–682, 2013.
  * [38] Xiao Li, Songyao Zhang, Xi Jiang, Shu Zhang, Junwei Han, Lei Guo, and Tuo Zhang.  Cortical development coupling between surface area and sulcal depth on macaque brains.  Brain Structure and Function, pages 1–17, 2022.


# A principled approach for weighted multilayer network aggregation

Junyao Kuang, Caterina Scoglio Junyao Kuang is with the Department of
Electrical and Computer Engineering, Kansas State University, Manhattan, KS,
66502. E-mail<EMAIL_ADDRESS>Caterina Scoglio is with the Department of
Electrical and Computer Engineering, Kansas State University, Manhattan, KS,
66502. E-mail<EMAIL_ADDRESS>
###### Abstract

A multilayer network depicts different types of interactions among the same
set of nodes. For example, protease networks consist of five to seven layers,
where different layers represent distinct types of experimentally confirmed
molecule interactions among proteins. In a multilayer protease network, the
co-expression layer is obtained through the meta-analysis of transcriptomic
data from various sources and platforms. While in some researches the co-
expression layer is in turn represented as a multilayered network, a
fundamental problem is how to obtain a single-layer network from the
corresponding multilayered network. This process is called multilayer network
aggregation. In this work, we propose a _maximum a posteriori_ estimation-
based algorithm for multilayer network aggregation. The method allows to
aggregate a weighted multilayer network while conserving the core information
of the layers. We evaluate the method through an unweighted friendship network
and a multilayer gene co-expression network. We compare the aggregated gene
co-expression network with a network obtained from conflated datasets and a
network obtained from averaged weights. The Von Neumann entropy is adopted to
compare the mixedness of the three networks, and, together with other network
measurements, shows the effectiveness of the proposes method.

###### Index Terms:

multilayer network, gene co-expression network, network aggregation, Von
Neumann entropy

##  1 Introduction

Networks are widely used to describe the interactions between nodes within a
system. The interactions can be either represented by unweighted edges or
quantified by weights on edges. In terms of epidemic analysis, weighted
networks provide more accurate results than their unweighted counterparts
according to [1, 2]. In different types of networks, edges are detected in
various ways. For example, in a social network, an edge is detected if two
people recognize each other or meet up at the same places during a period [3,
4]. A weighted social network can be obtained by simply asking people how
often they meet their friends during a period [5, 6]. Similarly, in
constructing gene co-expression networks, the weights of edges are determined
through correlation coefficients, and an unweighted gene co-expression network
can be obtained by selecting edges with correlation coefficients greater than
a threshold [7]. In papers [7, 8, 9], the Pearson correlation is used to
determine edges for gene co-expression networks. The gene expression levels
are recorded after subjects receive different stimuli, the Pearson correlation
coefficients are calculated for every pair of genes and a fixed threshold is
used to select edges. For example, if two genes have a Pearson correlation
coefficient greater than 0.8, then an edge is detected between the two genes
[9].

However, when the data come from various sources and platforms, it is
difficult to conflate different datasets and find a single-layer network. For
example, to construct a gene co-expression network, authors in [7] construct a
network from a conflated dataset, where gene expression values are tested
under different platforms, such as RNA-sequencing and microarray. The gene
expression values are normalized through a median-shifted method, and the
Pearson correlation coefficient is used to detect edges. The drawback of the
method is that each expression value must be recorded more than once in each
condition. The correlation coefficients in the co-expression map are
significantly affected by the normalization method, and most of the
correlation coefficients are greater than 0.9, which means that most of the
genes in the network are correlated. To avoid mixing different experiments
together, the multilayer network is used to depict different types of
experiments in some circumstances [6, 8]. For examples, authors in [8, 9, 10]
use multilayer networks to represent the relations of genes. However, when
there are too many layers, it is difficult to analyze the network properties,
such as centralities and degree distribution [11, 12, 13]. Besides, in the co-
expression network, the layers obtained through different types of experiments
need to be aggregated as a single layer in a more complex multilayer network
[14]. In this work, we are interested in constructing a single layer network
by aggregating a multilayer network which is obtained from various sources.

Although various methods have been used to construct single layer network, but
the methods are either too difficult to implement or based on random criteria.
Paper [15] compares three methods for multilayer network aggregation and link
prediction, including a direct aggregation, a result aggregation, and a layer
co-analysis method. In all of the three methods, the networks are mapped to
vectors by training the _node2vec_ neural network [16]. The drawbacks of the
algorithm are obvious; it is necessary to have a large dataset to train the
neural network, and a high-performance processor is required. Authors in paper
[14] use the Von Neumann entropy to evaluate the reducibility of a multilayer
network. And the authors adopt a greedy agglomerative hierarchical clustering
algorithm to explore the space of partitions of different layers, but the
layers are aggregated directly. In paper [11, 12], the authors introduce a
selection method, which selects nodes or edges of interest with high
centralities to obtain the final single-layer network. Also, a method that
filters nodes and edges with lower centralities is introduced to obtain a
single-layer network. Paper [8] presents a method that constructs a network
for each dataset, and an edge is confirmed in if it is observed in more than
two datasets. To obtain a single layer gene co-expression network, a rank
aggregation method is introduced in [9]. In paper [17], the authors use a
statistical method for aggregating unweighted networks, in which the
probability of observing an edge is decided by the significance level obtained
by Bonferroni correction.

In this work, we propose a simple and principled method to aggregate
multilayer networks, which is the main contribution. The method introduced in
this paper is based on statistical distribution of edge weights across the
layers. Inspired by [3, 18, 19, 20], we adopt the _maximum a posteriori_ (MAP)
estimation and expectation-maximization algorithm to determine the weights of
edges in the target single-layer network. In the derivation of the proposed
method, we make the following assumptions. First, each layer could contain
part or all of the nodes of interest, and the edge weights of different layers
are non-negative. Second, we assume that the different edges in the multilayer
network are independent and identically distributed. Also, the expected weight
of each edge follows Poisson distribution. Further, the ensemble distribution
of all the weights of edges is assumed to follow the exponential distribution,
as introduced in [21].

The paper is organized as follows. In section two, we describe the model for
constructing the target layer with the MAP estimation. In Section three, to
validate the effectiveness of the proposed method, we implement the algorithm
to aggregate an _Anopheles Gambiae_ gene co-expression network and a
friendship multilayer network. For the co-expression network, we compare the
properties of the resulted network with a network constructed from conflated
datasets and a network aggregated directly, and it shows that the proposed
method can reduce experimental uncertainty. We conclude the article in section
four and discuss future works.

##  2 MAP Estimation-Based Network Aggregation

We denote the target network as $A$, and assume there are $N$ nodes across all
of the $D$ layers. Note that different layers could have different sets of
nodes. For every pair nodes $i$ and $j$ in the network $A$, we assume the
expected weight is $\lambda_{ij}$. Therefore, the problem can be simplified to
determining the parameter set $\lambda=\\{\lambda_{ij}\ |\ i\leq N,j\leq
N\\}$, given the parameters of multilayer network $D$. In this method, we use
a probabilistic model, the likelihood of the model parameters given $D$ is
$P(\lambda\ |\ D)$. Applying the Bayes theorem, we have

$\displaystyle P(\lambda\ |\ D)=\frac{P(D\ |\ \lambda)P(\lambda)}{P(D)},$ (1)

where $P(D\ |\ \lambda)$, $P(\lambda)$ and $P(D)$ are the likelihood of the
data sets $D$ under network constructed through parameter set $\lambda$, the
prior probability, and the marginal likelihood, respectively. In a particular
dataset, it is possible that only a part of the $N$ nodes are tested, so we
can only obtain the weights for some of the edges. To determine the expression
of $P(D\ |\ \lambda)$, we assume the relationship between $i$ and $j$ is
detected $k_{ij}$ times in the $D$ layers, and the weight is denoted as
$\xi_{ij}^{\alpha}$ in a particular layer $\alpha$. For the node pair $i$ and
$j$, we assume the weight of the edge follows a Poisson distribution, so the
probability of $\xi_{ij}^{\alpha}$ can be expressed as

$\displaystyle P(\xi_{ij}^{\alpha}\ |\
\lambda_{ij})=\frac{e^{-\lambda_{ij}}\lambda_{ij}^{\xi_{ij}^{\alpha}}}{\xi_{ij}^{\alpha}!}.$
(2)

Further, the likelihood function of $\lambda_{ij}$ is then

$\displaystyle L(\lambda_{ij}\ |\ \xi_{ij}^{1},...,\xi_{ij}^{k_{ij}})=$
$\displaystyle\prod_{\alpha=1}^{k_{ij}}\frac{e^{-\lambda_{ij}}\lambda_{ij}^{\xi_{ij}^{\alpha}}}{\xi_{ij}^{\alpha}!}$
$\displaystyle=$
$\displaystyle\frac{e^{-k_{ij}\lambda_{ij}}\lambda_{ij}^{\sum_{\alpha}^{k_{ij}}\xi_{ij}^{\alpha}}}{\prod_{\alpha=1}^{k_{ij}}\xi_{ij}^{\alpha}!}.$
(3)

In maximum likelihood estimation (MLE), the prior probability $P(\lambda)$ is
regarded as a constant, and the prior probability will be leave out when
applying the expectation-maximization algorithm to find the optimal
parameters. On the contrary, the prior probability is a requirement in MAP
estimation. The prior probabilities are generally provided by pilot studies
[10]. However, the prior probabilities for each pair of nodes are unknown in
our assumptions. Instead, we assume a prior distribution for the parameter set
$\lambda$. In the target network, we assume the edge weights follow
exponential distribution [21], which is given by

$\displaystyle P(\lambda_{ij})=\theta e^{-\lambda_{ij}\theta},$ (4)

where $P(\lambda_{ij})$ and $\theta$ are the prior probability of
$\lambda_{ij}$ and the rate parameter of the exponential distribution,
respectively. The posterior probability of the target network is then

$\displaystyle P(\lambda\ |\
D)=\frac{1}{P(D)}\prod_{i,j}\frac{e^{-k_{ij}\lambda_{ij}}\lambda_{ij}^{\sum_{\alpha}^{k_{ij}}\xi_{ij}^{\alpha}}\theta
e^{-\lambda_{ij}\theta}}{\prod_{\alpha=1}^{k_{ij}}\xi_{ij}^{\alpha}!}.$ (5)

However, equation 5 is still intractable. To maximize the posterior
probability, we take the logarithm form of the expression. Leaving out
constant term $P(D)$ and the factorial terms. We have

$\displaystyle L(\lambda\ |\ D)=$
$\displaystyle\sum_{i,j}[-k_{ij}\lambda_{ij}+\sum_{\alpha}^{k_{ij}}\xi_{ij}^{\alpha}log(\lambda_{ij})$
$\displaystyle+$ $\displaystyle log(\theta)-\lambda_{ij}\theta].$ (6)

We use the expectation-maximization algorithm to find the optimal values of
$\\{\lambda_{ij}\\}$. Taking the derivation of equation 2 in terms of
$\lambda_{ij}$ and $\theta$, the estimates for $\lambda_{ij}$ and $\theta$ are

$\displaystyle\theta=\frac{N^{2}}{\sum_{i,j}\lambda_{ij}},$ (7)
$\displaystyle\lambda_{ij}=\frac{\sum_{\alpha}^{k_{ij}}\xi_{ij}^{\alpha}}{k_{ij}+\theta}.$
(8)

The optimal $\lambda_{ij}$ is calculated by updating $\theta$ and
$\lambda_{ij}$ alternatively. In our experiments, random initial values are
assigned to the two parameters.

Since the edge weight is assumed to follow Poisson distribution, so the
expected edge weight between node $i$ and node $j$ in the target network is
exactly $\lambda_{ij}$. The distribution of $\\{\lambda_{ij}\\}$ in the target
network is an exponential distribution with rate parameter $\theta$.

In aggregating multilayer networks, it is difficult to evaluate the target
network, since there is no ground truth network for comparison. To intuitively
compare how mixedness of different networks, we adopt the Von Neumann entropy
itntroduced in [14]. In quantum mechanics, a system can be described by a
semidefinite positive matrix with eigenvalues summing up to one, and the Von
Neumann entropy is used to evaluate the mixedness of the system. A larger Von
Neumann entropy represents a more mixed state. According to [14], the
statistical ensemble of a network can be measured by the Von Neumann entropy.
To evaluate the distinguishability between the multilayer networks and the
aggregated network, the authors use Von Neumann entropy to study the
reducibility of multilayer networks in [14].

The Laplacian matrix _L_ of a graph has non-negative eigenvalues. Therefore,
the Laplacian matrix can be used to measure the Von Neumann entropy. The Von
Neumann entropy is defined as

$\displaystyle h_{A}=-\sum_{i=1}^{N}\gamma_{i}log(\gamma_{i}).$ (9)

where $\gamma_{i}$ is the eigenvalue of the Laplacian matrix. The elements in
the Laplacian matrix are normalized by twice the sum of edge weights. In this
work, we use the Von Neumann entropy to evaluate the mixedness of different
single-layer networks.

##  3 Numerical Results

It is difficult to validate the proposed method in section 2, since there is
no ground truth network for comparison [3]. Firstly, we aggregate a single
layer gene co-expression network and compare the properties of the network
with two single-layer networks obtained through traditional methods. Secondly,
we verify that the method does not distinguish the layers when aggregating
unweighted multilayer networks, which is the expected result.

###  3.1 Results on weighted multilayer network

In the first case, we construct a single layer gene co-expression network
through the MAP method introduced in section 2 (for short, we call it MAPNet),
and the nodes in the network are genes from _Anopheles Gambiae_. To compare
the results, we construct a network (ConNet) from the conflated datasets, and
a network (AveNet) with edge weights averaged across all of the layers.

####  3.1.1 Data pre-process

In detecting gene co-expression networks, the genes are tested under different
conditions, and the experiments are conducted under various platforms [22, 23,
24, 25, 26, 27, 28, 29]. For example, to study how the genes functioning in
different development stages, some researchers [25, 22] measure the expression
levels of _Anopheles Gambiae_ genes in different embryonic and developmental
stages. Also, some researchers [23, 26] study how blood meal and plasmodium
falciparum infections affect the expression levels of genes, etc. In this
work, we apply the construction method to 290 genes of interest which come
from 11 experiments (or datasets). Each dataset contains 12 to 20 conditions,
and there are 156 conditions in total. However, in each dataset, some genes
could have missing values or even not tested. For example, the two nodes $x$
and $y$ have the following expression values:

$x=[\times,o,\times,\times,o,\times,\times,\times,\times,\times,\times,\times,\times,\times,\times,o,\times]$

$y=[\times,\times,\times,\times,o,\times,\times,\times,\times,\times,\times,\times,\times,o,\times,o,o]$.

Where there are 17 conditions in the dataset. $"\times"$ represents expression
value, and $"o"$ represents a missing value, the number of paired elements
between nodes $x$ and $y$ is 12. Therefore, the edge weight between the two
nodes can be calculated through the Pearson correlation. Note that only gene
pairs with at least five paired elements are considered, and negative
coefficients are set to zeros. So the range of edge weights is in range [0,
1]. In our case, we normalize the edge weights to integers in the range [0,
10]. Therefore, we can obtain an eleven-layer weighted network. As a
comparison, we construct a single-layer network, i.e., ConNet, from the 156
conflated conditions. Since the 156 conditions come from 11 datasets, we
normalize the expression values of each condition by the so-called z-score
normalization, i.e., $z_{i}=\frac{x_{i}-\bar{x}}{\sigma}$, where $\sigma$ is
the standard deviation of all the expression values in the condition.
Similarly, we compute the Pearson correlation coefficients of all the gene
pairs with at least five paired elements and set negative values to zeros.

####  3.1.2 Network aggregation and comparison

We aggregate the multilayer network through the MAP method introduced in
section 2. In Fig. 1, we show the cumulative distribution function (CDF) of
edge weights in the aggregated network, as well as the corresponding
exponential distribution. the coefficient of determination is 0.93 which means
that the exponential distribution fits well the edge weights distribution.

Figure 1: CDF of the edge weights distribution and the corresponding
exponential distribution. Figure 2: The relation between the number of
datasets and rate parameters and coefficients of determination. Results are
averaged over 10 combinations of the 11 datasets.

In Fig. 2, we show how the number of layers affects the model. The red curve
shows the coefficients of determination between the obtained exponential
distribution and the distribution of the estimated $\\{\lambda_{ij}\\}$. As we
can see, the coefficients of determination are above 0.9 when there are at
least four datasets used, which means the exponential distribution fits well
the edge weights distribution. The blue curve shows the variation of the rate
parameters with respect to the number of datasets. The rate parameters
converge when there are at least four datasets. It is because when there is
more evidence, i.e., more layers, the edge weights are confirmed by more
experimental results. The estimated edge weights are not affected by systemic
errors. A network aggregated with 1 to 3 layers is highly affected by some
unpredictable factors, such as experimental biases. The single layer network
aggregated from more layers can help to improve the reliability of the
constructed network and reduce uncertainty.

Figure 3: Comparison between Pearson correlation of conflated data set and
$\lambda_{ij}$. The $\lambda_{ij}$ are obtained with 11 datasets. Figure 4:
Comparison between averaged Pearson correlation coefficients and
$\lambda_{ij}$. The $\\{\lambda_{ij}\\}$ are obtained with 11 datasets.

The comparison between the ConNet and MAPNet is shown in Fig. 3. The edge
weights of the two networks are normalized by their corresponding maximums. As
shown in the figure, the edge weights obtained through the two methods are
heterogeneous. In the neighborhood of the vertical axis, there some edges have
high correlation coefficients in the ConNet but low $\lambda_{ij}$ values in
MAPNet. This is because, in some datasets, there are some genes with less than
five paired elements. But after conflating the conditions of different
datasets together, more genes are detected with at least five paired elements.
More importantly, genes from the same dataset tend to have similar systemic
biases, such as experimental manipulation and platform differences. This is
the reason that some researchers aggregate multilayer networks to obtain a
single-layer co-expression network[8, 9] rather than constructing a network
from conflated datasets.

Figure 5: Edge weights of the aggregated network and the corresponding
exponential distribution.

We then compare the $\\{\lambda_{ij}\\}$ with the averaged Pearson correlation
coefficients, i.e., the AveNet. In Fig. 4, it can be seen that there are also
some values close to the vertical axis. The phenomenon is due to edges with
low rate of presence are given lower edge weights. Also, we can see that the
points are situated in different radial lines since edges detected in more
datasets are given higher weights in the proposed method. Finally, we compare
the node strength distribution of the three networks. The CDFs of the node
strength distributions are shown in Fig. 5. It can be seen, in the MAPNet,
there is a significant number of nodes having low node strength, and this is
accord with the results in Fig. 3 and Fig. 4.

TABLE I: Properties of the three networks  Properties | MAPNet | ConNet | AveNet
---|---|---|---
No. of nodes | 290 | 290 | 290
Ave. node strength | 50.6 | 47.7 | 68.0
Min. node strength | 0.1 | 10.7 | 11.6
Max. node strength | 89.3 | 106.8 | 109.4
Ave. clustering coef. | 0.16 | 0.18 | 0.05
Eigenvector centrality | 0.052 | 0.053 | 0.056
Entropy | 9.6 | 9.8 | 10.1

Without setting thresholds for edges, we calculate the properties of the three
weighted networks. The statistics of the three networks are shown in table I.
As we can see, the MAPNet has a much smaller minimum node strength and maximum
node strength. The proposed method can significantly filter out edges with low
rate of presence across all of the layers. The ConNet and AveNet have higher
entropy since the two methods tend to conserve all the information of the
network as long as the edges are present in any of the datasets. On the
contrary, the proposed method tends to conserve the edges that have a high
rate of presence, and gives these edges higher weights as confirmed in Fig. 3
and Fig. 4.

###  3.2 Verification on unweighted multilayer network

In the second case, we implement the proposed method to an unweighted
multilayer network, and we do not distinguish the layers. In this case, the
edge weights in the aggregated network will be the rate of presence across the
layers. For example, if an edge is present in the first three layers, and
another edge is present in layers 3, 4, and 5, then we expect the weights of
the two edges are the same in the aggregated network.

The network we use to demonstrate the method is a friendship network that
consists of 5 layers. The nodes in the multilayer network are the employees of
the Department of Computer Science at Aarhus University, Denmark. The five
layers describe online and offline relationships between the employees, i.e.,
Facebook, Leisure, Work, Co-authorship, and Lunch. In Fig. 6, the blue curve
shows the CDF of the distribution of edge weights. The CDF of the
corresponding exponential distribution is the red curve in Fig. 6. As we can
see, there are five steps, which corresponding to the rate of presence of
edges in the five layers.

Figure 6: Edge weights of the aggregated network and the corresponding
exponential distribution.

##  4 Conclusion

In this paper, we proposed a _maximum a posteriori_ estimation-based algorithm
to aggregate mutlilayer networks. In constructing single-layer network from
various datasets, we first need to obtain a weighted network for each dataset.
The edge weights between node pairs are assumed to follow the Poisson
distribution across all the layers. Also, the edge weights in the target
network are thought to follow the exponential distribution. We apply the
expectation-maximization algorithm to find the optimized edge weights and rate
parameter of the exponential distribution.

In the weighted network example, the algorithm converges quickly. The ultimate
rate parameter and expected edge weights are affected by the number of
datasets. Both the coefficients of determination and rate parameters become
stable when there are at least 4 datasets. To evaluate the target network, we
adopted the Von Neumann entropy to measure the mixedness of the network. In
section 3, we show that the network obtained through our method conserves the
core information of the multilayer network. In the unweighted network example,
where the edges are not distinguishable, the edge weight in aggregated network
is exactly the rate of presence across all of the layers.

We compared the method to a network obtained through conflated datasets and a
network by averaging edge weights directly. The advantage of the proposed
method is that we are not required to normalize the expression values, and
thus the systemic errors that come from platforms can be reduced. However,
compared to the ConNet, the drawback of the proposed method is that there are
must be enough conditions in each dataset to find the edge weights in our
example. The experiments that contain less than 10 conditions are not taken
into consideration, which is the drawback of using a multilayer network.

## Acknowledgments

This works was supported by the National Institutes of Health under Grant No.
1R01AI140760-01A1.

## References

  * [1] Q. Yang, C. Yi, A. Vajdi, L.W. Cohnstaedt, H. Wu, X. Guo and C.M. Scoglio, ”Short-term forecasts and long-term mitigation evaluations for the COVID-19 epidemic in Hubei Province, China,” _Infectious Disease Modelling_ , vol. 5, pp. 563-574, 2020, doi.org/10.1016/j.idm.2020.08.001
  * [2] Q. Yang and H. Lu and C. Scoglio and M. de Jong and D. Gruenbacher, ”A network-based approach to model the development of city branding in China,” _Journal of Computers, Environment and Urban Systems_ , vol. 72, pp. 161-176, Nov. 2018, doi.org/10.1016/j.compenvurbsys.2018.06.007
  * [3] M. Newman, ”Network structure from rich but noisy data,” _Nature Physics_ , vol. 14, pp. 14, 2018, doi.org/10.1038/s41567-018-0076-1
  * [4] N. Eagle, A.S. Pentland, ”Reality mining: Sensing complex social systems,” _Personal Ubiquitous Comput._ , vol. 98, no. 36, pp. 255-268, May. 2006, doi.org/10.1007/s00779-005-0046-3
  * [5] N. Eagle, A.S. Pentland and D. Lazer, ”Inferring friendship network structure by using mobile phone data,” _Proceedings of the national academy of sciences_ , vol. 106, no. 36, pp. 15274-15278, Sept. 2009, doi.org/10.1073/pnas.0900282106
  * [6] M. Matteo, B. Micenkova and L. Rossi, ”Combinatorial analysis of multiple networks,” _arXiv preprint_ , 2013, arXiv:1303.4986.
  * [7] R.M. MacCallum, S.N. Redmond and G.K. Christophides, ”An expression map for Anopheles gambiae,” _BMC genomics_ , vol. 12, no. 1, pp. 1-16, Dec. 2011, doi.org/10.1186/1471-2164-12-620
  * [8] P. Adler, R. Kolde, M. Kull, A. Tkachenko, H. Peterson, J. Reimand, and J. Vilo, ”Mining for coexpression across hundreds of datasets using novel rank aggregation and visualization methods,” _Genome biology_ , vol. 10, pp. 1-11, Dec. 2009, doi.org/10.1186/gb-2009-10-12-r139
  * [9] H.K. Lee, A.K. Hsu, J. Sajdak, J. Qin, and P. Pavlidis, ”Coexpression analysis of human genes across many microarray data sets,” _BMC genomics_ , vol. 14, no. 6, pp. 1-16, Jun. 2004, doi: 10.1101/gr.1910904.
  * [10] J. Kuang and C. Scoglio, ”Layer reconstruction and missing link prediction of multilayer network with Maximum A Posteriori estimation,” _arXiv preprint_ , 2021, arXiv:2101.02733
  * [11] R. Interdonatoa, M. Magnanib, D. Pernac, A. Tagarellic, and D. Vega, ”Multilayer network simplification: Approaches, models and methods,” _Computer Science Review_ , vol. 36, pp. 100246, May. 2020, doi: 10.1016/j.cosrev.2020.100246.
  * [12] P. Bródka, K. Skibicki, P. Kazienko and K. Musiał, ”A degree centrality in multi-layered social network,” _International Conference on Computational Aspects of Social Networks (CASoN)_ , pp. 237-242, 2011, doi: 10.1109/CASON.2011.6085951.
  * [13] T. Chakraborty and R. Narayanam, ”Cross-layer betweenness centrality in multiplex networks with applications,” _IEEE 32nd International Conference on Data Engineering (ICDE’16)_ , pp. 397-408, 2016, doi: 10.1109/ICDE.2016.7498257.
  * [14] M. De Domenico, V. Nicosia, A. Arenas, and V. Latora, ”Structural reducibility of multilayer networks,” _Nature Communications_ , vol. 6, pp. 6864, Apr. 2015, doi: 10.1038/ncomms7864.
  * [15] W. Liu, P. Chen, S. Yeung, T. Suzumura and L. Chen, ”Principled multilayer network embedding,” _IEEE International Conference on Data Mining Workshops (ICDMW’17)_ , pp. 134-141, Nov. 2017, doi: 10.1109/ICDMW.2017.23.
  * [16] A. Grover and J. Leskovec, ”Structural reducibility of multilayer networks,” _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD’16)_ , pp. 855-864, Aug. 2016, doi:10.1145/2939672.2939754.
  * [17] K. Baltakys, J. Kanniainen, and F. Emmert-Streib, ”Multilayer aggregation with statistical validation: Application to investor networks,” _Scientific reports_ , vol. 8, pp. 8198, May. 2018, doi:10.1038/s41598-018-26575-2.
  * [18] M. Newman, ”Estimating network structure from unreliable measurements,” _Physical Review E_ , vol. 98, no. 6, pp. 062321, Dec. 2018, doi:10.1103/PhysRevE.98.062321.
  * [19] J. Young, G.T. Cantwell, and M. Newman, ”Robust Bayesian inference of network structure from unreliable data,” _arXiv preprint_ , 2020, arXiv:2008.03334.
  * [20] J. Young, F.S. Valdovinos, and M. Newman, ”Reconstruction of plant–pollinator networks from observational data,” _bioRxiv_ , 2019, doi: https://doi.org/10.1101/754077.
  * [21] V. Amati, A. Mol, T. Shafie, C. Hofman, and U. Brandes, ”A Framework for Reconstructing Archaeological Networks Using Exponential Random Graph Models,” _J Archaeol. Method Theory_ , vol. 27, pp. 192–219, 2020, doi:10.1007/s10816-019-09423-z.
  * [22] A.C. Koutsos, C. Blass, S. Meister, S. Schmidt, R.M. MacCallum, M.B. Soares, F.H. Collins, V. Benes, E. Zdobnov, F.C. Kafatos, and G.K. Christophides, ”Life cycle transcriptome of the malaria mosquito Anopheles gambiae and comparison with the fruitfly Drosophila melanogaster,” _Proceedings of the National Academy of Sciences_ , vol. 104, no. 27, pp. 11304-11309, July. 2007, doi: 10.1073/pnas.0703988104.
  * [23] O. Marinotti, E. Calvo, Q.K. Nguyen, S. Dissanayake, J.M.C. Ribeiro and A.A. James, ”Genome‐wide analysis of gene expression in adult Anopheles gambiae,” _Insect molecular biology_ , vol. 15, no. 1, pp. 1-12, Feb. 2007, doi.org/10.1111/j.1365-2583.2006.00610.x.
  * [24] B.J. Cassone, K. Mouline, M.W. Hahn, B.J. White, M. Pombi, F. Simard, C. Costantini, and N.J. Besansky, ”Differential gene expression in incipient species of Anopheles gambiae,” _Molecular ecology_ , vol. 17, no. 10, pp. 2491-2504, May. 2008, doi: 10.1111/j.1365-294X.2008.03774.x.
  * [25] Y. Goltsev, G.L. Rezende, K. Vranizan, G. Lanzaro, D. Valle, and M. Levine, ”Developmental and evolutionary basis for drought tolerance of the Anopheles gambiae embryo,” _Developmental biology_ , vol. 330, no. 2, pp. 462-470, Jun. 2009, doi: 10.1016/j.ydbio.2009.02.038.
  * [26] A.M. Mendes, P.H. Awono-Ambene, S.E. Nsango, A. Cohuet, D. Fontenille, F.C. Kafatos, G.K. Christophides, I. Morlais, and D. Vlachou, ”Infection intensity-dependent responses of Anopheles gambiae to the African malaria parasite Plasmodium falciparum,” _Infection and immunity_ , vol. 79, no. 11, pp. 4708-4715, Nov. 2011, doi:10.1128/IAI.05647-11.
  * [27] B.J. Cassone, M.J. Molloy, C. Cheng, J.C. Tan, M.W. Hahn, and N.J. Besansky, ”Divergent transcriptional response to thermal stress by Anopheles gambiae larvae carrying alternative arrangements of inversion 2La,” _Molecular ecology_ , vol. 20, no. 12, pp. 2567-2580, Nov. 2011, doi:10.1111/j.1365-294X.2011.05114.x.
  * [28] D.A. Baker, T. Nolan, B. Fischer, A. Pinder, A. Crisanti, and S. Russell , ”A comprehensive gene expression atlas of sex-and tissue-specificity in the malaria vector, Anopheles gambiae,” _BMC Genomics_ , vol. 12, pp. 296, Jun. 2011, doi:10.1186/1471-2164-12-296.
  * [29] S.S.C. Rund, T.Y. Hou, S.M. Ward, F.H. Collins, and G.E. Duffield, ”Genome-wide profiling of diel and circadian gene expression in the malaria vector Anopheles gambiae,” _Proceedings of the National Academy of Sciences_ , vol. 108, no. 32, pp. E421-E430, Aug. 2011, doi:10.1073/pnas.1100584108.