# install cell skipped during CI (deps preinstalled into system Python)

import h5py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import remfile
import seaborn as sns
from dandi.dandiapi import DandiAPIClient
from pynwb import NWBHDF5IO
from scipy import stats
from tqdm import tqdm

# Set plotting style to match paper
plt.style.use('default')
sns.set_palette("Set2")

def setup_figure_style():
    """Setup matplotlib parameters to match paper style"""
    plt.rcParams.update({
        'font.size': 8,
        'axes.titlesize': 10,
        'axes.labelsize': 9,
        'xtick.labelsize': 8,
        'ytick.labelsize': 8,
        'legend.fontsize': 8,
        'figure.titlesize': 12,
        'axes.linewidth': 0.8,
        'axes.spines.top': False,
        'axes.spines.right': False,
        'xtick.major.width': 0.8,
        'ytick.major.width': 0.8,
        'xtick.minor.width': 0.6,
        'ytick.minor.width': 0.6,
    })

setup_figure_style()
print("Libraries imported and plotting style configured")

Libraries imported and plotting style configured

# Session ID helpers for the per-mouse-day dandiset format.
# DANDI asset paths look like:
#   sub-Subject<date><cat><animal>/sub-..._ses-<cellType>++<state>++<pharm>++<geno>++<date>_<modality>.nwb
# The 5-token session_id encodes the experimental cohort (no figure number).

def get_session_id(asset_path):
    fname = asset_path.split('/')[-1]
    if '_ses-' not in fname:
        return ''
    return fname.split('_ses-', 1)[1].rsplit('_', 1)[0]

def parse_session_id(ses):
    t = ses.split('++')
    if len(t) < 5:
        return {}
    return {'cell_type': t[0], 'state': t[1], 'pharm': t[2], 'genotype': t[3], 'date': t[4]}

def get_modality(asset_path):
    fname = asset_path.split('/')[-1]
    if not fname.endswith('.nwb'):
        return ''
    return fname[:-len('.nwb')].rsplit('_', 1)[-1]

def is_f1_dspn_somatic(asset_path):
    """Figure 1 dSPN somatic excitability: per-mouse-day icephys files for dSPN cell type."""
    fields = parse_session_id(get_session_id(asset_path))
    if not fields: return False
    return get_modality(asset_path) in ('icephys', 'icephys+ophys') and fields['cell_type'] == 'dSPN'

def get_condition_label(condition_subfolder):
    """Pretty condition label from the IRT condition_subfolder column."""
    return condition_subfolder  # The condition column is already human-readable

print('Utility functions defined for new dandiset format.')

Utility functions defined for new dandiset format.

def count_action_potentials(
    voltage_trace_mV: np.ndarray, timestamps_s: np.ndarray, threshold_mV: float = 0.0
) -> int:
    """Threshold-crossing spike count within a 500 ms stimulus window.

    Matches approach in the original analysis; uses 200-700 ms from sweep start.
    """
    if voltage_trace_mV.size == 0:
        return 0
    if timestamps_s[-1] - timestamps_s[0] < 0.8:
        return 0
    start_t = timestamps_s[0] + 0.2
    end_t = timestamps_s[0] + 0.7
    i0 = np.searchsorted(timestamps_s, start_t)
    i1 = np.searchsorted(timestamps_s, end_t)
    if i0 >= i1 or i1 > voltage_trace_mV.size:
        i0 = voltage_trace_mV.size // 4
        i1 = 3 * voltage_trace_mV.size // 4
    x = voltage_trace_mV[i0:i1]
    if x.size == 0:
        return 0
    spikes = 0
    i = 0
    while i < x.size:
        if x[i] > threshold_mV:
            spikes += 1
            while i < x.size and x[i] > threshold_mV:
                i += 1
        else:
            i += 1
    return spikes

def calculate_rheobase(current_steps: list[float], spike_counts: list[int]) -> float:
    """Calculate rheobase as minimum current that elicits at least one spike."""
    for current, spikes in zip(current_steps, spike_counts):
        if spikes >= 1:
            return current
    return np.nan

def safe_mean_sem(values: np.ndarray) -> tuple[float, float]:
    """Calculate mean and standard error, handling edge cases."""
    if len(values) == 0:
        return np.nan, 0.0
    if len(values) == 1:
        return float(values[0]), 0.0
    return float(np.mean(values)), float(np.std(values, ddof=1) / np.sqrt(len(values)))

print("Analysis functions defined")

