Added code to process auditory output data and to visualize connectivity matrix of LSNM models

Author: Antonio Ulloa

analysis/compute_BOLD_balloon_model_auditory.py

@@ -75,9 +75,9 @@
 
 # define the name of the input files where the synaptic activities are stored
 # we have four files: TC-PSL, Tones-PSL, TC-DMS, Tones-DMS
-SYN_file_2 = 'output.TC_PSL/synaptic_in_ROI.npy'
+SYN_file_4 = 'output.TC_PSL/synaptic_in_ROI.npy'
 SYN_file_3 = 'output.Tones_PSL/synaptic_in_ROI.npy'
-SYN_file_4 = 'output.TC_DMS/synaptic_in_ROI.npy'
+SYN_file_2 = 'output.TC_DMS/synaptic_in_ROI.npy'
 SYN_file_1 = 'output.Tones_DMS/synaptic_in_ROI.npy'
 
 # define the name of the output file where the BOLD timeseries will be stored

analysis/compute_NAI_auditory.py

@@ -53,9 +53,23 @@
 # total number of trials in simulated experiment
 number_of_trials = 12
 
+# define the names of the input files where the MEG timeseries are contained, for four conditions:
+# TC-PSL, Tones-PSL, TC-DMS and Tones-DMS
+#MEG_TC_PSL_subj = np.array(['subject_4_with_feedback/output.TC_PSL/meg_source_activity.npy',
+#                             'subject_original_with_feedback/output.TC_PSL/meg_source_activity.npy'])
+#MEG_Tones_PSL_subj = np.array(['subject_4_with_feedback/output.Tones_PSL/meg_source_activity.npy',
+#                             'subject_original_with_feedback/output.Tones_PSL/meg_source_activity.npy'])
+#MEG_TC_DMS_subj = np.array(['subject_4_with_feedback/output.TC_DMS/meg_source_activity.npy',
+#                             'subject_original_with_feedback/output.TC_DMS/meg_source_activity.npy'])
+#MEG_Tones_DMS_subj = np.array(['subject_4_with_feedback/output.Tones_DMS/meg_source_activity.npy',
+#                                'subject_original_with_feedback/output.Tones_DMS/meg_source_activity.npy'])
+
 # name of the input file where MEG source activity is stored
 MEG_source_file = 'meg_source_activity.npy'
 
+# name of the ouput file where the NAI timecourses will be stored
+NAI_file = 'NAI_timecourse.npy'
+
 # Load MEG source activity data files
 MEG_source_activity = np.load(MEG_source_file)
 
@@ -120,10 +134,15 @@
 
 print 'Modulation Index:', MI
 
+NAI_timecourse = aud_MEG_source_activity
+
+# now, save all NAI timecourses to a single file 
+np.save(NAI_file, NAI_timecourse)
+
 # Set up figure to plot MEG source dynamics averaged across trials
 fig = plt.figure(1)
 
-plt.suptitle('MEG SOURCE DYNAMICS AVERAGED ACROSS TRIALS')
+plt.suptitle('NEURAL ACTIVITY INDEX TIME COURSE AVERAGED ACROSS TRIALS')
 
 # Plot MEG signal
 aud_plot_S1 = plt.plot(aud_MEG_source_activity[0], label='S1')

analysis/compute_PSC_auditory.py

@@ -83,7 +83,7 @@
 num_of_syn_blocks = 1              # we have more synaptic blocks than fmri blocks
                                     # because we get rid of blocks in BOLD timeseries
 
-num_of_subjects = 1
+num_of_subjects = 4
 
 scans_removed = 0             # number of scans removed from BOLD computation
 synaptic_steps_removed = 0    # number of synaptic steps removed from synaptic
@@ -102,10 +102,25 @@
 
 # define the names of the input files where the BOLD timeseries are contained, for four conditions:
 # TC-PSL, Tones-PSL, TC-DMS and Tones-DMS
