Neural population dynamics of direct electrical stimulation of neocortex.

Code ↔ Paper

7 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 7 matches

[1] § STAR★METHODS › QUANTIFICATION AND STATISTICAL ANALYSIS › Evoked potential analysis › Aligning and averaging the evoked potential ↔ figure2.ipynb, lines 805–922 · score 0.62 · CCF distance, peak aligned, normalized distance, evoked potential, boundaries, amplitude
[2] § STAR★METHODS › METHOD DETAILS › Visual and electrical stimulation behavioral detection assay › Software and hardware implementation: ↔ dev.py, lines 242–287 · score 0.59 · spout positioning, Pyglet, solenoid, lick, monitors, delivery
[3] § STAR★METHODS › METHOD DETAILS › Visual and electrical stimulation behavioral detection assay › Task protocol: ↔ plotter.py, lines 54–179 · score 0.59 · water delivery, daily, days, weight, curves, reward
[4] § STAR★METHODS › QUANTIFICATION AND STATISTICAL ANALYSIS › Anisotropic volumetric analyses › Convex-hull comparison to an isotropic sphere hull ↔ figure3.ipynb, lines 1754–1827 · score 0.54 · convex hull, sphere hull, tessellated, Figure 3
[5] § RESULTS › Evoked potential spatial extent increases sub-linearly with amplitude ↔ figure2.ipynb, lines 78–136 · score 0.54 · radial distance, point source model, meters, evoked potentials, Figure 2
[6] § RESULTS › Anisotropic volumetric profile of the induced electrical field ↔ figure3.ipynb, lines 891–980 · score 0.50 · Euclidean distance, fitted spheres, mass, contact, radius, Figure 3
[7] § STAR★METHODS › METHOD DETAILS › Visual and electrical stimulation behavioral detection assay › Trial structure: ↔ params.py, lines 26–77 · score 0.50 · variable, lapse, alarm, licking, Catch, 1.5 s

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 paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Jupyter notebook · 2,260 lines · 88 KB · GPL-3.0 · 2 matches

