OSCR

Neuronal excitability and parameter variability in the Hodgkin-Huxley model.

Code ↔ Paper

14 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 14 matches
  1. [1] § Methods › Spatial discretization and boundary conditions ↔ Figures 4-7/core_simulation.py, lines 28–131 · score 0.86 · steady state gating, Newton step, initial guesses, damped, iteration, smaller
  2. [2] § Methods › Spatial discretization and boundary conditions ↔ Figure 8/sim_core.py, lines 20–122 · score 0.83 · steady state gating, Newton step, initial guesses, damped, iteration, derivative
  3. [3] § Results ↔ Figures 4-7/core_analysis.py, lines 536–602 · score 0.75 · propagation speed, action potential propagated, distal site, Successful propagation, injection site, distance
  4. [4] § Results ↔ Figure 8/analysis.py, lines 17–166 · score 0.67 · upward threshold crossing, stimulus period, stimulus onset, APs, traces
  5. [5] § Results ↔ Figures 4-7/core_analysis.py, lines 292–419 · score 0.64 · upward threshold crossing, firing frequencies, duration, APs, histograms, onset
  6. [6] § Methods › Numerical integration and implementation ↔ Figures 4-7/core_simulation.py, lines 28–131 · score 0.63 · steady state, Newton, damped, iteration, enforcing, derivative
  7. [7] § Methods › Numerical integration and implementation ↔ Figures 4-7/core_simulation.py, lines 233–318 · score 0.61 · reversal potentials, gating variables, clipping, solver, interval, membrane
  8. [8] § Methods › Bootstrap parameter estimation ↔ Figure 1/analyze_ah.ipynb, lines 27–78 · score 0.61 · curve fit, bootstrap resampling, seed, Parameter uncertainties, guesses, replacement
  9. [9] § Methods › Bootstrap parameter estimation ↔ Figure 1/analyze_am.ipynb, lines 27–84 · score 0.61 · curve fit, bootstrap resampling, seed, Parameter uncertainties, guesses, replacement
  10. [10] § Methods › Numerical integration and implementation ↔ Figure 8/sim_core.py, lines 20–122 · score 0.61 · steady state, Newton, damped, iteration, enforcing, derivative
  11. [11] § Methods › Spatial discretization and boundary conditions ↔ Figures 4-7/core_simulation.py, lines 188–231 · score 0.58 · Sealed end, boundary, finite, dynamics, derivative, coupled
  12. [12] § Methods › Numerical integration and implementation ↔ Figure 8/sim_core.py, lines 125–239 · score 0.56 · Runge Kutta, Diffrax, adaptive, numerical, JAX, simulations
  13. [13] § Methods › Biophysical model of an excitable cable ↔ Figures 4-7/core_simulation.py, lines 233–318 · score 0.55 · Iext, axial, externally, ionic, radius, Cm
  14. [14] § Methods › Numerical integration and implementation ↔ Figure 8/sim_core.py, lines 125–239 · score 0.52 · Model parameters, safeguards, clipping, solver, numerical, variables

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 · 582 lines · 23 KB · no license · 5 matches

  1. """
  2. Core simulation module for the JAX-optimized spatial Hodgkin-Huxley cable model.
  4. Contains the JAX-compiled simulation functions for running the cable model.
  5. """
  7. # Configure XLA environment BEFORE importing JAX
  8. from config import configure_xla_environment
  9. configure_xla_environment()
  11. import jax
  12. import jax.numpy as jnp
  13. from jax import jit, vmap, lax, pmap
  14. import diffrax
  15. from functools import partial
  17. # Import configuration
  18. from config import (
  19. N_SEGMENTS, DX, LAMBDA_CABLE, INJECTION_SEGMENT, RECORDING_SEGMENT,
  20. T_START, T_END, RAMP_T, A_RADIUS
  21. )
  23. # Set JAX to use CPU and enable 64-bit precision
  24. jax.config.update('jax_platform_name', 'cpu')
  25. jax.config.update("jax_enable_x64", True)
  28. def find_equilibrium_jax(params):
  29. """
  30. Solve for an equilibrium (resting) potential V_rest by forcing net ionic current = 0
  31. with gating variables at their steady-state values.
  32. Uses Newton's method implemented in JAX for JIT compilation.
  34. For spatial cable model, returns same initial conditions for all N segments.
  35. """
  36. # Unpack parameters
  37. (gbar_Na, E_Na, A_am, B_am, C_am, A_bm, D_bm, A_ah, D_ah,
  38. E_bh, F_bh, G_bh,
  39. gbar_K, E_K, A_alpha, V_alpha, k_alpha, A_beta, tau_beta,
  40. C_m, E_l, G_l) = params
  42. # Define gating functions using JAX operations
  43. def alpha_m(V):
  44. denom = 1.0 - jnp.exp(-(V + B_am) / C_am)
  45. # Use where to handle division safely
  46. return jnp.where(jnp.abs(denom) > 1e-9, A_am * (V + B_am) / denom, A_am * C_am)
  48. def beta_m(V):
  49. V_safe = jnp.clip(V, -150.0, 100.0)
  50. return A_bm * jnp.exp(-V_safe / D_bm)
  52. def alpha_h(V):
  53. V_safe = jnp.clip(V, -150.0, 100.0)
  54. return A_ah * jnp.exp(-V_safe / D_ah)
  56. def beta_h(V):
  57. return E_bh / (1.0 + jnp.exp(-(V + F_bh) / G_bh))
  59. def alpha_n(V):
  60. V_safe = jnp.clip(V, -150.0, 100.0)
  61. denom = 1.0 - jnp.exp(-(V_alpha + V_safe) / k_alpha)
  62. return jnp.where(jnp.abs(denom) > 1e-9, A_alpha * (V_alpha + V_safe) / denom, A_alpha * k_alpha)
  64. def beta_n(V):
  65. V_safe = jnp.clip(V, -150.0, 100.0)
  66. return A_beta * jnp.exp(-V_safe / tau_beta)
  68. # Steady-state gating
  69. def m_inf(V):
  70. return alpha_m(V) / (alpha_m(V) + beta_m(V))
  72. def h_inf(V):
  73. return alpha_h(V) / (alpha_h(V) + beta_h(V))
  75. def n_inf(V):
  76. return alpha_n(V) / (alpha_n(V) + beta_n(V))
  78. # Net current at steady-state gating
  79. def net_current(V):
  80. mm = m_inf(V)
  81. hh = h_inf(V)
  82. nn = n_inf(V)
  83. I_Na = gbar_Na * mm**3 * hh * (V - E_Na)
  84. I_K = gbar_K * nn**4 * (V - E_K)
  85. I_L = G_l * (V - E_l)
  86. return I_Na + I_K + I_L
  88. # Robust Newton's method with bounds and damping
  89. V = -65.0 # Initial guess
  90. V_min = -100.0 # Lower bound
  91. V_max = 50.0 # Upper bound
  93. for i in range(100): # Increased iterations for better convergence
  94. f = net_current(V)
  96. # Numerical derivative with small step
  97. h_step = 1e-7
  98. df = (net_current(V + h_step) - f) / h_step
  100. # Prevent division by zero or very small derivatives
  101. df = jnp.where(jnp.abs(df) < 1e-10, 1e-10, df)
  103. # Newton step with adaptive damping
  104. dV = -f / df
  105. # Limit step size to prevent overshooting
  106. dV = jnp.clip(dV, -5.0, 5.0) # Smaller steps for stability
  108. V_new = V + 0.3 * dV # Stronger damping (0.3 instead of 0.5)
  110. # Enforce bounds
  111. V_new = jnp.clip(V_new, V_min, V_max)
  113. # Check for convergence
  114. change = jnp.abs(V_new - V)
  115. V = V_new
  117. # Early exit if converged (but keep iterating for JIT)
  118. V = jnp.where(change < 1e-8, V, V) # Tighter tolerance
  120. V_rest = V
  121. m0 = m_inf(V_rest)
  122. h0 = h_inf(V_rest)
  123. n0 = n_inf(V_rest)
  125. # Return initial conditions for all segments (uniform across axon)
  126. V0_all = jnp.tile(V_rest, N_SEGMENTS)
  127. m0_all = jnp.tile(m0, N_SEGMENTS)
  128. h0_all = jnp.tile(h0, N_SEGMENTS)
  129. n0_all = jnp.tile(n0, N_SEGMENTS)
  131. return V0_all, m0_all, h0_all, n0_all
  134. def simulate_HH_AP_jax(params, init_conditions, tmax=100.0, dt=0.005, I_amplitude=10.0):
  135. """
  136. Simulate the spatially-extended Hodgkin-Huxley cable equation using diffrax PDE solver.
  137. JAX-compatible and JIT-compiled for performance.
  139. Cable geometry: 10 cm length, 0.5 mm diameter, 80 segments
  140. Stimulus applied at injection site, voltage recorded at injection and distal sites
  141. Uses adaptive Tsit5 solver for accuracy with stiff cable equations
  143. **Memory-efficient version: Only returns voltage at 2 key locations**
  145. Parameters:
  146. -----------
  147. params : tuple
  148. Hodgkin-Huxley parameters
  149. init_conditions : tuple
  150. Initial conditions (V0_all, m0_all, h0_all, n0_all)
  151. tmax : float
  152. Simulation duration (ms)
  153. dt : float
  154. Initial timestep hint (ms)
  155. I_amplitude : float
  156. Total current injection amplitude in μA
  158. Returns:
  159. --------
  160. times : array
  161. Time vector
  162. V_injection : array
  163. Voltage at injection site
  164. V_distal : array
  165. Voltage at distal recording site
  166. """
  167. (gbar_Na, E_Na, A_am, B_am, C_am, A_bm, D_bm, A_ah, D_ah,
  168. E_bh, F_bh, G_bh,
  169. gbar_K, E_K, A_alpha, V_alpha, k_alpha, A_beta, tau_beta,
  170. C_m, E_l, G_l) = params
  172. # Gating functions for Na (m,h) using JAX operations
  173. def alpha_m(V):
  174. denom = 1.0 - jnp.exp(-(V + B_am) / C_am)
  175. return jnp.where(jnp.abs(denom) > 1e-9, A_am*(V + B_am)/denom, A_am*C_am)
  177. def beta_m(V):
  178. V_safe = jnp.clip(V, -150.0, 100.0)
  179. return A_bm * jnp.exp(-V_safe / D_bm)
  181. def alpha_h(V):
  182. V_safe = jnp.clip(V, -150.0, 100.0)
  183. return A_ah * jnp.exp(-V_safe / D_ah)
  185. def beta_h(V):
  186. return E_bh / (1.0 + jnp.exp(-(V + F_bh)/G_bh))
  188. # K gating
  189. def alpha_n(V):
  190. V_safe = jnp.clip(V, -150.0, 100.0) # Prevent exp overflow
  191. denom = 1.0 - jnp.exp(-(V_alpha + V_safe) / k_alpha)
  192. return jnp.where(jnp.abs(denom) > 1e-9, A_alpha*(V_alpha + V_safe)/denom, A_alpha*k_alpha)
  194. def beta_n(V):
  195. V_safe = jnp.clip(V, -150.0, 100.0) # Prevent exp overflow
  196. return A_beta * jnp.exp(-V_safe / tau_beta)
  198. # Cable equation PDE system for diffrax with spatial coupling
  199. def cable_hh_dynamics(t, y, args):
  200. I_amp = args
  201. N = N_SEGMENTS
  203. # Extract state variables for all segments
  204. # State vector: [V[0:N], m[N:2N], h[2N:3N], n[3N:4N]]
  205. V = y[0:N]
  206. m = y[N:2*N]
  207. h = y[2*N:3*N]
  208. n = y[3*N:4*N]
  210. # Clip state variables to prevent numerical instability
  211. # Clip V below E_Na to give the constraint enforcement room to work
  212. V = jnp.clip(V, -150.0, E_Na)
  213. m = jnp.clip(m, 0.0, 1.0)
  214. h = jnp.clip(h, 0.0, 1.0)
  215. n = jnp.clip(n, 0.0, 1.0)
  217. # Calculate spatial derivative (second derivative for diffusion)
  218. # d²V/dx² using finite differences
  219. d2Vdx2 = jnp.zeros(N)
  221. # Interior points: central difference
  222. d2Vdx2 = d2Vdx2.at[1:N-1].set(
  223. (V[2:N] - 2*V[1:N-1] + V[0:N-2]) / (DX**2)
  224. )
  226. # Boundary conditions: sealed ends (no-flux boundary, dV/dx = 0)
  227. # For no-flux BC, use ghost point approach: V[-1]=V[1] and V[N]=V[N-2]
  228. # Left boundary: d²V/dx²[0] = (V[1] - 2*V[0] + V[1]) / DX² = 2*(V[1] - V[0]) / DX²
  229. # Right boundary: d²V/dx²[N-1] = (V[N-2] - 2*V[N-1] + V[N-2]) / DX² = 2*(V[N-2] - V[N-1]) / DX²
  230. d2Vdx2 = d2Vdx2.at[0].set(2.0 * (V[1] - V[0]) / (DX**2))
  231. d2Vdx2 = d2Vdx2.at[N-1].set(2.0 * (V[N-2] - V[N-1]) / (DX**2))
  233. # Calculate ionic currents for all segments
  234. INa = gbar_Na * m**3 * h * (V - E_Na)
  235. IK = gbar_K * n**4 * (V - E_K)
  236. IL = G_l * (V - E_l)
  238. # External current injection with smooth ramping (1ms ramp up/down)
  239. # Convert total current (μA) to current density (μA/cm²) based on segment size
  240. # This ensures same total current regardless of spatial discretization
  241. Iext = jnp.zeros(N)
  242. ramp_up = jnp.clip((t - T_START) / RAMP_T, 0.0, 1.0)
  243. ramp_down = jnp.clip((T_END - t) / RAMP_T, 0.0, 1.0)
  244. # I_amp is total current (μA), convert to density by dividing by segment membrane area
  245. segment_membrane_area = 2.0 * jnp.pi * A_RADIUS * DX # cm²
  246. current_density = I_amp / segment_membrane_area # μA/cm²
  247. current_value = current_density * ramp_up * ramp_down
  248. Iext = Iext.at[INJECTION_SEGMENT].set(current_value)
  250. # Cable equation: C_m * dV/dt = I_axial_density - I_ion + I_ext
  251. # I_axial = LAMBDA_CABLE * d²V/dx² where LAMBDA_CABLE = πa²/R_a in S*cm
  252. # Units: (S*cm) * (mV/cm²) = mV*S/cm = mA/cm (current per unit length)
  253. # Convert to current density (μA/cm²) by dividing by membrane area per unit length (2πa):
  254. # Current density = (mA/cm) / (2πa cm) * 1000 = μA/cm²
  255. membrane_circumference = 2.0 * jnp.pi * A_RADIUS # cm
  256. spatial_current_density = (LAMBDA_CABLE * d2Vdx2 / membrane_circumference) * 1000.0 # μA/cm²
  257. dV = (spatial_current_density - INa - IK - IL + Iext) / C_m
  259. # Enforce reversal potential constraint: when V >= E_Na, dV/dt must be ≤ 0
  260. # This prevents V from exceeding E_Na regardless of injected current
  261. dV = jnp.where(V >= E_Na, jnp.minimum(dV, 0.0), dV)
  263. # Gating variable dynamics (same for all segments)
  264. dm = alpha_m(V) * (1.0 - m) - beta_m(V) * m
  265. dh = alpha_h(V) * (1.0 - h) - beta_h(V) * h
  266. dn = alpha_n(V) * (1.0 - n) - beta_n(V) * n
  268. # Return concatenated derivatives
  269. return jnp.concatenate([dV, dm, dh, dn])
  271. # Initial conditions - tuple of (V_all, m_all, h_all, n_all) each of shape (N_SEGMENTS,)
  272. V0_all, m0_all, h0_all, n0_all = init_conditions
  274. # Clamp to valid ranges
  275. V0_all = jnp.clip(V0_all, -150.0, 100.0)
  276. m0_all = jnp.clip(m0_all, 0.0, 1.0)
  277. h0_all = jnp.clip(h0_all, 0.0, 1.0)
  278. n0_all = jnp.clip(n0_all, 0.0, 1.0)
  280. # Concatenate all initial conditions
  281. y0 = jnp.concatenate([V0_all, m0_all, h0_all, n0_all])
  283. # Setup ODE solver - use Tsit5 (5th order) with adaptive stepping for better efficiency
  284. term = diffrax.ODETerm(cable_hh_dynamics)
  285. solver = diffrax.Tsit5()
  286. stepsize_controller = diffrax.PIDController(rtol=1e-5, atol=1e-6)
  288. # Save at 0.1 ms intervals instead of every timestep
  289. # This reduces memory by ~70× while maintaining sufficient resolution for analysis
  290. save_interval = 0.1 # ms
  291. ts = jnp.arange(0.0, tmax + save_interval/2, save_interval)
  292. saveat = diffrax.SaveAt(ts=ts)
  294. # Solve cable PDE with adaptive-step Tsit5 solver
  295. solution = diffrax.diffeqsolve(
  296. term,
  297. solver,
  298. t0=0.0,
  299. t1=tmax,
  300. dt0=dt, # Initial step size hint
  301. y0=y0,
  302. args=I_amplitude,
  303. saveat=saveat,
  304. stepsize_controller=stepsize_controller,
  305. max_steps=1000000, # High limit for adaptive stepping (Tsit5 needs many steps)
  306. throw=False # Don't throw errors, return result with error code
  307. )
  309. # Memory-efficient: Extract only the two voltage traces we need
  310. # Instead of returning full solution (800 vars × 70k timesteps = 560M values)
  311. # Return only 2 voltage traces (2 × 70k timesteps = 140k values per simulation)
  312. V_solution = solution.ys[:, :N_SEGMENTS] # Shape: (n_times, N_SEGMENTS)
  314. # Enforce V ≤ E_Na constraint on the extracted segments
  315. V_injection = jnp.clip(V_solution[:, INJECTION_SEGMENT], -150.0, E_Na)
  316. V_distal = jnp.clip(V_solution[:, RECORDING_SEGMENT], -150.0, E_Na)
  318. return solution.ts, V_injection, V_distal
  321. @partial(jit, static_argnames=('tmax', 'dt'))
  322. def simulate_batch_jax(params_batch, tmax=100.0, dt=0.02, I_amplitude=10.0):
  323. """
  324. Simulate multiple parameter sets using pmap for multi-core parallelization.
  325. Uses spatially-extended cable model.
  326. Records voltage at injection site and distal site for each simulation.
  328. Strategy: Split batch across CPU cores using pmap, with lax.scan within each core.
  329. This balances parallelism with memory efficiency.
  331. Parameters:
  332. -----------
  333. params_batch : array
  334. Batch of parameter sets
  335. tmax : float
  336. Simulation duration (ms)
  337. dt : float
  338. Initial timestep hint (ms)
  339. I_amplitude : float
  340. Total current injection amplitude in μA
  342. Returns:
  343. --------
  344. times : array
  345. Time vector
  346. voltage_injection : array
  347. Voltage traces at injection site (n_sims, n_times)
  348. voltage_distal : array
  349. Voltage traces at distal site (n_sims, n_times)
  350. """
  351. # Vectorize equilibrium finding across batch
  352. batch_find_eq = vmap(find_equilibrium_jax)
  353. V_rest_batch, m0_batch, h0_batch, n0_batch = batch_find_eq(params_batch)
  355. # Stack initial conditions: shape (batch_size, 4, N_SEGMENTS)
  356. # We'll pass this through and unpack in simulate_HH_AP_jax
  357. init_conds_batch = jnp.stack([V_rest_batch, m0_batch, h0_batch, n0_batch], axis=1)
  359. # Get number of devices
  360. n_devices = len(jax.devices())
  361. batch_size = params_batch.shape[0]
  363. # Pad batch to be divisible by n_devices if needed
  364. remainder = batch_size % n_devices
  365. if remainder != 0:
  366. pad_size = n_devices - remainder
  367. params_batch = jnp.concatenate([
  368. params_batch,
  369. jnp.repeat(params_batch[-1:], pad_size, axis=0)
  370. ], axis=0)
  371. init_conds_batch = jnp.concatenate([
  372. init_conds_batch,
  373. jnp.repeat(init_conds_batch[-1:], pad_size, axis=0)
  374. ], axis=0)
  376. # Reshape to (n_devices, sims_per_device, ...)
  377. per_device = params_batch.shape[0] // n_devices
  378. params_split = params_batch.reshape(n_devices, per_device, -1)
  379. # For init_conds, preserve the (4, N_SEGMENTS) structure
  380. init_conds_split = init_conds_batch.reshape(n_devices, per_device, 4, N_SEGMENTS)
  382. # Define function to run on each device using lax.scan
  383. def process_device_batch(params_device, init_conds_device):
  384. def scan_fn(carry, inputs):
  385. params, init_conds = inputs
  386. # Unpack init_conds from stacked array to tuple
  387. init_tuple = (init_conds[0], init_conds[1], init_conds[2], init_conds[3])
  388. times, V_inj, V_dist = simulate_HH_AP_jax(params, init_tuple, tmax, dt, I_amplitude)
  389. return carry, (V_inj, V_dist)
  391. _, (V_inj_batch, V_dist_batch) = lax.scan(
  392. scan_fn,
  393. None,
  394. (params_device, init_conds_device)
  395. )
  396. return V_inj_batch, V_dist_batch
  398. # Use pmap to parallelize across devices
  399. pmap_fn = pmap(process_device_batch)
  400. voltage_injection_split, voltage_distal_split = pmap_fn(params_split, init_conds_split)
  402. # Reshape back to (total_batch_size, n_times)
  403. voltage_injection = voltage_injection_split.reshape(-1, voltage_injection_split.shape[-1])
  404. voltage_distal = voltage_distal_split.reshape(-1, voltage_distal_split.shape[-1])
  406. # Remove padding if we added any
  407. if remainder != 0:
  408. voltage_injection = voltage_injection[:batch_size]
  409. voltage_distal = voltage_distal[:batch_size]
  411. # Get times from first simulation (all identical)
  412. # Unpack first init_conds to tuple
  413. init_tuple_0 = (init_conds_batch[0, 0], init_conds_batch[0, 1], init_conds_batch[0, 2], init_conds_batch[0, 3])
  414. times, _, _ = simulate_HH_AP_jax(params_batch[0], init_tuple_0, tmax, dt, I_amplitude)
  416. return times, voltage_injection, voltage_distal
  419. def incremental_simulation_jax(param_values, time_test, dt, nT, chunk_size=5000, save_path='temp_results_jax', I_amplitude=10.0):
  420. """
  421. Run spatially-extended cable model simulations in chunks using JAX vectorization.
  422. Records voltage at injection site and distal site for statistical analysis.
  423. Much faster than multiprocessing approach due to JIT compilation and vectorization.
  425. Parameters:
  426. -----------
  427. param_values : array
  428. Parameter sets from sampling (NumPy array)
  429. time_test : array
  430. Time vector
  431. dt : float
  432. Time step
  433. nT : int
  434. Number of time points
  435. chunk_size : int
  436. Number of simulations to run before saving
  437. save_path : str
  438. Directory to save temporary results
  439. I_amplitude : float
  440. Total current injection amplitude in μA
  442. Returns:
  443. --------
  444. save_path : str
  445. Path to directory containing chunk files
  446. num_chunks : int
  447. Number of chunk files created
  448. """
  449. import os
  450. import time
  451. import numpy as np
  453. # Create temporary directory
  454. if not os.path.exists(save_path):
  455. os.makedirs(save_path)
  457. total_sims = len(param_values)
  459. # Warm-up compilation strategy:
  460. # - First chunk uses small size (50) to safely compile JIT code
  461. # - Subsequent chunks use larger size (500) to reuse compiled code for 10× speedup
  462. warmup_chunk_size = 50 # Safe for compilation
  463. production_chunk_size = chunk_size # Use requested size after warm-up
  465. # Calculate chunks: first chunk is warmup_size, rest are production_size
  466. warmup_sims = min(warmup_chunk_size, total_sims)
  467. remaining_sims = total_sims - warmup_sims
  468. num_production_chunks = int(np.ceil(remaining_sims / production_chunk_size)) if remaining_sims > 0 else 0
  469. num_chunks = 1 + num_production_chunks # 1 warmup + N production
  471. print(f"Running JAX simulations with warm-up compilation strategy:")
  472. print(f" Chunk 0 (warm-up): {warmup_chunk_size} simulations - triggers JIT compilation")
  473. print(f" Chunks 1-{num_chunks-1}: {production_chunk_size} simulations each - reuses compiled code")
  474. print(f" Total chunks: {num_chunks}, Total simulations: {total_sims}")
  475. print(f"\nExpected performance:")
  476. print(f" Chunk 0: ~30-50 seconds (compilation + execution)")
  477. print(f" Chunks 1+: ~5-10 seconds each (execution only, 10× larger batches!)\n")
  479. # Pre-compute static values to avoid recompilation
  480. tmax_static = float(time_test[-1])
  481. dt_static = float(dt)
  483. for chunk_idx in range(num_chunks):
  484. # Determine chunk boundaries based on warm-up vs production
  485. if chunk_idx == 0:
  486. # Warm-up chunk: small size for safe compilation
  487. chunk_start = 0
  488. chunk_end = warmup_sims
  489. current_chunk_size = warmup_chunk_size
  490. chunk_type = "WARM-UP (JIT compilation)"
  491. else:
  492. # Production chunk: large size using compiled code
  493. prod_idx = chunk_idx - 1
  494. chunk_start = warmup_sims + prod_idx * production_chunk_size
  495. chunk_end = min(warmup_sims + (prod_idx + 1) * production_chunk_size, total_sims)
  496. current_chunk_size = production_chunk_size
  497. chunk_type = "PRODUCTION (reusing compiled code)"
  499. actual_chunk_size = chunk_end - chunk_start
  501. print(f"\nProcessing JAX chunk {chunk_idx + 1}/{num_chunks} - {chunk_type}")
  502. print(f" Samples {chunk_start}-{chunk_end} ({actual_chunk_size} simulations)")
  503. chunk_start_time = time.time()
  505. # Get parameter subset for this chunk
  506. params_chunk = param_values[chunk_start:chunk_end]
  508. # Pad chunk to match expected size to avoid recompilation
  509. if actual_chunk_size < current_chunk_size:
  510. # Repeat last parameter set to fill chunk
  511. padding_size = current_chunk_size - actual_chunk_size
  512. padding = np.tile(params_chunk[-1:], (padding_size, 1))
  513. params_chunk = np.vstack([params_chunk, padding])
  514. print(f" Note: Padded from {actual_chunk_size} to {current_chunk_size} samples")
  516. # Convert to JAX array
  517. params_jax = jnp.array(params_chunk)
  519. try:
  520. # Run vectorized simulation
  521. if chunk_idx == 0:
  522. print(f" Compiling + running {current_chunk_size} simulations (I={I_amplitude} μA)...")
  523. else:
  524. print(f" Running {current_chunk_size} simulations in parallel (I={I_amplitude} μA)...")
  525. sim_start = time.time()
  526. _, voltage_injection, voltage_distal = simulate_batch_jax(params_jax, tmax_static, dt_static, I_amplitude)
  527. sim_time = time.time() - sim_start
  529. # Trim padding from results if needed
  530. if actual_chunk_size < current_chunk_size:
  531. voltage_injection = voltage_injection[:actual_chunk_size]
  532. voltage_distal = voltage_distal[:actual_chunk_size]
  534. # Convert back to NumPy for saving
  535. chunk_V_injection = np.array(voltage_injection, dtype=np.float32)
  536. chunk_V_distal = np.array(voltage_distal, dtype=np.float32)
  538. print(f" Simulations completed in {sim_time:.2f}s ({actual_chunk_size/sim_time:.1f} sims/sec)")
  540. # Save chunks to disk (separate files for each recording site)
  541. chunk_file_injection = os.path.join(save_path, f'chunk_{chunk_idx}_injection.npy')
  542. chunk_file_distal = os.path.join(save_path, f'chunk_{chunk_idx}_distal.npy')
  543. estimated_size_mb = (chunk_V_injection.nbytes + chunk_V_distal.nbytes) / (1024**2)
  544. print(f" Saving chunk {chunk_idx + 1}/{num_chunks} (~{estimated_size_mb:.1f} MB)...")
  546. np.save(chunk_file_injection, chunk_V_injection)
  547. np.save(chunk_file_distal, chunk_V_distal)
  549. # Verify save
  550. if (os.path.exists(chunk_file_injection) and os.path.getsize(chunk_file_injection) > 1000 and
  551. os.path.exists(chunk_file_distal) and os.path.getsize(chunk_file_distal) > 1000):
  552. chunk_total_time = time.time() - chunk_start_time
  553. print(f" ✓ Chunk {chunk_idx + 1}/{num_chunks} saved successfully (total: {chunk_total_time:.2f}s)")
  554. else:
  555. print(f" WARNING: Chunk {chunk_idx + 1} may not have saved correctly")
  557. except Exception as e:
  558. print(f" ERROR during chunk {chunk_idx + 1}: {e}")
  559. print(f" Attempting to save partial results...")
  561. # Try saving smaller sub-batches
  562. try:
  563. sub_batch_size = min(1000, actual_chunk_size // 4)
  564. for sub_idx in range(0, actual_chunk_size, sub_batch_size):
  565. sub_end = min(sub_idx + sub_batch_size, actual_chunk_size)
  566. sub_params = params_jax[sub_idx:sub_end]
  567. _, sub_injection, sub_distal = simulate_batch_jax(sub_params, tmax_static, dt_static, I_amplitude)
  568. sub_chunk_inj = np.array(sub_injection, dtype=np.float32)
  569. sub_chunk_dist = np.array(sub_distal, dtype=np.float32)
  570. sub_file_inj = os.path.join(save_path, f'chunk_{chunk_idx}_sub_{sub_idx//sub_batch_size}_injection.npy')
  571. sub_file_dist = os.path.join(save_path, f'chunk_{chunk_idx}_sub_{sub_idx//sub_batch_size}_distal.npy')
  572. np.save(sub_file_inj, sub_chunk_inj)
  573. np.save(sub_file_dist, sub_chunk_dist)
  574. print(f" Chunk {chunk_idx + 1} saved as sub-chunks")
  575. except Exception as e2:
  576. print(f" CRITICAL ERROR: Could not save chunk {chunk_idx + 1}: {e2}")
  577. raise
  579. # Free memory
  580. del chunk_V_injection, chunk_V_distal
  582. return save_path, num_chunks

core_simulation.py at commit cc19f51, no license · at the source

Overview

Authors: Alon Korngreen1,2
ORCID iDs: Alon Korngreen
  1. The Leslie and Susan Gonda Interdisciplinary Brain Research Center, Bar-Ilan University, Ramat Gan, Israel
  2. The Mina and Everard Goodman Faculty of Life Sciences, Bar-Ilan University, Ramat Gan, Israel
Institutions: Bar-Ilan University (Israel)
Journal: PLoS computational biology, volume 22, issue 6, article e1014458
Dates: received 12 February 2026; accepted 15 June 2026; published online 29 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pcbi.1014458 · PMID 42371975 · PMCID PMC13336477 · OpenAlex W7166522723
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: computational modeling (no new data) (modality), none (in silico) (organism), computational (subfield)
Methods: Connectivity, Statistics
MeSH: Action Potentials*, Models, Neurological*, Neurons*, Animals, Axons, Computational Biology, Computer Simulation, Decapodiformes, Monte Carlo Method, Sodium, Sodium Channels (* major topic)
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Citations: not cited yet (Europe PMC); 73 references in the paper

Abstract

Biophysically detailed neuron models are often built as a one-way pipeline in which voltage-clamp data are reduced to a single set of best-fit channel parameters, which are then combined into a deterministic spiking model. This practice discards experimentally observed scatter and fitting uncertainty, obscuring the mechanisms by which robustness and degeneracy arise in excitable systems. Here, I reintroduce fitted-parameter uncertainty into the Hodgkin-Huxley model and embed uncertainty and global sensitivity analysis into model construction. I digitized sodium and potassium rate-constant data from the original Hodgkin and Huxley figures and used bootstrap resampling to estimate uncertainty in the fitted voltage-dependent kinetic parameters. I then propagated these uncertainty estimates through a spatially extended squid axon cable model using large-scale Monte Carlo simulations, in which each sample defined a complete set of kinetic, conductance, passive, and structural parameters. At the channel level, first-order Sobol sensitivity indices revealed that all kinetic parameters contribute to output variance in a strongly time-dependent manner, with distinct parameters controlling transient and steady-state behavior for potassium and sodium conductances. At the level of neuronal excitability, the simulations produced a heterogeneous population of firing behaviors, including non-firing, phasic, regular, and spontaneous activity. Across stimulus amplitudes, the dominant firing mode was a single spike at stimulus onset, consistent with the physiological role of the squid giant axon in rapid escape behavior. The canonical 1952 Hodgkin-Huxley parameter set fell within the regularly firing minority subpopulation, rather than representing a unique or dominant solution. In the phasic subpopulation, action potential propagation and conduction velocity varied widely yet remained within experimental ranges. Finally, global sensitivity analysis during spiking showed uniformly small first-order Sobol indices but large total-order indices, indicating that excitability is governed primarily by strong interactions among all parameters rather than by any subset. Together, these results support reframing the Hodgkin-Huxley model as an experimentally constrained ensemble of behaviors rather than a single privileged parameter set, with physiologically relevant firing patterns emerging from structured regions of the parameter space.

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 14 matches between paragraphs and lines of code.

alon67/HodgkinHuxleyVariability

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: cc19f51beb9cfc884ff0768ccec6833cfb6f763d, 30 December 2025
Languages: Python (15), Jupyter (8)
Size: 39 files, 23 scripts
Software Heritage: not archived
Found in: the text, “Sensitivity analysis”
Holds: README, environment (Figure 8/requirements.txt, Figures 4-7/requirements.txt), 8 notebooks
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (21 files), Matplotlib (11 files), pandas (9 files), SciPy (9 files), JAX (3 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
24 files
At the source: github.com/alon67/HodgkinHuxleyVariability

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 23 scripts, each with its path and the digest of its content;
  • 14 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 simulations and analyses are publicly available on GitHub at https://github.com/alon67/HodgkinHuxleyVariability. No experimental data were generated for this study.

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

Versions

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

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 1 author, 11 MeSH terms, 67 references.

Cite

This paper

Korngreen, A. (2026). Neuronal excitability and parameter variability in the Hodgkin-Huxley model. PLoS computational biology, 22(6), e1014458. https://doi.org/10.1371/journal.pcbi.1014458

BibTeX

@article{korngreen2026neuronal,
author = {Korngreen, Alon},
title = {{Neuronal excitability and parameter variability in the Hodgkin-Huxley model}},
journal = {PLoS computational biology},
year = {2026},
month = jun,
volume = {22},
number = {6},
pages = {e1014458},
publisher = {PLOS},
issn = {1553-734X},
doi = {10.1371/journal.pcbi.1014458},
url = {https://doi.org/10.1371/journal.pcbi.1014458},
pmid = {42371975},
pmcid = {PMC13336477}
}

RIS

TY - JOUR
AU - Korngreen, Alon
TI - Neuronal excitability and parameter variability in the Hodgkin-Huxley model
T2 - PLoS computational biology
J2 - PLoS Comput Biol
PY - 2026
DA - 2026/06/29
VL - 22
IS - 6
SP - e1014458
SN - 1553-734X
PB - PLOS
DO - 10.1371/journal.pcbi.1014458
UR - https://doi.org/10.1371/journal.pcbi.1014458
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pcbi.1014458",
"type": "article-journal",
"title": "Neuronal excitability and parameter variability in the Hodgkin-Huxley model",
"container-title": "PLoS computational biology",
"author": [
{
"family": "Korngreen",
"given": "Alon"
}
],
"container-title-short": "PLoS Comput Biol",
"volume": "22",
"issue": "6",
"page": "e1014458",
"DOI": "10.1371/journal.pcbi.1014458",
"PMID": "42371975",
"PMCID": "PMC13336477",
"ISSN": "1553-734X",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pcbi.1014458",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
29
]
]
}
}

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.1371/journal.pcbi.1014337 [code]
Fast reconstruction of degenerate populations of conductance-based neuron models from spike times.
Journal: PLoS computational biology
In common: pandas, SciPy, Matplotlib, 1 other tool, computational, computational modeling (no new data), 7 references
[2] doi:10.1093/bioinformatics/btag328 [code]
eFEL: electrophysiology feature extraction library.
Journal: Bioinformatics (Oxford, England)
In common: pandas, SciPy, Matplotlib, 1 other tool, none (in silico), 4 references
[3] doi:10.1371/journal.pone.0354021 [code]
A genetic algorithm for self-supervised models of oscillatory neurodynamics.
Journal: PloS one
In common: JAX, SciPy, Matplotlib, 1 other tool, computational modeling (no new data), 3 references
[4] doi:10.1371/journal.pcbi.1014177 [code]
Activity-dependent neuromodulation and calcium homeostasis cooperate to produce robust and modulable neuronal function.
Journal: PLoS computational biology
In common: none (in silico), computational modeling (no new data), 5 references
[5] doi:10.1007/s10827-026-00936-7 [code]
When can neuronal activity-dependent homeostatic plasticity maintain circuit-level properties?
Journal: Journal of computational neuroscience
In common: Matplotlib, NumPy, none (in silico), computational modeling (no new data), 3 references
[6] doi:10.1073/pnas.2533168123 [code]
Dendritic morphology and synaptic nonlinearities enhance functional complexity in human cortical neurons.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: pandas, SciPy, Matplotlib, 1 other tool, 3 references
[7] doi:10.1371/journal.pcbi.1014617 [code]
An in silico framework for dissecting the mechanistic origins of in vivo recorded neuronal activity.
Journal: PLoS computational biology
In common: none (in silico), computational, computational modeling (no new data), 3 references
[8] doi:10.1371/journal.pcbi.1014283 [code]
Spatial richness of neural magnetic fields.
Journal: PLoS computational biology
In common: SciPy, Matplotlib, NumPy, computational, computational modeling (no new data), 2 references
[9] doi:10.3389/fncom.2026.1771884 [code]
NEURONpyxl: fast, flexible, Python-integrated simulation of biophysical neural networks with complex plastic synapses.
Journal: Frontiers in computational neuroscience
In common: pandas, SciPy, Matplotlib, 1 other tool, computational modeling (no new data), 2 references
[10] doi:10.1007/s11571-026-10522-3 [code]
Acetylcholine enhances deviance detection in Hodgkin-Huxley neuronal networks.
Journal: Cognitive neurodynamics
In common: none (in silico), computational, computational modeling (no new data), 3 references

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.