-BOLD_TC_PSL_subj = np.array(['subject_original_with_feedback/output.TC_PSL/lsnm_bold_balloon.npy'])
-BOLD_Tones_PSL_subj = np.array(['subject_original_with_feedback/output.Tones_PSL/lsnm_bold_balloon.npy'])
-BOLD_TC_DMS_subj = np.array(['subject_original_with_feedback/output.TC_DMS/lsnm_bold_balloon.npy'])
-BOLD_Tones_DMS_subj = np.array(['subject_original_with_feedback/output.Tones_DMS/lsnm_bold_balloon.npy'])
+BOLD_TC_PSL_subj = np.array(['subject_original_with_feedback/output.TC_PSL/lsnm_bold_balloon.npy',
+                             'subject_4_with_feedback/output.TC_PSL/lsnm_bold_balloon.npy',
+                             'subject_5_with_feedback/output.TC_PSL/lsnm_bold_balloon.npy',
+                             'subject_6_with_feedback/output.TC_PSL/lsnm_bold_balloon.npy'])
+
+BOLD_Tones_PSL_subj = np.array(['subject_original_with_feedback/output.Tones_PSL/lsnm_bold_balloon.npy',
+                             'subject_4_with_feedback/output.Tones_PSL/lsnm_bold_balloon.npy',
+                             'subject_5_with_feedback/output.Tones_PSL/lsnm_bold_balloon.npy',
+                             'subject_6_with_feedback/output.Tones_PSL/lsnm_bold_balloon.npy'])
+
+BOLD_TC_DMS_subj = np.array(['subject_original_with_feedback/output.TC_DMS/lsnm_bold_balloon.npy',
+                             'subject_4_with_feedback/output.TC_DMS/lsnm_bold_balloon.npy',
+                             'subject_5_with_feedback/output.TC_DMS/lsnm_bold_balloon.npy',
+                             'subject_6_with_feedback/output.TC_DMS/lsnm_bold_balloon.npy'])
+
+BOLD_Tones_DMS_subj = np.array(['subject_original_with_feedback/output.Tones_DMS/lsnm_bold_balloon.npy',
+                                'subject_4_with_feedback/output.Tones_DMS/lsnm_bold_balloon.npy',
+                                'subject_5_with_feedback/output.Tones_DMS/lsnm_bold_balloon.npy',
+                                'subject_6_with_feedback/output.Tones_DMS/lsnm_bold_balloon.npy'])
 
 # set matplot lib parameters to produce visually appealing plots
 mpl.style.use('ggplot')
@@ -155,10 +170,11 @@
 print BOLD_TC_PSL.shape
 
 # Perform time-course normalization by converting each timeseries to percent signal change
-# for each subject and for each module (DMS relative to the mean of the PSL condition).
+# for each subject and for each module (TC and Tones DMS relative to the mean of the Tones PSL
+# condition, as per Husain et al (2004), pp 1711, "Simulating PET and fMRI".
 for s in range(0, num_of_subjects):
     for m in range(0, num_of_modules):
-        timecourse_mean = np.mean(BOLD_TC_PSL[s,m])
+        timecourse_mean = np.mean(BOLD_Tones_DMS[s,m])
         BOLD_TC_DMS[s,m] = BOLD_TC_DMS[s,m] / timecourse_mean * 100. - 100.
 
 for s in range(0, num_of_subjects):
@@ -167,14 +183,14 @@
         BOLD_Tones_DMS[s,m] = BOLD_Tones_DMS[s,m] / timecourse_mean * 100. - 100.
 
 # now calculate PSC of TC-DMS with respect ot Tones-DMS
-BOLD_DMS = BOLD_TC_DMS - BOLD_Tones_DMS
+BOLD_DMS = BOLD_TC_DMS #- BOLD_Tones_DMS
 
 mean_PSC_dms = np.zeros((num_of_subjects, num_of_modules))
 
-# now, perform a mean of PSCs of chosen scans
+# now, perform a mean of PSCs
 for s in range(0, num_of_subjects):
     for m in range(0, num_of_modules):
-        mean_PSC_dms[s,m] = np.mean(BOLD_DMS[s,m,8:17])
+        mean_PSC_dms[s,m] = np.mean(BOLD_DMS[s,m][5:10])
 
 # now, calculate the average PSC across subjects for each brain region and for each condition:
 BOLD_a1_dms_avg = np.mean(mean_PSC_dms[:,0]) 
@@ -303,7 +319,7 @@
 
 fig, ax = plt.subplots()
 
-#ax.set_ylim([0,3.5])
+#ax.set_ylim([0,16])
 
 # now, group the values to be plotted by brain module and by task condition
 

analysis/compute_avg_NAI_auditory.py

