Minimal model of interictal and ictal discharges “Epileptor-2” (Chizhov et al 2018)

 Download zip file 
Help downloading and running models
Accession:263074
"Seizures occur in a recurrent manner with intermittent states of interictal and ictal discharges (IIDs and IDs). The transitions to and from IDs are determined by a set of processes, including synaptic interaction and ionic dynamics. Although mathematical models of separate types of epileptic discharges have been developed, modeling the transitions between states remains a challenge. A simple generic mathematical model of seizure dynamics (Epileptor) has recently been proposed by Jirsa et al. (2014); however, it is formulated in terms of abstract variables. In this paper, a minimal population-type model of IIDs and IDs is proposed that is as simple to use as the Epileptor, but the suggested model attributes physical meaning to the variables. The model is expressed in ordinary differential equations for extracellular potassium and intracellular sodium concentrations, membrane potential, and short-term synaptic depression variables. A quadratic integrate-and-fire model driven by the population input current is used to reproduce spike trains in a representative neuron. ..."
Reference:
1 . Chizhov AV, Zefirov AV, Amakhin DV, Smirnova EY, Zaitsev AV (2018) Minimal model of interictal and ictal discharges "Epileptor-2". PLoS Comput Biol 14:e1006186 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Extracellular; Realistic Network; Synapse;
Brain Region(s)/Organism:
Cell Type(s): Abstract integrate-and-fire leaky neuron;
Channel(s): I Potassium; I Sodium; Na/K pump;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: Python; Pascal/Delphi; Mathematica; Javascript;
Model Concept(s): Synaptic Plasticity; Depression; Bifurcation; Oscillations; Epilepsy; Activity Patterns; Short-term Synaptic Plasticity;
Implementer(s): Chizhov, Anton [anton.chizhov at mail.ioffe.ru];
Search NeuronDB for information about:  I Sodium; I Potassium; Na/K pump;
unit Equations2;

interface
uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, TeEngine, Series, ExtCtrls, TeeProcs, Chart, Math,
  DDSpinEdit,
  Unit1,Unit2,Stimulation;

var
    { Variables }
    t,U,w,K,Na,nu,xD,m,V,dVdt
                                                :double;
    { Parameters }
    dt,Cm,tau_w,gI_gE,tau_K,ro_pump,tau_xD,dxD_reset,tau_m,tau_Na,
    g_U,gE,gL,C,
    I_stim,I_a,k1,Gain,u0,uu,nu_max,
    UT,U_reset,dw_reset,dNa_reset,dK_reset,roK_roNa,
    U1,U2,V0, NoiseAmpl,INaKpump,
    Nai_alpha,Ko_alpha,VK,VK0,VKnorm,K_bath,K0,Na0,DB
                                                :double;
    nt,nt_end,n_Write,n_Draw
                                                :integer;
    aaa,bbb                                     :text;
    Stim                                        :TStim;


procedure Cleaning;
procedure PlotSigmoid;
procedure Integrate;
function dexp(x :double) :double;
function Sigmoid(x :double) :double;
function Sigmoid2(x :double) :double;
function Sigmoid3(x :double) :double;
function Theta(x :double) :double;
function step(a,b :double) :double;    { step = a**b }
function max(x,y:double):double;

implementation


function dexp(x :double) :double;
begin
  if      x<-20 then dexp:=0{exp(-20)}
  else if x> 20 then dexp:=exp( 20)
  else               dexp:=exp(x);
end;

function Sigmoid(x :double) :double;
begin    Sigmoid:=1/(1+dexp(-x))  end;

function Sigmoid2(x :double) :double;
begin    Sigmoid2:=2/(1+dexp(-2*x)) - 1;  end;

function Sigmoid3(x :double) :double;
begin
  if      x<0 then Sigmoid3:=0
  else if x<1 then Sigmoid3:=x
  else             Sigmoid3:=1;
end;

function Theta(x :double) :double;
begin if x>=0 then Theta:=1 else Theta:=0; end;

function step(a,b :double) :double;    { step = a**b }
begin
  step:=exp(b*ln(a));
end;

function max(x,y:double):double;
begin  if x>y then max:=x  else max:=y;  end;

