# 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 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.
# Acetylcholine biosensor (Fig 5) files have modality 'ophys' and cellType 'pan'
# (pan-cellular GRABACh3.0 imaging, not cell-specific).
# State token values: CTRL | PD | OFF | ON (and their pretty names per the paper).

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_f5_achfp(asset_path):
    """Figure 5 acetylcholine biosensor: ophys modality with pan cellType."""
    fields = parse_session_id(get_session_id(asset_path))
    if not fields: return False
    return get_modality(asset_path) == 'ophys' and fields['cell_type'] == 'pan'

def get_condition_label(session_id):
    """Map state token to paper condition label."""
    fields = parse_session_id(session_id)
    if not fields: return 'unknown'
    return {
        'CTRL': 'UL control',
        'PD':   '6-OHDA',
        'OFF':  'LID off-state',
        'ON':   'LID on-state',
    }.get(fields['state'], fields['state'])

def process_acetylcholine_trial(
    timestamps: np.ndarray,
    fluorescence: np.ndarray,
    stim_time: float = 10.0,
) -> dict:
    """Process a single acetylcholine biosensor trial."""
    baseline_window = (stim_time - 1.0, stim_time)
    response_window = (stim_time, stim_time + 2.0)

    time_interval = timestamps[1] - timestamps[0]
    idx_baseline_start = np.argmin(np.abs(timestamps - baseline_window[0]))
    idx_baseline_end = np.argmin(np.abs(timestamps - baseline_window[1]))
    idx_response_start = np.argmin(np.abs(timestamps - response_window[0]))
    idx_response_end = np.argmin(np.abs(timestamps - response_window[1]))

    F0 = float(np.mean(fluorescence[idx_baseline_start:idx_baseline_end]))

    dF_over_F0 = (fluorescence - F0) / F0 if F0 > 0 else np.zeros_like(fluorescence)

    auc = float(np.sum(dF_over_F0[idx_response_start:idx_response_end]) * time_interval)
    peak_response = float(np.max(dF_over_F0[idx_response_start:idx_response_end]))

    return {
        "F0": F0,
        "dF_over_F0": dF_over_F0,
        "auc": auc,
        "peak_response": peak_response,
        "timestamps": timestamps,
        "fluorescence": fluorescence,
    }

def get_calibration_values(trials_data: list[dict]) -> tuple[float, float]:
    """Extract Fmax and Fmin using full-trace calibration averages with background subtraction."""
    ach_trials = [t for t in trials_data if t["treatment"] == "ACh_calibration"]
    ttx_trials = [t for t in trials_data if t["treatment"] == "TTX_calibration"]

    if ach_trials and ttx_trials:
        bg_pmt = 162.0
        ach_vals = []
        for t in ach_trials:
            ach_vals.extend((t["fluorescence"] - bg_pmt).tolist())
        ttx_vals = []
        for t in ttx_trials:
            ttx_vals.extend((t["fluorescence"] - bg_pmt).tolist())
        if ach_vals and ttx_vals:
            return float(np.mean(ach_vals)), float(np.mean(ttx_vals))

    # Fallback constants (from original script)
    return 493.47, 212.39

def normalize_fluorescence(trial_data: dict, Fmax: float, Fmin: float, stim_time: float = 10.0) -> dict:
    """Normalize fluorescence; compute AUC via trapezoidal integration in [stim, stim+2s]."""
    FI = Fmax - Fmin
    dF_over_FI = (trial_data["fluorescence"] - trial_data["F0"]) / FI if FI > 0 else np.zeros_like(trial_data["fluorescence"])
    trial_data["dF_over_FI"] = dF_over_FI

    idx_start = np.argmin(np.abs(trial_data["timestamps"] - stim_time))
    idx_end = np.argmin(np.abs(trial_data["timestamps"] - (stim_time + 2.0)))
    rs = dF_over_FI[idx_start:idx_end]
    ts = trial_data["timestamps"][idx_start:idx_end]

    if len(rs) > 1:
        trial_data["auc_normalized"] = float(np.trapezoid(rs, ts))
    else:
        trial_data["auc_normalized"] = 0.0

    trial_data["peak_normalized"] = float(np.max(rs)) if len(rs) > 0 else 0.0
    return trial_data

print("Analysis functions defined (trapz + full-trace calibration)")

