Extracellular Resistance as a 3D Grid - Bidomain Model

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

Extracellular Resistance as a 3D Grid - Bidomain Model

Post by shailesh » Sun Jul 15, 2012 11:08 am

Hi,

I assume it best to explain my objective first so that you could suggest a better approach to tackle my problem. I am attempting to develop a bidomain model (one grid representing cells, coupled via gap junctions, and the other grid their adjacent extracellular space). This would be a passive system and 3-dimensional with num(cells)>>>100. I have developed the first grid with gap junctions implemented as point processes via NMODL and it seems to be working fine. My task is now to build the second grid.

To start off, I am trying to keep it simple and just develop a 2 cell model (by a method which I can easily and efficiently scale for my larger 3D models). From what I have read on the forum, it means I should rule out using LinearMechanism. Do correct me if I am wrong here.

The most natural approach seemed to be to utilize 'extracellular' mechanism in NEURON. I would by default obtain the connection between the two grids via the cell membrane. I believe, I would only need one layer of extracellular, i.e. nlayer=1. The documentation stated that:
If other than 2 extracellular layers is desired, you may recompile the program by changing the nrn/src/nrnoc/options.h line #define EXTRACELLULAR 2 to the number of layers desired. Be sure to recompile both nrnoc and nrniv as well as any user defined .mod files that use the ELECTRODE_CURRENT statement.
I made the change to options.h file. But I am not clear on what was meant by "recompile both nrnoc and nrniv". I assumed it meant to run mknrndll on 'nrn71\src\nrnoc' folder. I found 'nrniv' in 'nrn71\bin' and ran mknrndll on that as well (I always assumed mknrndll was only for mod files, does it also operate on dll files?). But even then accessing Tools>Distributed Mechanisms>Viewers>Shape Name>(selecting section) I could still see parameters for the 2nd extracellular layer. How can I verify whether nlayer = 1 or 2?

