AdEx Model in Python from Essentials of Neuroscience with MATLAB: Module 3

adex-model-in-python.py

# This is a draft of a Python version of the code of this video: https://www.youtube.com/watch?v=8lDF0acgRdI&list=PLn0OLiymPak1b2aYULx6hDVU7wSGEUJqw&index=17
# Disclaimer: The output looks similar-enough to me. Certainly there are off-by-one differences and there might be other differences too
import math
import numpy as np
import matplotlib.pyplot as plt

# Part 1. Setup

SIMULATION_DURATION = 500  # Simulation time in ms
dt = 0.01  # descretization in ms

IMax = 500  # pA


def _t(t):
    return int(t / dt)


timevec = np.arange(0, SIMULATION_DURATION, dt)

Iinput = np.zeros(_t(SIMULATION_DURATION))
Iinput[_t(100) : _t(400)] = IMax

# Part 2. Constants

C = 200  # Capacitance
gl = 10  # Leak conductance
El = -58  # Leak voltage
vt = -50  # Resting membrane potential
delT = 2  # Spike width factor
a = 2  # Resonator/integrator variable
tauW = 120  # Adaptation decay time
b = 100  # Adaptation jump (update for w)
vreset = -46  # Voltage reset

# Part 3. Simulation

membPoten = np.zeros(_t(SIMULATION_DURATION))
vpeak = 0  # spike threshold

for timei in range(_t(SIMULATION_DURATION)):
    if timei == 0:
        v = vt
        w = 0
    else:
        Iin = Iinput[timei]
        v += dt * (-gl * (v - El) + gl * delT * math.exp((v - vt) / delT) + Iin - w) / C
        w += dt * (a * (membPoten[timei - 1] - El) - w) / tauW

    if v >= vpeak:
        v = vreset
        w = w + b
    else:
        membPoten[timei] = v

# Part 4. Plot

fig, ax1 = plt.subplots()

l2 = ax1.plot(timevec, membPoten, color="b", label="membPoten")
ax1.set_ylabel("Membrane Voltage (mV)")

ax2 = ax1.twinx()
l1 = ax2.plot(timevec, Iinput, color="g", label="Iinput")
plt.ylim(0, 5000)
ax2.set_ylabel("Iinput Current (pA)")


plt.xlabel("Time (ms)")
plt.title("Iinput vs Time and MembPoten vs Time")

ax1.legend(handles=l1 + l2)