Gap Junctions as Resistor and Capacitor in Parallel

[shailesh]

Hi,

Posting this as a separate query as it was on a different note to the one I posted earlier. I have been trying to extend the NMODL code for implementing a gap junction (as given in the NEURON book) to have both a resistive and capacitive component (i.e. a gap junction represented as a resistor and capacitor in parallel). For starters, I went through ModelDB for gap junctions implemented via NMODL but it turned out they were all purely resisitive. The NEURON forum also failed to return any results on the same and so I decided to have a hand at it.

At the very basic, I believe the significant change (to the gap junction code in the NEURON book) would be change the current calculation equation from:

Code: Select all

i = (v - vgap)/r

to

Code: Select all

i = ((v - vgap)/r) + ( (d(v - vgap)/dt) * Cj)

The second term will basically represent the capacitive current, i.e. ic = c * dv/dt where v=(v-vgap). I hope I am correct till here.

Now I was having a bit of trouble trying to calculate dv/dt in NMODL. Then I came across this post: http://www.neuron.yale.edu/phpbb/viewto ... =16&t=2146. I had my doubts whether I could use that exact same approach of calculating dv/dt = 1000*i_cap/cm. But for lack of a better alternative, I went ahead with it and this what I ended up with:

Code: Select all

NEURON {
	POINT_PROCESS Gap
	POINTER vgap, icap1, icap2		: 1 = same cell, 2 = other cell
	RANGE r, i, dv1dt, dv2dt, cm, Cj	: 1 = same cell, 2 = other cell
	NONSPECIFIC_CURRENT i
}

PARAMETER {
	r = 1e10 (megohm)			: Gap Junctional Coupling Resistance
	Cj = 1000 (picofarad)			: Gap Junctional Coupling Capacitance
}

ASSIGNED {
	v (millivolt)
	vgap (millivolt)
	i (nanoamp)
	icap1 (milliamp/cm2)			: 1 = same cell
	icap2 (milliamp/cm2)			: 2 = other cell
	cm (microfarad/cm2)
	dv1dt (millivolt/ms)			: 1 = same cell
	dv2dt (millivolt/ms)			: 2 = other cell
}

INITIAL {
	dv1dt = (1000)*icap1/cm			: 1 = same cell
	dv2dt = (1000)*icap2/cm			: 2 = other cell
}

BREAKPOINT {
	dv1dt = (1000)*icap1/cm			: 1 = same cell
	dv2dt = (1000)*icap2/cm			: 2 = other cell	
	i = ((v - vgap)/r) + ((dv1dt-dv2dt)*Cj)
}

UNITSON

The hoc code to go along with this is:

Code: Select all

create cella, cellb
access cella

objref g[2]

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

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

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

setpointer g[0].icap1, cella.i_cap(0.9999)
setpointer g[1].icap1, cellb.i_cap(0.0001)

setpointer g[0].icap2, cellb.i_cap(0.0001)
setpointer g[1].icap2, cella.i_cap(0.9999)

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

There were no compile errors etc. But the output was wayward and most certainly incorrect. As expected, setting Cj=0 does make it work properly as a purely resisitive gap junction.
If I were to jot down my questions:
1> Is the basic equation that I stated at the outset correct?
2> Would it be right to implement dv/dt = 1000*i_cap/cm here? I am a bit doubtful as, I believe, i_cap is dependent on only cm and Cj is not coming into the picture anywhere in its calculation.
3> On a broader level, what would be the correct/better method to implement this in NMODL as a point process? Any tips or pointers?

Thanks,
Shailesh Appukuttan

[ted]

The best way to implement a gap junction that involves capacitance is with the LinearMechanism class. This can be done with hoc statements, but the quickest way is with the Linear Circuit Builder (LCB henceforth). If you are not familiar with the LCB, please work through this tutorial
http://www.neuron.yale.edu/neuron/stati ... ncir1.html

Using the LCB to set up a gap junction that involves capacitance

1. Click on the LCB's Arrange radio button. Drag two instances of the "cell" icon (the one that has just a single electrode) onto the Linear Circuit Builder's canvas.

2. Click on the LCB's Label button, then on "Change". Click on one cell icon. This brings up a "Select a Section" tool that displays a shape plot. If the shapes of the two cells overlap, It may be necessary to use 3D Rotate to separate them (click on the square "menu box" in the left upper corner of the canvas and scroll down to 3D Rotate, then release the mouse button and click on the canvas above the shape of the model cells and click on the canvas, then drag down while keeping the mouse button depressed; the blue dot shows where the electrode is attached, and you can change its location by clicking on one of the sections; the section name and "range" of the selected point will be displayed above the shape plot; click on Accept to dismiss this tool).
At this point you should see two cell icons, each showing a different name.

