Skip to content
This repository was archived by the owner on Jan 27, 2023. It is now read-only.

Commit a4d3e72

Browse files
author
Antonio Ulloa
committed
updated simulation scripts
1 parent da6b4d4 commit a4d3e72

11 files changed

+421
-70
lines changed

‎analysis/avg_FC_TVB_ROIs_across_subjects.py‎

Lines changed: 169 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -86,6 +86,8 @@
8686

8787
from scipy.stats import t
8888

89+
from scipy.stats import itemfreq
90+
8991
from scipy import stats
9092

9193
import math as m
@@ -94,6 +96,8 @@
9496

9597
from scipy.stats import skew
9698

99+
import scipy.io
100+
97101
from matplotlib import cm as CM
98102

99103
from mne.viz import circular_layout, plot_connectivity_circle
@@ -130,77 +134,81 @@ def adj_to_list(input_filename,output_filename,delimiter):
130134
# construct array of subjects to be considered
131135
subjects = np.arange(10)
132136

137+
# define name of input file where Hagmann empirical data is stored (matlab file given to us
138+
# by Olaf Sporns and Chris Honey
139+
hagmann_data = 'DSI_release2_2011.mat'
140+
133141
# define output file where means, standard deviations, and variances will be stored
134142
RS_FC_avg_file = 'tvb_only_rs_fc_avg.npy'
135143
TB_FC_avg_file = 'tb_fc_avg.npy'
136144

