/*-------------------------------------------------------------------------- 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: 2005-08-19 --------------------------------------------------------------------------*/ #ifndef CN_LPANEURON_CC #define CN_LPANEURON_CC #include "CN_neuron.cc" LPAneuron::LPAneuron(int inlabel, double *the_p= LPA_p): neuron(inlabel, LPA_IVARNO, LPANEURON, the_p, LPA_PNO) { } LPAneuron::LPAneuron(int inlabel, tnvector inpos, double *the_p= LPA_p): neuron(inlabel, LPA_IVARNO, LPANEURON, inpos, the_p, LPA_PNO) { } inline double LPAneuron::E(double *x) { return x[idx]; } #define efunc(X,Y,Z) (1.0/(1.0+expl(((X)+(Y))/(Z)))) #define E x[idx] #define Capacit p[0] #define tau_Ca p[1] #define f p[2] #define C_Ca_0 p[3] #define E_Na p[4] #define E_K p[5] #define E_H p[6] #define E_leak p[7] #define g_H p[8] #define g_leak p[9] #define g_Na p[10] #define g_CaT p[11] #define g_CaS p[12] #define g_A p[13] #define g_KCa p[14] #define g_Kd p[15] #define Area p[16] void LPAneuron::currents(ostream &os, double *x) { _E_Ca=(12.5*log(3000.0/x[idx+12])); I_Na= g_Na*ipower(x[idx+1],3)*x[idx+2]*(E-E_Na)*Area; I_CaT= g_CaT*ipower(x[idx+3],3)*x[idx+4]*(E-_E_Ca)*Area; I_CaS= g_CaS*ipower(x[idx+5],3)*x[idx+6]*(E-_E_Ca)*Area; I_A= g_A*ipower(x[idx+7],3)*x[idx+8]*(E-E_K)*Area; I_KCa= g_KCa*ipower(x[idx+9],4)*(E-E_K)*Area; I_Kd= g_Kd*ipower(x[idx+10],4)*(E-E_K)*Area; I_H= g_H*x[idx+11]*(E-E_H)*Area; I_leak= g_leak*(E-E_leak)*Area; os << I_Na << " "; os << I_CaT << " "; os << I_CaS << " "; os << I_A << " "; os << I_KCa << " "; os << I_Kd << " "; os << I_H << " "; os << I_leak << " "; } void LPAneuron::derivative(double *x, double *dx) { Isyn= 0.0; forall (den_it) { Isyn+= den_it->c_value()->Isyn(x); } _E_Ca=(12.5*log(3000.0/x[idx+12])); I_Na= g_Na*ipower(x[idx+1],3)*x[idx+2]*(E-E_Na); I_CaT= g_CaT*ipower(x[idx+3],3)*x[idx+4]*(E-_E_Ca); I_CaS= g_CaS*ipower(x[idx+5],3)*x[idx+6]*(E-_E_Ca); I_A= g_A*ipower(x[idx+7],3)*x[idx+8]*(E-E_K); I_KCa= g_KCa*ipower(x[idx+9],4)*(E-E_K); I_Kd= g_Kd*ipower(x[idx+10],4)*(E-E_K); I_H= g_H*x[idx+11]*(E-E_H); I_leak= g_leak*(E-E_leak); dx[idx]= Isyn/Capacit -(I_Na + I_Kd + I_CaT + I_CaS + I_A + I_KCa + I_H + I_leak)*Area/Capacit; // differential eqn for mNa _minf=efunc(E,25.5,-5.29); _taum=2.64-2.52*efunc(E,120.0,-25.0); dx[idx+1]=(_minf-x[idx+1])/_taum; // differential eqn for hNa _hinf= efunc(E,48.9,5.18); _tauh= (1.34*efunc(E,62.9,-10.0))*(1.5+efunc(E,34.9,3.6)); if (_tauh < 1e-2) _tauh= 1e-2; // this is a horrible numerical instability // if the neuron is hyperpolarized, _tauh goes to 0 amplifying the already // clear disadvantageous rounding errors in (order of 1 - order of 1) // in mNa *a lot*. This leads to an explosion and eventually nan ... // cerr << _hinf << " " << _tauh << " "; dx[idx+2]=(_hinf-x[idx+2])/_tauh; // differential eqn for mCa1 _minf= efunc(E,27.1,-7.2); _taum= 43.4-42.6*efunc(E,68.1,-20.5); dx[idx+3]=(_minf-x[idx+3])/_taum; // differential eqn for hCa1 _hinf= efunc(E,32.1,5.5); _tauh= 210.0-179.6*efunc(E,55.0,-16.9); dx[idx+4]=(_hinf-x[idx+4])/_tauh; // differential eqn for mCa2 _minf= efunc(E,33.0,-8.1); _taum= 2.8+(14.0/(exp((E+27.0)/10.0)+exp((E+70.0)/-13.0))); dx[idx+5]=(_minf-x[idx+5])/_taum; // differential eqn for hCa2 _hinf= efunc(E,60,6.2); _tauh= 120+(300.0/(exp((E+55.0)/9.0)+exp((E+65.0)/-16))); dx[idx+6]=(_hinf-x[idx+6])/_tauh; // differential eqn for mA _minf= efunc(E,27.2,-8.7); _taum= 23.2-20.8*efunc(E,32.9,-15.2); dx[idx+7]=(_minf-x[idx+7])/_taum; // differential eqn for hA _hinf= efunc(E,56.9,4.9); _tauh= 77.2-58.4*efunc(E,38.9,-26.5); dx[idx+8]=(_hinf-x[idx+8])/_tauh; // differential eqn for mKCa _minf= (x[idx+12]/(x[idx+12]+3.0))*efunc(E,28.3,-12.6); _taum= 180.6-150.2*efunc(E,46.0,-22.7); dx[idx+9]=(_minf-x[idx+9])/_taum; // differential eqn for mKd _minf= efunc(E,12.3,-11.8); _taum= 14.4-12.8*efunc(E,28.3,-19.2); dx[idx+10]=(_minf-x[idx+10])/_taum; // differential eqn for mh _minf= efunc(E,75.0,5.5); _taum= 2.0/(exp((E+169.7)/-11.6)+exp((E-26.7)/14.3)); if (_taum < 1e-2) _taum= 1e-2; // same instability as for _tauh_Na ... dx[idx+11]=(_minf-x[idx+11])/_taum; // differential eqn for o, prob. of Ca mediated K channel activation dx[idx+12]=(-f*(I_CaT+I_CaS)-x[idx+12]+C_Ca_0)/tau_Ca; // the Isyn is calculated elsewhere dx[idx+13]= 0.0; } #undef efunc #undef E #undef Capacit #undef tau_Ca #undef f #undef C_Ca_0 #undef E_Na #undef E_K #undef E_H #undef E_leak #undef g_H #undef g_leak #undef g_Na #undef g_CaT #undef g_CaS #undef g_A #undef g_KCa #undef g_Kd #undef Area #endif