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

Commit 00dc2a1

Browse files
author
Antonio Ulloa
committed
Debugged analysis and visualization scripts
1 parent 50150c4 commit 00dc2a1

8 files changed

+389
-112
lines changed

analysis/average_BOLD_across_subjects.py

Lines changed: 20 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -87,26 +87,26 @@
8787

8888
# define the name of the input files where the BOLD and synaptic timeseries are
8989
# stored:
90-
syn_subj = ['subject_11/output.36trials/synaptic_in_ROI.npy',
91-
'subject_12/output.36trials/synaptic_in_ROI.npy',
92-
'subject_13/output.36trials/synaptic_in_ROI.npy',
93-
'subject_14/output.36trials/synaptic_in_ROI.npy',
94-
'subject_15/output.36trials/synaptic_in_ROI.npy',
95-
'subject_16/output.36trials/synaptic_in_ROI.npy',
96-
'subject_17/output.36trials/synaptic_in_ROI.npy',
97-
'subject_18/output.36trials/synaptic_in_ROI.npy',
98-
'subject_19/output.36trials/synaptic_in_ROI.npy',
99-
'subject_20/output.36trials/synaptic_in_ROI.npy']
100-
BOLD_subj = ['subject_11/output.36trials/lsnm_bold_balloon.npy',
101-
'subject_12/output.36trials/lsnm_bold_balloon.npy',
102-
'subject_13/output.36trials/lsnm_bold_balloon.npy',
103-
'subject_14/output.36trials/lsnm_bold_balloon.npy',
104-
'subject_15/output.36trials/lsnm_bold_balloon.npy',
105-
'subject_16/output.36trials/lsnm_bold_balloon.npy',
106-
'subject_17/output.36trials/lsnm_bold_balloon.npy',
107-
'subject_18/output.36trials/lsnm_bold_balloon.npy',
108-
'subject_19/output.36trials/lsnm_bold_balloon.npy',
109-
'subject_20/output.36trials/lsnm_bold_balloon.npy']
90+
syn_subj = ['subject_11/output.36trials.with_feedback/synaptic_in_TVB_ROI.npy',
91+
'subject_12/output.36trials.with_feedback/synaptic_in_TVB_ROI.npy',
92+
'subject_13/output.36trials.with_feedback/synaptic_in_TVB_ROI.npy',
93+
'subject_14/output.36trials.with_feedback/synaptic_in_TVB_ROI.npy',
94+
'subject_15/output.36trials.with_feedback/synaptic_in_TVB_ROI.npy',
95+
'subject_16/output.36trials.with_feedback/synaptic_in_TVB_ROI.npy',
96+
'subject_17/output.36trials.with_feedback/synaptic_in_TVB_ROI.npy',
97+
'subject_18/output.36trials.with_feedback/synaptic_in_TVB_ROI.npy',
98+
'subject_19/output.36trials.with_feedback/synaptic_in_TVB_ROI.npy',
99+
'subject_20/output.36trials.with_feedback/synaptic_in_TVB_ROI.npy']
100+
BOLD_subj = ['subject_11/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
101+
'subject_12/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
102+
'subject_13/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
103+
'subject_14/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
104+
'subject_15/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
105+
'subject_16/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
106+
'subject_17/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
107+
'subject_18/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
108+
'subject_19/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
109+
'subject_20/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy']
110110

111111
# open files that contain synaptic and fMRI BOLD timeseries
112112
lsnm_syn = np.zeros((num_of_subjects, num_of_modules, experiment_length))

analysis/avg_func_conn_across_subjects.py

Lines changed: 17 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -305,13 +305,23 @@
305305