Analysis functions defined

from collections import Counter

dandiset_id = "001538"
client = DandiAPIClient()
dandiset = client.get_dandiset(dandiset_id, "0.260527.1302")
assets_list = list(dandiset.get_assets())

# Filter for F1 dSPN somatic excitability: per-mouse-day icephys files with dSPN cell type
f1_somexc_assets = [a for a in assets_list if is_f1_dspn_somatic(a.path)]
print(f'Found {len(f1_somexc_assets)} F1 dSPN per-mouse-day files')

# Breakdown by condition (state + pharm)
cohort = Counter()
for a in f1_somexc_assets:
    f = parse_session_id(get_session_id(a.path))
    cohort[(f['state'], f['pharm'])] += 1
for k, n in sorted(cohort.items()):
    state, pharm = k
    pretty = {'OffState': 'LID off-state', 'OnState': 'LID on-state'}.get(state, state)
    if pharm != 'none':
        pretty += f' with {pharm}'
    print(f'  {pretty}: {n} mouse-days')

Found 18 F1 dSPN per-mouse-day files
  LID off-state: 8 mouse-days
  LID on-state with D1RaSch: 5 mouse-days
  LID on-state: 5 mouse-days

# Per-cell spike counting using Pattern B slicing.
# Each per-mouse-day file has MULTIPLE cells; we group by (subject, cell_id) and
# compute AP counts per current step using index_start + count from the IRT.

all_recording_data = []   # one row per (cell, current_step) sweep observation
cell_statistics = []      # per-cell aggregates (for rheobase)

print(f'Processing {len(f1_somexc_assets)} per-mouse-day files...')

for i, asset in enumerate(tqdm(f1_somexc_assets, desc='Streaming files')):
    s3_url = asset.get_content_url(follow_redirects=1, strip_query=False)
    file_system = remfile.File(s3_url, disk_cache=remfile.DiskCache("nwb-cache"))
    file = h5py.File(file_system, mode='r')
    io = NWBHDF5IO(file=file, load_namespaces=True)
    nwbfile = io.read()

    subject_id = nwbfile.subject.subject_id
    rec_df = nwbfile.intracellular_recordings.to_dataframe()

    # Filter to somatic sweeps only (each file has both somatic + dendritic)
    soma_mask = rec_df[('intracellular_recordings', 'dendrite_type')] == 'Soma'
    soma_rows = rec_df[soma_mask]

    # Group by cell_id within this mouse-day
    per_cell_data = {}  # (cell_id, condition) -> {step_pA: [ap_counts]}
    for _, row in soma_rows.iterrows():
        cell_id = int(row[('intracellular_recordings', 'cell_id')])
        condition = row[('intracellular_recordings', 'condition_subfolder')]
        current_pA = float(row[('intracellular_recordings', 'stimulus_current_pA')])

        # Pattern B slice into the consolidated CurrentClampSeries
        resp_ref = row[('responses', 'response')]
        ts = resp_ref.timeseries
        index_start = int(resp_ref.idx_start)
        count = int(resp_ref.count)
        voltage_v = ts.data[index_start: index_start + count]
        voltage_mV = voltage_v * 1000.0

        rate = float(ts.rate)
        timestamps_s = np.arange(count) / rate

        # Count APs in the 200-700 ms step window
        spike_count = count_action_potentials(voltage_mV, timestamps_s, threshold_mV=0.0)

        key = (cell_id, condition)
        if key not in per_cell_data:
            per_cell_data[key] = {}
        per_cell_data[key].setdefault(current_pA, []).append(spike_count)

        all_recording_data.append({
            'subject_id': subject_id,
            'cell_id': cell_id,
            'condition': condition,
            'current_pA': current_pA,
            'spike_count': spike_count,
        })

    # Per-cell aggregates: max APs across steps, rheobase (min step with >=1 AP)
    for (cell_id, condition), by_step in per_cell_data.items():
        steps = sorted(by_step.keys())
        means = {s: np.mean(by_step[s]) for s in steps}
        rheo = next((s for s in steps if means[s] >= 1), None)
        cell_statistics.append({
            'subject_id': subject_id,
            'cell_id': cell_id,
            'condition': condition,
            'max_aps': max(means.values()) if means else 0,
            'rheobase_pA': rheo,
        })

df_recordings = pd.DataFrame(all_recording_data)
df_cells = pd.DataFrame(cell_statistics)
print(f'\nProcessed {len(df_cells)} unique (subject, cell, condition) combinations')
print(f'Total sweep observations: {len(df_recordings)}')
print(df_cells.head())

