OSCR

HCN1 is a primary HCN Pacemaker Channel in Neurons.

Code ↔ Paper

3 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 3 matches
  1. [1] § Methods › Simulations of the effects of HCN1 channels in an SCN neuron ↔ HCN_rest-singlePo.py, lines 77–187 · score 0.71 · HCN1 open probability, Brian2 simulation, neuron, leak, trace
  2. [2] § Methods › Simulations of the effects of HCN1 channels in an SCN neuron ↔ HCN_rest-singlePo.py, lines 26–75 · score 0.71 · potassium reversal potential, steady state, membrane potential, dt, HCN1
  3. [3] § Methods › Simulations of the effects of HCN1 channels in an SCN neuron ↔ HCN_rest-singlePo.py, lines 77–187 · score 0.57 · voltage ramp, membrane potential, duration, Simulations, HCN, injected

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 905 lines · 36 KB · no license · 3 matches

  1. # Copyright (c) 2026 Christoph Schmidt-Hieber
  2. # SPDX-License-Identifier: MIT
  3. # Simulation of HCN channel effects on membrane potential during voltage ramps
  4. # Code accompanying Enke et al., 2026
  5. # v2026-01-30
  6. # Simplified code to compute capacitance from specified membrane time constant (tau_m)
  7. # v2026-01-21
  8. # LLM assistance was used to:
  9. # - document function arguments and return values
  10. # - assist with plotting
  11. # - assist with exporting the data to Excel
  12. # - produce data structures for the Excel export
  13. # - produce data structures to store the results of the simulations
  14. # - make the code more readable
  16. from brian2 import *
  17. import numpy as np
  18. import pandas as pd
  19. import matplotlib.pyplot as plt
  20. import matplotlib.gridspec as gridspec
  21. from scipy.interpolate import interp1d
  22. from itertools import product
  23. from pathlib import Path
  24. import pickle
  25. print("\n".join([f"{name}: {version}" for name, version in [(name, globals()[name].__version__ if hasattr(globals()[name], '__version__') else 'N/A') for name in ['brian2', 'numpy', 'pandas', 'matplotlib', 'scipy', 'itertools', 'pathlib', 'pickle']] if name in globals()]))
  26. def compute_conductances(V_rest, E_K_param, R_in, scaling, p_rest_rep):
  27. """
  28. Compute leak and HCN conductances from the biophysical parameters.
  30. Args:
  31. V_rest: Resting membrane potential (quantity).
  32. E_K_param: Potassium reversal potential (quantity).
  33. R_in: Input resistance (quantity).
  34. scaling: Dimensionless scaling factor for HCNX.
  35. p_rest_rep: Resting open probability (float).
  37. Returns:
  38. (g_leak, g_HCN_full, g_HCNX, G_total) as siemens quantities.
  39. """
  40. # Remove units and convert to float
  41. V_r = float(V_rest / volt)
  42. E_K_val = float(E_K_param / volt)
  43. E_HCN_val = float(E_HCN / volt)
  44. G_total = float(1 / R_in / siemens)
  46. # At the beginning of the ramp we assume that the system is at
  47. # a steady state and therefore we have two conditions:
  48. # 1) dV/dt is zero, so the net sum of currents is zero.
  49. # 2) Membrane potential is at V_rest
  50. # The total conductance of the cell is given by G_total = 1/R_in
  51. # and composed of g_leak + g_HCN1_rest + g_HCNX, where:
  52. # g_HCN1_rest = p_rest_rep * g_HCN_full
  53. # g_HCNX = scaling * g_HCN_full
  54. # Therefore:
  55. # g_leak + (p_rest_rep + scaling) * g_HCN_full = G_total [eq 1]
  56. # Net sum of currents at steady state is zero, therefore:
  57. # g_leak * (V_rest - E_K) + g_HCN1_rest * (V_rest - E_HCN) + g_HCNX * (V_rest - E_HCN) = 0
  58. # We can ignore capacitive currents because dV/dt is zero at steady state, hence C*dV/dt = 0.
  59. # Rearranging yields:
  60. # g_leak * (V_rest - E_K) + (p_rest_rep + scaling) * g_HCN_full * (V_rest - E_HCN) = 0 [eq 2]
  61. # Solving equation 1 for g_leak: g_leak = G_total - (p_rest_rep + scaling) * g_HCN_full,
  62. # and substituting into equation 2 yields:
  63. # G_total * (V_rest - E_K) + (p_rest_rep + scaling) * g_HCN_full * (E_K - E_HCN) = 0
  64. # Rearranging to isolate g_HCN_full gives:
  65. # G_total * (V_rest - E_K) = (p_rest_rep + scaling) * g_HCN_full * (E_HCN - E_K)
  66. # Solving for g_HCN_full yields:
  67. # g_HCN_full = G_total * (V_rest - E_K) / ((p_rest_rep + scaling) * (E_HCN_val - E_K_val))
  68. denom = (p_rest_rep + scaling) * (E_HCN_val - E_K_val)
  69. g_HCN_full = G_total * (V_r - E_K_val) / denom
  71. g_HCNX = scaling * g_HCN_full
  72. g_HCN1_rest = p_rest_rep * g_HCN_full
  73. g_leak = G_total - g_HCN1_rest - g_HCNX
  75. return g_leak * siemens, g_HCN_full * siemens, g_HCNX * siemens, G_total * siemens
  77. def run_sim(V_rest, E_K_param, R_in, scaling,
  78. P_HCN1_timed, t_interp, sim_duration, t_start_ramp, dVdt_target,
  79. p_rest_rep,
  80. inject_ramp=False, label="sim", rep_idx=None):
  81. """
  82. Run a single Brian2 simulation and store summary/trace outputs.
  84. Args:
  85. V_rest: Resting membrane potential (quantity).
  86. E_K_param: Potassium reversal potential (quantity).
  87. R_in: Input resistance (quantity).
  88. scaling: Dimensionless scaling factor for HCNX.
  89. P_HCN1_timed: TimedArray of HCN1 open probability.
  90. t_interp: Time vector with units for ramp/current construction.
  91. sim_duration: Total simulation duration (quantity).
  92. t_start_ramp: Time when the ramp starts (quantity).
  93. dVdt_target: Target voltage ramp slope (quantity).
  94. p_rest_rep: Resting open probability (float).
  95. inject_ramp: Whether to inject the compensatory ramp current.
  96. label: Key for storing trace outputs.
  97. rep_idx: Replicate identifier for summary output.
  98. """
  99. start_scope()
  101. g_leak_run, g_HCN_full_run, g_HCNX_run, G_total_run = compute_conductances(
  102. V_rest, E_K_param, R_in, scaling, p_rest_rep
  103. )
  104. # Compute total membrane capacitance from specified membrane time constant:
  105. # C_total = tau_m / R_in (since tau = R * C -> C = tau / R)
  106. C_total = tau_m / R_in
  108. if inject_ramp:
  109. V_ramp = V_start + dVdt_target * t_interp
  110. n_hold = int(hold_duration * ms / defaultclock.dt)
  112. I_hold = np.zeros(n_hold) * amp
  114. t_ramp = t_interp[n_hold:]
  115. V_ramp_dyn = V_ramp[n_hold:]
  116. P_HCN_dyn = P_HCN1_timed(t_ramp)
  118. I_ramp = -(
  119. C_total * dVdt_target
  120. - g_leak_run * (V_ramp_dyn - E_K_param)
  121. - g_HCN_full_run * P_HCN_dyn * (V_ramp_dyn - E_HCN)
  122. - g_HCNX_run * (V_ramp_dyn - E_HCN)
  123. ) # Quantity in amp
  125. # --- FIX: concatenate unitless magnitudes, then reattach units ---
  126. I_inj_array = np.concatenate([
  127. np.asarray(I_hold / amp),
  128. np.asarray(I_ramp / amp),
  129. ]) * amp
  131. else:
  132. I_inj_array = np.zeros(len(t_interp)) * amp
  134. I_inj_timed = TimedArray(I_inj_array, dt=defaultclock.dt)
  136. eqs = '''
  137. dv/dt = (I_leak + I_HCN1 + I_HCNX + I_inj) / C : volt
  138. I_leak = -g_leak * (v - E_K) : amp
  139. I_HCN1 = -g_HCN_full * P_HCN1_timed(t) * (v - E_HCN) : amp
  140. I_HCNX = -g_HCNX * (v - E_HCN) : amp
  141. I_inj = I_inj_timed(t) : amp
  142. C : farad (constant)
  143. E_K : volt (constant)
  144. g_leak : siemens (constant)
  145. g_HCN_full : siemens (constant)
  146. g_HCNX : siemens (constant)
  147. '''
  149. neuron = NeuronGroup(1, eqs, method='euler')
  150. neuron.v = V_rest
  151. neuron.C = C_total
  152. neuron.E_K = E_K_param
  153. neuron.g_leak = g_leak_run
  154. neuron.g_HCN_full = g_HCN_full_run
  155. neuron.g_HCNX = g_HCNX_run
  157. mon = StateMonitor(neuron, ['v', 'I_HCN1', 'I_leak', 'I_inj'], record=0)
  158. run(sim_duration)
  160. time = mon.t / ms
  161. v = mon.v[0] / mV
  162. i_hcn1 = mon.I_HCN1[0] / nA
  163. i_leak = mon.I_leak[0] / nA
  164. i_inj = mon.I_inj[0] / nA
  165. p_open = P_HCN1_timed(mon.t)
  167. i_start = int(t_start_ramp / defaultclock.dt)
  168. V_max = float(np.max(v[i_start:]))
  170. summary_data.append({
  171. "replicate": rep_idx,
  172. "E_K_mV": float(E_K_param / mV),
  173. "R_in_Mohm": float(R_in / ohm) / 1e6,
  174. "scaling": float(scaling),
  175. "V_max_mV": V_max,
  176. "injected": bool(inject_ramp),
  177. "label": label
  178. })
  180. time_series_data[label] = {
  181. "time_ms": np.asarray(time),
  182. "v_mV": np.asarray(v),
  183. "i_hcn1_nA": np.asarray(i_hcn1),
  184. "i_leak_nA": np.asarray(i_leak),
  185. "i_inj_nA": np.asarray(i_inj),
  186. "p_open": np.asarray(p_open)
  187. }
  189. def closest(val, arr):
  190. """
  191. Return the element of arr closest to val.
  193. Args:
  194. val: Target value (float-like).
  195. arr: Iterable of numeric values.
  196. """
  197. arr = np.asarray(list(arr), dtype=float)
  198. return float(arr[np.argmin(np.abs(arr - float(val)))])
  200. def labels_for_noinj_condition(summary_df, E_K_mV, R_in_Mohm, scaling):
  201. """
  202. Get labels matching a non-injected condition in summary_df.
  204. Args:
  205. summary_df: DataFrame with summary data and labels.
  206. E_K_mV: Potassium reversal potential in mV (float).
  207. R_in_Mohm: Input resistance in MΩ (float).
  208. scaling: Scaling value (float).
  209. """
  210. m = (
  211. (summary_df["injected"] == False) &
  212. (summary_df["E_K_mV"] == float(E_K_mV)) &
  213. (summary_df["R_in_Mohm"] == float(R_in_Mohm)) &
  214. (summary_df["scaling"] == float(scaling))
  215. )
  216. labels = summary_df.loc[m, "label"].dropna().tolist()
  217. labels = [l for l in labels if l in time_series_data]
  218. return labels
  220. def labels_for_injected(summary_df, time_series_data):
  221. """
  222. Get labels for injected traces that exist in time_series_data.
  224. Args:
  225. summary_df: DataFrame with summary data and labels.
  226. time_series_data: Dict of trace data keyed by label.
  227. """
  228. labels = summary_df.loc[summary_df["injected"] == True, "label"].dropna().tolist()
  229. labels = [l for l in labels if l in time_series_data]
  230. if len(labels) > 0:
  231. return labels
  232. raise RuntimeError("No injected traces found in summary_df['label'] that exist in time_series_data.")
  234. def mean_trace(time_series_data, labels, key):
  235. """
  236. Compute mean trace by interpolating all traces onto the shortest common time grid
  237. over the overlapping time window.
  238. Returns (t_ref, y_mean).
  240. Args:
  241. time_series_data: Dict of trace data keyed by label.
  242. labels: Labels to include in the mean.
  243. key: Trace field to average (e.g., "v_mV").
  244. """
  245. t_list = [np.asarray(time_series_data[l]["time_ms"], dtype=float) for l in labels]
  246. t_min = max(t[0] for t in t_list)
  247. t_max = min(t[-1] for t in t_list)
  248. if not (t_max > t_min):
  249. raise RuntimeError("No overlapping time window across traces for averaging.")
  251. masks = [(t >= t_min) & (t <= t_max) for t in t_list]
  252. lengths = [int(np.sum(m)) for m in masks]
  253. idx_ref = int(np.argmin(lengths))
  254. t_ref = t_list[idx_ref][masks[idx_ref]]
  256. Y = []
  257. for lbl, t, m in zip(labels, t_list, masks):
  258. y = np.asarray(time_series_data[lbl][key], dtype=float)
  259. Y.append(np.interp(t_ref, t[m], y[m]))
  260. Y = np.vstack(Y)
  261. return t_ref, np.mean(Y, axis=0)
  263. def plot_cloud_and_mean(ax, labels, xkey, ykey,
  264. color="#000000", alpha_cloud=0.25, lw_mean=2.0,
  265. y_transform=None, label_mean=None):
  266. """
  267. Plot all individual traces with the SAME color as the mean (transparent),
  268. then plot the mean trace on top (opaque).
  270. Args:
  271. ax: Matplotlib axis to plot on.
  272. labels: Labels of traces to plot.
  273. xkey: X-axis field name (e.g., "time_ms").
  274. ykey: Y-axis field name (e.g., "v_mV").
  275. color: Trace/mean color.
  276. alpha_cloud: Alpha for individual traces.
  277. lw_mean: Line width for mean trace.
  278. y_transform: Optional transform applied to y-values.
  279. label_mean: Legend label for the mean trace.
  280. """
  281. # Cloud
  282. for lbl in labels:
  283. ts = time_series_data[lbl]
  284. x = np.asarray(ts[xkey], dtype=float)
  285. y = np.asarray(ts[ykey], dtype=float)
  286. if y_transform is not None:
  287. y = y_transform(y)
  288. ax.plot(x, y, color=color, alpha=alpha_cloud, linewidth=1.0)
  290. # Mean
  291. t_mean, y_mean = mean_trace(time_series_data, labels, ykey)
  292. if y_transform is not None:
  293. y_mean = y_transform(y_mean)
  294. ax.plot(t_mean, y_mean, color=color, alpha=1.0, linewidth=lw_mean,
  295. label=label_mean)
  297. def plot_dist(ax, x_vals, df_subset, x_col, xlabel, ylabel=None, semilogx=False,
  298. point_size=22, point_alpha=0.55, color=None):
  299. """
  300. Mean±SD and replicate points. Color optionally set (e.g. to match panel b by E_K).
  302. Args:
  303. ax: Matplotlib axis to plot on.
  304. x_vals: Ordered x-axis values.
  305. df_subset: DataFrame subset with summary values.
  306. x_col: Column name for x values in df_subset.
  307. xlabel: X-axis label.
  308. ylabel: Y-axis label (or None to skip).
  309. semilogx: Whether to use a log-scaled x-axis.
  310. point_size: Marker size for replicate points.
  311. point_alpha: Alpha for replicate points.
  312. color: Color for mean±SD and replicate points. If None, use black / gray.
  313. """
  314. line_color = color if color is not None else "black"
  315. point_color = color if color is not None else "0.35"
  317. g = (df_subset.groupby(x_col, as_index=False)
  318. .agg(V_mean=("V_max_mV", "mean"),
  319. V_sd=("V_max_mV", "std"),
  320. n=("V_max_mV", "count")))
  321. g = g.set_index(x_col).reindex(x_vals).reset_index()
  323. x = g[x_col].to_numpy(dtype=float)
  324. y = g["V_mean"].to_numpy(dtype=float)
  325. yerr = g["V_sd"].to_numpy(dtype=float)
  327. if semilogx:
  328. ax.errorbar(x, y, yerr=yerr, fmt='o-', capsize=3, color=line_color,
  329. markerfacecolor=line_color, markeredgecolor=line_color, linewidth=1.6)
  330. ax.set_xscale("log")
  331. else:
  332. ax.errorbar(x, y, yerr=yerr, fmt='o-', capsize=3, color=line_color,
  333. markerfacecolor=line_color, markeredgecolor=line_color, linewidth=1.6)
  335. # Replicate points: slightly jitter only for semilogx readability; otherwise no jitter is fine here
  336. reps = df_subset.copy()
  337. if semilogx:
  338. jitter = np.exp(np.random.normal(loc=0.0, scale=0.03, size=len(reps)))
  339. x_rep = reps[x_col].to_numpy(dtype=float) * jitter
  340. else:
  341. x_rep = reps[x_col].to_numpy(dtype=float)
  343. ax.scatter(x_rep, reps["V_max_mV"].to_numpy(dtype=float),
  344. s=point_size, color=point_color, alpha=point_alpha, edgecolors="none")
  346. ax.set_xlabel(xlabel)
  347. if ylabel is not None:
  348. ax.set_ylabel(ylabel)
  350. def plot_scaling_curves_by_EK(ax, df, rin_fixed, ek_values, scaling_order,
  351. ek_color_map,
  352. xlabel="HCN2-4 scaling", ylabel=None,
  353. point_size=28, point_alpha=0.65):
  354. """
  355. Vmax vs scaling at fixed R_in, with one curve per E_K.
  356. Replicate points:
  357. - colored by E_K
  358. - plotted at exact scaling x positions, with a tiny EK-dependent x-offset
  359. so different EK groups at the same scaling don't overlap.
  361. Args:
  362. ax: Matplotlib axis to plot on.
  363. df: DataFrame with summary data.
  364. rin_fixed: Fixed R_in value (MΩ).
  365. ek_values: Ordered E_K values (mV).
  366. scaling_order: Ordered scaling values.
  367. ek_color_map: Mapping from E_K to color.
  368. xlabel: X-axis label.
  369. ylabel: Y-axis label (or None to skip).
  370. point_size: Marker size for replicate points.
  371. point_alpha: Alpha for replicate points.
  372. """
  373. df0 = df[df["R_in_Mohm"] == float(rin_fixed)].copy()
  374. if df0.empty:
  375. raise RuntimeError(f"No data for R_in_Mohm == {rin_fixed}")
  377. # Tiny, fixed offsets per EK (in x-axis units of 'scaling')
  378. # Chosen to be visually separable but negligible relative to tick spacing.
  379. offsets = np.linspace(-0.035, 0.035, num=max(1, len(ek_values)))
  380. ek_to_dx = {float(ek): float(dx) for ek, dx in zip(ek_values, offsets)}
  382. for ek in ek_values:
  383. d = df0[df0["E_K_mV"] == float(ek)].copy()
  384. if d.empty:
  385. continue
  387. col = ek_color_map[float(ek)]
  388. dx = ek_to_dx[float(ek)]
  390. g = (d.groupby("scaling", as_index=False)
  391. .agg(V_mean=("V_max_mV", "mean"),
  392. V_sd=("V_max_mV", "std"),
  393. n=("V_max_mV", "count")))
  394. g = g.set_index("scaling").reindex(scaling_order).reset_index()
  396. x = g["scaling"].to_numpy(float) + dx
  397. y = g["V_mean"].to_numpy(float)
  398. yerr = g["V_sd"].to_numpy(float)
  400. # Mean ± SD (also offset slightly to match the replicate cloud)
  401. ax.errorbar(x, y, yerr=yerr, fmt='o-', capsize=3, linewidth=1.8,
  402. color=col, markerfacecolor=col, markeredgecolor=col,
  403. label=rf"$E_{{\mathrm{{K}}}}={ek:.1f}\,\mathrm{{mV}}$")
  405. # Replicate points: same scaling within EK, but EK-dependent dx
  406. ax.scatter(d["scaling"].to_numpy(float) + dx,
  407. d["V_max_mV"].to_numpy(float),
  408. s=point_size, color=col, alpha=point_alpha, edgecolors="none")
  410. ax.set_xlabel(xlabel)
  411. if ylabel is not None:
  412. ax.set_ylabel(ylabel)
  414. # === Constants ===
  415. # Membrane time constant (tau_m) used to compute total capacitance C_total = tau_m / R_in
  416. tau_m = 30 * ms
  417. E_HCN = 0 * mV
  418. V_start = -81 * mV
  419. V_end = -53 * mV
  420. hold_duration = 50 # ms
  421. sim_dt = 10.0 # µs
  422. defaultclock.dt = sim_dt * us
  424. # === Plotting config ===
  425. # Fixed parameter defaults used in plotting and summaries
  426. EK_DEFAULT = -95.0 # mV
  427. RIN_DEFAULT = 100.0 # MΩ
  428. SC_DEFAULT = 0.0
  430. # Colors (Matplotlib default cycle colors by explicit hex)
  431. COL_V = "#1f77b4" # blue
  432. COL_P = "#ff7f0e" # orange
  433. COL_HCN = "#2ca02c" # green
  434. COL_LEAK = "#d62728" # red
  435. COL_INJ = "#9467bd" # purple
  437. # For panel b (curves by EK): use 3 distinct colors, consistent for mean+replicates
  438. EK_COLOR_MAP = None # filled after we know EK values
  440. # === Simulation cache ===
  441. USE_SIM_CACHE = True
  442. SIM_CACHE_PATH = Path("sim_cache.pkl")
  444. # === Load HCN1 open probability data (time + 11 p-columns) ===
  445. file_name = "mHCN1_single_20260203.xlsx"
  446. df = pd.read_excel(file_name)
  448. iend = 1338
  449. # time column (shared)
  450. time_col = pd.to_numeric(df.iloc[:iend, 0], errors="coerce")
  452. replicates = []
  453. eps = 0.0 # threshold to define "non-zero"
  454. N_rest = 10 # first N non-zero points used to estimate p_rest per replicate
  456. for col_idx in range(1, 12): # columns 1..11
  457. p_col = pd.to_numeric(df.iloc[:iend, col_idx], errors="coerce")
  459. valid = (~time_col.isna()) & (~p_col.isna())
  460. t_vals = time_col[valid].to_numpy(dtype=float)
  461. p_vals = p_col[valid].to_numpy(dtype=float)
  463. # Replicate-specific p_rest from first N non-zero points
  464. n_take = min(N_rest, len(p_vals))
  465. p_rest_rep = float(np.mean(p_vals[:n_take])) # or np.median(p_vals[:n_take]) for robustness
  467. # Re-zero time so that dynamic starts at hold_duration cleanly
  468. t_vals = t_vals - t_vals[0]
  470. replicates.append((t_vals, p_vals, col_idx-1, p_rest_rep))
  474. # === Parameter grid derived from explicit value lists ===
  475. V_rest_values = [-81 * mV]
  476. E_K_values = [-95 * mV, -92.5 * mV, -90 * mV]
  477. R_in_values = [1e7 * ohm, 1e8 * ohm, 1e9 * ohm]
  478. scaling_values = [0.0, 0.5, 1.0, 1.5, 2.0]
  480. param_grid = [
  481. {"V_rest": V_rest, "E_K_param": E_K, "R_in": R_in, "scaling": scaling}
  482. for (V_rest, E_K, R_in, scaling) in product(V_rest_values, E_K_values, R_in_values, scaling_values)
  483. ]
  486. # === Run and store simulation results ===
  487. summary_data = []
  488. time_series_data = {}
  490. if USE_SIM_CACHE and SIM_CACHE_PATH.exists():
  491. with SIM_CACHE_PATH.open("rb") as f:
  492. cached = pickle.load(f)
  493. summary_data = cached.get("summary_data", [])
  494. time_series_data = cached.get("time_series_data", {})
  495. else:
  496. # === Run all simulations for each replicate ===
  497. for t_vals, p_vals, rep_idx, p_rest_rep in replicates:
  498. # Interpolate HCN open probability for this replicate
  499. t_hold = np.arange(0, hold_duration, sim_dt / 1000.0)
  500. p_hold = np.full_like(t_hold, p_rest_rep)
  502. t_dynamic = t_vals + hold_duration
  503. p_combined = np.concatenate((p_hold, p_vals))
  504. t_combined = np.concatenate((t_hold, t_dynamic))
  506. # float time bases in ms
  507. t_interp_ms = np.arange(0, t_combined[-1] + sim_dt/1000.0, sim_dt/1000.0)
  508. p_interp = interp1d(
  509. t_combined, p_combined, kind="linear", fill_value="extrapolate"
  510. )(t_interp_ms)
  512. # TimedArray expects values sampled at defaultclock.dt; dt already set in us
  513. P_HCN1_timed = TimedArray(p_interp, dt=defaultclock.dt)
  515. # Keep a quantity time vector for Brian2 computations that expect units
  516. t_interp = t_interp_ms * ms
  518. t_start_ramp = hold_duration * ms
  519. t_end_ramp = t_dynamic[-1] * ms
  520. ramp_duration = t_end_ramp - t_start_ramp
  521. dVdt_target = (V_end - V_start) / ramp_duration
  522. sim_duration = t_end_ramp
  524. # Baseline runs: exhaustive parameter grid
  525. for i, params in enumerate(param_grid):
  526. run_sim(
  527. **params,
  528. P_HCN1_timed=P_HCN1_timed,
  529. t_interp=t_interp,
  530. sim_duration=sim_duration,
  531. t_start_ramp=t_start_ramp,
  532. dVdt_target=dVdt_target,
  533. p_rest_rep=p_rest_rep,
  534. inject_ramp=False,
  535. label=f"rep{rep_idx}_grid{i}",
  536. rep_idx=rep_idx
  537. )
  539. # Injected ramp run
  540. run_sim(
  541. V_rest=-81 * mV, E_K_param=EK_DEFAULT * mV, R_in=1e8 * ohm, scaling=SC_DEFAULT,
  542. P_HCN1_timed=P_HCN1_timed,
  543. t_interp=t_interp,
  544. sim_duration=sim_duration,
  545. t_start_ramp=t_start_ramp,
  546. dVdt_target=dVdt_target,
  547. p_rest_rep=p_rest_rep,
  548. inject_ramp=True,
  549. label=f"rep{rep_idx}_injected",
  550. rep_idx=rep_idx
  551. )
  553. with SIM_CACHE_PATH.open("wb") as f:
  554. pickle.dump(
  555. {"summary_data": summary_data, "time_series_data": time_series_data},
  556. f,
  557. protocol=pickle.HIGHEST_PROTOCOL,
  558. )
  560. summary_df = pd.DataFrame(summary_data)
  562. # Normalize legacy column names if needed
  563. rename_map = {}
  564. if "E_K" in summary_df.columns and "E_K_mV" not in summary_df.columns:
  565. rename_map["E_K"] = "E_K_mV"
  566. if "R_in" in summary_df.columns and "R_in_Mohm" not in summary_df.columns:
  567. rename_map["R_in"] = "R_in_Mohm"
  568. if "V_max" in summary_df.columns and "V_max_mV" not in summary_df.columns:
  569. rename_map["V_max"] = "V_max_mV"
  570. summary_df = summary_df.rename(columns=rename_map)
  572. df_no_inj = summary_df[summary_df["injected"] == False].copy()
  573. df_inj = summary_df[summary_df["injected"] == True].copy()
  575. # -----------------------------------------------------------------------------
  576. # Collect labels for traces
  577. # -----------------------------------------------------------------------------
  578. # Panels a and d use (E_K, R_in, scaling) = (EK_DEFAULT, RIN_DEFAULT, SC_DEFAULT)
  579. # (= -95 mV, 100 MΩ, 0 by default).
  580. labels_noinj = labels_for_noinj_condition(summary_df, EK_DEFAULT, RIN_DEFAULT, SC_DEFAULT)
  581. if len(labels_noinj) == 0:
  582. raise RuntimeError(
  583. "Requested non-injected trace condition not found in summary_data/time_series_data:\n"
  584. f"E_K={EK_DEFAULT} mV, R_in={RIN_DEFAULT} MΩ, scaling={SC_DEFAULT}"
  585. )
  587. labels_inj = labels_for_injected(summary_df, time_series_data)
  589. # -----------------------------------------------------------------------------
  590. # Reference values for parameter sweeps (panels b and c)
  591. # -----------------------------------------------------------------------------
  592. E_K_vals = sorted(df_no_inj["E_K_mV"].unique())
  593. R_in_vals = sorted(df_no_inj["R_in_Mohm"].unique())
  594. sc_vals = sorted(df_no_inj["scaling"].unique())
  596. E_K_ref_for_rin = closest(EK_DEFAULT, E_K_vals)
  597. sc_ref_for_rin = closest(SC_DEFAULT, sc_vals)
  599. df_rin = df_no_inj[(df_no_inj["scaling"] == sc_ref_for_rin) & (df_no_inj["E_K_mV"] == E_K_ref_for_rin)].copy()
  600. rin_order = sorted(df_rin["R_in_Mohm"].unique())
  602. ek_values = sorted(df_no_inj["E_K_mV"].unique())
  603. scaling_order = sorted(df_no_inj["scaling"].unique())
  605. # Map EK->color (3 curves)
  606. # Use Matplotlib default qualitative colors for clean separation
  607. ek_palette = ["#1f77b4", "#ff7f0e", "#2ca02c"] # blue, orange, green
  608. EK_COLOR_MAP = {float(ek): ek_palette[i % len(ek_palette)] for i, ek in enumerate(ek_values)}
  610. # -----------------------------------------------------------------------------
  611. # Figure layout:
  612. # Top row: a (left) and d (right), aligned, each is a 3-row trace stack
  613. # Bottom row: b (left, scaling curves by E_K) and c (right, R_in sweep)
  614. # -----------------------------------------------------------------------------
  615. fig = plt.figure(figsize=(6, 10)) # wider, less tall (since a and c are side-by-side)
  617. gs_outer = gridspec.GridSpec(
  618. 2, 2, figure=fig,
  619. height_ratios=[3.2, 1.55],
  620. width_ratios=[1.0, 1.0],
  621. hspace=0.55, wspace=0.28
  622. )
  624. # Subgrids for trace stacks
  625. gs_a = gridspec.GridSpecFromSubplotSpec(
  626. 3, 1, subplot_spec=gs_outer[0, 0],
  627. height_ratios=[1.35, 1.15, 1.05],
  628. hspace=0.14
  629. )
  630. gs_d = gridspec.GridSpecFromSubplotSpec(
  631. 3, 1, subplot_spec=gs_outer[0, 1],
  632. height_ratios=[1.35, 1.15, 1.05],
  633. hspace=0.14
  634. )
  636. # Panels b and c in bottom row
  637. gs_bc = gridspec.GridSpecFromSubplotSpec(
  638. 1, 2, subplot_spec=gs_outer[1, :],
  639. width_ratios=[2.0, 1.0],
  640. wspace=0.35
  641. )
  643. # Axes creation
  644. ax_a1 = fig.add_subplot(gs_a[0, 0])
  645. ax_a2 = fig.add_subplot(gs_a[1, 0], sharex=ax_a1)
  646. ax_a3 = fig.add_subplot(gs_a[2, 0], sharex=ax_a1)
  648. ax_d1 = fig.add_subplot(gs_d[0, 0], sharex=ax_a1, sharey=ax_a1)
  649. ax_d2 = fig.add_subplot(gs_d[1, 0], sharex=ax_a1, sharey=ax_a2)
  650. ax_d3 = fig.add_subplot(gs_d[2, 0], sharex=ax_a1, sharey=ax_a3)
  652. ax_b1 = fig.add_subplot(gs_bc[0, 0])
  653. ax_c1 = fig.add_subplot(gs_bc[0, 1], sharey=ax_b1)
  655. fig.subplots_adjust(left=0.08, right=0.98, top=0.96, bottom=0.10)
  657. # Hide x tick labels on upper axes in stacks
  658. for ax in (ax_a1, ax_a2, ax_d1, ax_d2):
  659. ax.tick_params(bottom=False, labelbottom=False)
  661. # -----------------------------------------------------------------------------
  662. # Panel a: Non-injected traces (colored clouds + colored means)
  663. # -----------------------------------------------------------------------------
  664. plot_cloud_and_mean(ax_a1, labels_noinj, "time_ms", "p_open",
  665. color=COL_P, alpha_cloud=0.20, lw_mean=2.3)
  666. ax_a1.set_ylabel(r"$P_{\mathrm{open}}$")
  668. plot_cloud_and_mean(ax_a2, labels_noinj, "time_ms", "v_mV",
  669. color=COL_V, alpha_cloud=0.20, lw_mean=2.3)
  670. ax_a2.set_ylabel(r"$V$ (mV)")
  672. # Currents: two colors (no patterns)
  673. plot_cloud_and_mean(ax_a3, labels_noinj, "time_ms", "i_hcn1_nA",
  674. color=COL_HCN, alpha_cloud=0.18, lw_mean=2.2,
  675. y_transform=lambda y: -y, label_mean="HCN1 (mean)")
  676. plot_cloud_and_mean(ax_a3, labels_noinj, "time_ms", "i_leak_nA",
  677. color=COL_LEAK, alpha_cloud=0.18, lw_mean=2.2,
  678. y_transform=lambda y: -y, label_mean="Leak (mean)")
  680. ax_a3.set_ylabel("Current (nA)")
  681. ax_a3.set_xlabel("Time (ms)")
  682. ax_a3.legend(frameon=False, fontsize=9, loc="upper right")
  684. # -----------------------------------------------------------------------------
  685. # Panel d: Injected traces (colored clouds + colored means)
  686. # -----------------------------------------------------------------------------
  687. plot_cloud_and_mean(ax_d1, labels_inj, "time_ms", "p_open",
  688. color=COL_P, alpha_cloud=0.20, lw_mean=2.3)
  689. ax_d1.set_ylabel(r"$P_{\mathrm{open}}$")
  691. plot_cloud_and_mean(ax_d2, labels_inj, "time_ms", "v_mV",
  692. color=COL_V, alpha_cloud=0.20, lw_mean=2.3)
  693. ax_d2.set_ylabel(r"$V$ (mV)")
  695. plot_cloud_and_mean(ax_d3, labels_inj, "time_ms", "i_hcn1_nA",
  696. color=COL_HCN, alpha_cloud=0.18, lw_mean=2.2,
  697. y_transform=lambda y: -y, label_mean="HCN1 (mean)")
  698. plot_cloud_and_mean(ax_d3, labels_inj, "time_ms", "i_leak_nA",
  699. color=COL_LEAK, alpha_cloud=0.18, lw_mean=2.2,
  700. y_transform=lambda y: -y, label_mean="Leak (mean)")
  701. plot_cloud_and_mean(ax_d3, labels_inj, "time_ms", "i_inj_nA",
  702. color=COL_INJ, alpha_cloud=0.18, lw_mean=2.2,
  703. y_transform=lambda y: -y, label_mean="Injected (mean)")
  705. ax_d3.set_ylabel("Current (nA)")
  706. ax_d3.set_xlabel("Time (ms)")
  707. ax_d3.legend(frameon=False, fontsize=9, loc="upper right")
  709. # -----------------------------------------------------------------------------
  710. # Panel b: Scaling curves by E_K
  711. # -----------------------------------------------------------------------------
  712. plot_scaling_curves_by_EK(
  713. ax_b1,
  714. df=df_no_inj,
  715. rin_fixed=RIN_DEFAULT,
  716. ek_values=ek_values,
  717. scaling_order=scaling_order,
  718. ek_color_map=EK_COLOR_MAP,
  719. xlabel="HCN2-4 scaling",
  720. ylabel=r"Max $V$ (mV)",
  721. point_size=34, # larger
  722. point_alpha=0.60 # less transparent
  723. )
  725. # -----------------------------------------------------------------------------
  726. # Panel c: R_in parameter sweep (same E_K, scaling as panels a/d → match panel b color)
  727. # -----------------------------------------------------------------------------
  728. col_c = EK_COLOR_MAP[float(E_K_ref_for_rin)]
  729. plot_dist(
  730. ax_c1, rin_order, df_rin, "R_in_Mohm",
  731. xlabel=r"$R_{\mathrm{in}}$ (M$\Omega$)",
  732. ylabel=None,
  733. semilogx=True,
  734. point_size=34,
  735. point_alpha=0.65,
  736. color=col_c
  737. )
  738. ax_c1.tick_params(labelleft=False)
  740. ax_b1.set_title(rf"$R_{{\mathrm{{in}}}} = {RIN_DEFAULT:.0f}\,\mathrm{{M\Omega}}$", fontsize=10, pad=4)
  741. ax_b1.legend(frameon=False, fontsize=9, loc="best")
  742. # Remove x-axis line and ticks to emphasize categorical nature of scaling values
  743. # ax_b1.spines["bottom"].set_visible(False)
  744. # ax_b1.tick_params(bottom=False, labelbottom=True) # Keep labels, remove ticks
  745. ax_c1.set_title("mean ± SD (replicates)", fontsize=10, pad=4)
  747. # -----------------------------------------------------------------------------
  748. # Styling & panel labels
  749. # -----------------------------------------------------------------------------
  750. for ax in (ax_a1, ax_a2, ax_a3, ax_d1, ax_d2, ax_d3, ax_c1, ax_b1):
  751. ax.spines["top"].set_visible(False)
  752. ax.spines["right"].set_visible(False)
  753. ax.grid(False)
  755. # Panel letters aligned to the top of each block
  756. bbox_a = ax_a1.get_position()
  757. bbox_d = ax_d1.get_position()
  758. bbox_b = ax_b1.get_position()
  759. bbox_c = ax_c1.get_position()
  761. fig.text(bbox_a.x0 - 0.06, bbox_a.y1, "a", fontsize=13, fontweight="bold", va="top", ha="left")
  762. fig.text(bbox_d.x0 - 0.06, bbox_d.y1, "d", fontsize=13, fontweight="bold", va="top", ha="left")
  763. fig.text(bbox_b.x0 - 0.06, bbox_b.y1, "b", fontsize=13, fontweight="bold", va="top", ha="left")
  764. fig.text(bbox_c.x0 - 0.06, bbox_c.y1, "c", fontsize=13, fontweight="bold", va="top", ha="left")
  766. # -----------------------------------------------------------------------------
  767. # Export data to Excel
  768. # -----------------------------------------------------------------------------
  769. excel_file = "figure_data.xlsx"
  770. with pd.ExcelWriter(excel_file, engine='openpyxl') as writer:
  772. # Panel a: Non-injected traces
  773. # Get mean traces
  774. t_mean_a, v_mean_a = mean_trace(time_series_data, labels_noinj, "v_mV")
  775. _, p_mean_a = mean_trace(time_series_data, labels_noinj, "p_open")
  776. _, i_hcn1_mean_a = mean_trace(time_series_data, labels_noinj, "i_hcn1_nA")
  777. _, i_leak_mean_a = mean_trace(time_series_data, labels_noinj, "i_leak_nA")
  779. # Create DataFrame with mean traces
  780. df_panel_a = pd.DataFrame({
  781. "Time_ms": t_mean_a,
  782. "Voltage_mean_mV": v_mean_a,
  783. "P_open_mean": p_mean_a,
  784. "I_HCN1_mean_nA": -i_hcn1_mean_a, # Apply negative transform as in plot
  785. "I_leak_mean_nA": -i_leak_mean_a # Apply negative transform as in plot
  786. })
  788. # Add individual replicate traces (interpolated to mean time grid)
  789. for i, label in enumerate(labels_noinj):
  790. ts = time_series_data[label]
  791. t_orig = np.asarray(ts["time_ms"], dtype=float)
  792. df_panel_a[f"Voltage_rep{i+1}_mV"] = np.interp(t_mean_a, t_orig, np.asarray(ts["v_mV"], dtype=float))
  793. df_panel_a[f"P_open_rep{i+1}"] = np.interp(t_mean_a, t_orig, np.asarray(ts["p_open"], dtype=float))
  794. df_panel_a[f"I_HCN1_rep{i+1}_nA"] = -np.interp(t_mean_a, t_orig, np.asarray(ts["i_hcn1_nA"], dtype=float))
  795. df_panel_a[f"I_leak_rep{i+1}_nA"] = -np.interp(t_mean_a, t_orig, np.asarray(ts["i_leak_nA"], dtype=float))
  797. df_panel_a.to_excel(writer, sheet_name="Panel a", index=False)
  799. # Panel b: Scaling curves by E_K
  800. # Find max number of replicates across all E_K and scaling combinations
  801. df_panel_b_temp = df_no_inj[(df_no_inj["R_in_Mohm"] == float(RIN_DEFAULT))].copy()
  802. max_reps_b = int(df_panel_b_temp.groupby(["E_K_mV", "scaling"])["V_max_mV"].count().max()) if not df_panel_b_temp.empty else 0
  804. df_panel_b_list = []
  805. for ek in ek_values:
  806. d = df_no_inj[(df_no_inj["E_K_mV"] == float(ek)) &
  807. (df_no_inj["R_in_Mohm"] == float(RIN_DEFAULT))].copy()
  808. if d.empty:
  809. continue
  811. g_b = (d.groupby("scaling", as_index=False)
  812. .agg(V_max_mean_mV=("V_max_mV", "mean"),
  813. V_max_sd_mV=("V_max_mV", "std"),
  814. n_replicates=("V_max_mV", "count")))
  815. g_b = g_b.set_index("scaling").reindex(scaling_order).reset_index()
  816. g_b["E_K_mV"] = ek
  818. # Create replicate columns
  819. for rep_idx in range(1, max_reps_b + 1):
  820. col_name = f"V_max_rep{rep_idx}_mV"
  821. g_b[col_name] = np.nan
  823. # Add individual replicate values
  824. for sc_val in scaling_order:
  825. reps = d[d["scaling"] == sc_val]["V_max_mV"].values
  826. for j, rep_val in enumerate(reps):
  827. col_name = f"V_max_rep{j+1}_mV"
  828. g_b.loc[g_b["scaling"] == sc_val, col_name] = rep_val
  830. df_panel_b_list.append(g_b)
  832. if df_panel_b_list:
  833. df_panel_b = pd.concat(df_panel_b_list, ignore_index=True)
  834. # Reorder columns: E_K, scaling, then statistics and replicates
  835. cols = ["E_K_mV", "scaling", "V_max_mean_mV", "V_max_sd_mV", "n_replicates"]
  836. rep_cols = [c for c in df_panel_b.columns if c.startswith("V_max_rep")]
  837. df_panel_b = df_panel_b[cols + rep_cols]
  838. df_panel_b.to_excel(writer, sheet_name="Panel b", index=False)
  840. # Panel c: R_in parameter sweep
  841. g_c = (df_rin.groupby("R_in_Mohm", as_index=False)
  842. .agg(V_max_mean_mV=("V_max_mV", "mean"),
  843. V_max_sd_mV=("V_max_mV", "std"),
  844. n_replicates=("V_max_mV", "count")))
  845. g_c = g_c.set_index("R_in_Mohm").reindex(rin_order).reset_index()
  847. # Find max number of replicates to create all columns
  848. max_reps = int(df_rin.groupby("R_in_Mohm")["V_max_mV"].count().max())
  850. # Add individual replicate values
  851. df_panel_c = g_c.copy()
  852. for rep_idx in range(1, max_reps + 1):
  853. col_name = f"V_max_rep{rep_idx}_mV"
  854. df_panel_c[col_name] = np.nan
  856. for rin_val in rin_order:
  857. reps = df_rin[df_rin["R_in_Mohm"] == rin_val]["V_max_mV"].values
  858. for j, rep_val in enumerate(reps):
  859. col_name = f"V_max_rep{j+1}_mV"
  860. df_panel_c.loc[df_panel_c["R_in_Mohm"] == rin_val, col_name] = rep_val
  862. df_panel_c.to_excel(writer, sheet_name="Panel c", index=False)
  864. # Panel d: Injected traces
  865. # Get mean traces
  866. t_mean_d, v_mean_d = mean_trace(time_series_data, labels_inj, "v_mV")
  867. _, p_mean_d = mean_trace(time_series_data, labels_inj, "p_open")
  868. _, i_hcn1_mean_d = mean_trace(time_series_data, labels_inj, "i_hcn1_nA")
  869. _, i_leak_mean_d = mean_trace(time_series_data, labels_inj, "i_leak_nA")
  870. _, i_inj_mean_d = mean_trace(time_series_data, labels_inj, "i_inj_nA")
  872. # Create DataFrame with mean traces
  873. df_panel_d = pd.DataFrame({
  874. "Time_ms": t_mean_d,
  875. "Voltage_mean_mV": v_mean_d,
  876. "P_open_mean": p_mean_d,
  877. "I_HCN1_mean_nA": -i_hcn1_mean_d, # Apply negative transform as in plot
  878. "I_leak_mean_nA": -i_leak_mean_d, # Apply negative transform as in plot
  879. "I_injected_mean_nA": -i_inj_mean_d # Apply negative transform as in plot
  880. })
  882. # Add individual replicate traces (interpolated to mean time grid)
  883. for i, label in enumerate(labels_inj):
  884. ts = time_series_data[label]
  885. t_orig = np.asarray(ts["time_ms"], dtype=float)
  886. df_panel_d[f"Voltage_rep{i+1}_mV"] = np.interp(t_mean_d, t_orig, np.asarray(ts["v_mV"], dtype=float))
  887. df_panel_d[f"P_open_rep{i+1}"] = np.interp(t_mean_d, t_orig, np.asarray(ts["p_open"], dtype=float))
  888. df_panel_d[f"I_HCN1_rep{i+1}_nA"] = -np.interp(t_mean_d, t_orig, np.asarray(ts["i_hcn1_nA"], dtype=float))
  889. df_panel_d[f"I_leak_rep{i+1}_nA"] = -np.interp(t_mean_d, t_orig, np.asarray(ts["i_leak_nA"], dtype=float))
  890. df_panel_d[f"I_injected_rep{i+1}_nA"] = -np.interp(t_mean_d, t_orig, np.asarray(ts["i_inj_nA"], dtype=float))
  892. df_panel_d.to_excel(writer, sheet_name="Panel d", index=False)
  894. print(f"Data exported to {excel_file}")
  896. # -----------------------------------------------------------------------------
  897. # Save
  898. # -----------------------------------------------------------------------------
  899. fig.savefig("figure_combined.svg", format="svg", bbox_inches="tight")
  900. fig.savefig("figure_combined.pdf", format="pdf", bbox_inches="tight")
  901. plt.show()

