Limit-Cycle Proliferation Under Parametric Delayed Feedback in a Conductance-Based Neuron: Bifurcation Landscape, Orbit Catalog, and Capacity Analysis

Code ↔ Paper

36 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 36 matches

[1] § 3. Experimental Pipeline › 3.3. PS1: Write Protocol ↔ experiments/PS1_WriteProtocol.ipynb, lines 775–869 · score 0.90 · orbits achieve lock, Gate PS G1, orbit switch lock, median settling, ISI error, G1d
[2] § 3. Experimental Pipeline › 3.4. PS2: Read Protocol and Noise Robustness ↔ experiments/PS2_ReadProtocol.ipynb, lines 764–876 · score 0.88 · Gate PS G2, noisy accuracy, jitter accuracy, Clean accuracy, observation window, G2d
[3] § 2. Mathematical Framework › 2.2. Stability and the Orbit Landscape ↔ experiments/PS0b_FloquetValidation.ipynb, lines 468–548 · score 0.81 · trivial Floquet multipliers, QR iteration, phase shift, periodic orbit, history, delay
[4] § 4. Results › 4.1. PS0: The HH–DFC Orbit Catalog ↔ experiments/PS0_OrbitCatalogue.ipynb, lines 522–667 · score 0.76 · hierarchical clustering, qualitative categories, ISI separation, distinct orbit, linkage, map
[5] § 3. Experimental Pipeline › 3.2. PS0: Dense Orbit Catalog ↔ experiments/PS0_OrbitCatalogue.ipynb, lines 669–739 · score 0.76 · pairwise ISI separability, Gate PS G0, qualitative categories, distinct orbit, Catalog, PS0
[6] § 2. Mathematical Framework › 2.1. Hodgkin–Huxley Neuron with Delayed Feedback Control ↔ experiments/PS0b_FloquetValidation.ipynb, lines 715–843 · score 0.74 · Floquet multiplier validation, Floquet convergence, Period quality, ps0b floquetvalidation, trivial multiplier, empirical stability
[7] § 3. Experimental Pipeline › 3.5. PS3: Full Write–Read–Erase Demonstration ↔ experiments/PS3_FullDemo_n_100.ipynb, lines 786–876 · score 0.73 · erase verification, Gate PS G3, G3b, G3a, G3e, G3c
[8] § 3. Experimental Pipeline › 3.5. PS3: Full Write–Read–Erase Demonstration ↔ experiments/PS3_FullDemo_n_20.ipynb, lines 786–863 · score 0.73 · erase verification, Gate PS G3, G3b, G3a, G3e, G3c
[9] § 3. Experimental Pipeline › 3.2. PS0: Dense Orbit Catalog ↔ experiments/PS0_OrbitCatalogue.ipynb, lines 320–403 · score 0.72 · fingerprint vectors, quasi periodic, ISI CV, silent, cluster, chaotic
[10] § 3. Experimental Pipeline › 3.1. Gate Threshold Rationale ↔ experiments/PS0_OrbitCatalogue.ipynb, lines 669–739 · score 0.72 · ISI separability, PS G0, pairwise ISI separation, distinct orbit, catalog, PS0
[11] § 4. Results › 4.3. PS2: Read Protocol and Noise Robustness ↔ experiments/PS2_ReadProtocol.ipynb, lines 764–876 · score 0.70 · Gate PS G2, jitter accuracy, clean accuracy, G2d, G2a, G2b
[12] § 2. Mathematical Framework › 2.1. Hodgkin–Huxley Neuron with Delayed Feedback Control ↔ experiments/PS0_OrbitCatalogue.ipynb, lines 44–108 · score 0.70 · potassium activation, sodium activation, gating variables
[13] § 3. Experimental Pipeline › 3.7. PS5: Maximum Capacity from All 207 Orbit Types ↔ experiments/PS5_MaxCapacity.ipynb, lines 573–653 · score 0.70 · pairwise confusion matrix, greedy maximum subset, capacity curve, PS5, windows, selection
[14] § 4. Results › 4.2. PS1: Write Protocol Performance ↔ experiments/PS1_WriteProtocol.ipynb, lines 775–869 · score 0.67 · Median settling, orbits achieve, Lock rate, target ISI, G1c, G1a
[15] § 4. Results › 4.2. PS1: Write Protocol Performance ↔ experiments/PS1_WriteProtocol.ipynb, lines 871–987 · score 0.66 · ISI relative error, Switch lock rate, Median settling, switching matrix, PS1, Protocol
[16] § 4. Results › 4.1. PS0: The HH–DFC Orbit Catalog ↔ experiments/PS0b_FloquetValidation.ipynb, lines 468–548 · score 0.65 · QR iteration, period quality, Floquet multipliers, trivial, Validation, transient
[17] § 3. Experimental Pipeline › 3.4. PS2: Read Protocol and Noise Robustness ↔ experiments/PS2_ReadProtocol.ipynb, lines 1083–1183 · score 0.65 · Noise robustness, observation window, fingerprint templates, calibrated, jitter, transient
[18] § 4. Results › 4.6. PS5: Maximum Capacity ↔ experiments/PS5_MaxCapacity.ipynb, lines 278–351 · score 0.62 · pre simulation, ISI windows, confusion matrix, caching, PS5, Capacity
[19] § 4. Results › 4.6. PS5: Maximum Capacity ↔ experiments/PS5_MaxCapacity.ipynb, lines 573–653 · score 0.62 · Gate PS G5, greedy maximum subset, Capacity curve, PS5, accuracy, orbits
[20] § 3. Experimental Pipeline › 3.5. PS3: Full Write–Read–Erase Demonstration ↔ experiments/PS3_FullDemo_n_100.ipynb, lines 536–630 · score 0.58 · PS1 switching matrix, clean simulation, alphabet, sequence, transient, symbol
[21] § 3. Experimental Pipeline › 3.5. PS3: Full Write–Read–Erase Demonstration ↔ experiments/PS3_FullDemo_n_20.ipynb, lines 536–630 · score 0.58 · PS1 switching matrix, clean simulation, alphabet, sequence, transient, symbol
[22] § 4. Results › 4.5. PS4: Rate Coding Baseline Comparison ↔ experiments/PS4_RateBaseline.ipynb, lines 346–413 · score 0.57 · HH Rate capacity, HH Rate baseline, Capacity comparison, PS4, bar, HH DFC
[23] § 4. Results › 4.5. PS4: Rate Coding Baseline Comparison ↔ experiments/PS4_RateBaseline.ipynb, lines 438–534 · score 0.57 · Voltage traces, stable limit cycles, ISI series, PS4, Baseline, cm2
[24] § 4. Results › 4.3. PS2: Read Protocol and Noise Robustness ↔ experiments/PS2_ReadProtocol.ipynb, lines 878–924 · score 0.56 · max band, observation window, classification accuracy, calibration, PS2, min
[25] § 3. Experimental Pipeline › 3.2. PS0: Dense Orbit Catalog ↔ experiments/PS0_SecondaryIbias.ipynb, lines 626–752 · score 0.56 · 5.9–56.9 ms, 5.9 ms, dense, space, 200 ms, catalog
[26] § 4. Results › 4.5. PS4: Rate Coding Baseline Comparison ↔ experiments/PS5_MaxCapacity.ipynb, lines 498–571 · score 0.56 · capacity curve, Capacity comparison, HH Rate, bar, HH DFC, PS4
[27] § 2. Mathematical Framework › 2.1. Hodgkin–Huxley Neuron with Delayed Feedback Control ↔ experiments/PS0_OrbitCatalogue.ipynb, lines 110–232 · score 0.56 · warm started, delay buffer, RC3, voltage, transients, neuron
[28] § 4. Results › 4.5. PS4: Rate Coding Baseline Comparison ↔ experiments/PS4_RateBaseline.ipynb, lines 303–344 · score 0.55 · Gate PS G4, topological diversity, ISI range, DFC orbits, PS4, Baseline
[29] § 3. Experimental Pipeline › 3.5. PS3: Full Write–Read–Erase Demonstration ↔ experiments/PS5_MaxCapacity.ipynb, lines 278–351 · score 0.55 · evenly spacing, ISI window, capacity, accuracy, Phase, orbits
[30] § 3. Experimental Pipeline › 3.3. PS1: Write Protocol ↔ experiments/PS1_WriteProtocol.ipynb, lines 316–358 · score 0.55 · minimum pairwise ISI, maximizes, greedy, PS1, protocol, cataloged
[31] § 2. Mathematical Framework › 2.3. Orbit Library and Memory Formalism ↔ experiments/PS0_OrbitCatalogue.ipynb, lines 522–667 · score 0.54 · burst_p, qualitative categories, triplet, doublet, tonic, ISIs
[32] § 4. Results › 4.5. PS4: Rate Coding Baseline Comparison ↔ experiments/PS4_RateBaseline.ipynb, lines 346–413 · score 0.54 · Firing rate, HH rate, ISI CV, Hz, PS4, threshold
[33] § 2. Mathematical Framework › 2.3. Orbit Library and Memory Formalism ↔ experiments/PS0_SecondaryIbias.ipynb, lines 382–458 · score 0.53 · burst_p, qualitative categories, triplet, doublet, tonic, ISIs
[34] § Appendix A. HH Rate Functions and Parameters ↔ experiments/PS0_OrbitCatalogue.ipynb, lines 110–232 · score 0.53 · warm start, delay buffer, RC3, voltage, neuron, DFC
[35] § 2. Mathematical Framework › 2.4. ISI-Based Orbit Classification: POLD ↔ experiments/PS2_ReadProtocol.ipynb, lines 411–512 · score 0.51 · ISI sequence, floored, segment, correlation, zero, score
[36] § 4. Results › 4.2. PS1: Write Protocol Performance ↔ experiments/PS1_WriteProtocol.ipynb, lines 1009–1109 · score 0.50 · Gate PS G1, orbit switching matrix, settling, PS1, Protocol, delay

