diff --git a/examples/lif_fi_curve.py b/examples/lif_fi_curve.py new file mode 100644 index 0000000..b05af77 --- /dev/null +++ b/examples/lif_fi_curve.py @@ -0,0 +1,153 @@ +# Copyright (c) 2017 The University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" +LIF neuron F-I curve demonstration. + +Five IF_curr_exp (LIF) neurons are each driven by a different constant DC +current. The example records membrane voltage and spikes then produces: + + 1. Voltage traces — showing integrate-to-threshold, spike, and reset for + each current level. + 2. F-I curve — measured firing rate vs. injected current overlaid with + the analytical LIF prediction. + +LIF parameters used (defaults for IF_curr_exp): + v_rest = -65 mV v_thresh = -50 mV v_reset = -65 mV + tau_m = 20 ms tau_refrac = 2 ms cm = 1 nF + +Threshold current: I_thresh = (v_thresh - v_rest) * cm / tau_m = 0.75 nA + - neurons driven below this level stay at rest (subthreshold) + - neurons driven above fire at a rate that increases with current +""" + +import numpy as np +import matplotlib.pyplot as plt +import pyNN.spiNNaker as sim + +# ── Simulation parameters ──────────────────────────────────────────────────── +SIMTIME = 1000.0 # ms — total simulation duration +TIMESTEP = 1.0 # ms +STIM_START = 100.0 # ms — DC current onset +STIM_STOP = 900.0 # ms — DC current offset + +# DC current levels to sweep (nA). Two below threshold, three above. +CURRENTS = [0.3, 0.6, 0.9, 1.2, 1.5] +N = len(CURRENTS) + +# LIF neuron parameters (explicit for clarity) +LIF_PARAMS = dict( + v_rest=-65.0, v_reset=-65.0, v_thresh=-50.0, + tau_m=20.0, tau_refrac=2.0, cm=1.0, + tau_syn_E=5.0, tau_syn_I=5.0, +) + +# ── Analytical F-I curve ───────────────────────────────────────────────────── + + +def lif_fi_theory(currents, v_rest=-65.0, v_reset=-65.0, v_thresh=-50.0, + tau_m=20.0, tau_refrac=2.0, cm=1.0): + """Return theoretical LIF firing rates (Hz) for each current (nA).""" + R_m = tau_m / cm # MΩ (tau_m in ms, cm in nF → mV/nA) + rates = [] + for I in currents: + V_inf = v_rest + R_m * I # steady-state voltage (mV) + if V_inf <= v_thresh: + rates.append(0.0) + else: + # time from v_reset to v_thresh + t_spike = tau_m * np.log( + (V_inf - v_reset) / (V_inf - v_thresh)) # ms + rates.append(1000.0 / (tau_refrac + t_spike)) # Hz + return np.array(rates) + + +# ── Build and run network ──────────────────────────────────────────────────── +sim.setup(timestep=TIMESTEP, min_delay=TIMESTEP) +sim.set_number_of_neurons_per_core(sim.IF_curr_exp, 100) + +neurons = sim.Population( + N, sim.IF_curr_exp(**LIF_PARAMS), label="LIF neurons") + +for cell, amp in zip(neurons, CURRENTS): + cell.inject(sim.DCSource(amplitude=amp, start=STIM_START, stop=STIM_STOP)) + +neurons.record(["v", "spikes"]) + +sim.run(SIMTIME) + +neo = neurons.get_data(variables=["v", "spikes"]) +v_data = neo.segments[0].filter(name="v")[0] # shape: (time, N) +spike_trains = neo.segments[0].spiketrains # list of N SpikeTrains + +sim.end() + +# ── Compute measured firing rates ──────────────────────────────────────────── +stim_duration_s = (STIM_STOP - STIM_START) / 1000.0 # seconds +measured_rates = [] +for train in spike_trains: + in_window = [s for s in train + if STIM_START < float(s) < STIM_STOP] + measured_rates.append(len(in_window) / stim_duration_s) + +theory_currents = np.linspace(0, max(CURRENTS) + 0.2, 200) +theory_rates = lif_fi_theory(theory_currents, **LIF_PARAMS) + +# ── Plot ───────────────────────────────────────────────────────────────────── +fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8)) +fig.suptitle("LIF Neuron — F-I Curve Demonstration", fontsize=13) + +colors = plt.cm.viridis(np.linspace(0.15, 0.85, N)) + +# Panel 1: voltage traces +for i in range(N): + ax1.plot(v_data.times, v_data[:, i], + color=colors[i], label=f"I = {CURRENTS[i]:.1f} nA") +ax1.axhline(LIF_PARAMS["v_thresh"], color="red", + linestyle="--", linewidth=1.0, label=f"Threshold ({LIF_PARAMS['v_thresh']} mV)") +ax1.axhline(LIF_PARAMS["v_rest"], color="gray", + linestyle=":", linewidth=0.8, label=f"V_rest ({LIF_PARAMS['v_rest']} mV)") +ax1.axvspan(STIM_START, STIM_STOP, alpha=0.08, color="yellow", label="Stimulus window") +ax1.set_xlabel("Time (ms)") +ax1.set_ylabel("Membrane potential (mV)") +ax1.set_title("Membrane Voltage Traces") +ax1.legend(loc="upper right", fontsize=8, ncol=2) +ax1.set_xlim(0, SIMTIME) + +# Panel 2: F-I curve +ax2.plot(theory_currents, theory_rates, + color="tomato", linewidth=1.5, linestyle="--", label="Analytical LIF") +ax2.plot(CURRENTS, measured_rates, + "o", color="steelblue", markersize=9, zorder=5, label="Simulated") +ax2.axvline(0.75, color="gray", linestyle=":", linewidth=1.0, + label=r"$I_{thresh}$ = 0.75 nA") +ax2.set_xlabel("Injected current (nA)") +ax2.set_ylabel("Firing rate (Hz)") +ax2.set_title("F-I Curve (Firing Rate vs. Injected Current)") +ax2.legend(fontsize=9) +ax2.set_xlim(0, max(CURRENTS) + 0.2) +ax2.set_ylim(bottom=0) +ax2.grid(True, alpha=0.3) + +plt.tight_layout() +plt.savefig("lif_fi_curve.png", dpi=150, bbox_inches="tight") +plt.show() + +# ── Console summary ────────────────────────────────────────────────────────── +print(f"\n{'Current (nA)':<14} {'Measured (Hz)':<16} {'Theory (Hz)':<14} Status") +print("-" * 58) +theory_at_points = lif_fi_theory(CURRENTS, **LIF_PARAMS) +for I, meas, theo in zip(CURRENTS, measured_rates, theory_at_points): + status = "subthreshold" if meas == 0 else "firing" + print(f"{I:<14.1f} {meas:<16.1f} {theo:<14.1f} {status}")