Numerical Integration of Izhikevich and HH model neurons (Stewart and Bair 2009)

 Download zip file 
Help downloading and running models
Accession:117361
The Parker-Sochacki method is a new technique for the numerical integration of differential equations applicable to many neuronal models. Using this method, the solution order can be adapted according to the local conditions at each time step, enabling adaptive error control without changing the integration timestep. We apply the Parker-Sochacki method to the Izhikevich ‘simple’ model and a Hodgkin-Huxley type neuron, comparing the results with those obtained using the Runge-Kutta and Bulirsch-Stoer methods.
Reference:
1 . Stewart RD, Bair W (2009) Spiking neural network simulation: numerical integration with the Parker-Sochacki method. J Comput Neurosci 27:115-33 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism:
Cell Type(s): Hodgkin-Huxley neuron;
Channel(s):
Gap Junctions:
Receptor(s): AMPA; Glutamate;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: C or C++ program; MATLAB;
Model Concept(s): Simplified Models; Detailed Neuronal Models; Methods;
Implementer(s): Stewart, Robert [Robert.Stewart at pharm.ox.ac.uk];
Search NeuronDB for information about:  AMPA; Glutamate; Gaba; Glutamate;
/*Utilities for numerical integration of the modified Traub-Miles neuron*/
/*Written by Dr Robert Stewart for Stewart & Bair, 2009*/

/*Derivs routine*/
void tm_derivs(double *y, double *dydx, double *fp){
	double alpha_n, beta_n, alpha_m, beta_m, alpha_h, beta_h;
	/*Parameters in nS for cell of 20,000 micometer surface area (2e-4 cm^2)*/
	double gNa = 20000, gK = 6000, gL = 10;
	double C = 200; /*pF*/
 	double ENa = 50; /*Sodium reversal potential in mV*/
	double EK = -90; /*Potassium reversal potential in mV*/
	double EL = -65; /*Leak reversal potential in mV*/
	double E_ampa = 0; /*AMPA reversal potential in mV*/
	double E_gaba = -80; /*GABA reversal potential in mV*/
	double Vr = -63;
	double E_alpha_n = Vr+15;
	double E_beta_n = Vr+10;
	double E_alpha_m = Vr+13;
	double E_beta_m = Vr+40;
	double E_alpha_h = Vr+17;
	double E_beta_h = Vr+40;

	double v = y[0];
	double n = y[1];
	double m = y[2];
	double h = y[3];
	double g_ampa = y[4];
	double g_gaba = y[5];
	double n4 = n*n*n*n; 
	double m3h = m*m*m*h;
	double dt = fp[0];
	double co_g_ampa = fp[1]; 
	double co_g_gaba = fp[2];
	double I_inj = fp[8];

	dydx[0] = dt*(I_inj - gK*n4*(v-EK) - gNa*m3h*(v-ENa) - gL*(v-EL) - g_ampa*(v-E_ampa) - g_gaba*(v-E_gaba))/C;
	if(v == E_alpha_n) alpha_n = 0.032*5; /*protect against div by zero*/
	else alpha_n = 0.032 * (E_alpha_n-v)/(exp((E_alpha_n-v)/5) - 1);
	beta_n = 0.5*exp((E_beta_n-v)/40);      
	dydx[1] = dt*(alpha_n*(1-n) - beta_n*n);
	if(v == E_alpha_m) alpha_m = 0.32*4; /*protect against div by zero*/
	else alpha_m = 0.32 * (E_alpha_m-v)/(exp((E_alpha_m-v)/4) - 1);
	if(v == E_beta_m) beta_m = 0.28*5; /*protect against div by zero*/
	else beta_m = 0.28 * (v-E_beta_m)/(exp((v-E_beta_m)/5) - 1);
	dydx[2] = dt*(alpha_m*(1-m) - beta_m*m);
	alpha_h = 0.128 * exp((E_alpha_h-v)/18); 
	beta_h = 4 / (exp((E_beta_h-v)/5)+1);
	dydx[3] = dt*(alpha_h*(1-h) - beta_h*h);
	dydx[4] = co_g_ampa*g_ampa;
	dydx[5] = co_g_gaba*g_gaba;
}