Analysis functions defined (trapz + full-trace calibration)

from collections import Counter

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

f5_achfp_assets = [a for a in assets_list if is_f5_achfp(a.path)]
print(f'Found {len(f5_achfp_assets)} F5 acetylcholine biosensor files')

cohort = Counter(get_condition_label(get_session_id(a.path)) for a in f5_achfp_assets)
print('\nBreakdown by condition:')
for cond, n in sorted(cohort.items()):
    print(f'  {cond}: {n} files')

Found 11 F5 acetylcholine biosensor files

Breakdown by condition:
  6-OHDA: 3 files
  LID off-state: 5 files
  UL control: 3 files

# Initialize data collection
all_trial_data = []  # All individual trials for analysis
stimulation_type = "single_pulse"  # Focus on single pulse stimulation

print(f"Processing Figure 5 acetylcholine biosensor files for {stimulation_type} stimulation...\n")

for i, asset in enumerate(tqdm(f5_achfp_assets, desc="Processing F5 AChFP files")):
    session_id = get_session_id(asset.path)
    condition = get_condition_label(session_id)

    tqdm.write(f"  {i+1}/{len(f5_achfp_assets)}: {condition} - {session_id}")

    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()

    trials_df = nwbfile.trials.to_dataframe()
    fluorescence_module = nwbfile.processing["ophys"]["Fluorescence"]

    trials_data = []
    for _, trial in trials_df.iterrows():
        series_name = trial['roi_series_name']
        series = fluorescence_module.roi_response_series[series_name]

        if "GRABCh" not in series_name:
            continue

        treatment_mapping = {
            'control': 'control',
            '50nM_dopamine': 'dopamine',
            'quinpirole': 'quinpirole',
            'sulpiride': 'sulpiride',
            'TTX_calibration': 'TTX_calibration',
            'acetylcholine_calibration': 'ACh_calibration'
        }

        treatment = treatment_mapping.get(trial['treatment'], trial['treatment'])
        stimulation = 'calibration' if 'calibration' in treatment else trial['stimulation']

        # Series data shape may be (T, n_roi); squeeze to 1D for single-ROI analysis
        fluorescence = np.squeeze(series.data[:])
        if fluorescence.ndim != 1:
            fluorescence = fluorescence[:, 0]  # take first ROI if still 2D
        timestamps = series.get_timestamps()
        timestamps = timestamps - timestamps[0]

        abs_stim_time = trial['stimulus_start_time']
        orig_start = series.get_timestamps()[0]
        stim_time = abs_stim_time - orig_start

        trials_data.append({
            "condition": condition,
            "treatment": treatment,
            "stimulation": stimulation,
            "subject_id": nwbfile.subject.subject_id,
            "session_id": nwbfile.identifier,
            "series_name": series_name,
            "fluorescence": fluorescence,
            "timestamps": timestamps,
            "stim_time": stim_time,
            "file": asset.path.split('/')[-1]
        })

    Fmax, Fmin = get_calibration_values(trials_data)

    for trial in trials_data:
        if trial["treatment"] in ["ACh_calibration", "TTX_calibration"]:
            continue
        if trial["stimulation"] != stimulation_type:
            continue

        trial_data = process_acetylcholine_trial(trial["timestamps"], trial["fluorescence"], trial["stim_time"])
        trial_data = normalize_fluorescence(trial_data, Fmax, Fmin, trial["stim_time"])

        all_trial_data.append({
            "condition": trial["condition"],
            "treatment": trial["treatment"],
            "subject_id": trial["subject_id"],
            "session_id": trial["session_id"],
            "series_name": trial["series_name"],
            "F0": trial_data["F0"],
            "auc_unnormalized": trial_data["auc"],
            "auc_normalized": trial_data["auc_normalized"],
            "peak_unnormalized": trial_data["peak_response"],
            "peak_normalized": trial_data["peak_normalized"],
            "file": trial["file"],
        })

    tqdm.write(f"    Processed {len(trials_data)} trials")

    io.close()
    file.close()

df = pd.DataFrame(all_trial_data)

print(f"\nData processing complete:")
print(f"  Total measurements: {len(df)}")
print(f"  Experimental conditions: {df['condition'].nunique()}")