# %%
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os,glob
import seaborn as sns
from open_ephys.analysis import Session
from pynwb import NWBHDF5IO, NWBFile
from pynwb.file import Subject
from datetime import datetime
from dateutil.tz import tzlocal
from jlh_ephys.analysis import Analysis
from jlh_ephys.utils import choose_stim_parameter, OE
from jlh_ephys.psth_raster import psth_arr
#from jlh_ephys import spatial_tools as st ## compatability issues maybe need to install from allensdk.core.mouse_connectivity_cache import MouseConnectivityCache in its own env
from jlh_ephys.raw import find_artifact_start, align_data, raw_heatmap
from ccf_3D.tools import herbs_processing as hp
from ccf_3D.tools.metrics import distance
import matplotlib
matplotlib.rcParams['pdf.fonttype'] = 42
matplotlib.rcParams['ps.fonttype'] = 42
from scipy.ndimage import gaussian_filter1d
# %%
jlh31 = Analysis('jlh31', '2023-01-24') # all probes, many params and sites
jlh32 = Analysis('jlh32', '2023-02-18') # all probes, many params and sites
jlh33 = Analysis('jlh33', '2023-02-22') # probe B and C, many params and trains in diff recordings
jlh34 = Analysis('jlh34', '2023-05-15') # all 3 many params
jlh39 = Analysis('jlh39', '2023-08-16') # 3 probes for bi, no A for mono
jlh40 = Analysis('jlh40', '2023-08-21') # probeB and probe C only: monopolar -- with diff pulse width
jlh48 = Analysis('jlh48', '2024-11-05') # bb test, sorted only no bb portion (full bi and mono)
jlh49 = Analysis('jlh49v1', '2024-07-11') # bb test, just no bb portion (full bi only)
jlh53 = Analysis('jlh53concat', '2024-09-03') # regular, all 3 probes, probe A broke sometime.
recordings = [jlh31, jlh32, jlh33, jlh34, jlh39, jlh40, jlh48, jlh49, jlh53]
for r in recordings:
r.get_brain_regs()
r.get_electrode_coords()
nn_contacts_labels = [6,11,3,14,1,16,2,15,5,12,4,13,7,10,8,9]
EP_path = r"C:\Users\jordan\Documents\combined_figs\evoked_potentials"
amp_alphas = {5: 0.3, 25: 0.5, 50: 0.7, 100: 1.0}
bipolar_color = '#82A69D'
cathodal_color = '#161943'
anodal_color = '#D99551'
# %%
distances = {}
distances_norm = {}
for r in recordings:
contact = r.trials.contact_negative[0]
distances[r.mouse] = r.get_dists(contact)
temp_norm_dist = {}
for probe, dists in r.distances.items():
# normalize min dist to 0
min = np.min(dists)
dists_norm = dists - min
# make prior to crossing point negative and after positive
cross_index = np.where(dists_norm == 0)[0][0] # this gives the first index where data is non-positive
dists_norm[cross_index:] = -np.abs(dists_norm[cross_index:])
temp_norm_dist[probe] = dists_norm
dists[cross_index:] = -np.abs(dists[cross_index:])
r.distances[probe] = dists
distances_norm[r.mouse] = temp_norm_dist
# %% [markdown]
# # point source modeling (generating a null hypothesis) Figure 2a-b
# %%
rho = 5.56
currents_uA = [5, 25, 50, 100]
currents_A = [i * 1e-6 for i in currents_uA]
# Radial distances (meters)
r_min = 0.5e-4 # 0.05 mm or 50 microns
r_max = 2e-3 # 1 mm
N_points = 600
r_values = np.linspace(r_min, r_max, N_points)
fig = plt.figure(figsize=(3.54, 3.54))
grays = ['0.8', '0.5', '0.3', '0.1']
half_maxs_dist = [] # duh they're all the same LOL ignore this
full_maxs_dists = []
full_max_threshold = 0.1
for i, I in enumerate(currents_A):
# Compute potential (in Volts)
V = (rho * I) / (4 * np.pi * r_values)
# Distance in mm for plotting
r_mm = r_values * 1e3
hundred_microns = np.abs(r_mm - 0.100).argmin()
max = V[hundred_microns]
half_max = max / 2
half_max_dist = np.abs(V - half_max).argmin()
half_maxs_dist.append(r_mm[half_max_dist])
full_max_dist = np.abs(V - full_max_threshold).argmin()
full_maxs_dists.append(r_mm[full_max_dist])
# Plot: x is potential (V), y is distance (mm)
label_str = f"I = {I*1e6:.0f} µA"
plt.plot(V, r_mm, linewidth=2, color=grays[i], label=label_str)
plt.plot(V, -r_mm, linewidth=2, color=grays[i]) # mirror around 0 for ±r
plt.xlabel('Potential V (V)', fontsize=10)
plt.ylabel('Distance r (mm)', fontsize=10)
plt.title('Point-Source Potential in a Homogeneous Medium', fontsize=10)
# Match the original axis ranges:
# Potential from 0 to 1, distance from -1 mm to 1 mm
plt.xlim(0, 1)
plt.ylim(-1, 1)
plt.legend(fontsize=8)
plt.tight_layout()
save_str = r'C:\Users\jordan\Documents\combined_figs\evoked_potentials\point_source_modeling\potential_vs_distance.png'
#plt.savefig(save_str)
save_str = save_str = r'C:\Users\jordan\Documents\combined_figs\evoked_potentials\point_source_modeling\potential_vs_distance.pdf'
#plt.savefig(save_str)
full_maxs_dists
# %%
# point-source spatial extent vs current using analytic threshold crossing
rho = 5.56 # resistivity (Ω·m)
currents_uA = [5, 25, 50, 100] # µA
currents_A = np.array(currents_uA) * 1e-6
V_thresh = 0.05 # volts; "near zero" threshold (try 0.02–0.1)
r_thresh_m = (rho * currents_A) / (4 * np.pi * V_thresh)
r_thresh_mm = r_thresh_m * 1e3
fig, ax = plt.subplots(figsize=(2, 2.5))
ax.plot(currents_uA, r_thresh_mm, '-', lw=2, color = 'k')
ax.scatter(currents_uA, r_thresh_mm, color = 'k', s=28)
ax.set_xlabel('current (µA)')
ax.set_ylabel(r'$r_{\mathrm{th}}$ (mm)')
ax.set_title(f'spatial extent at V_th = {V_thresh:.3f} V')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.grid(False)
ax.set_yticks([0, 0.5, 1.0])
ax.set_xticks([0, 25, 50, 100])
plt.tight_layout()
save_str = r'C:\Users\jordan\Documents\combined_figs\evoked_potentials\point_source_modeling\potential_vs_distance_spatialextent.pdf'
#plt.savefig(save_str)
# %%
# %% [markdown]
# # algorithm for calculating spatial boundaries of evoked potentials (load data from csv of pre-calculated)
# %%
## note this is will only work if you have raw data which is available upon request
# thus this is mainly demonstrative -- I provide a curated dataframe of the results where all figures were calculated
from scipy.ndimage import gaussian_filter1d
pre = -0.3
post = 5
snapshots = [1] # ms
snapshot_idxs = [int(pre / 1000 * 30_000) + int(s / 1000 * 30_000) for s in snapshots]
snap_alphas = [1]
data_list = []
for r in recordings:
for probe in r.probes:
if r.mouse in ['jlh33', 'jlh40']: # this fixes the issue with recorded vs geometric probes
if probe == 'probeB':
raw_probe = 'probeA'
elif probe == 'probeC':
raw_probe = 'probeB'
else:
raw_probe = probe
channels = len(r.probe_coords[probe]) + 20
for run in r.trials.run.unique():
fig, axs = plt.subplots(1, 2, figsize=(11, 8))
stim_times = r.trials[r.trials.run == run].start_time.values
title = f'{probe}_Run{run} {r.parameters[run]}'
data = r.raw.get_chunk(raw_probe, stim_times, pre=pre + 3, post=post + 3, chs=np.arange(channels))
try:
aligned_data = align_data(data, pre, post, channels, threshold=500, median_subtraction=False)
except:
print(f'Error in recording {r.mouse} probe {probe} amp {amp}')
aligned_data = data
dists = distances[r.mouse][probe]
sub_chs = np.arange(len(dists), len(dists) + 15)
aligned_sub_data = aligned_data - np.median(aligned_data[:,:, sub_chs], axis=2)[:,:, np.newaxis]
# Use the axs object for the current amp
raw_heatmap(aligned_sub_data, pre=pre, post=post, dists=dists, vmin=-500, vmax=500,
save=False, title=f'{r.parameters[run]}', ax=axs[0])
mean_aligned_sub_data = np.mean(aligned_sub_data, axis=0)
for snap_idx, (snapshot_time, idx) in enumerate(zip(snapshots, snapshot_idxs)):
data_snap = aligned_data[:, idx, :]
data_snap = np.abs(data_snap - np.median(data_snap[:, sub_chs], axis=1)[:, np.newaxis])
data_snap[:,191] = data_snap[:,192] # ch 191 seems to be the reference?
# algo for slope detection
smoothed_gauss = gaussian_filter1d(data_snap, sigma=2, axis=1)
mean_data = np.mean(data_snap, axis=0)
mean_gauss = np.mean(smoothed_gauss, axis=0)
sd_gauss = np.std(smoothed_gauss, axis=0)
baseline_mean = np.mean(mean_gauss[-15:])
baseline_std = np.std(mean_gauss[-15:])
max_idx = np.argmax(mean_gauss[30:])+30
max_volt = mean_gauss[max_idx]
max_threshold = 40
min_voltage = 20
half_max = max_volt / 2
if max_volt > max_threshold:
uppers = mean_gauss[max_idx:]
lowers = mean_gauss[:max_idx] # doesn't allow bottom 30 channels to be considered
upper_chs_outside_potential = np.where(uppers < min_voltage)[0]
if len(upper_chs_outside_potential) > 0:
upper_bound = (upper_chs_outside_potential[0] - 2) + max_idx #minus two to get channel safely within range
half_max_upper = np.where(uppers > half_max)[0][-1] + max_idx
else:
upper_bound = -1
half_max_upper = -1
lower_chs_outside_potential = np.where(lowers < min_voltage)[0]
if len(lower_chs_outside_potential) > 0:
lower_bound = lower_chs_outside_potential[-1] + 2 #plus 2 to get channel safely within range
half_max_lower = np.where(lowers < half_max)[0][-1]
else:
lower_bound = -1
half_max_lower = -1
else:
lower_bound, upper_bound = -1,-1
axs[1].plot(mean_data, np.arange(0,channels), color='k', label='Raw Mean')
axs[1].plot(mean_gauss, np.arange(0,channels), color='blue', label='Gaussian Mean')
axs[1].fill_betweenx(np.arange(0,channels), mean_gauss - sd_gauss, mean_gauss + sd_gauss, color='blue', alpha=0.3, label='SD')
axs[1].set_ylabel('Channel')
axs[1].set_xlabel('Voltage')
axs[1].set_title('snapshots')
axs[1].set_yticks(np.arange(0,channels,10))
axs[1].legend()
if max_volt > max_threshold:
axs[1].axhline(lower_bound, color='cyan', linestyle='--', label='Slope Lower Bound')
axs[1].axhline(upper_bound, color='cyan', linestyle='--', label='Slope Upper Bound')
axs[1].axhline(half_max_lower, color = 'orange', linestyle = '--', label = 'Half Max Lower')
axs[1].axhline(half_max_upper, color = 'orange', linestyle = '--', label = 'Half Max Upper')
axs[1].axhline(max_idx, color='black', linestyle='--', label='Max')
axs[0].axhline(lower_bound, color='cyan', linestyle='--')
axs[0].axhline(upper_bound, color='cyan', linestyle='--')
axs[0].axhline(half_max_lower, color = 'orange', linestyle = '--')
axs[0].axhline(half_max_upper, color = 'orange', linestyle = '--',)
axs[0].axhline(max_idx, color='black', linestyle='--')
fig.suptitle(title, size=20)
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
base_path = os.path.join(EP_path, "ap_potential_w_voltageline")
sub_folder = os.path.join(base_path, "all_runs_med_sub_spatialdetection", 'abs_volt_threshold_from_max40_min20_halfmax_sigma3', '1ms', r.mouse)
if not os.path.exists(sub_folder):
os.makedirs(sub_folder)
#plt.savefig(os.path.join(sub_folder, f'{r.mouse}_{probe}_run{run}.pdf'))
#plt.savefig(os.path.join(sub_folder, f'{r.mouse}_{probe}_run{run}.png'))
plt.close()
# Append all relevant details to the list
data_list.append({
'recording': r.mouse,
'probe': probe,
'run': run,
'stim_parameters': r.parameters[run],
'max_idx': max_idx,
'max_volt': max_volt,
'half_max': half_max,
'half_max_upper_idx': half_max_upper,
'half_max_lower_idx': half_max_lower,
'lower_bound_idx': lower_bound,
'upper_bound_idx': upper_bound,
'mean_gauss': mean_gauss.tolist(),
'std_gauss': sd_gauss.tolist(),
'mean_aligned_sub_data': mean_aligned_sub_data
})
# Convert the list to a dataframe
df_results_w_array = pd.DataFrame(data_list)
# %%
# load the results from the EP boundaries algorithm with some intermediaries
import pickle
df_results = pd.read_pickle(os.path.join('intermediates', 'EP_boundaries_w_array_distances_curated.pkl'))
recording_names = ['jlh31', 'jlh32', 'jlh33', 'jlh34', 'jlh39', 'jlh40', 'jlh48', 'jlh49v1', 'jlh53concat']
recordings = [jlh31, jlh32, jlh33, jlh34, jlh39, jlh40, jlh48, jlh49, jlh53]
# Create the recording map
recording_map = {name: recording for name, recording in zip(recording_names, recordings)}
with open(r'intermediates\mean_aligned_sub_data.pkl', 'rb') as f:
array_list = pickle.load(f)
df_results['mean_aligned_sub_data'] = array_list
# %% [markdown]
# # average EP heatmaps (Figure 2 e)
# %%
# note, requires raw data to run this code, which is available upon request
from scipy.ndimage import gaussian_filter
# average fields by polarity and amplitude
polarities = ['bipolar', 'monopolar']
fig, axs = plt.subplots(4, 4, figsize=(15, 12)) # 3 rows (probes), 4 columns (amplitudes)
fig.suptitle(f'1ms aligned spatial response by polarities', fontsize=14)
probes = ['probeA', 'probeB', 'probeC']
amps = [-5, -25, -50, -100, 5, 25, 50, 100]
pad = 400
mean_array_bycondition = {}
pre = -0.3
post = 5
df_results_w_array = df_results.copy()
for polarity in polarities:
for i, amp in enumerate(amps):
if amp < 0 and polarity == 'bipolar':
row = 0
elif amp > 0 and polarity == 'bipolar':
row = 1
elif amp < 0 and polarity == 'monopolar':
row = 2
elif amp > 0 and polarity == 'monopolar':
row = 3
if np.abs(amp) == 5:
col = 0
elif np.abs(amp) == 25:
col = 1
elif np.abs(amp) == 50:
col = 2
elif np.abs(amp) == 100:
col = 3
ax = axs[row, col]
# Filter data for the current probe and amplitude
df = df_results_w_array[
(df_results_w_array['max_volt'] > 40) &
(df_results_w_array['amplitude'] == amp) &
(df_results_w_array['polarity'] == polarity) &
(df_results_w_array['half_max_lower_idx'] != -1) &
(df_results_w_array['half_max_upper_idx'] != -1) &
(df_results_w_array['lower_bound_idx'] != -1) &
(df_results_w_array['upper_bound_idx'] != -1)
].reset_index()
# Check if the dataframe is empty
if df.empty:
ax.set_title(f'{amp}uA {probe} (No Data)')
ax.axis('off')
continue
combined_array = np.zeros((len(df), df.iloc[0]['mean_aligned_sub_data'].shape[0], 800)) #800 for padded
for _, dfrow in df.iterrows():
array = dfrow['mean_aligned_sub_data']
gauss_array = gaussian_filter1d(array, sigma=1, axis=1)
## pad it
max_idx = dfrow['max_idx']
upper_gauss = gauss_array[:, max_idx:]
lower_gauss = gauss_array[:, :max_idx]
if upper_gauss.shape[1] < pad:
upper_pad = [np.nan] * (pad - upper_gauss.shape[1])
padded_upper_gauss = np.zeros((upper_gauss.shape[0], pad))
for i in range(upper_gauss.shape[0]):
padded_upper_gauss[i] = np.concatenate([upper_gauss[i], upper_pad])
if lower_gauss.shape[1] < pad:
lower_pad = [np.nan] * (pad - lower_gauss.shape[1])
padded_lower_gauss = np.zeros((lower_gauss.shape[0], pad))
for i in range(lower_gauss.shape[0]):
padded_lower_gauss[i] = np.concatenate([lower_pad, lower_gauss[i]])
padded_gauss = np.concatenate([padded_lower_gauss, padded_upper_gauss], axis=1)
combined_array[_, :, :] = padded_gauss
mean_array_bycondition[f'{polarity}_{amp}'] = combined_array #3D n (probe, run, etc) x samples x channels
data_to_plot = np.nanmean(combined_array, axis=0).T
data_to_plot = gaussian_filter(data_to_plot, sigma=2)
#data_to_plot = gaussian_filter1d(data_to_plot, sigma=2, axis=1)
#data_to_plot = gaussian_filter1d(data_to_plot, sigma = 2, axis=0)
time_ms = np.linspace(-pre, post, data_to_plot.shape[0])
ax.imshow(data_to_plot, aspect='auto', vmax=200, vmin=-200, origin='lower', cmap='vlag',
extent=[time_ms[0], time_ms[-1], 0, data_to_plot.shape[0]])
ax.set_ylim(300, 500)
ax.set_title(f'{amp}uA {polarity}')
ax.set_xlim(-pre, 3)
if row == 3:
ax.set_xlabel('Time (ms)')
else:
ax.set_xticks([])
if col == 0:
ax.set_ylabel('Channel')
ax.set_yticklabels(np.arange(-100, 101, 50))
else:
ax.set_yticks([])
#path = r'C:\Users\jordan\Documents\combined_figs\evoked_potentials'
#plt.savefig(os.path.join(path, 'avg_heatmaps_by_polarity_amp_mean_smoothed_3ms_100ch_updated.pdf'))
#plt.savefig(os.path.join(path, 'avg_heatmaps_by_polarity_amp_mean_smoothed_3ms_100ch_updated.png'))
# %% [markdown]
# # average EP at 1 ms across distances (Figure 2f and Supp)
# %%
polarities = ['bipolar', 'monopolar']
fig, axs = plt.subplots(4, 4, figsize=(7, 7)) # 3 rows (probes), 4 columns (amplitudes)
fig.suptitle(f'1ms aligned spatial response by polarities', fontsize=14)
probes = ['probeA', 'probeB', 'probeC']
amps = [-5, -25, -50, -100, 5, 25, 50, 100]
mean_boundaries_bycondition = {}
for polarity in polarities:
for i, amp in enumerate(amps):
if amp < 0 and polarity == 'bipolar':
row = 0
elif amp > 0 and polarity == 'bipolar':
row = 1
elif amp < 0 and polarity == 'monopolar':
row = 2
elif amp > 0 and polarity == 'monopolar':
row = 3
if np.abs(amp) == 5:
col = 0
elif np.abs(amp) == 25:
col = 1
elif np.abs(amp) == 50:
col = 2
elif np.abs(amp) == 100:
col = 3
ax = axs[row, col]
# Filter data for the current probe and amplitude
df = df_results[
(df_results['max_volt'] > 40) &
(df_results['amplitude'] == amp) &
(df_results['polarity'] == polarity) &
(df_results['half_max_lower_idx'] != -1) &
(df_results['half_max_upper_idx'] != -1) &
(df_results['lower_bound_idx'] != -1) &
(df_results['upper_bound_idx'] != -1)
]
# Check if the dataframe is empty
if df.empty:
ax.set_title(f'{amp}uA {probe} (No Data)')
ax.axis('off')
continue
# Create a matrix of aligned responses
aligned_matrix = np.array(df['aligned_mean_gauss'].tolist())
mean_aligned_matrix = np.nanmean(aligned_matrix, axis=0)
std_aligned_matrix = np.nanstd(aligned_matrix, axis=0)
sem_aligned_matrix = std_aligned_matrix / np.sqrt(np.sum(~np.isnan(aligned_matrix), axis=0)) # SEM
normalized_chs = np.arange(-400, 400, 1)
normalized_distances = np.arange(-4000, 4000, 10)
upper_bound_idx = int(np.mean(df['upper_bound_idx_aligned'].values))
lower_bound_idx = int(np.mean(df['lower_bound_idx_aligned'].values))
half_max_upper_idx = int(np.mean(df['half_max_upper_idx_aligned'].values))
half_max_lower_idx = int(np.mean(df['half_max_lower_idx_aligned'].values))
mean_boundaries_bycondition[f'{polarity}_{amp}'] = {
'upper_bound': normalized_distances[upper_bound_idx],
'lower_bound': normalized_distances[lower_bound_idx],
'half_max_upper': normalized_distances[half_max_upper_idx],
'half_max_lower': normalized_distances[half_max_lower_idx],
'mean_max': mean_aligned_matrix[400],
'mean_AUC_fullEP': np.trapz(mean_aligned_matrix[lower_bound_idx:upper_bound_idx]),
'mean_AUC_halfEP': np.trapz(mean_aligned_matrix[half_max_lower_idx:half_max_upper_idx])
}
ax.axhline(normalized_distances[lower_bound_idx], color='blue', linestyle='--', label='Slope Lower Bound')
ax.axhline(normalized_distances[upper_bound_idx], color='blue', linestyle='--', label='Slope Upper Bound')
ax.axhline(normalized_distances[half_max_lower_idx], color = 'orange', linestyle = '--', label = 'Half Max Lower')
ax.axhline(normalized_distances[half_max_upper_idx], color = 'orange', linestyle = '--', label = 'Half Max Upper')
# Plot the SEM error cloud
ax.fill_betweenx(
normalized_distances,
mean_aligned_matrix - sem_aligned_matrix,
mean_aligned_matrix + sem_aligned_matrix,
color='gray', alpha=0.3, label='SEM'
)
# Plot the mean aligned matrix
ax.plot(mean_aligned_matrix, normalized_distances, color='black', linewidth=2, label='Mean')
# Set titles and labels
if row == 0:
ax.set_title(f'{amp}uA')
if col == 0:
ax.set_ylabel(f'{polarity}', fontsize=11)
else:
ax.set_yticks([]) # Hide y-ticks for all columns except the first
if row == 3:
ax.set_xlabel('uV', fontsize=8)
else:
ax.set_xticks([]) # Hide x-ticks for all rows except the last
# Add legend only to the last subplot
if row == 3 and col == 3:
ax.legend(loc='upper right', fontsize=8)
ax.set_xlim([0, 800])
ax.set_ylim([-1500, 1500])
plt.tight_layout()
path = os.path.join(r'C:\Users\jordan\Documents\combined_figs\evoked_potentials')
#plt.savefig(os.path.join(path, 'spatial_response_by_polarity.png'))
#plt.savefig(os.path.join(path, 'spatial_response_by_polarity.pdf'))
# %%
mean_boundaries_bycondition = {}
fig, axs = plt.subplots(1, 3, figsize=(8, 3)) # 2 rows for polarities, 2 columns for amplitude sign
fig.suptitle('1ms aligned spatial response by polarities', fontsize=12)
colors = [cathodal_color, anodal_color, bipolar_color]
amps = [-5, -25, -50, -100]
for i, pol in enumerate(['cathodal', 'anodal', 'bipolar']):
ax = axs[i]
if pol == 'bipolar':
polarity = 'bipolar'
else:
polarity = 'monopolar'
for k, amp in enumerate(amps):
if pol == 'anodal':
amp = np.abs(amp) # Make positive for anodal
df = df_results[
(df_results['max_volt'] > 40) &
(df_results['amplitude'] == amp) &
(df_results['polarity'] == polarity) &
(df_results['half_max_lower_idx'] != -1) &
(df_results['half_max_upper_idx'] != -1) &
(df_results['lower_bound_idx'] != -1) &
(df_results['upper_bound_idx'] != -1)
]
if df.empty:
print('skipping empty dataframe for', amp, pol)
continue # Skip this iteration if the dataframe is empty
aligned_matrix = np.array(df['aligned_mean_gauss'].tolist())
mean_aligned_matrix = np.nanmean(aligned_matrix, axis=0)
sem_aligned_matrix = np.nanstd(aligned_matrix, axis=0) / np.sqrt(np.sum(~np.isnan(aligned_matrix), axis=0))
normalized_distances = np.arange(-4000, 4000, 10)
# Plot SEM error cloud
ax.fill_betweenx(
normalized_distances,
mean_aligned_matrix - sem_aligned_matrix,
mean_aligned_matrix + sem_aligned_matrix,
color='gray', alpha=0.2 # Use different colors for each amplitude
)
# Plot the mean response
ax.plot(mean_aligned_matrix, normalized_distances, label=f'{amp}uA', color=colors[i], alpha = 0.3 + (k * 0.2))
ax.set_xlim([0, 800])
ax.set_ylim([-1000, 1000])
ax.set_title(f'{pol}')
ax.set_xlabel('uV')
ax.set_ylabel('Distance (um)')
ax.legend(loc='upper right', fontsize=8)
# Adjust layout and save figures
plt.tight_layout()
#path = r'C:\Users\jordan\Documents\combined_figs\evoked_potentials'
#plt.savefig(os.path.join(path, 'spatial_response_by_polarity_and_sign_100ch.png'))
#plt.savefig(os.path.join(path, 'spatial_response_by_polarity_and_sign_100ch.pdf'))
# %% [markdown]
# # spatial extents superficial and deep (Figure 2g and Supp CCF vs peak-aligned)
# %%
# ccf distances and peak aligned distances plots
from scipy.stats import sem
from scipy.stats import linregress
for dist_type in ['ccf', 'peak_aligned']:
fig, axs = plt.subplots(1, 3, figsize=(8, 3)) # 2 rows for polarities, 2 columns for amplitude sign
fig.suptitle(f'{dist_type} distances 1D spatial bounds (up vs down) w fits', fontsize=12)
colors = [cathodal_color, anodal_color, bipolar_color]
amps = [-5, -25, -50, -100]
for i, pol in enumerate(['cathodal', 'anodal', 'bipolar']):
ax = axs[i]
ax.set_title(pol, color = colors[i])
if pol == 'bipolar':
polarity = 'bipolar'
else:
polarity = 'monopolar'
all_upper_dists = []
all_lower_dists = []
all_amps = []
for k, amp in enumerate(amps):
if pol == 'anodal':
amp = np.abs(amp) # Make positive for anodal
ax.set_xlabel('Amplitude (uA)')
ax.set_ylabel('Distance (units)')
ax.invert_yaxis()
df = df_results[
(df_results['max_volt'] > 40) &
(df_results['amplitude'] == amp) &
(df_results['polarity'] == polarity) &
(df_results['half_max_lower_idx'] != -1) &
(df_results['half_max_upper_idx'] != -1) &
(df_results['lower_bound_idx'] != -1) &
(df_results['upper_bound_idx'] != -1)
]
if dist_type == 'ccf':
upper_dists = df['upper_bound_dist'].values
lower_dists = df['lower_bound_dist'].values
elif dist_type == 'peak_aligned':
normalized_distances = np.arange(-4000, 4000, 10)[::-1]
upper_bound_idxs = df['upper_bound_idx_aligned'].values
lower_bound_idxs = df['lower_bound_idx_aligned'].values
upper_dists = [normalized_distances[idx] for idx in upper_bound_idxs]
lower_dists = [normalized_distances[idx] for idx in lower_bound_idxs]
# Collect data for scatter plot
all_upper_dists.extend(upper_dists)
all_lower_dists.extend(lower_dists)
all_amps.extend([amp] * len(upper_dists))
# Scatter individual data points
ax.scatter([amp] * len(upper_dists), upper_dists, color=colors[i], s = 3, alpha=0.8, zorder=2)
ax.scatter([amp] * len(lower_dists), lower_dists, color=colors[i], s = 3, alpha=0.8, zorder=2)
# print mean and std of upper and lower distances
mean_upper = np.mean(all_upper_dists)
mean_lower = np.mean(all_lower_dists)
std_upper = np.std(all_upper_dists)
std_lower = np.std(all_lower_dists)
print(f'{pol} {dist_type} {amp}uA: Upper Mean: {mean_upper:.2f}, Std: {std_upper:.2f}; Lower Mean: {mean_lower:.2f}, Std: {std_lower:.2f}')
# Fit polynomial regression to the data
if len(all_amps) > 1: # Ensure we have enough data points
upper_coeffs = np.polyfit(all_amps, all_upper_dists, deg=2)
lower_coeffs = np.polyfit(all_amps, all_lower_dists, deg=2)
upper_poly = np.poly1d(upper_coeffs)
lower_poly = np.poly1d(lower_coeffs)
# Generate smooth curves for plotting
amp_range = np.linspace(np.min(all_amps), np.max(all_amps), 100)
upper_fit = upper_poly(amp_range)
lower_fit = lower_poly(amp_range)
# Calculate standard error of the mean (SEM) for the fits
upper_fit_err = sem(all_upper_dists)
lower_fit_err = sem(all_lower_dists)
# Plot polynomial regression lines
ax.plot(amp_range, upper_fit, color='darkblue', label='Upper Bound Fit', zorder=4)
ax.plot(amp_range, lower_fit, color='darkred', label='Lower Bound Fit', zorder=4)
# Add error bands
ax.fill_between(amp_range, upper_fit - upper_fit_err, upper_fit + upper_fit_err, color='blue', alpha=0.2, label='Upper Fit Error', zorder=1)
ax.fill_between(amp_range, lower_fit - lower_fit_err, lower_fit + lower_fit_err, color='red', alpha=0.2, label='Lower Fit Error', zorder=1)
x = all_amps
y_upper = all_upper_dists
y_lower = all_lower_dists
# Linear fit (Upper)
slope_u, intercept_u, r_u, pval_u, stderr_u = linregress(x, y_upper)
r2_linear_upper = r_u ** 2
print(f"Linear Fit (Upper): R² = {r2_linear_upper:.2f}, Slope = {slope_u:.2f}, Intercept = {intercept_u:.2f}, p = {pval_u:.2e}")
# Quadratic fit (Upper)
coeffs_quad_upper = np.polyfit(x, y_upper, 2)
p_quad_upper = np.poly1d(coeffs_quad_upper)
y_fit_quad_upper = p_quad_upper(x)
r2_quad_upper = 1 - (np.sum((y_upper - y_fit_quad_upper)**2) / np.sum((y_upper - np.mean(y_upper))**2))
print(f"Quadratic Fit (Upper): R² = {r2_quad_upper:.2f}, Coeffs = {coeffs_quad_upper}")
# Linear fit (Lower)
slope_l, intercept_l, r_l, pval_l, stderr_l = linregress(x, y_lower)
r2_linear_lower = r_l ** 2
print(f"Linear Fit (Lower): R² = {r2_linear_lower:.2f}, Slope = {slope_l:.2f}, Intercept = {intercept_l:.2f}, p = {pval_l:.2e}")
# Quadratic fit (Lower)
coeffs_quad_lower = np.polyfit(x, y_lower, 2)
p_quad_lower = np.poly1d(coeffs_quad_lower)
y_fit_quad_lower = p_quad_lower(x)
r2_quad_lower = 1 - (np.sum((y_lower - y_fit_quad_lower)**2) / np.sum((y_lower - np.mean(y_lower))**2))
print(f"Quadratic Fit (Lower): R² = {r2_quad_lower:.2f}, Coeffs = {coeffs_quad_lower}")
# Add color-coded text annotations for R² values
ax.text(
0.6, 0.90,
f"Quad: R²: {r2_quad_upper:.2f}",
color='blue', fontsize=5, verticalalignment='top', transform=ax.transAxes
)
ax.text(
0.6, 0.80,
f"Lin: R²: {r2_linear_upper:.2f}",
color='darkblue', fontsize=5, verticalalignment='top', transform=ax.transAxes
)
ax.text(
0.6, 0.10,
f"Quad: R²: {r2_quad_lower:.2f}",
color='red', fontsize=5, verticalalignment='top', transform=ax.transAxes)
ax.text(
0.6, 0.20,
f"Lin: R²: {r2_linear_lower:.2f}",
color='darkred', fontsize=5, verticalalignment='top', transform=ax.transAxes)
if pol == 'anodal':
ax.set_xticks(np.abs(amps))
else:
ax.set_xticks(amps)
ax.invert_xaxis() # Invert x-axis for negative amplitudes
ax.invert_yaxis()
# Adjust layout and save the figure
plt.tight_layout(rect=[0, 0, 1, 0.96])
path = r'C:\Users\jordan\Documents\combined_figs\evoked_potentials'
#plt.savefig(os.path.join(path, f'distances_{dist_type}_with_fit_updated.png'))
#plt.savefig(os.path.join(path, f'distances_{dist_type}_with_fit_updated.pdf'))
# %%
# ccf distances and peak aligned distances plots
from scipy.stats import sem, linregress
for dist_type in ['ccf', 'peak_aligned']:
fig, axs = plt.subplots(1, 3, figsize=(8, 3))
fig.suptitle(f'{dist_type} distances 1D spatial bounds (up vs down) — linear fits', fontsize=12)
colors = [cathodal_color, anodal_color, bipolar_color]
amps = [-5, -25, -50, -100]
for i, pol in enumerate(['cathodal', 'anodal', 'bipolar']):
ax = axs[i]
ax.set_title(pol, color=colors[i])
polarity = 'bipolar' if pol == 'bipolar' else 'monopolar'
all_upper_dists = []
all_lower_dists = []
all_amps = []
mean_upper_per_amp = []
mean_lower_per_amp = []
sem_upper_per_amp = []
sem_lower_per_amp = []
for k, amp in enumerate(amps):
if pol == 'anodal':
amp = np.abs(amp)
df = df_results[
(df_results['max_volt'] > 40) &
(df_results['amplitude'] == amp) &
(df_results['polarity'] == polarity) &
(df_results['half_max_lower_idx'] != -1) &
(df_results['half_max_upper_idx'] != -1) &
(df_results['lower_bound_idx'] != -1) &
(df_results['upper_bound_idx'] != -1)
]
if dist_type == 'ccf':
upper_dists = df['upper_bound_dist'].values
lower_dists = df['lower_bound_dist'].values
else:
normalized_distances = np.arange(-4000, 4000, 10)[::-1]
upper_dists = [normalized_distances[idx] for idx in df['upper_bound_idx_aligned'].values]
lower_dists = [normalized_distances[idx] for idx in df['lower_bound_idx_aligned'].values]
# scatter points
ax.scatter([amp] * len(upper_dists), upper_dists, color=colors[i], s=3, alpha=0.6, zorder=2)
ax.scatter([amp] * len(lower_dists), lower_dists, color=colors[i], s=3, alpha=0.6, zorder=2)
# accumulate
all_upper_dists.extend(upper_dists)
all_lower_dists.extend(lower_dists)
all_amps.extend([amp] * len(upper_dists))
# store mean ± sem for plotting mean dots
mean_upper_per_amp.append(np.mean(upper_dists))
mean_lower_per_amp.append(np.mean(lower_dists))
sem_upper_per_amp.append(sem(upper_dists))
sem_lower_per_amp.append(sem(lower_dists))
print(f'{pol} {dist_type} {amp}uA: upper mean={np.mean(upper_dists):.2f}, lower mean={np.mean(lower_dists):.2f}')
# plot means with error bars
ax.errorbar(amps if pol != 'anodal' else np.abs(amps),
mean_upper_per_amp, yerr=sem_upper_per_amp,
fmt='o', color='lightgray', markersize=4, label='upper mean', zorder=4)
ax.errorbar(amps if pol != 'anodal' else np.abs(amps),
mean_lower_per_amp, yerr=sem_lower_per_amp,
fmt='o', color='black', markersize=4, label='lower mean', zorder=4)
# linear fits
if len(all_amps) > 1:
x = np.asarray(all_amps)
y_upper = np.asarray(all_upper_dists)
y_lower = np.asarray(all_lower_dists)
# linear regressions
slope_u, intercept_u, r_u, p_u, _ = linregress(x, y_upper)
slope_l, intercept_l, r_l, p_l, _ = linregress(x, y_lower)
r2_u = r_u ** 2
r2_l = r_l ** 2
x_line = np.linspace(np.min(x), np.max(x), 100)
ax.plot(x_line, slope_u * x_line + intercept_u, color='lightgray', linewidth=1.2)
ax.plot(x_line, slope_l * x_line + intercept_l, color='black', linewidth=1.2)
print(f"{pol} upper linear: R²={r2_u:.2f}, slope={slope_u:.2f}, p={p_u:.2e}")
print(f"{pol} lower linear: R²={r2_l:.2f}, slope={slope_l:.2f}, p={p_l:.2e}")
# annotate
ax.text(0.60, 0.85, f"Upper R²={r2_u:.2f}", color='lightgray', fontsize=6, va='top', transform=ax.transAxes)
ax.text(0.60, 0.20, f"Lower R²={r2_l:.2f}", color='black', fontsize=6, va='top', transform=ax.transAxes)
# aesthetics
if pol == 'anodal':
ax.set_xticks(np.abs(amps))
else:
ax.set_xticks(amps)
ax.invert_xaxis()
ax.set_xlabel('amplitude (uA)')
ax.set_ylabel('boundary distance')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.invert_yaxis()
ax.set_ylim([1500, -1500])
ax.set_yticks([1500, 1000, 500, 0, -500, -1000, -1500], ['1500', '', '', '0', '', '', '-1500'])
plt.tight_layout(rect=[0, 0, 1, 0.96])
path = r'C:\Users\jordan\Documents\combined_figs\evoked_potentials'
# plt.savefig(os.path.join(path, f'distances_{dist_type}_linear_means.png'))
# plt.savefig(os.path.join(path, f'distances_{dist_type}_linear_means.pdf'))
# %%
# full version: signed means ± SEM (clouds) and absolute linear fits (offset scatters)
from scipy.stats import sem, linregress
for dist_type in ['ccf', 'peak_aligned']:
colors = [cathodal_color, anodal_color, bipolar_color]
amps_base = [-5, -25, -50, -100]
# ========== PLOT 1: signed means ± SEM as error "clouds" with individual points ==========
fig1, axs1 = plt.subplots(1, 3, figsize=(8, 3))
fig1.suptitle(f'{dist_type} — signed boundary distances (mean ± SEM clouds)', fontsize=12)
for i, pol in enumerate(['cathodal', 'anodal', 'bipolar']):
ax = axs1[i]
ax.set_title(pol, color=colors[i])
polarity = 'bipolar' if pol == 'bipolar' else 'monopolar'
amps = np.abs(amps_base) if pol == 'anodal' else amps_base
mean_upper, mean_lower = [], []
sem_upper, sem_lower = [], []
for amp_in in amps_base:
amp = abs(amp_in) if pol == 'anodal' else amp_in
df = df_results[
(df_results['max_volt'] > 40) &
(df_results['amplitude'] == amp) &
(df_results['polarity'] == polarity) &
(df_results['half_max_lower_idx'] != -1) &
(df_results['half_max_upper_idx'] != -1) &
(df_results['lower_bound_idx'] != -1) &
(df_results['upper_bound_idx'] != -1)
]
if dist_type == 'ccf':
upper_dists = df['upper_bound_dist'].values
lower_dists = df['lower_bound_dist'].values
else:
normalized_distances = np.arange(-4000, 4000, 10)[::-1]
upper_dists = [normalized_distances[idx] for idx in df['upper_bound_idx_aligned'].values]
lower_dists = [normalized_distances[idx] for idx in df['lower_bound_idx_aligned'].values]
ax.scatter([amp]*len(upper_dists), upper_dists, color=colors[i], s=4, alpha=0.4, zorder=1)
ax.scatter([amp]*len(lower_dists), lower_dists, color=colors[i], s=4, alpha=0.4, zorder=1)
mean_upper.append(np.mean(upper_dists))
mean_lower.append(np.mean(lower_dists))
sem_upper.append(sem(upper_dists))
sem_lower.append(sem(lower_dists))
# shaded SEM "clouds"
ax.plot(amps, mean_upper, color='gray', lw=1.2)
ax.fill_between(amps,
np.array(mean_upper)-np.array(sem_upper),
np.array(mean_upper)+np.array(sem_upper),
color='gray', alpha=0.25, label='upper ± sem')
ax.plot(amps, mean_lower, color='black', lw=1.2)
ax.fill_between(amps,
np.array(mean_lower)-np.array(sem_lower),
np.array(mean_lower)+np.array(sem_lower),
color='black', alpha=0.25, label='lower ± sem')
# formatting
if pol == 'anodal':
ax.set_xticks(np.abs(amps_base))
else:
ax.set_xticks(amps_base)
ax.invert_xaxis()
ax.set_xlabel('amplitude (µA)')
if i == 0:
ax.set_ylabel('boundary distance (µm)')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.invert_yaxis()
ax.set_ylim([1500, -1500])
ax.set_yticks([1500, 1000, 500, 0, -500, -1000, -1500],
['1500', '', '', '0', '', '', '-1500'])
if i == 0:
ax.legend(frameon=False, fontsize=7, loc='lower left')
plt.tight_layout(rect=[0, 0, 1, 0.95])
#plt.savefig(os.path.join(path, f'{dist_type}_means_sem_cloud.png'))
#plt.savefig(os.path.join(path, f'{dist_type}_means_sem_cloud.pdf'))
# ========== PLOT 2: absolute distances with linear fits + offset scatters ==========
fig2, axs2 = plt.subplots(1, 3, figsize=(8, 3))
fig2.suptitle(f'{dist_type} — absolute boundary distances (linear fits + offset scatters)', fontsize=12)
offset = 3.0 # µA offset between upper and lower scatters
for i, pol in enumerate(['cathodal', 'anodal', 'bipolar']):
ax = axs2[i]
ax.set_title(pol, color=colors[i])
polarity = 'bipolar' if pol == 'bipolar' else 'monopolar'
amps = np.abs(amps_base) if pol == 'anodal' else amps_base
x_upper, y_upper = [], []
x_lower, y_lower = [], []
for amp_in in amps_base:
amp = abs(amp_in) if pol == 'anodal' else amp_in
df = df_results[
(df_results['max_volt'] > 40) &
(df_results['amplitude'] == amp) &
(df_results['polarity'] == polarity) &
(df_results['half_max_lower_idx'] != -1) &
(df_results['half_max_upper_idx'] != -1) &
(df_results['lower_bound_idx'] != -1) &
(df_results['upper_bound_idx'] != -1)
]
if dist_type == 'ccf':
upper_dists = df['upper_bound_dist'].values
lower_dists = df['lower_bound_dist'].values
else:
normalized_distances = np.arange(-4000, 4000, 10)[::-1]
upper_dists = [normalized_distances[idx] for idx in df['upper_bound_idx_aligned'].values]
lower_dists = [normalized_distances[idx] for idx in df['lower_bound_idx_aligned'].values]
au = np.abs(upper_dists)
al = np.abs(lower_dists)
x_upper.extend([amp]*len(au))
y_upper.extend(au)
x_lower.extend([amp]*len(al))
y_lower.extend(al)
# offset scatter positions
ax.scatter(np.array([amp]*len(au)) + offset, au,
marker='o', s=24, color='gray', alpha=0.6, edgecolor='none', zorder=2, label=None)
ax.scatter(np.array([amp]*len(al)) - offset, al,
marker='s', s=24, color='black', alpha=0.6, edgecolor='none', zorder=2, label=None)
# regression fits
if len(x_upper) > 1 and len(x_lower) > 1:
xu, yu = np.array(x_upper, float), np.array(y_upper, float)
xl, yl = np.array(x_lower, float), np.array(y_lower, float)
su, bu, ru, pu, _ = linregress(xu, yu)
sl, bl, rl, pl, _ = linregress(xl, yl)
x_line = np.linspace(np.min(amps), np.max(amps), 100)
ax.plot(x_line, sl * x_line + bl, color='black', lw=1.5, ls='-', label='|lower| fit')
ax.plot(x_line, su * x_line + bu, color='gray', lw=1.5, ls='--', label='|upper| fit')
print(f"{pol} {dist_type}: upper R²={ru**2:.2f}, lower R²={rl**2:.2f}")
if pol == 'anodal':
ax.set_xticks(np.abs(amps_base))
else:
ax.set_xticks(amps_base)
ax.invert_xaxis()
ax.set_xlabel('amplitude (µA)')
if i == 0:
ax.set_ylabel('absolute boundary distance (µm)')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set_ylim(0, 1500)
ax.set_yticks([0, 500, 1000, 1500])
if i == 0:
ax.legend(frameon=False, fontsize=7, loc='upper left')
plt.tight_layout(rect=[0, 0, 1, 0.95])
path = r'C:\Users\jordan\Documents\combined_figs\evoked_potentials'
#plt.savefig(os.path.join(path, f'{dist_type}_abs_linear_fit_offset.png'))
#plt.savefig(os.path.join(path, f'{dist_type}_abs_linear_fit_offset.pdf'))
# %%
# config
dist_type = 'ccf' # 'ccf' or 'peak_aligned'
title = f'50 µA radius across pols'
# helper: pull combined absolute distances (upper + lower) at 50 µA for one polarity
def _pull_combined_50(df_results, dist_type, pol_name):
# map to df polarity/amp
if pol_name == 'bipolar':
polarity = 'bipolar'; amp = -50
elif pol_name == 'anodal':
polarity = 'monopolar'; amp = 50
else:
polarity = 'monopolar'; amp = -50 # cathodal
df = df_results[
(df_results['max_volt'] > 40) &
(df_results['amplitude'] == amp) &
(df_results['polarity'] == polarity) &
(df_results['half_max_lower_idx'] != -1) &
(df_results['half_max_upper_idx'] != -1) &
(df_results['lower_bound_idx'] != -1) &
(df_results['upper_bound_idx'] != -1)
].copy()
# get distances per side
if dist_type == 'ccf':
upper = df['upper_bound_dist'].to_numpy(dtype=float)
lower = df['lower_bound_dist'].to_numpy(dtype=float)
else: # 'peak_aligned'
# aligned distances are in 10 µm steps from -4000..+3990 reversed
normalized_distances = np.arange(-4000, 4000, 10)[::-1]
upper = np.array([normalized_distances[idx] for idx in df['upper_bound_idx_aligned'].to_numpy()], dtype=float)
lower = np.array([normalized_distances[idx] for idx in df['lower_bound_idx_aligned'].to_numpy()], dtype=float)
# combine as magnitudes; n doubles (upper+lower treated as separate observations)
pooled = np.concatenate([np.abs(upper), np.abs(lower)])
return pooled
# build tidy df
rows = []
for pol in ['cathodal', 'anodal', 'bipolar']:
vals = _pull_combined_50(df_results, dist_type, pol)
rows += [{'polarity': pol, 'dist_abs': v} for v in vals]
plot_df = pd.DataFrame(rows)
# plot: strip + mean±sem in your stim colors
fig, ax = plt.subplots(figsize=(2, 3.0))
ax.set_title(title, fontsize=9)
ax.set_xlabel('polarity (50 µA overall)')
ax.set_ylabel('distance (µm)')
for pol, col in [('cathodal', cathodal_color), ('anodal', anodal_color), ('bipolar', bipolar_color)]:
sub = plot_df[plot_df['polarity'] == pol]
sns.stripplot(
data=sub, x='polarity', y='dist_abs',
jitter=0.25, alpha=0.65, size=3,
color=col, ax=ax
)
# overlay mean ± sem
g = plot_df.groupby('polarity')['dist_abs']
xcats = ['cathodal', 'anodal', 'bipolar']
xlocs = np.arange(len(xcats))
means = [g.get_group(cat).mean() for cat in xcats]
sems = [g.get_group(cat).sem() for cat in xcats]
for xi, cat, m, s in zip(xlocs, xcats, means, sems):
ax.errorbar(
xi, m, yerr=s, fmt='D-', lw=1.2, ms=5,
color={'cathodal': cathodal_color, 'anodal': anodal_color, 'bipolar': bipolar_color}[cat],
capsize=3, zorder=3
)
ax.set_yticks([0, 1000, 2000])
# styling
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
# optional fixed ylim for comparability:
# ax.set_ylim(0, 1500)
plt.tight_layout()
path = r'C:\Users\jordan\Documents\combined_figs\evoked_potentials'
#plt.savefig(os.path.join(path, f'ccf_dist_50ua_avgdist.png'))
#plt.savefig(os.path.join(path, f'ccf_dist_50ua_avgdist.pdf'))
# %%
import pandas as pd
import numpy as np
from statsmodels.formula.api import ols
import statsmodels.api as sm
# Filter for valid ccf entries
df_ccf = df_results[
(df_results['max_volt'] > 40) &
(df_results['half_max_lower_idx'] != -1) &
(df_results['half_max_upper_idx'] != -1) &
(df_results['lower_bound_idx'] != -1) &
(df_results['upper_bound_idx'] != -1) &
(df_results['amplitude'] != 10) &
(df_results['amplitude'] != -10) # Exclude 10 uA
]
# Build long-form dataframe
records = []
for _, row in df_ccf.iterrows():
for bound_type, dist_col in zip(['upper', 'lower'], ['upper_bound_dist', 'lower_bound_dist']):
records.append({
'distance': row[dist_col],
'amplitude': np.abs(row['amplitude']),
'polarity': (
'bipolar' if row['polarity'] == 'bipolar'
else 'anodal' if row['amplitude'] > 0
else 'cathodal'
),
'bound_type': bound_type
})
df_long = pd.DataFrame(records)
df_long['amplitude'] = df_long['amplitude'].astype('category')
df_long['polarity'] = df_long['polarity'].astype('category')
df_long['bound_type'] = df_long['bound_type'].astype('category')
for bound in ['upper', 'lower']:
print(f"\n---- {bound.upper()} BOUND ----")
df_bound = df_long[df_long['bound_type'] == bound]
model = ols('distance ~ C(polarity) * amplitude', data=df_bound).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)
model = ols('distance ~ C(polarity) * amplitude * C(bound_type)', data=df_long).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)
# %%
# Treat amplitude as continuous — make sure it's numeric
df_cathodal['amplitude'] = df_cathodal['amplitude'].astype(float)
# Run model with continuous amplitude and interaction with bound_type
model = ols('distance ~ amplitude * C(bound_type)', data=df_cathodal).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)
# %% [markdown]
# # half spatial extents superficial and deep (Figure 2g and Supp CCF vs peak-aligned)
# %%
# ccf distances and peak aligned distances plots
from scipy.stats import sem
for dist_type in ['ccf', 'peak_aligned']:
fig, axs = plt.subplots(1, 3, figsize=(8, 3)) # 2 rows for polarities, 2 columns for amplitude sign
fig.suptitle(f'{dist_type} distances 1D half spatial bounds (up vs down) w fits', fontsize=12)
colors = [cathodal_color, anodal_color, bipolar_color]
for i, pol in enumerate(['cathodal', 'anodal', 'bipolar']):
ax = axs[i]
ax.set_title(pol, color = colors[i])
if pol == 'bipolar':
polarity = 'bipolar'
else:
polarity = 'monopolar'
all_upper_dists = []
all_lower_dists = []
all_amps = []
all_mean_uppers = []
all_mean_lowers = []
all_std_uppers = []
all_std_lowers = []
all_sem_uppers = []
all_sem_lowers = []
amps = [-5, -25, -50, -100]
for k, amp in enumerate(amps):
if pol == 'anodal':
amp = np.abs(amp) # Make positive for anodal
ax.set_xlabel('Amplitude (uA)')
ax.set_ylabel('Distance (units)')
#ax.invert_yaxis()
df = df_results[
(df_results['max_volt'] > 40) &
(df_results['amplitude'] == amp) &
(df_results['polarity'] == polarity) &
(df_results['half_max_lower_idx'] != -1) &
(df_results['half_max_upper_idx'] != -1) &
(df_results['lower_bound_idx'] != -1) &
(df_results['upper_bound_idx'] != -1)
]
if dist_type == 'ccf':
upper_dists = df['half_max_upper_dist'].values
lower_dists = df['half_max_lower_dist'].values
elif dist_type == 'peak_aligned':
normalized_distances = np.arange(-4000, 4000, 10)[::-1]
upper_bound_idxs = df['half_max_upper_idx_aligned'].values
lower_bound_idxs = df['half_max_lower_idx_aligned'].values
upper_dists = [normalized_distances[idx] for idx in upper_bound_idxs]
lower_dists = [normalized_distances[idx] for idx in lower_bound_idxs]
# Collect data for scatter plot
all_upper_dists.extend(upper_dists)
all_lower_dists.extend(lower_dists)
all_amps.extend([amp] * len(upper_dists))
# Scatter individual data points
ax.scatter([amp] * len(upper_dists), upper_dists, color=colors[i], s = 3, alpha=0.8, zorder=2)
ax.scatter([amp] * len(lower_dists), lower_dists, color=colors[i], s = 3, alpha=0.8, zorder=2)
# print mean and std of upper and lower distances
mean_upper = np.mean(all_upper_dists)
mean_lower = np.mean(all_lower_dists)
std_upper = np.std(all_upper_dists)
std_lower = np.std(all_lower_dists)
sem_upper = sem(all_upper_dists)
sem_lower = sem(all_lower_dists)
all_mean_uppers.append(mean_upper)
all_mean_lowers.append(mean_lower)
all_std_uppers.append(std_upper)
all_std_lowers.append(std_lower)
all_sem_uppers.append(sem_upper)
all_sem_lowers.append(sem_lower)
print(f'{pol} {dist_type} {amp}uA: Upper Mean: {mean_upper:.2f}, Std: {std_upper:.2f}; Lower Mean: {mean_lower:.2f}, Std: {std_lower:.2f}')
if pol == 'anodal':
amps = [5, 25, 50, 100]
else:
amps = [-5, -25, -50, -100]
ax.plot(amps, all_mean_uppers, color = colors[i], marker = 'o', markersize = 4, linewidth = 1)
ax.plot(amps, all_mean_lowers, color = colors[i], marker = 'o', markersize = 4, linewidth = 1)
ax.fill_between(amps,
np.array(all_mean_uppers) - np.array(all_sem_uppers),
np.array(all_mean_uppers) + np.array(all_sem_uppers),
color=colors[i], alpha=0.2, label='Upper Bound SEM', zorder=1)
ax.fill_between(amps,
np.array(all_mean_lowers) - np.array(all_sem_lowers),
np.array(all_mean_lowers) + np.array(all_sem_lowers),
color=colors[i], alpha=0.2, label='Lower Bound SEM', zorder=1)
# Fit polynomial regression to the data
if len(all_amps) > 1: # Ensure we have enough data points
upper_coeffs = np.polyfit(all_amps, all_upper_dists, deg=2)
lower_coeffs = np.polyfit(all_amps, all_lower_dists, deg=2)
upper_poly = np.poly1d(upper_coeffs)
lower_poly = np.poly1d(lower_coeffs)
# Generate smooth curves for plotting
amp_range = np.linspace(np.min(all_amps), np.max(all_amps), 100)
upper_fit = upper_poly(amp_range)
lower_fit = lower_poly(amp_range)
# Calculate standard error of the mean (SEM) for the fits
upper_fit_err = sem(all_upper_dists)
lower_fit_err = sem(all_lower_dists)
# Plot polynomial regression lines
#ax.plot(amp_range, upper_fit, color='darkblue', label='Upper Bound Fit', zorder=4)
#ax.plot(amp_range, lower_fit, color='darkred', label='Lower Bound Fit', zorder=4)
# Add error bands
#ax.fill_between(amp_range, upper_fit - upper_fit_err, upper_fit + upper_fit_err, color='blue', alpha=0.2, label='Upper Fit Error', zorder=1)
#ax.fill_between(amp_range, lower_fit - lower_fit_err, lower_fit + lower_fit_err, color='red', alpha=0.2, label='Lower Fit Error', zorder=1)
x = all_amps
y_upper = all_upper_dists
y_lower = all_lower_dists
# Fit linear regression for upper bounds
coeffs_linear_upper = np.polyfit(x, y_upper, 1)
p_linear_upper = np.poly1d(coeffs_linear_upper)
y_fit_linear_upper = p_linear_upper(x)
r2_linear_upper = 1 - (np.sum((y_upper - y_fit_linear_upper)**2) / np.sum((y_upper - np.mean(y_upper))**2))
# Fit quadratic regression for upper bounds
coeffs_quad_upper = np.polyfit(x, y_upper, 2)
p_quad_upper = np.poly1d(coeffs_quad_upper)
y_fit_quad_upper = p_quad_upper(x)
r2_quad_upper = 1 - (np.sum((y_upper - y_fit_quad_upper)**2) / np.sum((y_upper - np.mean(y_upper))**2))
# Fit linear regression for lower bounds
coeffs_linear_lower = np.polyfit(x, y_lower, 1)
p_linear_lower = np.poly1d(coeffs_linear_lower)
y_fit_linear_lower = p_linear_lower(x)
r2_linear_lower = 1 - (np.sum((y_lower - y_fit_linear_lower)**2) / np.sum((y_lower - np.mean(y_lower))**2))
# Fit quadratic regression for lower bounds
coeffs_quad_lower = np.polyfit(x, y_lower, 2)
p_quad_lower = np.poly1d(coeffs_quad_lower)
y_fit_quad_lower = p_quad_lower(x)
r2_quad_lower = 1 - (np.sum((y_lower - y_fit_quad_lower)**2) / np.sum((y_lower - np.mean(y_lower))**2))
# Add color-coded text annotations for R² values
# Upper bounds
#ax.text(
# 0.6, 0.90,
# f"Quad: R²: {r2_quad_upper:.2f}",
# color='blue', fontsize=5, verticalalignment='top', transform=ax.transAxes
#)
#ax.text(
# 0.6, 0.80,
# f"Lin: R²: {r2_linear_upper:.2f}",
# color='darkblue', fontsize=5, verticalalignment='top', transform=ax.transAxes
#)
#ax.text(
# 0.6, 0.10,
# f"Quad: R²: {r2_quad_lower:.2f}",
# color='red', fontsize=5, verticalalignment='top', transform=ax.transAxes)
#ax.text(
# 0.6, 0.20,
# f"Lin: R²: {r2_linear_lower:.2f}",
# color='darkred', fontsize=5, verticalalignment='top', transform=ax.transAxes)
if pol == 'anodal':
ax.set_xticks(np.abs(amps))
else:
ax.set_xticks(amps)
ax.invert_xaxis() # Invert x-axis for negative amplitudes
#ax.invert_yaxis()
ax.set_ylim([-1000, 1000])
ax.set_yticks([1000, 500, 0, -500, -1000], ['1', '', '0', '', '-1'])
ax.invert_yaxis()
# Adjust layout and save the figure
plt.tight_layout(rect=[0, 0, 1, 0.96])
path = r'C:\Users\jordan\Documents\combined_figs\evoked_potentials'
#plt.savefig(os.path.join(path, f'half_distances_{dist_type}_updated.png'))
#plt.savefig(os.path.join(path, f'half_distances_{dist_type}_updated.pdf'))
# %%
import pandas as pd
import numpy as np
from statsmodels.formula.api import ols
import statsmodels.api as sm
# Filter for valid ccf entries
df_ccf = df_results[
(df_results['max_volt'] > 40) &
(df_results['half_max_lower_idx'] != -1) &
(df_results['half_max_upper_idx'] != -1) &
(df_results['lower_bound_idx'] != -1) &
(df_results['upper_bound_idx'] != -1) &
(df_results['amplitude'] != 10) &
(df_results['amplitude'] != -10) # Exclude 10 uA
]
# Build long-form dataframe
records = []
for _, row in df_ccf.iterrows():
for bound_type, dist_col in zip(['upper', 'lower'], ['half_max_upper_idx_aligned', 'half_max_lower_idx_aligned']):
records.append({
'distance': normalized_distances[row[dist_col]],
'amplitude': np.abs(row['amplitude']),
'polarity': (
'bipolar' if row['polarity'] == 'bipolar'
else 'anodal' if row['amplitude'] > 0
else 'cathodal'
),
'bound_type': bound_type
})
df_long = pd.DataFrame(records)
df_long['amplitude'] = df_long['amplitude'].astype('category')
df_long['polarity'] = df_long['polarity'].astype('category')
df_long['bound_type'] = df_long['bound_type'].astype('category')
for bound in ['upper', 'lower']:
print(f"\n---- {bound.upper()} BOUND ----")
df_bound = df_long[df_long['bound_type'] == bound]
model = ols('distance ~ C(polarity) * amplitude', data=df_bound).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)
model = ols('distance ~ C(polarity) * amplitude * C(bound_type)', data=df_long).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)
# %%
6.040685e-03
# %%
df.columns
# %%
# ccf distances and peak aligned distances plots
from scipy.stats import sem
for dist_type in ['ccf', 'peak_aligned']:
fig, axs = plt.subplots(1, 3, figsize=(8, 3)) # 2 rows for polarities, 2 columns for amplitude sign
fig.suptitle(f'{dist_type} distances 1D half spatial bounds (up vs down) w fits', fontsize=12)
colors = [cathodal_color, anodal_color, bipolar_color]
for i, pol in enumerate(['cathodal', 'anodal', 'bipolar']):
ax = axs[i]
ax.set_title(pol, color = colors[i])
if pol == 'bipolar':
polarity = 'bipolar'
else:
polarity = 'monopolar'
all_upper_dists = []
all_lower_dists = []
all_amps = []
all_mean_uppers = []
all_mean_lowers = []
all_std_uppers = []
all_std_lowers = []
all_sem_uppers = []
all_sem_lowers = []
amps = [-5, -25, -50, -100]
for k, amp in enumerate(amps):
if pol == 'anodal':
amp = np.abs(amp) # Make positive for anodal
ax.set_xlabel('Amplitude (uA)')
ax.set_ylabel('Distance (units)')
#ax.invert_yaxis()
df = df_results[
(df_results['max_volt'] > 40) &
(df_results['amplitude'] == amp) &
(df_results['polarity'] == polarity) &
(df_results['half_max_lower_idx'] != -1) &
(df_results['half_max_upper_idx'] != -1) &
(df_results['lower_bound_idx'] != -1) &
(df_results['upper_bound_idx'] != -1)
]
if dist_type == 'ccf':
upper_dists = df['upper_bound_dist'].values
lower_dists = df['lower_bound_dist'].values
elif dist_type == 'peak_aligned':
normalized_distances = np.arange(-4000, 4000, 10)[::-1]
upper_bound_idxs = df['upper_bound_idx_aligned'].values
lower_bound_idxs = df['lower_bound_idx_aligned'].values
upper_dists = [normalized_distances[idx] for idx in upper_bound_idxs]
lower_dists = [normalized_distances[idx] for idx in lower_bound_idxs]
# Collect data for scatter plot
all_upper_dists.extend(upper_dists)
all_lower_dists.extend(lower_dists)
all_amps.extend([amp] * len(upper_dists))
# Scatter individual data points
ax.scatter([amp] * len(upper_dists), upper_dists, color=colors[i], s = 3, alpha=0.8, zorder=2)
ax.scatter([amp] * len(lower_dists), lower_dists, color=colors[i], s = 3, alpha=0.8, zorder=2)
# print mean and std of upper and lower distances
mean_upper = np.mean(all_upper_dists)
mean_lower = np.mean(all_lower_dists)
std_upper = np.std(all_upper_dists)
std_lower = np.std(all_lower_dists)
sem_upper = sem(all_upper_dists)
sem_lower = sem(all_lower_dists)
all_mean_uppers.append(mean_upper)
all_mean_lowers.append(mean_lower)
all_std_uppers.append(std_upper)
all_std_lowers.append(std_lower)
all_sem_uppers.append(sem_upper)
all_sem_lowers.append(sem_lower)
print(f'{pol} {dist_type} {amp}uA: Upper Mean: {mean_upper:.2f}, Std: {std_upper:.2f}; Lower Mean: {mean_lower:.2f}, Std: {std_lower:.2f}')
if pol == 'anodal':
amps = [5, 25, 50, 100]
else:
amps = [-5, -25, -50, -100]
ax.plot(amps, all_mean_uppers, color = colors[i], marker = 'o', markersize = 4, linewidth = 1)
ax.plot(amps, all_mean_lowers, color = colors[i], marker = 'o', markersize = 4, linewidth = 1)
ax.fill_between(amps,
np.array(all_mean_uppers) - np.array(all_sem_uppers),
np.array(all_mean_uppers) + np.array(all_sem_uppers),
color=colors[i], alpha=0.2, label='Upper Bound SEM', zorder=1)
ax.fill_between(amps,
np.array(all_mean_lowers) - np.array(all_sem_lowers),
np.array(all_mean_lowers) + np.array(all_sem_lowers),
color=colors[i], alpha=0.2, label='Lower Bound SEM', zorder=1)
# Fit polynomial regression to the data
if len(all_amps) > 1: # Ensure we have enough data points
upper_coeffs = np.polyfit(all_amps, all_upper_dists, deg=2)
lower_coeffs = np.polyfit(all_amps, all_lower_dists, deg=2)
upper_poly = np.poly1d(upper_coeffs)
lower_poly = np.poly1d(lower_coeffs)
# Generate smooth curves for plotting
amp_range = np.linspace(np.min(all_amps), np.max(all_amps), 100)
upper_fit = upper_poly(amp_range)
lower_fit = lower_poly(amp_range)
# Calculate standard error of the mean (SEM) for the fits
upper_fit_err = sem(all_upper_dists)
lower_fit_err = sem(all_lower_dists)
# Plot polynomial regression lines
#ax.plot(amp_range, upper_fit, color='darkblue', label='Upper Bound Fit', zorder=4)
#ax.plot(amp_range, lower_fit, color='darkred', label='Lower Bound Fit', zorder=4)
# Add error bands
#ax.fill_between(amp_range, upper_fit - upper_fit_err, upper_fit + upper_fit_err, color='blue', alpha=0.2, label='Upper Fit Error', zorder=1)
#ax.fill_between(amp_range, lower_fit - lower_fit_err, lower_fit + lower_fit_err, color='red', alpha=0.2, label='Lower Fit Error', zorder=1)
x = all_amps
y_upper = all_upper_dists
y_lower = all_lower_dists
# Fit linear regression for upper bounds
coeffs_linear_upper = np.polyfit(x, y_upper, 1)
p_linear_upper = np.poly1d(coeffs_linear_upper)
y_fit_linear_upper = p_linear_upper(x)
r2_linear_upper = 1 - (np.sum((y_upper - y_fit_linear_upper)**2) / np.sum((y_upper - np.mean(y_upper))**2))
# Fit quadratic regression for upper bounds
coeffs_quad_upper = np.polyfit(x, y_upper, 2)
p_quad_upper = np.poly1d(coeffs_quad_upper)
y_fit_quad_upper = p_quad_upper(x)
r2_quad_upper = 1 - (np.sum((y_upper - y_fit_quad_upper)**2) / np.sum((y_upper - np.mean(y_upper))**2))
# Fit linear regression for lower bounds
coeffs_linear_lower = np.polyfit(x, y_lower, 1)
p_linear_lower = np.poly1d(coeffs_linear_lower)
y_fit_linear_lower = p_linear_lower(x)
r2_linear_lower = 1 - (np.sum((y_lower - y_fit_linear_lower)**2) / np.sum((y_lower - np.mean(y_lower))**2))
# Fit quadratic regression for lower bounds
coeffs_quad_lower = np.polyfit(x, y_lower, 2)
p_quad_lower = np.poly1d(coeffs_quad_lower)
y_fit_quad_lower = p_quad_lower(x)
r2_quad_lower = 1 - (np.sum((y_lower - y_fit_quad_lower)**2) / np.sum((y_lower - np.mean(y_lower))**2))
# Add color-coded text annotations for R² values
# Upper bounds
#ax.text(
# 0.6, 0.90,
# f"Quad: R²: {r2_quad_upper:.2f}",
# color='blue', fontsize=5, verticalalignment='top', transform=ax.transAxes
#)
#ax.text(
# 0.6, 0.80,
# f"Lin: R²: {r2_linear_upper:.2f}",
# color='darkblue', fontsize=5, verticalalignment='top', transform=ax.transAxes
#)
#ax.text(
# 0.6, 0.10,
# f"Quad: R²: {r2_quad_lower:.2f}",
# color='red', fontsize=5, verticalalignment='top', transform=ax.transAxes)
#ax.text(
# 0.6, 0.20,
# f"Lin: R²: {r2_linear_lower:.2f}",
# color='darkred', fontsize=5, verticalalignment='top', transform=ax.transAxes)
if pol == 'anodal':
ax.set_xticks(np.abs(amps))
else:
ax.set_xticks(amps)
ax.invert_xaxis() # Invert x-axis for negative amplitudes
#ax.invert_yaxis()
ax.set_ylim([-1500, 1500])
ax.set_yticks([1500, 1000, 500, 0, -500, -1000,-1500], ['1.5', '', '', '0', '', '', '-1.5'])
ax.invert_yaxis()
# Adjust layout and save the figure
plt.tight_layout(rect=[0, 0, 1, 0.96])
path = r'C:\Users\jordan\Documents\combined_figs\evoked_potentials'
#plt.savefig(os.path.join(path, f'distances_{dist_type}_updated.png'))
#plt.savefig(os.path.join(path, f'distances_{dist_type}_updated.pdf'))
# %%
min_ch_list = []
for _, row in df_results.iterrows():
recording = row['recording']
probe = row['probe']
r = recording_map[recording]
dists = distances[recording][probe]
arg_min = np.argmin(np.abs(dists))
min_ch_list.append(arg_min)
df_results['closest_ch_ccf'] = min_ch_list
# %%
df = df_results.copy()
df = df[
(df['max_volt'] > 40) &
(df['half_max_lower_idx'] != -1) &
(df['half_max_upper_idx'] != -1) &
(df['lower_bound_idx'] != -1) &
(df['upper_bound_idx'] != -1)
].reset_index()
df['normalized_closest_ch'] = df['closest_ch_ccf'] - df['max_idx']
fig = plt.figure(figsize=(3, 3))
swarm_plot = sns.swarmplot(data=df, x='recording', y='normalized_closest_ch', hue = 'recording', palette='viridis', size = 2)
# Set labels with specified font size
plt.xlabel('Recording', fontsize=10)
plt.xticks(rotation=45)
plt.ylabel('closest ccf channel - max_response ch', fontsize=10)
plt.title('Normalized Closest Channel by Recording', fontsize=12)
#plt.savefig(os.path.join(path, 'normalized_closest_ch_by_recording.png'))
#plt.savefig(os.path.join(path, 'normalized_closest_ch_by_recording.pdf'))
from scipy.stats import f_oneway
# Make sure amplitude is treated as a group
grouped = df.groupby('recording')['normalized_closest_ch']
# Extract list of arrays (one per amplitude)
groups = [group.values for _, group in grouped]
# Run ANOVA
fval, pval = f_oneway(*groups)
print(f"One-way ANOVA: F = {fval:.4f}, p = {pval:.4e}")
from statsmodels.stats.multicomp import pairwise_tukeyhsd
tukey = pairwise_tukeyhsd(endog=df['normalized_closest_ch'],
groups=df['recording'],
alpha=0.05)
print(tukey.summary())
# %%
df = df_results.copy()
df = df[
(df['max_volt'] > 40) &
(df['half_max_lower_idx'] != -1) &
(df['half_max_upper_idx'] != -1) &
(df['lower_bound_idx'] != -1) &
(df['upper_bound_idx'] != -1) &
(df['amplitude'] != 10) &
(df['amplitude'] != -10)
].reset_index()
df['normalized_closest_ch'] = df['closest_ch_ccf'] - df['max_idx']
fig = plt.figure(figsize = (3, 3))
swarm_plot = sns.swarmplot(data=df, x='amplitude', y='normalized_closest_ch', palette='viridis', s = 2)
# Set labels with specified font size
plt.xlabel('amplitude', fontsize=10)
plt.xticks(rotation=45)
plt.ylabel('closest ccf channel - max_response ch', fontsize=10)
plt.title('Normalized Closest Channel by amplitude', fontsize=12)
from scipy.stats import f_oneway
# Make sure amplitude is treated as a group
grouped = df.groupby('amplitude')['normalized_closest_ch']
# Extract list of arrays (one per amplitude)
groups = [group.values for _, group in grouped]
# Run ANOVA
fval, pval = f_oneway(*groups)
print(f"One-way ANOVA: F = {fval:.4f}, p = {pval:.4e}")
# %%
df = df_results.copy()
df = df[
(df['max_volt'] > 40) &
(df['half_max_lower_idx'] != -1) &
(df['half_max_upper_idx'] != -1) &
(df['lower_bound_idx'] != -1) &
(df['upper_bound_idx'] != -1)
].reset_index()
df['normalized_closest_ch'] = df['closest_ch_ccf'] - df['max_idx']
fig = plt.figure(figsize = (2, 3))
df['group'] = 'all'
swarm_plot = sns.swarmplot(data=df, x='group', y='normalized_closest_ch', hue = 'group', palette='viridis', s = 2)
# Set labels with specified font size
plt.xlabel('amplitude', fontsize=10)
plt.xticks(rotation=45)
plt.ylabel('closest ccf channel - max_response ch', fontsize=10)
plt.title('Normalized Closest Channel by amplitude', fontsize=12)
# print the mean and std and sem of the normalized closest channel
mean_closest_ch = df['normalized_closest_ch'].mean()
std_closest_ch = df['normalized_closest_ch'].std()
sem_closest_ch = df['normalized_closest_ch'].sem()
print(f'Mean normalized closest channel: {mean_closest_ch:.2f}, Std: {std_closest_ch:.2f} SEM: {sem_closest_ch:.2f}')
plt.axhline(0, color='black', linestyle='--', linewidth=0.5)
path = r'C:\Users\jordan\Documents\combined_figs\evoked_potentials'
#plt.savefig(os.path.join(path, 'normalized_closest_ch_all.png'))
#plt.savefig(os.path.join(path, 'normalized_closest_ch_all.pdf'))
# %% [markdown]
# # peak voltage, half and full AUC regressions for amplitudes
# %%
df.head()
# %%
#scatters
fig, axs = plt.subplots(3, 3, figsize=(4, 4))
amplitudes = [-5, -25, -50, -100]
for col, pol in enumerate(['cathodal', 'anodal', 'bipolar']):
if pol == 'bipolar':
polarity = 'bipolar'
else:
polarity = 'monopolar'
if pol == 'anodal':
amps = np.abs(amplitudes)
else:
amps = amplitudes
df = df_results[
(df_results['max_volt'] > 40) &
(df_results['amplitude'].isin(amps)) &
(df_results['polarity'] == polarity) &
(df_results['half_max_lower_idx'] != -1) &
(df_results['half_max_upper_idx'] != -1) &
(df_results['lower_bound_idx'] != -1) &
(df_results['upper_bound_idx'] != -1)
]
variables = ['max_volt', 'half_max_AUC', 'full_AUC']
titles = ['Max Voltage by Amplitude', 'Half Max AUC by Amplitude', 'Full AUC by Amplitude']
for i, (var, title) in enumerate(zip(variables, titles)):
row = i
ax = axs[row, col]
# Scatter plot of individual points
sns.scatterplot(
data=df, x='amplitude', y=var, color = colors[col], s=6, ax=ax, legend=False, alpha=1
)
# fit first and second order regression lines style
sns.regplot(
data=df, x='amplitude', y=var, scatter=False, ci=68, order = 2, line_kws={'color': 'black', 'linewidth': 0.5}, ax=ax
)
sns.regplot(
data=df, x='amplitude', y=var, scatter=False, ci=68, order = 1, line_kws={'color': 'blue', 'linewidth': 0.5}, ax=ax
)
# print mean and std of the variables
for amp in [-5, -25, -50, -100]:
if pol == 'anodal':
amp = np.abs(amp)
amp_df = df[df['amplitude'] == amp]
mean_value = amp_df[var].mean()
std_value = amp_df[var].std()
print(f'{pol}_{amp}: {var} mean ± std: {mean_value:.2f} ± {std_value:.2f}')
#if row == 0:
#ax.set_title(title, fontsize=6)
if row == 2:
ax.set_xlabel('Amplitude (uA)', fontsize=6)
ax.set_xticks(amps, amps, fontsize = 6)
else:
ax.set_xlabel('')
ax.set_xticks(amps, ['', '', '', ''], fontsize = 6)
filtered_df = df[['amplitude', var]].dropna()
x = filtered_df['amplitude']
y = filtered_df[var]
# Fit linear regression (order=1)
coeffs_linear = np.polyfit(x, y, 1)
p_linear = np.poly1d(coeffs_linear)
y_fit_linear = p_linear(x)
r2_linear = 1 - (np.sum((y - y_fit_linear)**2) / np.sum((y - np.mean(y))**2))
# Fit quadratic regression (order=2)
coeffs_quad = np.polyfit(x, y, 2)
p_quad = np.poly1d(coeffs_quad)
y_fit_quad = p_quad(x)
r2_quad = 1 - (np.sum((y - y_fit_quad)**2) / np.sum((y - np.mean(y))**2))
ax.text(0.05, 0.95, f"Linear R² = {r2_linear:.2f}\nQuadratic R²= {r2_quad:.2f}",
transform=ax.transAxes, fontsize=4, verticalalignment='top')
if pol != 'anodal':
ax.invert_xaxis()
if var == 'max_volt':
ax.set_ylim([0, 2000])
if col == 0:
ax.set_ylabel('Max Voltage (uV)', fontsize=6)
ax.set_yticks([0, 1000, 2000], [0, 1, 2], fontsize = 6)
else:
ax.set_ylabel('')
ax.set_yticks([0, 1000, 2000], ['', '', ''], fontsize = 6)
elif var == 'half_max_AUC':
ax.set_ylim([0, 80000])
if col == 0:
ax.set_ylabel('Half Max AUC (mV*ms)', fontsize=6)
ax.set_yticks([0, 40000, 80000], [0, 4, 8], fontsize = 6)
else:
ax.set_ylabel('')
ax.set_yticks([0, 40000, 80000], ['', '', ''], fontsize = 6)
elif var == 'full_AUC':
ax.set_ylim([0, 120000])
if col == 0:
ax.set_yticks([0, 60000, 120000], [0, 6, 12], fontsize = 6)
ax.set_ylabel('Full AUC (mV*ms)', fontsize=6)
else:
ax.set_yticks([0, 60000, 120000], ['', '', ''], fontsize = 6)
ax.set_ylabel('')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
# Adjust layout and save the figure
#plt.suptitle()
plt.tight_layout()
#plt.savefig(os.path.join(path, f'max_auc_scatter_with_means__poly_updated.png'))
#plt.savefig(os.path.join(path, f'max_auc_scatter_with_means__poly_updated.pdf'))
# %% [markdown]
# # comparing LFP and AP
# %%
r = jlh31
probe = 'probeB'
amp = -100
channels = 300
stim_times = choose_stim_parameter(r.trials, amp=amp, pulse_number=1, pulse_duration=100, polarity='monopolar')[0:75]
fig, axes = plt.subplots(1,2)
## AP plot
pre = -0.3
post = 5
data = r.raw.get_chunk(probe, stim_times, pre=pre + 3, post=post + 3, chs=np.arange(channels))
aligned_data = align_data(data, pre, post, channels, threshold=500, median_subtraction=False)
dists = distances[r.mouse][probe]
# Use the axs object for the current amp
sub_chs = np.arange(len(dists), len(dists) + 20)
aligned_sub_data = aligned_data - np.median(aligned_data[:,:, sub_chs], axis=2)[:,:, np.newaxis]
ax = axes[0]
raw_heatmap(aligned_sub_data, pre=pre, post=post, dists=dists, vmin=-300, vmax=300,
save=False, title=f'{amp} uA', ax=ax)
# LFP
pre = 0
post = 100
data = r.raw.get_chunk(probe, stim_times, pre=pre, post=post, band = 'lfp', chs=np.arange(channels))
aligned_data = data
dists = distances[r.mouse][probe]
# Use the axs object for the current amp
sub_chs = np.arange(len(dists), len(dists) + 40)
aligned_sub_data = aligned_data - np.median(aligned_data[:,:, sub_chs], axis=2)[:,:, np.newaxis]
ax = axes[1]
raw_heatmap(aligned_sub_data, pre=pre, post=post, dists=dists, vmin=-800, vmax=800,
save=False, title=f'{amp} uA', ax=ax)
plt.tight_layout()
# %% [markdown]
# # volumes
# %%
from scipy.spatial import ConvexHull
import numpy as np
def calculate_volume(points):
"""
Calculate the volume of a 3D shape defined by 6 points using a convex hull.
Parameters:
points (ndarray): A (6, 3) numpy array where each row represents a 3D point (x, y, z).
Returns:
float: The volume of the convex hull enclosing the points.
"""
# Check if the input is valid
if points.shape != (6, 3):
raise ValueError("Input must be a (6, 3) array of 3D points.")
# Compute the convex hull
hull = ConvexHull(points) # points should be (6,3) array
volume = hull.volume
# Return the volume of the convex hull
return hull.volume
# %%
#3D stuff... this requires measurements for all 3 probes...
polarities = ['bipolar', 'monopolar']
probes = ['probeA', 'probeB', 'probeC']
amps = [-5, -25, -50, -100, 5, 25, 50, 100]
volume_dict = {}
points_dict = {}
points_dict_rec = {}
rec_dict = {}
for polarity in polarities:
for i, amp in enumerate(amps):
# Filter data for the current probe and amplitude
df = df_results[
(df_results['max_volt'] > 40) &
(df_results['amplitude'] == amp) &
(df_results['polarity'] == polarity) &
(df_results['lower_bound_idx'] != -1) &
(df_results['upper_bound_idx'] != -1)
]
#print(df['recording'].unique())
for record in df['recording'].unique():
probes = df[df['recording'] == record]['probe'].unique()
if len(probes) == 3:
key = f'{polarity}_{amp}'
r = recording_map[record]
a_df = df[(df['recording'] == record) & (df['probe'] == 'probeA')]
a_upper_ch = a_df['upper_bound_idx'].values[0]
a_lower_ch = a_df['lower_bound_idx'].values[0]
try:
a_up_ccf = r.probe_coords['probeA'][a_upper_ch]
except:
a_up_cff = r.probe_coords['probeA'][-1]
a_low_ccf = r.probe_coords['probeA'][a_lower_ch]
b_df = df[(df['recording'] == record) & (df['probe'] == 'probeB')]
b_upper_ch = b_df['upper_bound_idx'].values[0]
b_lower_ch = b_df['lower_bound_idx'].values[0]
try:
b_up_ccf = r.probe_coords['probeB'][b_upper_ch]
except:
b_up_ccf = r.probe_coords['probeB'][-1]
b_low_ccf = r.probe_coords['probeB'][b_lower_ch]
c_df = df[(df['recording'] == record) & (df['probe'] == 'probeC')]
c_upper_ch = c_df['upper_bound_idx'].values[0]
c_lower_ch = c_df['lower_bound_idx'].values[0]
try:
c_up_ccf = r.probe_coords['probeC'][c_upper_ch]
except:
c_up_ccf = r.probe_coords['probeC'][-1]
c_low_ccf = r.probe_coords['probeC'][c_lower_ch]
points = np.array([a_up_ccf, a_low_ccf, b_up_ccf, b_low_ccf, c_up_ccf, c_low_ccf])
hull = ConvexHull(points)
volume = hull.volume
if key not in volume_dict:
volume_dict[key] = []
volume_dict[key].append(volume)
if key not in points_dict:
points_dict[key] = []
points_dict[key].append(points)
key_rec = f'{record}_{polarity}_{amp}'
if key_rec not in points_dict_rec:
points_dict_rec[key_rec] = []
points_dict_rec[key_rec].append(points)
if key not in rec_dict:
rec_dict[key] = []
rec_dict[key].append(record)
# %%
volume_dict_mm3 = {key: np.array(values) / 1e9 for key, values in volume_dict.items()}
mc_25 = volume_dict_mm3['monopolar_-25']
mc_50 = volume_dict_mm3['monopolar_-50']
mc_100 = volume_dict_mm3['monopolar_-100']
mc_x = [-25, -50, -100]
mc_y = [np.mean(mc_25), np.mean(mc_50), np.mean(mc_100)]
mc_err = [sem(mc_25), sem(mc_50), sem(mc_100)]
# print mean +/- std of monopolar volumes
print(f'Cathodal Mean Volumes: -25uA: {np.mean(mc_25):.2f}+/-{np.std(mc_25):.2f} mm³, -50uA: {np.mean(mc_50):.2f}+/-{np.std(mc_50):.2f}mm³, -100uA: {np.mean(mc_100):.2f}+/-{np.std(mc_100):.2f} mm³')
fig, ax = plt.subplots(1, 1, figsize=(2, 3.5))
# Line and error fill for monopolar (cathodal)
ax.plot(mc_x, mc_y, label='Cathodal', color=cathodal_color, marker='o')
ax.fill_between(mc_x, np.array(mc_y) - np.array(mc_err), np.array(mc_y) + np.array(mc_err),
color=cathodal_color, alpha=0.2)
# Plot individual values with jitter
jitter_strength = 1.1
for x, values in zip(mc_x, [mc_25, mc_50, mc_100]):
jittered_x = np.random.normal(x, scale=jitter_strength, size=len(values))
ax.scatter(jittered_x, values, color=cathodal_color, s=5, alpha=0.7, zorder=3)
# Labels and formatting
ax.set_title('Mean Volume')
ax.set_xlabel('Amplitude (uA)')
ax.set_ylabel('Volume (mm³)')
ax.invert_xaxis()
ax.set_xticks([-25, -50, -100])
ax.set_ylim(0, 1)
ax.set_yticks([0, 0.5, 1])
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.tight_layout()
path = r'C:\Users\jordan\Documents\combined_figs\evoked_potentials'
#plt.savefig(os.path.join(path, 'monopolar_EPvolume.png'))
#plt.savefig(os.path.join(path, 'monopolar_EPvolume.pdf'))
# %%
from scipy.stats import linregress
import numpy as np
# Stack all individual volumes
all_volumes = np.concatenate([mc_25, mc_50, mc_100])
all_amplitudes = np.concatenate([
np.full_like(mc_25, -25),
np.full_like(mc_50, -50),
np.full_like(mc_100, -100)
])
# Run linear regression
slope, intercept, r_value, p_value, std_err = linregress(np.abs(all_amplitudes), all_volumes)
# Print regression results
print(f"Linear Regression: Volume = {slope:.4f} * Amplitude + {intercept:.4f}")
print(f"R² = {r_value**2:.4f}, p = {p_value:.4e}")
# %%
9.1036e-02
# %%
# Quadratic fit
quad_coeffs = np.polyfit(np.abs(all_amplitudes), all_volumes, deg=2)
quad_model = np.poly1d(quad_coeffs)
# Predicted values
y_quad_fit = quad_model(np.abs(all_amplitudes))
# R² for quadratic
ss_res_quad = np.sum((all_volumes - y_quad_fit) ** 2)
ss_tot = np.sum((all_volumes - np.mean(all_volumes)) ** 2)
r2_quad = 1 - (ss_res_quad / ss_tot)
print(f"Quadratic Fit: Volume = {quad_coeffs[0]:.4e} * A² + {quad_coeffs[1]:.4e} * A + {quad_coeffs[2]:.4e}")
print(f"R² (Quadratic) = {r2_quad:.4f}")
# %%
from scipy.stats import sem
# Convert µm³ to mm³ (1 mm³ = 10^9 µm³)
volume_dict_mm3 = {key: np.array(values) / 1e9 for key, values in volume_dict.items()}
bc_25 = volume_dict_mm3['bipolar_-25']
bc_50 = volume_dict_mm3['bipolar_-50']
bc_100 = volume_dict_mm3['bipolar_-100']
mc_25 = volume_dict_mm3['monopolar_-25']
mc_50 = volume_dict_mm3['monopolar_-50']
mc_100 = volume_dict_mm3['monopolar_-100']
ma_25 = volume_dict_mm3['monopolar_25']
ma_50 = volume_dict_mm3['monopolar_50']
ma_100 = volume_dict_mm3['monopolar_100']
bc_x = [-25, -50, -100]
bc_y = [np.mean(bc_25), np.mean(bc_50), np.mean(bc_100)]
bc_err = [sem(bc_25), sem(bc_50), sem(bc_100)]
mc_x = [-25, -50, -100]
mc_y = [np.mean(mc_25), np.mean(mc_50), np.mean(mc_100)]
mc_err = [sem(mc_25), sem(mc_50), sem(mc_100)]
ma_x = [25, 50, 100]
ma_y = [np.mean(ma_25), np.mean(ma_50), np.mean(ma_100)]
ma_err = [sem(ma_25), sem(ma_50), sem(ma_100)]
print(f'Cathodal Mean Volumes: -25uA: {np.mean(mc_25):.2f}+/-{np.std(mc_25):.2f} mm³, -50uA: {np.mean(mc_50):.2f}+/-{np.std(mc_50):.2f}mm³, -100uA: {np.mean(mc_100):.2f}+/-{np.std(mc_100):.2f} mm³')
print(f'Anodal Mean Volumes: 25uA: {np.mean(ma_25):.2f}+/-{np.std(ma_25):.2f} mm³, 50uA: {np.mean(ma_50):.2f}+/-{np.std(ma_50):.2f}mm³, 100uA: {np.mean(ma_100):.2f}+/-{np.std(ma_100):.2f} mm³')
print(f'Bipolar Mean Volumes: -25uA: {np.mean(bc_25):.2f}+/-{np.std(bc_25):.2f} mm³, -50uA: {np.mean(bc_50):.2f}+/-{np.std(bc_50):.2f}mm³, -100uA: {np.mean(bc_100):.2f}+/-{np.std(bc_100):.2f} mm³')
fig, ax = plt.subplots(1, 1, figsize=(3.5, 3.5))
# Bipolar
ax.plot(bc_x, bc_y, label='Bipolar', color=bipolar_color, marker='o')
ax.fill_between(bc_x, np.array(bc_y) - np.array(bc_err), np.array(bc_y) + np.array(bc_err), color=bipolar_color, alpha=0.2)
# Cathodal
ax.plot(mc_x, mc_y, label='Cathodal', color=cathodal_color, marker='o')
ax.fill_between(mc_x, np.array(mc_y) - np.array(mc_err), np.array(mc_y) + np.array(mc_err), color=cathodal_color, alpha=0.2)
# Anodal
ax.plot(ma_x, ma_y, label='Anodal', color=anodal_color, marker='o')
ax.fill_between(ma_x, np.array(ma_y) - np.array(ma_err), np.array(ma_y) + np.array(ma_err), color=anodal_color, alpha=0.2)
# Labels and legend
ax.set_title('Mean Volume')
ax.set_xlabel('Amplitude (uA)')
ax.set_ylabel('Volume (mm³)')
#ax.legend(loc='upper right')
ax.set_xticks([-100, -50, -25, 25, 50, 100])
ax.set_ylim(0, 1)
ax.set_yticks([0, 0.5, 1])
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.tight_layout()
#plt.savefig(os.path.join(path, 'volume_by_amplitude.png'))
#plt.savefig(os.path.join(path, 'volume_by_amplitude.pdf'))
# %%
import statsmodels.api as sm
from statsmodels.formula.api import ols
model = ols("volume ~ C(polarity) * amplitude", data=df_vol).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)
# %%

