diff --git a/code/spikes_analysis.py b/code/spikes_analysis.py index df95c5e..ea11eb8 100644 --- a/code/spikes_analysis.py +++ b/code/spikes_analysis.py @@ -1,5 +1,6 @@ import matplotlib.pyplot as plt import numpy as np +import scipy.stats as ss from read_chirp_data import * from utility import * from IPython import embed @@ -9,14 +10,15 @@ sampling_rate = 40 #kHz data_dir = "../data" dataset = "2018-11-09-ad-invivo-1" # parameters for binning, smoothing and plotting +cut_window = 60 chirp_size = 14 #ms neuronal_delay = 5 #ms -chirp_start = int((-chirp_size/2+neuronal_delay+50)*sampling_rate) -chirp_end = int((chirp_size/2+neuronal_delay+51)*sampling_rate) +chirp_start = int((-chirp_size/2+neuronal_delay+cut_window)*sampling_rate) +chirp_end = int((chirp_size/2+neuronal_delay+cut_window+1)*sampling_rate) num_bin = 12 window = 1 #ms -time_axis = np.arange(-50, 50, 1/sampling_rate) #steps -spike_bins = np.arange(-50, 51) #ms +time_axis = np.arange(-cut_window, cut_window, 1/sampling_rate) #steps +spike_bins = np.arange(-cut_window, cut_window+1) #ms # read data from files spikes = read_chirp_spikes(os.path.join(data_dir, dataset)) @@ -24,17 +26,11 @@ eod = read_chirp_eod(os.path.join(data_dir, dataset)) chirp_times = read_chirp_times(os.path.join(data_dir, dataset)) # make a delta f map for the quite more complicated keys -df_map = {} -for k in spikes.keys(): - df = k[1] - if df in df_map.keys(): - df_map[df].append(k) - else: - df_map[df] = [k] +df_map = map_keys(spikes) # differentiate between phases phase_vec = np.arange(0, 1+1/num_bin, 1/num_bin) -cut_range = np.arange(-50*sampling_rate, 50*sampling_rate, 1) +cut_range = np.arange(-cut_window*sampling_rate, cut_window*sampling_rate, 1) # make dictionaries for spiketimes df_phase_time = {} @@ -52,7 +48,7 @@ for deltaf in df_map.keys(): # get spikes between 50 ms befor and after the chirp spikes_to_cut = np.asarray(spikes[rep][phase]) - spikes_cut = spikes_to_cut[(spikes_to_cut > -50) & (spikes_to_cut < 50)] + spikes_cut = spikes_to_cut[(spikes_to_cut > -cut_window) & (spikes_to_cut < cut_window)] spikes_idx = np.round(spikes_cut*sampling_rate) # also save as binary, 0 no spike, 1 spike binary_spikes = np.isin(cut_range, spikes_idx)*1 @@ -65,28 +61,47 @@ for deltaf in df_map.keys(): df_phase_time[deltaf][idx] = [spikes_cut] df_phase_binary[deltaf][idx] = binary_spikes + # for plotting and calculating iterate over delta f and phases for df in df_phase_time.keys(): + beat_duration = int(abs(1/df*1000)) #ms + beat_window = 0 + while beat_window+beat_duration <= cut_window: + beat_window = beat_window+beat_duration for phase in df_phase_time[df].keys(): trials_binary = df_phase_binary[df][phase] - sr_chirp = np.zeros(len(trials_binary)) - sr_beat = np.zeros(len(trials_binary)) + train_chirp = [] + train_beat = [] + spikerate_chirp = np.zeros(len(trials_binary)) + spikerate_beat = np.zeros(len(trials_binary)) for i, trial in enumerate(trials_binary): smoothed_trial = smooth(trial, window, 1/sampling_rate) - sr_chirp[i] = np.mean(smoothed_trial[chirp_start:chirp_end]) - sr_beat[i] = np.mean(smoothed_trial[0:chirp_start]) - - for rate_chirp in sr_chirp: - for rate_beat in sr_beat: - r = np.corrcoef(rate_chirp, rate_beat) - print(r) - - embed() - exit() - - #csi = (spikerate_chirp-spikerate_befor)/(spikerate_chirp+spikerate_befor) + train_chirp.append(smoothed_trial[chirp_start:chirp_end]) + train_beat.append(smoothed_trial[chirp_start-beat_window:chirp_start]) + spikerate_chirp[i] = np.mean(smoothed_trial[chirp_start:chirp_end]) + spikerate_beat[i] = np.mean(smoothed_trial[chirp_start-beat_window:chirp_start]) + + rcs = [] + rbs = [] + for i, train in enumerate(train_chirp): + for j, train2 in enumerate(train_chirp): + if np.array_equal(train, train2): + continue + else: + rc, _ = ss.pearsonr(train, train2) + rb, _ = ss.pearsonr(train_beat[i], train_beat[j]) + rcs.append(rc) + rbs.append(rb) + + r_train_chirp = np.mean(rcs) + r_train_beat = np.mean(rbs) + + csi_train = (r_train_chirp - r_train_beat) / (r_train_chirp + r_train_beat) + print(csi_train) + csi_rate = (np.std(spikerate_chirp) - np.std(spikerate_beat)) / (np.std(spikerate_chirp) + np.std(spikerate_beat)) + print(csi_rate) # plot #plot_trials = df_phase_time[df][phase]