void tm_first(double **y, double **co, double *fp){
	double v,n,m,h,a,b,c,d,g_ampa,g_gaba,alpha_n,alpha_m,beta_m,n4,m3h,v1;
	double e,f,g,q,r,s,v_alpha_n,v_alpha_m,v_beta_m,co_alpha_n,co_alpha_m,co_beta_m;
	double chi,psi,xi,co_K,co_Na,I_inj;
	double gNa = 20000, gK = 6000, gL = 10;
 	double ENa = 50; /*Sodium reversal potential in mV*/
	double EK = -90; /*Potassium reversal potential in mV*/
	double EL = -65; /*Leak reversal potential in mV*/
	double E_ampa = 0; /*AMPA reversal potential in mV*/
	double E_gaba = -80; /*GABA reversal potential in mV*/
	
	alpha_n = fp[3]; alpha_m = fp[4]; co_K = fp[5]; co_Na = fp[6];
	I_inj = fp[8];
	v=y[0][0]; n=y[1][0]; m=y[2][0]; h=y[3][0]; g_ampa=y[4][0]; g_gaba=y[5][0];
	a=y[6][0]; b=y[7][0]; c=y[8][0]; d=y[9][0]; 
	e=y[10][0]; f=y[11][0];	g=y[12][0]; q=y[13][0]; r=y[14][0]; s = y[15][0];
	chi = y[16][0]; psi = y[17][0]; xi = y[18][0];
	
	v1 = co[0][0]*(v*chi+co_K*EK+co_Na*ENa+gL*EL+g_ampa*E_ampa+g_gaba*E_gaba+I_inj);
	y[0][1] = v1;
	y[1][1] = co[1][0]*(n*psi + alpha_n);/*n*/
	y[2][1] = co[2][0]*(m*xi + alpha_m);/*m*/
	y[3][1] = co[3][0]*(b*(1-h) - 4*d);/*h*/
	y[4][1] = co[4][0]*g_ampa;
	y[5][1] = co[5][0]*g_gaba;
	y[6][1] = co[6][0]*v1*a;/*a*/
	y[7][1] = co[7][0]*v1*b;/*b*/
	y[8][1] = co[8][0]*v1*c;/*c*/
	y[9][1] = (y[3][1] - d*y[8][1])/(c+1); /*d = h/(c+1)*/
	y[10][1] = co[10][0]*v1*e;/*e*/
	y[11][1] = co[11][0]*v1*f;/*f*/
	y[12][1] = co[12][0]*v1*g;/*g*/
	y[13][1] = (-v1 - q*y[10][1])/(e-1);/*q*/
	y[14][1] = (-v1 - r*y[11][1])/(f-1);/*r*/
	y[15][1] = (v1 - s*y[12][1])/(g-1);/*s*/
}

