GluN2A Enables Noradrenergic Control of Prefrontal Oscillations and Cognitive Flexibility.
The 2 matches
- [1] § Materials and Methods › Multielectrode array (MEA) recordings ↔ MEA LFP _ Open OnDemand_LC_mPFC.ipynb, lines 375–439 · score 0.69 · power spectral density, Absolute power, Multitaper, windows, band, gamma
- [2] § Materials and Methods › Multielectrode array (MEA) recordings ↔ MEA LFP _ Open OnDemand_NE_Pharmacology.ipynb, lines 374–438 · score 0.69 · power spectral density, Absolute power, Multitaper, windows, band, gamma
Paper
Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC
The paper is loaded when this pane is shown.
The authors' code
Jupyter notebook · 1,495 lines · 52 KB · no license · 1 match
- # %%
- from scipy.signal import butter, lfilter
- from scipy import signal
- import numpy as np
- import matplotlib.pyplot as plt
- import seaborn as sns
- import sklearn as sk
- import h5py
- import os.path
- import pandas as pd
- from sklearn.decomposition import PCA
- from sklearn.cluster import KMeans
- from sklearn.cluster import MiniBatchKMeans
- from sklearn.metrics import silhouette_score
- from sklearn.cluster import DBSCAN
- from sklearn import metrics
- import librosa.display
- # from kneed import KneeLocator
- import os
- import urllib
- import numpy as np
- from scipy.io import loadmat
- from tensorpac import Pac, EventRelatedPac, PreferredPhase
- from tensorpac.utils import PeakLockedTF, PSD, ITC, BinAmplitude
- from tensorpac.signals import pac_signals_wavelet
- import matplotlib.pyplot as plt
- plt.style.use('seaborn-poster')
- # %%
- save_path = "/home/hosseini/Desktop/2024/07.08.2024/Preprocessed/Output"
- save_pathh = "/home/hosseini/Desktop/2024/07.08.2024/Preprocessed/Output/mPFC_Slice2_Multitaper"
- save_path_coherence = "/home/hosseini/Desktop/2024/07.08.2024/Preprocessed/Output/mPFC_coherence_Slice2"
- name = '\_07082024_mPFC_Slice2'
- name2 = str('_07082024_mPFC_Slice2')
- strain = 'WT'
- # filenames
- data_path = "/home/hosseini/Desktop/2024/07.08.2024/Preprocessed"
- file_xtrain_vehicle = data_path + r"/2024-07-08mPFC_Slice2_LC_to_mPFC_CAG_ChR2_eYFP_10Second_Stim_A_Baseline__0002.h5" #### Vehicle
- file_xtrain_GFC = data_path + r"/2024-07-08mPFC_Slice2_LC_to_mPFC_CAG_ChR2_eYFP_10Second_Stim_A_Stimulation_.h5" #### NE
- file_xtrain_GFC_DOB = data_path + r"/2024-07-08mPFC_Slice2_LC_to_mPFC_CAG_ChR2_eYFP_10Second_Stim_A_Yohimbine+Stimulation__0002.h5" #### NE + DRUG
- # file_xtrain_vehicle2 = data_path + r"/2024-06-17__Slice2_WT_mPFC_NE2__0001.h5" #### NE2
- ################################################################### Baseline
- with h5py.File(file_xtrain_vehicle, 'r') as hdf:
- dt1 = hdf.get('Data/Recording_0/AnalogStream/Stream_0')
- dt1_items = list(dt1.items())
- ChannelData_vehicle = np.array(dt1.get('ChannelData'))
- print(ChannelData_vehicle.shape)
- ################################################################### Kanaite
- with h5py.File(file_xtrain_GFC, 'r') as hdf:
- dt1 = hdf.get('Data/Recording_0/AnalogStream/Stream_0')
- dt1_items = list(dt1.items())
- ChannelData_GFC = np.array(dt1.get('ChannelData'))
- print(ChannelData_GFC.shape)
- ################################################################### Kanaite + Bicuculline
- with h5py.File(file_xtrain_GFC_DOB, 'r') as hdf:
- dt1 = hdf.get('Data/Recording_0/AnalogStream/Stream_0')
- dt1_items = list(dt1.items())
- ChannelData_GFC_DOB = np.array(dt1.get('ChannelData'))
- print(ChannelData_GFC_DOB.shape)
- sf = 20000 # sampling requency
- print(f'time recorded in min: {((len(ChannelData_vehicle[1])/sf)/60)}')
- # %%
- # dim_x=60 # number of electrodes
- # dim_y=2401980 # time in seconds with fs 20k
- data_concat = np.concatenate((ChannelData_vehicle[:,:],
- ChannelData_GFC[:,:],
- ChannelData_GFC_DOB[:,:],
- # ChannelData_vehicle2[:,:],
- ), axis=1)
- np.shape(data_concat)
- # %% [markdown]
- # # Lowpass filter & downsample the data
- # %%
- def filter_data(data, low, high, sf, order=2):
- # Determine Nyquist frequency
- nyq = sf/2
- # Set bands
- low = low/nyq
- high = high/nyq
- # Calculate coefficients
- b, a = butter(order, [low, high], btype='bandpass',)
- # Filter signal
- filtered_data = lfilter(b, a, data)
- return filtered_data
- def downsample_data(data, sf, target_sf):
- # Calculate the initial downsampling factor
- factor = int(sf / target_sf)
- # Check if downsampling is necessary
- if factor < 1:
- raise ValueError("Target sampling frequency must be less than original sampling frequency.")
- # Downsample using multiple steps if necessary
- data_down = data
- while factor > 1:
- # Apply decimation in chunks to avoid excessive decimation
- step_factor = min(factor, 10) # Decimate in steps of max 10 to avoid aliasing
- data_down = signal.decimate(data_down, step_factor, ftype='fir')
- sf = sf / step_factor # Update the sampling rate after each decimation
- factor = int(sf / target_sf) # Recalculate the factor for the next iteration
- return data_down, int(sf)
- # %% [markdown]
- # # Apply the eprevious steps (notch 60, filter, downsample) to Combined data
- # %%
- from scipy import signal
- from scipy.signal import butter, lfilter
- sf = 20000 # sampling frequency
- target_sf = 2000 # target frequency for downsampling
- f0 = 60.0 # Frequency to be removed from signal (Hz)
- w0 = f0/(target_sf/2) # Normalized Frequency
- Q = 30 # Quality factor
- # Design notch filter
- b, a = signal.iirnotch(w0, Q)
- start = 0
- end = int(len(data_concat[0])/sf)
- ################################################### all data comined
- clean_data = []
- lfp_power = []
- # clean_data_power = []
- for chanlNo in range(len(data_concat)):
- data = data_concat[chanlNo].ravel()
- # compare the raw data with the filtered signal
- spike_data = filter_data(data[start*sf:end*sf], low=4, high=400, sf=sf) # Band pass filter the data for LFP analysis
- # Now we use the above function to downsample the signal.
- lfp_data, sf_lfp = downsample_data(spike_data, sf=sf, target_sf=target_sf)
- lfp_power.append(signal.welch(lfp_data, fs=sf_lfp, window='hamming', nperseg=1024, scaling='spectrum'))
- clean_data.append(signal.lfilter(b, a, lfp_data))
- # %% [markdown]
- # # Multitaper Spectrogram
- # %%
- %run /home/hosseini/Desktop/multitaper_spectrogram.py
- # %%
- # from multitaper_spectrogram_python import multitaper_spectrogram # import multitaper_spectrogram function from the multitaper_spectrogram_python.py file
- # Set spectrogram params
- fs = 2000 # Sampling Frequency
- frequency_range = [1, 250] # Limit frequencies from 0 to 25 Hz
- time_bandwidth = 5 # Set time-half bandwidth
- num_tapers = 8 # Set number of tapers (optimal is time_bandwidth*2 - 1)
- window_params = [4, 1] # Window size is 4s with step size of 1s
- min_nfft = 0 # No minimum nfft
- detrend_opt = 'constant' # detrend each window by subtracting the average
- multiprocess = True # use multiprocessing
- n_jobs = -1 # use 3 cores in multiprocessing
- weighting = 'unity' # weight each taper at 1
- plot_on = True # plot spectrogram
- clim_scale = False # do not auto-scale colormap
- verbose = False # print extra info
- xyflip = False # do not transpose spect output matrix
- # %%
- import warnings
- warnings.filterwarnings('ignore')
- import os.path
- #first we filter the raw data and then downsample the raw data
- from matplotlib.backends.backend_pdf import PdfPages
- pdf_pages = PdfPages(os.path.join(save_pathh,'Multitaper.pdf'))
- for i in range (len(clean_data)):
- spect, stimes, sfreqs = multitaper_spectrogram(clean_data[i],
- fs,
- frequency_range,
- time_bandwidth,
- num_tapers,
- window_params, min_nfft, detrend_opt,
- multiprocess, n_jobs,
- weighting,
- # plot_on,
- clim_scale,
- verbose, xyflip)
- # For saving the graphs
- spect_data = spect
- clim = np.percentile(spect_data, [5, 95]) # Scale colormap from 5th percentile to 95th
- # freqs, times, spectrogram = signal.spectrogram(data)#######
- fig = plt.figure(figsize=(10, 5))
- librosa.display.specshow(librosa.amplitude_to_db(spect, ref=np.max), x_axis='time', y_axis='linear',
- x_coords=stimes, y_coords=sfreqs, shading='auto',
- cmap= "jet"
- )
- plt.colorbar(label='Power [dB]')
- plt.xlabel("Time [MM:SS]")
- plt.ylabel("Frequency [Hz]")
- if clim_scale:
- plt.clim(clim) # actually change colorbar scale
- plt.title(f'Electrod No. : {i}')
- plt.axhline(y = 4, color = 'w', linestyle = '--', linewidth= 1.0 )
- plt.axhline(y = 10, color = 'w', linestyle = '--', linewidth= 1.0 )
- plt.axhline(y = 30, color = 'w', linestyle = '--', linewidth= 1.0 )
- plt.axvline(x = 5*60, color = 'w', linestyle = '--', linewidth= 2.0 ) # 5 mins
- plt.axvline(x = 10*60, color = 'w', linestyle = '--', linewidth= 2.0 )
- # plt.axvline(x = 15*60, color = 'w', linestyle = '--', linewidth= 2.0 )
- pdf_pages.savefig(fig)
- ######## Saving
- # plt.savefig(os.path.join(save_pathh , 'ElectrodeNo{}.svg'.format(i)),bbox_inches = 'tight')
- pdf_pages.close()
- # %% [markdown]
- # # STOP! DELETE UNWANTED CHANNELS
- # %%
- unwanted_num = [0,1,2,3,4,5,6,7,10,11,12,13,15,19,20,23,27,32,35,38,51,54,
- ]# channels to be removed
- # List of numbers with corresponding indices as shown in the image
- numbers_dict = {
- 54: 0, 56: 1, 55: 2, 46: 3, 45: 4, 44: 5, 36: 6, 35: 7, 34: 8, 26: 9,
- 16: 10, 25: 11, 15: 12, 24: 13, 14: 14, 13: 15, 23: 16, 12: 17, 22: 18,
- 11: 19, 21: 20, 33: 21, 32: 22, 31: 23, 43: 24, 42: 25, 41: 26, 52: 27,
- 51: 28, 53: 29, 63: 30, 61: 31, 62: 32, 71: 33, 72: 34, 73: 35, 81: 36,
- 82: 37, 83: 38, 91: 39, 101: 40, 92: 41, 102: 42, 93: 43, 103: 44,
- 104: 45, 94: 46, 105: 47, 95: 48, 106: 49, 96: 50, 84: 51, 85: 52,
- 86: 53, 74: 54, 75: 55, 76: 56, 65: 57, 66: 58, 64: 59
- }
- # Function to get corresponding indices, excluding unwanted numbers
- def get_corresponding_values(selected_nums):
- filtered_values = [num for num in selected_nums if num in unwanted_num]
- return filtered_values
- # L2/3
- L_2_3 = get_corresponding_values([
- 19, 20, 23, 26, 28, 31, 33, 36, 39, 40,
- # 10, 9, 6, 3, 1, 58, 56, 53, 50, 49
- ])
- # L5
- L_5 = get_corresponding_values([
- # 19, 20, 23, 26, 28, 31, 33, 36, 39, 40,
- 17, 18, 22, 25, 27, 32, 34, 41, 42,
- 15, 16, 24, 29, 30, 35, 38, 43, 44,
- 13, 8, 0, 59, 54, 51, 46, 45, 12,
- 12, 11, 7, 4, 2, 57, 55, 52, 48, 47,
- 10, 9, 6, 3, 1, 58, 56, 53, 50, 49
- ]
- )
- print(L_2_3, L_5)
- # %%
- layer ='Layer2_3'
- import warnings
- warnings.filterwarnings('ignore')
- import os.path
- #first we filter the raw data and then downsample the raw data
- from matplotlib.backends.backend_pdf import PdfPages
- pdf_pages = PdfPages(os.path.join(save_pathh,'Multitaper_Layer2_3.pdf'))
- for i in L_2_3:
- spect, stimes, sfreqs = multitaper_spectrogram(clean_data[i],
- fs,
- frequency_range,
- time_bandwidth,
- num_tapers,
- window_params, min_nfft, detrend_opt,
- multiprocess, n_jobs,
- weighting,
- # plot_on,
- clim_scale,
- verbose, xyflip)
- # For saving the graphs
- spect_data = spect
- clim = np.percentile(spect_data, [5, 95]) # Scale colormap from 5th percentile to 95th
- # freqs, times, spectrogram = signal.spectrogram(data)#######
- fig = plt.figure(figsize=(10, 5))
- librosa.display.specshow(librosa.amplitude_to_db(spect, ref=np.max), x_axis='time', y_axis='linear',
- x_coords=stimes, y_coords=sfreqs, shading='auto',
- cmap= "jet"
- )
- # plt.colorbar(label='Power [dB]')
- plt.xlabel("Time [MM:SS]")
- plt.ylabel("Frequency [Hz]")
- plt.title(f'Electrod No. : {i}')
- plt.axhline(y = 4, color = 'w', linestyle = '--', linewidth= 2.0 )
- plt.axhline(y = 10, color = 'w', linestyle = '--', linewidth= 2.0 )
- plt.axhline(y = 25, color = 'w', linestyle = '--', linewidth= 2.0 )
- plt.axvline(x = 5*60 ,color = 'w', linestyle = '--', linewidth= 2.0 ) #
- plt.axvline(x = 10*60 ,color = 'w', linestyle = '--', linewidth= 2.0 ) #
- pdf_pages.savefig(fig)
- ######## Saving
- plt.savefig(os.path.join(save_pathh , 'ElectrodeNo{}_Layer2_3.svg'.format(i)),bbox_inches = 'tight')
- pdf_pages.close()
- # %%
- layer ='Layer5'
- import warnings
- warnings.filterwarnings('ignore')
- import os.path
- #first we filter the raw data and then downsample the raw data
- from matplotlib.backends.backend_pdf import PdfPages
- pdf_pages = PdfPages(os.path.join(save_pathh,'Multitaper_Layer5.pdf'))
- for i in L_5:
- spect, stimes, sfreqs = multitaper_spectrogram(clean_data[i],
- fs,
- frequency_range,
- time_bandwidth,
- num_tapers,
- window_params, min_nfft, detrend_opt,
- multiprocess, n_jobs,
- weighting,
- # plot_on,
- clim_scale,
- verbose, xyflip)
- # For saving the graphs
- spect_data = spect
- clim = np.percentile(spect_data, [5, 95]) # Scale colormap from 5th percentile to 95th
- # freqs, times, spectrogram = signal.spectrogram(data)#######
- fig = plt.figure(figsize=(10, 5))
- librosa.display.specshow(librosa.amplitude_to_db(spect, ref=np.max), x_axis='time', y_axis='linear',
- x_coords=stimes, y_coords=sfreqs, shading='auto',
- cmap= "jet"
- )
- # plt.colorbar(label='Power [dB]')
- plt.xlabel("Time [MM:SS]")
- plt.ylabel("Frequency [Hz]")
- plt.title(f'Electrod No. : {i}')
- plt.axhline(y = 4, color = 'w', linestyle = '--', linewidth= 2.0 )
- plt.axhline(y = 10, color = 'w', linestyle = '--', linewidth= 2.0 )
- plt.axhline(y = 25, color = 'w', linestyle = '--', linewidth= 2.0 )
- plt.axvline(x = 5*60 ,color = 'w', linestyle = '--', linewidth= 2.0 ) #
- plt.axvline(x = 10*60 ,color = 'w', linestyle = '--', linewidth= 2.0 ) #
- pdf_pages.savefig(fig)
- ######## Saving
- plt.savefig(os.path.join(save_pathh , 'ElectrodeNo{}_Layer5.svg'.format(i)),bbox_inches = 'tight')
- pdf_pages.close()
- # %%
- # %%
- # delta = [1, 4]
- # theta = [4, 10]
- # beta = [10, 30]
- # gamma = [30, 80]
- # ref = https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4009705/
- def bandpower(data, sf, band, method='welch', window_sec=None, relative=False):
- """Compute the average power of the signal x in a specific frequency band.
- Requires MNE-Python >= 0.14.
- Parameters
- ----------
- data : 1d-array
- Input signal in the time-domain.
- sf : float
- Sampling frequency of the data.
- band : list
- Lower and upper frequencies of the band of interest.
- method : string
- Periodogram method: 'welch' or 'multitaper'
- window_sec : float
- Length of each window in seconds. Useful only if method == 'welch'.
- If None, window_sec = (1 / min(band)) * 2.
- relative : boolean
- If True, return the relative power (= divided by the total power of the signal).
- If False (default), return the absolute power.
- Return
- ------
- bp : float
- Absolute or relative band power.
- """
- from scipy.signal import welch
- from scipy.integrate import simps
- from mne.time_frequency import psd_array_multitaper
- band = np.asarray(band)
- low, high = band
- # Compute the modified periodogram (Welch)
- if method == 'welch':
- if window_sec is not None:
- nperseg = window_sec * sf
- else:
- nperseg = (2 / low) * sf
- freqs, psd = welch(data, sf, nperseg=nperseg)
- elif method == 'multitaper':
- psd, freqs = psd_array_multitaper(data, sf, adaptive=False,
- normalization='full', verbose=0, n_jobs=-1)
- # Frequency resolution
- freq_res = freqs[1] - freqs[0]
- # Find index of band in frequency vector
- idx_band = np.logical_and(freqs >= low, freqs <= high)
- # Integral approximation of the spectrum using parabola (Simpson's rule)
- bp = simps(psd[idx_band], dx=freq_res)
- if relative:
- bp /= simps(psd, dx=freq_res)
- return bp
- # %% [markdown]
- # # Layer 2/3
- # %%
- from scipy import signal
- from scipy.signal import butter, lfilter
- sf = 20000 # sampling frequency
- target_sf = 1000 # target frequency for downsampling
- f0 = 60.0 # Frequency to be removed from signal (Hz)
- w0 = f0/(target_sf/2) # Normalized Frequency
- Q = 30 # Quality factor
- # Design notch filter
- b, a = signal.iirnotch(w0, Q)
- start = 0
- end = int(len(ChannelData_vehicle[0])/sf)
- ################################################### baseline
- clean_data_vehicle = []
- for chanlNo in (L_2_3):
- data = ChannelData_vehicle[chanlNo].ravel()
- spike_data = filter_data(data[start*sf:end*sf], low=1, high= 500, sf=sf) # Band pass filter the data for LFP analysis
- # Now we use the above function to downsample the signal.
- lfp_data, sf_lfp = downsample_data(spike_data, sf=sf, target_sf=target_sf)
- clean_data_vehicle.append(signal.lfilter(b, a, lfp_data))
- ################################################### kainate
- end = int(len(ChannelData_GFC[0])/sf)
- clean_data_GFC = []
- for chanlNo in (L_2_3):
- data = ChannelData_GFC[chanlNo].ravel()
- spike_data = filter_data(data[start*sf:end*sf], low=1, high= 500, sf=sf) # Band pass filter the data for LFP analysis
- # Now we use the above function to downsample the signal.
- lfp_data, sf_lfp = downsample_data(spike_data, sf=sf, target_sf=target_sf)
- clean_data_GFC.append(signal.lfilter(b, a, lfp_data))
- # ################################################### kainate + bicuculine
- end = int(len(ChannelData_GFC_DOB[0])/sf)
- clean_data_GFC_DOB = []
- for chanlNo in (L_2_3):
- data = ChannelData_GFC_DOB[chanlNo].ravel()
- spike_data = filter_data(data[start*sf:end*sf], low=1, high= 500, sf=sf) # Band pass filter the data for LFP analysis
- # Now we use the above function to downsample the signal.
- lfp_data, sf_lfp = downsample_data(spike_data, sf=sf, target_sf=target_sf)
- clean_data_GFC_DOB.append(signal.lfilter(b, a, lfp_data))
- # %% [markdown]
- # # stop ! Change the NAME of the drug
- # %%
- import warnings
- warnings.filterwarnings('ignore')
- import numpy as np
- import pandas as pd
- import os
- layer = 'L_2_3'
- # Frequency bands and data groupings
- freq_bands = {
- 'delta': [1, 4],
- 'theta': [4, 10],
- 'beta': [10, 25],
- 'lowgamma': [30, 60],
- 'highgamma': [60, 100],
- 'gamma': [30, 100],
- 'ripple': [100, 200]
- }
- data_groups = {
- 'vehicle': clean_data_vehicle,
- 'NE ': clean_data_GFC,
- 'NE_YHM': clean_data_GFC_DOB,
- }
- # Time segments
- time_segments = {
- '_immediate': (38, 43),
- '_delayed': (38, 98),
- }
- # Function to compute band power for each time segment
- def compute_band_power(data, fs, freq_band, start, end):
- return bandpower(data[start * fs:end * fs], fs, freq_band, 'welch')
- # Initialize results dictionary for band power per time segment
- band_results = {band: {group: {segment: [] for segment in time_segments} for group in data_groups} for band in freq_bands}
- # Main loop for bandpower calculation with time segments
- for i in range(len(L_2_3)):
- for band_name, band_range in freq_bands.items():
- for group_name, group_data in data_groups.items():
- for segment_name, (start, end) in time_segments.items():
- # Compute band power for the specific time segment
- bp = compute_band_power(group_data[i], target_sf, band_range, start, end)
- band_results[band_name][group_name][segment_name].append(bp)
- # Prepare data for export, now considering time segments
- export_data = {}
- for band in freq_bands:
- concatenated_data = []
- for group in data_groups:
- for segment in time_segments:
- segment_data = band_results[band][group][segment]
- # Ensure the segment data is not empty and contains valid arrays
- if len(segment_data) > 0:
- # Convert any zero-dimensional arrays into one-dimensional ones
- segment_data = [np.atleast_1d(data) for data in segment_data]
- concatenated_data.append(np.concatenate(segment_data))
- # Only concatenate if we have non-empty arrays
- if len(concatenated_data) > 0:
- export_data[band] = np.concatenate(concatenated_data)
- else:
- export_data[band] = np.array([]) # Handle case where there's no data
- # Create conditions list, including the data groups and time segments
- conditions = np.concatenate([
- len(band_results['delta'][group][segment]) * [f"{group}_{segment}"]
- for group in data_groups for segment in time_segments if len(band_results['delta'][group][segment]) > 0
- ])
- # Calculate average band power for each condition
- average_cols = {}
- for band in freq_bands:
- if len(export_data[band]) > 0:
- # Split data for the band into separate arrays by condition (group + segment)
- data_per_condition = np.array_split(export_data[band], len(data_groups) * len(time_segments))
- avg_per_condition = [np.mean(segment) for segment in data_per_condition]
- # Create a column repeating the average value for alignment with the DataFrame length
- average_col = np.concatenate([[avg] * len(segment) for avg, segment in zip(avg_per_condition, data_per_condition)])
- average_cols[f'Average {band}'] = average_col
- else:
- average_cols[f'Average {band}'] = np.array([]) # Handle empty case
- # Metadata (modify or add actual values as needed)
- layer_col = [layer] * len(conditions) # Replace with actual data if available
- id_col = [name] * len(conditions) # Replace with actual data if available
- strain_col = [strain] * len(conditions) # Replace with actual data if available
- # Create DataFrame for export, now including time segments
- df_export = pd.DataFrame({
- **export_data,
- 'condition': conditions,
- **average_cols,
- 'Layer': layer_col,
- 'ID': id_col,
- 'strain': strain_col
- })
- # Export to Excel
- file_name = f"{name}_BandPower_{layer}.xlsx"
- save_path_full = os.path.join(save_path, file_name)
- with pd.ExcelWriter(save_path_full, engine='xlsxwriter') as writer:
- df_export.to_excel(writer, sheet_name='BandPower', index=False)
- print("Data exported successfully.")
- # %%
- import seaborn as sns
- import matplotlib.pyplot as plt
- # Load the exported data from the previous code
- df = pd.read_excel(os.path.join(save_path, file_name))
- # Frequency bands and titles
- bands = ['delta', 'theta', 'beta', 'lowgamma', 'highgamma', 'ripple']
- titles = [
- r'[$\delta$] Power', r'[$\theta$] Power', r'[$\beta$] Power',
- r'low [$\gamma$] Power', r'high [$\gamma$] Power', r'Cortical Ripple Power'
- ]
- # Updated condition labels based on time segments and groups
- xticklabels = [
- 'Vehicle_immediate', 'Vehicle_delayed',
- 'NE_immediate', 'NE_delayed',
- 'NE_YHM_immediate', 'NE_YHM_delayed'
- ]
- # Create subplots for each frequency band
- fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(15, 5))
- # Loop over each frequency band and create a plot
- for i, band in enumerate(bands):
- row, col = divmod(i, 3) # Determine row and column for each subplot
- # Create boxplot for each condition and frequency band
- sns.boxplot(data=df, x='condition', y=band, ax=ax[row, col], showfliers=False)
- # Create swarmplot (or stripplot) to show individual points on top of the boxplot
- sns.stripplot(data=df, x='condition', y=band, color='k', ax=ax[row, col], dodge=True)
- # Set the x-axis labels to the updated condition names
- ax[row, col].set_xticklabels(xticklabels)
- ax[row, col].tick_params(axis='x', rotation=90, size= 2) # Rotate x-axis labels for readability
- # Set title for each subplot based on the band
- ax[row, col].set_title(titles[i])
- # Set y-axis to log scale for better visualization of power values
- ax[row, col].set_yscale("log")
- # Remove top and right spines to make the plot look cleaner
- ax[row, col].spines['top'].set_visible(False)
- ax[row, col].spines['right'].set_visible(False)
- # Set y-axis label for the first subplot
- if band == 'delta':
- ax[row, col].set_ylabel('Power [\u03bc$V^2$]')
- # Adjust the layout for better spacing
- plt.tight_layout()
- # Save the figure as a PDF
- plt.savefig(os.path.join(save_path, name + 'BandPower_L2_3.pdf'), bbox_inches='tight')
- # Show the plot
- plt.show()
- # %%
- import numpy as np
- from scipy import signal
- import matplotlib.pyplot as plt
- import os
- import seaborn as sns
- # Helper function to append suffix to filenames
- def save_figure_with_suffix(filename_base, suffix='_L2_3', save_path='.', ext='svg'):
- filename = f"{filename_base}{suffix}.{ext}"
- plt.savefig(os.path.join(save_path, filename), bbox_inches='tight')
- # Coherence calculation for data within layer 2/3
- def calculate_coherence(data_list, sf, freq_range=None):
- num_channels = len(data_list)
- coherence_matrix = np.zeros((num_channels, num_channels))
- # Loop through all pairs of channels
- for i in range(num_channels):
- for j in range(i + 1, num_channels):
- f, Cxy = signal.coherence(data_list[i], data_list[j], fs=sf, nperseg=1024)
- if freq_range:
- # Find indices within the desired frequency range
- freq_indices = np.where((f >= freq_range[0]) & (f <= freq_range[1]))
- coherence_matrix[i, j] = np.mean(Cxy[freq_indices])
- else:
- coherence_matrix[i, j] = np.mean(Cxy)
- return coherence_matrix
- def plot_coherence_matrix(matrix, title='', vmin=None, vmax=None):
- plt.figure(figsize=(7, 5))
- ax = sns.heatmap(matrix, annot=False, cmap='hot', vmin=vmin, vmax=vmax, cbar_kws={'label': 'Coherence'})
- ax.set_title(title)
- ax.set_xlabel('Channel')
- ax.set_ylabel('Channel')
- # Define start and end times for base and stim32 separately
- intervals = [(38, 98)]
- fs = target_sf
- # Extract data segments for baseline and stimulation
- clean_data_Vehicle = [data[start * fs:end * fs] for data in clean_data_vehicle for start, end in intervals]
- clean_data_NE = [data[start * fs:end * fs] for data in clean_data_GFC for start, end in intervals]
- clean_data_NE_YHM = [data[start * fs:end * fs] for data in clean_data_GFC_DOB for start, end in intervals]
- # Calculate coherence matrices
- coherence_matrix_Vehicle = calculate_coherence(clean_data_Vehicle, sf=target_sf, freq_range=(30, 80))
- coherence_matrix_NE = calculate_coherence(clean_data_NE, sf=target_sf, freq_range=(30, 80))
- coherence_matrix_NE_YHM = calculate_coherence(clean_data_NE_YHM, sf=target_sf, freq_range=(30, 80))
- # Find global min and max values for the color scale
- min_coherence = min(coherence_matrix_Vehicle.min(), coherence_matrix_NE.min(),
- coherence_matrix_NE_YHM.min(), )
- max_coherence = max(coherence_matrix_Vehicle.max(), coherence_matrix_NE.max(),
- coherence_matrix_NE_YHM.max(), )
- # Plot and save coherence matrices
- for matrix, label in zip([coherence_matrix_Vehicle, coherence_matrix_NE, coherence_matrix_NE_YHM],
- ["Vehicle", "NE", "NE_YHM"]):
- plot_coherence_matrix(matrix, title=f"Coherence Matrix - {label}", vmin=min_coherence, vmax=max_coherence)
- save_figure_with_suffix(f'coherence_matrix_{label}', save_path=save_path_coherence)
- plt.show()
- # Calculate and print average coherence for quantification
- coherence_avgs = {
- "Vehicle": np.mean(coherence_matrix_Vehicle[coherence_matrix_Vehicle > 0]),
- "NE": np.mean(coherence_matrix_NE[coherence_matrix_NE > 0]),
- "NE_YHM": np.mean(coherence_matrix_NE_YHM[coherence_matrix_NE_YHM > 0]),
- }
- for label, avg in coherence_avgs.items():
- print(f"Average Coherence - {label}:", avg)
- # %%
- data_vehicle= []
- data_GFC= []
- data_GFC_DOB= []
- start = 38
- end = 98
- fs = target_sf
- for i in range(len(clean_data_vehicle)):
- data_vehicle.append(downsample_data(clean_data_vehicle[i][start * fs:end * fs], sf, target_sf)[0]) # [::10] is used for downsampling and faster process
- data_GFC.append(downsample_data(clean_data_GFC[i][start * fs:end * fs], sf, target_sf)[0])
- data_GFC_DOB.append(downsample_data(clean_data_GFC_DOB[i][start * fs:end * fs], sf, target_sf)[0])
- # %%
- sf = target_sf
- p_obj = Pac(idpac=(4, 0, 0), f_pha=(1, 30, 1, .2), f_amp=(1, 200, 1, 2)) # 1-30 DTAB
- pha = p_obj.filter(sf, np.array(data_vehicle), ftype='phase', n_jobs=-1)
- amp = p_obj.filter(sf, np.array(data_vehicle), ftype='amplitude', n_jobs=-1)
- pac_vehicle = p_obj.fit(pha, amp)
- pha = p_obj.filter(sf, np.array(data_GFC), ftype='phase', n_jobs=-1)
- amp = p_obj.filter(sf, np.array(data_GFC), ftype='amplitude', n_jobs=-1)
- pac_GFC = p_obj.fit(pha, amp)
- pha = p_obj.filter(sf, np.array(data_GFC_DOB), ftype='phase', n_jobs=-1)
- amp = p_obj.filter(sf, np.array(data_GFC_DOB), ftype='amplitude', n_jobs=-1)
- pac_GFC_DOB = p_obj.fit(pha, amp)
- # ################# saving ndarrays
- import numpy as np
- from pathlib import Path
- path = Path(save_path).expanduser()
- path.mkdir(parents=True, exist_ok=True)
- np.save(path/(name + layer + '_pac_vehicle'), pac_vehicle)
- np.save(path/(name + layer + '_pac_NE'), pac_GFC)
- np.save(path/(name + layer + '_pac_NE_YHM'), pac_GFC_DOB)
- # %%
- import matplotlib.pyplot as plt
- import numpy as np
- condition= ['Vehicle', 'NE',
- 'NE_YHM',
- ]
- pac = [pac_vehicle, pac_GFC,
- pac_GFC_DOB,
- ]
- # Set up the figure size for the plot
- plt.figure(figsize=(12, 9))
- # Initialize variable to store the global vmax
- global_vmax = -np.inf # Start with a very low value
- # First loop to determine global vmax
- for k in pac:
- pacc = k.mean(axis=1) # Compute the mean across time points/trials
- global_vmax = max(global_vmax, pacc.max()) # Update global vmax with the max value across all conditions
- # Loop through each condition and its corresponding PAC
- for i, (k, cond) in enumerate(zip(pac, condition)):
- pacc = k.mean(axis=-1)
- # Check the shape of pacc after mean calculation
- print(f"Shape of pacc for {cond}: {pacc.shape}")
- # Set up the color scaling using the global vmax
- kw = dict(vmax=global_vmax, vmin=0, cmap='jet')
- # Plotting the comodulogram
- plt.subplot(2, 2, i + 1)
- p_obj.comodulogram(pacc, title=f"PAC {cond}", **kw)
- # Adding vertical lines at specific points for reference
- plt.axvline(x=4, color='k', linestyle='--', linewidth=2.0)
- plt.axvline(x=10, color='k', linestyle='--', linewidth=2.0)
- # Adjust layout for better spacing
- plt.tight_layout()
- ####### Saving
- import os.path
- fileName = name2
- plt.savefig(os.path.join(save_path , fileName + layer + '_PAC.svg'),bbox_inches = 'tight')
- plt.savefig(os.path.join(save_path , fileName + layer + '_PAC.pdf'), bbox_inches = 'tight')
- # %%
- # define the preferred phase object
- sf = target_sf
- # pp_obj = PreferredPhase(f_pha=[1, 4], f_amp=(30, 100, 2, 2)) # alpha (1-4 Hz)
- pp_obj = PreferredPhase(f_pha=[4,8], f_amp=(30, 200, 2, 2)) # THETA (4-8 Hz)
- n_jobs = -1
- # compute the preferred phase (reuse the amplitude computed above)
- pha_vehicle = pp_obj.filter(sf, np.array(data_vehicle), ftype='phase', n_jobs=n_jobs)
- amp_vehicle = pp_obj.filter(sf, np.array(data_vehicle), ftype='amplitude', n_jobs=n_jobs)
- pha_GFC = pp_obj.filter(sf, np.array(data_GFC), ftype='phase', n_jobs=n_jobs)
- amp_GFC = pp_obj.filter(sf, np.array(data_GFC), ftype='amplitude', n_jobs=n_jobs)
- pha_GFC_DOB = pp_obj.filter(sf, np.array(data_GFC_DOB), ftype='phase', n_jobs=n_jobs)
- amp_GFC_DOB = pp_obj.filter(sf, np.array(data_GFC_DOB), ftype='amplitude', n_jobs=n_jobs)
- ampbin_vehicle, _, vecbin = pp_obj.fit(pha_vehicle, amp_vehicle, n_bins=72)
- ampbin_GFC, _, vecbin = pp_obj.fit(pha_GFC, amp_GFC, n_bins=72)
- ampbin_GFC_DOB, _, vecbin = pp_obj.fit(pha_GFC_DOB, amp_GFC_DOB, n_bins=72)
- # mean binned amplitude across trials
- ampbin_vehicle = np.squeeze(ampbin_vehicle).mean(-1).T
- ampbin_GFC = np.squeeze(ampbin_GFC).mean(-1).T
- ampbin_GFC_DOB = np.squeeze(ampbin_GFC_DOB).mean(-1).T
- ################# saving ndarrays theta
- import numpy as np
- from pathlib import Path
- path.mkdir(parents=True, exist_ok=True)
- np.save(path/(name2 + layer + '_ampbin_vehicle'), ampbin_vehicle)
- np.save(path/(name2 + layer + '_ampbin_NE'), ampbin_GFC)
- np.save(path/(name2 + layer + '_ampbin_NE_YHM'), ampbin_GFC_DOB)
- # %%
- import numpy as np
- import matplotlib.pyplot as plt
- # List of all the amplitude bins (all frequency conditions)
- ampbins = [ampbin_vehicle, ampbin_GFC, ampbin_GFC_DOB, ]
- # Find the global minimum and maximum values across all amplitude bins
- global_min = np.min([np.min(amp) for amp in ampbins])
- global_max = np.max([np.max(amp) for amp in ampbins])
- # Now set vmin and vmax to these global values
- vmin_value = global_min
- vmax_value = global_max
- # Define the plotting dictionary with consistent color range
- kw_plt = dict(cmap='RdBu_r', interp=0.1, cblabel='Bins Amplitude',
- colorbar=True, y=1.3, fz_title=20, plotas='contour',
- vmin=vmin_value, vmax=vmax_value)
- # Create the figure
- plt.figure(figsize=(27, 23))
- # Plot the data for each frequency condition using the updated kw_plt dictionary
- pp_obj.polar(ampbin_vehicle, vecbin, pp_obj.yvec, subplot=331, title='vehicle', **kw_plt)
- pp_obj.polar(ampbin_GFC, vecbin, pp_obj.yvec, subplot=332, title='NE', **kw_plt)
- pp_obj.polar(ampbin_GFC_DOB, vecbin, pp_obj.yvec, subplot=333, title='NE_YHM', **kw_plt)
- # Adjust layout to prevent overlapping elements
- plt.tight_layout()
- ########### Saving
- import os.path
- fileName = name2
- plt.savefig(os.path.join(save_path , fileName + layer + '_PreferredPhase_ThetaGamma.pdf'), bbox_inches = 'tight')
- # %%
- # %% [markdown]
- # # Layer 5
- # %%
- from scipy import signal
- from scipy.signal import butter, lfilter
- sf = 20000 # sampling frequency
- target_sf = 1000 # target frequency for downsampling
- f0 = 60.0 # Frequency to be removed from signal (Hz)
- w0 = f0/(target_sf/2) # Normalized Frequency
- Q = 30 # Quality factor
- # Design notch filter
- b, a = signal.iirnotch(w0, Q)
- start = 0
- end = int(len(ChannelData_vehicle[0])/sf)
- ################################################### baseline
- clean_data_vehicle = []
- for chanlNo in (L_5):
- data = ChannelData_vehicle[chanlNo].ravel()
- spike_data = filter_data(data[start*sf:end*sf], low=1, high= 500, sf=sf) # Band pass filter the data for LFP analysis
- # Now we use the above function to downsample the signal.
- lfp_data, sf_lfp = downsample_data(spike_data, sf=sf, target_sf=target_sf)
- clean_data_vehicle.append(signal.lfilter(b, a, lfp_data))
- ################################################### kainate
- end = int(len(ChannelData_GFC[0])/sf)
- clean_data_GFC = []
- for chanlNo in (L_5):
- data = ChannelData_GFC[chanlNo].ravel()
- spike_data = filter_data(data[start*sf:end*sf], low=1, high= 500, sf=sf) # Band pass filter the data for LFP analysis
- # Now we use the above function to downsample the signal.
- lfp_data, sf_lfp = downsample_data(spike_data, sf=sf, target_sf=target_sf)
- clean_data_GFC.append(signal.lfilter(b, a, lfp_data))
- # ################################################### kainate + bicuculine
- end = int(len(ChannelData_GFC_DOB[0])/sf)
- clean_data_GFC_DOB = []
- for chanlNo in (L_5):
- data = ChannelData_GFC_DOB[chanlNo].ravel()
- spike_data = filter_data(data[start*sf:end*sf], low=1, high= 500, sf=sf) # Band pass filter the data for LFP analysis
- # Now we use the above function to downsample the signal.
- lfp_data, sf_lfp = downsample_data(spike_data, sf=sf, target_sf=target_sf)
- clean_data_GFC_DOB.append(signal.lfilter(b, a, lfp_data))
- # %%
- import warnings
- warnings.filterwarnings('ignore')
- import numpy as np
- import pandas as pd
- import os
- layer = 'L_5'
- # Frequency bands and data groupings
- freq_bands = {
- 'delta': [1, 4],
- 'theta': [4, 10],
- 'beta': [10, 25],
- 'lowgamma': [30, 60],
- 'highgamma': [60, 100],
- 'gamma': [30, 100],
- 'ripple': [100, 200]
- }
- data_groups = {
- 'vehicle': clean_data_vehicle,
- 'NE ': clean_data_GFC,
- 'NE_YHM': clean_data_GFC_DOB,
- }
- # Time segments
- time_segments = {
- '_immediate': (38, 43),
- '_delayed': (38, 98),
- }
- # Function to compute band power for each time segment
- def compute_band_power(data, fs, freq_band, start, end):
- return bandpower(data[start * fs:end * fs], fs, freq_band, 'welch')
- # Initialize results dictionary for band power per time segment
- band_results = {band: {group: {segment: [] for segment in time_segments} for group in data_groups} for band in freq_bands}
- # Main loop for bandpower calculation with time segments
- for i in range(len(L_5)):
- for band_name, band_range in freq_bands.items():
- for group_name, group_data in data_groups.items():
- for segment_name, (start, end) in time_segments.items():
- # Compute band power for the specific time segment
- bp = compute_band_power(group_data[i], target_sf, band_range, start, end)
- band_results[band_name][group_name][segment_name].append(bp)
- # Prepare data for export, now considering time segments
- export_data = {}
- for band in freq_bands:
- concatenated_data = []
- for group in data_groups:
- for segment in time_segments:
- segment_data = band_results[band][group][segment]
- # Ensure the segment data is not empty and contains valid arrays
- if len(segment_data) > 0:
- # Convert any zero-dimensional arrays into one-dimensional ones
- segment_data = [np.atleast_1d(data) for data in segment_data]
- concatenated_data.append(np.concatenate(segment_data))
- # Only concatenate if we have non-empty arrays
- if len(concatenated_data) > 0:
- export_data[band] = np.concatenate(concatenated_data)
- else:
- export_data[band] = np.array([]) # Handle case where there's no data
- # Create conditions list, including the data groups and time segments
- conditions = np.concatenate([
- len(band_results['delta'][group][segment]) * [f"{group}_{segment}"]
- for group in data_groups for segment in time_segments if len(band_results['delta'][group][segment]) > 0
- ])
- # Calculate average band power for each condition
- average_cols = {}
- for band in freq_bands:
- if len(export_data[band]) > 0:
- # Split data for the band into separate arrays by condition (group + segment)
- data_per_condition = np.array_split(export_data[band], len(data_groups) * len(time_segments))
- avg_per_condition = [np.mean(segment) for segment in data_per_condition]
- # Create a column repeating the average value for alignment with the DataFrame length
- average_col = np.concatenate([[avg] * len(segment) for avg, segment in zip(avg_per_condition, data_per_condition)])
- average_cols[f'Average {band}'] = average_col
- else:
- average_cols[f'Average {band}'] = np.array([]) # Handle empty case
- # Metadata (modify or add actual values as needed)
- layer_col = [layer] * len(conditions) # Replace with actual data if available
- id_col = [name] * len(conditions) # Replace with actual data if available
- strain_col = [strain] * len(conditions) # Replace with actual data if available
- # Create DataFrame for export, now including time segments
- df_export = pd.DataFrame({
- **export_data,
- 'condition': conditions,
- **average_cols,
- 'Layer': layer_col,
- 'ID': id_col,
- 'strain': strain_col
- })
- # Export to Excel
- file_name = f"{name}_BandPower_{layer}.xlsx"
- save_path_full = os.path.join(save_path, file_name)
- with pd.ExcelWriter(save_path_full, engine='xlsxwriter') as writer:
- df_export.to_excel(writer, sheet_name='BandPower', index=False)
- print("Data exported successfully.")
- # %%
- import seaborn as sns
- import matplotlib.pyplot as plt
- # Load the exported data from the previous code
- df = pd.read_excel(os.path.join(save_path, file_name))
- # Frequency bands and titles
- bands = ['delta', 'theta', 'beta', 'lowgamma', 'highgamma', 'ripple']
- titles = [
- r'[$\delta$] Power', r'[$\theta$] Power', r'[$\beta$] Power',
- r'low [$\gamma$] Power', r'high [$\gamma$] Power', r'Cortical Ripple Power'
- ]
- # Updated condition labels based on time segments and groups
- xticklabels = [
- 'Vehicle_immediate', 'Vehicle_delayed',
- 'NE_immediate', 'NE_delayed',
- 'NE_YHM_immediate', 'NE_YHM_delayed'
- ]
- # Create subplots for each frequency band
- fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(15, 5))
- # Loop over each frequency band and create a plot
- for i, band in enumerate(bands):
- row, col = divmod(i, 3) # Determine row and column for each subplot
- # Create boxplot for each condition and frequency band
- sns.boxplot(data=df, x='condition', y=band, ax=ax[row, col], showfliers=False)
- # Create swarmplot (or stripplot) to show individual points on top of the boxplot
- sns.stripplot(data=df, x='condition', y=band, color='k', ax=ax[row, col], dodge=True)
- # Set the x-axis labels to the updated condition names
- ax[row, col].set_xticklabels(xticklabels)
- ax[row, col].tick_params(axis='x', rotation=90, size= 2) # Rotate x-axis labels for readability
- # Set title for each subplot based on the band
- ax[row, col].set_title(titles[i])
- # Set y-axis to log scale for better visualization of power values
- ax[row, col].set_yscale("log")
- # Remove top and right spines to make the plot look cleaner
- ax[row, col].spines['top'].set_visible(False)
- ax[row, col].spines['right'].set_visible(False)
- # Set y-axis label for the first subplot
- if band == 'delta':
- ax[row, col].set_ylabel('Power [\u03bc$V^2$]')
- # Adjust the layout for better spacing
- plt.tight_layout()
- # Save the figure as a PDF
- plt.savefig(os.path.join(save_path, name + 'BandPower_L5.pdf'), bbox_inches='tight')
- # Show the plot
- plt.show()
- # %%
- import numpy as np
- from scipy import signal
- import matplotlib.pyplot as plt
- import os
- import seaborn as sns
- # Helper function to append suffix to filenames
- def save_figure_with_suffix(filename_base, suffix='_L5', save_path='.', ext='svg'):
- filename = f"{filename_base}{suffix}.{ext}"
- plt.savefig(os.path.join(save_path, filename), bbox_inches='tight')
- # Coherence calculation for data within layer 2/3
- def calculate_coherence(data_list, sf, freq_range=None):
- num_channels = len(data_list)
- coherence_matrix = np.zeros((num_channels, num_channels))
- # Loop through all pairs of channels
- for i in range(num_channels):
- for j in range(i + 1, num_channels):
- f, Cxy = signal.coherence(data_list[i], data_list[j], fs=sf, nperseg=1024)
- if freq_range:
- # Find indices within the desired frequency range
- freq_indices = np.where((f >= freq_range[0]) & (f <= freq_range[1]))
- coherence_matrix[i, j] = np.mean(Cxy[freq_indices])
- else:
- coherence_matrix[i, j] = np.mean(Cxy)
- return coherence_matrix
- def plot_coherence_matrix(matrix, title='', vmin=None, vmax=None):
- plt.figure(figsize=(7, 5))
- ax = sns.heatmap(matrix, annot=False, cmap='hot', vmin=vmin, vmax=vmax, cbar_kws={'label': 'Coherence'})
- ax.set_title(title)
- ax.set_xlabel('Channel')
- ax.set_ylabel('Channel')
- # Define start and end times for base and stim32 separately
- intervals = [(0, -1)]
- fs = target_sf
- # Extract data segments for baseline and stimulation
- clean_data_Vehicle = [data[start * fs:end * fs] for data in clean_data_vehicle for start, end in intervals]
- clean_data_NE = [data[start * fs:end * fs] for data in clean_data_GFC for start, end in intervals]
- clean_data_NE_YHM = [data[start * fs:end * fs] for data in clean_data_GFC_DOB for start, end in intervals]
- # Calculate coherence matrices
- coherence_matrix_Vehicle = calculate_coherence(clean_data_Vehicle, sf=target_sf, freq_range=(30, 80))
- coherence_matrix_NE = calculate_coherence(clean_data_NE, sf=target_sf, freq_range=(30, 80))
- coherence_matrix_NE_YHM = calculate_coherence(clean_data_NE_YHM, sf=target_sf, freq_range=(30, 80))
- # Find global min and max values for the color scale
- min_coherence = min(coherence_matrix_Vehicle.min(), coherence_matrix_NE.min(),
- coherence_matrix_NE_YHM.min(), )
- max_coherence = max(coherence_matrix_Vehicle.max(), coherence_matrix_NE.max(),
- coherence_matrix_NE_YHM.max())
- # Plot and save coherence matrices
- for matrix, label in zip([coherence_matrix_Vehicle, coherence_matrix_NE, coherence_matrix_NE_YHM, ],
- ["Vehicle", "NE", "NE_YHM", ]):
- plot_coherence_matrix(matrix, title=f"Coherence Matrix - {label}", vmin=min_coherence, vmax=max_coherence)
- save_figure_with_suffix(f'coherence_matrix_{label}', save_path=save_path_coherence)
- plt.show()
- # Calculate and print average coherence for quantification
- coherence_avgs = {
- "Vehicle": np.mean(coherence_matrix_Vehicle[coherence_matrix_Vehicle > 0]),
- "NE": np.mean(coherence_matrix_NE[coherence_matrix_NE > 0]),
- "NE_YHM": np.mean(coherence_matrix_NE_YHM[coherence_matrix_NE_YHM > 0]),
- }
- for label, avg in coherence_avgs.items():
- print(f"Average Coherence - {label}:", avg)
- # %% [markdown]
- # # Preferred Phase
- # %%
- # The preferred phase is given by the phase bin at which the amplitude is maximum. Said differently,
- # the gamma amplitude is binned according to the alpha phase. The preferred-phase is computed as the
- # maximum of this histogram. To identify this preferred-phase, we propose here a polar plotting method to
- # visualize retrieve the previous result using the histogram method
- # %%
- data_vehicle= []
- data_GFC= []
- data_GFC_DOB= []
- start = 38
- end = 98
- fs = target_sf
- for i in range(len(clean_data_vehicle)):
- data_vehicle.append(downsample_data(clean_data_vehicle[i][start * fs:end * fs], sf, target_sf)[0]) # [::10] is used for downsampling and faster process
- data_GFC.append(downsample_data(clean_data_GFC[i][start * fs:end * fs], sf, target_sf)[0])
- data_GFC_DOB.append(downsample_data(clean_data_GFC_DOB[i][start * fs:end * fs], sf, target_sf)[0])
- # %% [markdown]
- # # PAC (Phase-amplitude coupling) Theta-Gamma
- # %% [markdown]
- # # stop ! Change the NAME of the drug
- # %%
- sf = target_sf
- p_obj = Pac(idpac=(4, 0, 0), f_pha=(1, 30, 1, .2), f_amp=(1, 200, 1, 2)) # 1-30 DTAB
- pha = p_obj.filter(sf, np.array(data_vehicle), ftype='phase', n_jobs=-1)
- amp = p_obj.filter(sf, np.array(data_vehicle), ftype='amplitude', n_jobs=-1)
- pac_vehicle = p_obj.fit(pha, amp)
- pha = p_obj.filter(sf, np.array(data_GFC), ftype='phase', n_jobs=-1)
- amp = p_obj.filter(sf, np.array(data_GFC), ftype='amplitude', n_jobs=-1)
- pac_GFC = p_obj.fit(pha, amp)
- pha = p_obj.filter(sf, np.array(data_GFC_DOB), ftype='phase', n_jobs=-1)
- amp = p_obj.filter(sf, np.array(data_GFC_DOB), ftype='amplitude', n_jobs=-1)
- pac_GFC_DOB = p_obj.fit(pha, amp)
- # ################# saving ndarrays
- import numpy as np
- from pathlib import Path
- path = Path(save_path).expanduser()
- path.mkdir(parents=True, exist_ok=True)
- np.save(path/(name + layer + '_pac_vehicle'), pac_vehicle)
- np.save(path/(name + layer + '_pac_NE'), pac_GFC)
- np.save(path/(name + layer + '_pac_NE_YHM'), pac_GFC_DOB)
- # %%
- import matplotlib.pyplot as plt
- import numpy as np
- condition= ['Vehicle', 'NE',
- 'NE_YHM',
- ]
- pac = [pac_vehicle, pac_GFC,
- pac_GFC_DOB,
- ]
- # Set up the figure size for the plot
- plt.figure(figsize=(12, 9))
- # Initialize variable to store the global vmax
- global_vmax = -np.inf # Start with a very low value
- # First loop to determine global vmax
- for k in pac:
- pacc = k.mean(axis=1) # Compute the mean across time points/trials
- global_vmax = max(global_vmax, pacc.max()) # Update global vmax with the max value across all conditions
- # Loop through each condition and its corresponding PAC
- for i, (k, cond) in enumerate(zip(pac, condition)):
- pacc = k.mean(axis=-1)
- # Check the shape of pacc after mean calculation
- print(f"Shape of pacc for {cond}: {pacc.shape}")
- # Set up the color scaling using the global vmax
- kw = dict(vmax=global_vmax, vmin=0, cmap='jet')
- # Plotting the comodulogram
- plt.subplot(2, 2, i + 1)
- p_obj.comodulogram(pacc, title=f"PAC {cond}", **kw)
- # Adding vertical lines at specific points for reference
- plt.axvline(x=4, color='k', linestyle='--', linewidth=2.0)
- plt.axvline(x=10, color='k', linestyle='--', linewidth=2.0)
- # Adjust layout for better spacing
- plt.tight_layout()
- ####### Saving
- import os.path
- fileName = name2
- plt.savefig(os.path.join(save_path , fileName + layer + '_PAC.svg'),bbox_inches = 'tight')
- plt.savefig(os.path.join(save_path , fileName + layer + '_PAC.pdf'), bbox_inches = 'tight')
- # %%
- # %% [markdown]
- # # Preferred Phase Theta-Gamma
- # %% [markdown]
- # # stop ! Change the NAME of the drug
- # %%
- # define the preferred phase object
- sf = target_sf
- # pp_obj = PreferredPhase(f_pha=[1, 4], f_amp=(30, 100, 2, 2)) # alpha (1-4 Hz)
- pp_obj = PreferredPhase(f_pha=[4,8], f_amp=(30, 200, 2, 2)) # THETA (4-8 Hz)
- n_jobs = -1
- # compute the preferred phase (reuse the amplitude computed above)
- pha_vehicle = pp_obj.filter(sf, np.array(data_vehicle), ftype='phase', n_jobs=n_jobs)
- amp_vehicle = pp_obj.filter(sf, np.array(data_vehicle), ftype='amplitude', n_jobs=n_jobs)
- pha_GFC = pp_obj.filter(sf, np.array(data_GFC), ftype='phase', n_jobs=n_jobs)
- amp_GFC = pp_obj.filter(sf, np.array(data_GFC), ftype='amplitude', n_jobs=n_jobs)
- pha_GFC_DOB = pp_obj.filter(sf, np.array(data_GFC_DOB), ftype='phase', n_jobs=n_jobs)
- amp_GFC_DOB = pp_obj.filter(sf, np.array(data_GFC_DOB), ftype='amplitude', n_jobs=n_jobs)
- ampbin_vehicle, _, vecbin = pp_obj.fit(pha_vehicle, amp_vehicle, n_bins=72)
- ampbin_GFC, _, vecbin = pp_obj.fit(pha_GFC, amp_GFC, n_bins=72)
- ampbin_GFC_DOB, _, vecbin = pp_obj.fit(pha_GFC_DOB, amp_GFC_DOB, n_bins=72)
- # mean binned amplitude across trials
- ampbin_vehicle = np.squeeze(ampbin_vehicle).mean(-1).T
- ampbin_GFC = np.squeeze(ampbin_GFC).mean(-1).T
- ampbin_GFC_DOB = np.squeeze(ampbin_GFC_DOB).mean(-1).T
- ################# saving ndarrays theta
- import numpy as np
- from pathlib import Path
- path.mkdir(parents=True, exist_ok=True)
- np.save(path/(name2 + layer + '_ampbin_vehicle'), ampbin_vehicle)
- np.save(path/(name2 + layer + '_ampbin_NE'), ampbin_GFC)
- np.save(path/(name2 + layer + '_ampbin_NE_YHM'), ampbin_GFC_DOB)
- # %%
- import numpy as np
- import matplotlib.pyplot as plt
- # List of all the amplitude bins (all frequency conditions)
- ampbins = [ampbin_vehicle, ampbin_GFC, ampbin_GFC_DOB, ]
- # Find the global minimum and maximum values across all amplitude bins
- global_min = np.min([np.min(amp) for amp in ampbins])
- global_max = np.max([np.max(amp) for amp in ampbins])
- # Now set vmin and vmax to these global values
- vmin_value = global_min
- vmax_value = global_max
- # Define the plotting dictionary with consistent color range
- kw_plt = dict(cmap='RdBu_r', interp=0.1, cblabel='Bins Amplitude',
- colorbar=True, y=1.3, fz_title=20, plotas='contour',
- vmin=vmin_value, vmax=vmax_value)
- # Create the figure
- plt.figure(figsize=(27, 23))
- # Plot the data for each frequency condition using the updated kw_plt dictionary
- pp_obj.polar(ampbin_vehicle, vecbin, pp_obj.yvec, subplot=331, title='vehicle', **kw_plt)
- pp_obj.polar(ampbin_GFC, vecbin, pp_obj.yvec, subplot=332, title='NE', **kw_plt)
- pp_obj.polar(ampbin_GFC_DOB, vecbin, pp_obj.yvec, subplot=333, title='NE_YHM', **kw_plt)
- # Adjust layout to prevent overlapping elements
- plt.tight_layout()
- # ########### Saving
- # import os.path
- # fileName = name2
- # plt.savefig(os.path.join(save_path , fileName + layer + '_PreferredPhase_ThetaGamma.pdf'), bbox_inches = 'tight')
- # %%
- # %%
- # %%
- # %%
- # %%
- # %%
- # %%
- # %%
- # %%
- # %%
- # %%
- # %%
MEA LFP _ Open OnDemand_LC_mPFC.ipynb at commit 29e3811, no license · at the source
Overview
- Department of Pharmacology, University of Michigan Medical School, Ann Arbor, Michigan 48109
- Neuroscience Graduate Program, University of Michigan Medical School, Ann Arbor, Michigan 48109
Abstract
Cognitive flexibility—the ability to adapt behavior when contingencies change—is impaired in psychiatric disorders involving prefrontal dysfunction. The medial prefrontal cortex (mPFC) relies on noradrenergic input from the locus ceruleus (LC), yet the molecular mechanisms enabling this neuromodulatory control remain unclear. Here we show that GluN2A-containing NMDA receptors are required for LC–mPFC regulation of network dynamics and reversal learning in male mice. Optogenetic activation of LC→mPFC projections enhanced reversal learning in wild-type (WT) and heterozygous mice but not in global Grin2a knock-outs, whereas LC inhibition impaired performance only in WT animals. In slices, norepinephrine (NE) and LC stimulation induced gamma and high-frequency oscillations in WT mPFC that were blocked by α2-adrenergic antagonism, but these oscillatory responses were undetectable in Grin2a mutants. Grin2a mutants also exhibited increased LC axonal density and elevated NE transporter expression in the prelimbic cortex, consistent with enhanced noradrenergic clearance capacity. Together, these findings identify GluN2A as a key determinant of LC–prefrontal circuit function supporting cognitive flexibility. They suggest that functional deficits in these mutants should be interpreted within the context of compensatory structural hyperinnervation resulting from global GluN2A deficiency, which may reflect a developmental adaptation rather than acute signaling loss. Furthermore, these results promote α2-adrenergic pathways as potential entry points for restoring prefrontal network coordination.
Reproduced under the paper's license (CC BY), from the paper cited above.
Repositories
Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.
NeuroDataa/PatchClampAnalysis
Availability: 1 check, the latest on 27 September 2026: the link is dead
- 27 September 2026: the link is dead
NeuroDataa/Grin2a_LC_mPFC_pMEA_Analysis
29e3811b47f2151df71c5dd3445ada824d04a819, 13 May 2025Availability: 1 check, the latest on 27 September 2026: the link answers
- 27 September 2026: the link answers
2 files
- MEA LFP _ Open OnDemand_LC_mPFC.ipynb, Jupyter, 1,495 lines, 1 match
- MEA LFP _ Open OnDemand_NE_Pharmacology
.ipynb , Jupyter, 1,376 lines, 1 match
NeuroData/PatchClampAnalysis
Availability: 1 check, the latest on 27 September 2026: the link is dead
- 27 September 2026: the link is dead
NeuroDataa/MEA_Analysis
Availability: 1 check, the latest on 27 September 2026: the link answers
- 27 September 2026: the link answers
The paper's code and data availability statement is in the Data section.
Tracing map
Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.
What the map holds:
- 4 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
- 2 scripts, each with its path and the digest of its content;
- 2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
- neither the text of the paper nor the code itself.
Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.
Data
No dataset and no data link were found in the paper.
Data and code availability
All original analysis code generated for this study has been deposited in publicly accessible GitHub repositories:
https://
This study did not generate new or unique reagents.
Reproduced under the paper's license (CC BY), from the paper cited above.
Versions
The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.
Version 1, 27 September 2026: the first record
Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 6 keywords, 12 MeSH terms, 1 funder, 33 references.
Cite
This paper
Hosseini, H., Evans-Martin, S., Bogomilsky, E., & Jones, K. S. (2026). GluN2A Enables Noradrenergic Control of Prefrontal Oscillations and Cognitive Flexibility. eNeuro, 13(8), ENEURO.0050-26.2026. https://
BibTeX
@article{hosseini2026glu
author = {Hosseini, Hassan and Evans-Martin, Sky and Bogomilsky, Emma and Jones, Kevin S},
title = {{GluN2A Enables Noradrenergic Control of Prefrontal Oscillations and Cognitive Flexibility}},
journal = {eNeuro},
year = {2026},
month = aug,
volume = {13},
number = {8},
pages = {ENEURO.0050--26.2026},
publisher = {Society for Neuroscience},
issn = {2373-2822},
doi = {10.1523/
url = {https://
pmid = {42498666},
pmcid = {PMC13456957}
}
RIS
TY - JOUR
AU - Hosseini, Hassan
AU - Evans-Martin, Sky
AU - Bogomilsky, Emma
AU - Jones, Kevin S
TI - GluN2A Enables Noradrenergic Control of Prefrontal Oscillations and Cognitive Flexibility
T2 - eNeuro
J2 - eNeuro
PY - 2026
DA - 2026/
VL - 13
IS - 8
SP - ENEURO.0050
EP - 26.2026
SN - 2373-2822
PB - Society for Neuroscience
DO - 10.1523/
UR - https://
LA - en
ER -
CSL-JSON
{
"id": "10.1523/
"type": "article-journal",
"title": "GluN2A Enables Noradrenergic Control of Prefrontal Oscillations and Cognitive Flexibility",
"container-title": "eNeuro",
"author": [
{
"family": "Hosseini",
"given": "Hassan"
},
{
"family": "Evans-Martin",
"given": "Sky"
},
{
"family": "Bogomilsky",
"given": "Emma"
},
{
"family": "Jones",
"given": "Kevin S"
}
],
"container-title-short":
"volume": "13",
"issue": "8",
"page": "ENEURO.0050-26.2026",
"DOI": "10.1523/
"PMID": "42498666",
"PMCID": "PMC13456957",
"ISSN": "2373-2822",
"publisher": "Society for Neuroscience",
"URL": "https://
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
7
]
]
}
}
The tracing map gets a citation of its own once an author has validated it and it has a DOI.
Similar papers
The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.
- [1] doi:10.1038/s41593-026-02258-4 [code]
- Laminar organization of cellular microcircuits modulating human interictal epileptiform discharges.Journal: Nature neuroscienceIn common: MNE-Python, seaborn, scikit-learn, 4 other tools, 1 reference
- [2] doi:10.3389/fpsyg.2026.1774068 [code]
- Analysis of cognitive mechanisms in phoneme perception and pronunciation errors among Korean language learners.Journal: Frontiers in psychologyIn common: MNE-Python, h5py, seaborn, 5 other tools
- [3] doi:10.3389/fncom.2026.1786996 [code]
- Schumann-anchored golden ratio organization of human neural oscillations.Journal: Frontiers in computational neuroscienceIn common: MNE-Python, h5py, seaborn, 5 other tools
- [4] doi:10.1038/s41598-026-50946-9 [code]
- Fixation-related potentials reveal that confusing program code elicits a late frontal positivity.Journal: Scientific reportsIn common: MNE-Python, h5py, seaborn, 5 other tools
- [5] doi:10.1038/s41593-026-02285-1 [code]
- Fixation duration on natural scenes is explained by memory encoding not processing demand.Journal: Nature neuroscienceIn common: MNE-Python, h5py, seaborn, 5 other tools
- [6] doi:10.1162/imag.a.1227 [code]
- Large language models reveal the neural tracking of linguistic context in attended and unattended multi-talker speech.Journal: Imaging neuroscience (Cambridge, Mass.)In common: MNE-Python, h5py, seaborn, 5 other tools
- [7] doi:10.1038/s41531-026-01372-1 [code]
- Varying patterns of association between cortical large-scale networks and subthalamic nucleus activity in Parkinson's disease.Journal: NPJ Parkinson's diseaseIn common: MNE-Python, h5py, seaborn, 5 other tools
- [8] doi:10.7554/elife.106543 [code]
- Stimulus dependencies-rather than next-word prediction-can explain pre-onset brain encoding in naturalistic listening designs.Journal: eLifeIn common: MNE-Python, h5py, seaborn, 5 other tools
- [9] doi:10.1117/1.nph.13.2.025001 [code]
- Surface-based image reconstruction optimization for high-density functional near-infrared spectroscopy.Journal: NeurophotonicsIn common: MNE-Python, h5py, seaborn, 5 other tools
- [10] doi:10.1038/s41597-025-05174-7 [code]
- A large-scale MEG and EEG dataset for object recognition in naturalistic scenesJournal: n/aIn common: MNE-Python, h5py, seaborn, 5 other tools
Contribute
The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.
Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.
Claim this paper
Correct its record
Say what each link of this record is, remove the ones that are not the paper's, add the ones that are missing. The correction becomes a new version of the record, in its Versions section.
Validate its tracing map
You validate the map as this page shows it: 4 repositories of the authors' code, each at its verified commit and with its license, 2 scripts, and 2 matches between paragraphs and code (see the Code and Map sections). It then receives a DOI on Zenodo, with you (your ORCID iD) and OSCR as its creators; the code itself is not deposited.
The map's fingerprint: sha256:c090dbdb6194630b…
Add the badge to its README
The badge links the code to this page. Copy one of these into the README of the paper's code: only you decide where it goes, and nothing is changed for you.
Markdown
[, paste the snippet at the top, then “Commit changes…” and, to review it first, “Create a new branch and start a pull request”. You open the pull request; OSCR asks for no permission.
Request its removal
To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).
Discussion, reproductions, activity
Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.
Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.
Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.