Paper

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

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Jupyter notebook · 1,038 lines · 35 KB · MIT · 8 matches

# %% [markdown]
# # PS0 — Dense Orbit Catalogue
# ## Option C: HH Delay-Directed Orbit Selection
# ### HHSMC Project — Corrected Pipeline — February 2026
#
# **Purpose:** Build comprehensive catalogue of all stable periodic orbits accessible via (K, τ) parameter switching at fixed I_bias.
#
# **Prerequisite:** ExpF0v3 → Gate G0 FAILED → Track B (Option C)
# - Strong chaos found (λ₁ up to 0.218 ms⁻¹) but only 2–4 ISI clusters
# - 75.1% periodic (5,711/7,600 pts) is the opportunity
#
# **Gate PS-G0:** ≥15 distinct orbit types, ≥3 qualitative categories
# %% [markdown]
# ## CELL 1 — Setup and Imports
# %%
import numpy as np
from numba import njit
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from matplotlib.colors import ListedColormap
import json, os, time, warnings
from datetime import datetime
from scipy.cluster.hierarchy import fcluster, linkage
from scipy.spatial.distance import pdist, squareform
from collections import Counter
# --- Google Drive Mount ---
try:
from google.colab import drive
drive.mount('/content/drive')
OUTPUT_DIR = '/content/drive/My Drive/HHSMC/full_study/PS0_orbit_catalogue'
ON_COLAB = True
except ImportError:
OUTPUT_DIR = './PS0_results'
ON_COLAB = False
os.makedirs(OUTPUT_DIR, exist_ok=True)
print(f"Output directory: {OUTPUT_DIR}")
print(f"Timestamp: {datetime.now().isoformat()}")
print(f"Phase PS0 — Dense Orbit Catalogue for Option C")
# %% [markdown]
# ## CELL 2 — HH Model (V-shifted convention, rest = 0)
# %%
# --- Fixed biophysical parameters ---
C_M = 1.0 # μF/cm²
G_NA = 120.0 # mS/cm²
G_K = 36.0 # mS/cm²
G_L = 0.3 # mS/cm²
E_NA = 115.0 # mV (V-shifted)
E_K = -12.0 # mV (V-shifted)
E_L = 10.6 # mV (V-shifted)
DT = 0.01 # ms integration step
@njit
def alpha_m(V):
"""Sodium activation rate with L'Hôpital limit."""
x = 25.0 - V
if abs(x) < 1e-7:
return 1.0
return 0.1 * x / (np.exp(x / 10.0) - 1.0)
@njit
def beta_m(V):
return 4.0 * np.exp(-V / 18.0)
@njit
def alpha_h(V):
return 0.07 * np.exp(-V / 20.0)
@njit
def beta_h(V):
return 1.0 / (np.exp((30.0 - V) / 10.0) + 1.0)
@njit
def alpha_n(V):
"""Potassium activation rate with L'Hôpital limit."""
x = 10.0 - V
if abs(x) < 1e-7:
return 0.1
return 0.01 * x / (np.exp(x / 10.0) - 1.0)
@njit
def beta_n(V):
return 0.125 * np.exp(-V / 80.0)
@njit
def hh_rhs(V, m, h, n, I_total):
"""Right-hand side of HH equations."""
I_Na = G_NA * m*m*m * h * (V - E_NA)
I_K = G_K * n*n*n*n * (V - E_K)
I_L = G_L * (V - E_L)
dV = (I_total - I_Na - I_K - I_L) / C_M
dm = alpha_m(V) * (1.0 - m) - beta_m(V) * m
dh = alpha_h(V) * (1.0 - h) - beta_h(V) * h
dn = alpha_n(V) * (1.0 - n) - beta_n(V) * n
return dV, dm, dh, dn
@njit
def hh_steady_state(V):
"""Steady-state gating variables at voltage V."""
am = alpha_m(V); bm = beta_m(V)
ah = alpha_h(V); bh = beta_h(V)
an = alpha_n(V); bn = beta_n(V)
return am/(am+bm), ah/(ah+bh), an/(an+bn)
# %% [markdown]
# ## CELL 3 — Simulation Engine (HH + DFC with RK4)
# %%
@njit
def simulate_hh_dfc_full(I_bias, K, tau_ms, T_total_ms, dt=0.01,
T_transient_ms=0.0):
"""
Simulate HH neuron with Pyragas DFC. Returns spike times, ISIs,
and the full voltage trace for pattern analysis.
Corrections applied:
- RC3: Warm-started delay buffer (free-run for τ ms)
- DFC recomputed at each RK4 substep
- L'Hôpital limits in all rate functions
"""
n_steps = int(T_total_ms / dt)
buf_size = max(int(tau_ms / dt), 1)
# Initial conditions at rest
V = 0.0
m, h, n = hh_steady_state(V)
# --- Warm-start delay buffer (RC3 fix) ---
V_buf = np.zeros(buf_size)
for ws in range(buf_size):
I_total = I_bias # No DFC during warmup
dV1, dm1, dh1, dn1 = hh_rhs(V, m, h, n, I_total)
V2 = V + 0.5*dt*dV1; m2 = m + 0.5*dt*dm1
h2 = h + 0.5*dt*dh1; n2 = n + 0.5*dt*dn1
dV2, dm2, dh2, dn2 = hh_rhs(V2, m2, h2, n2, I_total)
V3 = V + 0.5*dt*dV2; m3 = m + 0.5*dt*dm2
h3 = h + 0.5*dt*dh2; n3 = n + 0.5*dt*dn2
dV3, dm3, dh3, dn3 = hh_rhs(V3, m3, h3, n3, I_total)
V4 = V + dt*dV3; m4 = m + dt*dm3
h4 = h + dt*dh3; n4 = n + dt*dn3
dV4, dm4, dh4, dn4 = hh_rhs(V4, m4, h4, n4, I_total)
V = V + (dt/6.0)*(dV1+2*dV2+2*dV3+dV4)
m = min(max(m + (dt/6.0)*(dm1+2*dm2+2*dm3+dm4), 0.0), 1.0)
h = min(max(h + (dt/6.0)*(dh1+2*dh2+2*dh3+dh4), 0.0), 1.0)
n = min(max(n + (dt/6.0)*(dn1+2*dn2+2*dn3+dn4), 0.0), 1.0)
V_buf[ws % buf_size] = V
buf_idx = 0
# --- Main simulation with DFC ---
max_spikes = int(T_total_ms / 2) + 100
spike_times_raw = np.empty(max_spikes)
n_spikes_raw = 0
V_prev = V
# Store voltage trace (subsampled: every 10 steps = 0.1 ms)
subsample = 10
n_trace = n_steps // subsample + 1
V_trace = np.empty(n_trace)
trace_idx = 0
for step in range(n_steps):
# Read delayed voltage
V_delayed = V_buf[buf_idx]
# DFC control
I_ctrl = K * (V_delayed - V)
I_total = I_bias + I_ctrl
# RK4 with DFC recomputed at each substep
dV1, dm1, dh1, dn1 = hh_rhs(V, m, h, n, I_total)
Vk2 = V + 0.5*dt*dV1
I_ctrl_2 = K * (V_delayed - Vk2)
mk2 = m + 0.5*dt*dm1; hk2 = h + 0.5*dt*dh1; nk2 = n + 0.5*dt*dn1
dV2, dm2, dh2, dn2 = hh_rhs(Vk2, mk2, hk2, nk2, I_bias + I_ctrl_2)
Vk3 = V + 0.5*dt*dV2
I_ctrl_3 = K * (V_delayed - Vk3)
mk3 = m + 0.5*dt*dm2; hk3 = h + 0.5*dt*dh2; nk3 = n + 0.5*dt*dn2
dV3, dm3, dh3, dn3 = hh_rhs(Vk3, mk3, hk3, nk3, I_bias + I_ctrl_3)
Vk4 = V + dt*dV3
I_ctrl_4 = K * (V_delayed - Vk4)
mk4 = m + dt*dm3; hk4 = h + dt*dh3; nk4 = n + dt*dn3
dV4, dm4, dh4, dn4 = hh_rhs(Vk4, mk4, hk4, nk4, I_bias + I_ctrl_4)
V_new = V + (dt/6.0)*(dV1 + 2*dV2 + 2*dV3 + dV4)
m_new = min(max(m + (dt/6.0)*(dm1 + 2*dm2 + 2*dm3 + dm4), 0.0), 1.0)
h_new = min(max(h + (dt/6.0)*(dh1 + 2*dh2 + 2*dh3 + dh4), 0.0), 1.0)
n_new = min(max(n + (dt/6.0)*(dn1 + 2*dn2 + 2*dn3 + dn4), 0.0), 1.0)
# Update delay buffer
V_buf[buf_idx] = V_new
buf_idx = (buf_idx + 1) % buf_size
# Spike detection: rising threshold crossing at V = 0 mV
if V_prev <= 0.0 and V_new > 0.0:
if n_spikes_raw < max_spikes:
spike_times_raw[n_spikes_raw] = step * dt
n_spikes_raw += 1
# Subsample voltage
if step % subsample == 0 and trace_idx < n_trace:
V_trace[trace_idx] = V_new
trace_idx += 1
V_prev = V_new
V = V_new; m = m_new; h = h_new; n = n_new
# Extract spikes in analysis window
spike_times = spike_times_raw[:n_spikes_raw]
V_trace = V_trace[:trace_idx]
# Filter to analysis window
if T_transient_ms > 0:
mask = spike_times >= T_transient_ms
spike_times_analysis = spike_times[mask]
else:
spike_times_analysis = spike_times
if len(spike_times_analysis) >= 2:
isi_array = np.diff(spike_times_analysis)
else:
isi_array = np.empty(0)
return spike_times_analysis, isi_array, V_trace
# %% [markdown]
# ## CELL 4 — ISI Pattern Analysis Functions
# %%
@njit
def isi_cv(isi_array):
"""Coefficient of variation of ISI series."""
if len(isi_array) < 3:
return -1.0
mu = np.mean(isi_array)
if mu < 1e-10:
return -1.0
return np.std(isi_array) / mu
def detect_pattern_period(isis, max_period_spikes=12):
"""
Detect repeating ISI pattern via autocorrelation.
For periodic orbits with complex patterns (e.g., doublets, triplets),
the ISI series has a repeating pattern: [a, b, a, b, ...] for doublets,
[a, b, c, a, b, c, ...] for triplets, etc.
Returns: (pattern_length, pattern_isis, confidence)
pattern_length: number of ISIs in one pattern repeat (1=tonic)
pattern_isis: the repeating ISI sequence
confidence: normalized autocorrelation at pattern period
"""
if len(isis) < 6:
return 1, isis[:1] if len(isis) > 0 else np.array([0.0]), 0.0
isis = np.array(isis, dtype=np.float64)
n = len(isis)
# Try pattern lengths 1 to max_period_spikes
best_len = 1
best_conf = 0.0
best_pattern = isis[:1]
for p in range(1, min(max_period_spikes + 1, n // 3 + 1)):
# Check if ISI[i] ≈ ISI[i+p] for all i
n_compare = min(n - p, 3 * p) # Compare at least 3 repetitions
if n_compare < p:
continue
diffs = np.zeros(n_compare)
for i in range(n_compare):
diffs[i] = abs(isis[i] - isis[i + p])
# Median absolute deviation from pattern repeat
med_diff = np.median(diffs)
med_isi = np.median(isis)
if med_isi > 0:
relative_error = med_diff / med_isi
else:
continue
# Confidence: 1 - relative_error (capped at 0)
conf = max(0.0, 1.0 - relative_error * 10.0)
if conf > best_conf and conf > 0.8:
best_len = p
best_conf = conf
# Extract pattern: average over repetitions
pattern = np.zeros(p)
count = 0
for rep in range(n // p):
start = rep * p
if start + p <= n:
pattern += isis[start:start+p]
count += 1
if count > 0:
pattern /= count
best_pattern = pattern
# If best_len == 1, verify it's truly tonic (not just defaulting)
if best_len == 1 and len(isis) >= 3:
cv = np.std(isis) / np.mean(isis) if np.mean(isis) > 0 else 999
if cv < 0.02:
best_conf = 1.0
best_pattern = np.array([np.mean(isis)])
return best_len, best_pattern, best_conf
def classify_dynamics(isis, spike_times, V_trace, T_analysis_ms):
"""
Classify the dynamics at a (K, τ) point.
Returns: dict with classification and fingerprint
"""
result = {
'class': 'unknown',
'n_spikes': len(spike_times),
'isi_cv': -1.0,
'isi_mean': 0.0,
'isi_std': 0.0,
'firing_rate': 0.0,
'pattern_length': 0,
'pattern_isis': [],
'pattern_period_ms': 0.0,
'fingerprint': np.zeros(6),
}
# Silent
if len(spike_times) < 5:
result['class'] = 'silent'
return result
# Check for depolarization block (V stays high)
if V_trace is not None and len(V_trace) > 100:
last_quarter = V_trace[3*len(V_trace)//4:]
if np.mean(last_quarter) > 30.0 and np.std(last_quarter) < 5.0:
result['class'] = 'depol_block'
return result
cv = isi_cv(isis)
result['isi_cv'] = cv
result['isi_mean'] = float(np.mean(isis))
result['isi_std'] = float(np.std(isis))
result['firing_rate'] = len(spike_times) / (T_analysis_ms / 1000.0)
if cv < 0:
result['class'] = 'insufficient'
return result
# Classify by ISI variability
if cv < 0.02:
# Low variability: tonic or complex periodic
p_len, p_isis, p_conf = detect_pattern_period(isis)
result['pattern_length'] = p_len
result['pattern_isis'] = p_isis.tolist()
result['pattern_period_ms'] = float(np.sum(p_isis))
if p_len == 1:
result['class'] = 'tonic'
else:
result['class'] = f'periodic_p{p_len}'
elif cv < 0.10:
# Medium variability: could be complex periodic with jitter
p_len, p_isis, p_conf = detect_pattern_period(isis)
result['pattern_length'] = p_len
result['pattern_isis'] = p_isis.tolist()
result['pattern_period_ms'] = float(np.sum(p_isis))
if p_conf > 0.7:
result['class'] = f'periodic_p{p_len}'
else:
result['class'] = 'quasi_periodic'
elif cv < 0.15:
result['class'] = 'quasi_periodic'
else:
result['class'] = 'chaotic'
# Build fingerprint vector for clustering
# [mean_ISI, ISI_CV, pattern_length, pattern_period, firing_rate, ISI_range]
isi_range = float(np.max(isis) - np.min(isis)) if len(isis) > 0 else 0.0
result['fingerprint'] = np.array([
result['isi_mean'],
result['isi_cv'],
float(result['pattern_length']),
result['pattern_period_ms'],
result['firing_rate'],
isi_range,
])
return result
# %% [markdown]
# ## CELL 5 — Dense Orbit Sweep
# %%
def run_dense_sweep(I_bias, K_range, tau_range,
T_sim=3000.0, T_transient=500.0):
"""
Dense sweep of (K, τ) plane at fixed I_bias.
Parameters
----------
I_bias : float — Fixed bias current (μA/cm²)
K_range : array — DFC gain values
tau_range : array — DFC delay values (ms)
T_sim : float — Total simulation time (ms)
T_transient : float — Transient to discard (ms)
Returns
-------
results_grid : 2D array of classification dicts
summary : dict with aggregate statistics
"""
n_K = len(K_range)
n_tau = len(tau_range)
total = n_K * n_tau
T_analysis = T_sim - T_transient
print(f"Dense Sweep: {n_K} × {n_tau} = {total} grid points")
print(f" I_bias = {I_bias:.1f} μA/cm²")
print(f" K: [{K_range[0]:.3f}, {K_range[-1]:.3f}] ({n_K} values)")
print(f" τ: [{tau_range[0]:.1f}, {tau_range[-1]:.1f}] ms ({n_tau} values)")
print(f" T_sim = {T_sim:.0f} ms, T_transient = {T_transient:.0f} ms")
# JIT warmup
print("JIT warmup...", end=" ", flush=True)
_ = simulate_hh_dfc_full(10.0, 0.5, 50.0, 200.0, DT, 100.0)
print("done.")
# Storage
results_grid = [[None]*n_tau for _ in range(n_K)]
class_grid = np.full((n_K, n_tau), -1, dtype=np.int32)
isi_mean_grid = np.full((n_K, n_tau), np.nan)
cv_grid = np.full((n_K, n_tau), np.nan)
fr_grid = np.full((n_K, n_tau), np.nan)
# Class encoding
class_map = {
'silent': 0, 'depol_block': 1, 'tonic': 2,
'quasi_periodic': 3, 'chaotic': 4, 'insufficient': 5,
'unknown': 6,
}
# Periodic patterns: 10 + pattern_length
# e.g., periodic_p2 → 12, periodic_p3 → 13
t0 = time.time()
done = 0
for ik, K_val in enumerate(K_range):
for it, tau_val in enumerate(tau_range):
spk, isis, V_trace = simulate_hh_dfc_full(
I_bias, K_val, tau_val, T_sim, DT, T_transient)
result = classify_dynamics(isis, spk, V_trace, T_analysis)
result['K'] = float(K_val)
result['tau_ms'] = float(tau_val)
results_grid[ik][it] = result
# Encode class for grid
cls = result['class']
if cls in class_map:
class_grid[ik, it] = class_map[cls]
elif cls.startswith('periodic_p'):
p = int(cls.split('_p')[1])
class_grid[ik, it] = 10 + p
else:
class_grid[ik, it] = 6
isi_mean_grid[ik, it] = result['isi_mean']
cv_grid[ik, it] = result['isi_cv']
fr_grid[ik, it] = result['firing_rate']
done += 1
if done % 500 == 0:
elapsed = time.time() - t0
rate = done / elapsed
eta = (total - done) / rate
print(f" [{done}/{total}] {elapsed:.0f}s elapsed, "
f"~{eta:.0f}s remaining", flush=True)
elapsed = time.time() - t0
print(f"Dense sweep complete in {elapsed:.1f}s ({elapsed/60:.1f} min)")
# --- Summary statistics ---
class_counts = Counter()
for ik in range(n_K):
for it in range(n_tau):
cls = results_grid[ik][it]['class']
class_counts[cls] += 1
print(f"\nClassification summary:")
for cls, count in sorted(class_counts.items(), key=lambda x: -x[1]):
pct = 100.0 * count / total
print(f" {cls:20s}: {count:5d} ({pct:5.1f}%)")
summary = {
'I_bias': I_bias,
'n_K': n_K,
'n_tau': n_tau,
'total_points': total,
'T_sim': T_sim,
'T_transient': T_transient,
'class_counts': dict(class_counts),
'elapsed_s': elapsed,
}
return results_grid, class_grid, isi_mean_grid, cv_grid, fr_grid, summary
# %% [markdown]
# ## CELL 6 — Orbit Fingerprint Clustering
# %%
def cluster_orbits(results_grid, K_range, tau_range,
min_isi_separation=2.0, linkage_method='average'):
"""
Cluster periodic orbits by ISI fingerprint.
For all periodic points, compute pairwise distances in fingerprint
space and cluster using agglomerative hierarchical clustering.
Parameters
----------
min_isi_separation : float — minimum ISI mean difference (ms) between
distinct orbit types
Returns
-------
orbit_types : list of dicts, one per distinct orbit type
orbit_assignments : dict mapping (ik, it) → type_id
"""
n_K = len(K_range)
n_tau = len(tau_range)
# Collect all periodic points
periodic_points = []
for ik in range(n_K):
for it in range(n_tau):
r = results_grid[ik][it]
cls = r['class']
if cls == 'tonic' or cls.startswith('periodic_p'):
periodic_points.append({
'ik': ik, 'it': it,
'K': r['K'], 'tau_ms': r['tau_ms'],
'class': cls,
'isi_mean': r['isi_mean'],
'isi_cv': r['isi_cv'],
'pattern_length': r['pattern_length'],
'pattern_period_ms': r['pattern_period_ms'],
'pattern_isis': r['pattern_isis'],
'firing_rate': r['firing_rate'],
'fingerprint': r['fingerprint'],
})
n_periodic = len(periodic_points)
print(f"\nOrbit clustering: {n_periodic} periodic points")
if n_periodic == 0:
print(" No periodic points found!")
return [], {}
# Build fingerprint matrix
# Use: [isi_mean, pattern_period, pattern_length]
# Weight ISI mean most heavily (primary distinguishing feature)
fp_matrix = np.zeros((n_periodic, 3))
for i, p in enumerate(periodic_points):
fp_matrix[i, 0] = p['isi_mean'] # weight 1.0
fp_matrix[i, 1] = p['pattern_period_ms'] * 0.5 # weight 0.5
fp_matrix[i, 2] = p['pattern_length'] * 5.0 # weight 5.0 (categorical)
# Pairwise distances
if n_periodic > 1:
dists = pdist(fp_matrix, metric='euclidean')
Z = linkage(dists, method=linkage_method)
# Cut at threshold such that clusters differ by at least min_isi_separation
labels = fcluster(Z, t=min_isi_separation, criterion='distance')
else:
labels = np.array([1])
n_clusters = len(set(labels))
print(f" Found {n_clusters} orbit clusters (threshold={min_isi_separation} ms)")
# Build orbit type catalogue
orbit_types = []
orbit_assignments = {}
for cid in sorted(set(labels)):
members = [periodic_points[i] for i in range(n_periodic)
if labels[i] == cid]
# Representative: point closest to cluster centroid in ISI mean
isi_means = [m['isi_mean'] for m in members]
centroid_isi = np.mean(isi_means)
best_idx = np.argmin([abs(m['isi_mean'] - centroid_isi) for m in members])
rep = members[best_idx]
# Determine qualitative category
p_len = rep['pattern_length']
if p_len == 1:
category = 'tonic'
elif p_len == 2:
category = 'doublet'
elif p_len == 3:
category = 'triplet'
else:
category = f'burst_p{p_len}'
orbit_type = {
'type_id': int(cid),
'category': category,
'pattern_length': p_len,
'isi_mean': float(centroid_isi),
'isi_std': float(np.std(isi_means)),
'pattern_isis': rep['pattern_isis'],
'pattern_period_ms': rep['pattern_period_ms'],
'firing_rate': rep['firing_rate'],
'n_members': len(members),
'representative_K': rep['K'],
'representative_tau': rep['tau_ms'],
'K_range': [float(min(m['K'] for m in members)),
float(max(m['K'] for m in members))],
'tau_range': [float(min(m['tau_ms'] for m in members)),
float(max(m['tau_ms'] for m in members))],
}
orbit_types.append(orbit_type)
# Record assignments
for m in members:
orbit_assignments[(m['ik'], m['it'])] = int(cid)
# Sort by ISI mean
orbit_types.sort(key=lambda x: x['isi_mean'])
# Print catalogue
print(f"\n{'='*80}")
print(f"ORBIT TYPE CATALOGUE — {len(orbit_types)} distinct types")
print(f"{'='*80}")
print(f"{'ID':>4} {'Category':>10} {'PLen':>5} {'ISI_mean':>9} "
f"{'ISI_std':>8} {'Period':>8} {'FR(Hz)':>8} {'Members':>8} "
f"{'K*':>6} {'τ*':>8}")
print("-"*80)
for ot in orbit_types:
print(f"{ot['type_id']:>4} {ot['category']:>10} "
f"{ot['pattern_length']:>5} "
f"{ot['isi_mean']:>9.2f} {ot['isi_std']:>8.3f} "
f"{ot['pattern_period_ms']:>8.2f} {ot['firing_rate']:>8.1f} "
f"{ot['n_members']:>8} {ot['representative_K']:>6.3f} "
f"{ot['representative_tau']:>8.1f}")
# Count categories
cats = Counter(ot['category'] for ot in orbit_types)
print(f"\nQualitative categories: {dict(cats)}")
return orbit_types, orbit_assignments
# %% [markdown]
# ## CELL 7 — Gate PS-G0 Evaluation
# %%
def evaluate_gate_PS_G0(orbit_types, min_types=15, min_categories=3):
"""
Gate PS-G0: Enough orbit diversity for meaningful memory?
Requirements:
- ≥15 distinct orbit types
- Pairwise ISI mean difference > 2 ms for distinguishable pairs
- At least 3 qualitative categories
"""
print("\n" + "="*70)
print("GATE PS-G0 — ORBIT DIVERSITY ASSESSMENT")
print("="*70)
n_types = len(orbit_types)
categories = set(ot['category'] for ot in orbit_types)
n_categories = len(categories)
# Pairwise ISI separability
n_separable_pairs = 0
n_total_pairs = 0
min_separation = float('inf')
for i in range(n_types):
for j in range(i+1, n_types):
diff = abs(orbit_types[i]['isi_mean'] - orbit_types[j]['isi_mean'])
n_total_pairs += 1
if diff > 2.0:
n_separable_pairs += 1
min_separation = min(min_separation, diff)
# Gate checks
gate_types = n_types >= min_types
gate_categories = n_categories >= min_categories
gate_all = gate_types and gate_categories
print(f"\n Orbit types found: {n_types} (gate: ≥{min_types})"
f" {'✓ PASS' if gate_types else '✗ FAIL'}")
print(f" Categories: {n_categories} {dict(Counter(ot['category'] for ot in orbit_types))}"
f" (gate: ≥{min_categories})"
f" {'✓ PASS' if gate_categories else '✗ FAIL'}")
print(f" Separable pairs (Δ ISI > 2ms): {n_separable_pairs}/{n_total_pairs}")
print(f" Minimum ISI separation: {min_separation:.3f} ms")
if gate_all:
print(f"\n ✓✓✓ GATE PS-G0: PASS — Proceed to Phase PS1 ✓✓✓")
print(f" Orbit library has sufficient diversity for orbit-coded memory.")
else:
if not gate_types:
print(f"\n ✗ Insufficient orbit types. Try denser grid or secondary I_bias values.")
if not gate_categories:
print(f"\n ✗ Insufficient qualitative categories.")
print(f"\n ✗✗✗ GATE PS-G0: FAIL — Consider adjustments ✗✗✗")
gate_result = {
'decision': 'PASS' if gate_all else 'FAIL',
'n_types': n_types,
'n_categories': n_categories,
'categories': list(categories),
'gate_types': gate_types,
'gate_categories': gate_categories,
'n_separable_pairs': n_separable_pairs,
'n_total_pairs': n_total_pairs,
'min_separation_ms': float(min_separation) if min_separation < float('inf') else 0.0,
}
print("="*70)
return gate_result
# %% [markdown]
# ## CELL 8 — Visualization
# %%
def plot_bifurcation_map(class_grid, K_range, tau_range, I_bias, save_dir=None):
"""Bifurcation map: dynamics classification across (K, τ) plane."""
fig, ax = plt.subplots(1, 1, figsize=(14, 8))
# Custom colormap
# 0:silent(black), 1:depol(gray), 2:tonic(blue), 3:quasi(yellow),
# 4:chaotic(red), 5:insufficient(white), 6:unknown(white),
# 10+p: periodic_p (greens)
unique_classes = sorted(set(class_grid.flatten()))
colors_dict = {
0: '#1a1a2e', 1: '#4a4a4a', 2: '#2196F3', 3: '#FFC107',
4: '#F44336', 5: '#EEEEEE', 6: '#EEEEEE',
}
# Periodic patterns: gradient from green to purple
periodic_colors = ['#4CAF50', '#66BB6A', '#81C784', '#A5D6A7',
'#00BCD4', '#26C6DA', '#4DD0E1', '#80DEEA',
'#7E57C2', '#9575CD', '#B39DDB', '#CE93D8']
# Map class codes to sequential indices for colormapping
class_to_idx = {}
color_list = []
idx = 0
for c in unique_classes:
class_to_idx[c] = idx
if c in colors_dict:
color_list.append(colors_dict[c])
elif c >= 10:
p = c - 10
color_list.append(periodic_colors[min(p - 1, len(periodic_colors) - 1)])
else:
color_list.append('#CCCCCC')
idx += 1
# Remap grid
plot_grid = np.zeros_like(class_grid, dtype=float)
for ik in range(len(K_range)):
for it in range(len(tau_range)):
plot_grid[ik, it] = class_to_idx.get(class_grid[ik, it], 0)
cmap = ListedColormap(color_list)
im = ax.pcolormesh(tau_range, K_range, plot_grid,
cmap=cmap, shading='auto',
vmin=-0.5, vmax=len(color_list)-0.5)
ax.set_xlabel('τ (ms)', fontsize=13)
ax.set_ylabel('K (DFC gain)', fontsize=13)
ax.set_title(f'Bifurcation Map — HH + DFC at I_bias = {I_bias:.1f} μA/cm²',
fontsize=14)
ax.set_xscale('log')
# Legend
import matplotlib.patches as mpatches
class_names = {
0: 'Silent', 1: 'Depol. Block', 2: 'Tonic', 3: 'Quasi-periodic',
4: 'Chaotic', 5: 'Insufficient',
}
handles = []
for c in unique_classes:
if c in class_names:
name = class_names[c]
elif c >= 10:
name = f'Periodic p{c-10}'
else:
continue
color = color_list[class_to_idx[c]]
handles.append(mpatches.Patch(color=color, label=name))
ax.legend(handles=handles, loc='upper right', fontsize=9,
ncol=2, framealpha=0.9)
plt.tight_layout()
if save_dir:
plt.savefig(os.path.join(save_dir, 'PS0_bifurcation_map.png'),
dpi=300, bbox_inches='tight')
plt.show()
def plot_isi_landscape(isi_mean_grid, cv_grid, K_range, tau_range,
I_bias, save_dir=None):
"""ISI mean and CV landscape plots."""
fig, axes = plt.subplots(1, 2, figsize=(18, 7))
# ISI mean
ax = axes[0]
isi_plot = np.copy(isi_mean_grid)
isi_plot[isi_plot <= 0] = np.nan
im = ax.pcolormesh(tau_range, K_range, isi_plot,
cmap='viridis', shading='auto')
ax.set_xlabel('τ (ms)', fontsize=12)
ax.set_ylabel('K', fontsize=12)
ax.set_title(f'Mean ISI (ms) — I_bias={I_bias:.1f}', fontsize=13)
ax.set_xscale('log')
plt.colorbar(im, ax=ax, label='Mean ISI (ms)')
# ISI CV
ax = axes[1]
cv_plot = np.copy(cv_grid)
cv_plot[cv_plot < 0] = np.nan
im = ax.pcolormesh(tau_range, K_range, cv_plot,
cmap='hot', shading='auto', vmin=0, vmax=0.5)
ax.set_xlabel('τ (ms)', fontsize=12)
ax.set_ylabel('K', fontsize=12)
ax.set_title(f'ISI CV — I_bias={I_bias:.1f}', fontsize=13)
ax.set_xscale('log')
plt.colorbar(im, ax=ax, label='ISI CV')
plt.tight_layout()
if save_dir:
plt.savefig(os.path.join(save_dir, 'PS0_isi_landscape.png'),
dpi=300, bbox_inches='tight')
plt.show()
def plot_orbit_catalogue(orbit_types, save_dir=None):
"""Visualize the orbit type catalogue."""
if not orbit_types:
print("No orbit types to plot.")
return
n = len(orbit_types)
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
# 1. ISI mean distribution
ax = axes[0]
isi_means = [ot['isi_mean'] for ot in orbit_types]
categories = [ot['category'] for ot in orbit_types]
cat_colors = {'tonic': '#2196F3', 'doublet': '#4CAF50',
'triplet': '#FF9800', 'burst_p4': '#9C27B0'}
colors = [cat_colors.get(c, '#607D8B') for c in categories]
ax.barh(range(n), isi_means, color=colors, edgecolor='black', alpha=0.8)
ax.set_yticks(range(n))
ax.set_yticklabels([f"T{ot['type_id']} ({ot['category']})"
for ot in orbit_types], fontsize=8)
ax.set_xlabel('Mean ISI (ms)', fontsize=11)
ax.set_title('Orbit Types by ISI Mean', fontsize=12)
ax.invert_yaxis()
# 2. Pairwise separation matrix
ax = axes[1]
sep_matrix = np.zeros((n, n))
for i in range(n):
for j in range(n):
sep_matrix[i, j] = abs(isi_means[i] - isi_means[j])
im = ax.imshow(sep_matrix, cmap='YlOrRd', aspect='auto')
ax.set_xlabel('Orbit Type', fontsize=11)
ax.set_ylabel('Orbit Type', fontsize=11)
ax.set_title('Pairwise ISI Separation (ms)', fontsize=12)
plt.colorbar(im, ax=ax, label='|Δ ISI| (ms)')
# 3. Parameter space coverage
ax = axes[2]
for ot in orbit_types:
color = cat_colors.get(ot['category'], '#607D8B')
ax.scatter(ot['representative_tau'], ot['representative_K'],
c=color, s=80, edgecolors='black', linewidth=0.5, zorder=5)
ax.annotate(f"T{ot['type_id']}", (ot['representative_tau'], ot['representative_K']),
fontsize=7, ha='center', va='bottom')
ax.set_xlabel('τ (ms)', fontsize=11)
ax.set_ylabel('K', fontsize=11)
ax.set_title('Representative Points in (K, τ) Space', fontsize=12)
ax.set_xscale('log')
plt.tight_layout()
if save_dir:
plt.savefig(os.path.join(save_dir, 'PS0_orbit_catalogue.png'),
dpi=300, bbox_inches='tight')
plt.show()
# %% [markdown]
# ## CELL 9 — JSON Serialization Helper
# %%
def clean_for_json(obj):
"""Convert numpy types for JSON serialization."""
if isinstance(obj, dict):
return {k: clean_for_json(v) for k, v in obj.items()}
elif isinstance(obj, list):
return [clean_for_json(v) for v in obj]
elif isinstance(obj, (np.integer,)):
return int(obj)
elif isinstance(obj, (np.floating,)):
return float(obj)
elif isinstance(obj, np.ndarray):
return obj.tolist()
elif isinstance(obj, (np.bool_,)):
return bool(obj)
return obj
# %% [markdown]
# ## CELL 10 — MAIN EXECUTION
# %%
if __name__ == '__main__' or True: # Always run in Colab
print("="*70)
print("PS0 — DENSE ORBIT CATALOGUE")
print("Option C: HH Delay-Directed Orbit Selection")
print("="*70)
# ---- Define parameter grid ----
# Fixed I_bias from ExpF0v3 pilot (standard HH spiking regime)
I_BIAS = 10.0 # μA/cm²
# Dense (K, τ) grid — 10× denser than ExpF0v3
K_range = np.arange(0.0, 2.02, 0.02) # 101 values: 0.00, 0.02, ..., 2.00
tau_range = np.geomspace(1.0, 200.0, 100) # 100 values, log-spaced
# Include K=0 as baseline (no DFC → natural tonic firing)
print(f"\nGrid: {len(K_range)} × {len(tau_range)} "
f"= {len(K_range)*len(tau_range)} points")
print(f"I_bias = {I_BIAS:.1f} μA/cm²\n")
# ==== DENSE SWEEP ====
print("\n" + "="*50)
print("PHASE PS0-A — Dense Parameter Sweep")
print("="*50)
(results_grid, class_grid, isi_mean_grid, cv_grid, fr_grid,
summary) = run_dense_sweep(I_BIAS, K_range, tau_range,
T_sim=3000.0, T_transient=500.0)
# Save raw results
np.savez_compressed(os.path.join(OUTPUT_DIR, 'PS0_grid_data.npz'),
class_grid=class_grid,
isi_mean_grid=isi_mean_grid,
cv_grid=cv_grid,
fr_grid=fr_grid,
K_range=K_range,
tau_range=tau_range)
print(f"Grid data saved to {OUTPUT_DIR}")
# ==== VISUALIZATION ====
print("\n" + "="*50)
print("PHASE PS0-B — Bifurcation Analysis")
print("="*50)
plot_bifurcation_map(class_grid, K_range, tau_range, I_BIAS,
save_dir=OUTPUT_DIR)
plot_isi_landscape(isi_mean_grid, cv_grid, K_range, tau_range,
I_BIAS, save_dir=OUTPUT_DIR)
# ==== ORBIT CLUSTERING ====
print("\n" + "="*50)
print("PHASE PS0-C — Orbit Fingerprint Clustering")
print("="*50)
orbit_types, orbit_assignments = cluster_orbits(
results_grid, K_range, tau_range,
min_isi_separation=2.0)
plot_orbit_catalogue(orbit_types, save_dir=OUTPUT_DIR)
# Save orbit catalogue
with open(os.path.join(OUTPUT_DIR, 'PS0_orbit_types.json'), 'w') as f:
json.dump(clean_for_json(orbit_types), f, indent=2)
with open(os.path.join(OUTPUT_DIR, 'PS0_summary.json'), 'w') as f:
json.dump(clean_for_json(summary), f, indent=2)
print(f"\nOrbit catalogue saved to {OUTPUT_DIR}")
# ==== GATE PS-G0 ====
print("\n")
gate_result = evaluate_gate_PS_G0(orbit_types, min_types=15, min_categories=3)
with open(os.path.join(OUTPUT_DIR, 'gate_PS_G0_result.json'), 'w') as f:
json.dump(clean_for_json(gate_result), f, indent=2)
# ==== SECONDARY I_BIAS SWEEPS (if needed) ====
if gate_result['decision'] == 'FAIL':
print("\n" + "="*50)
print("SECONDARY SWEEPS — Additional I_bias values")
print("="*50)
print(" Primary sweep at I_bias=10.0 insufficient.")
print(" Running secondary sweeps at I_bias = 7.0, 8.5, 12.0, 15.0")
print(" (Set SECONDARY_SWEEPS = True and re-run this cell)")
# ==== FINAL SUMMARY ====
print(f"\n{'='*70}")
print(f"PS0 COMPLETE — ALL RESULTS SAVED")
print(f"{'='*70}")
print(f"Output directory: {OUTPUT_DIR}")
print(f"Files:")
print(f" PS0_grid_data.npz — Raw classification grid")
print(f" PS0_orbit_types.json — Orbit type catalogue")
print(f" PS0_summary.json — Sweep summary statistics")
print(f" gate_PS_G0_result.json — Gate decision")
print(f" PS0_bifurcation_map.png — (K, τ) dynamics heatmap")
print(f" PS0_isi_landscape.png — ISI mean and CV landscapes")
print(f" PS0_orbit_catalogue.png — Orbit types visualization")
print(f"\nGate PS-G0: {gate_result['decision']}")
if gate_result['decision'] == 'PASS':
print(f" → PROCEED TO PHASE PS1 (Write Protocol & Settling Time)")
else:
print(f" → Run secondary I_bias sweeps or adjust clustering threshold")

PS0_OrbitCatalogue.ipynb at commit 465b302, under MIT · at the source

Overview

Authors: Mohammad O Alhawarat¹, Ayman J Alnsour¹, Mohammed A F Al-Husainy¹, Khalil M Abdelnaby^1,2

Faculty of Information Technology, Al-Ahliyya Amman University, Amman 19328, Jordan; (A.J.A.); (M.A.F.A.-H.); (K.M.A.)
Systems and Computers Engineering Department, Faculty of Engineering, Al-Azhar University, Nasr City, Cairo 11765, Egypt

Journal: Entropy (Basel, Switzerland), volume 28, issue 6, article 678

Dates: received 17 April 2026; accepted 4 June 2026; published online 11 June 2026; in print June 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI · PMCID PMC13298313

Status: code verified

Categories: computational modeling (no new data) (modality), none (in silico) (organism)

Keywords: Hodgkin–Huxley model, parametric delayed feedback, limit-cycle multiplicity, orbit-coded memory, bifurcation landscape, autaptic feedback, neuromorphic computing

Citations: not cited yet (Europe PMC); 89 references in the paper

Abstract

We show that a single Hodgkin–Huxley (HH) neuron with Pyragas-type delayed feedback control (DFC) can store multiple symbols as stable periodic orbits, where the specific orbit is selected by tuning the DFC gain K and time delay τ. Sweeping the (K,τ) parameter plane at fixed bias current Ibias = 10.0 μA/cm2 reveals 207 orbit types across 12 topological categories, with inter-spike interval (ISI) means from 5.9 to 56.9 ms. We establish: (i) a write protocol that reliably locks orbits with 13.9 ms median settling time; (ii) a novel Pattern-Oriented Limit-cycle Decoder (POLD) that reads orbits at 100% accuracy from only five observed ISIs (1200 trials across 12 orbits; Wilson 95% CI: 99.7–100%); (iii) a complete single-symbol write–read–erase (W–R–E) cycle with 100% read accuracy, 92% erase verification, and no decay over hold durations up to 50 s; and (iv) a fully validated 12-symbol memory capacity with a read-discriminable upper bound of 67 symbols (11.2× over rate coding; write viability confirmed only for the conservative 12-symbol subset). Reliable orbit addressing needs delay precision of ±2%, which constitutes a write-precision specification and not a fundamental capacity limit. These findings show that parametric delayed feedback is a viable mechanism for limit-cycle-based information storage in conductance-based spiking neurons. The biological interpretation is analogical, not direct: the ±2% delay-precision requirement exceeds what has been demonstrated for biological autaptic variability, and the orbit-coded memory framing is best understood as a computational proof-of-principle aimed at neuromorphic engineering, not as a claim about biological working memory.

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

malhawarat/HH-DFC-OrbitMemory

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 465b302c206d7fcb197fc7c26ab2b35914991471, 4 July 2026

Languages: Jupyter (9)

Size: 11 files, 9 scripts

Software Heritage: not archived

Found in: “Supplementary Materials”

Holds: README, license file, 9 notebooks

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

Tools: Matplotlib (9 files), Numba (9 files), NumPy (9 files), SciPy (2 files)

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

27 September 2026: the link answers

11 files

experiments/PS0_OrbitCatalogue.ipynb — Jupyter, 1,038 lines, 8 matches
experiments/PS0_SecondaryIbias.ipynb — Jupyter, 752 lines, 2 matches
experiments/PS0b_FloquetValidation.ipynb — Jupyter, 843 lines, 3 matches
experiments/PS1_WriteProtocol.ipynb — Jupyter, 1,109 lines, 5 matches
experiments/PS2_ReadProtocol.ipynb — Jupyter, 1,183 lines, 5 matches
experiments/PS3_FullDemo_n_100.ipynb — Jupyter, 1,161 lines, 2 matches
experiments/PS3_FullDemo_n_20.ipynb — Jupyter, 1,148 lines, 2 matches
experiments/PS4_RateBaseline.ipynb — Jupyter, 605 lines, 4 matches
experiments/PS5_MaxCapacity.ipynb — Jupyter, 653 lines, 5 matches
LICENSE — License, 21 lines
README.md — Text, 147 lines

At the source: github.com/malhawarat/HH-DFC-OrbitMemory

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

All simulation code (six Jupyter notebooks, PS0–PS5), processed output data, and the complete orbit catalog are available in a public GitHub repository at https://github.com/malhawarat/HH-DFC-OrbitMemory (accessed on 1 June 2026). Raw simulation outputs exceeding the repository size limit are available from the corresponding author upon reasonable request.

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, 4 authors, 7 keywords, 80 references.

Cite

This paper

Alhawarat, M. O., Alnsour, A. J., Al-Husainy, M. A. F., & Abdelnaby, K. M. (2026). Limit-Cycle Proliferation Under Parametric Delayed Feedback in a Conductance-Based Neuron: Bifurcation Landscape, Orbit Catalog, and Capacity Analysis. Entropy (Basel, Switzerland), 28(6), 678.

BibTeX

@article{alhawarat2026limit, author = {Alhawarat, Mohammad O and Alnsour, Ayman J and Al-Husainy, Mohammed A F and Abdelnaby, Khalil M}, title = {{Limit-Cycle Proliferation Under Parametric Delayed Feedback in a Conductance-Based Neuron: Bifurcation Landscape, Orbit Catalog, and Capacity Analysis}}, journal = {Entropy (Basel, Switzerland)}, year = {2026}, month = jun, volume = {28}, number = {6}, pages = {678}, publisher = {Multidisciplinary Digital Publishing Institute (MDPI)}, issn = {1099-4300}, pmcid = {PMC13298313} }

RIS

TY - JOUR AU - Alhawarat, Mohammad O AU - Alnsour, Ayman J AU - Al-Husainy, Mohammed A F AU - Abdelnaby, Khalil M TI - Limit-Cycle Proliferation Under Parametric Delayed Feedback in a Conductance-Based Neuron: Bifurcation Landscape, Orbit Catalog, and Capacity Analysis T2 - Entropy (Basel, Switzerland) J2 - Entropy (Basel) PY - 2026 DA - 2026/06/01 VL - 28 IS - 6 SP - 678 SN - 1099-4300 PB - Multidisciplinary Digital Publishing Institute (MDPI) LA - en ER -

CSL-JSON

{ "id": "pmcid:PMC13298313", "type": "article-journal", "title": "Limit-Cycle Proliferation Under Parametric Delayed Feedback in a Conductance-Based Neuron: Bifurcation Landscape, Orbit Catalog, and Capacity Analysis", "container-title": "Entropy (Basel, Switzerland)", "author": [ { "family": "Alhawarat", "given": "Mohammad O" }, { "family": "Alnsour", "given": "Ayman J" }, { "family": "Al-Husainy", "given": "Mohammed A F" }, { "family": "Abdelnaby", "given": "Khalil M" } ], "container-title-short": "Entropy (Basel)", "volume": "28", "issue": "6", "page": "678", "PMCID": "PMC13298313", "ISSN": "1099-4300", "publisher": "Multidisciplinary Digital Publishing Institute (MDPI)", "language": "en", "issued": { "date-parts": [ [ 2026, 6, 1 ] ] } }

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

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.