print("\nMeasurement breakdown by condition and treatment:")
for condition in df['condition'].unique():
    condition_data = df[df['condition'] == condition]
    print(f"  {condition}: {len(condition_data)} total measurements")
    for treatment in condition_data['treatment'].unique():
        treatment_data = condition_data[condition_data['treatment'] == treatment]
        print(f"    {treatment}: n={len(treatment_data)}")
    print()

Processing Figure 5 acetylcholine biosensor files for single_pulse stimulation...

Processing F5 AChFP files:   0%|          | 0/11 [00:00<?, ?it/s]

Processing F5 AChFP files:   0%|          | 0/11 [00:00<?, ?it/s]

  1/11: UL control - pan++CTRL++none++WT++20240402

Processing F5 AChFP files:   0%|          | 0/11 [00:04<?, ?it/s]

Processing F5 AChFP files:   9%|▉         | 1/11 [00:04<00:41,  4.17s/it]

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

# Collapse to per-file averages (first two trials per file×treatment)

plot_df = df[df["treatment"].isin(["control", "dopamine", "quinpirole",
"sulpiride"])].copy()
plot_df = plot_df.dropna(subset=["auc_normalized"]).reset_index().rename(columns={"index":
"_row_order"})

# Avoid GroupBy.apply deprecation by sorting then using head(2)

plot_df_sorted = plot_df.sort_values(["condition", "treatment", "file", "_row_order"])
first_two = (
    plot_df_sorted.groupby(["condition", "treatment", "file"], sort=False,
as_index=False)
    .head(2)
)

collapsed = (
    first_two.groupby(["condition", "treatment", "file"], as_index=False)
["auc_normalized"].mean()
)

# Create Figure 5F style plot (no function)

fig, ax = plt.subplots(1, 1, figsize=(10, 6))

# Conditions and labels (top)

conditions = ["UL control", "6-OHDA", "LID off-state"]
condition_labels = ["control", "6-OHDA", "LID off-state"]

# Treatments by condition (bottom)

treatments_by_condition = {
    "UL control": ["control", "quinpirole", "sulpiride"],
    "6-OHDA": ["control", "dopamine", "quinpirole", "sulpiride"],
    "LID off-state": ["control", "dopamine", "quinpirole", "sulpiride"],
}

# Display labels for treatments

display_label = {
    "control": "control",
    "dopamine": "+DA",
    "quinpirole": "+quinpirole",
    "sulpiride": "+sulpiride",
}

# Colors to match paper: black (UL control), red (6-OHDA), blue (LID off)

condition_colors = ["black", "red", "blue"]

# Build dynamic x-positions per group

group_gap = 1
current_pos = 1
x_positions = []
x_labels = []
group_bounds = []  # (start_pos, end_pos) per condition

all_data = []
all_colors = []

for i, condition in enumerate(conditions):
    cond_df = collapsed[collapsed["condition"] == condition]
    start_pos = current_pos
    treatments = treatments_by_condition[condition]

    for t in treatments:
        series = cond_df[cond_df["treatment"] == t]["auc_normalized"]
        all_data.append(series.values if not series.empty else np.array([]))
        all_colors.append(condition_colors[i])
        x_positions.append(current_pos)
        x_labels.append(display_label[t])
        current_pos += 1

    end_pos = current_pos - 1
    group_bounds.append((start_pos, end_pos))
    current_pos += group_gap  # gap before next condition

# Create box plots

box_plot = ax.boxplot(
    all_data,
    positions=x_positions,
    patch_artist=True,
    boxprops=dict(linewidth=1.5),
    whiskerprops=dict(linewidth=1.5),
    capprops=dict(linewidth=1.5),
    medianprops=dict(color="black", linewidth=2),
    flierprops=dict(marker="o", markersize=3, markeredgecolor="black", alpha=0.7),
    widths=0.6,
)

# Color boxes and add points

for i, (patch, color) in enumerate(zip(box_plot["boxes"], all_colors)):
    if color == "black":
        patch.set_facecolor("white")
    elif color == "red":
        patch.set_facecolor("#ffcccc")  # Light red
    else:
        patch.set_facecolor("#ccccff")  # Light blue
    patch.set_edgecolor(color)
    patch.set_linewidth(1.5)

    if len(all_data[i]) > 0:
        x_pos = np.random.normal(x_positions[i], 0.05, len(all_data[i]))
        ax.scatter(x_pos, all_data[i], color=color, alpha=0.7, s=20, zorder=3)