The next concern was to have the extracellular space of adjacent cells connected/coupled (the ECF at last segment of one cell connected to the ECF of first segment at the other cell). I believe, by default, the ECFs for different sections in NEURON are disconnected. I searched through the forum for any related articles, and one of the likely ones seemed to be regarding ephaptic coupling (viewtopic.php?f=12&t=1812&p=6503&hilit=ephaptic#p6503). But it seemed to suggest using LinearMechanism which again seems to be infeasible for the larger model here.

A seemingly weird thought that came up, considering that I had already coupled the cells via gap junctions, was to have a similar mechanism for coupling/connecting the extracellular resistances. Though I had my own doubts, I gave it a try. I made another NMODL file (pretty much identical to gap.mod):

Code: Select all

NEURON {
	POINT_PROCESS GapE
	POINTER vext,vgapext
	RANGE r, i
	NONSPECIFIC_CURRENT i
}

PARAMETER {
	r = 1e10 (megohm)
}

ASSIGNED {
	vext (millivolt)
	vgapext (millivolt)
	i (nanoamp)
}

BREAKPOINT {
	i = (vext - vgapext)/r
}
The relevant hoc code being:

Code: Select all

//Coupling Extracellular
objref ge[2]
for i = 0 ,1 {
	ge[i] = new GapE(0.5)
	ge[i].r = 10
}
cella ge[0].loc(0.9999)
cellb ge[1].loc(0.0001)
setpointer ge[0].vext, cella.vext[0](0.9999)
setpointer ge[1].vext, cellb.vext[0](0.0001)
setpointer ge[0].vgapext, cellb.vext[0](0.0001)
setpointer ge[1].vgapext, cella.vext[0](0.9999)
But I realized (should have thought of it in the first place) that the NONSPECIFIC_CURRENT i would only be injected in the corresponding intracellular space and not at the extracellular nodes. Is there any way to overcome the same?

Another tame effort was to change the BREAKPOINT block to:

Code: Select all

BREAKPOINT {
	vext = vext - ((vext - vgapext)/r)
}
But obviously this was also doomed to fail as vext would just get overwritten with values calculated here vs those for extracellular mechanism rather than each contributing to the final value. I am overdoing it by this stage, but can setting e_extracellular of first segment of cell2 = vext[0] of last segment of cell1 (and vice versa) offer any possibility? Apologies in advance if it doesn't make sense, I think I need to take my mind off it for a while.

I realize I am thinking along a single track, so do advice me on what would be a better approach to develop the extracellular grid of the bidomain model.

Thanks,
Shailesh Appukuttan

ted
Site Admin
Posts: 5592
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Extracellular Resistance as a 3D Grid - Bidomain Model

Post by ted » Sun Jul 15, 2012 1:02 pm

shailesh wrote:To start off, I am trying to keep it simple and just develop a 2 cell model (by a method which I can easily and efficiently scale for my larger 3D models).
Always a good way to start.
From what I have read on the forum, it means I should rule out using LinearMechanism. Do correct me if I am wrong here.
Done in a different post.
The most natural approach seemed to be to utilize 'extracellular' mechanism in NEURON. I would by default obtain the connection between the two grids via the cell membrane.
Yes. Use extracellular's i_membrane as current sources to drive the nodes in your extracellular grid, and couple the grid's local potentials back to the extracellular membrane via e_extracellular. The default values of extracellular's parameters are such that current through the xraxial and xc elements is vanishingly small, and voltage drop across the xg elements is negligible. Consequently, membrane current flows "radially" to ground through the e_extracellular voltage sources.
I believe, I would only need one layer of extracellular, i.e. nlayer=1.
True in theory, but why divert your effort to compiling from source code unless it's necessary--and you won't know if it's necessary until you run into space or speed issues. I'd advise sticking with plain vanilla NEURON for now.
The next concern was to have the extracellular space of adjacent cells connected/coupled
That is the purpose of your extracellular grid.
A seemingly weird thought that came up, considering that I had already coupled the cells via gap junctions, was to have a similar mechanism for coupling/connecting the extracellular resistances.
If you're thinking "how can I join the extracellular layers of one cell to those of another cell at the site of a gap junction," forget about it. Coupling between cells is to be done by your extracellular grid.

Some investigators have literally split the task into two parts, with NEURON simulating a cell and some other program simulating the surrounding conductive medium. Simulation execution then proceeds like so

Code: Select all

Initialize the model cell and model bulk conductor to t=0
REPEAT
  given the extracellular field and membrane currents at time t
    use NEURON to calculate the membrane currents at t+dt
    use a Laplace solver to calculate the field in the bulk conductor to t+dt
    exchange data between NEURON and the Laplace solver
      i.e. NEURON sends membrane currents at t+dt to the Laplace solver,
      and the Laplace solver sends the extracellular potentials at t+dt to NEURON
  update t to t+dt
UNTIL done
Could all of this be done in NEURON? In principle, yes. The chief problem is how to deal with extracellular space. The task is easiest if the conductive medium is resistive and isotropic, in which case a network model whose cells contain a total of N intracellular nodes (segments) requires an NxN matrix to represent the extracellular interactions between segments. The simplest implementation of a simulation run would then proceed like this:

Code: Select all

Initialize the cells and extracellular potentials at t=0
REPEAT
  given the e_extracellular values and membrane currents at time t
    advance all cells by dt
    from the new i_membrane values, calculate new e_extracellular values
      i.e. by e_ex = Rx * i_m where i_m is the vector of membrane currents,
        Rx is the matrix of transfer resistances between the model's segments,
        and e_ex is the vector of extracellular potentials
          that result from the membrane currents
  update t to t+dt
  update the e_extracellular values from the elements of e_ex
UNTIL done

shailesh
Posts: 104
Joined: Thu Mar 10, 2011 12:11 am

Re: Extracellular Resistance as a 3D Grid - Bidomain Model

Post by shailesh » Mon Jul 16, 2012 2:20 pm

A BIG thanks for the point by point response. I did see several other posts where extracellular was being modeled in another software (Matlab/Comsol/Maxwell etc). But I did prefer to stick to NEURON for the same and hence will try the second approach that you have suggested.

And yes I will go ahead with vanilla NEURON (nlayer=2) for the time being to keep things simple.

Thanks again!
Shailesh Appukuttan

shailesh
Posts: 104
Joined: Thu Mar 10, 2011 12:11 am

Connecting Extracellular Resistances Using LinearMechanism

Post by shailesh » Wed Jun 26, 2013 6:48 am

After a gap of around an year, I come back to this topic (last time around I found it a bit difficult to figure my way out and hence moved on to other things - but this time I don't have an option).

I had made up my mind to use NEURON + Python for doing this based on your suggestion; using Python for modeling the extracellular conductive medium. But I happened to come across this model on modeldb dealing with ephaptic coupling: http://senselab.med.yale.edu/modeldb/Sh ... \ephap.hoc. The basic concept, I believe, had to be common in both cases of connecting the extracellular spaces of two cells.

They modeled the same using 'LinearMechanism'. I am trying to adopt the same technique but my knowledge of 'LinearMechanism' is very patchy (picked up from the documentation, forum and scanning codes in modeldb) and hence would like your opinion. (I wonder whether this is equivalent to the procedure you suggested for getting everything done within NEURON)

As before, I started off simple with just two identical cells, named 'cella' and 'cellb', having 'extracellular' mechanisms. Below I have attached the relevant code as well as my interpretaion of what it is doing. Do let me know if I am wrong at any stage:

Code: Select all

Model:
              CELL A                                      CELL B
o/`--o--'\/\/`--o--'\/\/`--o--'\o       o/`--o--'\/\/`--o--'\/\/`--o--'\o vext + v
     |          |          |                 |          |          |
    ---        ---        ---               ---        ---        ---
   |   |      |   |      |   |             |   |      |   |      |   |
    ---        ---        ---               ---        ---        ---
     |          |          |                 |          |          |
     |          |          |                 |          |          |     i_membrane
     |          |          |                 |  xraxial |          |
 /`--o--'\/\/`--o--'\/\/`--o--'\#########/`--o--'\/\/`--o--'\/\/`--o--'vext
To be modeled: ######### (resistive link between the two extracellular spaces)

LinearMechanism is being used as follows:

Code: Select all

c*(dy/dt) + g*y = b
lm = new LinearMechanism(cmat, gmat, e, bvec, sl, xl, layer)

Here,
c = cmat = | 0  0 |           // As no capacitive element considered between extracellular spaces
           | 0  0 |

b= bvec = | 0 |               // Not sure when this is non-zero
          | 0 |

sl = | cella |                // Relevant sections
     | cellb |

xvec = xl = | 0.9999 |        // Location on sections; non-terminal extremities
            | 0.0001 |

layervec = layer = | 1 |        Hence, y = e = | cella.vext[0](0.9999) |
                   | 1 |                       | cellb.vext[0](0.0001) |


As c = cmat = 0 & b = 0, the equation comes down to:
g*y = 0
The thing I am most unsure about is the matrix for 'g' - the conductance matrix

From the following modeldb files:
> Migliore et. al 2005: http://senselab.med.yale.edu/modeldb/sh ... db\gap.hoc
> Bokil et. al 2001: http://senselab.med.yale.edu/modeldb/Sh ... \ephap.hoc

The general observation seems to be:

Code: Select all

g = gmat = |  ge   -ge |        where ge = extracellular conductance (more on that later)
           | -ge    ge |
I am not clear on how the above matrix for g arises. My understanding was that the non-diagonal elements should be 'ge' - as conductance = ge from cella.vext[0](0.9999) to cellb.vext[0](0.0001), and vice versa. Also, why do we have '-ge' for the diagonal elements here? Shouldn't they be just 0 as we do not want to add any conductance between cella.vext[0](0.9999) and itself; similary for cellb.vext[0](0.0001). Why -ge?

Suppose we did go ahead with the 'gmat' matrix as in above modeldb files, then expanding the equation:

Code: Select all

g*y = 0
|  ge   -ge |  x  | cella.vext[0](0.9999) | =  | 0 |
| -ge    ge |     | cellb.vext[0](0.0001) |    | 0 |

i.e.
 [ge*cella.vext[0](0.9999)] - [ge*cellb.vext[0](0.0001)] = 0
-[ge*cella.vext[0](0.9999)] + [ge*cellb.vext[0](0.0001)] = 0

Both the above equations essentially being the same, we get,
cella.vext[0](0.9999) = cellb.vext[0](0.0001)
>> This is certainly not what we want and neither what we could expect - as there is resistive pathway (=1/ge) between them which would cause a potential gradient.
At this point I am lost and can't figure our where I have gone wrong and messed up things.

The relevant hoc code is as follows:

Code: Select all

	objref gmat, cmat, bvec, e, xl, layer, sl, lm
	ncells = 2

	gmat = new Matrix(ncells, ncells, 2) // sparse
	cmat = new Matrix(ncells, ncells, 2) // sparse and 0
	bvec = new Vector(ncells) // also 0
	e = new Vector(ncells) // extracellular potential
	sl = new SectionList()
	xl = new Vector(ncells)
	layer = new Vector(ncells)
	layer.fill(1)

	cella sl.append()
	xl.x[0] = 0.9999
	cellb sl.append()
	xl.x[1] = 0.0001
	
	gmat.x[0][1] = 200e-1/area(0.9999)
	gmat.x[1][0] = 200e-1/area(0.9999)

	gmat.x[0][1] = -200e-1/area(0.9999)  //Not sure
	gmat.x[1][0] = -200e-1/area(0.9999)  //Not sure

	lm = new LinearMechanism(cmat, gmat, e, bvec, sl, xl, layer)
APPENDED LATER - The 'correct' approach?!
---------------------------------------------------
Approaching the current balance equations bottom-up (as should be the proper method I presume), for the model we have:

Let,
Iext_ab = Current flowing through extracellular resistive link from cella to cellb
Iext_ba = Current flowing through extracellular resistive link from cellb to cella
Va = cella.vext[0](0.9999)
Vb = cellb.vext[0](0.0001)
ge = conductance of extracellular resitive link

The current balance equations would be obtained as follows:

Code: Select all

Iext_ab = -Iext_ba
ge * (Va-Vb) = ge * (Vb-Va)
i.e.
geVa - geVb = geVb - geVa

Considering that b=0 in the modelDB files, I suppose are represented in matrix form as:
|  ge   -ge |  x  | Va | =  | 0 |
| -ge    ge |     | Vb |    | 0 |
I am not totally convinced with this derivation - especially the representation in the matrix form. Could you point out where I have gone wrong?

Also, a couple of things I wanted to confirm:
>> 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)

The primary question might be whether this LinearMechanism appraoch would work for connecting the extracellular spaces. If yes, then that seems best for me. But even if not, I am keen to understand how it works and how it has been used in the mentioned modeldb files.

I have tried to state the query as elaborately as possible so that it could also serve as a possible reference for future users of 'LinearMechanism' and/or these modeldb files.
Last edited by shailesh on Sun Aug 11, 2013 1:22 pm, edited 1 time in total.

shailesh
Posts: 104
Joined: Thu Mar 10, 2011 12:11 am

Re: Extracellular Resistance as a 3D Grid - Bidomain Model

Post by shailesh » Sun Aug 11, 2013 1:20 pm

For the last month or more, I have been trying my hand at using LinearMechanism. Despite having made some headway I am still some way off my goal. Currently, I am stuck with the following errors: "spFactor error: Zero Diagonal" & "spFactor error: Singular". But before going into those, I want to ensure that my understanding and implementation of LinearMechanism is correct. I had outlined my approach in the post above. Could someone have a look and let me know if I am headed down the right path? Thanks in advance!

p.s. the queries towards the end, regarding 'ge' and 'xraxial' units, have already been discussed here: viewtopic.php?f=12&t=549#p11541

ted
Site Admin
Posts: 5592
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Connecting Extracellular Resistances Using LinearMechani

Post by ted » Tue Aug 13, 2013 5:29 pm

Lacking time to reply to your questions in a single session, I decided to start by replying incrementally, expanding gradually until I address as many of your questions as I can.
shailesh wrote:

Code: Select all

b= bvec = | 0 |               // Not sure when this is non-zero
It will be nonzero if a particular equation
c*dy/dt + g*y = b
has a nonzero right hand side (think "driving function" or "signal source").
The general observation seems to be:

Code: Select all

g = gmat = |  ge   -ge |        where ge = extracellular conductance (more on that later)
           | -ge    ge |
I am not clear on how the above matrix for g arises.
Think about a circuit as being a graph in which nodes (locations at which v is calculated) are connected by edges (circuit elements through which current flows). Next consider a circuit in which two nodes are connected by a conductance G. Let's call them node 0 and node 1, so the voltages at these nodes will be called v0 and v1. The current that leaves node 0 through g is
(v0 - v1)*G
and the current that leaves node 1 through g is
(v1 - v0)*G
The v vector is
v0
v1
and the g matrix is
G -G
-G G

If these two nodes are part of a larger circuit, the g matrix for the entire circuit will of course have many more rows and columns, but the submatrix that represents the coupling of 0 and 1 by G will still be
G -G
-G G

shailesh
Posts: 104
Joined: Thu Mar 10, 2011 12:11 am

Re: Extracellular Resistance as a 3D Grid - Bidomain Model

Post by shailesh » Tue Aug 13, 2013 11:50 pm

Thanks Ted. So far so good - I follow you till that point. I attempted such a derivation in that post under the head:
APPENDED LATER - The 'correct' approach
Only query being what I have posted in the next part there...
Suppose we did go ahead with the 'gmat' matrix as in above modeldb files, then expanding the equation:
Looking forward to the next installment :)
p.s. Atleast I am convinced that I am headed in the right way... now keen to know why it was right!

