function [xdot,xinit,option] = chattering( t, x, flag)
% chattering neuron model by Aoyagi et al.

global	stim
global gsk
global gcat
switch flag
case ''
	% state variables
	v=x(1);		% membrane potential (mV)
	m=x(2);		% sodium current activation
	h=x(3);		% sodium current inactivation
	n=x(4);		% potassium current activation
	mca=x(5);	% calcium current activation
	msk=x(6);	% sk current activation
	mcat=x(7);	% cation current activation
	ca=x(8);	% calcium concentration
	oc=x(9);	% ratio of the buffer moecules occupied by Ca2+
	% input current (uA/cm^2)
	Ie = interp1(stim(:,1),stim(:,2),t);	% linear interpolation

	% membrane capacitance (uF/cm^2)
	C=1;
	% maximum conductances (mS/cm^2)
	gna = 130;		% sodium
	gk = 35;		% potassium
	gl = 0.13;		% leak
	% gsk = ;	     	% sk
	% gcat = ;	     	% cation
	% reversal potentials (mV)
	Ena = 60;		% sodium
	Ek = -97;		% potassium
	El = -68.8;		% leak
	Ecat = -42 ;    	% cation
	% other parameters
	Phi = 10;		% temperature factor
	F = 96.5;		% Faraday constant*10^-3 (C/mol)
	T = 293;		% temperature (K)
	R = 8.31;		% gas constant (J/K/mol)
	Pmax = 0.01;		% (uA/(uM mV cm^2))
	Ca_out = 2.1*1000;	% outside calcium concentration (uM)
	KdSK = 0.4;		% half-activation calcium concentration (uM)
	Kdcat = 15;		% half-activation calcium concentration (uM)
	Psi = 2.8;		% inverse of association rate constant (uM ms)
	km = 0.3;		% (1/ms)
	kp = 0.1;		% (1/ms)
	Kmpump = 0.75;		% (uM)
	gpump = 3.6;		% (uM/ms)
	B = 30;			% total concentration of buffer molecules (uM)
	eta = 0.027;

	% activation and inactivation rates
	am = 0.1*(v+25)/(1-exp(-(v+25)/10));
	bm = 4*exp(-(v+50)/12);
	ah = 0.07*exp(-(v+42)/10);
	bh = 1/(1+exp(-(v+12)/10));
	an = 0.01*(v+26)/(1-exp(-(v+26)/10));
	bn = 0.125*exp(-(v+36)/25);
	amca = 1.6/(1+exp(-0.072*(v-5)));
	bmca = 0.02*(v+8.69)/(exp((v+8.69)/5.36)-1);

	% other variables
	Ighk = Pmax*v*(ca-Ca_out*exp(-2*F*v/(R*T)))/(1-exp(-2*F*v/(R*T)));
	m_SK_inf = ca/(ca+KdSK);
	tau_SK = Psi/(ca+KdSK);
	m_cat_inf = ca/(ca+Kdcat);
	tau_cat = Psi/(ca+Kdcat);
	% xdot
	xdot = [ (Ie -gna*m^3*h*(v-Ena) -gk*n^4*(v-Ek) -gl*(v-El) -mca^2*Ighk -gsk*msk*(v-Ek) -gcat*mcat*(v-Ecat))/C;
		Phi*(am*(1-m)-bm*m);
		Phi*(ah*(1-h)-bh*h);
		Phi*(an*(1-n)-bn*n);
		Phi*(amca*(1-mca)-bmca*mca);
		(m_SK_inf-msk)/tau_SK;
		(m_cat_inf-mcat)/tau_cat;
		-eta*mca^2*Ighk+km*B*oc-kp*ca*B*(1-oc)-gpump*ca/(ca+Kmpump);
		-km*oc+kp*ca*(1-oc)];
case 'init'			% Initialization
	xdot = stim(1,1):0.1:stim(end,1);	% time span
	xinit = [ -75; 0.3; 0.7; 0.3; 0.3; 0.3; 0.3; 0.3; 0.3];	% initial values
	option = odeset('MaxStep',10);	% option
end