# ============================================================================ # # 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_PSC.py) was created on December 2, 2015. # # Author: Antonio Ulloa # # Last updated by Antonio Ulloa on February 4, 2016 # **************************************************************************/ # compute_PSC.py # # Calculate and plot the Percent Signal Change (PSC) of fMRI BOLD timeseries across # subjects for 2 conditions: DMS and passive viewing. # # The inputs are BOLD timeseries for each subject, one timeseries per brain module. # We take each timeseries and convert # each to PSC by dividing every time point by the baseline (the mean of the whole # time course for a given subject and brain area), then multiplying by 100 # (Source: BrainVoyager v20.0 User's Guide). That way, # we obtain a # whole-experiment timecourse in percent signal change, per brain area, per subject. # Then, we average together all time points for each condition, DMS and control, # separately, across subjects, per brain area. That gives us one PSC number per # brain area per condition, for all subjects. # # We also do a paired t-test to check for the statistical significance of the difference # between DMS and control conditions, per brain area. # import numpy as np import matplotlib.pyplot as plt import matplotlib as mpl import pandas as pd import math as math from scipy.stats import poisson from scipy.stats import t # define the length of both each trial and the whole experiment # in synaptic timesteps, as well as total number of trials experiment_length = 3960 trial_length = 110 num_of_trials = 36 # number of trials per subject num_of_fmri_blocks = 12 # how many blocks of trials in the experiment num_of_syn_blocks = 12 # we have more synaptic blocks than fmri blocks # because we get rid of blocks in BOLD timeseries num_of_subjects = 10 scans_removed = 0 # number of scans removed from BOLD computation synaptic_steps_removed = 0 # number of synaptic steps removed from synaptic # computation num_of_modules = 8 # V1, V4, IT, FS, D1, D2, FR, cIT (contralateral IT) # define total number of trials across subjects total_trials = num_of_trials * num_of_subjects # define total number of blocks across subjects total_fmri_blocks = num_of_fmri_blocks * num_of_subjects #total_syn_blocks = num_of_syn_blocks * num_of_subjects synaptic_timesteps = experiment_length - synaptic_steps_removed # define the names of the input files where the BOLD timeseries are contained BOLD_ts_subj=np.array(['subject_11/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy', 'subject_12/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy', 'subject_13/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy', 'subject_14/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy', 'subject_15/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy', 'subject_16/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy', 'subject_17/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy', 'subject_18/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy', 'subject_19/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy', 'subject_20/output.36trials.with_feedback/bold_balloon_TVB_ROI.npy']) # set matplot lib parameters to produce visually appealing plots #mpl.style.use('ggplot') # define output file where means, standard deviations, and variances will be stored PSC_stats_FILE = 'psc_stats_TVB_ROI.txt' # define neural synaptic time interval in seconds. The simulation data is collected # one data point at synaptic intervals (10 simulation timesteps). Every simulation # timestep is equivalent to 5 ms. Ti = 0.005 * 10 # Total time of scanning experiment in seconds (timesteps X 5) T = 198 # Time for one complete trial in milliseconds Ttrial = 5.5 # the scanning happened every Tr interval below (in milliseconds). This # is the time needed to sample hemodynamic activity to produce # each fMRI image. Tr = 2 num_of_scans = T / Tr - scans_removed # Construct a numpy array of time points time_in_seconds = np.arange(0, T, Tr) scanning_timescale = np.arange(0, synaptic_timesteps, synaptic_timesteps / (T/Tr)) # open files containing BOLD time-series BOLD_ts = np.zeros((num_of_subjects, num_of_modules, num_of_scans)) for idx in range(0, num_of_subjects): BOLD_ts[idx] = np.load(BOLD_ts_subj[idx]) # Perform time-course normalization by converting each timeseries to percent signal change # for each subject and for each module relative to the mean of each time course. for s in range(0, num_of_subjects): for m in range(0, num_of_modules): timecourse_mean = np.mean(BOLD_ts[s,m]) BOLD_ts[s,m] = BOLD_ts[s,m] / timecourse_mean * 100. - 100. mean_PSC_dms = np.zeros((num_of_subjects, num_of_modules)) mean_PSC_ctl = np.zeros((num_of_subjects, num_of_modules)) scans_in_each_block = int(round(num_of_scans / (num_of_fmri_blocks / 2.), 0)) print 'Scans in each block: ', scans_in_each_block # now, perform a mean of PSCs of scan 4, 5, 6 of each condition for s in range(0, num_of_subjects): for m in range(0, num_of_modules): for i in range(0, num_of_scans, scans_in_each_block): mean_PSC_dms[s,m] = (BOLD_ts[s,m,i+4] + BOLD_ts[s,m,i+5] + BOLD_ts[s,m,i+6]) / 3. mean_PSC_ctl[s,m] = (BOLD_ts[s,m,i+11]+ BOLD_ts[s,m,i+12]+ BOLD_ts[s,m,i+13]) / 3. # now, calculate the average PSC across subjects for each brain region and for each condition: BOLD_v1_dms_avg = np.mean(mean_PSC_dms[:,0]) BOLD_v4_dms_avg = np.mean(mean_PSC_dms[:,1]) BOLD_it_dms_avg = np.mean(mean_PSC_dms[:,2]) BOLD_fs_dms_avg = np.mean(mean_PSC_dms[:,3]) BOLD_d1_dms_avg = np.mean(mean_PSC_dms[:,4]) BOLD_d2_dms_avg = np.mean(mean_PSC_dms[:,5]) BOLD_fr_dms_avg = np.mean(mean_PSC_dms[:,6]) BOLD_v1_ctl_avg = np.mean(mean_PSC_ctl[:,0]) BOLD_v4_ctl_avg = np.mean(mean_PSC_ctl[:,1]) BOLD_it_ctl_avg = np.mean(mean_PSC_ctl[:,2]) BOLD_fs_ctl_avg = np.mean(mean_PSC_ctl[:,3]) BOLD_d1_ctl_avg = np.mean(mean_PSC_ctl[:,4]) BOLD_d2_ctl_avg = np.mean(mean_PSC_ctl[:,5]) BOLD_fr_ctl_avg = np.mean(mean_PSC_ctl[:,6]) # calculate the variance as well, across subjects, for each brain regions and for each condition: BOLD_v1_dms_var = np.var(mean_PSC_dms[:,0]) BOLD_v4_dms_var = np.var(mean_PSC_dms[:,1]) BOLD_it_dms_var = np.var(mean_PSC_dms[:,2]) BOLD_fs_dms_var = np.var(mean_PSC_dms[:,3]) BOLD_d1_dms_var = np.var(mean_PSC_dms[:,4]) BOLD_d2_dms_var = np.var(mean_PSC_dms[:,5]) BOLD_fr_dms_var = np.var(mean_PSC_dms[:,6]) BOLD_v1_ctl_var = np.var(mean_PSC_ctl[:,0]) BOLD_v4_ctl_var = np.var(mean_PSC_ctl[:,1]) BOLD_it_ctl_var = np.var(mean_PSC_ctl[:,2]) BOLD_fs_ctl_var = np.var(mean_PSC_ctl[:,3]) BOLD_d1_ctl_var = np.var(mean_PSC_ctl[:,4]) BOLD_d2_ctl_var = np.var(mean_PSC_ctl[:,5]) BOLD_fr_ctl_var = np.var(mean_PSC_ctl[:,6]) # calculate the standard deviation of the PSC across subjects, for each brain region and for each # condition BOLD_v1_dms_std = np.std(mean_PSC_dms[:,0]) BOLD_v4_dms_std = np.std(mean_PSC_dms[:,1]) BOLD_it_dms_std = np.std(mean_PSC_dms[:,2]) BOLD_fs_dms_std = np.std(mean_PSC_dms[:,3]) BOLD_d1_dms_std = np.std(mean_PSC_dms[:,4]) BOLD_d2_dms_std = np.std(mean_PSC_dms[:,5]) BOLD_fr_dms_std = np.std(mean_PSC_dms[:,6]) BOLD_v1_ctl_std = np.std(mean_PSC_ctl[:,0]) BOLD_v4_ctl_std = np.std(mean_PSC_ctl[:,1]) BOLD_it_ctl_std = np.std(mean_PSC_ctl[:,2]) BOLD_fs_ctl_std = np.std(mean_PSC_ctl[:,3]) BOLD_d1_ctl_std = np.std(mean_PSC_ctl[:,4]) BOLD_d2_ctl_std = np.std(mean_PSC_ctl[:,5]) BOLD_fr_ctl_std = np.std(mean_PSC_ctl[:,6]) # Display all PSCs for both DMS and CTL, which represent the average of # the Percent Signal Change per each brain area per each module print 'BOLD V1 DMS mean, std, var: [', BOLD_v1_dms_avg, BOLD_v1_dms_std, BOLD_v1_dms_var, ']' print 'BOLD V4 DMS mean, std, var: [', BOLD_v4_dms_avg, BOLD_v4_dms_std, BOLD_v4_dms_var, ']' print 'BOLD IT DMS mean, std, var: [', BOLD_it_dms_avg, BOLD_it_dms_std, BOLD_it_dms_var, ']' print 'BOLD FS DMS mean, std, var: [', BOLD_fs_dms_avg, BOLD_fs_dms_std, BOLD_fs_dms_var, ']' print 'BOLD D1 DMS mean, std, var: [', BOLD_d1_dms_avg, BOLD_d1_dms_std, BOLD_d1_dms_var, ']' print 'BOLD D2 DMS mean, std, var: [', BOLD_d2_dms_avg, BOLD_d2_dms_std, BOLD_d2_dms_var, ']' print 'BOLD FR DMS mean, std, var: [', BOLD_fr_dms_avg, BOLD_fr_dms_std, BOLD_fr_dms_var, ']' print 'BOLD V1 CTL mean, std, var: [', BOLD_v1_ctl_avg, BOLD_v1_ctl_std, BOLD_v1_ctl_var, ']' print 'BOLD V4 CTL mean, std, var: [', BOLD_v4_ctl_avg, BOLD_v4_ctl_std, BOLD_v4_ctl_var, ']' print 'BOLD IT CTL mean, std, var: [', BOLD_it_ctl_avg, BOLD_it_ctl_std, BOLD_it_ctl_var, ']' print 'BOLD FS CTL mean, std, var: [', BOLD_fs_ctl_avg, BOLD_fs_ctl_std, BOLD_fs_ctl_var, ']' print 'BOLD D1 CTL mean, std, var: [', BOLD_d1_ctl_avg, BOLD_d1_ctl_std, BOLD_d1_ctl_var, ']' print 'BOLD D2 CTL mean, std, var: [', BOLD_d2_ctl_avg, BOLD_d2_ctl_std, BOLD_d2_ctl_var, ']' print 'BOLD FR CTL mean, std, var: [', BOLD_fr_ctl_avg, BOLD_fr_ctl_std, BOLD_fr_ctl_var, ']' # Calculate the statistical significance by using a one-tailed paired t-test: # We are going to have one group of 10 subjects, doing both DMS and control task # STEPS: # (1) Set up hypotheses: # The NULL hypothesis is: # * The mean difference between paired observations (DMS and CTL) is zero # Our alternative hypothesis is: # * The mean difference between paired observations (DMS and CTL) is not zero # (2) Set a significance level: alpha = 0.05 # (3) What is the critical value and the rejection region? n = 10 - 1 # sample size minus 1 rejection_region = 1.833 # as found on t-test table for t and dof given, # values of t above rejection_region will be rejected # (4) compute the value of the test statistic # calculate differences between the pairs of data: d = mean_PSC_dms - mean_PSC_ctl # calculate the mean of those differences d_mean = np.mean(d, axis=0) # calculate the standard deviation of those differences d_std = np.std(d, axis=0) # calculate square root of sample size minus 1 sqrt_n = math.sqrt(n) # calculate standard error of the mean differences d_sem = d_std/sqrt_n # calculate the t statistic: t_star = d_mean / d_sem print 'Mean differences: ', d_mean print 't-values for PSC comparisons (V1, V4, IT, FS, D1, D2, FR, cIT): ', t_star print 'Dimensions of mean differences array', d_mean.shape print 'Dimensions of std of differences array', d_std.shape # increase font size prior to plotting plt.rcParams.update({'font.size': 15}) # define number of groups to plot N = 1 # create a list of x locations for each group index = np.arange(N) width = 0.2 # width of the bars fig, ax = plt.subplots() #ax.set_ylim([0,3.5]) # now, group the values to be plotted by brain module and by task condition rects_v1_dms = ax.bar(index, d_mean[0], width, color='yellow', label='V1', yerr=d_sem[0] ) rects_v4_dms = ax.bar(index + width, d_mean[1], width, color='green', label='V4', yerr=d_sem[1] ) rects_it_dms = ax.bar(index + width*2, d_mean[2], width, color='blue', label='IT', yerr=d_sem[2] ) rects_fs_dms = ax.bar(index + width*3, d_mean[3], width, color='orange', label='FS', yerr=d_sem[3] ) rects_d1_dms = ax.bar(index + width*4, d_mean[4], width, color='red', label='D1', yerr=d_sem[4] ) rects_d2_dms = ax.bar(index + width*5, d_mean[5], width, color='pink', label='D2', yerr=d_sem[5] ) rects_fr_dms = ax.bar(index + width*6, d_mean[6], width, color='purple', label='FR', yerr=d_sem[6] ) #ax.set_title('PSC ACROSS SUBJECTS IN ALL BRAIN REGIONS') # get rid of x axis ticks and labels ax.set_xticks([]) ax.set_ylabel('Signal change differences') # Shrink current axis by 10% to make space for legend box = ax.get_position() ax.set_position([box.x0, box.y0, box.width * 0.85, box.height]) # place a legend to the right of the figure plt.legend(loc='center left', bbox_to_anchor=(1.02, .5)) #set up figure to plot BOLD signal (normalized to PSC) fig=plt.figure(2) # plot V1 BOLD time-series in yellow ax = plt.subplot(7,1,1) ax.set_yticks([]) ax.set_xticks([]) ax.plot(BOLD_ts[0], linewidth=3.0, color='yellow') # plot V4 BOLD time-series in yellow ax = plt.subplot(7,1,2) ax.set_yticks([]) ax.set_xticks([]) ax.plot(BOLD_ts[1], linewidth=3.0, color='lime') # plot IT BOLD time-series in yellow ax = plt.subplot(7,1,3) ax.set_yticks([]) ax.set_xticks([]) ax.plot(BOLD_ts[2], linewidth=3.0, color='blue') # plot FS BOLD time-series in yellow ax = plt.subplot(7,1,4) ax.set_yticks([]) ax.set_xticks([]) ax.plot(BOLD_ts[3], linewidth=3.0, color='orange') # plot D1 BOLD time-series in yellow ax = plt.subplot(7,1,5) ax.set_yticks([]) ax.set_xticks([]) ax.plot(BOLD_ts[4], linewidth=3.0, color='red') # plot D2 BOLD time-series in yellow ax = plt.subplot(7,1,6) ax.set_yticks([]) ax.set_xticks([]) ax.plot(BOLD_ts[5], linewidth=3.0, color='pink') # plot FR BOLD time-series in yellow ax = plt.subplot(7,1,7) ax.set_yticks([]) ax.set_xticks([]) ax.plot(BOLD_ts[6], linewidth=3.0, color='darkorchid') # Show the plots on the screen plt.show()