Tutorial 2 - The Leaky Integrate-and-Fire Neuron — snntorch 0.9.1 documentation (2024)

2.1 Spiking Neurons: Intuition

In our brains, a neuron might be connected to 1,000 \(-\) 10,000other neurons. If one neuron spikes, all downhill neurons mightfeel it. But what determines whether a neuron spikes in the first place?The past century of experiments demonstrate that if a neuron experiencessufficient stimulus at its input, then it might become excited and fire its own spike.

Where does this stimulus come from? It could be from:

• the sensory periphery,

• an invasive electrode artificially stimulating the neuron, or in most cases,

• from other pre-synaptic neurons.

Given that these spikes are very short bursts of electrical activity, itis quite unlikely for all input spikes to arrive at the neuron body inprecise unison. This indicates the presence of temporal dynamics that‘sustain’ the input spikes, kind of like a delay.

2.2 The Passive Membrane

Like all cells, a neuron is surrounded by a thin membrane. This membraneis a lipid bilayer that insulates the conductive saline solution withinthe neuron from the extracellular medium. Electrically, the twoconductive solutions separated by an insulator act as a capacitor.

Another function of this membrane is to control what goes in and out ofthis cell (e.g., ions such as Na\(^+\)). The membrane is usuallyimpermeable to ions which blocks them from entering and exiting theneuron body. But there are specific channels in the membrane that aretriggered to open by injecting current into the neuron. This chargemovement is electrically modelled by a resistor.

The following block will derive the behaviour of a LIF neuron fromscratch. If you’d prefer to skip the math, then feel free to scroll onby; we’ll take a more hands-on approach to understanding the LIF neurondynamics after the derivation.

Optional: Derivation of LIF Neuron Model

Now say some arbitrary time-varying current \(I_{\rm in}(t)\) is injected into the neuron,be it via electrical stimulation or from other neurons. The total current in the circuit is conserved, so:

\[I_{\rm in}(t) = I_{R} + I_{C}\]

From Ohm’s Law, the membrane potential measured between the insideand outside of the neuron \(U_{\rm mem}\) is proportional tothe current through the resistor:

\[I_{R}(t) = \frac{U_{\rm mem}(t)}{R}\]

The capacitance is a proportionality constant between the chargestored on the capacitor \(Q\) and \(U_{\rm mem}(t)\):

\[Q = CU_{\rm mem}(t)\]

The rate of change of charge gives the capacitive current:

\[\frac{dQ}{dt}=I_C(t) = C\frac{dU_{\rm mem}(t)}{dt}\]

Therefore:

\[I_{\rm in}(t) = \frac{U_{\rm mem}(t)}{R} + C\frac{dU_{\rm mem}(t)}{dt}\]

\[\implies RC \frac{dU_{\rm mem}(t)}{dt} = -U_{\rm mem}(t) + RI_{\rm in}(t)\]

The right hand side of the equation is of units[Voltage]. On the left hand side of the equation,the term \(\frac{dU_{\rm mem}(t)}{dt}\) is of units[Voltage/Time]. To equate it to the left hand side (i.e., voltage),\(RC\) must be of unit [Time]. We refer to \(\tau = RC\) as the time constant of the circuit:

\[\tau \frac{dU_{\rm mem}(t)}{dt} = -U_{\rm mem}(t) + RI_{\rm in}(t)\]

The passive membrane is therefore described by a linear differential equation.

For a derivative of a function to be of the same form as the original function,i.e., \(\frac{dU_{\rm mem}(t)}{dt} \propto U_{\rm mem}(t)\), this impliesthe solution is exponential with a time constant \(\tau\).

Say the neuron starts at some value \(U_{0}\) with no further input,i.e., \(I_{\rm in}(t)=0.\) The solution of the linear differential equation is:

\[U_{\rm mem}(t) = U_0e^{-\frac{t}{\tau}}\]

The general solution is shown below.

Optional: Forward Euler Method to Solving the LIF Neuron Model

We managed to find the analytical solution to the LIF neuron, but it isunclear how this might be useful in a neural network. This time,let’s instead use the forward Euler method to solve the previous linearordinary differential equation (ODE). This approach might seemarduous, but it gives us a discrete, recurrent representation of the LIFneuron. Once we reach this solution, it can be applied directly to a neuralnetwork. As before, the linear ODE describing the RC circuit is:

\[\tau \frac{dU(t)}{dt} = -U(t) + RI_{\rm in}(t)\]

The subscript from \(U(t)\) is omitted for simplicity.

First, let’s solve this derivative without taking the limit\(\Delta t \rightarrow 0\):

\[\tau \frac{U(t+\Delta t)-U(t)}{\Delta t} = -U(t) + RI_{\rm in}(t)\]

For a small enough \(\Delta t\), this gives a good enoughapproximation of continuous-time integration. Isolating the membrane atthe following time step gives:

\[U(t+\Delta t) = U(t) + \frac{\Delta t}{\tau}\big(-U(t) + RI_{\rm in}(t)\big)\]

The following function represents this equation:

def leaky_integrate_neuron(U, time_step=1e-3, I=0, R=5e7, C=1e-10): tau = R*C U = U + (time_step/tau)*(-U + I*R) return U

The default values are set to \(R=50 M\Omega\) and\(C=100pF\) (i.e., \(\tau=5ms\)). These are quiterealistic with respect to biological neurons.

Now loop through this function, iterating one time step at a time.The membrane potential is initialized at \(U=0.9 V\), with the assumption thatthere is no injected input current, \(I_{\rm in}=0 A\).The simulation is performed with a millisecond precision\(\Delta t=1\times 10^{-3}\)s.

num_steps = 100U = 0.9U_trace = [] # keeps a record of U for plottingfor step in range(num_steps): U_trace.append(U) U = leaky_integrate_neuron(U) # solve next step of Uplot_mem(U_trace, "Leaky Neuron Model")

This exponential decay seems to match what we expected!

3.1 Lapicque: Without Stimulus

Instantiate Lapicque’s neuron using the following line of code.R & C are modified to simpler values, while keeping the previous timeconstant of \(\tau=5\times10^{-3}\)s.

time_step = 1e-3R = 5C = 1e-3# leaky integrate and fire neuron, tau=5e-3lif1 = snn.Lapicque(R=R, C=C, time_step=time_step)

The neuron model is now stored in lif1. To use this neuron:

Inputs

• cur_in: each element of \(I_{\rm in}\) is sequentially passed as an input (0 for now)

• mem: the membrane potential, previously \(U[t]\), is also passed as input. Initialize it arbitrarily as \(U[0] = 0.9~V\).

Outputs

• spk_out: output spike \(S_{\rm out}[t+\Delta t]\) at the next time step (‘1’ if there is a spike; ‘0’ if there is no spike)

• mem: membrane potential \(U_{\rm mem}[t+\Delta t]\) at the next time step

These all need to be of type torch.Tensor.

# Initialize membrane, input, and outputmem = torch.ones(1) * 0.9 # U=0.9 at t=0cur_in = torch.zeros(num_steps, 1) # I=0 for all tspk_out = torch.zeros(1) # initialize output spikes

These values are only for the initial time step \(t=0\).To analyze the evolution of mem over time, create a list mem_rec to record these values at every time step.

# A list to store a recording of membrane potentialmem_rec = [mem]

Now it’s time to run a simulation! At each time step, mem isupdated and stored in mem_rec:

# pass updated value of mem and cur_in[step]=0 at every time stepfor step in range(num_steps): spk_out, mem = lif1(cur_in[step], mem) # Store recordings of membrane potential mem_rec.append(mem)# convert the list of tensors into one tensormem_rec = torch.stack(mem_rec)# pre-defined plotting functionplot_mem(mem_rec, "Lapicque's Neuron Model Without Stimulus")

The membrane potential decays over time in the absence of any inputstimuli.

3.2 Lapicque: Step Input

Now apply a step current \(I_{\rm in}(t)\) that switches on at\(t=t_0\). Given the linear first-order differential equation:

\[\tau \frac{dU_{\rm mem}}{dt} = -U_{\rm mem} + RI_{\rm in}(t),\]

the general solution is:

\[U_{\rm mem}=I_{\rm in}(t)R + [U_0 - I_{\rm in}(t)R]e^{-\frac{t}{\tau}}\]

If the membrane potential is initialized to\(U_{\rm mem}(t=0) = 0 V\), then:

\[U_{\rm mem}(t)=I_{\rm in}(t)R [1 - e^{-\frac{t}{\tau}}]\]

Based on this explicit time-dependent form, we expect\(U_{\rm mem}\) to relax exponentially towards \(I_{\rm in}R\).Let’s visualize what this looks like by triggering a current pulse of\(I_{in}=100mA\) at \(t_0 = 10ms\).

# Initialize input current pulsecur_in = torch.cat((torch.zeros(10, 1), torch.ones(190, 1)*0.1), 0) # input current turns on at t=10# Initialize membrane, output and recordingsmem = torch.zeros(1) # membrane potential of 0 at t=0spk_out = torch.zeros(1) # neuron needs somewhere to sequentially dump its output spikesmem_rec = [mem]

This time, the new values of cur_in are passed to the neuron:

num_steps = 200# pass updated value of mem and cur_in[step] at every time stepfor step in range(num_steps): spk_out, mem = lif1(cur_in[step], mem) mem_rec.append(mem)# crunch -list- of tensors into one tensormem_rec = torch.stack(mem_rec)plot_step_current_response(cur_in, mem_rec, 10)

As \(t\rightarrow \infty\), the membrane potential\(U_{\rm mem}\) exponentially relaxes to \(I_{\rm in}R\):

>>> print(f"The calculated value of input pulse [A] x resistance [Ω] is: {cur_in[11]*lif1.R} V")>>> print(f"The simulated value of steady-state membrane potential is: {mem_rec[200][0]} V")The calculated value of input pulse [A] x resistance [Ω] is: 0.5 VThe simulated value of steady-state membrane potential is: 0.4999999403953552 V

Close enough!

3.3 Lapicque: Pulse Input

Now what if the step input was clipped at \(t=30ms\)?

# Initialize current pulse, membrane and outputscur_in1 = torch.cat((torch.zeros(10, 1), torch.ones(20, 1)*(0.1), torch.zeros(170, 1)), 0) # input turns on at t=10, off at t=30mem = torch.zeros(1)spk_out = torch.zeros(1)mem_rec1 = [mem]

# neuron simulationfor step in range(num_steps): spk_out, mem = lif1(cur_in1[step], mem) mem_rec1.append(mem)mem_rec1 = torch.stack(mem_rec1)plot_current_pulse_response(cur_in1, mem_rec1, "Lapicque's Neuron Model With Input Pulse", vline1=10, vline2=30)

\(U_{\rm mem}\) rises just as it did for the step input, but now itdecays with a time constant of \(\tau\) as in our first simulation.

Let’s deliver approximately the same amount of charge\(Q = I \times t\) to the circuit in half the time. This means theinput current amplitude must be increased by a little, and thetime window must be decreased.

# Increase amplitude of current pulse; half the time.cur_in2 = torch.cat((torch.zeros(10, 1), torch.ones(10, 1)*0.111, torch.zeros(180, 1)), 0) # input turns on at t=10, off at t=20mem = torch.zeros(1)spk_out = torch.zeros(1)mem_rec2 = [mem]# neuron simulationfor step in range(num_steps): spk_out, mem = lif1(cur_in2[step], mem) mem_rec2.append(mem)mem_rec2 = torch.stack(mem_rec2)plot_current_pulse_response(cur_in2, mem_rec2, "Lapicque's Neuron Model With Input Pulse: x1/2 pulse width", vline1=10, vline2=20)

Let’s do that again, but with an even faster input pulse and higheramplitude:

# Increase amplitude of current pulse; quarter the time.cur_in3 = torch.cat((torch.zeros(10, 1), torch.ones(5, 1)*0.147, torch.zeros(185, 1)), 0) # input turns on at t=10, off at t=15mem = torch.zeros(1)spk_out = torch.zeros(1)mem_rec3 = [mem]# neuron simulationfor step in range(num_steps): spk_out, mem = lif1(cur_in3[step], mem) mem_rec3.append(mem)mem_rec3 = torch.stack(mem_rec3)plot_current_pulse_response(cur_in3, mem_rec3, "Lapicque's Neuron Model With Input Pulse: x1/4 pulse width", vline1=10, vline2=15)

Now compare all three experiments on the same plot:

compare_plots(cur_in1, cur_in2, cur_in3, mem_rec1, mem_rec2, mem_rec3, 10, 15, 20, 30, "Lapicque's Neuron Model With Input Pulse: Varying inputs")

As the input current pulse amplitude increases, the rise time of themembrane potential speeds up. In the limit of the input current pulsewidth becoming infinitesimally small, \(T_W \rightarrow 0s\), themembrane potential will jump straight up in virtually zero rise time:

# Current spike inputcur_in4 = torch.cat((torch.zeros(10, 1), torch.ones(1, 1)*0.5, torch.zeros(189, 1)), 0) # input only on for 1 time stepmem = torch.zeros(1)spk_out = torch.zeros(1)mem_rec4 = [mem]# neuron simulationfor step in range(num_steps): spk_out, mem = lif1(cur_in4[step], mem) mem_rec4.append(mem)mem_rec4 = torch.stack(mem_rec4)plot_current_pulse_response(cur_in4, mem_rec4, "Lapicque's Neuron Model With Input Spike", vline1=10, ylim_max1=0.6)

The current pulse width is now so short, it effectively looks like aspike. That is to say, charge is delivered in an infinitely short periodof time, \(I_{\rm in}(t) = Q/t_0\) where \(t_0 \rightarrow 0\).More formally:

\[I_{\rm in}(t) = Q \delta (t-t_0),\]

where \(\delta (t-t_0)\) is the Dirac-Delta function. Physically, itis impossible to ‘instantaneously’ deposit charge. But integrating\(I_{\rm in}\) gives a result that makes physical sense, as we canobtain the charge delivered:

\[1 = \int^{t_0 + a}_{t_0 - a}\delta(t-t_0)dt\]

\[f(t_0) = \int^{t_0 + a}_{t_0 - a}f(t)\delta(t-t_0)dt\]

Here,\(f(t_0) = I_{\rm in}(t_0=10) = 0.5A \implies f(t) = Q = 0.5C\).

Hopefully you have a good feel of how the membrane potential leaks atrest, and integrates the input current. That covers the ‘leaky’ and‘integrate’ part of the neuron. How about the fire?

3.4 Lapicque: Firing

So far, we have only seen how a neuron will react to spikes at theinput. For a neuron to generate and emit its own spikes at the output,the passive membrane model must be combined with a threshold.

If the membrane potential exceeds this threshold, then a voltage spikewill be generated, external to the passive membrane model.

Modify the leaky_integrate_neuron function from before to adda spike response.

# R=5.1, C=5e-3 for illustrative purposesdef leaky_integrate_and_fire(mem, cur=0, threshold=1, time_step=1e-3, R=5.1, C=5e-3): tau_mem = R*C spk = (mem > threshold) # if membrane exceeds threshold, spk=1, else, 0 mem = mem + (time_step/tau_mem)*(-mem + cur*R) return mem, spk

Set threshold=1, and apply a step current to get this neuronspiking.

# Small step current inputcur_in = torch.cat((torch.zeros(10), torch.ones(190)*0.2), 0)mem = torch.zeros(1)mem_rec = []spk_rec = []# neuron simulationfor step in range(num_steps): mem, spk = leaky_integrate_and_fire(mem, cur_in[step]) mem_rec.append(mem) spk_rec.append(spk)# convert lists to tensorsmem_rec = torch.stack(mem_rec)spk_rec = torch.stack(spk_rec)plot_cur_mem_spk(cur_in, mem_rec, spk_rec, thr_line=1, vline=109, ylim_max2=1.3, title="LIF Neuron Model With Uncontrolled Spiking")

Oops - the output spikes have gone out of control! This isbecause we forgot to add a reset mechanism. In reality, each time aneuron fires, the membrane potential hyperpolarizes back to its restingpotential.

Implementing this reset mechanism into our neuron:

# LIF w/Reset mechanismdef leaky_integrate_and_fire(mem, cur=0, threshold=1, time_step=1e-3, R=5.1, C=5e-3): tau_mem = R*C spk = (mem > threshold) mem = mem + (time_step/tau_mem)*(-mem + cur*R) - spk*threshold # every time spk=1, subtract the threhsold return mem, spk