ted
Site Admin
Posts: 5592
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Extracellular Resistance as a 3D Grid - Bidomain Model

Post by ted » Wed Aug 14, 2013 11:28 am

shailesh wrote:Only query being what I have posted in the next part there...
Suppose we did go ahead with the 'gmat' matrix as in above modeldb files, then expanding the equation:
The LinearMechanism's g matrix is not the same as the g matrix in the model's system equations cy' + gy = b. The LinearMechanism's g matrix merely adds new off-diagonal terms to the system equations' g matrix. Which is to say that the elements in the LinearMechanism's g matrix affect the system equations that govern certain y variables (voltages), but they do not prescribe what any particular voltage will do unless they constitute the complete system equation for that particular node. And they don't, at least not if the node is connected to other nodes or signal sources.

shailesh
Posts: 104
Joined: Thu Mar 10, 2011 12:11 am

Re: Extracellular Resistance as a 3D Grid - Bidomain Model

Post by shailesh » Wed Aug 14, 2013 11:45 am

Oh ok... I think I see it now! The equation:
c*(dy/dt) + g*y = b
just involved scalar multiplications, and not the matrix multiplications that I mixed in. I suppose I just got into the flow of things seeing all the matrices, and when I had it on paper, matrix multiplication seemed the most natural way to go about it - especially since the order of the matrices were also feasible for doing so.