3. Click on Arrange, then drag a resistor and a capacitor onto the canvas and connect one end of each to one electrode, and the other ends to the other electrode. Now you are ready to specify parameter values, initial conditions, etc..
A suggestion: before proceeding to specify parameters etc., it would be a good idea to save the LCB to a session file so you can retrieve it later if necessary. Hint: make an init.hoc file that has the following contents--

Code: Select all

load_file("nrngui.hoc")
load_file("cellspec.hoc") // contains the code that creates your two model cells
load_file("lcb.ses") // session file to which you have saved just the linear circuit builder

Then you will be able to recreate everything by using NEURON to execute init.hoc. And after you have further customized your LCB, just save it to the same ses file so you'll be sure to get the most recent version of it.

[shailesh]

Hi Ted,

Thanks for the quick response. I had actually posted another post (viewtopic.php?f=8&t=2650&p=10559#p10559) around the same time in which I had used the Linear Circuit Builder for implementing a gap junction with a resistor and capacitor in parallel. But that was more with a different intention (to test the Impedance plot). But still thanks for the instructions on using LCB; I used to end up (going roundabout) changing access sections on xterm to be able to place different cells on the LCB, didn't realize that they were hidden from view and could be accessed by 3D rotate. So thats one more thing learnt!

The reason I am really keen on implementing the same via NMODL is:
> I am working on a network of interconnected cells where num(cells)>>>100 and thus coupling via LCB would be too tedious and unmanageable
> I am already having several such models where resistive coupling has been implemented via NMODL (gap.mod). To test the effect of a capacitive component in the coupling for these models, it would have been ideal to just edit gap.mod
> (adapting a line from the NEURON book) I am more of a inveterate programmer than a GUI user :)

It would be a huge help if you can direct me on how I should go about from here, especially in implementing the same via NMODL.

Thanks,
Shailesh Appukuttan

[ted]

The reason I am really keen on implementing the same via NMODL is:
> I am working on a network of interconnected cells where num(cells)>>>100 and thus coupling via LCB would be too tedious and unmanageable

Of course. But NMODL is not the way to implement gap junctions that involve capacitance. To quote my previous post

The best way to implement a gap junction that involves capacitance is with the LinearMechanism class. This can be done with hoc statements

So check out the documentation of LinearMechanism in the Programmer's Reference, then go to ModelDB
https://senselab.med.yale.edu/ModelDB/
and find some working examples by using its Google search to locate model entries that use LinearMechanism.

[shailesh]

I did go through ModelDB and other related articles on the Forum and located 2 hoc implementation of gap junctions via LinearMechanism:
> Olfactory bulb mitral cell: synchronization by gap junctions (Migliore et al 2005)
http://senselab.med.yale.edu/modeldb/Sh ... db\gap.hoc
>CA1 oriens alveus interneurons: signaling properties (Minneci et al. 2007)
http://senselab.med.yale.edu/modeldb/Sh ... nt\gap.hoc

I have been going through them, but facing a bit of difficulty in understanding LinearMechanism implemented via hoc code (the matrices for starters). I would resist asking any doubts at this stage, before I get to better grips with it. But do suggest, if any, starting material / course tutorial etc dealing with it that I could learn from. Any sections from the NEURON book?

One last thing: I am tempted to ask why the NMODL code (from the first post) turns out to be incorrect. I am pretty sure it is, but couldn't figure out why precisely.

Thanks again for your time and effort,
Shailesh Appukuttan

[ted]

I have been going through them, but facing a bit of difficulty in understanding LinearMechanism implemented via hoc code (the matrices for starters).

The documentation is almost obscurely concise. Suggest you make a toy model consisting of two single compartment cells coupled by a gap junction, for which you will know the analytical solution, make your best guess about the matrix syntax, see what happens, and if the result is other than you expected, revise and retry.

But do suggest, if any, starting material / course tutorial etc dealing with it that I could learn from. Any sections from the NEURON book?

Not yet.

One last thing: I am tempted to ask why the NMODL code (from the first post) turns out to be incorrect.

The discretized cable equation's ODEs are not exposed in a way that grants direct access to the derivatives of the state variables. Membrane potential can only be perturbed by affecting charge balance, i.e. by mechanisms that WRITE currents.

[shailesh]

Returning to this discussion after a long hiatus...

I took your advice and decided to model R-C gap junctions using LinearMechanism. Though by now I have some grip on using LinearMechanisms, I decided to make use of the Linear Circuit Builder (LCB) to start off with and take a look at the LinearMechanism code it produces. Here I had some basic queries:

When I have just a single section in NEURON (with no other elements in the LCB except that section), then on clicking "Print Matrix Info", we get the following:

Code: Select all

Matrix equations
LinearCircuitNode[0]
LinearCircuitNode[1]
   ElementBase[1] InsideCell[0] 0
