Introduction to Neural Networks

Programming Lab 2: The Integrate-and-Fire (LIF) Neuron

Getting started

As before, we start by importing the libraries we will use:

import numpy as np 
import matplotlib.pyplot as plt 
from itertools import islice # import this to slice time within the "for" loop

The LIF Neuron

The LIF neuron represents one of the simplest descriptions of spiking neurons. Developed in 1907 by Pierre Lapicque, the LIF model predated the modern understanding of the ionic basis of the action potential. While more complex models can better describe the spiking activity of single neurons (e.g. the Hodgkin-Huxley model, for which the Nobel prize was awarded in 1963), the LIF neuron has proven highly useful in describing the dynamics of large networks. Because neurons in cortex receive strong synaptic input from the surrounding network, the simplicity of the LIF neuron represents a key path to explain neuronal responses, going from the most basic current inputs to the complex fluctuations created by thousands of synaptic contacts in cortex.

The LIF model is defined by the equation:

\begin{equation} \tau_m \frac{dV}{dt} = E_L - V + R_m I_e \tag{1}\label{lif} \end{equation}

where \(\tau_m\) is the membrane time constant, \(V\) is the membrane potential, \(E_L\) is the resting potential, \(R_m\) is the membrane resistance, and \(C_m\) is the membrane capacitance. Remember that \(C_m\) determines how much current it will take to drive the membrane potential to change at a certain rate, while the membrane resistance \(R_m\) determines the current required to hold the neuron away from \(E_L\) at steady state. The relationship of the neuronal membrane, ion channels, and these parameters we discussed in class is depicted below for the potassium conductance:

Typical parameter values for neurons in cortex are:

# parameters
Rm     = 1e6    # resistance (ohm)
Cm     = 2e-8   # capacitance (farad)
taum   = Rm*Cm  # time constant (seconds)
Vr     = -.060  # resting membrane potential (volt)
Vreset = -.070  # membrane potential after spike (volt)
Vth    = -.050  # spike threshold (volt)
Vs     = .020   # spiking potential (volt)

Time axis

Next, we need to define a "time axis", the discrete points at which we will evaluate the simulation:

dt     = .001   # simulation time step (seconds)
T      = 1.0    # total time to simulate (seconds)
time   = np.linspace(dt,T,T/dt) # vector of timepoints we will simulate

Helper functions

Next, we need to define some functions that will help us with setting up (or initializing), running, and plotting the simulations:

def initialize_simulation():
    # zero-pad membrane potential vector 'V' and spike vector 'spikes'
    V      = np.zeros(time.size) # preallocate vector for simulated membrane potentials
    spikes = np.zeros(time.size) # vector to denote when spikes happen - spikes will be added after LIF simulation
    V[0]   = Vr # set first time point to resting potential
    return V,spikes

def logistic_map(a,x0,nsteps):
    # function to simulate logistic map:
    # x_{n+1} = a * x_n * (1-x_n)
    x    = np.zeros(nsteps)
    x[0] = x0
    for ii in range(1,nsteps):
        x[ii] = a * x[ii-1] * (1-x[ii-1])
    return x

def plot_potentials(time,V,timeSpikes):
    # plots membrane potential (V) against time (time), and marks spikes with red markers (timeSpikes)
    plt.show()
    plt.plot(time,V,'k',timeSpikes,np.ones(timeSpikes.size)*Vs,'ro')
    plt.ylabel('membrane potential (mV)')
    plt.xlabel('time (seconds)')

def check_solutions( result, solution_filename ):
    # check solutions against provided values
    solution = np.load( solution_filename )
    if ( np.linalg.norm( np.abs( result - solution ) ) < 0.1 ):
    	print( '\n\n ---- problem solved successfully ---- \n\n' )

Numerical integration

Given a differential equation of the form:

\begin{equation} x'(t) = f(x) \end{equation}

we can approximate its solution by the Euler method. To explain this, we can start with the finite difference approximation to the differential equation: \begin{equation} x'(t) \approx \frac{x(t+h) - x(t)}{h} \end{equation} which can be rearranged into the form: \begin{equation} x(t+h) \approx x(t) + hx'(t) \end{equation} and can then be numerically instantiated by the discrete recursion: \begin{equation} x_{n+1} = x_n + h f(x_n) \end{equation}

More elaborate, higher-order integration schemes then result from considering different strategies for the finite difference approximation to the original differential equation. Each integration strategy has benefits and drawbacks, which can also change the numerical result obtained.

Simulating the LIF neuron

Now, we want to write a function that integrates one time step of the LIF neuron, using the Euler method:

def integrate_and_fire( V, spikes, ii, Ie ):
    # function to integrate changes in local membrane potential and fire if threshold reached
    # V - vector of membrane potential
    # spikes - spike marker vector
    # i - index (applied to V and spikes) for current time step
    # Ie - input current at this time step (scalar of unit amp)
			
	# 1: calculate change in membrane potential (dV)
	# 2: integrate over given time step (Euler method)
	# 3: does the membrane potential exceed threshold (V > Vth)?
				
    return V,spikes # output the membrane potential vector and the {0,1} vector of spikes

Problem 1: Step current input

Simulate the LIF neuron with the given parameters from 0 to 1 second. Create a current input \(I_e\) that starts at 0 A, steps up to 15 nA at 200 ms, then steps down to 0 A at 800 ms.

def problem_1():
    #////////////////////////////////////
    #  problem 1 - step current input //
    #//////////////////////////////////
    #
    # Implement a leaky integrate and fire (LIF) neuron with parameters given 
    # above. 
    # 
    # Create a current input which:
    #       - starts at 0 A
    #       - steps up to 15 nA at stim_time[0]
    #       - steps down to 0 A at stim_time[1]
    #
    # Output:
    # Plot the resulting simulated membrane potential of the LIF neuron.
    #

    # problem-specific parameters
    stim_time = [.2,.8] # time (seconds) when current turns ON and turns OFF
    
    V,spikes = initialize_simulation() 	# initialize simulation

    ***ADD CODE HERE*** 		# iterate over each time step
    
    # add spikes to create membrane potential waveforms
    V[spikes==1] = Vs
    
    # plot membrane potential
    plot_potentials(time,V,time[spikes==1])
    plt.title('P1 : step current input')

    # output:
    V_prob1 = V
    return V_prob1

Output: Plot the results, save the membrane potential in a vector named "V_prob1", and run the checking function with problem1.npy:

check_solutions( V_prob1, 'problem1.npy' )