single resistor as gap junction and connecting cells

[pzhang]

What's the best way to model gap junctions?

I'm trying to model a fly heart, I have the single fly heart cells build using channel builder and inserted new mechs, and I want to couple them together with gap junctions. Right now I have done this by creating 20 inidividual cells with 20 cell builders and then using linear circuit builder, connecting them via a single resistor between each cells, in a line.

However, now I wanted to change some parameters for each of the cells, diam and L, and I was hoping to use hoc and a for-loop.

ii=1
for (ii=1; ii<21; ii+=1) access soma_ii
{diam=8 L=50 gnabar_Inaf=0.00078 gCaL_ICaL=0.0005
gKr_IKr= 0.000588 gKs_IKs=0.001 gto_Ito=0.0006}

all my cells are named soma_1, soma_2, and etc, I don't know how I can access each different cells as the above code won't let me do that, as it doesn't understand this soma_ii deal.

Is this the best way to model a gap junction in this case?

I stimulate the first cell using a Iclamp, and it'll generate an action potential, but that seems to die off after 5 or 6 cells, is this a problem with the single cell or is it because of the resistors. I was thinking if my resistor value is not big enough, after the current is delivered to about the 3rd cells, it would leak and not have enough to get it to the threshold value....

Thanks in advance!

[ted]

The GUI is great for rapid development, for managing model cells that have complex
anatomy and/or nonuniform biophysical properties, and for prototyping network models
(using the Network Builder to make small nets that can be mined for code that is
reused to build large network models, as described in chapter 11 of The NEURON
Book). Unfortunately the Network Builder does not yet offer a way to deal with gap
junctions.

For your particular use, it would be best to adopt a strategy similar to that which was
employed by
M. Migliore, M. Hines and GM Shepherd
The role of distal dendritic gap junctions in synchronization
of mitral cell axonal output, J. Comput. Neurosci., 18:151-161 (2005).
Source code is available from ModelDB
http://senselab.med.yale.edu/senselab/m ... odel=43039

Their model involved only two cells, but the approach should scale well. First, they
defined a cell class, which makes it easy to manage common properties and spawn
new instances of cells. You can also alter the properties of individual cell instances.
Second, they defined a gap junction class, which uses the LinearMechanism class to
set up connections between cells, but insulates you from the messy details of fiddling
with the LinearMechanism class at the hoc level.

Let me know if you have any other questions about this, or any help in adapting their
approach to your needs.

ii=1
for (ii=1; ii<21; ii+=1) access soma_ii
{diam=8 L=50 gnabar_Inaf=0.00078 gCaL_ICaL=0.0005
gKr_IKr= 0.000588 gKs_IKs=0.001 gto_Ito=0.0006}

all my cells are named soma_1, soma_2, and etc, I don't know how I can access each different cells as the above code won't let me do that, as it doesn't understand this soma_ii deal.

Two comments here.
1. There is never a reason to use more than one access statement. Use access just
once, to tell hoc which section in your model is "most important" (in your case, maybe
the one that is stimulated first). If there are multiple access statements, an unqualified
variable or parameter name, such as v(0.5), will refer to a property of a completely
different part of the cell depending on the most recently executed access statement.
hoc offers two perfectly good methods for specifying which section you mean: dot
notation i.e. sectionname.attribute, and the "temporary section stack" forms
sectionname statement
and
sectionname { stmt1 stmt2 . . . }
2. The "soma_ii thing"? The for loop doesn't iterate over your cells because soma_ii is
not the name of a cell--you called them soma_1 etc.. The for loop doesn't affect the ii
in the string soma_ii. A slightly awkward workaround would be to use sprint to build
strings that contain the proper names, then execute those strings.

But that's completely unnecessary if you use object oriented programming to manage
your cells. If your cells are instances of a particular class, you would create them with
syntax something like this (assuming the name of the class is FlyHeartCell)
NUMCELLS = 20
objref mc[NUMCELLS]
for ii=0,NUMCELLS-1 mc[ii] = new FlyHeartCell()
Notice that the ii in mc[ii] is a variable, not part of a string, so its value can be changed
by the for loop.

[pzhang]

The model you recommended looks really nice, but I'm having trouble understanding it. I couldn't access model view and I'm not sure how does the gap junction is inserted.
if I just use resister to connect each individual cells, it seems to create a ok connection between the cells, is this theoretically incorrect to model a gap junction? Or is it okay? because after connecting all 20 cells together with resistors, I find it works okay.

[ted]

if I just use resister to connect each individual cells, it seems to create a ok connection between the cells, is this theoretically incorrect to model a gap junction? Or is it okay? because after connecting all 20 cells together with resistors, I find it works okay.

It may work but you have to admit that it's cumbersome, and that makes it prone to user
errors. I am a strong advocate of using the GUI as much as possible, but only when it
makes things easier, not harder.