cmat
 -0       -0
 -0       -1
gmat
 1        0
 -0       -0
singular values
1       1
Circuit equations have 0 inconsistencies

Q1. What do the two rows/cols refer to here?
>> I suppose one (2nd row/col in the matrix) would map on to the specified section(location). But, wasn't sure what the first one refers to? Ground possibly? I tried to test that by connecting the section's terminal to ground via a resistor - and there were changes in gmat[1][0] and gmat[1][1], possibly confirming my assumption?!

Q2. Why is cmat[1][1] = -1?
Q3. Why is gmat[0][0] = 1?
>> I was expecting cmat and gmat to be zero in the absence of any other circuit element on the LCB

Q4. What do 'singular values' indicate? Do they correspond to any of the matrices supplied to LinearMechanism?

If we consider the actual scenario of two cells connected via a parallel R-C gap junction (Rj = 10 MOhm, Cj = 5 nF), then I get:

Code: Select all

Matrix equations
LinearCircuitNode[14]
LinearCircuitNode[8]
   ElementBase[3] InsideCell[0] 0
   ElementBase[4] Resistor[0] 0
   ElementBase[7] Capacitor[0] 0
LinearCircuitNode[13]
   ElementBase[7] Capacitor[0] 1
   ElementBase[4] Resistor[0] 1
   ElementBase[2] InsideCell[0] 0
cmat
 -0       -0       -0
 -0       -0.997   -0.00318
 -0       -0.00318 -0.997
gmat
 1        0        0
 -0       6.37e-05 -6.37e-05
 -0       -6.37e-05 6.37e-05
singular values
1       1       0.993506
Circuit equations have 0 inconsistencies

I have confirmed the values of Rj and Cj from the matrices and those specified in LCB:

** Verify value of gmat **
gmat = 6.37e-05 (S/cm2)
1 / (gmat * area(0.5) * 1e-8) = 9994031 (Ohm) = 9.994 (MOhm) = Rj

** Verify value of cmat **
cmat = 0.00318 (1000uF/cm2)
cmat * 1000 * area(0.5) * 1e-8 = 0.0049951323 (uF) = 4.9951323 (nF) = Cj