Processing 18 per-mouse-day files...

Streaming files:   0%|          | 0/18 [00:00<?, ?it/s]

Streaming files:   6%|▌         | 1/18 [00:01<00:25,  1.49s/it]

Streaming files:  11%|█         | 2/18 [00:04<00:33,  2.11s/it]

Streaming files:  17%|█▋        | 3/18 [00:06<00:34,  2.29s/it]

Streaming files:  22%|██▏       | 4/18 [00:08<00:29,  2.13s/it]

Streaming files:  28%|██▊       | 5/18 [00:10<00:28,  2.18s/it]

Streaming files:  33%|███▎      | 6/18 [00:13<00:27,  2.33s/it]

Streaming files:  39%|███▉      | 7/18 [00:16<00:28,  2.56s/it]

Streaming files:  44%|████▍     | 8/18 [00:17<00:20,  2.07s/it]

# Create F-I curves plot
fig, ax = plt.subplots(1, 1, figsize=(4.5, 3.5))  # Adjusted aspect ratio to match reference

# Define condition plotting styles to match reference - use only circles with larger size
condition_styles = {
    "LID off-state": {
        "color": "black", 
        "marker": "o", 
        "linestyle": "-", 
        "label": "off-state",
        "markerfacecolor": "white",
        "markeredgecolor": "black"
    },
    "LID on-state": {
        "color": "black", 
        "marker": "o",  # Changed from "s" to "o"
        "linestyle": "-", 
        "label": "on-state",
        "markerfacecolor": "black",
        "markeredgecolor": "black"
    },
    "LID on-state with SCH": {
        "color": "gray",  # Changed back to gray for antagonist
        "marker": "o",  # Changed from "^" to "o" 
        "linestyle": "-",  # Changed from "--" to "-"
        "label": "on-state+D1R\nantagonist",
        "markerfacecolor": "gray",
        "markeredgecolor": "gray"
    }
}

# Process each condition
for condition in ["LID off-state", "LID on-state", "LID on-state with SCH"]:
    if condition not in df_recordings["condition"].unique():
        print(f"Warning: {condition} not found in data")
        continue
    
    condition_data = df_recordings[df_recordings["condition"] == condition]
    
    # Calculate mean and SEM for each current step
    summary_data = []
    for current, group in condition_data.groupby("current_pA"):
        mean_spikes, sem_spikes = safe_mean_sem(group["spike_count"].values)
        summary_data.append({
            "current_pA": current, 
            "mean_spikes": mean_spikes, 
            "sem_spikes": sem_spikes
        })
    
    summary_df = pd.DataFrame(summary_data).sort_values("current_pA")
    
    # Filter to current range used in paper (0-300 pA)
    summary_df = summary_df[
        (summary_df["current_pA"] >= 0) & (summary_df["current_pA"] <= 300)
    ]
    
    if len(summary_df) == 0:
        print(f"Warning: No data in 0-300pA range for {condition}")
        continue
    
    # Plot with error bars - larger markers
    style = condition_styles[condition]
    ax.errorbar(
        summary_df["current_pA"], 
        summary_df["mean_spikes"], 
        yerr=summary_df["sem_spikes"],
        marker=style["marker"], 
        color=style["color"], 
        linestyle=style["linestyle"],
        linewidth=2, 
        markersize=8,  # Increased from 4 to 8
        capsize=3, 
        capthick=1, 
        label=style["label"],
        markerfacecolor=style["markerfacecolor"],
        markeredgecolor=style["markeredgecolor"], 
        markeredgewidth=1.5,  # Increased edge width
    )

# Add dotted horizontal line at y=0 to show baseline
ax.axhline(y=0, color='gray', linestyle=':', alpha=0.7, linewidth=1)

# Formatting to match paper style
ax.set_xlabel("current (pA)", fontsize=12)
ax.set_ylabel("number of APs", fontsize=12)
ax.set_xlim(0, 300)
ax.set_ylim(-1, 18)  # Extended y-limit to include negative values
ax.set_xticks([0, 100, 200, 300])
ax.set_yticks([0, 5, 10, 15])

# Style the axes
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_linewidth(1.5)
ax.spines['bottom'].set_linewidth(1.5)
ax.tick_params(axis='both', which='major', labelsize=10, width=1.5, length=5)

