From a9d3dc64d496c2f8ea3e97e3fc468b8885804978 Mon Sep 17 00:00:00 2001 From: dhruvm9 Date: Fri, 15 Sep 2023 16:39:40 -0400 Subject: [PATCH] Figure 5 --- Analysis/Figure5B.py | 83 +++++++++ Analysis/Figure5C.py | 116 ++++++++++++ Analysis/Figure5D.py | 172 ++++++++++++++++++ Analysis/Figure5E.py | 121 ++++++++++++ Analysis/Figure5F.py | 88 +++++++++ .../{detectDOWNs.py => detectDOWNs_HDC.py} | 0 Analysis/detectDOWNs_RSC.py | 74 ++++++++ 7 files changed, 654 insertions(+) create mode 100644 Analysis/Figure5B.py create mode 100644 Analysis/Figure5C.py create mode 100644 Analysis/Figure5D.py create mode 100644 Analysis/Figure5E.py create mode 100644 Analysis/Figure5F.py rename Analysis/{detectDOWNs.py => detectDOWNs_HDC.py} (100%) create mode 100644 Analysis/detectDOWNs_RSC.py diff --git a/Analysis/Figure5B.py b/Analysis/Figure5B.py new file mode 100644 index 0000000..f53816f --- /dev/null +++ b/Analysis/Figure5B.py @@ -0,0 +1,83 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Sep 15 16:05:14 2023 + +@author: dhruv +""" + +#import dependencies +import numpy as np +import pandas as pd +import scipy.io +import pynapple as nap +import os +import matplotlib.pyplot as plt +from scipy.stats import kendalltau + +#%% Data organization + +data_directory = '/media/DataDhruv/Dropbox (Peyrache Lab)/Peyrache Lab Team Folder/Data/AdrianPoSub/###AllPoSub' +datasets = np.genfromtxt(os.path.join(data_directory,'dataset_Hor_DM.list'), delimiter = '\n', dtype = str, comments = '#') + +readpath = '/media/DataDhruv/Dropbox (Peyrache Lab)/Peyrache Lab Team Folder/Projects/PoSub-UPstate/Data' +writepath = '/home/dhruv/Code/MehrotraLevenstein_2023/Analysis/OutputFiles' + +s = 'A3707-200317' +name = s.split('/')[-1] +path = os.path.join(data_directory, s) +rawpath = os.path.join(readpath,s) + +#%% Loading the data + +data = nap.load_session(rawpath, 'neurosuite') +data.load_neurosuite_xml(rawpath) +spikes = data.spikes + +#%% Load cell types + +filepath = os.path.join(path, 'Analysis') +listdir = os.listdir(filepath) +file = [f for f in listdir if 'CellTypes' in f] +celltype = scipy.io.loadmat(os.path.join(filepath,file[0])) + +pyr = [] +interneuron = [] +hd = [] + +for i in range(len(spikes)): + if celltype['ex'][i] == 1 and celltype['gd'][i] == 1: + pyr.append(i) + +for i in range(len(spikes)): + if celltype['fs'][i] == 1 and celltype['gd'][i] == 1: + interneuron.append(i) + +for i in range(len(spikes)): + if celltype['hd'][i] == 1 and celltype['gd'][i] == 1: + hd.append(i) + +#%% Load UP and DOWN states + +file = os.path.join(writepath, name +'.evt.py.dow') +down_ep = data.read_neuroscope_intervals(name = 'DOWN', path2file = file) + +file = os.path.join(writepath, name +'.evt.py.upp') +up_ep = data.read_neuroscope_intervals(name = 'UP', path2file = file) + +#%% Compute PETH + +cc2 = nap.compute_eventcorrelogram(spikes, nap.Tsd(up_ep['start'].values), binsize = 0.005, windowsize = 0.255, ep = up_ep, norm = True) +tmp = pd.DataFrame(cc2) +tmp = tmp.rolling(window=8, win_type='gaussian',center=True,min_periods=1).mean(std = 2) + +#%% Example units + +plt.figure() +plt.plot(tmp[pyr][23][-0.05:0.2]) +plt.plot(tmp[pyr][50][-0.05:0.2]) +plt.axhline(1, linestyle = '--', color = 'silver') +plt.axvline(0, color = 'k') +plt.yticks([1], ['mean rate']) +plt.xticks([0], ['DU']) +plt.gca().set_box_aspect(1) diff --git a/Analysis/Figure5C.py b/Analysis/Figure5C.py new file mode 100644 index 0000000..e707a93 --- /dev/null +++ b/Analysis/Figure5C.py @@ -0,0 +1,116 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Sep 15 10:01:54 2023 + +@author: dhruv +""" +#import dependencies +import numpy as np +import pandas as pd +import scipy.io +import pynapple as nap +import os +import matplotlib.pyplot as plt +from scipy.stats import kendalltau + +#%% Data organization + +data_directory = '/media/DataDhruv/Dropbox (Peyrache Lab)/Peyrache Lab Team Folder/Data/AdrianPoSub/###AllPoSub' +datasets = np.genfromtxt(os.path.join(data_directory,'dataset_Hor_DM.list'), delimiter = '\n', dtype = str, comments = '#') + +readpath = '/media/DataDhruv/Dropbox (Peyrache Lab)/Peyrache Lab Team Folder/Projects/PoSub-UPstate/Data' +writepath = '/home/dhruv/Code/MehrotraLevenstein_2023/Analysis/OutputFiles' + +uponset = [] +PMR = [] + +for s in datasets: + print(s) + name = s.split('/')[-1] + path = os.path.join(data_directory, s) + rawpath = os.path.join(readpath,s) + +###Loading the data + data = nap.load_session(rawpath, 'neurosuite') + data.load_neurosuite_xml(rawpath) + spikes = data.spikes + +###Load cell types + filepath = os.path.join(path, 'Analysis') + listdir = os.listdir(filepath) + file = [f for f in listdir if 'CellTypes' in f] + celltype = scipy.io.loadmat(os.path.join(filepath,file[0])) + + pyr = [] + interneuron = [] + hd = [] + + for i in range(len(spikes)): + if celltype['ex'][i] == 1 and celltype['gd'][i] == 1: + pyr.append(i) + + for i in range(len(spikes)): + if celltype['fs'][i] == 1 and celltype['gd'][i] == 1: + interneuron.append(i) + + for i in range(len(spikes)): + if celltype['hd'][i] == 1 and celltype['gd'][i] == 1: + hd.append(i) + +###Load UP and DOWN states + file = os.path.join(writepath, name +'.evt.py.dow') + down_ep = data.read_neuroscope_intervals(name = 'DOWN', path2file = file) + + file = os.path.join(writepath, name +'.evt.py.upp') + up_ep = data.read_neuroscope_intervals(name = 'UP', path2file = file) + +#%% Compute PETH + + cc2 = nap.compute_eventcorrelogram(spikes, nap.Tsd(up_ep['start'].values), binsize = 0.005, windowsize = 0.255, ep = up_ep, norm = True) + tmp = pd.DataFrame(cc2) + tmp = tmp.rolling(window=8, win_type='gaussian',center=True,min_periods=1).mean(std = 2) + dd2 = tmp[0:0.155] + + #Only EX cells + ee = dd2[pyr] + +#%% Compute peak-to-mean ratio (PMR) and UP onset delay + + for i in range(len(ee.columns)): + a = np.where(ee.iloc[:,i] > 0.5) + if len(a[0]) > 0: + PMR.append(ee.iloc[:,i].max()) + res = ee.iloc[:,i].index[a] + uponset.append(res[0]) + +#%% Plotting + +binsize = 0.005 #In seconds +(counts,onsetbins,peakbins) = np.histogram2d(uponset, PMR, + bins=[len(np.arange(0, 0.155, binsize))+1, len(np.arange(0, 0.155, binsize)) + 1], + range = [[-0.0025, 0.1575],[0.5, 3.6]]) + +masked_array = np.ma.masked_where(counts == 0, counts) +cmap = plt.cm.viridis +cmap.set_bad(color='white') + +plt.figure() +plt.imshow(masked_array.T, origin='lower', extent = [onsetbins[0], onsetbins[-1], peakbins[0], peakbins[-1]], + aspect = 'auto', cmap = cmap) +plt.colorbar(ticks = [min(counts.flatten()) + 1, max(counts.flatten())], label = '# cells') +plt.xlabel('UP onset delay (s)') +plt.ylabel('PMR') +plt.gca().set_box_aspect(1) + +y_est = np.zeros(len(uponset)) +m, b = np.polyfit(uponset, PMR, 1) +for i in range(len(uponset)): + y_est[i] = m*uponset[i] + +plt.plot(uponset, y_est + b, color = 'r', zorder = 3) +plt.axhline(1, color ='k', linestyle = '--') + +#%% Stats + +corr, p = kendalltau(uponset, PMR) diff --git a/Analysis/Figure5D.py b/Analysis/Figure5D.py new file mode 100644 index 0000000..0f5cfcc --- /dev/null +++ b/Analysis/Figure5D.py @@ -0,0 +1,172 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Sep 15 11:44:08 2023 + +@author: dhruv +""" + +#import dependencies +import numpy as np +import pandas as pd +import scipy.io +import pynapple as nap +import os +import matplotlib.pyplot as plt +from itertools import combinations +from scipy.stats import pearsonr +from sklearn.manifold import Isomap + +#%% Data organization + +data_directory = '/media/DataDhruv/Dropbox (Peyrache Lab)/Peyrache Lab Team Folder/Data/AdrianPoSub/###AllPoSub' +datasets = np.genfromtxt(os.path.join(data_directory,'dataset_Hor_DM.list'), delimiter = '\n', dtype = str, comments = '#') + +readpath = '/media/DataDhruv/Dropbox (Peyrache Lab)/Peyrache Lab Team Folder/Projects/PoSub-UPstate/Data' +writepath = '/home/dhruv/Code/MehrotraLevenstein_2023/Analysis/OutputFiles' + +allPETH = pd.DataFrame() +uponset = [] +PMR = [] + +for s in datasets: + print(s) + name = s.split('/')[-1] + path = os.path.join(data_directory, s) + rawpath = os.path.join(readpath,s) + +###Loading the data + data = nap.load_session(rawpath, 'neurosuite') + data.load_neurosuite_xml(rawpath) + spikes = data.spikes + +###Load cell types + filepath = os.path.join(path, 'Analysis') + listdir = os.listdir(filepath) + + file = [f for f in listdir if 'CellDepth' in f] + celldepth = scipy.io.loadmat(os.path.join(filepath,file[0])) + depth = celldepth['cellDep'] + + + file = [f for f in listdir if 'CellTypes' in f] + celltype = scipy.io.loadmat(os.path.join(filepath,file[0])) + + pyr = [] + interneuron = [] + hd = [] + + for i in range(len(spikes)): + if celltype['ex'][i] == 1 and celltype['gd'][i] == 1: + pyr.append(i) + + for i in range(len(spikes)): + if celltype['fs'][i] == 1 and celltype['gd'][i] == 1: + interneuron.append(i) + + for i in range(len(spikes)): + if celltype['hd'][i] == 1 and celltype['gd'][i] == 1: + hd.append(i) + +###Load UP and DOWN states + file = os.path.join(writepath, name +'.evt.py.dow') + down_ep = data.read_neuroscope_intervals(name = 'DOWN', path2file = file) + + file = os.path.join(writepath, name +'.evt.py.upp') + up_ep = data.read_neuroscope_intervals(name = 'UP', path2file = file) + +#%% Compute PETH + + cc2 = nap.compute_eventcorrelogram(spikes, nap.Tsd(up_ep['start'].values), binsize = 0.005, windowsize = 0.255, ep = up_ep, norm = True) + tmp = pd.DataFrame(cc2) + tmp = tmp.rolling(window=8, win_type='gaussian',center=True,min_periods=1).mean(std = 2) + dd2 = tmp[0:0.155] + + #Only EX cells + ee = dd2[pyr] + +#%% Concatenate all PETHs for Isomap + + if len(ee.columns) > 0: + + tokeep = [] + + for i in range(len(ee.columns)): + a = np.where(ee.iloc[:,i] > 0.5) + if len(a[0]) > 0: + tokeep.append(ee.columns[i]) + PMR.append(ee.iloc[:,i].max()) + res = ee.iloc[:,i].index[a] + uponset.append(res[0]) + + allPETH = pd.concat([allPETH, ee[tokeep]], axis = 1) + +#%% Isomap + +projection = Isomap(n_components = 2, n_neighbors = 50).fit_transform(allPETH.T.values) +H = PMR/(max(PMR)) + +norm_onset = uponset/max(uponset) +cmap = plt.cm.OrRd + +plt.figure(figsize = (8,8)) +plt.scatter(projection[:,0], projection[:,1], c = cmap(H)) +plt.gca().set_box_aspect(1) +plt.xlabel('dim 1') +plt.ylabel('dim 2') + +#%% Specific examples from Isomap + +allPETH.columns = range(allPETH.columns.size) + +summ = {} +summ['peth'] = allPETH +summ['p1'] = projection[:,0] +summ['p2'] = projection[:,1] + +examples = [2, 208, 226, 378, 528, 828, 1018, 1061, 1081] + +for i in examples: + + plt.figure() + plt.title('Example ' + str(i)) + plt.plot(summ['peth'][i]) + plt.axhline(y = 1, linestyle = '--', color = 'k') + plt.xlabel('Time from DU (ms)') + plt.ylabel('Norm. rate') + +#%% Plot of PMR as a function of Isomap dimension 1 + +r, p = pearsonr(projection[:,0], PMR) + +plt.figure() +plt.scatter(projection[:,0], PMR, label = 'r = ' + str(round(r,2))) +plt.gca().set_box_aspect(1) +plt.xlabel('dim 1') +plt.ylabel('PMR') +plt.legend(loc = 'upper right') + +#%% Gradient vector + +pairs = list(combinations(summ['peth'].columns, 2)) + +F_pmrr = [] + +for i,p in enumerate(pairs): + diff_pmrr = PMR[p[0]] - PMR[p[1]] + + dx = summ['p1'][p[0]] - summ['p1'][p[1]] + dy = summ['p2'][p[0]] - summ['p2'][p[1]] + + F_pmrr.append([diff_pmrr/dx, diff_pmrr/dy]) + +mean_pmrr = np.mean(F_pmrr, axis = 0) + +origin = np.array([[0, 0], [0, 0]]) +plt.figure() +plt.title('Gradient vector') +plt.xlim(0, 0.4) +plt.ylim(-0.4, 0) +plt.xlabel('dim 1') +plt.ylabel('dim 2') +plt.quiver(origin[0], origin[1], mean_pmrr[0] , mean_pmrr[1], angles = 'xy', scale_units = 'xy', scale = 1) \ No newline at end of file diff --git a/Analysis/Figure5E.py b/Analysis/Figure5E.py new file mode 100644 index 0000000..f9fc847 --- /dev/null +++ b/Analysis/Figure5E.py @@ -0,0 +1,121 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Sep 15 10:21:30 2023 + +@author: dhruv +""" +#import dependencies +import numpy as np +import pandas as pd +import scipy.io +import pynapple as nap +import os +import matplotlib.pyplot as plt +from scipy.stats import kendalltau + +#%% Data organization + +data_directory = '/media/DataDhruv/Recordings/WatsonBO' +datasets = np.genfromtxt(os.path.join(data_directory,'dataset_DM.list'), delimiter = '\n', dtype = str, comments = '#') + +uponset = [] +PMR = [] + +for s in datasets: + print(s) + name = s.split('/')[-1] + path = os.path.join(data_directory, s) + +#%% Load the relevant mat files + + listdir = os.listdir(path) + file = [f for f in listdir if 'spikes' in f] + spikedata = scipy.io.loadmat(os.path.join(path,file[0])) + + file = [f for f in listdir if 'events' in f] + events = scipy.io.loadmat(os.path.join(path,file[0])) + + file = [f for f in listdir if 'states' in f] + states = scipy.io.loadmat(os.path.join(path,file[0])) + + file = [f for f in listdir if 'CellClass' in f] + cellinfo = scipy.io.loadmat(os.path.join(path,file[0])) + +#%% Parsing mat files for variables of interest + +###Load EX cells + ex = cellinfo['CellClass'][0][0][1][0] + pyr = [] + + for i in range(len(ex)): + if ex[i] == 1: + pyr.append(i) + +###Load spikes (only for EX cells) + spks = spikedata['spikes'] + spk = {} + for i in range(len(spks[0][0][1][0])): + spk[i] = nap.Ts(spks[0][0][1][0][i]) + spikes = nap.TsGroup(spk) + spikes = spikes[pyr] + +###Load sleep states + sleepstate = states['SleepState'] + wake_ep = nap.IntervalSet(start = sleepstate[0][0][0][0][0][0][:,0], end = sleepstate[0][0][0][0][0][0][:,1]) + nrem_ep = nap.IntervalSet(start = sleepstate[0][0][0][0][0][1][:,0], end = sleepstate[0][0][0][0][0][1][:,1]) + rem_ep = nap.IntervalSet(start = sleepstate[0][0][0][0][0][2][:,0], end = sleepstate[0][0][0][0][0][2][:,1]) + +###Load UP and DOWN states + slowwaves = events['SlowWaves'] + up_ep = nap.IntervalSet( start = slowwaves[0][0][2][0][0][0][:,0], end = slowwaves[0][0][2][0][0][0][:,1]) + down_ep = nap.IntervalSet( start = slowwaves[0][0][2][0][0][1][:,0], end = slowwaves[0][0][2][0][0][1][:,1]) + +#%% Compute PETH + + cc = nap.compute_eventcorrelogram(spikes, nap.Tsd(up_ep['start'].values), binsize = 0.005, windowsize = 0.255, ep = up_ep, norm = True) + tmp = pd.DataFrame(cc) + tmp = tmp.rolling(window=8, win_type='gaussian',center=True,min_periods=1).mean(std = 2) + dd = tmp[0:0.155] + +#%% Compute peak-to-mean ratio (PMR) and UP onset delay + + for i in range(len(dd.columns)): + a = np.where(dd.iloc[:,i] > 0.5) + if len(a[0]) > 0: + PMR.append(dd.iloc[:,i].max()) + res = dd.iloc[:,i].index[a] + uponset.append(res[0]) + +#%% Plotting + +binsize = 0.005 #In seconds +(counts,onsetbins,peakbins) = np.histogram2d(uponset, PMR, + bins = [len(np.arange(0, 0.155, binsize)) + 1, len(np.arange(0, 0.155, binsize)) + 1], + range = [[-0.0025, 0.1575],[0.5, 3.6]]) + +masked_array = np.ma.masked_where(counts == 0, counts) +cmap = plt.cm.viridis +cmap.set_bad(color='white') + +plt.figure() +plt.imshow(masked_array.T, origin='lower', extent = [onsetbins[0],onsetbins[-1],peakbins[0],peakbins[-1]], + aspect='auto', cmap = cmap) +plt.colorbar(ticks = [min(counts.flatten()) + 1, max(counts.flatten())], label = '# cells') +plt.xlabel('UP onset delay (s)') +plt.ylabel('PMR') +plt.gca().set_box_aspect(1) + +y_est = np.zeros(len(uponset)) +m, b = np.polyfit(uponset, PMR, 1) +for i in range(len(uponset)): + y_est[i] = m*uponset[i] + +plt.plot(uponset, y_est + b, color = 'r', zorder = 3) +plt.axhline(1, color ='k', linestyle = '--') + +#%% Stats + +corr, p = kendalltau(uponset, PMR) + + \ No newline at end of file diff --git a/Analysis/Figure5F.py b/Analysis/Figure5F.py new file mode 100644 index 0000000..8f524dd --- /dev/null +++ b/Analysis/Figure5F.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Sep 15 10:45:53 2023 + +@author: dhruv +""" + +#import dependencies +import numpy as np +import pandas as pd +import pynapple as nap +import os +import matplotlib.pyplot as plt +from scipy.stats import kendalltau + +#%% Data organization + +data_directory = '/media/dhruv/LaCie1/A7800' +datasets = np.genfromtxt(os.path.join(data_directory,'dataset_DM.list'), delimiter = '\n', dtype = str, comments = '#') + +PMR = [] +uponset = [] + +for s in datasets: + print(s) + name = s.split('/')[-1] + path = os.path.join(data_directory, s) + +###Load the data + data = nap.load_session(path, 'neurosuite') + data.load_neurosuite_xml(path) + spikes = data.spikes + +###Load UP and DOWN states + file = os.path.join(path, name +'.evt.py.dow') + down_ep = data.read_neuroscope_intervals(name = 'DOWN', path2file = file) + + file = os.path.join(path, name +'.evt.py.upp') + up_ep = data.read_neuroscope_intervals(name = 'UP', path2file = file) + +#%% Compute PETH + + cc2 = nap.compute_eventcorrelogram(spikes, nap.Tsd(up_ep['start'].values), binsize = 0.005, windowsize = 0.255, ep = up_ep, norm = True) + tmp = pd.DataFrame(cc2) + tmp = tmp.rolling(window=8, win_type='gaussian',center=True,min_periods=1).mean(std = 2) + dd2 = tmp[0:0.155] + +#%% Compute peak-to-mean ratio (PMR) and UP onset delay + + for i in range(len(dd2.columns)): + a = np.where(dd2.iloc[:,i] > 0.5) + if len(a[0]) > 0: + PMR.append(dd2.iloc[:,i].max()) + res = dd2.iloc[:,i].index[a] + uponset.append(res[0]) + +#%% Plotting + +binsize = 0.005 #In seconds +(counts,onsetbins,peakbins) = np.histogram2d(uponset, PMR, + bins = [len(np.arange(0, 0.155, binsize)) + 1, len(np.arange(0, 0.155, binsize)) + 1], + range = [[-0.0025, 0.1575],[0.5, 3.6]]) + +masked_array = np.ma.masked_where(counts == 0, counts) +cmap = plt.cm.viridis +cmap.set_bad(color='white') + +plt.figure() +plt.imshow(masked_array.T, origin='lower', extent = [onsetbins[0],onsetbins[-1],peakbins[0],peakbins[-1]], + aspect='auto', cmap = cmap) +plt.colorbar(ticks = [min(counts.flatten()) + 1, max(counts.flatten())], label = '# cells') +plt.xlabel('UP onset delay (s)') +plt.ylabel('PMR') +plt.gca().set_box_aspect(1) + +y_est = np.zeros(len(uponset)) +m, b = np.polyfit(uponset, PMR, 1) +for i in range(len(uponset)): + y_est[i] = m*uponset[i] + +plt.plot(uponset, y_est + b, color = 'r', zorder = 3) +plt.axhline(1, color ='k', linestyle = '--') + +#%% Stats + +corr, p = kendalltau(uponset, PMR) + \ No newline at end of file diff --git a/Analysis/detectDOWNs.py b/Analysis/detectDOWNs_HDC.py similarity index 100% rename from Analysis/detectDOWNs.py rename to Analysis/detectDOWNs_HDC.py diff --git a/Analysis/detectDOWNs_RSC.py b/Analysis/detectDOWNs_RSC.py new file mode 100644 index 0000000..4984bcc --- /dev/null +++ b/Analysis/detectDOWNs_RSC.py @@ -0,0 +1,74 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Sep 15 11:16:31 2023 + +@author: dhruv +""" + +#import dependencies +import numpy as np +import os +import pynapple as nap +import matplotlib.pyplot as plt + +#%% Data organization + +data_directory = '/media/dhruv/LaCie1/A7800' +datasets = np.genfromtxt(os.path.join(data_directory,'dataset_DM.list'), delimiter = '\n', dtype = str, comments = '#') + +for s in datasets: + print(s) + name = s.split('/')[-1] + path = os.path.join(data_directory, s) + +###Loading the data + data = nap.load_session(path, 'neurosuite') + data.load_neurosuite_xml(path) + spikes = data.spikes + +#%% Fishing out wake and sleep epochs + + epochs = data.epochs + sleep_ep = epochs['sleep'] + wake_ep = epochs['wake'] + + file = os.path.join(path, name +'.sws.evt') + sws_ep = data.read_neuroscope_intervals(name = 'SWS', path2file = file) + +#%% Detection of UP and DOWN states + + bin_size = 0.01 #In seconds + smoothing_window = 0.02 #In seconds + + #Create binned spike times + rates = spikes.count(bin_size, sws_ep) + + #Smooth binned spike times + total2 = rates.as_dataframe().rolling(window = 100, win_type = 'gaussian', + center = True, min_periods = 1, + axis = 0).mean(std = int(smoothing_window/bin_size)) + + #Create population rate/multi-unit activity (MUA) + total2 = total2.sum(axis = 1) + total2 = nap.Tsd(total2) + + #Apply the threshold + idx = total2.threshold(np.percentile(total2.values,20),'below') + + #Exclusion criteria + down_ep = idx.time_support + down_ep = nap.IntervalSet(start = down_ep['start'], end = down_ep['end']) + down_ep = down_ep.drop_short_intervals(bin_size) + down_ep = down_ep.merge_close_intervals(bin_size*2) + down_ep = down_ep.drop_short_intervals(bin_size*3) + down_ep = down_ep.drop_long_intervals(bin_size*50) + + #Create the UP state + up_ep = nap.IntervalSet(down_ep['end'][0:-1], down_ep['start'][1:]) + down_ep = sws_ep.intersect(down_ep) + +#%% Write the data as an event file + + data.write_neuroscope_intervals(extension = '.evt.py.dow', isets = down_ep, name = 'DOWN') + data.write_neuroscope_intervals(extension = '.evt.py.upp', isets = up_ep, name = 'UP') \ No newline at end of file