void tm_iter(double **y, double **co, double *fp, int p){
	double v,n,m,h,g_ampa,g_gaba,a,b,c,d,e,f,g,q,r,s,chi,psi,xi;
	double v_chi,alpha_n,alpha_m,beta_m,co_K,co_Na; 
	double n_psi,m_xi,bh,n4,m3h,cd,eq,fr,gs,temp;
	int i,nv=6;
	double gNa = 20000, gK = 6000, gL = 10;
 	double ENa = 50; /*Sodium reversal potential in mV*/
	double EK = -90; /*Potassium reversal potential in mV*/
	double EL = -65; /*Leak reversal potential in mV*/
	double E_ampa = 0; /*AMPA reversal potential in mV*/
	double E_gaba = -80; /*GABA reversal potential in mV*/
	double co_alpha_n = 0.032;
	double co_alpha_m = 0.32;
	double co_beta_m = 0.28;
	
	v=y[0][0]; n=y[1][0]; m=y[2][0]; h=y[3][0]; g_ampa=y[4][0]; g_gaba=y[5][0];
	a=y[6][0]; b=y[7][0]; c=y[8][0]; d=y[9][0]; 
	e=y[10][0]; f=y[11][0];	g=y[12][0]; q=y[13][0]; r=y[14][0]; s = y[15][0];
	chi = y[16][0]; psi = y[17][0]; xi = y[18][0];
	
	y[21][0] = m*m*m; y[23][0] = n*n*n*n; y[19][0] = n*n; y[24][0] = n*n;
  
  /*psi and n times psi*/
  alpha_n = co_alpha_n*y[13][p];    		
  y[17][p] = -(alpha_n+y[6][p]); /*psi*/
  cauchy_prod(p,y[17],psi,y[1],n,&n_psi); /* n*psi */

  /*xi and m times xi*/
  alpha_m = co_alpha_m*y[14][p];
  beta_m = co_beta_m*y[15][p];
  y[18][p] = -(alpha_m+beta_m); /*xi*/
  cauchy_prod(p,y[18],xi,y[2],m,&m_xi); /* m*xi */

  /*powers*/
  cauchy_prod(p, y[1], n, y[1], n, &y[19][p]); /*n^2*/
  cauchy_prod(p, y[2], m, y[2], m, &y[20][p]); /*m^2*/ 
  cauchy_prod(p, y[20], m*m, y[2], m, &y[21][p]); /*m^3*/ 

  /*testing direct power algorithm from Knuth (p526)*/ 
  /*series_pow(p,y[1],n,y[23],n*n*n*n,4.0); /*n^4*/
  /*series_pow(p,y[23],n*n*n*n,y[24],n*n,0.5); /*(n^4)^.5*/
  /*mexPrintf("n2 error: %3.15f\n",y[24][p]-y[19][p]); /*square root works :-)*/

  cauchy_prod(p,y[7],b,y[3],h-1,&bh);    
  cauchy_prod(p, y[19], n*n, y[19], n*n, &n4); 
  cauchy_prod(p, y[21], m*m*m, y[3], h, &m3h);
  co_K = gK*n4; co_Na = gNa*m3h; 
  y[16][p] = -co_K - co_Na - y[4][p] - y[5][p]; /*chi*/
  cauchy_prod(p,y[16],chi,y[0],v, &v_chi); 

  /*v*/
  y[22][p] = co[0][0]*(v_chi + co_K*EK + co_Na*ENa + y[4][p]*E_ampa + y[5][p]*E_gaba);
  y[0][p+1] = y[22][p]/(p+1); /*v*/

  /*n,m,h,g_ampa,g_gaba*/
  y[1][p+1] = co[1][p]*(n_psi + alpha_n); /*n*/
  y[2][p+1] = co[2][p]*(m_xi + alpha_m); /*m*/
  y[3][p+1] = co[3][p]*(-bh - 4*y[9][p]); /*h*/
  y[4][p+1] = co[4][p]*y[4][p]; /*g_ampa*/
  y[5][p+1] = co[5][p]*y[5][p]; /*g_gaba*/

  /*a,b,c,d*/
  cauchy_prod(p,y[22],y[0][1],y[6],a,&y[6][p+1]); y[6][p+1]*=co[6][p];/*a*/ 
  cauchy_prod(p,y[22],y[0][1],y[7],b,&y[7][p+1]); y[7][p+1]*=co[7][p];/*b*/
  cauchy_prod(p,y[22],y[0][1],y[8],c,&y[8][p+1]); y[8][p+1]*=co[8][p];/*c*/
  series_div(p,y[3][p+1],y[8],(c+1),y[9],d); /*d*/

  cauchy_prod(p,y[22],y[0][1],y[10],e,&y[10][p+1]);y[10][p+1]*=co[10][p];/*e*/ 
  series_div(p,-y[0][p+1],y[10],(e-1),y[13],q); /*q*/

  cauchy_prod(p,y[22],y[0][1],y[11],f,&y[11][p+1]);y[11][p+1]*=co[11][p];/*f*/
  series_div(p,-y[0][p+1],y[11],(f-1),y[14],r); /*r*/

  cauchy_prod(p,y[22],y[0][1],y[12],g,&y[12][p+1]);y[12][p+1]*=co[12][p];/*g*/ 
  series_div(p,y[0][p+1],y[12],(g-1),y[15],s); /*g*/
}

Loading data, please wait...