# Add legend in top left corner instead of text labels on plot
ax.legend(loc="upper left", frameon=False, fontsize=9)

plt.tight_layout()
plt.show()

# Print summary statistics
print("=== F-I CURVE ANALYSIS SUMMARY ===\n")
for condition in df_recordings["condition"].unique():
    condition_data = df_recordings[df_recordings["condition"] == condition]
    n_recordings = len(condition_data)
    n_cells = len(df_cells[df_cells["condition"] == condition])
    
    # Calculate statistics at different current levels
    current_levels = [50, 100, 150, 200, 250]
    print(f"{condition}:")
    print(f"  Sample: {n_recordings} recordings from {n_cells} cells")
    
    for current in current_levels:
        current_data = condition_data[
            (condition_data["current_pA"] >= current - 10) & 
            (condition_data["current_pA"] <= current + 10)
        ]
        if len(current_data) > 0:
            mean_spikes = current_data["spike_count"].mean()
            sem_spikes = current_data["spike_count"].std() / np.sqrt(len(current_data))
            print(f"  {current}pA: {mean_spikes:.1f} ± {sem_spikes:.1f} APs (n={len(current_data)})")
    print()

=== F-I CURVE ANALYSIS SUMMARY ===

LID off-state:
  Sample: 358 recordings from 15 cells
  50pA: 0.0 ± 0.0 APs (n=31)
  100pA: 0.0 ± 0.0 APs (n=15)
  150pA: 0.1 ± 0.0 APs (n=29)
  200pA: 2.2 ± 0.8 APs (n=16)
  250pA: 5.5 ± 1.0 APs (n=28)

LID on-state:
  Sample: 309 recordings from 15 cells
  50pA: 0.1 ± 0.1 APs (n=30)
  100pA: 1.1 ± 0.6 APs (n=14)
  150pA: 5.6 ± 0.9 APs (n=31)
  200pA: 10.7 ± 1.1 APs (n=15)
  250pA: 13.6 ± 0.5 APs (n=29)

LID on-state with SCH:
  Sample: 210 recordings from 10 cells
  50pA: 0.0 ± 0.0 APs (n=20)
  100pA: 1.2 ± 0.8 APs (n=10)
  150pA: 6.3 ± 0.9 APs (n=20)
  200pA: 11.4 ± 0.8 APs (n=10)
  250pA: 14.7 ± 0.5 APs (n=20)

# Create rheobase comparison plot
fig, ax = plt.subplots(1, 1, figsize=(3.5, 4))  # Adjusted to match reference aspect ratio

# Filter out cells with invalid rheobase values
valid_cells = df_cells.dropna(subset=["rheobase_pA"])

# Prepare data for box plot
conditions_order = ["LID off-state", "LID on-state", "LID on-state with SCH"]
condition_labels = ["off", "on", "on + SCH"]  # Shorter labels like reference

# Get data for each condition
plot_data = []
actual_labels = []

for condition, label in zip(conditions_order, condition_labels):
    condition_data = valid_cells[valid_cells["condition"] == condition]["rheobase_pA"]
    if len(condition_data) > 0:
        plot_data.append(condition_data.values)
        actual_labels.append(label)
    else:
        print(f"Warning: No rheobase data for {condition}")

# Create box plot with styling to match reference
bp = ax.boxplot(
    plot_data, 
    labels=actual_labels, 
    patch_artist=True,
    boxprops=dict(facecolor="white", color="black", linewidth=1.5),
    whiskerprops=dict(color="black", linewidth=1.5), 
    capprops=dict(color="black", linewidth=1.5),
    medianprops=dict(color="black", linewidth=3),  # Thicker median line
    flierprops=dict(marker="o", markerfacecolor="gray", markersize=3, 
                  markeredgecolor="black", alpha=0.7),
    widths=0.6  # Slightly narrower boxes like reference
)

# Add individual data points with jitter - all gray
for i, data in enumerate(plot_data):
    x_vals = np.random.normal(i + 1, 0.05, len(data))
    ax.scatter(x_vals, data, color='gray', alpha=0.8, s=20, zorder=3, 
               edgecolors='black', linewidths=0.5)  # Larger markers (s=20)

# Add significance bracket and annotation like in reference
# Add bracket between first and second conditions
y_max = max([max(data) for data in plot_data]) + 20
bracket_height = y_max + 10
ax.plot([1, 1], [y_max, bracket_height], 'k-', linewidth=1)
ax.plot([2, 2], [y_max, bracket_height], 'k-', linewidth=1)
ax.plot([1, 2], [bracket_height, bracket_height], 'k-', linewidth=1)
ax.text(1.5, bracket_height + 10, '****', ha='center', va='bottom', fontsize=12, fontweight='bold')

