Spurious voltage changes

Moderator: wwlytton

Post Reply
peterbratby
Posts: 6
Joined: Fri Jul 20, 2018 1:00 am

Spurious voltage changes

Post by peterbratby »

Hi,

I am trying to resolve an issue using NEURON model gap junctions, in which spurious voltage changes occur in gap junction-connected sections.

The following Python code demonstrates the issue. There are three identical sections, each receiving an identical depolarizing current. The first two sections are connected to each other by a gap junction. I would expect all three sections to show the same change in voltage.

However, as can be seen in the voltage traces, the gap junction-connected sections show a smaller voltage change than the isolated section.

The two gap junction-connected sections have identical voltages throughout the simulation, so I would not expect a current to flow through the gap junction. Indeed, the values of gap_junction1._ref_i and gap_junction2._ref_i remain zero throughout.

I noticed that the issue depends on the step size h.dt. Decreasing the step size decreases the voltage difference.

I have also reproduced this issue using parallel transfer to model gap junctions.

Neuron version: VERSION 7.5 master (6b4c19f)

I hope someone can help.

Code: Select all

from neuron import h
import matplotlib.pyplot as plt

g = 10.

h.load_file("stdrun.hoc")
h.tstop = 100


section1 = h.Section()
stim1 = h.IClamp(section1(0.5))
stim1.delay = 0
stim1.amp = 1
stim1.dur = 100

v1 = h.Vector()
v1.record(section1(0.5)._ref_v)

section2 = h.Section()
stim2 = h.IClamp(section2(0.5))
stim2.delay = 0
stim2.amp = 1
stim2.dur = 100

v2 = h.Vector()
v2.record(section2(0.5)._ref_v)

gap_junction1 = h.gap_pointer(section1(0.5))
gap_junction2 = h.gap_pointer(section2(0.5))    
h.setpointer(section2(0.5)._ref_v, 'vgap', gap_junction1)
h.setpointer(section1(0.5)._ref_v, 'vgap', gap_junction2)
gap_junction1.g = gap_junction2.g = g

i1 = h.Vector()
i1.record(gap_junction1._ref_i)
i2 = h.Vector()
i2.record(gap_junction2._ref_i)

section3 = h.Section()
stim3 = h.IClamp(section3(0.5))
stim3.delay = 0
stim3.amp = 1
stim3.dur = 100

v3 = h.Vector()
v3.record(section3(0.5)._ref_v)

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

h.dt=2.5
h.finitialize(-60)

h.run()

plt.figure()
plt.plot(t, v1)
plt.plot(t, v2)
plt.plot(t, v3)

plt.figure()
plt.plot(t, i1)
plt.plot(t, i2)

Code: Select all

NEURON {
	POINT_PROCESS gap_pointer
	NONSPECIFIC_CURRENT i
	RANGE g, i
	POINTER vgap
}
PARAMETER {

	g = 2e-3 (microsiemen)
}
ASSIGNED {
	i (nanoamp)
	vgap (millivolt)
	v (millivolt)
}
BREAKPOINT {
	i = (v - vgap)*g
}
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Spurious voltage changes

Post by ted »

Excellent question. Very observant, especially your note that reducing dt makes the difference smaller. You have created an example that prominently reveals the kind of error that can occur when a gap junction is implemented in such a way that off-diagonal terms are omitted from the system's Jacobian matrix. Not only were your model's sections entirely capacitive (which is a bit of a perversity--I can't think of a single published model implemented with NEURON that has such sections, but then I'm not familiar with all published NEURON models), you also increased the time step by two orders of magnitude from NEURON's default value of 0.025 ms--conditions that emphasized the error. If your sections' membranes contained ion channels (e.g. the pas mechanism), instead of being purely capacitive, and the time step was 0.025 ms, then the trajectories of section1(0.5).v and section2(0.5).v would have differed only slightly from that of section3(0.5).v, and that small difference would have been visible only on close examination of the transient phase of the response to injected current (that is, times earlier than about 3 x membrane time constant).

Solution? When using NMODL-implemented gap junctions, try reducing dt by a factor of 5 or 10 and see if that has a significant qualitative effect on simulation results. If it does, try reducing dt further until you get to the point where you're satisfied that results have converged sufficiently.

Alternatively, you could try adaptive integration, which will automatically choose dt and integration order to satisfy your specified error tolerance. Or, implement gap junctions with the LinearMechanism class, which avoids the problem entirely by including off-diagonal terms. For example, your model implemented with the LinearCircuitBuilder shows identical membrane potential trajectories in the three sections. Of course this last suggestion is inapplicable if the coupled cells are on different hosts.
peterbratby
Posts: 6
Joined: Fri Jul 20, 2018 1:00 am

Re: Spurious voltage changes

Post by peterbratby »

Hi Ted,

Thanks for you reply. This makes a lot of sense.