figure2.ipynb at commit c4aaa60, under GPL-3.0 · at the source

OverviewAuthors: Jordan L Hickman1,2,3, Grant Hughes1,3, Eashan Sahai1, Moriah Miles1, Daniel J Denman1,4,5
ORCID iDs: Eashan Sahai, Daniel J Denman
Department of Biophysics and Physiology, University of Colorado Anschutz Medical Campus, Aurora, CO, USA
Medical Scientist Training Program, University of Colorado Anschutz Medical Campus, Aurora, CO, USA
Neuroscience Graduate Program, University of Colorado Anschutz Medical Campus, Aurora, CO, USA
Department of Ophthalmology, University of Colorado Anschutz Medical Campus, Aurora, CO, USA
Lead contact
Institutions: University of Colorado Anschutz (United States)
Journal: Cell reports, volume 45, issue 6, article 117420
Dates: published online 28 May 2026; in print 23 June 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1016/j.celrep.2026.117420 · PMID 42207642 · PMCID PMC13404931 · OpenAlex W4415660446
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.10.28.685195
Status: code verified
Categories: extracellular electrophysiology (units, LFP) (modality), mouse (organism), systems (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Single-unit activity, calcium imaging, Physiology & signal measures, Connectivity
Keywords: Visual cortex, Electrical stimulation, Neural Dynamics, Neuromodulation, Intracortical Microstimulation, Icms, Neuropixels, Cp: Neuroscience, Perceptual Detection
MeSH: Electric Stimulation*, Neocortex*, Neurons*, Action Potentials, Animals, Male, Mice, Mice, Inbred C57BL, Pyramidal Cells, Visual Cortex (* major topic)
Topic: Neuroscience and Neural Engineering (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex
Funding: National Institute of Neurological Disorders and Stroke (R01NS120850); National Eye Institute (R00EY028612); NINDS NIH HHS (R01 NS120850); NEI NIH HHS (R00 EY028612); National Institutes of Health
Citations: not cited yet (Europe PMC); 122 references in the paper
Research resources: electrophysiology data in NWB format RRID:SCR_015242
AbstractThe abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

Its files are read in the Code ↔ Paper reader above, with 7 matches between paragraphs and lines of code.

denmanlab/am4100_code

License: MIT

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: b1587c33836b7c8626e5b3816ff4d4753b978043, 22 May 2024

Languages: MATLAB (30), Python (2), Jupyter (1)

Size: 53 files, 33 scripts

Software Heritage: archived

Found in: the text, “Footnotes”

Holds: license file, 1 notebook

Not found: README, CITATION.cff, environment file, tests, continuous integration, documentation

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

34 files

am4100.py, Python, 153 lines
matlab_from_amsystems/ComConstants.m, MATLAB, 116 lines
matlab_from_amsystems/DrawBiphasic.m, MATLAB, 392 lines
matlab_from_amsystems/DrawMonophasic.m, MATLAB, 215 lines
matlab_from_amsystems/DrawRamp.m, MATLAB, 313 lines
matlab_from_amsystems/EXAMPLE_1_USB.m, MATLAB, 115 lines
matlab_from_amsystems/EXAMPLE_2_Ethernet.m, MATLAB, 106 lines
matlab_from_amsystems/IDnum.m, MATLAB, 158 lines
matlab_from_amsystems/LoadLibrary.m, MATLAB, 45 lines
matlab_from_amsystems/LoadWindow.m, MATLAB, 218 lines
matlab_from_amsystems/Plotit.m, MATLAB, 396 lines
matlab_from_amsystems/RetEamp1.m, MATLAB, 38 lines
matlab_from_amsystems/RetEamp3.m, MATLAB, 38 lines
matlab_from_amsystems/SetCom.m, MATLAB, 66 lines
matlab_from_amsystems/UpdateEvents.m, MATLAB, 35 lines
matlab_from_amsystems/VorIvalue.m, MATLAB, 59 lines
matlab_from_amsystems/ams4100_GUI.m, MATLAB, 534 lines
matlab_from_amsystems/ams4100_hClass.m, MATLAB, 1,553 lines
matlab_from_amsystems/checkTimes.m, MATLAB, 46 lines
matlab_from_amsystems/controlsGUI.m, MATLAB, 51 lines
matlab_from_amsystems/defaultData.m, MATLAB, 70 lines
matlab_from_amsystems/getOffsetAmps.m, MATLAB, 81 lines
matlab_from_amsystems/getTimeValues.m, MATLAB, 85 lines
matlab_from_amsystems/getValueNames.m, MATLAB, 71 lines
matlab_from_amsystems/nidaq_test.m, MATLAB, 13 lines
matlab_from_amsystems/processUserInput.m, MATLAB, 421 lines
matlab_from_amsystems/setTimeAMS.m, MATLAB, 45 lines
matlab_from_amsystems/timeNum.m, MATLAB, 52 lines
matlab_from_amsystems/trainNum.m, MATLAB, 38 lines
matlab_from_amsystems/untitled2.m, MATLAB, 2 lines
matlab_from_amsystems/untitled3.m, MATLAB, 1 line
read_serial.py, Python, 26 lines
workshopping_am4100class.ipynb, Jupyter, 190 lines
LICENSE, License, 21 lines

At the source: github.com/denmanlab/am4100_code

denmanlab/mouse_behavior

License: none: the authors keep all their rights

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: 59fb7046918c4f9cb361b8df584a525af5a0cc2e, 1 August 2025

Languages: Python (17), Jupyter (1)

Size: 60 files, 18 scripts

Software Heritage: not archived

Found in: the text, “Software and hardware implementation:”

Holds: README, 1 notebook

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: NumPy (10 files), pandas (10 files), Matplotlib (7 files), imageio (1 file), OpenCV (1 file), Pillow (1 file), seaborn (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

19 files

archive/arduino_controller.py, Python, 96 lines
archive/camera.py, Python, 304 lines
archive/camera2.py, Python, 103 lines
archive/mouse_detection.py, Python, 559 lines
archive/mouse_detection2.py, Python, 673 lines
archive/mouse_detection_contrast.py, Python, 937 lines
archive/mouse_detection_shaping.py, Python, 924 lines
archive/mouse_detection_timeout.py, Python, 744 lines
archive/solenoid_controller.py, Python, 39 lines
archive/sync_to_backups.py, Python, 187 lines
arduino_controller_jlh.py, Python, 55 lines
custom_timer.py, Python, 11 lines
dev.py, Python, 533 lines, 1 match
params.py, Python, 177 lines, 1 match
plotter.py, Python, 710 lines, 1 match
stimuli.py, Python, 29 lines
testing.py, Python, 28 lines
testing_backup_notebook.ipynb, Jupyter, 108 lines
README.md, Text, 1 line

At the source: github.com/denmanlab/mouse_behavior

denmanlab/estim_populations

License: GPL-3.0

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: c4aaa608314c8e8e91aee422baa69cd1b443a682, 29 October 2025

Languages: Python (12), Jupyter (7)

Size: 97 files, 19 scripts

Software Heritage: not archived

Found in: the text, “Dataset”

Holds: README, license file, environment (requirements.txt), 7 notebooks

Not found: CITATION.cff, tests, continuous integration, documentation

Tools: Matplotlib (3 files), NumPy (3 files), Neurodata Without Borders (PyNWB, MatNWB) (3 files), Open Ephys analysis tools (3 files), pandas (3 files), SciPy (3 files), seaborn (3 files), statsmodels (3 files), scikit-image (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

5 files

figure1.ipynb, Jupyter, 825 lines
figure2.ipynb, Jupyter, 2,260 lines, 2 matches
figure3.ipynb, Jupyter, 2,826 lines, 2 matches
repository limit reached (2,000 files or 30 MB): the rest is at the source (16 files)
LICENSE, License, 674 lines
README.md, Text, 77 lines

At the source: github.com/denmanlab/estim_populations

marcomusy/vedo

License: MIT

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: 88f8b2f9c5c6267b5dceda0e4ab0e15b3130910a, 4 August 2026

Languages: Python (559), Jupyter (11), Shell (10), JavaScript (2)

Size: 642 files, 582 scripts

Software Heritage: archived

Found in: the text, “3D plotting in anatomical brain space”

Holds: README, license file, environment (pyproject.toml), tests, continuous integration, documentation, 11 notebooks

Not found: CITATION.cff

Tools: NumPy (221 files), SciPy (29 files), Pillow (12 files), Matplotlib (8 files), pandas (4 files), scikit-image (2 files), imageio (1 file), JAX (1 file), napari (1 file), NetworkX (1 file), OpenCV (1 file), scikit-learn (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

584 files

cli_entry.py, Python, 17 lines
docs/assets/search-expand.js, JavaScript, 10 lines
docs/examples_db.js, JavaScript, 2,821 lines
docs/logos/embl_logo.py, Python, 44 lines
docs/logos/lab_logo_maker.py, Python, 21 lines
docs/logos/logo_vedo_simple.py, Python, 8 lines
examples/advanced/capping_mesh.py, Python, 56 lines
examples/advanced/contours2mesh.py, Python, 22 lines
examples/advanced/convex_hull.py, Python, 13 lines
examples/advanced/cut_and_cap.py, Python, 16 lines
examples/advanced/cut_with_mesh1.py, Python, 19 lines
examples/advanced/cut_with_points1.py, Python, 24 lines
examples/advanced/cut_with_points2.py, Python, 26 lines
examples/advanced/diffuse_data.py, Python, 67 lines
examples/advanced/fitline.py, Python, 28 lines
examples/advanced/fitplanes.py, Python, 25 lines
examples/advanced/fitspheres1.py, Python, 30 lines
examples/advanced/fitspheres2.py, Python, 31 lines
examples/advanced/geodesic_curve.py, Python, 17 lines
examples/advanced/geological_model.py, Python, 156 lines
examples/advanced/gyroid.py, Python, 22 lines
examples/advanced/interpolate_field.py, Python, 66 lines
examples/advanced/interpolate_scalar1.py, Python, 21 lines
examples/advanced/interpolate_scalar2.py, Python, 33 lines
examples/advanced/interpolate_scalar3.py, Python, 17 lines
examples/advanced/interpolate_scalar4.py, Python, 34 lines
examples/advanced/interpolate_scalar5.py, Python, 49 lines
examples/advanced/line2mesh_quads.py, Python, 31 lines
examples/advanced/line2mesh_tri.py, Python, 26 lines
examples/advanced/measure_curvature1.py, Python, 45 lines
examples/advanced/measure_curvature2.py, Python, 80 lines
examples/advanced/mesh_smoother1.py, Python, 21 lines
examples/advanced/mesh_smoother2.py, Python, 21 lines
examples/advanced/meshquality.py, Python, 33 lines
examples/advanced/moving_least_squares1D.py, Python, 38 lines
examples/advanced/moving_least_squares2D.py, Python, 52 lines
examples/advanced/multi_viewer2.py, Python, 64 lines
examples/advanced/recosurface.py, Python, 43 lines
examples/advanced/run_all.sh, Shell, 16 lines
examples/advanced/skeletonize.py, Python, 17 lines
examples/advanced/spline_draw1.py, Python, 16 lines
examples/advanced/spline_draw2.py, Python, 18 lines
examples/advanced/splitmesh.py, Python, 11 lines
examples/advanced/timer_callback0.py, Python, 22 lines
examples/advanced/timer_callback1.py, Python, 51 lines
examples/advanced/timer_callback2.py, Python, 63 lines
examples/advanced/timer_callback3.py, Python, 55 lines
examples/advanced/voronoi2.py, Python, 17 lines
examples/advanced/warp1.py, Python, 31 lines
examples/advanced/warp2.py, Python, 30 lines
examples/advanced/warp3.py, Python, 103 lines
examples/advanced/warp4a.py, Python, 160 lines
examples/advanced/warp4b.py, Python, 120 lines
examples/advanced/warp5.py, Python, 141 lines
examples/advanced/warp6.py, Python, 33 lines
examples/animation/airplane1.py, Python, 21 lines
examples/animation/airplane2.py, Python, 27 lines
examples/animation/aizawa_attractor.py, Python, 67 lines
examples/animation/aspring1.py, Python, 40 lines
examples/animation/aspring2_player.py, Python, 49 lines
examples/animation/brownian2d.py, Python, 130 lines
examples/animation/doubleslit.py, Python, 59 lines
examples/animation/drag_chain.py, Python, 32 lines
examples/animation/fourier_epicycles.py, Python, 77 lines
examples/animation/gas.py, Python, 134 lines
examples/animation/grayscott.py, Python, 83 lines
examples/animation/gyroscope1.py, Python, 62 lines
examples/animation/koch_fractal.py, Python, 34 lines
examples/animation/lorenz.py, Python, 45 lines
examples/animation/mag_field1.py, Python, 69 lines
examples/animation/multiple_pendulum.py, Python, 126 lines
examples/animation/optics_base.py, Python, 324 lines
examples/animation/optics_main1.py, Python, 108 lines
examples/animation/optics_main2.py, Python, 59 lines
examples/animation/optics_main3.py, Python, 34 lines
examples/animation/particle_simulator.py, Python, 123 lines
examples/animation/pendulum_3d.py, Python, 49 lines
examples/animation/pendulum_ode.py, Python, 66 lines
examples/animation/run_all.sh, Shell, 14 lines
examples/animation/self_org_maps2d.py, Python, 78 lines
examples/animation/spline_ease.py, Python, 27 lines
examples/animation/springs_fem.py, Python, 81 lines
examples/animation/trail.py, Python, 23 lines
examples/animation/tunnelling1.py, Python, 52 lines
examples/animation/tunnelling2.py, Python, 79 lines
examples/animation/value_iteration.py, Python, 123 lines
examples/animation/volterra.py, Python, 72 lines
examples/animation/wave_equation1d.py, Python, 117 lines
examples/animation/wave_equation2d.py, Python, 65 lines
examples/basic/align1.py, Python, 21 lines
examples/basic/align2.py, Python, 29 lines
examples/basic/align3.py, Python, 26 lines
examples/basic/align4.py, Python, 26 lines
examples/basic/align5.py, Python, 31 lines
examples/basic/align6.py, Python, 21 lines
examples/basic/background_image.py, Python, 24 lines
examples/basic/boolean.py, Python, 29 lines
examples/basic/boundaries.py, Python, 18 lines
examples/basic/buildmesh.py, Python, 19 lines
examples/basic/buttons1.py, Python, 31 lines
examples/basic/buttons2.py, Python, 52 lines
examples/basic/buttons3.py, Python, 44 lines
examples/basic/cartoony.py, Python, 37 lines
examples/basic/cells_within_bounds.py, Python, 30 lines
examples/basic/clustering.py, Python, 24 lines
examples/basic/color_mesh_cells1.py, Python, 25 lines
examples/basic/color_mesh_cells2.py, Python, 27 lines
examples/basic/colorcubes.py, Python, 69 lines
examples/basic/colorlines.py, Python, 39 lines
examples/basic/colormap_list.py, Python, 43 lines
examples/basic/colormaps.py, Python, 35 lines
examples/basic/connected_vtx.py, Python, 25 lines
examples/basic/cut_freehand.py, Python, 62 lines
examples/basic/cut_interactive.py, Python, 61 lines
examples/basic/delaunay2d.py, Python, 20 lines
examples/basic/delete_mesh_pts.py, Python, 30 lines
examples/basic/distance2mesh.py, Python, 19 lines
examples/basic/extrude1.py, Python, 17 lines
examples/basic/extrude2.py, Python, 25 lines
examples/basic/fillholes.py, Python, 14 lines
examples/basic/flatarrow.py, Python, 15 lines
examples/basic/glyphs1.py, Python, 45 lines
examples/basic/glyphs2.py, Python, 35 lines
examples/basic/glyphs3.py, Python, 21 lines
examples/basic/hover_legend.py, Python, 29 lines
examples/basic/input_box.py, Python, 39 lines
examples/basic/interaction_modes1.py, Python, 24 lines
examples/basic/interaction_modes2.py, Python, 31 lines
examples/basic/interaction_modes3.py, Python, 24 lines
examples/basic/interaction_modes4.py, Python, 40 lines
examples/basic/keypress.py, Python, 27 lines
examples/basic/largestregion.py, Python, 17 lines
examples/basic/legendbox.py, Python, 29 lines
examples/basic/light_sources.py, Python, 33 lines
examples/basic/lightings.py, Python, 15 lines
examples/basic/lin_interpolate.py, Python, 20 lines
examples/basic/manyspheres.py, Python, 21 lines
examples/basic/mesh_alphas.py, Python, 20 lines
examples/basic/mesh_coloring.py, Python, 36 lines
examples/basic/mesh_custom.py, Python, 27 lines
examples/basic/mesh_lut.py, Python, 48 lines
examples/basic/mesh_map2cell.py, Python, 22 lines
examples/basic/mesh_modify.py, Python, 15 lines
examples/basic/mesh_sharemap.py, Python, 16 lines
examples/basic/mesh_threshold.py, Python, 17 lines
examples/basic/mirror.py, Python, 14 lines
examples/basic/mouseclick1.py, Python, 28 lines
examples/basic/mouseclick2.py, Python, 33 lines
examples/basic/mousehighlight.py, Python, 30 lines
examples/basic/mousehover0.py, Python, 26 lines
examples/basic/mousehover1.py, Python, 43 lines
examples/basic/mousehover2.py, Python, 28 lines
examples/basic/mousehover3.py, Python, 38 lines
examples/basic/multirenderers.py, Python, 26 lines
examples/basic/multiwindows1.py, Python, 52 lines
examples/basic/multiwindows2.py, Python, 34 lines
examples/basic/pca_ellipse.py, Python, 21 lines
examples/basic/pca_ellipsoid.py, Python, 49 lines
examples/basic/record_play.py, Python, 13 lines
examples/basic/ribbon.py, Python, 19 lines
examples/basic/rotate_image.py, Python, 20 lines
examples/basic/run_all.sh, Shell, 10 lines
examples/basic/scalarbars.py, Python, 35 lines
examples/basic/shadow1.py, Python, 13 lines
examples/basic/shadow2.py, Python, 20 lines
examples/basic/shadow3.py, Python, 17 lines
examples/basic/shrink.py, Python, 11 lines
examples/basic/silhouette1.py, Python, 25 lines
examples/basic/silhouette2.py, Python, 50 lines
examples/basic/silhouette3.py, Python, 12 lines
examples/basic/skybox.py, Python, 15 lines
examples/basic/slider_browser.py, Python, 35 lines
examples/basic/sliders1.py, Python, 39 lines
examples/basic/sliders2.py, Python, 66 lines
examples/basic/sliders3d.py, Python, 30 lines
examples/basic/sliders_hsv.py, Python, 95 lines
examples/basic/sliders_range.py, Python, 50 lines
examples/basic/specular.py, Python, 33 lines
examples/basic/spline_tool.py, Python, 35 lines
examples/basic/ssao.py, Python, 23 lines
examples/basic/surf_intersect.py, Python, 17 lines
examples/basic/texture_coords.py, Python, 35 lines
examples/basic/texturecubes.py, Python, 21 lines
examples/basic/tube_radii.py, Python, 28 lines
examples/basic/voronoi1.py, Python, 19 lines
examples/extras/chemistry1.py, Python, 49 lines
examples/extras/chemistry2.py, Python, 35 lines
examples/extras/clone2d.py, Python, 23 lines
examples/extras/ellipt_fourier_desc.py, Python, 28 lines
examples/extras/export_numpy.py, Python, 22 lines
examples/extras/export_threejs.py, Python, 39 lines
examples/extras/export_x3d.py, Python, 22 lines
examples/extras/fast_simpl.py, Python, 40 lines
examples/extras/flag_labels1.py, Python, 30 lines
examples/extras/flag_labels2.py, Python, 17 lines
examples/extras/icon.py, Python, 17 lines
examples/extras/image_editor.py, Python, 19 lines
examples/extras/image_false_colors.py, Python, 27 lines
examples/extras/image_fft.py, Python, 24 lines
examples/extras/image_probe.py, Python, 37 lines
examples/extras/image_rgba.py, Python, 31 lines
examples/extras/image_to_mesh.py, Python, 35 lines
examples/extras/iminuit1.py, Python, 34 lines
examples/extras/iminuit2.py, Python, 48 lines
examples/extras/inset.py, Python, 49 lines
examples/extras/madcad1.py, Python, 79 lines
examples/extras/madcad2.py, Python, 34 lines
examples/extras/magic-class1.py, Python, 46 lines
examples/extras/make_video.py, Python, 54 lines
examples/extras/meshio_read.py, Python, 15 lines
examples/extras/meshlib1.py, Python, 143 lines
examples/extras/morphomatics_riemann.py, Python, 51 lines
examples/extras/morphomatics_tube.py, Python, 45 lines
examples/extras/napari1.py, Python, 28 lines
examples/extras/nelder-mead.py, Python, 29 lines
examples/extras/nevergrad_opt.py, Python, 43 lines
examples/extras/printc.py, Python, 36 lines
examples/extras/pygeodesic1.py, Python, 18 lines
examples/extras/pygmsh_cut.py, Python, 36 lines
examples/extras/pymeshlab1.py, Python, 45 lines
examples/extras/pymeshlab2.py, Python, 40 lines
examples/extras/pysr_regression.py, Python, 79 lines
examples/extras/qt_cutter.py, Python, 54 lines
examples/extras/qt_tabs.py, Python, 119 lines
examples/extras/qt_window1.py, Python, 57 lines
examples/extras/qt_window2.py, Python, 57 lines
examples/extras/qt_window3.py, Python, 57 lines
examples/extras/quaternions.py, Python, 105 lines
examples/extras/remesh_ACVD.py, Python, 46 lines
examples/extras/remesh_meshfix.py, Python, 48 lines
examples/extras/run_all.sh, Shell, 14 lines
examples/extras/spherical_harmonics1.py, Python, 91 lines
examples/extras/tensor_grid1.py, Python, 13 lines
examples/extras/tetgen1.py, Python, 60 lines
examples/extras/trame_ex1.py, Python, 51 lines
examples/extras/trame_ex2.py, Python, 54 lines
examples/extras/trame_ex3.py, Python, 97 lines
examples/extras/trimesh_nearest.py, Python, 35 lines
examples/extras/trimesh_ray.py, Python, 41 lines
examples/extras/trimesh_section.py, Python, 45 lines
examples/extras/trimesh_shortest.py, Python, 43 lines
examples/extras/wx_window.py, Python, 49 lines
examples/notebooks/align1.ipynb, Jupyter, 29 lines
examples/notebooks/distance2mesh.ipynb, Jupyter, 24 lines
examples/notebooks/interpolate_volume.ipynb, Jupyter, 36 lines
examples/notebooks/legosurface.ipynb, Jupyter, 23 lines
examples/notebooks/manipulate_camera.ipynb, Jupyter, 31 lines
examples/notebooks/numpy2volume.ipynb, Jupyter, 28 lines
examples/notebooks/pca.ipynb, Jupyter, 29 lines
examples/notebooks/shrink.ipynb, Jupyter, 15 lines
examples/notebooks/slider2d.ipynb, Jupyter, 62 lines
examples/notebooks/sphere.ipynb, Jupyter, 19 lines
examples/notebooks/test_types.ipynb, Jupyter, 67 lines
examples/pyplot/andrews_cluster.py, Python, 44 lines
examples/pyplot/anim_lines.py, Python, 67 lines
examples/pyplot/caption.py, Python, 25 lines
examples/pyplot/custom_axes1.py, Python, 54 lines
examples/pyplot/custom_axes2.py, Python, 39 lines
examples/pyplot/custom_axes3.py, Python, 30 lines
examples/pyplot/custom_axes4.py, Python, 43 lines
examples/pyplot/earthquake_browser.py, Python, 49 lines
examples/pyplot/embed_matplotlib1.py, Python, 19 lines
examples/pyplot/embed_matplotlib2.py, Python, 72 lines
examples/pyplot/explore5d.py, Python, 60 lines
examples/pyplot/fill_gap.py, Python, 36 lines
examples/pyplot/fit_circle.py, Python, 33 lines
examples/pyplot/fit_curve1.py, Python, 32 lines
examples/pyplot/fit_curve2.py, Python, 45 lines
examples/pyplot/fit_erf.py, Python, 51 lines
examples/pyplot/fit_polynomial1.py, Python, 44 lines
examples/pyplot/fit_polynomial2.py, Python, 69 lines
examples/pyplot/fonts3d.py, Python, 114 lines
examples/pyplot/goniometer.py, Python, 29 lines
examples/pyplot/graph_lineage.py, Python, 38 lines
examples/pyplot/graph_network.py, Python, 51 lines
examples/pyplot/histo_1d_a.py, Python, 8 lines
examples/pyplot/histo_1d_b.py, Python, 38 lines
examples/pyplot/histo_1d_c.py, Python, 33 lines
examples/pyplot/histo_1d_d.py, Python, 44 lines
examples/pyplot/histo_1d_e.py, Python, 25 lines
examples/pyplot/histo_2d_a.py, Python, 27 lines
examples/pyplot/histo_2d_b.py, Python, 29 lines
examples/pyplot/histo_3d.py, Python, 19 lines
examples/pyplot/histo_gauss.py, Python, 32 lines
examples/pyplot/histo_hexagonal.py, Python, 21 lines
examples/pyplot/histo_manual.py, Python, 88 lines
examples/pyplot/histo_pca.py, Python, 37 lines
examples/pyplot/histo_polar.py, Python, 46 lines
examples/pyplot/histo_spheric.py, Python, 15 lines
examples/pyplot/histo_violin.py, Python, 27 lines
examples/pyplot/intersect2d.py, Python, 34 lines
examples/pyplot/isolines.py, Python, 35 lines
examples/pyplot/latex.py, Python, 16 lines
examples/pyplot/lines_intersect.py, Python, 27 lines
examples/pyplot/markers.py, Python, 32 lines
examples/pyplot/markpoint.py, Python, 17 lines
examples/pyplot/np_matrix.py, Python, 24 lines
examples/pyplot/pie_chart.py, Python, 17 lines
examples/pyplot/plot_bars.py, Python, 33 lines
examples/pyplot/plot_density2d.py, Python, 26 lines
examples/pyplot/plot_density3d.py, Python, 23 lines
examples/pyplot/plot_density4d.py, Python, 42 lines
examples/pyplot/plot_empty.py, Python, 37 lines
examples/pyplot/plot_errband.py, Python, 47 lines
examples/pyplot/plot_errbars.py, Python, 46 lines
examples/pyplot/plot_extra_yaxis.py, Python, 63 lines
examples/pyplot/plot_fxy0.py, Python, 61 lines
examples/pyplot/plot_fxy1.py, Python, 58 lines
examples/pyplot/plot_fxy2.py, Python, 115 lines
examples/pyplot/plot_hexcells.py, Python, 25 lines
examples/pyplot/plot_multi.py, Python, 24 lines
examples/pyplot/plot_pip.py, Python, 45 lines
examples/pyplot/plot_polar.py, Python, 38 lines
examples/pyplot/plot_spheric.py, Python, 27 lines
examples/pyplot/plot_stream.py, Python, 36 lines
examples/pyplot/quiver.py, Python, 13 lines
examples/pyplot/run_all.sh, Shell, 14 lines
examples/pyplot/scatter1.py, Python, 21 lines
examples/pyplot/scatter2.py, Python, 57 lines
examples/pyplot/scatter3.py, Python, 52 lines
examples/pyplot/scatter_large.py, Python, 22 lines
examples/pyplot/triangulate2d.py, Python, 14 lines
examples/pyplot/whiskers.py, Python, 46 lines
examples/run_all.sh, Shell, 19 lines
examples/volumetric/app_isobrowser.py, Python, 11 lines
examples/volumetric/app_raycaster.py, Python, 13 lines
examples/volumetric/colorize_volume.py, Python, 34 lines
examples/volumetric/delaunay3d.py, Python, 28 lines
examples/volumetric/densifycloud.py, Python, 19 lines
examples/volumetric/earth_model.py, Python, 84 lines
examples/volumetric/erode_dilate.py, Python, 13 lines
examples/volumetric/euclidian_dist.py, Python, 12 lines
examples/volumetric/interpolate_volume.py, Python, 36 lines
examples/volumetric/isosurfaces1.py, Python, 19 lines
examples/volumetric/isosurfaces2.py, Python, 53 lines
examples/volumetric/legosurface.py, Python, 16 lines
examples/volumetric/lowpassfilter.py, Python, 18 lines
examples/volumetric/mesh2volume.py, Python, 17 lines
examples/volumetric/multiscalars.py, Python, 39 lines
examples/volumetric/numpy2volume0.py, Python, 13 lines
examples/volumetric/numpy2volume1.py, Python, 21 lines
examples/volumetric/numpy2volume2.py, Python, 16 lines
examples/volumetric/numpy_imread.py, Python, 25 lines
examples/volumetric/off_furniture.py, Python, 245 lines
examples/volumetric/office.py, Python, 23 lines
examples/volumetric/probe_line1.py, Python, 22 lines
examples/volumetric/probe_line2.py, Python, 34 lines
examples/volumetric/probe_points.py, Python, 19 lines
examples/volumetric/read_volume1.py, Python, 28 lines
examples/volumetric/read_volume2.py, Python, 32 lines
examples/volumetric/read_volume3.py, Python, 14 lines
examples/volumetric/rectl_grid1.py, Python, 21 lines
examples/volumetric/run_all.sh, Shell, 14 lines
examples/volumetric/slab_vol.py, Python, 22 lines
examples/volumetric/slice_mesh.py, Python, 18 lines
examples/volumetric/slice_plane1.py, Python, 33 lines
examples/volumetric/slice_plane2.py, Python, 24 lines
examples/volumetric/slice_plane3.py, Python, 39 lines
examples/volumetric/slicer1.py, Python, 22 lines
examples/volumetric/slicer2.py, Python, 45 lines
examples/volumetric/slicer_set_volume.py, Python, 65 lines
examples/volumetric/streamlines1.py, Python, 37 lines
examples/volumetric/streamlines2.py, Python, 15 lines
examples/volumetric/streamlines2_cylinderwidget.py, Python, 37 lines
examples/volumetric/streamlines2_linewidget.py, Python, 39 lines
examples/volumetric/streamlines2_pointwidget.py, Python, 35 lines
examples/volumetric/streamlines2_spherewidget.py, Python, 35 lines
examples/volumetric/streamlines3.py, Python, 23 lines
examples/volumetric/streamlines4.py, Python, 25 lines
examples/volumetric/struc_grid1.py, Python, 27 lines
examples/volumetric/tensors.py, Python, 31 lines
examples/volumetric/tet_astyle.py, Python, 18 lines
examples/volumetric/tet_build.py, Python, 44 lines
examples/volumetric/tet_cut1.py, Python, 19 lines
examples/volumetric/tet_cut2.py, Python, 38 lines
examples/volumetric/tet_explode.py, Python, 40 lines
examples/volumetric/tet_isos_slice.py, Python, 25 lines
examples/volumetric/tet_threshold.py, Python, 17 lines
examples/volumetric/tetralize_surface.py, Python, 31 lines
examples/volumetric/ugrid1.py, Python, 13 lines
examples/volumetric/ugrid2.py, Python, 12 lines
examples/volumetric/vol2points.py, Python, 13 lines
examples/volumetric/volume_from_mesh.py, Python, 16 lines
examples/volumetric/volume_operations.py, Python, 59 lines
examples/volumetric/volume_sharemap.py, Python, 27 lines
examples/volumetric/warp_scalars.py, Python, 15 lines
tests/common/run_all.sh, Shell, 14 lines
tests/common/test_0_imports.py, Python, 10 lines
tests/common/test_bins.py, Python, 28 lines
tests/common/test_chemistry_protein.py, Python, 59 lines
tests/common/test_cli.py, Python, 156 lines
tests/common/test_cylinder_widget.py, Python, 78 lines
tests/common/test_data_format.py, Python, 66 lines
tests/common/test_dataset_algorithms.py, Python, 182 lines
tests/common/test_explicit_grid.py, Python, 71 lines
tests/common/test_file_io_cube.py, Python, 28 lines
tests/common/test_file_io_opendx.py, Python, 73 lines
tests/common/test_file_io_writers.py, Python, 30 lines
tests/common/test_hist.py, Python, 26 lines
tests/common/test_keymap.py, Python, 100 lines
tests/common/test_lazy_exports.py, Python, 63 lines
tests/common/test_line_widget.py, Python, 90 lines
tests/common/test_logging_rich.py, Python, 50 lines
tests/common/test_mesh_core.py, Python, 177 lines
tests/common/test_numpy_hist.py, Python, 18 lines
tests/common/test_numpy_hist2.py, Python, 18 lines
tests/common/test_numpy_hist3.py, Python, 15 lines
tests/common/test_plotter_dry_run.py, Python, 59 lines
tests/common/test_plotter_modes.py, Python, 126 lines
tests/common/test_plotter_scene.py, Python, 31 lines
tests/common/test_point_widget.py, Python, 44 lines
tests/common/test_pointcloud_analyze.py, Python, 152 lines
tests/common/test_pointcloud_fits.py, Python, 139 lines
tests/common/test_pyplot.py, Python, 122 lines
tests/common/test_pyplot_functions.py, Python, 390 lines
tests/common/test_pyplot_stats.py, Python, 26 lines
tests/common/test_scalarbar.py, Python, 27 lines
tests/common/test_shapes.py, Python, 233 lines
tests/common/test_slicer3d_plotter.py, Python, 80 lines
tests/common/test_smoke_imports.py, Python, 20 lines
tests/common/test_sphere_widget.py, Python, 66 lines
tests/common/test_summary_rich.py, Python, 30 lines
tests/common/test_summary_rich_rollout.py, Python, 77 lines
tests/common/test_text2d_and_colors.py, Python, 140 lines
tests/common/test_trame_backend.py, Python, 40 lines
tests/common/test_transformations.py, Python, 113 lines
tests/common/test_transformations_quaternion.py, Python, 60 lines
tests/common/test_utils.py, Python, 34 lines
tests/common/test_volume_api.py, Python, 95 lines
tests/common/test_vtkclasses_slow_notice.py, Python, 72 lines
tests/conftest.py, Python, 31 lines
tests/issues/discussion_1190.py, Python, 37 lines
tests/issues/discussion_527.py, Python, 37 lines
tests/issues/discussion_716.py, Python, 16 lines
tests/issues/discussion_751.py, Python, 26 lines
tests/issues/discussion_784.py, Python, 13 lines
tests/issues/discussion_800.py, Python, 36 lines
tests/issues/discussion_942.py, Python, 68 lines
tests/issues/discussion_978.py, Python, 56 lines
tests/issues/issue_1025.py, Python, 43 lines
tests/issues/issue_1077.py, Python, 16 lines
tests/issues/issue_1109.py, Python, 46 lines
tests/issues/issue_1118.py, Python, 127 lines
tests/issues/issue_1146.py, Python, 59 lines
tests/issues/issue_1218.py, Python, 94 lines
tests/issues/issue_1221.py, Python, 169 lines
tests/issues/issue_1230.py, Python, 36 lines
tests/issues/issue_656.py, Python, 9 lines
tests/issues/issue_805.py, Python, 15 lines
tests/issues/issue_851.py, Python, 13 lines
tests/issues/issue_854.py, Python, 18 lines
tests/issues/issue_856.py, Python, 18 lines
tests/issues/issue_871a.py, Python, 123 lines
tests/issues/issue_871b.py, Python, 31 lines
tests/issues/issue_893.py, Python, 33 lines
tests/issues/issue_905.py, Python, 59 lines
tests/issues/issue_908.py, Python, 35 lines
tests/issues/issue_939.py, Python, 5 lines
tests/issues/issue_946.py, Python, 58 lines
tests/issues/issue_948.py, Python, 35 lines
tests/issues/issue_950.py, Python, 17 lines
tests/issues/issue_953.py, Python, 58 lines
tests/issues/issue_968.py, Python, 10 lines
tests/issues/run_all.sh, Shell, 11 lines
tests/issues/test_force_anim.py, Python, 76 lines
tests/issues/test_fxy_bessel1.py, Python, 41 lines
tests/issues/test_remove_objects.py, Python, 37 lines
tests/issues/test_sph_harm2.py, Python, 97 lines
tests/snippets/run_all.sh, Shell, 11 lines
tests/snippets/test_cell_types.py, Python, 177 lines
tests/snippets/test_closewindow.py, Python, 50 lines
tests/snippets/test_compare_fit1.py, Python, 81 lines
tests/snippets/test_discourse_1956.py, Python, 9 lines
tests/snippets/test_docs_sniplets.py, Python, 350 lines
tests/snippets/test_elastic_pendulum.py, Python, 70 lines
tests/snippets/test_ellipsoid_main_axes.py, Python, 41 lines
tests/snippets/test_fit_shapiro.py, Python, 67 lines
tests/snippets/test_interactive_plotxy1.py, Python, 26 lines
tests/snippets/test_interactive_plotxy2.py, Python, 66 lines
tests/snippets/test_project_point_on_variety.py, Python, 18 lines
vedo/__init__.py, Python, 525 lines
vedo/addons/__init__.py, Python, 47 lines
vedo/addons/axes.py, Python, 1,833 lines
vedo/addons/core.py, Python, 690 lines
vedo/addons/cutters.py, Python, 730 lines
vedo/addons/interaction.py, Python, 1,129 lines
vedo/addons/measure.py, Python, 581 lines
vedo/addons/sliders.py, Python, 463 lines
vedo/addons/ui.py, Python, 423 lines
vedo/addons/widgets.py, Python, 318 lines
vedo/applications/__init__.py, Python, 40 lines
vedo/applications/animation.py, Python, 757 lines
vedo/applications/browsers.py, Python, 345 lines
vedo/applications/chemistry.py, Python, 752 lines
vedo/applications/editing.py, Python, 607 lines
vedo/applications/morphing.py, Python, 827 lines
vedo/applications/slicing.py, Python, 1,089 lines
vedo/assembly.py, Python, 707 lines
vedo/backends.py, Python, 595 lines
vedo/cli.py, Python, 1,635 lines
vedo/colors.py, Python, 1,378 lines
vedo/core/__init__.py, Python, 22 lines
vedo/core/common.py, Python, 1,415 lines
vedo/core/data.py, Python, 377 lines
vedo/core/input.py, Python, 70 lines
vedo/core/points.py, Python, 295 lines
vedo/core/summary.py, Python, 122 lines
vedo/core/transformations.py, Python, 1,556 lines
vedo/core/volume.py, Python, 289 lines
vedo/external/__init__.py, Python, 29 lines
vedo/external/conversions.py, Python, 357 lines
vedo/file_io/__init__.py, Python, 50 lines
vedo/file_io/export_npy.py, Python, 26 lines
vedo/file_io/export_threejs.py, Python, 1,005 lines
vedo/file_io/export_x3d.py, Python, 155 lines
vedo/file_io/loaders.py, Python, 1,090 lines
vedo/file_io/network.py, Python, 107 lines
vedo/file_io/scene.py, Python, 588 lines
vedo/file_io/terminal.py, Python, 58 lines
vedo/file_io/video.py, Python, 245 lines
vedo/file_io/writers.py, Python, 200 lines
vedo/grids/__init__.py, Python, 23 lines
vedo/grids/explicit.py, Python, 486 lines
vedo/grids/image.py, Python, 1,584 lines
vedo/grids/rectilinear.py, Python, 526 lines
vedo/grids/structured.py, Python, 487 lines
vedo/grids/tetmesh.py, Python, 527 lines
vedo/grids/unstructured.py, Python, 969 lines
vedo/lazy_imports.py, Python, 47 lines
vedo/mesh/__init__.py, Python, 17 lines
vedo/mesh/core.py, Python, 2,884 lines
vedo/mesh/metrics.py, Python, 399 lines
vedo/plotter/__init__.py, Python, 16 lines
vedo/plotter/camera.py, Python, 270 lines
vedo/plotter/events.py, Python, 91 lines
vedo/plotter/interaction.py, Python, 563 lines
vedo/plotter/io.py, Python, 78 lines
vedo/plotter/keymap.py, Python, 1,092 lines
vedo/plotter/lifecycle.py, Python, 271 lines
vedo/plotter/modes.py, Python, 1,571 lines
vedo/plotter/runtime.py, Python, 2,583 lines
vedo/plotter/scene.py, Python, 289 lines
vedo/plotter/session.py, Python, 56 lines
vedo/pointcloud/__init__.py, Python, 33 lines
vedo/pointcloud/analyze.py, Python, 1,257 lines
vedo/pointcloud/core.py, Python, 500 lines
vedo/pointcloud/cut.py, Python, 658 lines
vedo/pointcloud/fits.py, Python, 589 lines
vedo/pointcloud/reconstruct.py, Python, 675 lines
vedo/pointcloud/transform.py, Python, 487 lines
vedo/pyplot/__init__.py, Python, 24 lines
vedo/pyplot/charts.py, Python, 1,263 lines
vedo/pyplot/figure.py, Python, 655 lines
vedo/pyplot/functions.py, Python, 1,695 lines
vedo/pyplot/graph.py, Python, 375 lines
vedo/pyplot/stats.py, Python, 585 lines
vedo/settings.py, Python, 816 lines
vedo/shapes/__init__.py, Python, 88 lines
vedo/shapes/analysis.py, Python, 95 lines
vedo/shapes/branding.py, Python, 47 lines
vedo/shapes/curves.py, Python, 54 lines
vedo/shapes/curves_core.py, Python, 1,538 lines
vedo/shapes/curves_extras.py, Python, 602 lines
vedo/shapes/glyphs.py, Python, 253 lines
vedo/shapes/latex.py, Python, 103 lines
vedo/shapes/markers.py, Python, 425 lines
vedo/shapes/primitives.py, Python, 64 lines
vedo/shapes/primitives_planar.py, Python, 693 lines
vedo/shapes/primitives_solids.py, Python, 756 lines
vedo/shapes/text.py, Python, 15 lines
vedo/shapes/text2d.py, Python, 345 lines
vedo/shapes/text3d.py, Python, 370 lines
vedo/shapes/text_utils.py, Python, 164 lines
vedo/transformations.py, Python, 7 lines
vedo/utils.py, Python, 2,883 lines
vedo/visual/__init__.py, Python, 29 lines
vedo/visual/mixins.py, Python, 1,061 lines
vedo/visual/runtime.py, Python, 2,099 lines
vedo/volume/__init__.py, Python, 17 lines
vedo/volume/core.py, Python, 1,513 lines
vedo/volume/slicing.py, Python, 273 lines
vedo/vtkclasses.py, Python, 1,065 lines
LICENSE, License, 21 lines
README.md, Text, 306 lines

At the source: github.com/marcomusy/vedo

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;
636 scripts, each with its path and the digest of its content;
7 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

Datasets cited

doi:10.48324/dandi.000774/0.260520.1753, at DANDI; found in the resources table

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 2, 28 September 2026

Publisher: n/a → Cell Press

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 9 keywords, 10 MeSH terms, 5 funders, 120 references, 1 RRID.

Cite

This paper

Hickman, J. L., Hughes, G., Sahai, E., Miles, M., & Denman, D. J. (2026). Neural population dynamics of direct electrical stimulation of neocortex. Cell reports, 45(6), 117420. https://doi.org/10.1016/j.celrep.2026.117420

BibTeX

@article{hickman2026neural, author = {Hickman, Jordan L and Hughes, Grant and Sahai, Eashan and Miles, Moriah and Denman, Daniel J}, title = {{Neural population dynamics of direct electrical stimulation of neocortex}}, journal = {Cell reports}, year = {2026}, month = may, volume = {45}, number = {6}, pages = {117420}, publisher = {Cell Press}, issn = {2211-1247}, doi = {10.1016/j.celrep.2026.117420}, url = {https://doi.org/10.1016/j.celrep.2026.117420}, pmid = {42207642}, pmcid = {PMC13404931} }

RIS

TY - JOUR AU - Hickman, Jordan L AU - Hughes, Grant AU - Sahai, Eashan AU - Miles, Moriah AU - Denman, Daniel J TI - Neural population dynamics of direct electrical stimulation of neocortex T2 - Cell reports J2 - Cell Rep PY - 2026 DA - 2026/05/28 VL - 45 IS - 6 SP - 117420 SN - 2211-1247 PB - Cell Press DO - 10.1016/j.celrep.2026.117420 UR - https://doi.org/10.1016/j.celrep.2026.117420 LA - en ER -

CSL-JSON

{ "id": "10.1016/j.celrep.2026.117420", "type": "article-journal", "title": "Neural population dynamics of direct electrical stimulation of neocortex", "container-title": "Cell reports", "author": [ { "family": "Hickman", "given": "Jordan L" }, { "family": "Hughes", "given": "Grant" }, { "family": "Sahai", "given": "Eashan" }, { "family": "Miles", "given": "Moriah" }, { "family": "Denman", "given": "Daniel J" } ], "container-title-short": "Cell Rep", "volume": "45", "issue": "6", "page": "117420", "DOI": "10.1016/j.celrep.2026.117420", "PMID": "42207642", "PMCID": "PMC13404931", "ISSN": "2211-1247", "publisher": "Cell Press", "URL": "https://doi.org/10.1016/j.celrep.2026.117420", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 28 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

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.

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.

Neural population dynamics of direct electrical stimulation of neocortex.

Code ↔ Paper

Paper

The authors' code

Overview

Abstract

Repositories

denmanlab/am4100_code

denmanlab/mouse_behavior

denmanlab/estim_populations

marcomusy/vedo

Tracing map

Data

Datasets cited

Versions

Version 2, 28 September 2026

Version 1, 28 September 2026: the first record

Cite

This paper

Similar papers

Contribute

Claim this paper

Correct its record

Validate its tracing map

Add the badge to its README

Request its removal

Discussion, reproductions, activity