/*--------------------------------------------------------------------------
Author: Thomas Nowotny
Institute: Institute for Nonlinear Dynamics
University of California San Diego
La Jolla, CA 92093-0402
email to: tnowotny@ucsd.edu
initial version: 2002-01-25
--------------------------------------------------------------------------*/
#ifndef LPRNEURON_CC
#define LPRNEURON_CC
#include "CN_neuron.cc"
#define efunc(X,Y,Z) (1.0/(1.0+exp(((X)-(Y))/(Z))))
#define Eaxon x[idx]
#define Esoma x[idx+15]
LPRneuron::LPRneuron(int inlabel, double *the_p= LPR_p):
neuron(inlabel, LPR_IVARNO, LPRNEURON, the_p, LPR_PNO)
{
}
LPRneuron::LPRneuron(int inlabel, vector<int> inpos, double *the_p= LPR_p):
neuron(inlabel, LPR_IVARNO, LPRNEURON, inpos, the_p, LPR_PNO)
{
}
inline double LPRneuron::E(double *x)
{
using namespace LPR5;
return x[idx+15]+p[V_shift];
}
void LPRneuron::currents(ostream &os, double *x)
{
using namespace LPR5;
static double INa, ICa, IoCa, Id, IA, Ih, Il, IM, IVV, tmp;
// differential eqn for the axon membrane potential
INa= pw3(x[idx+1])*x[idx+2]*p[g_Na]*(Eaxon-p[V_Na]);
Id= pw4(x[idx+8])*p[g_Kd]*(Esoma-p[V_K]);
IM= p[g_M]*x[idx+14]*(Eaxon-p[V_M]);
Il= p[g_leak]*(Eaxon-p[V_leak]);
IVV= p[g_VV]*(Esoma-Eaxon);
os << INa << " ";
os << Id << " ";
os << IM << " ";
os << Il << " ";
os << IVV << " ";
tmp= exp(Esoma/p[RTF]);
ICa= (x[idx+3]*x[idx+4]*p[g_CaT]+x[idx+5]*p[g_CaS])*
p[P_Ca]*Esoma*(x[idx+13]*tmp-p[Ca_out])/(tmp-1.0);
IoCa= x[idx+6]*x[idx+7]*p[g_KCa]*(Esoma-p[V_K]);
tmp= efunc(Esoma,p[V_aA],p[s_aA]);
IA= pw3(x[idx+9])*p[g_A]*(Esoma-p[V_K])*(tmp*x[idx+10]+(1.0-tmp)*x[idx+11]);
Ih= x[idx+12]*p[g_h]*(Esoma-p[V_h]);
Il= p[g_leaks]*(Esoma-p[V_leak]);
os << ICa << " ";
os << IoCa << " ";
os << IA << " ";
os << Ih << " ";
os << Il << endl;
}
void LPRneuron::derivative(double *x, double *dx)
{
using namespace LPR5;
static double a, b, hinf, kh, minf, km;
static double INa, Id, ICa, IoCa, IA, Ih, Il, IM, IVV;
static double Isyn;
Isyn= 0.0;
forall(den, den_it) {
Isyn+= (*den_it)->Isyn(x);
}
// differential eqn for the axon membrane potential
INa= pw3(x[idx+1])*x[idx+2]*p[g_Na]*(Eaxon-p[V_Na]);
Id= pw4(x[idx+8])*p[g_Kd]*(Eaxon-p[V_K]);
IM= p[g_M]*x[idx+14]*(Eaxon-p[V_M]);
Il= p[g_leak]*(Eaxon-p[V_leak]);
IVV= p[g_VV]*(Esoma-Eaxon);
dx[idx]= (-INa - Id - IM - Il + IVV)/p[C_a];
// soma compartment
a= exp(Esoma/p[RTF]);
ICa= (x[idx+3]*x[idx+4]*p[g_CaT]+x[idx+5]*p[g_CaS])
*p[P_Ca]*Esoma*(x[idx+13]*a-p[Ca_out])/(a-1.0);
IoCa= x[idx+6]*x[idx+7]*p[g_KCa]*(Esoma-p[V_K]);
a= efunc(Esoma,p[V_aA],p[s_aA]);
IA= pw3(x[idx+9])*p[g_A]*(Esoma-p[V_K])*(a*x[idx+10]+(1.0-a)*x[idx+11]);
Ih= x[idx+12]*p[g_h]*(Esoma-p[V_h]);
Il= p[g_leaks]*(Esoma-p[V_leak]);
dx[idx+15]= (-ICa - IoCa - IA - Ih - Il - IVV +
(p[I_DC]+Isyn)*p[I_scale])/p[C_s];
// differential eqn for mNa
a= (3.5+0.1*x[idx]) / (1.0-exp(-3.5-0.1*x[idx]));
b= 4.0*exp(-(x[idx]+60.0)/18.0);
dx[idx+1]= (a*(1.0-x[idx+1])-b*x[idx+1])*p[k_slow]*p[k_fast];
// differential eqn for hNa
a= 0.07*exp(-x[idx]/20.0-3.0);
b= 1.0 / (exp(-3.0-0.1*x[idx])+1.0);
dx[idx+2]= (a*(1.0-x[idx+2])-b*x[idx+2])*p[k_slow]*p[k_fast];
// differential eqn for md
a= (-0.5-0.01*x[idx]) / (exp(-5.0-0.1*x[idx])-1.0);
b= 0.125*exp(-(x[idx]+60.0)/80.0);
dx[idx+8]= (a*(1.0-x[idx+8])-b*x[idx+8])*p[k_slow]*p[k_fast];
// differential eqn for mCaT
minf= efunc(Esoma,p[V_mCaT],p[s_mCaT]);
km= p[k_mCaT]*p[k_slow];
dx[idx+3]= (minf-x[idx+3])*km;
// differential eqn for hCaT
hinf= efunc(Esoma,p[V_hCaT],p[s_hCaT]);
kh= p[k_hCaT]*efunc(Esoma,p[V_khCaT],p[s_khCaT])*p[k_slow];
dx[idx+4]= (hinf-x[idx+4])*kh;
// differential eqn for mCaS
minf= efunc(Esoma,p[V_mCaS],p[s_mCaS]);
km= p[k_mCaS]*p[k_slow];
dx[idx+5]= (minf-x[idx+5])*km;
// differential eqn for mKCa
a= p[V_mKCa1]-p[f]*x[idx+13];
b= p[V_mKca2]-p[f]*x[idx+13];
minf= efunc(Esoma,a,p[s_mKCa1])
*efunc(Esoma,b,p[s_mKca2])*(x[idx+13]/(p[c_mKCa]+x[idx+13]));
km= p[k_mKCa]*p[k_slow];
dx[idx+6]= (minf-x[idx+6])*km;
// differential eqn for hKCa
hinf= p[c_hKCa1]/(p[c_hKCa2]+x[idx+13]);
kh= p[k_hKCa]*p[k_slow];
dx[idx+7]= (hinf-x[idx+7])*kh;
// differential eqn for mA
minf= efunc(Esoma,p[V_mA],p[s_mA]);
km= p[k_mA]*p[k_slow];
dx[idx+9]= (minf-x[idx+9])*km;
// differential eqn for hA1
hinf= efunc(Esoma,p[V_hA],p[s_hA]);
kh= p[k_hA1]*p[k_slow];
dx[idx+10]= (hinf-x[idx+10])*kh;
// differential eqn for hA2
kh= p[k_hA2]*efunc(Esoma,p[V_khA2],p[s_khA2])*p[k_slow];
dx[idx+11]= (hinf-x[idx+11])*kh;
// differential eqn for mh
minf= efunc(Esoma,p[V_mh],p[s_mh]);
km= p[k_mh]/efunc(Esoma,p[V_kmh],p[s_kmh])*p[k_slow];
dx[idx+12]= (minf-x[idx+12])*km;
// differential eqn for Ca concentration
dx[idx+13]= -p[c_ICa]*ICa - p[k_Ca]*(x[idx+13]-p[Ca_0]);
// differential equn for M current activation var
minf= efunc(Eaxon,p[V_mM],p[s_mM]);
km= p[k_mM]*efunc(Eaxon,p[V_kmM],p[s_kmM])*p[k_slow];
dx[idx+14]= (minf-x[idx+14])*km;
}
#undef efunc
#undef Eaxon
#undef Esoma
#endif
