I have modified the FH mod file to fit the SEF parameters, henceforth referred to as the FHSEF model. I have gone over the code for the past couple of days and the units check out and the equations seem to be correct.
To test my model, I have constructed a hoc file, test.hoc that consists of a single node driven by current clamp. Behavior of the model for both HH and FH models seems normal. I get an action potential of correct shape and the state variables seem to rise and decay as they should, and the currents seem fine as well.
When I test the single node with the FHSEF model, the action potential does not look right. On a plot of "v+vrest", the action potential starts at zero, approaches its peak at 120 mV, but then does not decay back to zero. Instead, the potential stays near the Nerst potential of sodium, 71 mV. Likewise, the state variables do not decay back to their steady-state (inf) values either. m and n stay near 1 (remain active), and h remains near zero (inactive). Plot available in test.hoc
I have checked the steady-state values and time constants of the states by using the Grapher to plot them against v. They have been verified against the same functions coded in Matlab as correct.
If someone could provide some insight into the reason v+vrest hovers at the sodium equilibrium potential, I would much appreciate it. v should return back to its resting value, ie v+vrest should return back to zero.
Frijns' paper shows a comparison of the AP generated by FH and one generated by SEF. Although I can reproduce the FH AP, I cannot do the same for the SEF model. I would like to believe my code is error free =), but that may not be the case. Both fhsef.mod and test.hoc are listed below.
Many thanks for the input!
-----Begin test.hoc---------
Code: Select all
load_file("nrngui.hoc")
// Physical properties
create node
access node
node {
nseg=1
L=2.5
diam=15
insert fhsef
Ra = 70
nai=10 nao=142.0 ki=141.0 ko=4.2
}
// Set temperature and v_init
{celsius=37}
vrest = -84.6
v_init = vrest
// experimental tools
objref stim
stim = new IClamp(.5)
stim.del = .1
stim.dur = .1
stim.amp = 2
ts = 1
// graphical interface appearance
objref g
g = new Graph(0)
g.view(0, 0, 2, 100, 120, 110, 300, 200)
addplot(g,1)
proc fig1() {
g.erase_all()
g.size(0,ts,0,120)
g.addexpr("v-vrest",2,1)
tstop = ts
run()
}
proc fig2() {
g.erase_all()
g.size(0,ts,0,1)
g.addexpr("m_fhsef",2,1)
g.addexpr("h_fhsef",3,1)
g.addexpr("n_fhsef",4,1)
tstop=ts
run()
}
proc fig3() {
g.erase_all
g.size(0,ts,-5,5)
g.addexpr("ina",2,1)
g.addexpr("ik",3,1)
g.addexpr("il_fhsef",5,1)
tstop = ts
run()
}
proc fig4() {
g.erase_all
g.size(0,ts,-5,5)
g.addexpr("ina+ik+il_fhsef",2,1)
g.addexpr("i_cap",3,1)
g.addexpr("ina+ik+il_fhsef+i_cap",5,1)
tstop = ts
run()
}
proc menu() {
xpanel("Choose graph")
xbutton("Membrane Voltage", "fig1()")
xbutton("State variables", "fig2()")
xbutton("Currents", "fig3()")
xbutton("Ionic Currents", "fig4()")
xpanel(0,100)
}
// run simulation
menu()