Draft:Epileptor-2

Epileptor-2 is the simplified mathematical model of epileptic discharges observed in the brain, named after the previous model called Epileptor. It replicates brief interictal discharges—observed as clusters of action potential spikes in the activity of individual neurons—and longer ictal discharges, represented as clusters of these shorter discharges. According to the Epileptor-2 model, brief interictal discharges are characterized as stochastic oscillations of the membrane potential and synaptic resources, while ictal discharges emerge as oscillations in the extracellular concentration of potassium ions and the intracellular concentration of sodium ions. The both models, Epileptor and Epileptor-2, demonstrate that ionic dynamics play a decisive role in the generation of pathological activity.

Below, a modified Epileptor-2 model is described, where the second extracellular compartment has been added. Ionic dynamics is described with two extracellular concentrations of potassium ions, at the proximity of neurons, $[K]_{o1}$, and at a distance, $[K]_{o2}$, and the intracellular concentration of sodium ions, $[Na]_{i}$:

${d[K]_{o1} \over dt}={[K]_{o2}-[K]_{o1} \over \tau _{K1}}+\delta _{K}~\nu (t)-2\beta I_{pump}(t)$

${d[K]_{o2} \over dt}={[K]_{bath}-[K]_{o2} \over \tau _{K2}}+{[K]_{o1}-[K]_{o2} \over \tau _{K1}}$

${d[Na]_{i} \over dt}={{[Na]_{i}}^{0}-[Na]_{i} \over \tau _{Na}}+\delta _{Na}\nu (t)-3I_{pump}(t)$

According to these equations, the concentrations increase because of firing with the rate $\nu$, relax with the time scales $\tau _{K1},~\tau _{K2}$, and $\tau _{Na}$, and decrease due to the sodium-potassium pump $I_{pump}$.

The membrane polarization $V$ averaged across population and the synaptic resource $x$ are governed by

${dV \over dt}={u\left(t\right)-V \over \tau _{m}}$

${dx \over dt}={1-x \over \tau _{x}}-\delta _{x}x~\nu (t)$

where

$u(t)=g_{K}(V_{K}(t)-V_{K}^{0})+G_{syn}(x-0.5)~\nu (t)+\sigma ~\eta (t)$

The source term $u(t)$ summarises the current mediated by the changes ot the potassium concentration, the synaptic current and the noise $\eta (t)$ (white gaussian noise). The synaptic resource decreases because of firing with the rate $\nu$ and recovers with the time scale $\tau _{x}$.

The potassium reversal potential is

$V_{K}(t)=26.6~\ln([K]_{o2}/130)$

The sodium-potassium pump is

