-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathHMA.py
210 lines (156 loc) · 6.92 KB
/
HMA.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
# -*- coding: utf-8 -*-
"""
Created on Tue May 11 19:00:58 2021
Hierarchical modular partition of FC networks, Integration and Segregation related metrics
using eigenmode analysis [1,2]. The codes were adapted from
https://github.com/TobousRong/Hierarchical-module-analysis
[1] Wang, R., Lin, P., Liu, M., Wu, Y., Zhou, T., & Zhou, C. (2019).
Hierarchical connectome modes and critical state jointly maximize
human brain functional diversity. Physical review letters, 123(3),
038301.
[2] Wang, R., Liu, M., Cheng, X., Wu, Y., Hildebrandt, A., & Zhou, C. (2021).
Segregation, integration, and balance of large-scale resting brain networks
configure different cognitive abilities. Proceedings of the National
Academy of Sciences, 118(23).
[3] Wang, R., Fan, Y., Wu, Y., & Zhou, C. (2021). Heterogeneous aging trajectories
within resting-state brain networks predict distinct ADHD symptoms in adults.
arXiv preprint arXiv:2107.13219.
@author: Carlos Coronel
"""
import numpy as np
def Functional_HP(FC):
'''
This function uses an eigenmode-based analysis to detect the hierarchical
modules in FC networks.
Parameters
----------
FC : numpy array.
functional connectivity matrix.
Returns
-------
Clus_num : list.
number of modules found at each eigenmode level.
Clus_size : list of arrays.
number of nodes belonging to each module at each eigenmode level.
H_all : nested list of arrays
it contains all the assignments to each module within hierarchies.
'''
N = FC.shape[0] #number of ROIs
H_all= [] #all modules assignments
#This method requires the complete positive connectivity in FC matrix,
#that generates the global integration in the first level.
FC[FC < 0] = 0
FC = (FC + FC.T) / 2
#singular value decomposition, where 'u' and 'v' are the eigenvectors, and
#'s' the singular values (eigenvalues).
u, s, v = np.linalg.svd(FC)
s[s<0] = 0 #For avoiding negative eigenvalues
#Number of nodes in each module at 1st level.
H1_1 = np.argwhere(u[:,1] < 0)[:,0]
H1_2 = np.argwhere(u[:,1] >= 0)[:,0]
H1 = [] #modules of mode 1
H1.append(H1_1)
H1.append(H1_2)
H_all.append(H1)
Clus_num= [1] #The first level has one module and corresponds to the global integration.
Clus_size = [[N]] #The first level has one module with N total number of ROIs (nodes).
for mode in range(1,N-1):
exec('H%i = []'%(mode + 1))
x = np.argwhere(u[:,mode + 1] >= 0)[:,0]
y = np.argwhere(u[:,mode + 1] < 0)[:,0]
H = []
for j in range(0, 2 * Clus_num[mode-1]):
#assume the number of cluster in j-1 level is 2^(mode-1)
H.append(eval('H' + '%s_%s'%(mode,j+1)))
idx = np.array([len(index) for index in H]) #length of each cluster in H
#Delete the cluster with 0 size
H = [H[full] for full in range(0,len(idx)) if idx[full] != 0]
idx = [idx[full] for full in range(0,len(idx)) if idx[full] != 0]
Clus_size.append(idx)
Clus_num.append(len(H)) #number of cluster
k = 0
for j in range(0, 2 * Clus_num[mode], 2): #modular number
Positive_Node = np.intersect1d(H[k], x)
Negative_Node = np.intersect1d(H[k], y)
k = k + 1
exec('H' + '%s_%s'%((mode + 1), (j + 2)) + '=' + 'Positive_Node')
exec('H' + '%s_%s'%((mode + 1), (j + 1)) + '=' + 'Negative_Node')
exec('H%i.append(H%s_%s)'%((mode + 1),(mode + 1), (j + 2)))
exec('H%i.append(H%s_%s)'%((mode + 1),(mode + 1), (j + 1)))
exec('H_all.append(H%i)'%(mode + 1))
return([Clus_num,Clus_size,H_all])
def Balance(FC, Clus_num, Clus_size):
'''
This function calculates the integration and segregation components.
Parameters
----------
FC : numpy array.
functional connectivity matrix.
Clus_num : list.
number of modules found at each eigenmode level.
Clus_size: list of arrays.
number of nodes belonging to each module at each eigenmode level.
Returns
-------
Hin : float.
Integration component.
Hse : float.
Segregation component.
'''
N = FC.shape[0] #number of ROIs
#This method requires the complete positive connectivity in FC matrix,
#that generates the global integration in the first level.
FC[FC < 0] = 0
FC = (FC + FC.T) / 2
#singular value decomposition, where 'u' and 'v' are the eigenvectors, and
#'s' the singular values (eigenvalues).
u, s, v = np.linalg.svd(FC)
s[s<0] = 0
s = s ** 2 #using the squared Lambda
p = np.zeros(N-1)
#modular size correction
for i in range(0,len(Clus_num) - 1):
p[i] = np.sum(np.abs(np.array(Clus_size[i]) - N / Clus_num[i])) / N
HF = s[0:(N-1)] * np.array(Clus_num) *(1-p);
Hin = np.sum(HF[0]) / N**2 #integration component
Hse = np.sum(HF[1:(N-1)]) / N**2 #segregation component
return([Hin,Hse])
def nodal_measures(FC, Clus_num, Clus_size):
'''
This function calculates the nodal (regional) integration and segregation components.
Parameters
----------
FC : numpy array.
functional connectivity matrix.
Clus_num : list.
number of modules found at each eigenmode level.
Clus_size: list of arrays.
number of nodes belonging to each module at each eigenmode level.
Returns
-------
Hin_nodal : list, float.
Integration component of the N total number of nodes.
Hse_nodal : list, float.
Segregation component of the N total number of nodes.
'''
N = FC.shape[0] #number of ROIs
#This method requires the complete positive connectivity in FC matrix,
#that generates the global integration in the first level.
FC[FC < 0] = 0
FC = (FC + FC.T) / 2
#singular value decomposition, where 'u' and 'v' are the eigenvectors, and
#'s' the singular values (eigenvalues).
u, s, v = np.linalg.svd(FC)
s[s<0] = 0
s = s ** 2 #using the squared Lambda
p = np.zeros(N-1)
#modular size correction
for i in range(0,len(Clus_num) - 1):
p[i] = np.sum(np.abs(np.array(Clus_size[i]) - N / Clus_num[i])) / N
HF = s[0:(N-1)] * np.array(Clus_num) *(1-p)
#Nodal (regional) integration and segregation components
Hin_nodal = HF[0] / N * u[:,0]**2
Hse_nodal = np.zeros(N)
for i in range(1,N-1):
Hse_nodal += HF[i] / N * u[:,i]**2
return([Hin_nodal,Hse_nodal])