Convergence of e_extracellular and membrane currents
Moderator: wwlytton

 Posts: 14
 Joined: Fri Nov 25, 2016 12:02 pm
Convergence of e_extracellular and membrane currents
Dear all,
I am using NEURON to model ephaptic coupling between axons in a nerve. My model necessarily requires coupling NEURON to a Poisson solver (FEM, although for minimal working examples, I just use phi = 1/(4*pi*sigma) * sum(membr.currents/distance) for calculating the extracellular potential from the membrane currents of axons).
At every time step in a simulation, I need an iterative process that brings the system to convergence. Specifically, I am trying to make sure that the membrane currents and the extracellular potential field (phi) are consistent with each other before moving on to the next time step.
To this end, my strategy is the one described below. The process that goes on right before every time step of the simulations I perform has roughly the following structure:

for (every time step):
h.advance()
svstateobj.save()
convergence = False
while not convergence:
1. Calculate phi from the membrane currents and store it in an array.
2. Play the corresponding phi on _ref_e_extracellular on every segment of every axon (note that the time series for every segment includes only one value, since it is only meant .
3. h.finitialize() (otherwise, I think that playing the computed phi into the axons is worhless)
4. svstateobj.restore(1) # The last saved state. For the first iteration of the convergence process, this is the state right after the last h.advance()
5. h.fcurrent()
6. Save state (svstateobj.save()).
7. Check if the values of phi and the membrane currents have converged (I have a certain tolerance value). If so, convergence = True and exit the while loop to move on to the next time step.

(Every simulation I run is actually a chain of simulations of length dt (the length of the time step).)
My hope was that h.fcurrent() recalculated the membrane currents having into account the played new values of phi, but it seems to ignore them. I think it makes sense, since as said in the documentation, it takes into account the STATES (which I guess they are m, n and h for the HH model), when a played vector is not such thing, but an event (am I right or did I misunderstand?).
So my question is: is there a way to make fcurrent acknowledge the pressence of phi so that the new membrane currents are calculated according to the influence of phi on the axon membrane at that time?
As a last resource, I could do this by other means (maybe by translating phi into currents on every segment like McIntyre & al. did), although if there is a way to do it with NEURON, that would make things much easier.
I would appreciate if you had the time to shed some light to me about this.
All the best.
I am using NEURON to model ephaptic coupling between axons in a nerve. My model necessarily requires coupling NEURON to a Poisson solver (FEM, although for minimal working examples, I just use phi = 1/(4*pi*sigma) * sum(membr.currents/distance) for calculating the extracellular potential from the membrane currents of axons).
At every time step in a simulation, I need an iterative process that brings the system to convergence. Specifically, I am trying to make sure that the membrane currents and the extracellular potential field (phi) are consistent with each other before moving on to the next time step.
To this end, my strategy is the one described below. The process that goes on right before every time step of the simulations I perform has roughly the following structure:

for (every time step):
h.advance()
svstateobj.save()
convergence = False
while not convergence:
1. Calculate phi from the membrane currents and store it in an array.
2. Play the corresponding phi on _ref_e_extracellular on every segment of every axon (note that the time series for every segment includes only one value, since it is only meant .
3. h.finitialize() (otherwise, I think that playing the computed phi into the axons is worhless)
4. svstateobj.restore(1) # The last saved state. For the first iteration of the convergence process, this is the state right after the last h.advance()
5. h.fcurrent()
6. Save state (svstateobj.save()).
7. Check if the values of phi and the membrane currents have converged (I have a certain tolerance value). If so, convergence = True and exit the while loop to move on to the next time step.

(Every simulation I run is actually a chain of simulations of length dt (the length of the time step).)
My hope was that h.fcurrent() recalculated the membrane currents having into account the played new values of phi, but it seems to ignore them. I think it makes sense, since as said in the documentation, it takes into account the STATES (which I guess they are m, n and h for the HH model), when a played vector is not such thing, but an event (am I right or did I misunderstand?).
So my question is: is there a way to make fcurrent acknowledge the pressence of phi so that the new membrane currents are calculated according to the influence of phi on the axon membrane at that time?
As a last resource, I could do this by other means (maybe by translating phi into currents on every segment like McIntyre & al. did), although if there is a way to do it with NEURON, that would make things much easier.
I would appreciate if you had the time to shed some light to me about this.
All the best.
Last edited by mcapllonchj on Mon Feb 20, 2017 11:45 am, edited 1 time in total.
Re: Convergence of e_extracellular and membrane currents
Yes, one problem with the iterative method as presented is that
in the "while not convergece" loop, no states
are changing. You might try replacing h.fcurrent with h.fadvance()
(getting rid of the h.advance at the beginning of the time step)
and eliminate step 6
I am presuming that your "Play" is
implemented with PtrVector to scatter the (spatial) vector of phi to the
values of e_extracellular. In that case you should do it just before your
step 5 (which I suggest be h.fadvance).
This is not completely satisfactory as instead of iterating
y(t+dt)_(i+1) = y(t) + dt*f(y(t+dt)_i) to convergence, you would be
iterating
y(t+dt)_(i+1) = y(t) + dt*f(v(t+dt)_i, all_other_states(t))
to convergence but maybe that is good enough for your purposes.
I'm afraid it is not possible, without changing the internal implementation
of NEURON's fixed step integration, to get a completely consistent
implicit method by iterations, so that
y(t+dt)_(i+1) = y(t) + dt*f(y(t+dt)_i) converges (for every state
including things like v, m, h, n). If it were not
for the presence of extracellular which changes the ODE sytem to a
DifferentialAlgebraic Equation (DAE) system it would be convenient to
use python or hoc along with cvode.f to do the iterations directly.
I can see that it would be worthwhile to extend cvode.f to work for DAE
which in neuron is restricted to c*y = f(y) (note that c may have off
diagonal elements and some rows may be 0)
Alternatively, since e_extracellular is a forcing function in this part
of the problem, it would be possible (again with internal implementation
changes to that v_internal = v_membrane + e_force) and then, since model
is still ODE, scatter phi to e_force (note that there is a relatively
new feature called cvode.use_fast_imem that calculates i_membrane without
requiring e_extracellular (very useful for faster LFP calculations))
Cvode (in its DAE mode) does solve your problem if you
scatter phi to e_extracellular on every evaluation
(using a python callback from a POINT_PROCESS BEFORE BREAKPOINT block.
(I believe, although I'm not positive, that i_membrane would be from the
previous internal IDA solver iteration)
in the "while not convergece" loop, no states
are changing. You might try replacing h.fcurrent with h.fadvance()
(getting rid of the h.advance at the beginning of the time step)
and eliminate step 6
I am presuming that your "Play" is
implemented with PtrVector to scatter the (spatial) vector of phi to the
values of e_extracellular. In that case you should do it just before your
step 5 (which I suggest be h.fadvance).
This is not completely satisfactory as instead of iterating
y(t+dt)_(i+1) = y(t) + dt*f(y(t+dt)_i) to convergence, you would be
iterating
y(t+dt)_(i+1) = y(t) + dt*f(v(t+dt)_i, all_other_states(t))
to convergence but maybe that is good enough for your purposes.
I'm afraid it is not possible, without changing the internal implementation
of NEURON's fixed step integration, to get a completely consistent
implicit method by iterations, so that
y(t+dt)_(i+1) = y(t) + dt*f(y(t+dt)_i) converges (for every state
including things like v, m, h, n). If it were not
for the presence of extracellular which changes the ODE sytem to a
DifferentialAlgebraic Equation (DAE) system it would be convenient to
use python or hoc along with cvode.f to do the iterations directly.
I can see that it would be worthwhile to extend cvode.f to work for DAE
which in neuron is restricted to c*y = f(y) (note that c may have off
diagonal elements and some rows may be 0)
Alternatively, since e_extracellular is a forcing function in this part
of the problem, it would be possible (again with internal implementation
changes to that v_internal = v_membrane + e_force) and then, since model
is still ODE, scatter phi to e_force (note that there is a relatively
new feature called cvode.use_fast_imem that calculates i_membrane without
requiring e_extracellular (very useful for faster LFP calculations))
Cvode (in its DAE mode) does solve your problem if you
scatter phi to e_extracellular on every evaluation
(using a python callback from a POINT_PROCESS BEFORE BREAKPOINT block.
(I believe, although I'm not positive, that i_membrane would be from the
previous internal IDA solver iteration)

 Posts: 14
 Joined: Fri Nov 25, 2016 12:02 pm
Re: Convergence of e_extracellular and membrane currents
Dear Hines,
Thank you very much for your response.
Actually I had not discovered PtrVector before. So far I was using:
I have used PtrVector together with pv.scatter in my code now and simulations are so much faster. Thank you for “pointing” that out to me.
I am new in working with CVode and I am now getting familiar with it to try your suggestion. I think it will make a noticeable difference.
Having said that, I think that the convergence process I am trying to implement pursues a different type of convergence. What I am trying to reach is electromagnetic convergence, as in reaching a consistent state of the system in the following way:
At each iteration, the computed phi is different from that of the previous iteration. Assuming v_in does not change, the membrane potential is changed by the same amount as phi. Then, the membrane (ionic) currents are also changed (and I think phi itself would generate new currents too), but I want to leave the states (and hence, the conductances) constant for all iterations. The reason for this is that the time scale of electromagnetic changes is too small to change the states, and these states come from the (dm/dt, dn/dt, dh/dt) computed in the previous time step, when they were calculated using the “correct” v. So I want to hold time constant during the convergence process.
Still having doubts on how this would work in the very first iteration, I must say.
I am thinking about this approach of advancing t to t+dt until reaching convergence. Does it aim to take into account physiological effects of phi that have been not caught by the length of the time step? I can see this point, but I am not very sure of applying phi(t+dt) as a forcing back in t=t (the future modifying the present state, if I did not misunderstand). How can I better understand this process?
Again, I am very thankful for your help.
All the best.
Thank you very much for your response.
Actually I had not discovered PtrVector before. So far I was using:
Code: Select all
for j, seg in enumerate(seglist):
phi[j].play(seg._ref_e_extracellular, h.dt)
I am new in working with CVode and I am now getting familiar with it to try your suggestion. I think it will make a noticeable difference.
Having said that, I think that the convergence process I am trying to implement pursues a different type of convergence. What I am trying to reach is electromagnetic convergence, as in reaching a consistent state of the system in the following way:
At each iteration, the computed phi is different from that of the previous iteration. Assuming v_in does not change, the membrane potential is changed by the same amount as phi. Then, the membrane (ionic) currents are also changed (and I think phi itself would generate new currents too), but I want to leave the states (and hence, the conductances) constant for all iterations. The reason for this is that the time scale of electromagnetic changes is too small to change the states, and these states come from the (dm/dt, dn/dt, dh/dt) computed in the previous time step, when they were calculated using the “correct” v. So I want to hold time constant during the convergence process.
Still having doubts on how this would work in the very first iteration, I must say.
I am thinking about this approach of advancing t to t+dt until reaching convergence. Does it aim to take into account physiological effects of phi that have been not caught by the length of the time step? I can see this point, but I am not very sure of applying phi(t+dt) as a forcing back in t=t (the future modifying the present state, if I did not misunderstand). How can I better understand this process?
Again, I am very thankful for your help.
All the best.
Re: Convergence of e_extracellular and membrane currents
Using savestate.restore will put all NEURON states and voltages back to the saved time (t) and fadvance will put them all to t+dt consistent with the forcing function phi. But realize this is all only first order
correct with the implicit fixed step method. It would be nice to verify that if phi does not change between iterations, that one always gets the same result on each iteration for all the states.
I am guessing you want to recalculate phi on each iteration so that it is consistent with the i_membrane at the end of the previous iteration. Since i_membrane is calculated
after the voltages are updated at the end of fadvance, your phi(t+dt) will ultimately be consistent with v(t+dt) in an implicit solver sense.
correct with the implicit fixed step method. It would be nice to verify that if phi does not change between iterations, that one always gets the same result on each iteration for all the states.
I am guessing you want to recalculate phi on each iteration so that it is consistent with the i_membrane at the end of the previous iteration. Since i_membrane is calculated
after the voltages are updated at the end of fadvance, your phi(t+dt) will ultimately be consistent with v(t+dt) in an implicit solver sense.

 Posts: 14
 Joined: Fri Nov 25, 2016 12:02 pm
Re: Convergence of e_extracellular and membrane currents
Many thanks for this. I will have a think to fully understand it and then implement it to try.
Just in case and for the various experiments I am going to do both with keeping time constant and using fadvance (specifically for the former, actually), is there a way in NEURON to force variables like the membrane currents (which I believe is an output and can not be modified with a scatter) to have specific values of my choice without doing fadvance?
Just in case and for the various experiments I am going to do both with keeping time constant and using fadvance (specifically for the former, actually), is there a way in NEURON to force variables like the membrane currents (which I believe is an output and can not be modified with a scatter) to have specific values of my choice without doing fadvance?

 Site Admin
 Posts: 5750
 Joined: Wed May 18, 2005 4:50 pm
 Location: Yale University School of Medicine
 Contact:
Re: Convergence of e_extracellular and membrane currents
Ah, but v_in does change. Any change of extracellular potential is instantly coupled to the inside of the cell by membrane capacitance, with no attenuation. This must be the case, because membrane capacitance prevents membrane potential from changing instantaneously (short of something producing an infinitely large local membrane current, and that's not going to happen). So rapid fluctuations of extracellular potential produce rapid fluctuations of v_in, and the amplitude of the fluctuations of v_in approach those of extracellular potential as frequency rises above 1/(2*PI*local membrane time constant). Think "no significant attenuation for frequencies > 5/(2*PI*local membrane time constant)."At each iteration, the computed phi is different from that of the previous iteration. Assuming v_in does not change, the membrane potential is changed by the same amount as phi.

 Posts: 14
 Joined: Fri Nov 25, 2016 12:02 pm
Re: Convergence of e_extracellular and membrane currents
Dear Ted, many thanks for that reply (and sorry for this late reply). I see your point. I will come back to this question at some point later on and let you know if I am successful implementing the convergence process.