0

I have a questions with regards to the implementation of the leaky integrate and fire model in Python using ODE solvers.

I have an implementation using the Euler method with a for loop:

import numpy as np
import matplotlib.pyplot as plt

# Define parameters
tau_m = 10  # membrane time constant (ms)
v_rest = -70  # resting membrane potential (mV)
v_thresh = -55  # spike threshold (mV)
v_reset = -80  # reset potential (mV)
i_e = 3  # input current (nA)
dt = 0.1  # time step (ms)
t_max = 100  # simulation time (ms)

# Initialize variables
v = v_rest  # membrane potential (mV)


# Simulate LIF model
for t in np.arange(0, t_max, dt):
    dv_dt = (-(v - v_rest) + i_e) / tau_m  # membrane potential differential equation
    v += dv_dt * dt  # update membrane potential
    if v >= v_thresh:  # spike condition
        v = v_reset  # reset membrane potential

Now, the issue is that I have a much larger model implemented using solve_ivp() and I need to integrate the LIF model into it.

Defining the differential term of the membrane potential (dv_dt) is doable and I can do it as follows:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Define LIF model
def lif(t, y, tau_m, v_rest, v_thresh, v_reset, i_e):
    dv_dt = (-(y - v_rest) + i_e) / tau_m  # membrane potential differential equation
    return dv_dt

# Define parameters
tau_m = 10  # membrane time constant (ms)
v_rest = -70  # resting membrane potential (mV)
v_thresh = -55  # spike threshold (mV)
v_reset = -80  # reset potential (mV)
i_e = 3  # input current (nA)
t_span = (0, 100)  # simulation time span (ms)
y0 = v_rest  # initial membrane potential (mV)

# Simulate LIF model
sol = solve_ivp(lambda t, y: lif(t, y, tau_m, v_rest, v_thresh, v_reset, i_e), t_span, [y0])

# Extract results
t = sol.t
v = sol.y[0]

But I was not sure how to implement the hard change in membrane potential after depolarization.

v = v_reset if v >= v_thresh  # spike condition

since when you call the solver solve_ivp(), the only output of the function is the differential term dv_dt. I have thought about doing:


dv_dt = (v_reset - y) if y >= v_thresh else (-(y - v_rest) + i_e) / tau_m

So that I implement the change to y in the next iteration through dv_dt, but although this decreases the membrane potential, it does not reset it to v_reset.

Any ideas would be appreciated!

OMS
  • 57
  • 1
  • 9
  • 1
    There is no magic trick to do this triggered state change, you need to break the integration with a terminal event if the threshold condition is triggered, manually set the new state and restart the integration. See a discussion of the classical bouncing ball: https://stackoverflow.com/questions/70369268/implement-ball-bouncing-in-free-fall-using-scipy-ivp – Lutz Lehmann Apr 21 '23 at 12:20

1 Answers1

0

After exploring all possible options, I came up with two solutions that go around having to stop and restart the integration. I believe this is inefficient for every cell depolarization, especially if modelling a neural network.

The three other solutions

  1. Use the Euler method (As stated in the original question)

  2. Modify the integrate and fire model to depolarize and repolarize based on the membrane potential, but applying fast changes to dVm/dt, and not Vm directly:


def depo(Vm, Vth):
  s = 1000
  return s if Vm >= Vth else 0

def repo(Vm, Vmax):
  s = 10000
  return s if Vm >= Vmax else 0


# Define LIF model
def lif(t, y, tau_m, v_rest, v_thresh, v_reset, i_e):
    dv_dt = (-(y - v_rest) + i_e + depo(y[0], Vth) - repo(y[0], Vmax)) / tau_m  # membrane potential differential equation
    return dv_dt
  1. Use the Morris Lecar or any other simplified model that doesn't require direct changes to the solution of the differential equation.

I hope this helps someone if they have the same question I had.

OMS
  • 57
  • 1
  • 9