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.
