I_clamp calculated incorrectly when voltage-clamping a multi-compartment model

[luadam4c]

Just to remind you that there is still this error. I finally got to this error after using x.slice() to make xolotl consider the model multi-compartment

Error using X_cf2eeaab1609998b11ce5d1d973b3915 [xolotl] multi-compartment models cannot be integrated when something is clamped yet.

Error in xolotl/integrate (line 161) [results{1:n_outputs+1}] = f(arguments,self.I_ext',self.V_clamp');

[sg-s]

yes...need to think about how to do this.

[sg-s]

I don't think I know how to do this -- for multi-compartment models, I'm using the Crank-Nicholson scheme described in Appendix A of Chapter 6 in Dayan & Abbott.

If a single compartment was to be clamped, this has to be changed -- not sure how. if you have some pointers let me know

[alec-hoyland]

Yeah, the problem is that, while defining the coefficients (b' and d' in the text) is the same, when you back-propagate to solve for the change in voltage, any clamped voltage breaks the chain. Michael Hines clearly thinks it's possible though with either implicit Euler or Crank-Nicholson.

This paper provides some insights. https://link.springer.com/content/pdf/10.1023%2FA%3A1008945925865.pdf

[luadam4c]

Yeah, the problem is that, while defining the coefficients (b' and d' in the text) is the same, when you back-propagate to solve for the change in voltage, any clamped voltage breaks the chain. Michael Hines clearly thinks it's possible though with either implicit Euler or Crank-Nicholson.

I think in NEURON the single electrode voltage clamp is treated as a current source. So you integrate things the same way as you do for current clamp, but you add a current to wherever the clamp is to compensate for the voltage change. That's exactly how the voltage clamp works anyway.

See NEURON's svclmp.mod:

TITLE svclmp.mod
COMMENT
Single electrode Voltage clamp with three levels.
Clamp is on at time 0, and off at time
dur1+dur2+dur3. When clamp is off the injected current is 0.
The clamp levels are amp1, amp2, amp3.
i is the injected current, vc measures the control voltage)
Do not insert several instances of this model at the same location in order to
make level changes. That is equivalent to independent clamps and they will
have incompatible internal state values.
The electrical circuit for the clamp is exceedingly simple:
vc ---'\/\/`--- cell
        rs

Note that since this is an electrode current model v refers to the
internal potential which is equivalent to the membrane potential v when
there is no extracellular membrane mechanism present but is v+vext when
one is present.
Also since i is an electrode current,
positive values of i depolarize the cell. (Normally, positive membrane currents
are outward and thus hyperpolarize the cell)
ENDCOMMENT

INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}

DEFINE NSTEP 3

NEURON {
    POINT_PROCESS SEClamp
    ELECTRODE_CURRENT i
    RANGE dur1, amp1, dur2, amp2, dur3, amp3, rs, vc, i
}

UNITS {
    (nA) = (nanoamp)
    (mV) = (millivolt)
    (uS) = (microsiemens)
}

PARAMETER {
    rs = 1 (megohm) <1e-9, 1e9>
    dur1 (ms)     amp1 (mV)
    dur2 (ms) <0,1e9> amp2 (mV)
    dur3 (ms) <0,1e9> amp3 (mV)
}

ASSIGNED {
    v (mV)  : automatically v + vext when extracellular is present
    i (nA)
    vc (mV)
    tc2 (ms)
    tc3 (ms)
    on
}

INITIAL {
    tc2 = dur1 + dur2
    tc3 = tc2 + dur3
    on = 0
}

BREAKPOINT {
    SOLVE icur METHOD after_cvode
    vstim()
}

PROCEDURE icur() {
    if (on) {
        i = (vc - v)/rs
    }else{
        i = 0
    }
}