I developed the example in order to isolate an issue I was having with my full model, which consists of many thousands of gap junction-coupled cells. In the model each cell is identical, and receives the same stimulating current. At first, the cells fire in synchrony (as would be expected), but after a few hundred milliseconds, they become out of phase and start firing erratically. I am pretty sure this is due to the numerical errors you described, even though the issue does not disappear when I reduce the timestep. To confirm this, I plan to implement the gap junctions using the linearmechanism class to see if the issue goes away.

Unfortunately, this will not be a solution for the full model, which will need to be run on multiple nodes. Perhaps I could use linearmechanism only for cells on the same node?

Anyway, I will post on here if I find a satisfactory solution.
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Spurious voltage changes

Post by ted »

At first, the cells fire in synchrony (as would be expected)
For the sake of discussion, forget for the moment that no biological system will have multiple cells that have identical properties, and stay entirely in the world of mathematicians' models. Why is it expected that absolute synchrony would be a stable operating mode of your model? Could initial synchrony be simply an artifact of initialization--an unstable condition, like balancing on a knife edge, from which the cells start to diverge as soon as there is the slightest perturbation (even numerical roundoff error)? Has a dynamical systems analysis of your model proven absolute synchrony to be stable, even for a model with just two of your cells? It would be easy to use the Linear Circuit Builder to set up a model with just two of your model cells coupled by a resistor, and use that to address questions such as:
If they start synchronzed, do they stay synchronized, or is synchrony destroyed by perturbations that result from the finite precision of floating point arithmetic and numerical integration?
If they stay synchronized, how big a perturbation is required to break synchrony, or does it always restore itself?
If they are asynchronous at the start, do they spontaneously become synchronized, and is that synchrony characterized by a finite phase shift between them?
Is there a range of parameters over which the model shows chaotic behavior?
peterbratby
Posts: 6
Joined: Fri Jul 20, 2018 1:00 am

Re: Spurious voltage changes

Post by peterbratby »

Thanks again for your help.

I have implemented the gap junction using linearmechanism, which I understand should be equivalent to using the Linear Circuit Builder.

With just two cells, I am unable to generate any unstable or chaotic behaviour. Even if the cells start unsynchronised, they very quickly (within one cycle) become synchronised. This is true even for very low values of conductance.

When I do the same tests using parallel transfer, I also do not get unstable behaviour. However, if I increase the value of the conductance, I can observe a significant difference to the linearmechanism case. The cells take longer to start spiking, but then spike at a much higher frequency. As I decrease the timestep, the solution converges to the linearmechanism case.

Certainly it seems that a higher value for conductance results in the inaccurate behaviour.

I am continuing to investigate this, using more cells.
peterbratby
Posts: 6
Joined: Fri Jul 20, 2018 1:00 am

Re: Spurious voltage changes

Post by peterbratby »

I have found some more curious results in a gap junction connected network using the linearmechanism class.

The following link shows the results of running a network of 75 identical cells connected by gap junctions implemented by linearmechanism. Initial conditions were set by voltage clamping each neuron for a random length of time. Each set of three traces shows the results of the simulation with identical initial conditions, but decreasing values for the timestep, as indicated by the values shown in the centre. The only difference between the left and right set of simulations is the value of the initial conditions. The black trace shows the membrane potential of one cell, the grey traces the membrane potential of all directly connected cells.

https://drive.google.com/file/d/1FDbxRa ... sp=sharing

As can be seen, the network at first synchronises, but soon becomes unstable. The point at which instability occurs seems to depend on the initial conditions, but not on the length of the timestep. In fact, decreasing the timestep seems to suggest that the simulation is converging on a solution.

I am trying to determine if this behaviour is an authentic result of the simulation, caused by an error in my implementation, or due to some numerical instability in the integration of the differential equations. Below I reproduce the code which implements a gap junction using the linearmechanism class.

Code: Select all

def make_gap_junction_linear(cell1, n1, s1, cell2, n2, s2, g):

    sl = h.SectionList()
    sl.append(sec=cell1.dendrites[n1])
    sl.append(sec=cell2.dendrites[n2])
    xvec = h.Vector([s1, s2])
    cm = h.Matrix(2, 2, 2)   
    gm = h.Matrix(2, 2)
    y = h.Vector(2)
    b = h.Vector(2)
    a1 = h.area(s1, sec=cell1.dendrites[n1])
    a2 = h.area(s2, sec=cell2.dendrites[n2])
    gm.setval(0, 0, g * 100/a1)
    gm.setval(0, 1, -g * 100/a1)
    gm.setval(1, 1, g * 100/a2)
    gm.setval(1, 0, -g * 100/a2)
    lm = h.LinearMechanism(cm, gm, y, b, sl, xvec)
    return sl, xvec, cm, gm, y, b, lm
    
I hope someone can help shed some light on this.
Post Reply