@@ -0,0 +1,132 @@
+# ============================================================================
+#
+#                            PUBLIC DOMAIN NOTICE
+#
+#       National Institute on Deafness and Other Communication Disorders
+#
+# This software/database is a "United States Government Work" under the 
+# terms of the United States Copyright Act. It was written as part of 
+# the author's official duties as a United States Government employee and 
+# thus cannot be copyrighted. This software/database is freely available 
+# to the public for use. The NIDCD and the U.S. Government have not placed 
+# any restriction on its use or reproduction. 
+#
+# Although all reasonable efforts have been taken to ensure the accuracy 
+# and reliability of the software and data, the NIDCD and the U.S. Government 
+# do not and cannot warrant the performance or results that may be obtained 
+# by using this software or data. The NIDCD and the U.S. Government disclaim 
+# all warranties, express or implied, including warranties of performance, 
+# merchantability or fitness for any particular purpose.
+#
+# Please cite the author in any work or product based on this material.
+# 
+# ==========================================================================
+
+
+
+# ***************************************************************************
+#
+#   Large-Scale Neural Modeling software (LSNM)
+#
+#   Section on Brain Imaging and Modeling
+#   Voice, Speech and Language Branch
+#   National Institute on Deafness and Other Communication Disorders
+#   National Institutes of Health
+#
+#   This file (compute_avg_NAI_auditory.py) was created on April 19, 2016.
+#
+#
+#   Author: Antonio Ulloa. Last updated by Antonio Ulloa on April 19, 2016  
+# **************************************************************************/
+
+# compute_avg_NAI_auditory.py
+#
+# Compute average of Neural Activity Index (NAI) of auditory primary cortex using NAI's for a 
+# a number of subjects (previously stored in a python data file)
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+num_of_subjects = 2
+NAI_ts_length = 30
+
+# define the names of the input files where the NAI timecourses are contained, for four conditions:
+# TC-PSL, Tones-PSL, TC-DMS and Tones-DMS
+NAI_TC_PSL_subj = np.array(['subject_4_with_feedback/output.TC_PSL/NAI_timecourse.npy',
+                             'subject_original_with_feedback/output.TC_PSL/NAI_timecourse.npy'])
+NAI_Tones_PSL_subj = np.array(['subject_4_with_feedback/output.Tones_PSL/NAI_timecourse.npy',
+                             'subject_original_with_feedback/output.Tones_PSL/NAI_timecourse.npy'])
+NAI_TC_DMS_subj = np.array(['subject_4_with_feedback/output.TC_DMS/NAI_timecourse.npy',
+                             'subject_original_with_feedback/output.TC_DMS/NAI_timecourse.npy'])
+NAI_Tones_DMS_subj = np.array(['subject_4_with_feedback/output.Tones_DMS/NAI_timecourse.npy',
+                                'subject_original_with_feedback/output.Tones_DMS/NAI_timecourse.npy'])
+
+
+
+# Load NAI source activity data files
+# open files containing NAI time-series
+NAI_TC_PSL = np.zeros((num_of_subjects, 2, NAI_ts_length))
+NAI_Tones_PSL = np.zeros((num_of_subjects, 2, NAI_ts_length))
+NAI_TC_DMS = np.zeros((num_of_subjects, 2, NAI_ts_length))
+NAI_Tones_DMS = np.zeros((num_of_subjects, 2, NAI_ts_length))
+for idx in range(0, num_of_subjects):
+    NAI_TC_PSL[idx] = np.load(NAI_TC_PSL_subj[idx])    
+    NAI_Tones_PSL[idx] = np.load(NAI_Tones_PSL_subj[idx])    
+    NAI_TC_DMS[idx] = np.load(NAI_TC_DMS_subj[idx])    
+    NAI_Tones_DMS[idx] = np.load(NAI_Tones_DMS_subj[idx])    
+
+# average NAI across subjects
+for idx in range(0, num_of_subjects):
+    NAI_TC_PSL_avg = np.mean(NAI_TC_PSL, axis=0)
+    NAI_Tones_PSL_avg = np.mean(NAI_Tones_PSL, axis=0)
+    NAI_TC_DMS_avg = np.mean(NAI_TC_DMS, axis=0)
+    NAI_Tones_DMS_avg = np.mean(NAI_Tones_DMS, axis=0)
+
+
+
+# plot the average NAI's
+
+# Set up figure to plot MEG source dynamics averaged across trials
+fig = plt.figure(1)
+
+plt.suptitle('NAI average timecourse during TC DMS')
+
+# Plot NAI
+aud_plot_S1 = plt.plot(NAI_TC_DMS_avg[0], label='S1')
+aud_plot_S2 = plt.plot(NAI_TC_DMS_avg[1], label='S2')
+
+plt.legend()
+
+# Set up figure to plot MEG source dynamics averaged across trials
+fig = plt.figure(2)
+
+plt.suptitle('NAI average timecourse during Tones DMS')
+
+# Plot NAI
+aud_plot_S1 = plt.plot(NAI_Tones_DMS_avg[0], label='S1')
+aud_plot_S2 = plt.plot(NAI_Tones_DMS_avg[1], label='S2')
+
+plt.legend()
+# Set up figure to plot MEG source dynamics averaged across trials
+fig = plt.figure(3)
+
+plt.suptitle('NAI average timecourse during TC PSL')
+
+# Plot NAI
+aud_plot_S1 = plt.plot(NAI_TC_PSL_avg[0], label='S1')
+aud_plot_S2 = plt.plot(NAI_TC_PSL_avg[1], label='S2')
+
+plt.legend()
+# Set up figure to plot MEG source dynamics averaged across trials
+fig = plt.figure(4)
+
+plt.suptitle('NAI average timecourse during Tones PSL')
+
+# Plot NAI
+aud_plot_S1 = plt.plot(NAI_Tones_PSL_avg[0], label='S1')
+aud_plot_S2 = plt.plot(NAI_Tones_PSL_avg[1], label='S2')
+
+plt.legend()
+
+# Show the plot on the screen
+plt.show()