The code used by Migliore et al. isn't heavily commented, and I can see how it might be
difficult to figure out. But you don't need all of their code--you just need to know how to
use their Gap class in your own models. This takes a bit of detective work.

The essential clue is this bit in their forfig6-modeldb.hoc file--

Code: Select all

rel=0.95

gap[0] = new Gap()
mt[0].tuftden gap[0].src(rel)
mt[1].tuftden gap[0].target(rel)
gap[0].g(ggap)

which looks like they're connecting the 0.95 location on the tuftden section of mt[0] to
the 0.95 location on the tuftden section of mt[1]. Next we look at gap.hoc and discover,
eventually, that it has a public member called g, which is actually a func g(). Furthermore,
func g() calls proc set_gm(), the first line of which is

Code: Select all

proc set_gm() { local us, a// conductance in nanosiemens

so we guess that gap[0].g(ggap) will set the gap junction conductance in nanosiemens
to the numeric value of ggap.

So here's how I reused the Gap class to make a model of a bunch of cells connected
by gap junctions. First I wrote an init.hoc file with these contents:

Code: Select all

load_file("nrngui.hoc")
load_file("mc.hoc")  // defines MC cell class
load_file("gap.hoc")  // defines Gap class
  // for managing gap junctions implemented with the LinearMechanism class
load_file("modelspec.hoc")  // assembles a model of MCs connected by Gaps
load_file("rig.ses")  // user interface ("experimental rig")

Then I used a CellBuilder to specify the properties of my cells (in this case, a single
compartment with hh). When I was finised specifying these properties, I saved the
CellBuilder to a ses file (of course; I called it mc.ses), and then I made the CellBuilder
write a hoc file containing a template that defines the cell class--here's how:
Go to the CellBuilder's Management page.
Click on its Cell Type button.
Click on the Classname button and enter the name of the cell class you want to create
(my choice was MC).
Click on the "Save hoc code in file" button and specify the name of the hoc file to write
(my choice was mc.hoc). By the way, on my PC, I have to pull down the lower edge of
the CellBuilder in order to see this button.

Now I needed to write the hoc code that creates the instances of cells and gap junctions,
and assemble the whole thing into a model. Here's the hoc file, which I called
modelspec.hoc :

Code: Select all

NUMCELLS = 20
objref mc[NUMCELLS], gap[NUMCELLS-1]

// create the cells
for j = 0,NUMCELLS-1 mc[j] = new MC()
access mc[0].soma  // make sure there's a default section
// now we have 20 mc

// arrange them in space for convenience when using the GUI
YOFFSET = 50  // ok for 20 cells
for j = 0,NUMCELLS-1 {
  // y position jumps back to 0 after every 5th cell
  mc[j].position(j*mc[j].soma.L, -(j%5)*YOFFSET, 0)
}

// set up the gap junctions between them
// in this model, the cells and gap junctions are connected in line
// mc[0] gap[0] mc[1] gap[1] mc[2] . . . mc[NUMCELLS-2] gap[NUMCELLS-2] mc[NUMCELLS-1]
// for each adjacent pair mc[j], mc[j+1], gap[j] 
//   attaches mc[j].soma(1-epsilon) to mc[j].soma(epsilon)
//   where epsilon is a small positive number
// rationale:  Gap class returns "divide by 0" error if node area is 0

srcloc = 1
tgtloc = 0

epsilon = 1e-3
func cliploc() {
  if ($1<=0) {
    $1 = epsilon
  } else {
    if ($1>=1) {
      $1 = $1-epsilon
    }
  }
  return $1
}

for j = 0,NUMCELLS-2 {
  gap[j] = new Gap()
  srcloc = cliploc(srcloc)
  mc[j].soma gap[j].src(srcloc)
  tgtloc = cliploc(tgtloc)
  mc[j+1].soma gap[j].target(tgtloc)
}

// specify coupling
proc setg() { local j
  if ($1<0) $1=0
  for j=0,NUMCELLS-2 gap[j].g($1)
}

GAPG = 10 // nS
setg(GAPG)

Be sure to let me know if you have any questions about this.
The only tricky bit here is func cliploc(), which forces the site of the gap to an internal
node (i.e. anyplace except for the 0 and 1 ends of a section). That's because gap.hoc
contains a couple of statements (in proc set_gm()) which involve division by the area associated with the location where the gap junction is attached. At 0 and 1, the area()
function returns 0, causing a "divide by zero" error.

All that remains is to create an experimental rig for stimulating and recording from the
model. I used
a RunControl
a Voltage axis graph showing v(0.5) (belonging to mc[0].soma, the default section),
mc[1].soma.v(0.5), and mc[2].soma.v(0.5)
a PointProcessManager set up as an IClamp at MC[0].soma(0.5), delivering 10 nA
for 0.1 ms starting at 1 ms.

I'll send you a zip file via email that contains the whole thing.

One last item: I ran tests with a simpler model (just two passive cells) to verify that the
conductance needed by Gap's g() method is in nanosiemens, and that the arguments
to its src() and target() methods are indeed the relative locations along the two sections
that are being connected.

[pzhang]

Thank you very much!
I will try this and let you know how it works out.

[Newberry]

Greetings, I felt that in the interest of forum organization that would post in this thread as opposed to creating my own.

I am also creating a network model using the same type of gap junctions as described with a few differences:
1) A much simpler model; 2 cells each composed of a soma and a dendrite connect via gap junction
2) Instead of using 1 cell class, I am using two.