procedure Parameters;
{ Basic parameters are from Izhikevich, p. 295 }
begin
  dt :=Form1.DDSpinEdit1.Value;
  Cm :=Form1.DDSpinEdit2.Value;
  g_U:=Form1.DDSpinEdit3.Value;
  n_Write:=trunc(Form1.DDSpinEdit4.Value);
  n_Draw :=trunc(Form1.DDSpinEdit19.Value);
  gE :=Form1.DDSpinEdit5.Value;
  I_a      :=Form1.DDSpinEdit7.Value;
  UT       :=Form1.DDSpinEdit8.Value;
  U_reset  :=Form1.DDSpinEdit9.Value;
  dw_reset :=Form1.DDSpinEdit10.Value;
  dK_reset :=Form1.DDSpinEdit14.Value;
  k1       :=Form1.DDSpinEdit15.Value;
  dNa_reset:=Form1.DDSpinEdit16.Value;
  tau_w  :=Form1.DDSpinEdit11.Value;
  gI_gE  :=Form1.DDSpinEdit12.Value;
  tau_K  :=Form1.DDSpinEdit13.Value;
  ro_pump:=Form1.DDSpinEdit17.Value;
  DB     :=Form1.DDSpinEdit20.Value;
  tau_xD :=Form1.DDSpinEdit21.Value;
  dxD_reset:=Form1.DDSpinEdit22.Value;
  Gain   :=Form1.DDSpinEdit23.Value;
  u0     :=Form1.DDSpinEdit24.Value;
  NoiseAmpl:=Form1.DDSpinEdit25.Value;
  tau_m  :=Form1.DDSpinEdit26.Value;
  tau_Na :=Form1.DDSpinEdit27.Value;
  roK_roNa:=Form1.DDSpinEdit28.Value;
  K0     :=Form1.DDSpinEdit29.Value;
  gL     :=Form1.DDSpinEdit30.Value;
  C      :=Form1.DDSpinEdit31.Value;
//  U1:=-5/2/g_U-sqrt(5/4/g_U/g_U-140/g_U);
//  U2:=-5/2/g_U+sqrt(25/4/g_U/g_U-140/g_U);
  U1:=-60; {mV}
  U2:=-40; {mV}
  V0:=-70; {mV}
  Nai_alpha:=20 {mM};
  Ko_alpha:=2.5 {mM};
  K_bath:=Form1.DDSpinEdit18.Value; // 8.5 {mM};
  nu_max:=100 {Hz};
end;

procedure InitialConditions;
begin
  U:=V0;
  w:=0;
  K:=K0;
  VK0   :=26.6{mV}*ln({K}K0/130{mM});
  VKnorm:=26.6{mV}*ln(3{mM}/130{mM});
  Na0:=10;
  Na:=Na0;
  m:=0;
  xD:=1;
  nu:=0;
  V:=0;      dVdt:=0;
  Stim:=TStim.Create;
  { Randomize }
  RandSeed:=2;
end;

procedure Cleaning;
begin
 { Cleaning }
  Form1.Series1.Clear;
  Form1.Series2.Clear;
  Form1.Series3.Clear;
  Form1.Series4.Clear;
  Form1.Series5.Clear;
  Form1.Series6.Clear;
  Form1.Series7.Clear;
  Form1.LineSeries1.Clear;
  Form1.LineSeries2.Clear;
  Form1.LineSeries3.Clear;
  Form1.LineSeries5.Clear;
  Form1.Series7.Active:=(tau_m>0);
  Form2.LineSeries6.Clear;
end;

procedure PlotSigmoid;
var i :integer;
begin
  assignFile(bbb,'Sigmoid.dat'); rewrite(bbb);
  Form1.LineSeries4.Clear;
  Parameters;
  for i:=0 to 150 do begin
      uu:=i;
      if Form1.CheckBox1.Checked then nu:=nu_max*max(Sigmoid3((uu-u0)/Gain),0)
                                 else nu:=nu_max*max(Sigmoid2((uu-u0)/Gain),0);
      Form1.LineSeries4.AddXY(uu,nu);
      writeln(bbb,uu:8:3,' ',nu:8:3);
  end;
  Application.ProcessMessages;
  close(bbb);
end;

procedure Writing;
begin
      { Writing }
      {           1             2         3                                    }
      write  (aaa,nt*dt:8:3,' ',U:8:3,' ',w:8:3,' ');
      {           4          5                          6                      }
      write  (aaa,U2:8:3,' ',INaKpump*1e3{mM/s}:8:3,' ',uu:8:3,' ');
      {           7           8          9                                     }
      write  (aaa,nu :8:3,' ',V :8:3,' ',K :8:3,' ');
      {           10          11          12                                   }
      write  (aaa,Na :8:3,' ',xD :8:3,' ',Stim.Current :8:3,' ');
      {           13                 14                                          }
      writeln(aaa,nt*dt/1000:8:3,' ',Stim.Rate :8:3);