$I_{pump}(t)=\rho /\left(1+\exp(3.5-[K]_{o1})/(1+\exp((25-[Na]_{o1})/3)\right)$

The sigmoidal dependence of the firing rate on voltage is

$\nu (t)=100\left(2/\left(1+\exp((6-V)/10)\right)-1\right)$

Representative neuron is governed by the quadratic integrate-and-fire model for the membrane potential $U$. It receives the same input $u(t)$ as the population:

${dU \over dt}={g_{U} \over C_{U}}(U-U_{1})(U-U_{2})+{g_{L} \over C_{U}}u(t),~~~~~if~~~U>U^{th}~~then~~U=U_{reset}$

The parameters are as follows:

$\tau _{K1}=25s,~\tau _{K2}=250s,~\tau _{Na}=20s,~\tau _{x}=2s,~\tau _{m}=0.01s,$$~\delta _{K}=0.04~mM,~\delta _{Na}=0.03~mM,~\delta _{x}=0.01,~\rho =0.8,~\beta =10,~\sigma =8,$$~G_{syn}/g_{L}=2.5~mV\cdot s,~g_{K}/g_{L}=0.5,~g_{L}=5~nS,$$[K]_{o}^{0}=3~mM,~[K]_{bath}=8.5~mM,~[Na]_{i}^{0}=10~mM,$$g_{U}=0.4,~C_{U}=0.2,~U^{th}=25,~U_{reset}=-50,~U_{1}=-60,~U_{2}=-40,~U_{0}=-70$

Ictal discharges in experiment and Epileptor-2 model. Membrane potential and extracellular potassium concentration have been recorded in rat brain slice during induced epileptiform activity. Membrane potential and the extracellular potassium concentration obtained in simulations are compared to data. The ictal discharges are composed of short bursts each containing several action potentilals.

The model's implementation in Python is given below:

import numpy as np
import pylab as pl


def Sigmoid(V):
    scc = 2. / (1 + np.exp(2 * (Vth - V) / K_v)) - 1
    psc = np.clip(scc, a_min=0, a_max=None)
    return vmax * psc

def I_pump(K, Na):
    return rho / ((1 + np.exp(3.5 - K)) * (1 + np.exp((25 - Na) / 3)))

def VK(Kconc):
    return 26.6 * np.log(Kconc / 130.)

def u_over_gL(K, nu, x, y, Ipi):
    eta = np.random.normal(0, 1)
    vkk = VK(y)
    #********************************************************************************************
    return (g_K * (vkk - vkk0) + G_syn * nu * (x - 0.5)
            + sigma * eta / np.sqrt(dt) / np.sqrt(1000))
    #********************************************************************************************

def dots(states):
    K, Na, x, V, U, y = states
    Ip_ = I_pump(K, Na)
    nu_ = Sigmoid(V)
    uu_ = u_over_gL(K, nu_, x, y, Ip_)
    #*************************************************************
    dotK  = (y    - K ) / tau_K  + deltaK  * nu_ - 2 * beta * Ip_
    dotNa = (Na_0 - Na) / tau_Na + deltaNa * nu_ - 3 * Ip_
    dotV  = (uu_ - V) / tau_m
    dotx  = (1 - x) / tau_x - deltax * x * nu_
    doty  = (Kbath - y) / tau_y + (K-y)/tau_K
    dotU  = 1000 * ((gU / CU) * (U - U1) * (U - U2) + gL * uu_ / CU)
    #*************************************************************
    return np.array([dotK, dotNa, dotx, dotV, dotU, doty])

# in ms
tau_K = 100/4  # s
tau_Na = 20  # s
tau_m = 0.01  # s
tau_x = 2  # s
tau_y = 250

deltaK = 0.02*2         # switches from ID to IID
deltaNa = 0.03
deltax = 0.01
rho = 0.2*4             # mM/ms
beta = 10.
sigma = 25./3
G_syn = 5./2
g_K = 0.5
K_0 = 3.
Kbath1 = 3
Kbath2 = 8.5
Na_0 = 10.
vmax = 100
Vth = 25./4
K_v = 20.
gL = 1*5  # from their mathematica notebook
# parameters of the observer:
gU = 0.4
CU = 200.
Uth = 25.
Vr = -50.
U1 = -60.
U2 = -40.
U_0 = -70.

V_0 = 0
x_0 = 1.
y_0 = Kbath1
nu_0 = 0
Ip_0 = 0
vkk0 = VK(K_0)

dt = 0.01  # in second
u_0 = u_over_gL(K_0, nu_0, x_0, y_0, Ip_0)

tt = 400.  # in second
lngth = int(tt / dt)
xax = np.arange(0, tt, dt)

arbig = np.zeros((lngth, 9))
arbig[0] = [K_0, Na_0, x_0, V_0, U_0, y_0, nu_0, Ip_0, u_0]
Ki, Nai, xi, Vi, Ui, yi, nui, Ipi, ui = arbig[0]

for ii, it in enumerate(xax[1:], 1):
    if it <= 50:
        Kbath = Kbath1
    else:
        Kbath = Kbath2
    state = dots(arbig[ii - 1, :6])
    Ki, Nai, xi, Vi, Ui, yi = arbig[ii - 1, :6] + dt * state

    #***********************************
    Ipi = I_pump(Ki, Nai)
    nui = Sigmoid(Vi)
    ui = u_over_gL(Ki, nui, xi, yi, Ipi)
    #***********************************

    if Ui > Uth:
        Ui = Vr
    arbig[ii] = [Ki, Nai, xi, Vi, Ui, yi, nui, Ipi, ui]

# Plotting **********************************************************************************
KK, NNa, xx, VV, UU, yy, nuu, II, uu = arbig.T

pl.close('all')
pl.figure(figsize=(10, 7))

ax = pl.subplot(511)
pl.plot(xax, UU, 'k')
pl.ylabel(r'$U (mV)$')

pl.subplot(512, sharex=ax)
pl.plot(xax, VV, 'r')
pl.ylabel(r'$V(mV)$')

pl.subplot(513, sharex=ax)
pl.plot(xax, xx, color=(0.5, 0., 0.5))
pl.ylabel(r'$x^D$')

pl.subplot(514, sharex=ax)
pl.plot(xax, NNa, 'r')
pl.ylabel(r'$[Na]_i$ (mM)')

pl.subplot(515, sharex=ax)
pl.plot(xax, (yy+KK)/2, color=(0.2, 0.6, 0.2), linewidth=3)
pl.ylabel(r'$[K]_o$ (mM)')

pl.tight_layout()
pl.show()