I tried to make a "more realistic" SEClamp, by adding a rate in which the voltage is measured and the current updated (normally something like 20 kHz). See this code (actually only some addition in procedure icur was necessary) :
Code: Select all
NEURON {
POINT_PROCESS SEClamp2
ELECTRODE_CURRENT i
RANGE dur1, amp1, dur2, amp2, dur3, amp3, rs, vc, i, rate
}
UNITS {
(nA) = (nanoamp)
(mV) = (millivolt)
(uS) = (microsiemens)
(kHz) = (/ms)
}
PARAMETER {
rs = 1 (megohm) <1e-9, 1e9>
dur1 (ms) amp1 (mV)
dur2 (ms) <0,1e9> amp2 (mV)
dur3 (ms) <0,1e9> amp3 (mV)
rate = 20 (kHz)
}
ASSIGNED {
v (mV) : automatically v + vext when extracellular is present
i (nA)
vc (mV)
tc2 (ms)
tc3 (ms)
on
tlast (ms)
}
INITIAL {
tc2 = dur1 + dur2
tc3 = tc2 + dur3
on = 0
tlast = t
}
BREAKPOINT {
SOLVE icur METHOD after_cvode
vstim()
}
PROCEDURE icur() {
if (on) {
if ((t - tlast) >= 1 / rate) { : this is new
i = (vc - v)/rs
tlast = t : this is new
}
}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()
}
However, when I use cvode in my model and I use a play vector to command the actual amp for SEClamp, during some (not all) voltage steps with SEClamp, Neuron gives an error like
"CVode... At t = /shortly after actual vc-change/ and h = /some value/ the error test failed repeatedly or with |h| = hmin."
I guess an "at_time()" command is missing somewhere but I do not understand exactly at which time point cvode calculates what, so I do not know where to put that.
Thank you for helping!