# Formatting to match reference
ax.set_ylabel("rheobase (pA)", fontsize=12)
ax.set_ylim(0, 350)  # Extended to accommodate significance bracket
ax.set_yticks([0, 100, 200, 300])

# Style the axes to match reference
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_linewidth(1.5)
ax.spines['bottom'].set_linewidth(1.5)
ax.tick_params(axis='both', which='major', labelsize=11, width=1.5, length=5)
ax.tick_params(axis='x', which='major', length=0)  # Remove x-axis tick marks

plt.tight_layout()
plt.show()

# Statistical analysis
print("=== RHEOBASE STATISTICAL ANALYSIS ===\n")

for i, (condition, data) in enumerate(zip(conditions_order, plot_data)):
    if condition in valid_cells["condition"].unique():
        rheobase_data = valid_cells[valid_cells["condition"] == condition]["rheobase_pA"]
        n_cells = len(rheobase_data)
        mean_rheo = rheobase_data.mean()
        sem_rheo = rheobase_data.std() / np.sqrt(n_cells)
        median_rheo = rheobase_data.median()
        q25 = rheobase_data.quantile(0.25)
        q75 = rheobase_data.quantile(0.75)
        
        print(f"{condition} (n={n_cells}):")
        print(f"  Mean: {mean_rheo:.1f} ± {sem_rheo:.1f} pA")
        print(f"  Median: {median_rheo:.1f} pA")
        print(f"  IQR: {q25:.1f} - {q75:.1f} pA")
        print(f"  Range: {rheobase_data.min():.1f} - {rheobase_data.max():.1f} pA\n")

# Statistical comparisons
print("Statistical Comparisons:")

# Off-state vs On-state
if "LID off-state" in valid_cells["condition"].unique() and "LID on-state" in valid_cells["condition"].unique():
    off_data = valid_cells[valid_cells["condition"] == "LID off-state"]["rheobase_pA"]
    on_data = valid_cells[valid_cells["condition"] == "LID on-state"]["rheobase_pA"]
    
    # Mann-Whitney U test
    u_stat, u_p = stats.mannwhitneyu(off_data, on_data, alternative='two-sided')
    print(f"\nOff-state vs On-state:")
    print(f"  Mann-Whitney U: {u_stat:.2f}, p = {u_p:.4f}")
    print(f"  Significantly different: {'Yes' if u_p < 0.05 else 'No'}")
    
    # Effect size
    mean_diff = on_data.mean() - off_data.mean()
    print(f"  Mean difference: {mean_diff:.1f} pA")

/tmp/ipykernel_2286/3448983558.py:24: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11.
  bp = ax.boxplot(

=== RHEOBASE STATISTICAL ANALYSIS ===

LID off-state (n=15):
  Mean: 233.3 ± 15.6 pA
  Median: 220.0 pA
  IQR: 180.0 - 270.0 pA
  Range: 160.0 - 340.0 pA

LID on-state (n=15):
  Mean: 125.3 ± 8.4 pA
  Median: 140.0 pA
  IQR: 100.0 - 140.0 pA
  Range: 60.0 - 180.0 pA

LID on-state with SCH (n=10):
  Mean: 126.0 ± 8.5 pA
  Median: 140.0 pA
  IQR: 120.0 - 140.0 pA
  Range: 80.0 - 160.0 pA

Statistical Comparisons:

Off-state vs On-state:
  Mann-Whitney U: 219.50, p = 0.0000
  Significantly different: Yes
  Mean difference: -108.0 pA

Installing requirements¶

Reproducing Figure 1E: dSPN Somatic Excitability Analysis¶

Setup and Data Loading¶

Import Libraries and Configure Plotting Style¶

Session ID Parsing and Filtering Functions¶

Spike Detection and Analysis Functions¶

Action Potential Counting¶

Load DANDI Dataset¶

Data Processing and Analysis¶

Process All NWB Files¶

Figure 1E: Frequency-Intensity (F-I) Curves¶

Action Potential Frequency vs Injected Current¶

Figure 1E: Rheobase Comparison¶

Minimum Current Required for Action Potential Generation¶

Summary¶

Key Findings¶

Methodological Notes¶

Biological Significance¶