# Y-axis styling

ax.set_ylabel("AUC (normalized ΔF/F0*s)", fontsize=12, fontweight="bold")
ax.set_ylim(-0.05, 0.7)
ax.set_yticks([0, 0.2, 0.4, 0.6])

# Dotted baseline at y=0

ax.axhline(y=0, color="black", linestyle=":", alpha=0.8, linewidth=1.0)

# X-limits

first_pos = group_bounds[0][0]
last_pos = group_bounds[-1][1]
ax.set_xlim(first_pos - 0.5, last_pos + 0.5)

# Bottom: treatment labels

ax.set_xticks(x_positions)
ax.set_xticklabels(x_labels, rotation=45, ha="right", fontsize=9)

# Top: condition labels centered above groups

group_centers = [0.5 * (s + e) for (s, e) in group_bounds]
ax_top = ax.twiny()
ax_top.set_xlim(ax.get_xlim())
ax_top.set_xticks(group_centers)
ax_top.set_xticklabels(condition_labels, fontsize=11, fontweight="bold")
ax_top.tick_params(axis="x", which="major", pad=0)
for spine in ax_top.spines.values():
    spine.set_visible(False)

# Separators and shading

for (start, end) in group_bounds[:-1]:
    ax.axvline(x=end + 0.5, color="lightgray", linestyle="-", alpha=0.5, linewidth=1)

shade_colors = ["lightgray", "lightcoral", "lightblue"]
for (start, end), shade in zip(group_bounds, shade_colors):
    ax.axvspan(start - 0.5, end + 0.5, alpha=0.1, color=shade, zorder=0)

# Final styling

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)

plt.tight_layout()
plt.show()

# Summary (collapsed per-file, control only)

print("=== FIGURE 5F: ACETYLCHOLINE RELEASE ANALYSIS ===")
control_collapsed = collapsed[collapsed["treatment"] == "control"]
print("\nSummary by Condition (Control Treatment Only):")
for condition in ["UL control", "6-OHDA", "LID off-state"]:
    cd = control_collapsed[control_collapsed["condition"] == condition]
    auc_vals = cd["auc_normalized"]
    files = cd["file"].nunique()
    print(f"  {condition}: n={len(auc_vals)} measurements from {files} files")
    if len(auc_vals) > 0:
        print(f"    Mean ± SEM: {auc_vals.mean():.3f} ± {auc_vals.sem():.3f}")
        print(f"    Median: {auc_vals.median():.3f}")
        print(f"    Range: {auc_vals.min():.3f} - {auc_vals.max():.3f}")

print(f"\nTotal per-file datapoints (all treatments): {len(collapsed)}")
print(f"Conditions analyzed: {collapsed['condition'].nunique()}")
print(f"Unique files processed: {collapsed['file'].nunique()}")

=== FIGURE 5F: ACETYLCHOLINE RELEASE ANALYSIS ===

Summary by Condition (Control Treatment Only):
  UL control: n=3 measurements from 3 files
    Mean ± SEM: 0.071 ± 0.006
    Median: 0.074
    Range: 0.061 - 0.080
  6-OHDA: n=3 measurements from 3 files
    Mean ± SEM: 0.318 ± 0.147
    Median: 0.188
    Range: 0.154 - 0.611
  LID off-state: n=5 measurements from 5 files
    Mean ± SEM: 0.356 ± 0.055
    Median: 0.316
    Range: 0.200 - 0.482

Total per-file datapoints (all treatments): 41
Conditions analyzed: 3
Unique files processed: 11

Installing requirements¶

Reproducing Figure 5F: GRABACh3.0 Acetylcholine Biosensor Analysis¶

Setup and Data Loading¶

Import Libraries and Configure Plotting Style¶

Session ID Parsing and Filtering Functions¶

Fluorescence Analysis Functions¶

GRABACh3.0 Signal Processing¶

Load DANDI Dataset¶

Data Processing and Analysis¶

Process All NWB Files¶

Figure 5F: Box Plot Visualization¶

Acetylcholine Release Dynamics Across Conditions¶

Summary¶

Key Findings¶

Methodological Notes¶

Biological Significance¶