analysis/compute_syn_auditory.py

@@ -116,15 +116,15 @@
 
 # add all units WITHIN each region together across space to calculate
 # synaptic activity in EACH brain region
-a1_syn = np.sum(ea1u + ea1d + ia1u + ia1d, axis = 1) + np.sum(tvb_ea1+tvb_ia1, axis=1)
-a2_syn = np.sum(ea2u + ea2c + ea2d + ia2u + ia2c + ia2d, axis = 1) + np.sum(tvb_ea2+tvb_ia2, axis=1)
-st_syn = np.sum(estg + istg, axis = 1) + np.sum(tvb_est+tvb_ist, axis=1)
+a1_syn = np.sum(ea1u + ea1d + ia1u + ia1d, axis = 1) #+ np.sum(tvb_ea1+tvb_ia1, axis=1)
+a2_syn = np.sum(ea2u + ea2c + ea2d + ia2u + ia2c + ia2d, axis = 1) #+ np.sum(tvb_ea2+tvb_ia2, axis=1)
+st_syn = np.sum(estg + istg, axis = 1) #+ np.sum(tvb_est+tvb_ist, axis=1)
 d1_syn = np.sum(efd1 + ifd1, axis = 1)
 d2_syn = np.sum(efd2 + ifd2, axis = 1)
 fs_syn = np.sum(exfs + infs, axis = 1)
 fr_syn = np.sum(exfr + infr, axis = 1)
 
-pf_syn = d1_syn + d2_syn + fs_syn + fr_syn + np.sum(tvb_epf+tvb_ipf, axis=1)
+pf_syn = d1_syn + d2_syn + fs_syn + fr_syn #+ np.sum(tvb_epf+tvb_ipf, axis=1)
 
 # get rid of the first time point ('zero point') bc it could skew correlations later
 #v1_syn[0] = v1_syn[1]
@@ -136,16 +136,6 @@
 #fr_syn[0] = fr_syn[1]
 #lit_syn[0] = lit_syn[1]
 
-# ...and normalize the synaptic activities of each module (convert to percentage signal change)
-#v1_syn = v1_syn / np.mean(v1_syn) * 100. - 100.
-#v4_syn = v4_syn / np.mean(v4_syn) * 100. - 100.
-#it_syn = it_syn / np.mean(it_syn) * 100. - 100.
-#d1_syn = d1_syn / np.mean(d1_syn) * 100. - 100.
-#d2_syn = d2_syn / np.mean(d2_syn) * 100. - 100.
-#fs_syn = fs_syn / np.mean(fs_syn) * 100. - 100.
-#fr_syn = fr_syn / np.mean(fr_syn) * 100. - 100.
-#lit_syn= lit_syn/ np.mean(lit_syn)* 100. - 100.
-
 # create a numpy array of timeseries
 synaptic = np.array([a1_syn, a2_syn, st_syn, pf_syn])
 

auditory_model/subject_2_with_feedback/output.TC_PSL/log.txt

@@ -0,0 +1,2 @@
+Simulation Start Time: Sun Apr 17 13:21:45 2016
+Simulation End Time: Sun Apr 17 15:43:14 2016

auditory_model/subject_2_with_feedback/output.Tones_DMS/log.txt

@@ -1,2 +1,2 @@
-Simulation Start Time: Sun Apr 17 10:57:11 2016
-Simulation End Time: Sun Apr 17 11:02:08 2016
+Simulation Start Time: Sun Apr 17 11:13:38 2016
+Simulation End Time: Sun Apr 17 13:07:49 2016

auditory_model/subject_2_with_feedback/output.Tones_PSL/log.txt

@@ -0,0 +1,2 @@
+Simulation Start Time: Sun Apr 17 16:03:14 2016
+Simulation End Time: Sun Apr 17 18:07:24 2016

auditory_model/subject_3_with_feedback/neuralnet.json

@@ -0,0 +1,2 @@
+Simulation Start Time: Sun Apr 17 18:12:33 2016
+Simulation End Time: Sun Apr 17 20:20:30 2016