Leaky Integrate and Fire Neuron

In this notebook we will use Python to simulate the leaky integrate and fire (LIF) neuron model. We’ll investigate, in particular, how the spiking activity varies as we adjust the input current \(I\).

1 Preliminaries

Before beginning, let’s load in the Python package we’ll need:

2 Numerical solutions - Introduction

How do we compute a numerical solution to the integrate and fire model? The basic idea is to rearrange the differential equation to get \(V(t+1)\) on the left hand side, and \(V(t)\) on the right hand side. Then, if we know what’s happening at time \(t\), we can solve for what’s happening at time \(t+1\).

For example, consider the differential equation:

\[ \dfrac{dV}{dt} = -\dfrac{V-V_*}{\tau} \]

In words, we can think of:

\(dV\) as the “change in voltage V”,

\(dt\) as the “change in time t”.

Let’s consider the case that we record the voltage \(V\) in discrete time steps. So we observe:

\(V[0], V[1], V[2], \ldots\)

at times:

\(dt, \, 2*dt, \, 3*dt, \ldots\)

where \(dt\) is the time between our samples of \(V\).

We can now write the “change in voltage V” as:

\[ dV = V(t+1) - V(t) \]

Notice that the change in voltage is the difference in V between two sequential time samples. Now, let’s rewrite \(\dfrac{dV}{dt}\) as,

\[ \dfrac{dV}{dt} = \dfrac{ V(t+1) - V(t) }{ dt } \]

where we’ve replaced \(dV\). Now, let’s substitute this expression into the equation at the top of this file:

\[ \dfrac{ V(t+1) - V(t) }{ dt } = \dfrac{I}{C}. \]

Solving this equation for \(V(t+1)\) you’ll find that:

\[ V(t+1) = V(t) + dt*(-\dfrac{V(t)-V_*}{\tau}) \]

Notice that, in this expression, we use our current value of the voltage \(V(t)\) and the model (\(-\dfrac{V(t)-V_*}{\tau}\)) to determine the next value of the voltage V(t+1).

Now, let’s program this equation in Python. First, let’s set the values for the parameters \(V_*\) and \(\tau\).

We also need to set the value for \(dt\). This defines the time step for our model. We must choose it small enough so that we don’t miss anything interesting. We’ll choose:

Let’s assume the units of time are seconds. So, we step forward in time by \(0.01\) s.

The right hand side of our equation is nearly defined, but we’re still missing one thing, \(V(t)\).

Q: What value do we assign to \(V(t)\)?

So here’s an easier question: what initial value do we assign to \(V(t)\)?

To start, we’ll create an array of zeros to hold our results for \(V\):

This array V consists of 1000 rows and 1 column. We can think of each row entry as corresponding to a discrete step in time. Our goal is to fill-in the values of V (i.e., step forward in time), in a way consistent with our model.

Let’s choose an initial value for V of 0.2, which in our simple model we’ll assume represents the rest state.

Q: Given the initial state V[0]=0.2, calculate V[1]. Then calcualte V[2].

After the two calculations above, we’ve moved forward two time steps into the future, from \(t=0\) s to \(t=0.01\) s, and then from \(t=0.01\) s to \(t=0.02\) s. But what if we want to know \(V\) at \(t=10\) s? Then, this iteration-by-hand procedure becomes much too boring and error-prone. So, what do we do? Let’s make the computer do it …

3 Numerical solutions - Implementation

Let’s computerize this iteration-by-hand procedure to find V[999].

To do so, we’ll use a for-loop.

Here’s what it looks like:

Q: Does this loop make sense? Describe what’s happening here. What is something?

Q: Why does the range command end at 999?

Execute this for-loop and examine the results in vector V. To do so, let’s plot V:

Q: What happens to the voltage after 1000 steps?

This plot is informative, but not great. Really, we’d like to plot the voltage as a function of time, not steps or indices. To do so, we need to define a time axis:

Q: What’s happening in the command above? Does it make sense? (If not, trying printing or plotting t.)

Now, with time defined, let’s redo the plot of the voltage with the axes labeled appropriately.

Finally, let’s put it all together . . .

4 LIF CODE (version 1)

So far, we constructed parts of the I&F model in bits-and-pieces. Let’s now collect all the code, compute a numerical solution to the I&F model, and plot the results (with appropriate axes).

Q: Adjust the parameter Vstar. What happens to V if Vstar=0? Can you set Vstar so that V > 20 within 10 s?

5 Voltage threshold & reset

Notice, our model is missing something important: the reset.

Q: Without the reset, how does the voltage behave as \(t\rightarrow\infty\) (if \(Vstar>0\))?

Now, let’s update our model to include the reset. To do so, we’ll add two things to our code.

  • First, we’ll define the voltage threshold Vth, and reset voltage Vreset.
  • Second, we’ll check to see if V exceeds Vth using an if-statement; if it does, then we’ll set V equal to Vreset.

Q: How will you update the code below to include the reset?

6 I&F CODE (version 2)

Now, let’s put it all together to make a complete I&F model (with a thershold and reset), simulate it, and plot the result.

Q: Adjust the parameter I. What happens to V if I=10? If I=100?

Q: Adjust the parameter C. What happens to V if C=0.1? If C=10?

Q: What is “spiking” in this I&F model?

7 Discussion

  1. Describe the main components of the IF model. Draw it (in some way).

  2. Describe the main components of the LIF model. Draw it (in some way).

  3. Describe the differences and similarities between the IF and LIF models.

  4. The IF model is meant to mimic a neurons activity. What is realistic about the IF model? What is unrealistic?

  5. Describe the roles of the IF model parameters Vreset and Vthreshold.

  6. Consider the IF model. Sketch voltage (V) versus time (t) for a small input current, for a large input current.

    1. How does an increase in capacitance (C) impact the dynamics?
    2. Can you interpret this physically?

  7. Plot the f-I curve for the IF model.