Mismatch between NEURON results and Circuit Output?

Anything that doesn't fit elsewhere.
Post Reply
shailesh
Posts: 104
Joined: Thu Mar 10, 2011 12:11 am

Mismatch between NEURON results and Circuit Output?

Post by shailesh »

Hi,

I have been trying to develop a network model involving extracellular mechanism as well as constructs using LinearMechanism. Since neither of these are intuitive (and prone to user errors), I made a small dummy model just to test my understanding of what was happening (with the LinearMechanisms omitted to start with). I then implemented the electrical equivalent circuit on a circuit simulation platform (here: MultiSim) to compare the results (to be sure what I coded was what I meant... could possibly also have used the cell-builder, but wanted to explicitly define all circuit elements to be certain). The results didn't seem to provide a good match in the presence of extracellular mechanism. Here is a summary of what I did:

Model: Two cells, cellA & cellB, with nseg = 1 connected via gap junctions. cellA is stimulated via IClamp.

Code for developing model:

Code: Select all

create cella, cellb

forall {	
	L = 20000
	diam = 20/PI	//For rounded off values of area()

	Ra = 181
	cm = 1

	insert pas
	g_pas = 1 / 132500
	e_pas = -50

	nseg = 1	
}

//Setting up gap junctions between the two cells
objref g[2]

for i = 0 ,1 {
	g[i] = new Gap(0.5)
	g[i].r = 10	//10 MOhm
}

cella g[0].loc(0.9999)
cellb g[1].loc(0.0001)

setpointer g[0].vgap, cellb.v(0.0001)
setpointer g[1].vgap, cella.v(0.9999)


//Function to incorporate extracellular mechanism in both cells with identical properties - Invoked from GUI panel
proc setExtra() {
	forall {
		insert extracellular
		xc[0] = 0
		xc[1] = 0
		xg[0] = 1e-3
		xg[1] = 1e9		
	}
}


access cella

objref stim
cella stim = new IClamp(0.5)

stim.del = 5
stim.amp = 1
stim.dur = 4995

tstop = 5000
Electrical Equivalent Circuit (in presence of extracellular mechanism... as that is when the problem appears):
Image
Circuit parameters were converted into absolute units by multiplying with area(0.5)*1e-8 wherever required, i.e. cm, g_pas, xg[0], xg[1]
Here XMM represents voltmeters. XMM1 = cellA.v(0.5), XMM2 = cellA.vext[0](0.5) , XMM3 = cellB.v(0.5), XMM4 = cellB.vext[0](0.5)

Results:

Extracellular fields: ABSENT
-----------------------------
NEURON Output: cellA.v(0.5) = -31.2654 mV, cellB.v(0.5) = -35.6096 mV
Circuit Output: cellA.v(0.5) = -31.265 mV, cellB.v(0.5) = -35.609 mV
=> Perfect Match

Extracellular fields: PRESENT
-----------------------------
NEURON Output: cellA.v(0.5) = -31.6764 mV, cellA.vext[0](0.5) = 138.291 uV, cellB.v(0.5) = -36.0065 mV, cellB.vext[0](0.5) = 105.611 uV
Circuit Output: cellA.v(0.5) = -31.279 mV, cellA.vext[0](0.5) = 141.291 uV, cellB.v(0.5) = -35.595 mV, cellB.vext[0](0.5) = 108.718 uV
=> Differences in the two cases

I did verify the above simulations on another circuit simulation platform (CircuitLab) but ended up with a similar mismatch with NEURON in presence of the extracellular field. I have read that slight differences may arise because of differences in the integration method. But considering that there were near perfect matches in the absence of extracellular mechanism, I am not sure what is causing the larger variation in its presence.

The simulations however match when considering a single cell in the presence/absence of extracellular fields:

Single Cell with Extracellular field: ABSENT
------------------------------------------------
NEURON Output: cellA.v(0.5) = -16.875 mV
Circuit Output: cellA.v(0.5) = -16.876 mV
=> Perfect Match