# Small step current inputcur_in = torch.cat((torch.zeros(10), torch.ones(190)*0.2), 0)mem = torch.zeros(1)mem_rec = []spk_rec = []# neuron simulationfor step in range(num_steps): mem, spk = leaky_integrate_and_fire(mem, cur_in[step]) mem_rec.append(mem) spk_rec.append(spk)# convert lists to tensorsmem_rec = torch.stack(mem_rec)spk_rec = torch.stack(spk_rec)plot_cur_mem_spk(cur_in, mem_rec, spk_rec, thr_line=1, vline=109, ylim_max2=1.3, title="LIF Neuron Model With Reset")

Bam. We now have a functional leaky integrate-and-fire neuron model!

Note that if \(I_{\rm in}=0.2 A\) and \(R<5 \Omega\), then\(I\times R < 1 V\). If threshold=1, then no spiking wouldoccur. Feel free to go back up, change the values, and test it out.

As before, all of that code is condensed by calling the built-in Lapicque neuron model from snnTorch:

# Create the same neuron as before using snnTorchlif2 = snn.Lapicque(R=5.1, C=5e-3, time_step=1e-3)>>> print(f"Membrane potential time constant: {lif2.R * lif2.C:.3f}s")"Membrane potential time constant: 0.025s"

# Initialize inputs and outputscur_in = torch.cat((torch.zeros(10, 1), torch.ones(190, 1)*0.2), 0)mem = torch.zeros(1)spk_out = torch.zeros(1)mem_rec = [mem]spk_rec = [spk_out]# Simulation run across 100 time steps.for step in range(num_steps): spk_out, mem = lif2(cur_in[step], mem) mem_rec.append(mem) spk_rec.append(spk_out)# convert lists to tensorsmem_rec = torch.stack(mem_rec)spk_rec = torch.stack(spk_rec)plot_cur_mem_spk(cur_in, mem_rec, spk_rec, thr_line=1, vline=109, ylim_max2=1.3, title="Lapicque Neuron Model With Step Input")

The membrane potential exponentially rises and then hits the threshold,at which point it resets. We can roughly see this occurs between\(105ms < t_{\rm spk} < 115ms\). As a matter of curiousity, let’ssee what the spike recording actually consists of:

>>> print(spk_rec[105:115].view(-1))tensor([0., 0., 0., 0., 1., 0., 0., 0., 0., 0.])

The absence of a spike is represented by \(S_{\rm out}=0\), and theoccurrence of a spike is \(S_{\rm out}=1\). Here, the spike occursat \(S_{\rm out}[t=109]=1\). If you are wondering why each of these entries is stored as a tensor, itis because in future tutorials we will simulate large scale neuralnetworks. Each entry will contain the spike responses of many neurons,and tensors can be loaded into GPU memory to speed up the trainingprocess.

If \(I_{\rm in}\) is increased, then the membrane potentialapproaches the threshold \(U_{\rm thr}\) faster:

# Initialize inputs and outputscur_in = torch.cat((torch.zeros(10, 1), torch.ones(190, 1)*0.3), 0) # increased currentmem = torch.zeros(1)spk_out = torch.zeros(1)mem_rec = [mem]spk_rec = [spk_out]# neuron simulationfor step in range(num_steps): spk_out, mem = lif2(cur_in[step], mem) mem_rec.append(mem) spk_rec.append(spk_out)# convert lists to tensorsmem_rec = torch.stack(mem_rec)spk_rec = torch.stack(spk_rec)plot_cur_mem_spk(cur_in, mem_rec, spk_rec, thr_line=1, ylim_max2=1.3, title="Lapicque Neuron Model With Periodic Firing")

A similar increase in firing frequency can also be induced by decreasingthe threshold. This requires initializing a new neuron model, but therest of the code block is the exact same as above:

# neuron with halved thresholdlif3 = snn.Lapicque(R=5.1, C=5e-3, time_step=1e-3, threshold=0.5)# Initialize inputs and outputscur_in = torch.cat((torch.zeros(10, 1), torch.ones(190, 1)*0.3), 0)mem = torch.zeros(1)spk_out = torch.zeros(1)mem_rec = [mem]spk_rec = [spk_out]# Neuron simulationfor step in range(num_steps): spk_out, mem = lif3(cur_in[step], mem) mem_rec.append(mem) spk_rec.append(spk_out)# convert lists to tensorsmem_rec = torch.stack(mem_rec)spk_rec = torch.stack(spk_rec)plot_cur_mem_spk(cur_in, mem_rec, spk_rec, thr_line=0.5, ylim_max2=1.3, title="Lapicque Neuron Model With Lower Threshold")