HCN_rest-singlePo.py, no license · at the source

Overview

Authors: Uta Enke1, Andrea Schweinitz1, Debanjan Tewari1, Christian Sattler1, Ralf Schmauder1, Christoph Schmidt-Hieber2, Klaus Benndorf1
ORCID iDs: Uta Enke, Christian Sattler, Ralf Schmauder, Christoph Schmidt-Hieber, Klaus Benndorf
  1. Institut für Physiologie II, Universitätsklinikum Jena, Friedrich-Schiller-Universität Jena, Jena, Germany
  2. Institut für Physiologie I, Universitätsklinikum Jena, Friedrich-Schiller-Universität Jena, Jena, Germany
Institutions: Jena University Hospital (Germany); Friedrich Schiller University Jena (Germany)
Journal: Nature communications, volume 17, issue 1, article 3745
Dates: received 21 August 2025; accepted 10 April 2026; published online 23 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41467-026-72257-3 · PMID 42026039 · PMCID PMC13106851 · OpenAlex W7155356039
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: intracellular / patch clamp (modality), human (organism), cellular / molecular (subfield)
Keywords: Neurophysiology, Single-molecule biophysics, Single-channel recording, Patch clamp, Ion channels
MeSH: Biological Clocks*, Cyclic Nucleotide-Gated Cation Channels*, Hyperpolarization-Activated Cyclic Nucleotide-Gated Channels*, Neurons*, Potassium Channels*, Action Potentials, Animals, Humans, Ion Channel Gating, Patch-Clamp Techniques (* major topic)
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Deutsche Forschungsgemeinschaft (German Research Foundation) (RU 2518 P02)
Citations: not cited yet (Europe PMC); 49 references in the paper