I built the cell classes using the GUI.

I followed the above instructions and made the following init.hoc file

Code: Select all

load_file("nrngui.hoc")
load_file("A.hoc") // defines A class
load_file("B.hoc") // defines B Class
load_file("gap.hoc") // defines Gap class
load_file("modelspec.hoc")
load_file("rig.ses") // user interface (experimental rig)

I do not have rig.ses - but this is the least of my troubles at the moment.

Class A and B have the following appearance in the GUI (ASCII), Where O =soma, and ----- = dendrite

O-----

Thus ideally I want the network to look as follows;

O-----GAP-----O

The reason I have two classes despite the identical appearance is that I want to be able to alter the individual properties of the soma/dendrites to look at the effect on the network as a whole.

What I'm having difficulty doing is connecting these two cells. I am not too terribly experienced with NEURON programming as of yet.

I feel that the "modelspec.hoc" could easily be modified to encapsulate the network I am building.

Thanks in advance.

[Newberry]

Over the weekend I've tried connecting the cells by simply adding the following code to a test .hoc file after i've defined Cell Class A, B and used the template for gap junction as defined by Migliore et al.

Code: Select all

objref gap[0]
gap[0] = new Gap()
A[0].dend gap[0].src(1)
B[0].dend gap[0].target(1)
gap[0].g(ggap)

***Note ggap is defined in the template***

I believe I'm on the right track here because I don't need a complex template as my network is only composed of 2 cells.

Right now, presuming I'm on the right track, I am faced with the following error

nrniv: subscript < 1 gap
in C:/nrn59/Gapjunct/test.hoc near line 245
objref gap[0]
^

Is the code I am attempting to use viable for connecting the cells and how do I fix this error?

[ted]

objref foo
is sufficient to create a single objref.
objref foo[n]
where n > 0 tells hoc to create several objrefs whose names are
foo[0], foo[1] . . . foo[n-1]
So what do you think hoc would do if you tell it
objref foo[0]

[Newberry]

Thanks for the guidance there ted.

I had to do a little bit more tinkering to the above code to get it to work properly. It had quite a few errors in it. I'll be doing some testing to see if I actually got it functioning correctly.

Thanks again.

[Newberry]

In case anyone is interested:

Init.hoc

Code: Select all

load_file("nrngui.hoc")
load_file("A.hoc") // defines A class
load_file("B.hoc") // defines B Class
load_file("gap.hoc") // defines Gap class
load_file("modelspec.hoc") // connects A to B via gap junction

Where gap.hoc is from Migliore et al.

modelspec.hoc

Code: Select all

objref gap
gap = new Gap()
Adend gap.src(0.95)
Bdend gap.target(0.95)
gap.g(10)

[ted]

Thanks for the followup on that. I'm sure there will be others who want to use the Gap
class, and they will benefit from your experience.

[shailesh]

In context to one of the earlier posts on this chain dealing with:
M. Migliore, M. Hines and GM Shepherd
The role of distal dendritic gap junctions in synchronization of mitral cell axonal output, J. Comput. Neurosci., 18:151-161 (2005)
File: http://senselab.med.yale.edu/modeldb/sh ... db\gap.hoc

One last item: I ran tests with a simpler model (just two passive cells) to verify that the
conductance needed by Gap's g() method is in nanosiemens, and that the arguments
to its src() and target() methods are indeed the relative locations along the two sections
that are being connected.

I was having trouble working the dimensions out here. If conductance is in nanosiemens, then:
func g() -> sets g_ in nanosiemens
Then,
proc set_gm() -> sets matrix gmat elements as follows:
us = .001*g_ //not sure why the variable is named 'us' giving the impression it is in microSiemens(uS); but doesn't *.001 make it picoSiemens(pS)?
a = 100/area(xvec.x[0])
gm.x[0][0] = us*a

i.e. gm.x[0][0] = .001*g_ * 100/area(xvec.x[0])