I am comfortable with the matrix for gmat - the bottom-right sub-matrix referring to the resistive component of the gap junction.
But because of the earlier query (Q2 above) I have trouble understanding cmat - the bottom-right sub-matrix referring to the capacitive component of the gap junction. Why are the diagonal elements (-1 + Cj') - where does the -1 arise from?

[hines]

Though by now I have some grip on using LinearMechanisms, I decided to make use of the Linear Circuit Builder (LCB) to start off with and take a look at the LinearMechanism code it produces.

I am not certain that you are aware of the "Create class" button in the "Parameters" panel of the Linear Circuit builder. That creates a hoc file that wraps the linear circuit in a class that can be
instantiated any number of times. It contains a lot of comments and cmat is correctly printed. If you have questions after studying the hoc file, we can deal with them here.
I see that the 'Print Matrix info' button is printing not cmat but a related matrix that is modified to avoid the false imputation of inconsistency when cell node icons are a part of the circuit.
The first equation is always vground = 0 and that trivial equation and variable is removed when the circuit is passed into the LinearMechanism instance.
The second equation is the contribution of the current balance equation for a node and cmat is really 0 , 0. I need to fix the ouptut of the "Print Matrix Info" button.

If you use the "Create class" button for gap jucntion with c and r, then drag two internal node icons onto the canvas and conect them with a C and an R and label the nodes v1 and v2.
It is not necessary that the nodes be associated with sections.

[hines]

http://www.neuron.yale.edu/hg/neuron/nr ... ef66325387
conatins the fix for the Linear Circuit Builder "Print Matrix Info" button.

[shailesh]

Just after posting the earlier query, I was about to add that interestingly when we 'Create Class' from LCB, the LinearMechanism hoc code consists of only 2x2 matrices as opposed to 3x3 when we 'Print Matrix Info' - and why the difference? - But your post already answers that. I edited my "lincir1.hoc" file based on the additions/deletions listed in the fix you posted. And yes, everything seems clear now. Except for this:

I always (wrongly) thought that 'y' (in the LinearMechanism) could only refer to a sections internal potential or extracellular potentials (as determined by layer = 0/1/2). But your last line:

It is not necessary that the nodes be associated with sections.

got me thinking as to how I could create an arbitrary node in HOC, which I could use as a junction for connections.

Using LCB I created two cells (their internal potentials giving us two nodes) connected via two resistances in series (giving us one more node located between the two resistances), for a total of 3 nodes. I tried to interpret the output of 'Create Class' and this is what I found:

Parameter :: Size
----------------------
_c :: 3 x 3 (Matrix)
_g :: 3 x 3 (Matrix)
_y :: 3 (Vector)
_y0 :: 3 (Vector)
_b :: 3 (Vector)
_sl :: 2 (SectionList)
_xloc :: 2 (Vector)
_layer :: 2 (Vector)

I realized that I was wrong in believing that ALL the above must of same size and the documentation confirmed it:

c and g must be square matrices of the same rank as the y and b vectors.

(and not neccessarily _sl, _xloc and _layer).

But I am unsure as to how NEURON figures out which equations in gmat (_g) and cmat (_c) correspond to nodes determined by the entries in _sl, _xloc, _layer. Also, how can I use another LinearMechanism (if I wanted) to refer to the 'non-section' node that we have above, (e.g. if I plan to use it as a junction from which multiple connections may occur)?

[hines]

It is not necessary that the nodes be associated with sections.

Sorry. Unclear. I only meant that when you are using cell icons, they have to be associated with specifc locations on specific sections before the circuit is instantiated.
However when creating a class, it is not necessary for any sections to exist. In that case the label on all the cell icons will be "no location". But the class does require that, for each of the
cell icons, at least one of the ports has a label to disambiguate each of the LM.<label>_loc(x) methods. (if the vi and vext ports both have labels, the vi label will be used for the LM.*_loc(x) method.

But I am unsure as to how NEURON figures out which equations in gmat (_g) and cmat (_c) correspond to nodes determined by the entries in _sl, _xloc, _layer.

The first _sl.count rows of _c and _g are associated with the section current balance equations and the first _sl.count _y elements are associated with the voltage at those nodes.

[shailesh]

That clears up quite a bit. Just a few more queries on this:

1> Under what conditions do we get the following message on "Print Matrix Info" in LCB:

Circuit equations have X inconsistencies

X: number>=0
Am I correct in believing that it would come-up when we have an incomplete/invalid circuit?

2> But I could not make out why it gave me the message in the following case:
I tried connecting two cells via resistor using LCB without defining any sections. As you said, the labels on the cell icons were "no location". I labeled the nodes as 'Va' and 'Vb'. When I clicked "Print Matrix Info" it printed the matrices and then finished with the line:

Circuit equations have 1 inconsistencies

What is the inconsistency involved here? And interestingly, at this point if I create a section and then change the label for even one of the cell icons on LCB, the inconsistency goes away. What was the cause of the inconsistency and how did it go away just by changing one of the cell icons?

3> The final query is not something that is bothering me at the moment, but more out of curiosity:
Suppose I wanted to model 3 cells (cella, cellb, cellc) connected in the following manner:

Code: Select all

cellA ----- X ----- cellB
            |
            |
          cellC

where lines represent resistive pathways and X is a node (not part of any section).
How would we implement this using LinearMechanism directly or via LCB? Suppose we used LCB and the labels were cella.v(0.5): Va, cellb.v(0.5): Vb, cellc.v(0.5): Vc and V@X = Vx. If I wanted to split the connections into two parts (right now just to bring out my point more than anything else) whereby, I want to create one class for connecting: cellA ----- X ----- cellB and a second class for creating: X ----- cellC.

Class 1 would give me two methods: lm.Va_loc(x) and lm.Vb_loc(x)... and nothing for node X.
Class 2 would give me one method: lm.Vc_loc(x)... and again nothing for node X.

In such a scenario, how can I covey to LM/LCB that the node X I am referring to in both the cases is the same common node? For section nodes I am aware that we make use of _sl, _xloc and _layer to specify the appropriate node. But what can we do for non-section nodes?

I realize that the above problem can be (I suppose) transformed via star to delta circuit conversion... but as I said this is just a hypothetical problem made-up for the above scenario. Not sure whether I have complicated a simple issue in posing my question.

[hines]

1) The LCB matrix is not invertable. The number is the number of zeros of the singular value decomposition. The simple reasons are: no ground, two batteries in parallel, two current clamps in
series, or any other situation in which two current balance equations end up identical. In other words N variables but the N equations are not linearly independent. During Gaussian elimination one of
the rows will end up as 0.

2) Until the location is specified, the cell port may have not path to ground. If you give it a ground, such as connecting a battery or resistor, and then connection the other end to ground, there will
be no inconsistency in terms of "Print Matrix Info", (Surprise to me, it even seems works during a simulation But I wouldn't trust it without further testing).
Note that you can create hoc classes even of inconsistent circuits. eg a resistor with nothing connected to either node.

3) distinct instances of (same or distinct) classes cannot share internal nodes. x has to become a cell location (e.g. a single compartment with cm=0 and no mechanisms) so that it can be used
as a sharable connection point (sum of currents from all edges into the node = 0.)

[shailesh]

Thanks a lot! That takes care of all the queries I had at the moment.