#pragma rtGlobals=1		// Use modern global access method.

//	Procedures for performing Hodgkin Huxley simulations of excitable membrane.
//	The simulations are based on:
//		Hodgkin, A.L. & Huxley, A.F. (1952) "A quantitative description of membrane current and its 
//		application to conduction and excitation in nerve" Journal of Physiology, vol. 117, pp. 500-544
//	Created 25/9/2000 - JBF
//	For more information about these procedures, or to report bugs, please contact:
//		James.Fallon@ieee.org

//	Performs simulation of membrane voltage based on the Hodgkin Huxley equations.
//	The current waveform is supplied and should be scaled in 's' and 'A'. An appropriate form is:
//		make /o/n=1000 Current
//		setscale /p x 0, 2e-5, "s", Current
//		setscale d 0, 0, "A", Current
//		Current = ((x > 1e-3) - (x > 2e-3))*10e-6
//	which is equivalent to a current pulse of 10µA for 1 ms. 
//	The rate constants are adjusted according to the temperature, Temp, in ¡C. 
//	The method used determines what other wavesforms are calculated:
//	Method	1:	Voltage only (returned in current datafolder)
//			2:	Voltage and current wavesforms (returned in 'Hodgkin Huxley' data folder)
Function SimulateHodgkinHuxley(Current, Temp, Method)
	Wave Current
	Variable Temp, Method
	string currentfldr = getdatafolder(1)
	newdatafolder /o/s root:Packages
	newdatafolder /o/s 'Hodgkin Huxley'
	make /o w_para={1,120,-115,36,12,0.3,-10.613, Temp}
	//	Copy current wave ...
	duplicate /o Current, Voltage, w_Current, w_m, w_h, w_n
	setscale /p x 0,deltax(Current)*1000,"ms",Voltage, w_Current, w_m, w_h, w_n
	w_Current *= -1e6
	switch (Method)
		case 2:	//	Voltage and currents ...
			duplicate /o Current INa, IK, Il
			setscale /p x 0,deltax(Current)*1000,"ms", INa, IK, Il
			break
		case 1:
		default:
			break
	endswitch
	//	Initialise starting values ...
	w_m[0]=0.052971; w_h[0]= 0.59665; w_n[0]= 0.3179
	voltage[0]= 0
	//	Perform simulation ...
	integrateODE /q=1/u=(numpnts(Voltage)) Hodgkin_Huxley, w_para, {Voltage, w_m, w_h, w_n}
	if (v_flag == 1)
		return v_flag
	endif
	//	Return resulting voltage wave ...
	duplicate /o Voltage $currentfldr+"Voltage"
	wave Voltage = $currentfldr+"Voltage"
	Voltage = -(voltage+60)/1000; setscale d 0,0,"V",Voltage
	setscale /p x 0,deltax(Current),"s",Voltage			
	switch (Method)
		case 2:	//	Voltage and currents ...
			INa = w_para[1]*(w_m^3)*w_h*(voltage-w_para[2])
			IK = w_para[3]*(w_n^4)*(voltage-w_para[4])
			Il = w_para[5]*(voltage-w_para[6])
			break
		case 1:
		default:
			break	
	endswitch
	setdatafolder $currentfldr
	return 0
End	//	SimulateHodgkinHuxley(Current, Temp, Method)

//	Function used to solve the Hodgkin Huxley equations for an excitable membrane.
//	The function calculates the derivatives of V, m, h, and n. The curent value is passed
//	in w_Y and the resulting derivative is returned in w_dY. The order in w_Y is:
//		w_Y[0] = V
//		w_Y[1] = m
//		w_Y[2] = h
//		w_Y[3] = n
//	Extra parameters can be passed in w_Para and theses are:
//		w_para[0] = Cm
//		w_para[1] = GNa
//		w_para[2] = VNa
//		w_para[3] = GK
//		w_para[4] = VK
//		w_para[5] = Gl
//		w_para[6] = Vl
//		w_para[7] = Temp
Function Hodgkin_Huxley(w_para, t, w_Y, w_dY)
	Wave w_para, w_Y, w_dY
	variable t
	//	Load local variables ...
	variable V = w_Y[0], m= w_Y[1], h= w_Y[2], n= w_Y[3]
	//	Reference current wave ...
	wave w_Current = w_Current
	//	Calculate the derivatives ...
	//	dV/dt
	w_dY[0] = (w_Current(t)- w_para[5]*(V-w_para[6]) - w_para[1]*m^3*h*(V-w_para[2]) - w_para[3]*n^4*(V-w_para[4]))/w_para[0]
	//	dm/dt
	w_dY[1] = (0.1*(V+25) / (exp((V+25)/10)-1) * (1-m) - 4*exp(V/18)*m)*3^((w_para[7]-6.3)/10)
	//	dh/dt
	w_dY[2] = (0.07*exp(V/20)*(1-h) - h/(exp((V+30)/10)+1))*3^((w_para[7]-6.3)/10)
	//	dn/dt
	w_dY[3] = (0.01*(V+10)/(exp((V+10)/10)-1)*(1-n) - 1/8*exp(V/80)*n)*3^((w_para[7]-6.3)/10)
End	//	Hodgkin_Huxley(w_para, t, w_Y, w_dY)