end;

procedure Integrate;
var
    dUdt,dwdt,dKdt,dNadt,dxDdt,dmdt,I_ex,I_in,y
                                              :double;
begin
  Parameters;
  InitialConditions;
  Cleaning;
  PlotSigmoid;
  assignFile(aaa,'t_U_w_U2_pump_uu_nu_m_K_Na_xD.dat'); rewrite(aaa);
  nt:=0;
  repeat nt:=nt+1;  t:=nt*dt;
    Parameters;
    VK:=26.6{mV}*ln(K/130{mM});
    { Na-K pump, [Wei 2014] }
    INaKpump:=ro_pump/(1+dexp(3.5-K))/(1+dexp((25-Na)/3));    {M/s}
    { Stimulation }
    Stim.Stimulate(t,dt);
    I_stim  :=(VK-VKnorm)*k1 + NoiseAmpl/sqrt(dt)*randG(0,1) + Stim.Current; // sqrt is introduced on 26.09.2017
    I_ex:=       gE*m*xD;
    I_in:=-gI_gE*gE*m;  if I_stim+I_ex+I_in>DB then I_in:=0;
    { Total current }
    uu      :=I_stim + I_ex + I_in{gE*m*(xD-0.5)} - INaKpump;     // !!!  Na-K-pump is new since 22.05.2018
    { Izhikevich *************************************}
    dUdt :=1/Cm    *(g_U*(U-U1)*(U-U2) - w + uu + I_a);
    dwdt :=1/tau_w *(0-w);
    { Sigmoid ****************************************************************}
    y:=uu;
    if Form1.CheckBox3.Checked { with V-equation } then y:=V;
    if Form1.CheckBox1.Checked then  nu:=nu_max*    Sigmoid3((y-u0)/Gain)
                               else  nu:=nu_max*max(Sigmoid2((y-u0)/Gain),0);
    { Main eqs. **************************************************************}
    dKdt :=1/tau_K *(K_bath-K) - 2*roK_roNa*INaKpump + dK_reset *(nu{+Stim.Rate})/1000;
    dNadt:=1/tau_Na*(Na0-Na)   - 3*         INaKpump + dNa_reset*(nu+Stim.Rate)/1000;
    if tau_m=0 then m:=nu else
    dmdt :=1/tau_m *(nu-m);
    dxDdt:=1/tau_xD*(1-xD)                  - dxD_reset*nu*xD/1000;
    dVdt :=1/C*(-gL*V + uu);
    {*************************************************************************}
    U :=U +dt*dUdt;
    w :=w +dt*dwdt;
    K :=K +dt*dKdt;     K:=max(0,K);
    Na:=Na+dt*dNadt;
    m :=m +dt*dmdt;
    xD:=xD+dt*dxDdt;
    V :=V +dt*dVdt;
    {****************************************************************}
    { Drawing }
    if (U>=UT) then begin
       n_Write:=1;  n_Draw:=1;
    end;
    if nt mod n_Draw = 0 then begin
       Form1.Series1.AddXY(nt*dt,U);
       Form1.Series2.AddXY(nt*dt,w);
       Form1.Series3.AddXY(nt*dt,U2);
       Form1.Series4.AddXY(nt*dt,INaKpump);
       Form1.Series5.AddXY(nt*dt,uu);
       Form1.Series6.AddXY(nt*dt,nu);
       if tau_m>0 then
       Form1.Series7.AddXY(nt*dt,m);
       Form1.LineSeries1.AddXY(nt*dt,K);
       Form1.LineSeries2.AddXY(nt*dt,Na);
       Form1.LineSeries3.AddXY(nt*dt,xD);
       Form1.LineSeries5.AddXY(nt*dt,V);
       Form2.LineSeries6.AddXY(nt*dt,max(Stim.Current,Stim.Rate));
    end;
    if nt mod 1000*n_Draw = 0 then  Application.ProcessMessages;
    if (nt mod n_Write = 0) or (Stim.Status='Pulse') then  Writing;
    {****************************************************************}
    { Reset }
    if (U>=UT) then begin
       U :=U_reset;
       w :=w  + dw_reset;
    end;
    nt_end:=trunc(Form1.DDSpinEdit6.Value/dt);
  until nt>=nt_end;
  Application.ProcessMessages;
  closeFile(aaa);
end;

end.

Loading data, please wait...