COMMENT
The SOLVE of icur() in the BREAKPOINT block is necessary to compute
i=(vc - v(t))/rs instead of i=(vc - v(t-dt))/rs
This is important for time varying vc because the actual i used in
the implicit method is equivalent to (vc - v(t)/rs due to the
calculation of di/dv from the BREAKPOINT block.
The reason this works is because the SOLVE statement in the BREAKPOINT block
is executed after the membrane potential is advanced.

It is a shame that vstim has to be called twice but putting the call
in a SOLVE block would cause playing a Vector into vc to be off by one
time step.
ENDCOMMENT

PROCEDURE vstim() {
    on = 1
    if (dur1) {at_time(dur1)}
    if (dur2) {at_time(tc2)}
    if (dur3) {at_time(tc3)}
    if (t < dur1) {
        vc = amp1
    }else if (t < tc2) {
        vc = amp2
    }else if (t < tc3) {
        vc = amp3
    }else {
        vc = 0
        on = 0
    }
    icur()
}

[alec-hoyland]

Right -- they use an imperfect voltage clamp where there is a load resistance on the electrode. Then you either have to treat is as a conductance in the compartment that just happens to cancel out the rest of the currents, or treat it as another compartment that then vastly overpowers the clamped compartment.

The alternative is to break up the integration step into solving several submatrices.

The advantage to the former method is that it's simpler to program and imagine. It's no faster than the submatrix route, though CN integration tends to cause spurious oscillations under certain VC conditions.

[luadam4c]

Right -- they use an imperfect voltage clamp where there is a load resistance on the electrode. Then you either have to treat is as a conductance in the compartment that just happens to cancel out the rest of the currents, or treat it as another compartment that then vastly overpowers the clamped compartment.

The alternative is to break up the integration step into solving several submatrices.

The advantage to the former method is that it's simpler to program and imagine. It's no faster than the submatrix route, though CN integration tends to cause spurious oscillations under certain VC conditions.

I think the fact that there is a series resistance for the electrode is actually very realistic. We often measure that series resistance when doing experiments anyway, so it will be helpful to be able to make that a parameter in the simulations. I'm not sure if a two electrode voltage clamp should be different though, as it is indeed a more perfect voltage clamp. But almost all mammalian electrophysiology deals with single electrode voltage clamp, so I think it is worthwhile to implement at least the former version.

[sg-s]

there is some experimental support for voltage clamping a compartment in a multi-compartment neuron model, check out demo_multi_compartment_clamp

[luadam4c]

there is some experimental support for voltage clamping a compartment in a multi-compartment neuron model, check out demo_multi_compartment_clamp

Thanks for making voltage clamp work!

However, when I tried the voltage clamp with my simple model, I get the same current response in the soma whether or not the compartments are connected by setting the tree_idx to 1 in the soma. In fact, when I tried to use the recorded current in my voltage clamp as the current trace in current clamp mode, I can't reproduce the voltage I started with:

load('xolotl_test_before_simulations.mat')

m3ha.V_clamp(3) = -60;
iResponse = m3ha.integrate;
iVClamp = iResponse(:, 3);
m3ha.I_ext = zeros(size(iResponse));
m3ha.I_ext(:, 3) = iVClamp;
vResponse = m3ha.integrate;
m3ha.plot

I'm expecting to see the soma being held at -60 mV

Code is attached

test_m3ha_voltage_clamp.zip

[sg-s]

this code snippet

m3ha.V_clamp(3) = -60;
iResponse = m3ha.integrate;
iVClamp = iResponse(:, 3);
m3ha.I_ext = zeros(size(iResponse));
m3ha.I_ext(:, 3) = iVClamp;
vResponse = m3ha.integrate;
m3ha.plot

does not perform a voltage clamp. it instead operates in current injection mode. you cannot simultaneously inject current and voltage clamp (for now). i'm also curious why you would ever want to do that (as far as i know, this is also unstable dynamically, and you would not do this in a real neuron)

[sg-s]

the model you uploaded behaves exactly as i would expect it to. here, i'm clamping the soma, and plotting the voltages in the other two compartments:

[luadam4c]

this code snippet

m3ha.V_clamp(3) = -60;
iResponse = m3ha.integrate;
iVClamp = iResponse(:, 3);
m3ha.I_ext = zeros(size(iResponse));
m3ha.I_ext(:, 3) = iVClamp;
vResponse = m3ha.integrate;
m3ha.plot

does not perform a voltage clamp. it instead operates in current injection mode. you cannot simultaneously inject current and voltage clamp (for now). i'm also curious why you would ever want to do that (as far as i know, this is also unstable dynamically, and you would not do this in a real neuron)

This code snippet is performing two experiments sequentially, not simultaneously.

%% First experiment: Voltage clamp
% Clamp the soma at -60 mV
m3ha.V_clamp(3) = -60;

% Record the current response in all compartments
iResponse = m3ha.integrate;

% Extract the current response of the soma
iVClamp = iResponse(:, 3);

%% Second experiment: Current clamp
% Initialize the current injection
m3ha.I_ext = zeros(size(iResponse));

% Use the current that supposedly went through the electrode at the soma during the last experiment for this experiment
m3ha.I_ext(:, 3) = iVClamp;

% Record the voltage response in all compartments
vResponse = m3ha.integrate;

% We expect the soma to be clamped at -60 mV, but it's not
m3ha.plot

[luadam4c]

the model you uploaded behaves exactly as i would expect it to. here, i'm clamping the soma, and plotting the voltages in the other two compartments:

What about the current recorded? How do you check if it was correct? Wouldn't you try injecting the same current back in current clamp mode and see if you can get the same voltage response as what you trying to clamp at in the voltage clamp experiment? That's exactly what I have been trying to do.

[sg-s]

ah i see -- interesting. that's a good check.

if you're doing this, you should use closed_loop = false, otherwise you will not be starting from the same initial conditions

[luadam4c]

ah i see -- interesting. that's a good check.

if you're doing this, you should use closed_loop = false, otherwise you will not be starting from the same initial conditions

closed_loop is already set to false in the xolotl object I attached.

Any thoughts on why I can't get the initial voltage back with this method?

[sg-s]

have you tried the exact same experiment with a single compartment? i want to make sure that this is valid, or if the way voltage clamp is done in the single compartment is OK

[luadam4c]

@sg-s Yes I have. It works great with the single compartment model.

[sg-s]

OK, thanks! that narrows it down -- i bet there is a mistake in how the clamping currents are computed. in fact, come to think of it, i suspect it doesn't take into account the extra conductances from coupling to other parts of the neuron.

but the voltage clamping itself should be fine -- what I mean is, clamping a compartment should be correct, and the error is limited to the estimation of the clamping current.

[luadam4c]

OK, thanks! that narrows it down -- i bet there is a mistake in how the clamping currents are computed. in fact, come to think of it, i suspect it doesn't take into account the extra conductances from coupling to other parts of the neuron.

but the voltage clamping itself should be fine -- what I mean is, clamping a compartment should be correct, and the error is limited to the estimation of the clamping current.

The model I attached only has a Leak channel in 3 compartments but no complicated active mechanisms. Perhaps you can start with such a minimal model? The current that flows through the electrode should be equal to the current that flows out of the somatic Leak channel plus the current that flows from the soma to the first dendritic section.

electrode current = leak current + axial current

[sg-s]

it is not true that in a single compartment model, voltage clamping a neuron and then injecting the clamping current reproduces the exact same clamping voltage:

these "escapes" are probably because the clamp isn't fast enough.

[sg-s]

indeed, when we block NaV, they are the same:

[sg-s]

and in multi-compartment models, there is indeed a mismatch, because the clamping current isn't computed correctly:

[sg-s]

added the contribution from axial currents, but it still doesn't match. not sure where the mismatch is coming from.

[alec-hoyland]

@sg-s, can you send me the code you used to make these last few figures? I'll take a stab at this.

[sg-s]

@alec-hoyland demo_multi_compartment_clamp reproduces the last figure and illustrates the problem

[alec-hoyland]

There's even a small error if you voltage clamp to a very simple waveform (e.g. V = 0). I'm going to poke around in the C++ code and see if I find anything. Is there any reason we're not convinced that the CN method is just this error-prone?

[alec-hoyland]

Possible things to look into:

• Should the upstream/downstream conductances include injected current? If the conductance should be the inverse of the input resistance, maybe they should?
• Is the area of the final compartment in a multi-compartment chain being computed correctly? Right now, the area is set to the cylindrical part (surface without top and bottom) for all compartments in multi-compartment models.

[sg-s]

i don't know. i think a good first step is do this in NEURON and see if we get the same thing. if NEURON gets it right and we don't, there is clearly a problem with what we're doing

[alec-hoyland]

I'll do it.

[sg-s]

maybe start with a single compartment and try to reproduce this first?

[alec-hoyland]

I give up. I've spent a week working on this, and I can't reproduce it in NEURON because I can't get the clamps to work with numpy. I cannot even get a minimum working example going. I've had two post-docs who use NEURON look at it, and they don't know either. I gave up when NEURON started complaining that h.run() isn't a recognized command.

Here is my test code

from neuron import h
import numpy as numpy
import matplotlib.pyplot as pyplot

sec     = h.Section(name='sec')
sec.insert('hh')
sec.L = sec.diam = 3 # cm

clamp = h.SEClamp(sec(0.5))
clamp.dur1 = 5 # s
clamp.rs = 0.01

x = numpy.linspace(0, 2 * numpy.pi, 50)
voltage = h.Vector(numpy.sin(x))
voltage.play(clamp._ref_amp1, h.dt)

t = h.Vector()
V = h.Vector()
I = h.Vector()
t.record(h._ref_t)
V.record(sec(0.5)._ref_v)
I.record(clamp._ref_i)

# simulate the model
h.finitialize(-65)
h.dt        = 0.1
h.t         = h.dt
h.tstop     = 5000
h.run()

pyplot.figure()
pyplot.plot(t, V)
pyplot.plot(t, I)
pyplot.show()

Here is my full model

## create a multi-compartment model of a soma and neurite
## voltage clamp the soma, record the rectifying current
## add the rectifying current back into the soma
## to attempt to reproduce voltage-clamped neuron activity

from neuron import h, gui
import matplotlib.pyplot as pyplot
import numpy as np

## Instantiate neuronal compartments

soma = h.Section(name='soma')
dend = h.Section(name='dend')

# properties for the soma
soma.L      = 0.05 / 10     # mm -> cm
soma.diam   = 0.025 / 5     # mm radius -> cm diameter
soma.cm     = 10 / 10       # nF/mm^2 -> μF/cm^2
soma.Ra     = 1e-3 * 1e5    # MΩ⋅mm -> Ω⋅cm
soma.nseg   = 1

# properties for the neurite compartments
dend.L      = 0.35
dend.diam   = 0.1 / 10
dend.cm     = 10 / 10
dend.Ra     = 1e-3 * 1e5
dend.nseg   = 10

## Add conductances

conds       = ['acurrent', 'cas', 'cat', 'hcurrent', 'kca', 'kd', 'na']
gbars       = np.array([104, 11.76, 4.7, 0, 390, 250, 0]) * 1e-4  # S/cm^2
Erevs       = [-80, 30, 30, -20, -80, -80, 50]              # mV

# add all conductances to the Section objects
# add calcium mechanism
for sec in [soma, dend]:
    sec.insert('cad')
    for cond in conds:
        sec.insert(cond)

# update the parameters of the conductances
for sec in [soma, dend]:
    for seg in sec:
        seg.acurrent.gbar = gbars[0]
        seg.cas.gbar = gbars[1]
        seg.cat.gbar = gbars[2]
        seg.hcurrent.gbar = gbars[3]
        seg.kca.gbar = gbars[4]
        seg.kd.gbar = gbars[5]
        seg.na.gbar = gbars[6]
        seg.acurrent.Erev = Erevs[0]
        # seg.cas.Erev = Erevs[1]
        # seg.cat.Erev = Erevs[2]
        seg.hcurrent.Erev = Erevs[3]
        seg.kca.Erev = Erevs[4]
        seg.kd.Erev = Erevs[5]
        seg.na.Erev = gbars[6]

## Generate the voltage waveform

# time parameters
h.dt        = 0.1
h.t         = h.dt
h.tstop     = 5000

# instantiate a single-electrode voltage clamp
t           = h.Vector(1e-3 * np.linspace(h.dt, h.tstop, h.tstop/h.dt))
V_clamp     = h.Vector(np.linspace(0, 30, len(t)) + 30 * np.sin(1e-3 * np.arange(0, len(t)) - 50))
myVClamp    = h.SEClamp(soma(0.5))
V_clamp.play(myVClamp._ref_amp1, t, True)

## Simulate the model while voltage-clamping

# acquire the clamping current
t_vec_v     = h.Vector()
v_vec_v     = h.Vector()
I_clamp     = h.Vector()
t_vec_v.record(h._ref_t)
v_vec_v.record(soma(0.5)._ref_v)
I_clamp.record(myVClamp._ref_i)

# simulate the model
h.finitialize(-65)
h.dt        = 0.1
h.t         = h.dt
h.tstop     = 5000
h.run()

pyplot.figure()
pyplot.plot(t_vec_v, v_vec_v)
pyplot.show()

## Simulate the model while current clamping

# storage vectors
t_vec_i     = h.Vector()
v_vec_i     = h.Vector()
t_vec_i.record(h._ref_t)
v_vec_i.record(soma(0.5)._ref_v)

# set up the current clamp
myIClamp    = h.IClamp(soma(0.5))
I_clamp.play(myIClamp._ref_i, t_vec_i, True)

# simulate the model
h.finitialize(-65)
h.dt        = 0.1
h.t         = h.dt
h.tstop     = 5000
h.run()

## Visualize the results

h.dt        = 0.1
h.t         = h.dt
h.tstop     = 5000

pyplot.figure
pyplot.plot(t_vec_v, v_vec_v)

# pyplot.figure(figsize=(8,4))
# pyplot.plot(t_vec_v, v_vec_v);
# pyplot.plot(t_vec_i, v_vec_i)
# pyplot.xlabel('time (ms)')
# pyplot.ylabel('mV')
# pyplot.show()

[sg-s]

fair enough. i'll see if i can take a crack at it

[alec-hoyland]

I pored over the code in compartment.hpp and rewrote several sections for clarity, including the computeClampingCurrent() function that is responsible for returning the clamping current of a multi-compartment model to directly follow Dayan & Abbott 2005. My results are here: #370 .

Here is the result at dt = 0.1

If you take the demo_multi_compartment_clamp.m script and turn sim_dt and dt way down (to about 0.001 ms is fine), you will see the spurious oscillations which are common in Crank-Nicolson simulations, however, you will also see, in the long-time limit, that the solutions converge.

Here is an example of the pretty part.

This does not mean that the CN method is numerically unstable. The CN method is formally stable, but spurious oscillations and other artifacts are known to appear, especially with voltage clamp (see Hines & Carnevale 1995) and the Moore et al. papers.

[sg-s]

this looks good! i think this makes sense. @luadam4c what do you think?

[sg-s]

@alec-hoyland : when you changed the sim_dt you accidentally also changed the frequency of the sine wave, and that's what causes your initial transient. here's the same script, but using the same sine wave. note there is no transient when the solution is off (it's always exactly correct) i used sim_dt=.01

[sg-s]

fixed in master