That’s what happens for a constant current injection. But in both deepneural networks and in the biological brain, most neurons will beconnected to other neurons. They are more likely to receive spikes,rather than injections of constant current.

3.5 Lapicque: Spike Inputs

Let’s harness some of the skills we learnt in Tutorial1,and use the snntorch.spikegen module to create some randomlygenerated input spikes.

# Create a 1-D random spike train. Each element has a probability of 40% of firing.spk_in = spikegen.rate_conv(torch.ones((num_steps,1)) * 0.40)

Run the following code block to see how many spikes have been generated.

>>> print(f"There are {int(sum(spk_in))} total spikes out of {len(spk_in)} time steps.")There are 85 total spikes out of 200 time steps.

fig = plt.figure(facecolor="w", figsize=(8, 1))ax = fig.add_subplot(111)splt.raster(spk_in.reshape(num_steps, -1), ax, s=100, c="black", marker="|")plt.title("Input Spikes")plt.xlabel("Time step")plt.yticks([])plt.show()

# Initialize inputs and outputsmem = torch.ones(1)*0.5spk_out = torch.zeros(1)mem_rec = [mem]spk_rec = [spk_out]# Neuron simulationfor step in range(num_steps): spk_out, mem = lif3(spk_in[step], mem) spk_rec.append(spk_out) mem_rec.append(mem)# convert lists to tensorsmem_rec = torch.stack(mem_rec)spk_rec = torch.stack(spk_rec)plot_spk_mem_spk(spk_in, mem_rec, spk_out, "Lapicque's Neuron Model With Input Spikes")

3.6 Lapicque: Reset Mechanisms

We already implemented a reset mechanism from scratch, but let’s dive alittle deeper. This sharp drop of membrane potential promotes areduction of spike generation, which supplements part of the theory onhow brains are so power efficient. Biologically, this drop of membranepotential is known as ‘hyperpolarization’. Following that, it ismomentarily more difficult to elicit another spike from the neuron.Here, we use a reset mechanism to model hyperpolarization.

There are two ways to implement the reset mechanism:

1. reset by subtraction (default) \(-\) subtract the thresholdfrom the membrane potential each time a spike is generated;

2. reset to zero \(-\) force the membrane potential to zero eachtime a spike is generated.

3. no reset \(-\) do nothing, and let the firing go potentially uncontrolled.

Instantiate another neuron model to demonstrate how to alternatebetween reset mechanisms. By default, snnTorch neuron models use reset_mechanism = "subtract".This can be explicitly overridden by passing the argumentreset_mechanism = "zero".

# Neuron with reset_mechanism set to "zero"lif4 = snn.Lapicque(R=5.1, C=5e-3, time_step=1e-3, threshold=0.5, reset_mechanism="zero")# Initialize inputs and outputsspk_in = spikegen.rate_conv(torch.ones((num_steps, 1)) * 0.40)mem = torch.ones(1)*0.5spk_out = torch.zeros(1)mem_rec0 = [mem]spk_rec0 = [spk_out]# Neuron simulationfor step in range(num_steps): spk_out, mem = lif4(spk_in[step], mem) spk_rec0.append(spk_out) mem_rec0.append(mem)# convert lists to tensorsmem_rec0 = torch.stack(mem_rec0)spk_rec0 = torch.stack(spk_rec0)plot_reset_comparison(spk_in, mem_rec, spk_rec, mem_rec0, spk_rec0)

Pay close attention to the evolution of the membrane potential,especially in the moments after it reaches the threshold. You may noticethat for “Reset to Zero”, the membrane potential is forced back to zeroafter each spike.

So which one is better? Applying "subtract" (the default value inreset_mechanism) is less lossy, because it does not ignore how muchthe membrane exceeds the threshold by.

On the other hand, applying a hard reset with "zero" promotessparsity and potentially less power consumption when running ondedicated neuromorphic hardware. Both options are available for you toexperiment with.

That covers the basics of a LIF neuron model!