/*
	Demonstration of the physiological consequences of the 
	cyclic allosteric gating scheme for Ih mediated by HCN2 in 
	thalamocortical relay cells.
	
	This simulation will reproduce	the current-clamp 
	experiment to examine sag and rebound responses to a 
	current step and depolarizing responses to a cAMP pulse 
	described in figure 7 of Wang et al. 2002.  The 
	thalamocortical cell properties are modified from Destexhe 
	et al 1996, with Ih removed and replaced with a kinetic 
	model of HCN2 and the T-type Ca current also removed to 
	focus on the role of the H-current.

	For further details email Matt Nolan at mfnolan@fido.cpmc.
	columbia.edu.
	
	References

	Destexhe, A., Bal, T., McCormick, D.A. and Sejnowski, T.J. 
	(1996). Ionic mechanisms underlying synchronized 
	oscillations and propagating waves in a model of ferret 
	thalamic slices. Journal of Neurophysiology 76, 2049-2070.
	
	Wang, J., Chen, S., Nolan, M.F. and Siegelbaum, S.A. (2002). 
	Activity-dependent regulation of HCN pacemaker channels by 
	cyclicAMP: signalling through dynamic allosteric coupling. 
	Neuron 36(3), 451-461.
*/

load_file("nrngui.hoc")


//  create TC cell

load_file("TCmod.tem")		// read geometry file

objectvar TC[1]
TC[0] = new sTC()

// setup TC cells in resting mode (no spontaneous oscillation)

TC[0].kl.gmax = 0.003
TC[0].soma.gbar_hcn2 = 0.015		// insert hcn2

//  setup simulation parameters

dt = 1		// one is very close to the converged solution at dt = 0.1 however
		// dt=1 gives similar enough results and is good for demonstrating the model
steps_per_ms=1	// unless this is changed from the default of 40 neuron will change
		// dt back to 0.025

delay = 10000	// delay before step or pulse, reduced from 20000 in published figure to improve simulation speed
tstop = 30000	// simulation stop time, changed from 40000
gstart = 8000	// graph start, 18000 in published fig
gstop = 24000	// graph end, 34000 in published fig
celsius = 36
v_init = -77
conc = 0		// cAMP concentration
variable_domain(&conc, 0, 1e9)

//  IClamps, use stim[1] to adjust DC

objref stim[2]

for i=0,1	{
	TC[0].soma stim[i] = new IClamp(0.5)
	stim[i].del=delay
	stim[i].dur=tstop
	stim[i].amp=0
}

//  cAMPclamp to adjust cAMP

TC[0].soma	{insert cAMPclamp}
	
proc resetcAMP()	{
	TC[0].soma	{
		conc1_cAMPclamp=0
		conc2_cAMPclamp=0
		del_cAMPclamp=0
		dur_cAMPclamp=0
	}
}
	
resetcAMP()

// setup recording variables

objref vm[1], state[4]

vm[0] = new Vector (40000)
for i=0,3	{state[i] = new Vector (40000)}
vm[0].record(&TC[0].soma.v(0.5),1)
state[0].record(&TC[0].soma.c_hcn2(0.5),1)	// record closed state
state[1].record(&TC[0].soma.o_hcn2(0.5),1)	// record open state
state[2].record(&TC[0].soma.cao_hcn2(0.5),1)	// record open cAMP bound state
state[3].record(&TC[0].soma.cac_hcn2(0.5),1)	// record closed cAMP bound state

// Set up graph for Em and I step with buttons for running with different cAMP concentrations

objref b, g[2]

b = new VBox()
b.intercept(1)
for i=0,1	g[i] = new Graph()
g[0].size(gstart,gstop,-85,-65)
g[1].size(gstart,gstop,0,1.0)
g[1].color(1)
g[1].label("C")
g[1].color(2)
g[1].label("O")
g[1].color(3)
g[1].label("AO")
g[1].color(4)
g[1].label("AC")
xpanel("")
xbutton("currentsteps", "experiment1()")
xbutton("cAMP pulse", "experiment2()")
xvalue("cAMP concentration (microM)", "conc", 1)
xpanel()
b.intercept(0)
b.map("", 300, 100, 450, 350)

// procedure to plot vectors
proc plotdata ()	{
	vm[0].line(g[0])
	for i=0,3	state[i].line(g[1], i+1, 1)	
}

//experiment 1, single step experiment at various step durations
proc experiment1 ()	{
	stim[1].amp=0
	resetcAMP()
	IClamp[0].del=delay
	IClamp[0].amp=-0.05
	TC[0].soma	{conc1_cAMPclamp=conc*1e-6}
	stim[0].dur=125
	init()
	run()
	plotdata()
	stim[0].dur=1000
	init()
	run()
	plotdata()
	stim[0].dur=4000
	init()
	run()
	plotdata()
}

//experiment 2, look at response to brief step increases in cAMP
proc experiment2 ()	{local i
	for i=0,1	IClamp[i].amp=0
	TC[0].soma {
		conc1_cAMPclamp=0
		del_cAMPclamp=delay
		dur_cAMPclamp=1000
		conc2_cAMPclamp = conc*1e-6
		}
	init()
	run()
	plotdata()
}