Problem encountered using Vector.play and somatic traces

[selfdestructo]

Hi,

I'm encountering an issue obtaining the results I expect using the Vector.play method. As in Pettersen & Einevoll (Biophys J, 94, 784-802, 2008), considering some multi-compartment neuron with passive channels on all dendrites and a single active somatic compartment with somatic stimuli, the somatic membrane potential fluctuations will give rise to somatic membrane- and axial currents.

If the somatic potential obtained is played into a new, purely passive but otherwise similar model, the currents that arise should indeed be similar, any somatic membrane voltage fluctuation give rise to currents, the same currents would give equivalent membrane voltage fluctuations, which could be shown analytically for a ball-n-stick, I presume.

With the numerics in NEURON however, this is not the case. For reasonable dt's the outcome of the model appear rather low-pass filtered, evident executing the attached script, without modifying soma.cm.
I've read the documentation for Vector.play() but the means of action I need to take to get agreement is not clear to me. I've found that using ridiculously small dt's, or setting the soma.cm to some near infinite value (say 1E6) reduce the problem, but the traces will always lag behind the trace I played back in the first place.

Any explanation and solution of the issue would be very welcome indeed.

Espen.

Code: Select all

#!/usr/bin/env python

import pylab as pl
from neuron import h 

def initialize():
    h.finitialize(v_init)
    h.fcurrent()

def integrate():
    while h.t < tstop:
        h.fadvance()

def go():
    initialize()
    integrate()

'''Main'''
pl.close('all')

h.dt = 0.05
tstop = 5
v_init = -65


'''Def some ball-n-stick with active soma, synaptic stimulus, record soma Vm'''
soma = h.Section()
dend = h.Section()

dend.connect(soma, 1, 0)

soma.L = 30
soma.nseg = 1
soma.diam = 30

dend.L = 1000
dend.nseg = 100
dend.diam = 3

for sec in h.allsec():
    sec.insert('pas')
    sec.Ra = 100
    sec.cm = 1
    sec.g_pas = 1./20000
    sec.e_pas = -65
    sec.insert('extracellular')
    
soma.insert('hh')

syn = h.AlphaSynapse(0.5, sec=soma)
syn.onset = 0.5
syn.gmax = 0.05
syn.e = 0

tvec = h.Vector()
tvec.record((h._ref_t))

vsoma_0 = h.Vector()
vsoma_0.record(soma(0.5)._ref_v)

isoma_0 = h.Vector()
isoma_0.record(soma(0.5)._ref_i_membrane)

go()


'''Define similar all-passive cell, playback of prev. obtained somapotential in soma,
should yield same somatic response and membrane currents'''
soma = h.Section()
dend = h.Section()

dend.connect(soma, 1, 0)

soma.L = 30
soma.nseg = 1
soma.diam = 30

dend.L = 1000
dend.nseg = 100
dend.diam = 3

for sec in h.allsec():
    sec.insert('pas')
    sec.Ra = 100
    sec.cm = 1
    sec.g_pas = 1./20000
    sec.e_pas = -65
    sec.insert('extracellular')

#soma.cm = 1E6

vsoma_0.play(soma(0.5)._ref_v, h.dt)

vsoma_1 = h.Vector()
vsoma_1.record(soma(0.5)._ref_v)

isoma_1 = h.Vector()
isoma_1.record(soma(0.5)._ref_i_membrane)

go()

'''Plottin'''
pl.figure(1)
pl.plot(tvec, vsoma_0, tvec, vsoma_1)
pl.axis('tight')
pl.xlabel('t')
pl.ylabel('soma potential')
pl.figure(2)
pl.plot(tvec, isoma_0, tvec, isoma_1)
pl.axis('tight')
pl.xlabel('t')
pl.ylabel('soma current')

[ted]

Suggestion and hint: instead of trying to drive soma.v, attach an SEClamp to the soma, set its series resistance parameter to something very small (0.01 or 0.001 megohm would be good), set its dur1 to 1e9, and drive its amp1 parameter. What happens now?

[selfdestructo]

Suggestion and hint: instead of trying to drive soma.v, attach an SEClamp to the soma, set its series resistance parameter to something very small (0.01 or 0.001 megohm would be good), set its dur1 to 1e9, and drive its amp1 parameter. What happens now?

Hi Ted,
I'm aware that SEClamp can be used for this purpose, but for my application it simply does not give any difference compared to setting the cm of soma to "infinity". I wanted to avoid the electrode implementations in general, as I want to extract the membrane currents that arise for further use (I'm working with extracellular potentials), and whenever those electrodes are used, I have to add the electrode current to the membrane currents anyway. Also the response of the latter cell is trailing one time-step anyway, which is ok. I'm still a bit puzzled though that playing straight into soma._ref_v doesn't work though for realistic values of cm.

Code: Select all

...
#soma.cm = 1E6

clamp = h.SEClamp(0.5, sec=soma)
clamp.rs = 0.001
clamp.dur1 = 1E9
vsoma_0.play(clamp._ref_amp1, h.dt)

iclamp = h.Vector()
iclamp.record(clamp._ref_i)

#vsoma_0.play(soma(0.5)._ref_v, h.dt)

vsoma_1 = h.Vector()
vsoma_1.record(soma(0.5)._ref_v)

isoma_1 = h.Vector()
isoma_1.record(soma(0.5)._ref_i_membrane)

go()

'''Plottin'''
pl.figure(1)
pl.plot(tvec, vsoma_0, tvec, vsoma_1)
pl.axis('tight')
pl.xlabel('t')
pl.ylabel('soma potential')
pl.figure(2)
pl.plot(tvec, pl.array(isoma_0) * h.area(0.5, section=soma) / 1E2,
        tvec, pl.array(isoma_1) * h.area(0.5, section=soma) / 1E2,
        tvec, iclamp,
        tvec, pl.array(isoma_1) * h.area(0.5, section=soma) / 1E2 - pl.array(iclamp))
pl.legend(['soma 1 i', 'soma 2 i', 'SEClamp i', 'Net i'], loc='best')
pl.axis('tight')
pl.xlabel('t')
pl.ylabel('soma current')

[ted]

I should have examined the code in your first post more closely. You were thinking that the only difference between your "active" model and your "passive" model was that the latter lacks the voltage-gated currents that are present in the former. However, the way you made the latter passive was not by changing gnabar_hh and gkbar_hh to 0, but by not inserting hh. Which means that the passive model lacked not only the hh mechanism's voltage-gated gna and gk, but also the hh mechanism's huge gleak--huge because it's 0.001 S/cm2. Leave that out, and the soma is a perfect capacitor. Pour charge into it, and it has nowhere to go but into the dendrite. Once in the passive dendrite, it's not going to slosh back into the soma--it's going to have to slowly trickle out the relatively sparse passive channels that are scattered along that long, narrow pipe. That changes the dynamics of the system tremendously.

(It's easy to forget that hh gives you not only gna and gk but also gleak. I do it myself occasionally, "just for practice." But without gleak, hh would't be the complete Hodgkin-Huxley mechanism.)

This would have been revealed by browsing the properties of the cell with the ModelView tool (does it work under Python?), or by using whatever is the Python equivalent of
forall psection()
--probably something like

Code: Select all

forsec in h.allsec():
  h.psection()