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

 Download zip file 
Help downloading and running models
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.
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;
Gap Junctions:
Receptor(s): AMPA; Glutamate;
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];
Search NeuronDB for information about:  AMPA; Glutamate; Gaba; Glutamate;
/*Utilities for numerical integration of the Izhikevich model with COBA synapses*/
/*Written by Dr Robert Stewart for Stewart & Bair, 2009*/

/*Calculates the derivatives*/
void iz_derivs(double *y, double *dydx, double *fp){
	double v,u,g_ampa,g_gaba,I,k,l,E_ampa,E_gaba,E,a,b,co_g_ampa,co_g_gaba;
	v = y[0]; u = y[1]; g_ampa = y[2]; g_gaba = y[3]; /*extract variables*/
	I = fp[0]; k = fp[1]; l = fp[2]; E_ampa = fp[3]; E_gaba = fp[4];
	E = fp[5]; a = fp[6]; b = fp[7]; co_g_ampa = fp[8]; co_g_gaba = fp[9];
	dydx[0] = E*(k*v*v + l*v - u + I - g_ampa*(v-E_ampa) - g_gaba*(v-E_gaba));
	dydx[1] = a*(b*v - u);
	dydx[2] = co_g_ampa*g_ampa;
	dydx[3] = co_g_gaba*g_gaba;

/*PS - first term*/
void iz_first(double **y, double **co, double *fp){ 
	double v,u,g_ampa,g_gaba,I,k,l,E_ampa,E_gaba,E,a,b,co_g_ampa,co_g_gaba,chi;
	v = y[0][0]; u = y[1][0]; g_ampa = y[2][0]; g_gaba = y[3][0]; chi = y[4][0];
	I = fp[0]; k = fp[1]; E_ampa = fp[3]; E_gaba = fp[4];
	E = fp[5]; a = fp[6]; b = fp[7]; co_g_ampa = fp[8]; co_g_gaba = fp[9];
	y[0][1] = E*(v*chi - u + E_ampa*g_ampa + E_gaba*g_gaba + I);
	y[1][1] = a*(b*v - u);	
	y[2][1] = co_g_ampa*g_ampa;
	y[3][1] = co_g_gaba*g_gaba; 
	y[4][1] = k*y[0][1] - y[2][1] - y[3][1];

/*PS - iteration function for higher order terms*/
void iz_iter(double **y, double **co, double *fp, int p){ 
	double v,k,l,E_ampa,E_gaba,b,chi,vchi; int i;
	k = fp[1]; E_ampa = fp[3]; E_gaba = fp[4]; b = fp[7];
	vchi = y[0][0]*y[4][p] + y[4][0]*y[0][p];
	for(i = 1; i < p; i++){vchi += y[0][i]*y[4][p-i];}
	y[0][p+1] = co[0][p]*(vchi - y[1][p] + E_ampa*y[2][p] + E_gaba*y[3][p]);
	y[1][p+1] = co[1][p]*(b*y[0][p] - y[1][p]);
	y[2][p+1] = co[2][p]*y[2][p];
	y[3][p+1] = co[3][p]*y[3][p];
	y[4][p+1] = k*y[0][p+1] - y[2][p+1] - y[3][p+1];

Loading data, please wait...