The integration method and results of one compartment with only HH channel

[gaoyy18]

I constructed a single section (one single compartment) with HH channel. I tried to verify the implicit method when updating the voltage and voltage-gated variables m, n and h.

Code: Select all

# single compartment with HH channel, and use implicit euler method by setting h.secondorder=0

h.celsius=6.3
h.secondorder = 0
class BallAndStick:
    def __init__(self, gid):
        self._gid = gid
        self._setup_morphology()
        self._setup_biophysics()

    def _setup_morphology(self):
        self.soma = h.Section(name="soma", cell=self)
        self.all = self.soma.wholetree()
        self.soma.L = 100 * μm
        self.soma.diam = 10 * µm

    def _setup_biophysics(self):
        for sec in self.all:
            print(sec)
            sec.Ra = 100  # Axial resistance in Ohm * cm
            sec.cm = 1  # Membrane capacitance in micro Farads / cm^2
        self.soma.insert("hh")
        for seg in self.soma:
            seg.hh.gnabar = 0.12  # Sodium conductance in S/cm2
            seg.hh.gkbar = 0.036  # Potassium conductance in S/cm2
            seg.hh.gl = 0.0003  # Leak conductance in S/cm2
            seg.hh.el = -54.3 * mV  # Reversal potential

    def __repr__(self):
        return "BallAndStick[{}]".format(self._gid)

# del my_cell
my_cell = BallAndStick(0)

Then I run the simulation and record related data at each time step to verify the method used when updating voltage and gating variables m, n and h.

Code: Select all

# stim
stim1 = h.IClamp(my_cell.soma(0.5))
stim1.get_segment()
stim1.delay = 0
stim1.dur = 4
stim1.amp = 0.3

# record
soma_v = h.Vector().record(my_cell.soma(0.5)._ref_v)
t = h.Vector().record(h._ref_t)

mvec = h.Vector().record(my_cell.soma(0.5).hh._ref_m)
nvec = h.Vector().record(my_cell.soma(0.5).hh._ref_n)
hvec = h.Vector().record(my_cell.soma(0.5).hh._ref_h)

gkvec = h.Vector().record(my_cell.soma(0.5).hh._ref_gk)
gnavec = h.Vector().record(my_cell.soma(0.5).hh._ref_gna)

ikvec = h.Vector().record(my_cell.soma(0.5)._ref_ik)
inavec = h.Vector().record(my_cell.soma(0.5)._ref_ina)
ilvec = h.Vector().record(my_cell.soma(0.5).hh._ref_il)

mtauvec = h.Vector().record(my_cell.soma(0.5).hh._ref_mtau)
minfvec = h.Vector().record(my_cell.soma(0.5).hh._ref_minf)

# start simulation
h.finitialize(-65 * mV)
h.continuerun(25 * ms)


# convert h.Vector to np.array
s1 = np.array(soma_v.to_python())
t = np.array(t.to_python())

mvalue = np.array(mvec.to_python())
nvalue = np.array(nvec.to_python())
hvalue = np.array(hvec.to_python())

gkvalue = np.array(gkvec.to_python())
gnavalue = np.array(gnavec.to_python())

ikvalue = np.array(ikvec.to_python())
inavalue = np.array(inavec.to_python())
ilvalue = np.array(ilvec.to_python())

mtauvalue = np.array(mtauvec.to_python())
minfvalue = np.array(minfvec.to_python())

Then I tested some relationships between values acquired above.
1.Based on the following test, I verify that gna[t+1] = gnabar*m[t]*m[t]*m[t]*h[t]

Code: Select all

print(gnavalue[0])
print(gnavalue[1])
print(gnavalue[2])
print("**")
print(np.power(mvalue[0], 3)* hvalue[0]*0.12)
print(np.power(mvalue[1], 3)* hvalue[1]*0.12)

I got the following results, and then I found matches mentioned above:

Code: Select all

1.0609192838829855e-05
1.0609192838829855e-05
1.069770674641954e-05
**
1.0609192838829853e-05
1.069770674641954e-05

2. And based on the following test, I verify that the voltage v[t+1] is updated from v[t] with implicit method, and gna[t+1] gk[t+1] are used when updating.
it means c1/dt*(v[t+1]-v[t]) + ihh = i1 (i1 is a constant external input)
--> ihh = i1 - c1/dt*(v[t+1]-v[t]) = itest (itest is a newly defined variable for simplicity in the following codes)
--> ihh = gna[t+1]*(v[t+1]-Ena) + gk[t+1]*(v[t+1]-Ek) + gl*(v[t+1]-El) = i1 - c1/dt*(v[t+1]-v[t]) = itest

Code: Select all

# constant parameters used above are listed here for test
area1 = my_cell.soma(0.5).area() * 1e-8     # in cm^2
i1 = stim1.amp * 1e-9                       # constant external input
c1 = area1 * my_cell.soma.cm * 1e-6         # capacitance

dt = 0.025 * 1e-3
v0 = -65 * 1e-3

# for i=0 to i=1
print("for i=0 to i=1:")
itest=i1-c1/dt*(s1[1]-s1[0])*1e-3
print(itest)

print("**")

ihh=((gnavalue[1]*(s1[1]-50))+(gkvalue[1]*(s1[1]+77))+0.0003*(s1[1]+54.3))*1e-3*area1
print(ihh)

I got the following results, and then I found matches mentioned above:

Code: Select all

for i=0 to i=1:
4.058046768907709e-12
**
4.058046768874686e-12

3. However, I can not find similar implicit differential equations used when updating the voltage-gated variables m, n and h
I thought the implicit method is:
(m[t+1]-m[t])/dt = (minf[t+1]-m[t+1])/mtau[t+1]
however it is not right, as the left and right sides of the above equation do not match:

Code: Select all

# for i=0 to i=1
left = (mvalue[1]-mvalue[0])/dt *1e-3    
print(left)
right = (minfvalue[1]-mvalue[1])/mtauvalue[1]
print(right)

I got not matched results:

Code: Select all

0.005900810542734114
0.0055982019471185465

My question is: what's the correct equation when updating the gating variables m, n and h, from m to m[i+1](taking m as example)?
It really troubled me a lot. Thank you for your guidance in advance.

[gaoyy18]

I realized the problem, please look at this: viewtopic.php?t=4567