Single Cell with Extracellular field: PRESENT
------------------------------------------------
NEURON Output: cellA.v(0.5) = -16.875 mV, cellA.vext[0](0.5) = 250 uV
Circuit Output: cellA.v(0.5) = -16.876 mV, cellA.vext[0](0.5) = 249.996 uV
=> Perfect Match

So we basically have scenarios where:
> Gap junctions working fine in the absence of extracellular fields
> Extracellular fields working fine for single cells
> Problems apparently when both gap junctions and extracellular fields are present

One other query: With nseg = 1 and cells connected only via gap junctions (i.e. no "connect" statement), is it correct to consider axial resistivity, both internal as well as external, (i.e. Ra and xraxial) being irrelevant with regards to the electrical equivalent circuit of the cell? Simulations seemed to suggest so, but wanted to confirm.

Thanks in advance for any inputs!
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: Mismatch between NEURON results and Circuit Output?

Post by hines »

Excellent test! It cannot be sufficiently emphasized how important it is to compare simulation results using a simple model with known solution. In this case:
oc>r1 = 33.125 + 0.25
oc>r2 = 10 + r1
oc>r = r2*r1/(r1 + r2)
oc>r-50
-31.138233
or from your circuit output
oc>-31.279 + 141.291*1e-3
-31.137709
The problem was the hoc code statements

Code: Select all

setpointer g[0].vgap, cellb.v(0.0001)
setpointer g[1].vgap, cella.v(0.9999)
where vgap is supposed to be the internal potential but sec.v... is the membrane potential.
I must admit this did not occur to me before your test focused on the discrepancy. The solution is to introduce vi in the gap.mod file so that

Code: Select all

 BREAKPOINT {
+	vi = v : for use by other side of gap (not vm)
 	i = (vgap-v)/r
 }
and then use

Code: Select all

setpointer g[0].vgap, g[1].vi
setpointer g[1].vgap, g[0].vi
With these changes, my simulation results for gap + extracellular are
oc>cella.v(0.5)
-31.27952
oc>cellb.v(0.5)
-35.59548
oc>cella.vext[0](0.5)
0.14128664
oc>cellb.vext[0](0.5)
0.10871336

Remark; One could equivalently compute vi in a BEFORE BREAKPOINT block but in this case there would be no numerical difference.
shailesh
Posts: 104
Joined: Thu Mar 10, 2011 12:11 am

Re: Mismatch between NEURON results and Circuit Output?

Post by shailesh »

Thanks Michael! That did it. It did not occur to me to think about how the extracellular mechanism could have been interfering with the regular gap junction implementation (felt it was more likely to do with an incorrect circuit implemenation or some worrying number crunching). Glad that you were able to direct me towards course correction before I was more at sea with my larger models :)
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: Mismatch between NEURON results and Circuit Output?

Post by hines »

At a certain point, "larger models', is going to mean parallel cluster (or multithread) simulations which will require a change from the POINTER style for vgap to
the parallel transfer style
http://www.neuron.yale.edu/neuron/stati ... source_var
Unfortunately, 11 months ago, there was a change
http://www.neuron.yale.edu/hg/neuron/nr ... ac2b0cef43
having to do with cache efficient recalculation of source_var pointers that restricted the source_var to being a segment membrane potential.
I will look into the possibility of checking each of the relevant segments to see if extracellular is inserted, and, if so, make sure the transfer is
v + vext[0]
shailesh
Posts: 104
Joined: Thu Mar 10, 2011 12:11 am

Re: Mismatch between NEURON results and Circuit Output?

Post by shailesh »

Yes I can see that would become an issue in parallel implementation of gap junctions (in presence of extracellular mechanism). Is a workaround possible in the current NEURON framework for parallel simulations? I haven't done any serious work with Parallel sims yet and so don't have much clue. I did read that I would have to move on to implementing gap junctions using Parallel Transfer method but haven't got down to doing that yet (maybe in a couple of weeks time). Thanks for the heads up.
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: Mismatch between NEURON results and Circuit Output?