306306
# convert to Pandas dataframe, using the transpose to convert to a format where the names
307307
# of the modules are the labels for each time-series
308-
d_mean = pd.DataFrame(np.array([d_syn_mean,
309-
d_fmri_mean]),
310-
columns=np.array(['V1', 'V4', 'FS', 'D1', 'D2', 'FR', 'cIT']),
308+
d_mean = pd.DataFrame(np.array([d_syn_mean[:-1],
309+
d_fmri_mean[:-1]]),
310+
columns=np.array(['V1',
311+
'V4',
312+
'FS',
313+
'D1',
314+
'D2',
315+
'FR']), #, 'cIT']),
311316
index=np.array(['ISA', 'fMRI']))
312-
d_std = pd.DataFrame(np.array([d_syn_sem,
313-
d_fmri_sem]),
314-
columns=np.array(['V1', 'V4', 'FS', 'D1', 'D2', 'FR', 'cIT']),
317+
d_std = pd.DataFrame(np.array([d_syn_sem[:-1],
318+
d_fmri_sem[:-1]]),
319+
columns=np.array(['V1',
320+
'V4',
321+
'FS',
322+
'D1',
323+
'D2',
324+
'FR']), #, 'cIT']),
315325
index=np.array(['ISA', 'fMRI']))
316326

317327
# now, plot means and std's using 'pandas framework...
@@ -324,7 +334,7 @@
324334
ax = plt.gca() # get hold of the axes
325335

326336
bars=d_mean.plot(yerr=d_std, ax=ax, kind='bar',
327-
color=['yellow', 'green', 'orange', 'red', 'pink', 'purple', 'lightblue'],
337+
color=['yellow', 'green', 'orange', 'red', 'pink', 'purple'], #, 'lightblue'],
328338
ylim=[-0.5, 0.4])
329339

330340
# change the location of the legend

analysis/compute_BOLD_balloon_model.py

Lines changed: 11 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -76,10 +76,10 @@
7676
from scipy.integrate import odeint
7777

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

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

8484
# define balloon model parameters...
8585
tau_s = 1.5 # rate constant of vasodilatory signal decay in seconds
@@ -154,8 +154,8 @@ def balloon_function(y, t, syn):
154154
Ti = 0.005 * 10
155155

156156
# Total time of scanning experiment in seconds (timesteps X 5)
157-
#T = 198
158-
T = 110
157+
T = 198
158+
#T = 110
159159

160160
# Time for one complete trial in milliseconds
161161
Ttrial = 5.5
@@ -319,10 +319,12 @@ def balloon_function(y, t, syn):
319319
# Thus, we need to multiply the datapoint times 50 ms...
320320
# ... and divide by 1000 to convert to seconds
321321
#t = np.linspace(0, 659*50./1000., num=660)
322-
t = np.linspace(0, synaptic_timesteps * 50.0 / 1000., num=synaptic_timesteps)
322+
t = np.linspace(0, synaptic_timesteps+1 * 50.0 / 1000., num=synaptic_timesteps+1)
323323

324324
# downsample the BOLD signal arrays to produce scan rate of 2 per second
325-
scanning_timescale = np.arange(0, synaptic_timesteps, synaptic_timesteps / (T/Tr))
325+
scanning_timescale = np.linspace(0, synaptic_timesteps, num=T/Tr)
326+
scanning_timescale = np.round(scanning_timescale)
327+
scanning_timescale = scanning_timescale.astype(int)
326328
mr_time = t[scanning_timescale]
327329
v1_BOLD = v1_BOLD[scanning_timescale]
328330
v4_BOLD = v4_BOLD[scanning_timescale]
@@ -374,9 +376,9 @@ def balloon_function(y, t, syn):
374376
plt.suptitle('SIMULATED SYNAPTIC ACTIVITY')
375377

376378
# Plot synaptic activities
377-
plt.plot(t, v1_syn, linewidth=3.0, color='yellow')
378-
plt.plot(t, it_syn, linewidth=3.0, color='blue')
379-
plt.plot(t, d1_syn, linewidth=3.0, color='red')
379+
plt.plot(t, syn[0], linewidth=3.0, color='yellow')
380+
plt.plot(t, syn[1], linewidth=3.0, color='blue')
381+
plt.plot(t, syn[2], linewidth=3.0, color='red')
380382
plt.gca().set_axis_bgcolor('black')
381383

382384
# Set up separate figures to plot fMRI BOLD signal

analysis/compute_PSC.py

Lines changed: 11 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -101,22 +101,22 @@
101101
synaptic_timesteps = experiment_length - synaptic_steps_removed
102102

103103
# define the names of the input files where the BOLD timeseries are contained
104-
BOLD_ts_subj=np.array(['subject_11/output.36trials.with_feedback/lsnm_bold_balloon.npy',
105-
'subject_12/output.36trials.with_feedback/lsnm_bold_balloon.npy',
106-
'subject_13/output.36trials.with_feedback/lsnm_bold_balloon.npy',
107-
'subject_14/output.36trials.with_feedback/lsnm_bold_balloon.npy',
108-
'subject_15/output.36trials.with_feedback/lsnm_bold_balloon.npy',
109-
'subject_16/output.36trials.with_feedback/lsnm_bold_balloon.npy',
110-
'subject_17/output.36trials.with_feedback/lsnm_bold_balloon.npy',
111-
'subject_18/output.36trials.with_feedback/lsnm_bold_balloon.npy',
112-
'subject_19/output.36trials.with_feedback/lsnm_bold_balloon.npy',
113-
'subject_20/output.36trials.with_feedback/lsnm_bold_balloon.npy'])
104+
BOLD_ts_subj=np.array(['subject_11/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
105+
'subject_12/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
106+
'subject_13/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
107+
'subject_14/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
108+
'subject_15/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
109+
'subject_16/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
110+
'subject_17/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
111+
'subject_18/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
112+
'subject_19/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy',
113+
'subject_20/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy'])
114114

115115
# set matplot lib parameters to produce visually appealing plots
116116
#mpl.style.use('ggplot')
117117

118118
# define output file where means, standard deviations, and variances will be stored
119-
PSC_stats_FILE = 'psc_stats.txt'
119+
PSC_stats_FILE = 'psc_stats_TVB_ROI.txt'
120120

121121
# define neural synaptic time interval in seconds. The simulation data is collected
122122
# one data point at synaptic intervals (10 simulation timesteps). Every simulation

analysis/compute_syn_visual_TVB.py

Lines changed: 198 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,198 @@
1+
# ============================================================================
2+
#
3+
# PUBLIC DOMAIN NOTICE
4+
#
5+
# National Institute on Deafness and Other Communication Disorders
6+
#
7+
# This software/database is a "United States Government Work" under the
8+
# terms of the United States Copyright Act. It was written as part of
9+
# the author's official duties as a United States Government employee and
10+
# thus cannot be copyrighted. This software/database is freely available
11+
# to the public for use. The NIDCD and the U.S. Government have not placed
12+
# any restriction on its use or reproduction.
13+
#
14+
# Although all reasonable efforts have been taken to ensure the accuracy
15+
# and reliability of the software and data, the NIDCD and the U.S. Government
16+
# do not and cannot warrant the performance or results that may be obtained
17+
# by using this software or data. The NIDCD and the U.S. Government disclaim
18+
# all warranties, express or implied, including warranties of performance,
19+
# merchantability or fitness for any particular purpose.
20+
#
21+
# Please cite the author in any work or product based on this material.
22+
#
23+
# ==========================================================================
24+
25+
26+
27+
# ***************************************************************************
28+
#
29+
# Large-Scale Neural Modeling software (LSNM)
30+
#
31+
# Section on Brain Imaging and Modeling
32+
# Voice, Speech and Language Branch
33+
# National Institute on Deafness and Other Communication Disorders
34+
# National Institutes of Health
35+
#
36+
# This file (compute_syn_visual_TVB.py) was created on June 4, 2016.
37+
#
38+
#
39+
# Author: Antonio Ulloa
40+
#
41+
# Last updated by Antonio Ulloa on June 4, 2016
42+
#
43+
# Based on computer code originally developed by Barry Horwitz et al
44+
# **************************************************************************/
45+
46+
# compute_syn_visual_TVB.py
47+
#
48+
# Calculate and plot simulated synaptic activities in given ROIs that include
49+
# only TVB connectome nodes
50+
#
51+
# ... using data from visual delay-match-to-sample simulation.
52+
# It also saves the synaptic activities for each and all modules in a python data file
53+
# (*.npy)
54+
# The data is saved in a numpy array with columns in the following order:
55+
#
56+
# V1 ROI (right hemisphere, contains only TVB nodes)
57+
# V4 ROI (right hemisphere, contains only TVB nodes)
58+
# IT ROI (right hemisphere, contains only TVB nodes)
59+
# FS ROI (right hemisphere, contains only TVB nodes)
60+
# D1 ROI (right hemisphere, contains only TVB nodes)
61+
# D2 ROI (right hemisphere, contains only TVB nodes)
62+
# FR ROI (right hemisphere, contains only TVB nodes)
63+
# IT ROI (left hemisphere, contains only TVB nodes)
64+
65+
import numpy as np
66+
67+
import matplotlib.pyplot as plt
68+
69+
import matplotlib as mpl
70+
71+
# set matplot lib parameters to produce visually appealing plots
72+
mpl.style.use('ggplot')
73+
74+
# define the name of the output file where the integrated synaptic activity will be stored
75+
syn_file = 'synaptic_in_TVB_ROI.npy'
76+
77+
# the following ranges define the location of the nodes within a given ROI in Hagmann's brain.
78+
# They were taken from the excel document:
79+
# "Hagmann's Talairach Coordinates (obtained from TVB).xlsx"
80+
# Extracted from The Virtual Brain Demo Data Sets
81+
# Please note that arrays in Python start from zero so one does need to account for that and shift
82+
# indices given by the above document by one location.
83+
# Use 6 nodes within rPCAL including host node 345
84+
v1_loc = range(344, 350) # Hagmann's brain nodes included within V1 ROI
85+
86+
# Use 6 nodes within rFUS including host node 393
87+
v4_loc = range(390, 396) # Hagmann's brain nodes included within V4 ROI
88+
89+
# Use 6 nodes within rPARH including host node 413
90+
it_loc = range(412, 418) # Hagmann's brain nodes included within IT ROI
91+
92+
# Use 6 nodes within rRMF including host node 74
93+
d1_loc = range(73, 79) # Hagmann's brain nodes included within D1 ROI
94+
95+
# Use 6 nodes within rPTRI including host node 41
96+
d2_loc = range(39, 45) # Hagmann's brain nodes included within D2 ROI
97+
98+
# Use 6 nodes within rPOPE including host node 47
99+
fs_loc = range(47, 53) # Hagmann's brain nodes included within FS ROI
100+
101+
# Use 6 nodes within rCMF including host node 125
102+
fr_loc = range(125, 131) # Hagmann's brain nodes included within FR ROI
103+
104+
# Use 6 nodes within lPARH
105+
lit_loc= range(911, 917) # Hagmann's brain nodes included within left IT ROI
106+
107+
# Load TVB nodes synaptic activity
108+
tvb_synaptic = np.load("tvb_synaptic.npy")
109+
110+
# Load TVB host node synaptic activities into separate numpy arrays
111+
tvb_ev1 = tvb_synaptic[:, 0, v1_loc[0]:v1_loc[-1]+1, 0]
112+
tvb_ev4 = tvb_synaptic[:, 0, v4_loc[0]:v4_loc[-1]+1, 0]
113+
tvb_eit = tvb_synaptic[:, 0, it_loc[0]:it_loc[-1]+1, 0]
114+
tvb_ed1 = tvb_synaptic[:, 0, d1_loc[0]:d1_loc[-1]+1, 0]
115+
tvb_ed2 = tvb_synaptic[:, 0, d2_loc[0]:d2_loc[-1]+1, 0]
116+
tvb_efs = tvb_synaptic[:, 0, fs_loc[0]:fs_loc[-1]+1, 0]
117+
tvb_efr = tvb_synaptic[:, 0, fr_loc[0]:fr_loc[-1]+1, 0]
118+
tvb_iv1 = tvb_synaptic[:, 1, v1_loc[0]:v1_loc[-1]+1, 0]
119+
tvb_iv4 = tvb_synaptic[:, 1, v4_loc[0]:v4_loc[-1]+1, 0]
120+
tvb_iit = tvb_synaptic[:, 1, it_loc[0]:it_loc[-1]+1, 0]
121+
tvb_id1 = tvb_synaptic[:, 1, d1_loc[0]:d1_loc[-1]+1, 0]
122+
tvb_id2 = tvb_synaptic[:, 1, d2_loc[0]:d2_loc[-1]+1, 0]
123+
tvb_ifs = tvb_synaptic[:, 1, fs_loc[0]:fs_loc[-1]+1, 0]
124+
tvb_ifr = tvb_synaptic[:, 1, fr_loc[0]:fr_loc[-1]+1, 0]
125+
126+
# now extract synaptic activity in the contralateral IT
127+
tvb_elit = tvb_synaptic[:, 0, lit_loc[0]:lit_loc[-1]+1, 0]
128+
tvb_ilit = tvb_synaptic[:, 1, lit_loc[0]:lit_loc[-1]+1, 0]
129+
130+
# add all units WITHIN each region together across space to calculate
131+
# synaptic activity in EACH brain region
132+
v1_syn = np.sum(tvb_ev1+tvb_iv1, axis=1)
133+
v4_syn = np.sum(tvb_ev4+tvb_iv4, axis=1)
134+
it_syn = np.sum(tvb_eit+tvb_iit, axis=1)
135+
d1_syn = np.sum(tvb_ed1+tvb_id1, axis=1)
136+
d2_syn = np.sum(tvb_ed2+tvb_id2, axis=1)
137+
fs_syn = np.sum(tvb_efs+tvb_ifs, axis=1)
138+
fr_syn = np.sum(tvb_efr+tvb_ifr, axis=1)
139+
140+
# now, add unit across space in the contralateral IT
141+
lit_syn = np.sum(tvb_elit + tvb_ilit, axis=1)
142+
143+
# create a numpy array of timeseries
144+
synaptic = np.array([v1_syn, v4_syn, it_syn, fs_syn, d1_syn, d2_syn, fr_syn, lit_syn])
145+
146+
# now, save all synaptic timeseries to a single file
147+
np.save(syn_file, synaptic)
148+
149+
# Extract number of timesteps from one of the matrices
150+
timesteps = v1_syn.shape[0]
151+
print 'Timesteps = ', timesteps
152+
153+
# Construct a numpy array of timesteps (data points provided in data file)
154+
# to convert from timesteps to time in seconds we do the following:
155+
# Each simulation time-step equals 5 milliseconds
156+
# However, we are recording only once every 10 time-steps
157+
# Therefore, each data point in the output files represents 50 milliseconds.
158+
# Thus, we need to multiply the datapoint times 50 ms...
159+
# ... and divide by 1000 to convert to seconds
160+
#t = np.linspace(0, 659*50./1000., num=660)
161+
t = np.linspace(0, timesteps * 50.0 / 1000., num=timesteps)
162+
163+
164+
# Set up figures to plot synaptic activity
165+
plt.figure()
166+
plt.suptitle('SIMULATED SYNAPTIC ACTIVITY IN V1')
167+
plt.plot(t, v1_syn)
168+
# Set up figures to plot synaptic activity
169+
plt.figure()
170+
plt.suptitle('SIMULATED SYNAPTIC ACTIVITY IN V4')
171+
plt.plot(v4_syn)
172+
# Set up figures to plot synaptic activity
173+
plt.figure()
174+
plt.suptitle('SIMULATED SYNAPTIC ACTIVITY IN IT')
175+
plt.plot(it_syn)
176+
# Set up figures to plot synaptic activity
177+
plt.figure()
178+
plt.suptitle('SIMULATED SYNAPTIC ACTIVITY IN FS')
179+
plt.plot(fs_syn)
180+
# Set up figures to plot synaptic activity
181+
plt.figure()
182+
plt.suptitle('SIMULATED SYNAPTIC ACTIVITY IN D1')
183+
plt.plot(d1_syn)
184+
# Set up figures to plot synaptic activity
185+
plt.figure()
186+
plt.suptitle('SIMULATED SYNAPTIC ACTIVITY IN D2')
187+
plt.plot(d2_syn)
188+
# Set up figures to plot synaptic activity
189+
plt.figure()
190+
plt.suptitle('SIMULATED SYNAPTIC ACTIVITY IN FR')
191+
plt.plot(fr_syn)
192+
# Set up figures to plot synaptic activity
193+
plt.figure()
194+
plt.suptitle('SIMULATED SYNAPTIC ACTIVITY IN LEFT IT')
195+
plt.plot(lit_syn)
196+
197+
# Show the plots on the screen
198+
plt.show()

0 commit comments

Comments
 (0)