(An example of how to complicate when things are actually simple.)
Thanks for your time!

shailesh
Posts: 104
Joined: Thu Mar 10, 2011 12:11 am

Re: Extracellular Resistance as a 3D Grid - Bidomain Model

Post by shailesh » Sun Aug 18, 2013 2:22 am

One more query to follow-up on:
consider a circuit in which two nodes are connected by a conductance G. Let's call them node 0 and node 1, so the voltages at these nodes will be called v0 and v1. The current that leaves node 0 through g is
(v0 - v1)*G
and the current that leaves node 1 through g is
(v1 - v0)*G
The v vector is
v0
v1
and the g matrix is
G -G
-G G
Matrix b is defined as
b= bvec = | 0 |
............| 0 |

So, the equation gy = b would turn out to be:
| G -G | * | v0 | = | 0 |
| -G G | * | v1 | = | 0 |

Is this correct? Because I am not sure how my earlier comment about scalar multiplication will apply here; and matrix multiplication would result in:
(v0 - v1) * G = 0
(v1 - v0) * G = 0
The above equations would be incorrect as their RHS = 0
And solving the above equations would result in v0 = v1

Sorry for the repeated queries on this, but I am still not able to understand.

ted
Site Admin
Posts: 5592
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Extracellular Resistance as a 3D Grid - Bidomain Model