137145
# declare ROI labels
138-
labels = ['rLOF',
139-
'rPORB',
140-
'rFP' ,
141-
'rMOF' ,
142-
'rPTRI',
143-
'rPOPE',
144-
'rRMF' ,
145-
'rSF' ,
146-
'rCMF' ,
147-
'rPREC',
148-
'rPARC',
149-
'rRAC' ,
150-
'rCAC' ,
151-
'rPC' ,
152-
'rISTC',
153-
'rPSTC',
154-
'rSMAR',
155-
'rSP' ,
156-
'rIP' ,
157-
'rPCUN',
158-
'rCUN' ,
159-
'rPCAL',
160-
'rLOCC',
161-
'rLING',
162-
'rFUS' ,
163-
'rPARH',
164-
'rENT' ,
165-
'rTP' ,
166-
'rIT' ,
167-
'rMT' ,
168-
'rBSTS',
169-
'rST' ,
170-
'rTT',
171-
'lLOF' ,
172-
'lPORB',
173-
'lFP' ,
174-
'lMOF' ,
175-
'lPTRI',
176-
'lPOPE',
177-
'LRMF' ,
178-
'lSF' ,
179-
'lCMF' ,
180-
'lPREC',
181-
'lPARC',
182-
'lRAC' ,
183-
'lCAC' ,
184-
'lPC' ,
185-
'lISTC',
186-
'lPSTC',
187-
'lSMAR',
188-
'lSP' ,
189-
'lIP' ,
190-
'lPCUN',
191-
'lCUN' ,
192-
'lPCAL',
193-
'lLOC' ,
194-
'lLING',
195-
'lFUS' ,
196-
'lPARH',
197-
'lENT' ,
198-
'lTP' ,
199-
'lIT' ,
200-
'lMT' ,
201-
'lBSTS',
202-
'lST' ,
203-
'lTT'
146+
labels = [' rLOF',
147+
'rPORB',
148+
' rFP',
149+
' rMOF',
150+
'rPTRI',
151+
'rPOPE',
152+
' rRMF',
153+
' rSF',
154+
' rCMF',
155+
'rPREC',
156+
'rPARC',
157+
' rRAC',
158+
' rCAC',
159+
' rPC',
160+
'rISTC',
161+
'rPSTC',
162+
'rSMAR',
163+
' rSP',
164+
' rIP',
165+
'rPCUN',
166+
' rCUN',
167+
'rPCAL',
168+
'rLOCC',
169+
'rLING',
170+
' rFUS',
171+
'rPARH',
172+
' rENT',
173+
' rTP',
174+
' rIT',
175+
' rMT',
176+
'rBSTS',
177+
' rST',
178+
' rTT',
179+
' lLOF',
180+
'lPORB',
181+
' lFP',
182+
' lMOF',
183+
'lPTRI',
184+
'lPOPE',
185+
' lRMF',
186+
' lSF',
187+
' lCMF',
188+
'lPREC',
189+
'lPARC',
190+
' lRAC',
191+
' lCAC',
192+
' lPC',
193+
'lISTC',
194+
'lPSTC',
195+
'lSMAR',
196+
' lSP',
197+
' lIP',
198+
'lPCUN',
199+
' lCUN',
200+
'lPCAL',
201+
'lLOCC',
202+
'lLING',
203+
' lFUS',
204+
'lPARH',
205+
' lENT',
206+
' lTP',
207+
' lIT',
208+
' lMT',
209+
'lBSTS',
210+
' lST',
211+
' lTT'
204212
]
205213

206214

@@ -249,6 +257,8 @@ def adj_to_list(input_filename,output_filename,delimiter):
249257
TVB_LSNM_DMS_subj9 = 'subject_19/output.DMSTask/xcorr_matrix_66_regions.npy'
250258
TVB_LSNM_DMS_subj10 = 'subject_20/output.DMSTask/xcorr_matrix_66_regions.npy'
251259

260+
# open matlab file that contains hagmann empirical data
261+
hagmann_empirical = scipy.io.loadmat(hagmann_data)
252262

253263
# open files that contain correlation coefficients
254264
tvb_rs_subj1 = np.load(TVB_RS_subj1)
@@ -408,7 +418,7 @@ def adj_to_list(input_filename,output_filename,delimiter):
408418
# plot the correlation coefficients for rPC
409419
fig = plt.figure('Correlation coefficients using rPC as seed')
410420
ax = fig.add_subplot(111)
411-
rPC_loc = labels.index('rPC')
421+
rPC_loc = labels.index(' rPC')
412422
y_pos = np.arange(len(labels))
413423
ax.barh(y_pos, tvb_rs_mean[rPC_loc])
414424
ax.set_yticks(y_pos)
@@ -1063,5 +1073,97 @@ def adj_to_list(input_filename,output_filename,delimiter):
10631073
plt.text(1, 22, 'r=' + '{:.2f}'.format(r))
10641074

10651075

1076+
###################################################################################################
1077+
# plot correlation coefficients of the simulated ROIs vs empirical ROIs
1078+
# first, find average of correlation coefficients within each lo-res ROI by averaging across hi-res
1079+
# ROIs located within each one of the 66 lo-res ROIs
1080+
##################################################################################################
1081+
1082+
# we need to apply a Fisher Z transformation to the correlation coefficients,
1083+
# prior to averaging.
1084+
empirical_fc_mat_Z = np.arctanh(hagmann_empirical['COR_fMRI_average'])
1085+
1086+
# initialize 66x66 matrix
1087+
empirical_fc_lowres_Z = np.zeros([66,66])
1088+
1089+
# subtract 1 from labels array bc numpy arrays start with zero
1090+
hagmann_empirical['roi_lbls'][0] = hagmann_empirical['roi_lbls'][0] - 1
1091+
1092+
# compress 998x998 FC matrix to 66x66 FC matrix by averaging
1093+
for i in range(0,998):
1094+
for j in range(0,998):
1095+
1096+
# extract low-res coordinates from hi-res labels matrix
1097+
x = hagmann_empirical['roi_lbls'][0][i]
1098+
y = hagmann_empirical['roi_lbls'][0][j]
1099+
1100+
empirical_fc_lowres_Z[x, y] += empirical_fc_mat_Z[i, j]
1101+
1102+
# count the number of times each lowres label appears in the hires matrix
1103+
freq_array = itemfreq(hagmann_empirical['roi_lbls'][0])
1104+
1105+
# divide each sum by the number of hires ROIs within each lowres ROI
1106+
for i in range(0,66):
1107+
for j in range(0,66):
1108+
total_freq = freq_array[i][1] * freq_array[j][1]
1109+
empirical_fc_lowres_Z[i,j] = empirical_fc_lowres_Z[i,j] / total_freq
1110+
1111+
# now, convert back to from Z to R correlation coefficients
1112+
empirical_fc_lowres = np.tanh(empirical_fc_lowres_Z)
1113+
1114+
# initialize figure to plot simulated FC
1115+
fig=plt.figure('Functional Connectivity Matrix of empirical BOLD (66 ROIs)')
1116+
ax = fig.add_subplot(111)
1117+
cmap = CM.get_cmap('jet', 10)
1118+
empirical_fc_lowres = np.asarray(empirical_fc_lowres)
1119+
cax = ax.imshow(empirical_fc_lowres, vmin=-0.4, vmax=1.0, interpolation='nearest', cmap=cmap)
1120+
ax.grid(False)
1121+
color_bar=plt.colorbar(cax)
1122+
1123+
# simulated RS FC needs to be rearranged prior to scatter plotting
1124+
new_TVB_RS_FC = np.zeros([66, 66])
1125+
for i in range(0, 66):
1126+
for j in range(0, 66):
1127+
1128+
# extract labels of current ROI label from simulated labels list of simulated FC
1129+
label_i = labels[i]
1130+
label_j = labels[j]
1131+
# extract index of corresponding ROI label from empirical FC labels list
1132+
emp_i = np.where(hagmann_empirical['anat_lbls'] == label_i)[0][0]
1133+
emp_j = np.where(hagmann_empirical['anat_lbls'] == label_j)[0][0]
1134+
new_TVB_RS_FC[emp_i, emp_j] = tvb_rs_mean[i, j]
1135+
1136+
# initialize figure to plot correlations between empirical and simulated FC
1137+
fig=plt.figure('Empirical vs Simulated FC')
1138+
1139+
# apply mask to get rid of upper triangle, including main diagonal
1140+
mask = np.tri(new_TVB_RS_FC.shape[0], k=0)
1141+
mask = np.transpose(mask)
1142+
new_TVB_RS_FC = np.ma.array(new_TVB_RS_FC, mask=mask) # mask out upper triangle
1143+
empirical_fc_lowres = np.ma.array(empirical_fc_lowres, mask=mask) # mask out upper triangle
1144+
1145+
# flatten the numpy cross-correlation matrix
1146+
corr_mat_sim_TVB_RS_FC = np.ma.ravel(new_TVB_RS_FC)
1147+
corr_mat_emp_FC = np.ma.ravel(empirical_fc_lowres)
1148+
1149+
# remove masked elements from cross-correlation matrix
1150+
corr_mat_sim_TVB_RS_FC = np.ma.compressed(corr_mat_sim_TVB_RS_FC)
1151+
corr_mat_emp_FC = np.ma.compressed(corr_mat_emp_FC)
1152+
1153+
# scatter plot
1154+
plt.scatter(corr_mat_emp_FC, corr_mat_sim_TVB_RS_FC)
1155+
plt.xlabel('Empirical FC')
1156+
plt.ylabel('Model FC')
1157+
1158+
print corr_mat_emp_FC, corr_mat_sim_TVB_RS_FC
1159+
1160+
# fit scatter plot with np.polyfit
1161+
m, b = np.polyfit(corr_mat_emp_FC, corr_mat_sim_TVB_RS_FC, 1)
1162+
plt.plot(corr_mat_emp_FC, m*corr_mat_emp_FC + b, '-', color='red')
1163+
1164+
# calculate correlation coefficient and display it on plot
1165+
r = np.corrcoef(corr_mat_emp_FC, corr_mat_sim_TVB_RS_FC)[1,0]
1166+
plt.text(0.5, 0.3, 'r=' + '{:.2f}'.format(r))
1167+
10661168
# Show the plots on the screen
10671169
plt.show()

‎analysis/compute_BOLD_balloon_66_regions.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -76,14 +76,14 @@
7676

7777

7878
# define the name of the input file where the synaptic activities are stored
79-
SYN_file = 'synaptic_in_66_ROIs.npy'
79+
SYN_file = 'synaptic_in_66_ROIs_incl_LSNM.npy'
8080

8181
# define the name of the output file where the BOLD timeseries will be stored
82-
BOLD_file = 'bold_balloon_66_regions.npy'
82+
BOLD_file = 'bold_balloon_66_regions_incl_LSNM.npy'
8383

8484
# define the name of the output file where the cross-correlation matrix will
8585
# be stored
86-
xcorr_file = 'xcorr_matrix_66_regions.npy'
86+
xcorr_file = 'xcorr_matrix_66_regions_incl_LSNM.npy'
8787

8888
# define balloon model parameters...
8989
tau_s = 1.5 # rate constant of vasodilatory signal decay in seconds

0 commit comments

Comments
 (0)