If we work out the dimensions, we get:
gm.x[0][0] = .001*nanoSiemens * 100/(square microns)
-as area() returns the area in square microns (documentation)
So,
gm.x[0][0] = .001*10e-9*Siemens * 100/(10e-8 square cm)
gm.x[0][0] = 0.01 * S/cm2

From the documentation for 'LinearMechanism', I read that:

Note that current balance equations of sections when 0 < x < 1 have dimensions of milliamp/cm2 and positive terms are outward. Thus... g elements have dimensions of S/cm2

So, isn't there a problem here? Isn't the conductance matrix having incorrect units?

[ted]

First a short summary of units in NEURON.
time is in ms
voltage in mV
conductance and current in S/cm2 and mA/cm2 for density mechanisms
but uS and nA for "point processes" ("point sources") such as synapses and gap junctions.

The authors wanted to specify gap junction conductances in nS, but NEURON needs point processes to use uS, so what scale factor should be applied?
We know that
1000 nS = 1 uS
so
1 = 0.001 uS/nS
This means that, if we have a parameter with numerical value p and units of nS,
the way to convert to units of uS is to multiply p by 0.001 uS/nS.
i.e.
p nS = 0.001*p uS
And that's why func set_gm() contains the assignment statement
us = .001*g
(I strongly advise against writing fractional numbers between -1 and 1 without a 0 in front of the decimal point--it's way too easy not to see that little dot, with potentially dire consequences)

From the documentation for 'LinearMechanism', I read that:

Note that current balance equations of sections when 0 < x < 1 have dimensions of milliamp/cm2 and positive terms are outward. Thus... g elements have dimensions of S/cm2

True. It just happens that conductance density in S/cm2 and total conductance in uS are numerically identical when surface area is 100 um2*, and that's why the factor a is calculated as 100/segment area in um2.

*--To see why, suppose you have a cell with 100 um2 surface area and total conductance of p uS. The assertion is that this implies a conductance density of p S/cm2. Here's how to show that this assertion is correct.
Clearly the conductance density is
0.01*p uS/um2
We need to convert this to units of S/cm2.
First let's convert um2 to cm2.
1 cm2 = 1e8 um2
so
1 = 1e8 um2/cm2
and we have
0.01*p uS/um2 * 1 = 0.01*p uS/um2 * 1e8 um2/cm2
= 1e6*p uS/cm2
Now let's convert uS to S.
1e6 uS = 1 S
so
1 = 1e-6 S/uS
and we have
1e6*p uS/cm2 * 1e-6 S/uS
= p S/cm2

[shailesh]

Thanks Ted for clarifying that (I actually feel a bit stupid for not having figured out the first part); but the second part was something I had not given thought to. Makes sense now. Further to this, I had two more queries dealing with units (within NEURON). I had posted it elsewhere (viewtopic.php?f=8&t=2654&p=11485#p11485) but pasting the relevant part here (you could take your time in going through the rest of that post - its long!):

>> I read that 'ge' needs to be provided in terms of 'S/cm2' for LinearMechanism when intended for non-terminal nodes. So if I need the resistive link to be (say) 5 MOhm = 200 nS, then is it correct to set ge = (200e-9)/(area(0.9999)*1e-8) = 200e-1/area(0.9999)

>> xraxial is required in units of MOhm/cm. Again, if (say) I need a 5 MOhm resistance between vext[0] nodes of two adjacent segments, then is it correct to set xraxial = 5/(0.0001*L/nseg)

Thanks in advance!

[ted]

Thanks Ted for clarifying that (I actually feel a bit stupid for not having figured out the first part)

Reconciling different units is endlessly confusing. It seems more a topic for accountants and bookkeepers than for mathematicians or physicists. The challenge of deriving conversion factors, especially on the spot in front of an audience, can induce mind-numbing tedium or sheer panic.

I read that 'ge' needs to be provided in terms of 'S/cm2' for LinearMechanism when intended for non-terminal nodes. So if I need the resistive link to be (say) 5 MOhm = 200 nS, then is it correct to set ge = (200e-9)/(area(0.9999)*1e-8) = 200e-1/area(0.9999)

That's exactly as I would do it on paper. My thought process without paper would have been something like this:
200 nS = 0.2 uS = 0.2 S/cm2 on 100 um2
If numeric value of area in um2 is some multiple of 100, conductance should be proportionally smaller, therefore
conductance density = (0.2 S/cm2)*(100/area) = (20/area) S/cm2

xraxial is required in units of MOhm/cm. Again, if (say) I need a 5 MOhm resistance between vext[0] nodes of two adjacent segments, then is it correct to set xraxial = 5/(0.0001*L/nseg)

Right again. You know R in megohms, and the distance is L/nseg but that's in um, so
R/(L/nseg) needs to be multiplied by (1e4 um / cm).