codexintegrate-and-fire-current-steps.py
integrate-and-fire-current-steps.py
"""Exact LIF current-step trajectories; output basename.png in cwd.
Tested with Python 3.14.4, Matplotlib 3.10.7, NumPy 2.3.5.
Units in calculations: ms, mV, pA, MOhm. Parameters are illustrative.
Events use analytic crossing times, not time-grid threshold detection.
"""
import math
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
REST_MV = -65.0
RESET_MV = -65.0
THRESHOLD_MV = -50.0
TAU_MS = 20.0
RESISTANCE_MOHM = 100.0
REFRACTORY_MS = 2.0
STEP_START_MS = 20.0
STEP_END_MS = 180.0
END_MS = 220.0
def voltage(elapsed_ms, initial_mv, equilibrium_mv):
return equilibrium_mv + (initial_mv - equilibrium_mv) * np.exp(-elapsed_ms / TAU_MS)
def equilibrium(current_pa):
# 1 MOhm * 1 pA = 0.001 mV.
return REST_MV + RESISTANCE_MOHM * current_pa * 0.001
def active_segments(current_pa):
"""Return continuous analytic segments, event times and refractory holds."""
limit_mv = equilibrium(current_pa)
segments, events, holds = [], [], []
start_ms, initial_mv = STEP_START_MS, REST_MV
while start_ms < STEP_END_MS:
crossing_ms = math.inf
if limit_mv > THRESHOLD_MV:
crossing_ms = start_ms + TAU_MS * math.log(
(limit_mv - initial_mv) / (limit_mv - THRESHOLD_MV))
end_ms = min(crossing_ms, STEP_END_MS)
times = np.linspace(start_ms, end_ms, 300)
values = voltage(times - start_ms, initial_mv, limit_mv)
segments.append((times, values))
if crossing_ms >= STEP_END_MS:
return segments, events, holds, float(values[-1])
events.append(crossing_ms)
hold_end_ms = min(crossing_ms + REFRACTORY_MS, STEP_END_MS)
holds.append((crossing_ms, hold_end_ms))
start_ms, initial_mv = hold_end_ms, RESET_MV
return segments, events, holds, RESET_MV
def main():
high_limit = equilibrium(200.0)
latency_ms = TAU_MS * math.log((high_limit - RESET_MV) / (high_limit - THRESHOLD_MV))
period_ms = latency_ms + REFRACTORY_MS
rate_hz = 1000.0 / period_ms
low = active_segments(100.0)
high = active_segments(200.0)
assert equilibrium(100.0) == -55.0 and high_limit == -45.0
assert not low[1] and len(high[1]) == 5
assert math.isclose(latency_ms, 20.0 * math.log(4.0), abs_tol=1e-12)
assert math.isclose(high[1][0] - STEP_START_MS, latency_ms, abs_tol=1e-12)
assert np.allclose(np.diff(high[1]), period_ms, rtol=0, atol=1e-12)
assert math.isclose(rate_hz, 33.6407, rel_tol=0, abs_tol=0.00005)
for times, values in high[0][:-1]:
assert math.isclose(float(values[-1]), THRESHOLD_MV, abs_tol=1e-12)
assert all(np.max(values) < THRESHOLD_MV for _, values in low[0])
plt.rcParams.update({"font.family": "DejaVu Sans", "font.size": 10})
fig, axes = plt.subplots(3, 1, figsize=(9.0, 7.0), sharex=True,
gridspec_kw={"height_ratios": [0.8, 1.25, 1.65]})
fig.set_facecolor("white")
blue, orange = "#256f95", "#ba5418"
for ax in axes:
ax.set_facecolor("white")
ax.spines[["top", "right"]].set_visible(False)
ax.grid(axis="y", color="#e7e7e7", linewidth=0.7)
ax.set_xlim(0, END_MS)
for current, color in ((100.0, blue), (200.0, orange)):
axes[0].plot([0, STEP_START_MS, STEP_START_MS, STEP_END_MS, STEP_END_MS, END_MS],
[0, 0, current, current, 0, 0], color=color, lw=2,
label=f"{current:.0f} pA")
axes[0].set(ylabel="Input I (pA)", ylim=(-12, 245), yticks=[0, 100, 200])
axes[0].legend(loc="upper right", frameon=False, ncol=2)
axes[0].set_title("Current applied from 20 to 180 ms", loc="left", fontsize=11)
for ax, current, color, data in zip(axes[1:], (100.0, 200.0), (blue, orange), (low, high)):
segments, events, holds, final_mv = data
ax.plot([0, STEP_START_MS], [REST_MV, REST_MV], color=color, lw=2)
for times, values in segments:
ax.plot(times, values, color=color, lw=2)
for event in events:
ax.plot([event, event], [THRESHOLD_MV, RESET_MV], color=color, lw=1, ls=":")
ax.plot(event, THRESHOLD_MV + 2.0, marker="|", ms=10, color=color, mew=2)
for start_ms, end_ms in holds:
ax.axvspan(start_ms, end_ms, color="#ece1d8", alpha=1.0, zorder=0)
ax.plot([start_ms, end_ms], [RESET_MV, RESET_MV], color=color, lw=2)
times = np.linspace(STEP_END_MS, END_MS, 250)
ax.plot(times, voltage(times - STEP_END_MS, final_mv, REST_MV), color=color, lw=2)
ax.axhline(THRESHOLD_MV, color="#555555", ls="--", lw=1)
ax.hlines(equilibrium(current), STEP_START_MS, STEP_END_MS, color=color, ls="--", lw=1)
ax.set(ylabel="V (mV)", ylim=(-68, -42), yticks=[-65, -55, -50, -45])
axes[1].set_title("100 pA: equilibrium −55 mV, no events", loc="left", fontsize=11)
axes[1].text(214, THRESHOLD_MV + 0.7, "threshold −50", ha="right", fontsize=9)
axes[2].set_title("200 pA: event ticks, dotted resets and shaded refractory holds", loc="left", fontsize=11)
axes[2].text(174, -44.2, "passive equilibrium −45 mV", ha="right", fontsize=9, color=orange)
first_event = high[1][0]
axes[2].annotate("first event\n27.7259 ms after onset", xy=(first_event, -50),
xytext=(2, -45.0), fontsize=9, va="top",
arrowprops={"arrowstyle": "->", "color": "#444444"})
axes[2].annotate("2 ms hold", xy=(first_event + 1, RESET_MV), xytext=(28, -61),
ha="center", fontsize=9,
arrowprops={"arrowstyle": "->", "color": "#444444"})
axes[2].set_xlabel("Time (ms)")
fig.suptitle("Leaky integrate-and-fire: two exact current-step responses", fontsize=14)
fig.text(0.5, 0.062,
"EL = reset = −65 mV • threshold = −50 mV • τ = 20 ms • R = 100 MΩ",
ha="center", fontsize=10)
fig.text(0.5, 0.034, f"Sustained 200 pA: period = {period_ms:.4f} ms, rate = {rate_hz:.4f} Hz",
ha="center", fontsize=10)
fig.text(0.5, 0.010, "Illustrative model parameters; analytic trajectories; event ticks are not spike waveforms",
ha="center", fontsize=9)
fig.tight_layout(rect=(0.02, 0.085, 0.99, 0.95), h_pad=1.1)
fig.savefig(Path(__file__).stem + ".png", dpi=130, facecolor="white", transparent=False)
plt.close(fig)
print(f"Equilibria: {equilibrium(100.0):.1f}, {high_limit:.1f} mV")
print(f"Latency: {latency_ms:.10f} ms; period: {period_ms:.10f} ms; rate: {rate_hz:.10f} Hz")
print("200 pA event times (ms):", ", ".join(f"{t:.6f}" for t in high[1]))
print("100 pA events: 0; C = 200 pF; gL = 10 nS; analytic checks passed")
if __name__ == "__main__":
main()