Abstract

Rhythmic activity of specialized pacemaker neurons in the brain is necessary to control alertness and circadian timing. Four HCN channels have been identified to generate the pacemaker current Ih or Iq, differing in activation speed, voltage dependence, single-channel conductance, and cAMP sensitivity. Here we show the time-resolved operation of single HCN1, HCN2 and HCN4 channels during the pacemaker depolarization using a dynamic neuronal action potential clamp at femtosiemens resolution. All channels produce a relevant open probability during pacemaker depolarization. However, only mHCN1 channels are significantly activated and deactivated in action potential cycles whereas the gating in mHCN2 and mHCN4 channels is at best barely resolvable and too slow. Simulations suggest that the role of HCN1 channels is to trigger the initial neuronal pacemaker depolarization before other depolarizing conductances take over this role. In conclusion, mHCN1 channels are the primary HCN pacemaker channels that operate as trigger channels for pacemaking.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

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

OSF 7g5ch

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Languages: Python (1)
Size: 4 files, 1 script
Software Heritage: not checked
Found in: “Code availability”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Brian 2 (1 file), Matplotlib (1 file), NumPy (1 file), pandas (1 file), SciPy (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
3 files
At the source: osf.io/7g5ch

Code availability

The code of the simulations (HCN_rest-singlePo.py) is available at OSF [osf.io/7g5ch]. No manuscript-specific software was used for the analysis of measured data.

Reproduced under the paper's license (CC BY), from the paper cited above.

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:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 1 script, each with its path and the digest of its content;
  • 3 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

All data included in the manuscript are deposited in the Open Science Framework [osf.io/7g5ch].

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 5 keywords, 10 MeSH terms, 1 funder, 48 references.

Cite

This paper

Enke, U., Schweinitz, A., Tewari, D., Sattler, C., Schmauder, R., Schmidt-Hieber, C., & Benndorf, K. (2026). HCN1 is a primary HCN Pacemaker Channel in Neurons. Nature communications, 17(1), 3745. https://doi.org/10.1038/s41467-026-72257-3

BibTeX

@article{enke2026hcn1,
author = {Enke, Uta and Schweinitz, Andrea and Tewari, Debanjan and Sattler, Christian and Schmauder, Ralf and Schmidt-Hieber, Christoph and Benndorf, Klaus},
title = {{HCN1 is a primary HCN Pacemaker Channel in Neurons}},
journal = {Nature communications},
year = {2026},
month = apr,
volume = {17},
number = {1},
pages = {3745},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-72257-3},
url = {https://doi.org/10.1038/s41467-026-72257-3},
pmid = {42026039},
pmcid = {PMC13106851}
}

RIS

TY - JOUR
AU - Enke, Uta
AU - Schweinitz, Andrea
AU - Tewari, Debanjan
AU - Sattler, Christian
AU - Schmauder, Ralf
AU - Schmidt-Hieber, Christoph
AU - Benndorf, Klaus
TI - HCN1 is a primary HCN Pacemaker Channel in Neurons
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/04/23
VL - 17
IS - 1
SP - 3745
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-72257-3
UR - https://doi.org/10.1038/s41467-026-72257-3
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-72257-3",
"type": "article-journal",
"title": "HCN1 is a primary HCN Pacemaker Channel in Neurons",
"container-title": "Nature communications",
"author": [
{
"family": "Enke",
"given": "Uta"
},
{
"family": "Schweinitz",
"given": "Andrea"
},
{
"family": "Tewari",
"given": "Debanjan"
},
{
"family": "Sattler",
"given": "Christian"
},
{
"family": "Schmauder",
"given": "Ralf"
},
{
"family": "Schmidt-Hieber",
"given": "Christoph"
},
{
"family": "Benndorf",
"given": "Klaus"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "3745",
"DOI": "10.1038/s41467-026-72257-3",
"PMID": "42026039",
"PMCID": "PMC13106851",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-72257-3",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
23
]
]
}
}

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

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1113/ep093939 [code]
Temperature dependence of HCN channel kinetics: A systematic comparison across mammalian isoforms and species.
Journal: Experimental physiology
In common: intracellular / patch clamp, cellular / molecular, 4 references
[2] doi:10.1016/j.isci.2026.116930
HCAR1-mediated lactate signaling modulates motor behavior and regulates spontaneous firing in Purkinje cells.
Journal: iScience
In common: cellular / molecular, 4 references
[3] doi:10.1523/jneurosci.1506-25.2026 [code]
Controlling Spatio-Temporal Sequences of Neural Activity by Local Synaptic Changes.
Journal: The Journal of neuroscience : the official journal of the Society for Neuroscience
In common: Brian 2, pandas, SciPy, 2 other tools, cellular / molecular, 1 reference
[4] doi:10.1002/hipo.70089 [code]
The Role of Plasticity in Replay: Stability Through Anti-Hebbian Rules.
Journal: Hippocampus
In common: Brian 2, pandas, SciPy, 2 other tools, cellular / molecular, 1 reference
[5] doi:10.1371/journal.pcbi.1014752 [code]
Hierarchical feature binding in a spiking neural network model of the primate ventral visual pathway.
Journal: PLoS computational biology
In common: Brian 2, pandas, SciPy, 2 other tools, 1 reference
[6] doi:10.1126/sciadv.aee9425 [code]
Probabilistic inference of homonymous and heteronymous recurrent inhibition in human muscles from large-scale motor neuron recordings.
Journal: Science advances
In common: Brian 2, pandas, SciPy, 2 other tools, 1 reference
[7] doi:10.1038/s41467-026-75705-2 [code]
Redundant prefrontal hemispheres adapt storage strategy to working memory demands.
Journal: Nature communications
In common: Brian 2, pandas, SciPy, 2 other tools, 1 reference
[8] doi:10.1523/jneurosci.0912-25.2026 [code]
Hierarchical Afferent Connectivity Drives Population-Wide Bursting Dynamics in a Computational Model of Human-Derived Excitatory Neuronal Networks.
Journal: The Journal of neuroscience : the official journal of the Society for Neuroscience
In common: Brian 2, pandas, SciPy, 2 other tools, 1 reference
[9] doi:10.1038/s41467-026-75704-3 [code]
A minimal model of working memory in neural systems and neuromorphic circuits.
Journal: Nature communications
In common: Brian 2, SciPy, Matplotlib, 1 other tool, 1 reference
[10] doi:10.1016/j.celrep.2026.117793 [code]
Clustered inputs engage dendritic nonlinearities and calcium signaling to support efficient place-field formation in CA1 pyramidal neurons.
Journal: Cell reports
In common: Brian 2, pandas, SciPy, 2 other tools, cellular / molecular

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.