Post by hines »

http://www.neuron.yale.edu/hg/neuron/nr ... 4dac40175f
handles the transfer of v + vext now for
section pc.source_var(&v(x), sid)
The correct currently accessed section is now required. ie.
pc.source_var(&section.v(x), sid)
would probably return an error (even if there is no extracellular involved).
shailesh
Posts: 104
Joined: Thu Mar 10, 2011 12:11 am

Re: Mismatch between NEURON results and Circuit Output?

Post by shailesh »

Thanks Michael, that was real quick! I was eager to test it out but having some trouble in compiling the latest code from Mercurial. I believe the query would be more apt under a different head and hence have posted it here:
viewtopic.php?f=5&t=3094
Would be great if you could help me out!
shailesh
Posts: 104
Joined: Thu Mar 10, 2011 12:11 am

Re: Mismatch between NEURON results and Circuit Output?

Post by shailesh »

Remark; One could equivalently compute vi in a BEFORE BREAKPOINT block but in this case there would be no numerical difference.
I was wondering about this. I tested it on two cells (nseg = 1) connected via gap junctions for a simulation duration of 0.1 ms.

vi in BEFORE BREAKPOINT:

Code: Select all

oc>cella.v(0.5)
        -64.963742
oc>cellb.v(0.5)
        -64.988668
oc>g[0].i
        -0.0018707837
oc>g[1].i
        0.0018707837
vi in BREAKPOINT:

Code: Select all

oc>cella.v(0.5)
        -64.963748
oc>cellb.v(0.5)
        -64.988668
oc>g[0].i
        -0.002153892
oc>g[1].i
        0.0018704306
So there was a small difference. cella appears to always read the old (one time step back) value of cellb.v(0.5), whereas cellb always got the current value of cella.v(0.5). I suppose this is down to the order of execution:
cella: BEFORE BREAKPOINT
cellb: BEFORE BREAKPOINT
cella: BREAKPOINT
cella: BREAKPOINT
cellb: BREAKPOINT
cellb: BREAKPOINT

Considering this would it be more appropriate to have vi computed in BEFORE BREAKPOINT? Or am I doing something wrong in testing this?

Also, given a choice where it is equivalent, is it better to move as many statements from BREAKPOINT block to BEFORE BREAKPOINT block for faster simulations... considering the fact that BREAKPOINT would unnecessarily execute all those statements twice. But I suppose for simple assignment statements it won't make a noticeable difference.
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: Mismatch between NEURON results and Circuit Output?

Post by hines »

Considering this would it be more appropriate to have vi computed in BEFORE BREAKPOINT?
Yes. It is not good for results to depend on the order of calculation (at least when the order is not under control)
Also the situation is more deterministic, or perhaps I should say isotropic, when one imagines a two or three dimensional array of cells so that each has 4 or 6 gap junctions where the neighbors are all pointing to the same vi at the center.
Note that this problem will not occur with the ParallelContext.source_var approach where the results will be the same as the BEFORE BREAKPOINT with POINTER style.
Also, given a choice where it is equivalent, is it better to move as many statements from BREAKPOINT block to BEFORE BREAKPOINT block for faster simulations
The only downside I would see would be some loss of clarity. It would makes the ordering issue between communicating mod files even more difficult to analyze.
A minor issue is that things like BEFORE and AFTER tend to get ported ot new architectures, e.g. GPU, last.
shailesh
Posts: 104
Joined: Thu Mar 10, 2011 12:11 am

Re: Mismatch between NEURON results and Circuit Output?

Post by shailesh »

Thanks that clears everything for now. Maybe more (hopefully not) when I start shifting to parallel implementation of gap junctions!
Post Reply