Post by ted » Mon Aug 19, 2013 10:49 am

Your interpretation of the gap junction's g matrix would be correct if the g matrix were the system matrix, i.e. if the model consisted of two nodes and no sources. But that is not the case. There is an entire model with lots of nodes and many sources, and it has a large system matrix. The g matrix merely specifies values that are added to elements of the system matrix.

shailesh
Posts: 104
Joined: Thu Mar 10, 2011 12:11 am

Re: Extracellular Resistance as a 3D Grid - Bidomain Model

Post by shailesh » Tue Aug 20, 2013 12:39 pm

Ok...... I think I got it this time. A more general query:
I was just wondering that there are so many models on modeldb, but I have rarely come across models where:
> there are two or more cells and
> all these cells have extracellular fields (ecf) and
> the ecfs of these cells should be linked (i.e. it is required that the ecf of one cell should influence the ecf of another cell (and possibly, in turn, its vm))

Is this really so rarely undertaken in computational modeling? Or am I missing out on several modeldb entries...

The only related modeldb entry I found was this:
> Ephaptic interactions in olfactory nerve
http://senselab.med.yale.edu/modeldb/Sh ... Cephap.hoc

ted
Site Admin
Posts: 5592
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Extracellular Resistance as a 3D Grid - Bidomain Model

Post by ted » Tue Aug 20, 2013 2:18 pm

shailesh wrote:Is this really so rarely undertaken in computational modeling?
That is the case. Until recently most models have ignored such effects.

Post Reply