Extracellular mechanism with variable e_extracellular

[Keivan]

I am trying to make a model of cardiac action potential propagation. Ephaptic coupling has a major role in the propagation of action potentials in these cells. In a simple model of ephaptic coupling, the potential in the cleft between two cells should be calculated. In order to do that I have made a mod file that I include it hereunder:

Code: Select all

COMMENT
Potential in the Intercalated disk
This mechanism should be inserted in the pre-junctional compartment

phi_c = A*Rc*Ro*(-dv_pre/dx+dv_post/dx)

Keivan Moradi, 2013

ENDCOMMENT

NEURON {
    POINT_PROCESS vDisk
    POINTER v2, vp1, vp2
    RANGE vc
}

UNITS {
	(nA) = (nanoamp)
	(mV) = (millivolt)
	(mS) = (millisiemens)
	(um) = (micrometer)
	(kO) = (kiloohm)
}

PARAMETER {
	dx = 2.75  (um)    : the length of pre and post junctional compartments that 
			           : should be a fixed value in the pre-junctional and junctional compartments
	A  = 380   (um2)   : The area of the cell cross section
	Rc = 1e4   (kO)    : 1e3 to 1e5 kilo-ohm the cleft-to-ground resistance
	Ro = 6.67  (mS/cm) : mS/cm the cytosolic conductivity 
}
    
ASSIGNED {
    v    (mV)
	v1   (mV) : voltage in the first compartment of the pre-junctional myocyte
	v2   (mV) : voltage in the parent compartment of the pre-junctional myocyte
	vp1  (mV) : voltage in the first compartment of the post-junctional myocyte
	vp2  (mV) : voltage in the parent compartment of the post-junctional myocyte
	vc   (mV) : the voltage in the cleft
}
 
BREAKPOINT {
	v1 = v
	vc = (0.0001)*(Ro*Rc*A*(-v1+v2+vp1-vp2)/dx)
}

now I want to connect the vc of the above mechanism to the e_extracellular of my pre and post-junctional compartments. Is it possible? If not, is it possible to hack NEURON's source code to do that?

I am going to test the effect of different ion channels in the cleft membrane to test their role in the propagation of action potential in a 100*100 grid of cardiac myocytes. This is the hoc code that I have written so far.

Code: Select all

load_file("nrngui.hoc")

begintemplate Myocyte
	objref disk, branch, Vdisk
	create body, RU, RD, LU, LD, RUd, RDd, LUd, LDd
	public body, RU, RD, LU, LD, RUd, RDd, LUd, LDd, conn, disk, branch, Vdisk
	
	proc init() {
		// connect section(0or1), parent(x)
		connect LU(1), body(0) 
		connect LD(1), body(0)
		connect RU(0), body(1)
		connect RD(0), body(1)
		connect LUd(1), LU(0)
		connect LDd(1), LD(0)
		connect RUd(0), RU(1)
		connect RDd(0), RD(1)
		
		disk = new SectionList()
		LUd disk.append()
		LDd disk.append()
		RUd disk.append()
		RDd disk.append()
		
		branch = new SectionList()
		LU branch.append()
		LD branch.append()
		RU branch.append()
		RD branch.append()
		
		body {
			nseg = 1
			diam = 22
			L = 67		// micrometer
			Ra = 123.0
			insert pas
			g_pas = 0.0001666
			e_pas = -60.0
		}
		
		forsec branch {
			nseg = 3
			// the branch diameter is half the body
			diam = 11
			// at L = diam/4 the area of the compartment cross-section is equal 
			// to the lateral area of the compartment
			L    = diam * nseg / 4
			Ra   = 123
			insert pas
			g_pas = 0.0001666
			e_pas = -60.0
		}
		
		forsec disk {
			nseg = 1
			diam = 11
			L    = diam / 4
			insert extracellular
		}
	}
	proc conn() {
		RUd Vdisk = new vDisk(0.5)
		setpointer Vdisk.v2, RU.v(0.75)
		setpointer Vdisk.vp1, $o1.LDd.v(0.5)
		setpointer Vdisk.vp2, $o1.LD.v(0.75)
		// find a way to connect Vdisk.vc to e_extracellular of RUd and $o1.LDd ???
	}
endtemplate Myocyte

objref myocyte[2]
for i=0,1 myocyte[i] = new Myocyte()
for i=1,1 {
	myocyte[i-1].conn (myocyte[i])
}

[ted]

I want to connect the vc of the above mechanism to the e_extracellular of my pre and post-junctional compartments.

You're off to a good start. Declare vc to be a POINTER, not a RANGE variable. Link it to the e_extracellular of your pre and post-junctional compartments with setpointer statements.

[Keivan]

Thank you Ted.
I have changed my code the way you suggested.

Code: Select all

COMMENT
Potential in the Intercalated disk
This mechanism should be inserted in the pre-junctional compartment

phi_c = A*Rc*Ro*(-dv_pre/dx+dv_post/dx)

Keivan Moradi, 2013

ENDCOMMENT

NEURON {
    POINT_PROCESS vDisk
    POINTER v2, vp1, vp2, vc
}

UNITS {
	(nA) = (nanoamp)
	(mV) = (millivolt)
	(mS) = (millisiemens)
	(um) = (micrometer)
	(kO) = (kiloohm)
}

PARAMETER {
	dx = 2.75  (um)    : the length of pre and post junctional compartments that 
			           : should be a fixed value in the pre-junctional and junctional compartments
	A  = 380   (um2)   : The area of the cell cross section
	Rc = 5.5   (kO)    : 1e3 to 1e5 kilo-ohm the cleft-to-ground resistance
	Ro = 6.67  (mS/cm) : the cytosolic conductivity 
}
    
ASSIGNED {
    v    (mV)
	v1   (mV) : voltage in the first compartment of the pre-junctional myocyte
	v2   (mV) : voltage in the parent compartment of the pre-junctional myocyte
	vp1  (mV) : voltage in the first compartment of the post-junctional myocyte
	vp2  (mV) : voltage in the parent compartment of the post-junctional myocyte
	vc   (mV) : the voltage in the cleft
}

INITIAL {
	vc = 0
}
 
BREAKPOINT {
	v1 = v
	vc = (0.0001)*(Ro*Rc*A*(-v1+v2+vp1-vp2)/dx)
}

Code: Select all

load_file("nrngui.hoc")

begintemplate Myocyte
	objref disk, branch, Vdisk
	create body, RU, RD, LU, LD, RUd, RDd, LUd, LDd
	public body, RU, RD, LU, LD, RUd, RDd, LUd, LDd, conn, disk, branch, Vdisk
	
	proc init() {
		// connect section(0or1), parent(x)
		connect LU(1), body(0) 
		connect LD(1), body(0)
		connect RU(0), body(1)
		connect RD(0), body(1)
		connect LUd(1), LU(0)
		connect LDd(1), LD(0)
		connect RUd(0), RU(1)
		connect RDd(0), RD(1)
		
		disk = new SectionList()
		LUd disk.append()
		LDd disk.append()
		RUd disk.append()
		RDd disk.append()
		
		branch = new SectionList()
		LU branch.append()
		LD branch.append()
		RU branch.append()
		RD branch.append()
		
		body {
			nseg = 1
			diam = 22
			L = 67		// micrometer
			Ra = 50
			insert pas
			g_pas = 0.0001666
			e_pas = -60.0
		}
		
		forsec branch {
			nseg = 3
			// the branch diameter is half the body
			diam = 11
			// at L = diam/4 the area of the compartment cross-section is equal 
			// to the lateral area of the compartment
			L    = diam * nseg / 4
			Ra   = 50
			insert pas
			g_pas = 0.0001666
			e_pas = -60.0
		}
		
		forsec disk {
			nseg = 1
			diam = 11
			L    = diam / 4
			Ra   = 50
			insert extracellular
				e_extracellular = 0
				xraxial[0] = 1e9
				xg[0] = 1e9
				xc[0] = 0
				xraxial[1] = 1e9
				xg[1] = 1e9
				xc[1] = 0
		}
		
		forall insert hh
	}
	proc conn() {
		RUd Vdisk = new vDisk(0.5)
		setpointer Vdisk.v2, RU.v(0.75)
		setpointer Vdisk.vp1, $o1.LDd.v(0.5)
		setpointer Vdisk.vp2, $o1.LD.v(0.75)
		setpointer Vdisk.vc, RUd.e_extracellular(0.5)
		setpointer Vdisk.vc, $o1.LDd.e_extracellular(0.5)
	}
endtemplate Myocyte

// create the main dendrite
// create Mdend
// Mdend{
	// nseg = 5
	// diam = 318
	// L = 60
	// Ra = 123
	// insert pas
	// g_pas = 0.0001666
	// e_pas = -60.0
// }
// objref ORdend
// Mdend ORdend = new SectionRef()

objref myocyte[2]
for i=0,1 myocyte[i] = new Myocyte()
for i=1,1 {
	myocyte[i-1].conn(myocyte[i])
}

// myocyte.disk.printnames()
// topology()
// objref s
 // s = new Shape()

// create an electrode in the soma
objref stim
myocyte[0].body stim = new IClamp(0.5)
stim.del = 100
stim.dur = 100
stim.amp = 1

load_file("test.ses")
tstop = 500
v_init = -64
// run()

test.ses:

Code: Select all

{load_file("nrngui.hoc")}
objectvar save_window_, rvp_
objectvar scene_vector_[3]
objectvar ocbox_, ocbox_list_, scene_, scene_list_
{ocbox_list_ = new List()  scene_list_ = new List()}
{pwman_place(0,0,0)}
{
save_window_ = new Graph(0)
save_window_.size(0,500,-80,40)
scene_vector_[2] = save_window_
{save_window_.view(0, -80, 500, 120, 451, 66, 432, 246.7)}
graphList[0].append(save_window_)
save_window_.save_name("graphList[0].")
save_window_.addvar("myocyte[1].LDd.v( 0.5 )", 2, 1, 0.8, 0.9, 2)
save_window_.addvar("myocyte[0].RUd.v( 0.5 )", 1, 1, 0.8, 0.9, 2)
save_window_.addvar("myocyte[0].body.v( 0.5 )", 3, 1, 0.8, 0.9, 2)
}
{
xpanel("RunControl", 0)
v_init = -64
xvalue("Init","v_init", 1,"stdinit()", 1, 1 )
xbutton("Init & Run","run()")
xbutton("Stop","stoprun=1")
runStopAt = 5
xvalue("Continue til","runStopAt", 1,"{continuerun(runStopAt) stoprun=1}", 1, 1 )
runStopIn = 1
xvalue("Continue for","runStopIn", 1,"{continuerun(t + runStopIn) stoprun=1}", 1, 1 )
xbutton("Single Step","steprun()")
t = 500
xvalue("t","t", 2 )
tstop = 500
xvalue("Tstop","tstop", 1,"tstop_changed()", 0, 1 )
dt = 0.025
xvalue("dt","dt", 1,"setdt()", 0, 1 )
steps_per_ms = 40
xvalue("Points plotted/ms","steps_per_ms", 1,"setdt()", 0, 1 )
screen_update_invl = 0.05
xvalue("Scrn update invl","screen_update_invl", 1,"", 0, 1 )
realtime = 1.2
xvalue("Real Time","realtime", 0,"", 0, 1 )
xpanel(1080,26)
}
objectvar scene_vector_[1]
{doNotify()}

Do you mean that if I use the above command then by changing Vdisk.vc the value of RUd.e_extracellular(0.5) would change? This is in contradiction of what I knew about the NEURON. Are you sure about this? If this is true, then I should be able to change v, vext, t and most of the internal variables of the neuron. Is this true?

I have tried hard not to involve myself with the scary LinearMechanism of NEURON, but it seems that it is inevitable. Could you please help me to understand it? I do not understand the meaning of c and g matrices and their relation with the traditional RC circuit we solve in every compartment. I have consulted with some of my electrical engineer friends, but these matrices seemed different with their knowledge of the matrix solution in electrical circuits. If it is hard to explain in this forum, would you please refer me to a good and simple documentation or tutorial?
For example, I want to learn how to wire the extracelluar part of two compartments to each other? Or, how can I connect a node in the extracellular mechanism to the other node in the other extracellular mechanism? How can I make my own extracellular mechanism with it? and something like this.

[shailesh]

Hi Keivan,

Happened to come across your post as I was also working on modeling extracellular fields and found the title interesting.
I also ended up having to deal with LinearMechanism.... and yeah it seemed best avoided :p (atleast at first)

You might find the following discussion useful:
http://www.neuron.yale.edu/phpBB/viewto ... f=8&t=2654
Especially Ted's responses!
There seem to be no tutorials around for LinearMechanism.

Cheers,
Shailesh Appukuttan

[ted]

Do you mean that if I use the above command then by changing Vdisk.vc the value of RUd.e_extracellular(0.5) would change?

To quote the Programmer's Reference entry for setpointer:

Connects pointer variables in membrane mechanisms to the address of var. eg. If $NEURONHOME/examples/nmodl/synpre.mod in linked into NEURON, then:

Code: Select all

    soma1 syn1=new synp(.5)
    setpointer syn1.vpre, axon2.v(1)

would enable the synapse in soma1 to observe the axon2 membrane potential.

I use this feature for two purposes in xtra.mod. First, it is used to couple an extracellular stimulus current to local extracellular potential. Here are the relevant parts of xtra.mod

Code: Select all

NEURON {
	SUFFIX xtra
	RANGE rx : transfer resistance between stim/rec electrode and axon
  . . .
	GLOBAL is : stimulus current
	POINTER ex : link this to local e_extracellular
}
PARAMETER {
	rx = 1 (megohm)
 . . .
INITIAL {
	ex = is*rx*(1e6)
  . . .
BEFORE BREAKPOINT { : before each cy' = f(y,t) setup
  ex = is*rx*(1e6)
}

and here is the hoc code that sets up the link from the hoc-specified stimulus current variable to

Code: Select all

forall {
  if (ismembrane("xtra")) {
    for (x, 0) {
      . . .
      setpointer ex_xtra(x), e_extracellular(x) : ex_xtra drives e_extracellular
    }
  }
}

Second, it is used to couple each compartment's local membrane current back to the cell-generated potential that would be observed at the stimulating/recording electrode. Here are the relevant parts of xtra

Code: Select all

NEURON {
	SUFFIX xtra
	RANGE er
  . . .
	POINTER im
}
INITIAL {
  . . .
	er = (10)*rx*im*area
}
AFTER SOLVE { : after each solution step
  er = (10)*rx*im*area
}

and the hoc statements that link membrane current to xtra's im variable.

Code: Select all

forall {
  if (ismembrane("xtra")) {
    for (x, 0) {
      setpointer im_xtra(x), i_membrane(x)
      . . .
    }
  }
}

If this is true, then I should be able to change v, vext, t and most of the internal variables of the neuron.

Yes and no, or more properly, yes as long as you don't introduce a conflict or force an abrupt state variable change during adaptive integration. For example, changing t is done in the "initialize to steady state" custom initialization, but one wouldn't do it in the middle of a simulation run, nor can I imagine a case in which it would be useful to do that with a POINTER variable. With regard to changing state variables, that is possible but you have to be aware of the fact that NEURON might just overwrite the value you assigned. Also, if you change a state variable in the middle of a simulation that uses adaptive integration, you need to call cvode.re_init() because the change is equivalent to starting a new initial value problem.

With regard to representing intercalated disks, I think the LinearMechanism has big advantages because it would allow a much simpler model specification. From your NMODL code I would guess that a disk has an equivalent circuit something like this

Code: Select all

A                     B
o--- ra ---+--- rb ---o
           |
           rc
           |
          gnd

where the A terminal is attached to an end of one myofibril, and the B terminal is attached to an end of another myofibril. Is that correct?

[Keivan]

Thank you Ted, I would test the suggested code and reply soon.
Thank you Shailesh for pointing me to that topic.

From your NMODL code I would guess that a disk has an equivalent circuit something like this

Code: Select all
A B
o--- ra ---+--- rb ---o
|
rc
|
gnd

where the A terminal is attached to an end of one myofibril, and the B terminal is attached to an end of another myofibril. Is that correct?

The complete equivalent circuit is shown in Fig. 1 of Hand and Griffith, 2010. My nmodl model, however, tries to deal with the ephaptic coupling only.
I have modified that figure to make a new figure that shows the intercalated disk between two cardiac myocytes in my view. The A and B parts of my figure are normal shape of the equivalent circuit in hand and Griffith, 2010, but I want the one I show in the C part. Technically, B and C are the same, but to me, C is modeled with a special section of neuron that is connected to another section both in the intracellular and the extracellular parts. The connection in the intracellular part models the gap junction coupling and the connection in the extracellular part models the ephaptic coupling. This way, I can not only add and remove whatever ion channel I like to the intercalated disk with the normal insert command of NEURON, but also I have access to vext (phi_c or vc in my nmodl model) that I need it in my future models. In addition, I want to be able to model both the gap junction and the ephaptic couplings separate from each other.

[ted]

It is very easy to use LinearMechanism to implement gap junctions. Here's an example:

Code: Select all

load_file("nrngui.hoc")

create a, b // will be two separate cells
a { // use pt3d specification of geometry to control where the cells are drawn in the GUI
  pt3dclear()
  pt3dadd(0,10,0,1) // a will be at top left
  pt3dadd(100,10,0,1.5)
}
b {
  pt3dclear()
  pt3dadd(100,-10,0,1) // b will be at bottom right
  pt3dadd(200,-10,0,1) // different diameter to make it interesting
}
forall {
  Ra = 100 // ohm cm
  insert pas
  g_pas = 1e-4
  insert extracellular
}
a nseg = 3 // since its diam is bigger
b nseg = 5

load_file("gap.hoc") // from ModelDB entry 43039, defines Gap class
objref gap
gap = new Gap()
aloc = 0.999 // where gap junction will be connected
bloc = 0.001
ggap = 15 // gap junction conductance in nS
a gap.src(aloc)
b gap.target(bloc)
gap.g(ggap)

load_file("gaprig.ses") // RunControl, IClamp 0.03 nA x 1 ms attached to a(0)
// v axis graph that shows v at a(0.001), a(1), b(0.001), and b(1)

You'll have to get gap.hoc from ModelDB entry 43039, and a ses file called gaprig.ses that sets up
* a RunControl
* an IClamp attached to a(0) that delivers 0.03 nA x 1 ms
* a v axis graph that shows v at a(0.001), a(1), b(0.001), and b(1)

Note that Gap.g expects an argument that specifies the gap junction conductance in nanosiemens.

I'll have to think about how to make a new class that uses LinearMechanism to implement the ephaptic coupling you want.

[Keivan]

As far as I know, ephaptic coupling by itself is the sufficient condition for the propagation of action potential from one myocyte to another myocyte. My ephaptic coupling mechanism, however, is unstable. The suggested modifications, as in the xtra mechanism, did not solved the problem. Do you have any idea why this model is so unstable?
mod file:

Code: Select all

COMMENT
Potential in the Intercalated disk
This mechanism should be inserted in the pre-junctional compartment

phi_c = A*Rc*Ro*(-dv_pre/dx+dv_post/dx)

Keivan Moradi, 2013

ENDCOMMENT

NEURON {
    SUFFIX vDisk
    POINTER v1, v2, vp1, vp2, vc
    GLOBAL Rc
}

UNITS {
	(nA) = (nanoamp)
	(mV) = (millivolt)
	(mS) = (millisiemens)
	(um) = (micrometer)
	(kO) = (kiloohm)
}

PARAMETER {
	dx = 2.75  (um)    : the length of pre and post junctional compartments that 
			           : should be a fixed value in the pre-junctional and junctional compartments
	A  = 380   (um2)   : The area of the cell cross section
	Rc = 1e4   (kO)    : 1e3 to 1e5 kilo-ohm the cleft-to-ground resistance
	Ro = 6.67  (mS/cm) : the cytosolic conductivity 
}
    
ASSIGNED {
    v    (mV)
	v1   (mV) : voltage in the first compartment of the pre-junctional myocyte
	v2   (mV) : voltage in the parent compartment of the pre-junctional myocyte
	vp1  (mV) : voltage in the first compartment of the post-junctional myocyte
	vp2  (mV) : voltage in the parent compartment of the post-junctional myocyte
	vc   (mV) : the voltage in the cleft
}

INITIAL {
	vc = 0
}
 
BREAKPOINT {
	vc = (0.0001)*(Ro*Rc*A*(v1-v2-vp1+vp2)/dx)
}

hoc file

Code: Select all

load_file("nrngui.hoc")

begintemplate Myocyte
	objref disk, branch
	create body, RU, RD, LU, LD, RUd, RDd, LUd, LDd
	public body, RU, RD, LU, LD, RUd, RDd, LUd, LDd, conn, disk, branch
	
	proc init() {
		// connect section(0or1), parent(x)
		connect LU(1), body(0) 
		connect LD(1), body(0)
		connect RU(0), body(1)
		connect RD(0), body(1)
		connect LUd(1), LU(0)
		connect LDd(1), LD(0)
		connect RUd(0), RU(1)
		connect RDd(0), RD(1)
		
		disk = new SectionList()
		LUd disk.append()
		LDd disk.append()
		RUd disk.append()
		RDd disk.append()
		
		branch = new SectionList()
		LU branch.append()
		LD branch.append()
		RU branch.append()
		RD branch.append()
		
		body {
			nseg = 1
			diam = 22
			L = 67		// micrometer
			Ra = 150	// ohm-cm
			insert pas
			g_pas = 0.0001666
			e_pas = -60.0
			insert hh
		}
		
		forsec branch {
			nseg = 3
			// the branch diameter is half the body
			diam = 11
			// at L = diam/4 the area of the compartment cross-section is equal 
			// to the lateral area of the compartment
			L    = nseg * diam / 4
			Ra   = 150
			insert pas
			g_pas = 0.0001666
			e_pas = -60.0
			insert hh
		}
		
		forsec disk {
			nseg = 1
			diam = 11
			L    = diam / 4
			Ra   = 150
			insert extracellular
				e_extracellular = 0
				xraxial[0] = 1e9
				xg[0] = 1e9
				xc[0] = 0
				xraxial[1] = 1e9
				xg[1] = 1e9
				xc[1] = 0
		}
	}
	proc conn() {
		RUd { 
			insert vDisk
				Rc_vDisk = 10.94
				setpointer v1_vDisk(0.5), RUd.v(0.5)
				setpointer v2_vDisk(0.5), RU.v(0.75)
				setpointer vp1_vDisk(.5), $o1.LDd.v(0.5)
				setpointer vp2_vDisk(.5), $o1.LD.v(0.75)
			
				setpointer vc_vDisk(.5), RUd.e_extracellular(0.5)
				setpointer vc_vDisk(.5), $o1.LDd.e_extracellular(0.5)
			
				// setpointer vc_vDisk(.5), RUd.vext[0](0.5)
				// setpointer vc_vDisk(.5), $o1.LDd.vext[0](0.5)
		}
	}
endtemplate Myocyte

objref myocyte[2]
for i=0,1 myocyte[i] = new Myocyte()
myocyte[0].conn(myocyte[1])

// create an electrode in the cell body
objref stim
myocyte[0].body stim = new IClamp(0.5)
	stim.del = 100
	stim.dur = 10
	stim.amp = 1

load_file("test.ses")
tstop = 500
v_init = -64
run()

ses file

Code: Select all

{load_file("nrngui.hoc")}
objectvar save_window_, rvp_
objectvar scene_vector_[3]
objectvar ocbox_, ocbox_list_, scene_, scene_list_
{ocbox_list_ = new List()  scene_list_ = new List()}
{pwman_place(0,0,0)}
{
save_window_ = new Graph(0)
save_window_.size(0,500,-80,40)
scene_vector_[2] = save_window_
{save_window_.view(0, -80, 500, 120, 450, 66, 432, 246.7)}
graphList[0].append(save_window_)
save_window_.save_name("graphList[0].")
save_window_.addvar("myocyte[1].LDd.v( 0.5 )", 2, 1, 0.8, 0.9, 2)
save_window_.addvar("myocyte[0].RUd.v( 0.5 )", 1, 1, 0.8, 0.9, 2)
save_window_.addvar("myocyte[0].body.v( 0.5 )", 3, 1, 0.8, 0.9, 2)
save_window_.addvar("myocyte[1].LD.v( 0.833333 )", 7, 1, 0.8, 0.9, 2)
}
{
xpanel("RunControl", 0)
v_init = -64
xvalue("Init","v_init", 1,"stdinit()", 1, 1 )
xbutton("Init & Run","run()")
xbutton("Stop","stoprun=1")
runStopAt = 5
xvalue("Continue til","runStopAt", 1,"{continuerun(runStopAt) stoprun=1}", 1, 1 )
runStopIn = 1
xvalue("Continue for","runStopIn", 1,"{continuerun(t + runStopIn) stoprun=1}", 1, 1 )
xbutton("Single Step","steprun()")
t = 500
xvalue("t","t", 2 )
tstop = 500
xvalue("Tstop","tstop", 1,"tstop_changed()", 0, 1 )
dt = 0.025
xvalue("dt","dt", 1,"setdt()", 0, 1 )
steps_per_ms = 40
xvalue("Points plotted/ms","steps_per_ms", 1,"setdt()", 0, 1 )
screen_update_invl = 0.05
xvalue("Scrn update invl","screen_update_invl", 1,"", 0, 1 )
realtime = 1.16
xvalue("Real Time","realtime", 0,"", 0, 1 )
xpanel(1080,25)
}
objectvar scene_vector_[1]
{doNotify()}

[Keivan]

Today, I spend some time to understand the LinearMechanism by reading the topic that was suggested by Shailesh. It was really enlightening. I wish that information was explained in the documentation. Anyway, as far as I understood, using LinearMechanism, we can add conductance and capacitance between two nodes in the model. I think my million dollar question is "how we can add a new shared node between two sections?"

[Keivan]

Since it is easy to wire two nodes in the model to each other, I made my ephaptic coupling mechanism.

Code: Select all

// This mechanism is very simple. 
// It connects two extracellular voltage nodes to each other
// with a low resistance wire so that they become a single node (theoretically)
// then we connect this voltage node with a resistor (xg) to earth (the extracellular solution)
// the extracelluar mechanisms with xc[0] = xc[1] = 0 do this

// originally by Bokil et. al. (2001), J. Neurosci. 21:RC173(1-5)
// 2013 : Keivan Moradi:
//		1. Modified to make it a template
//		2. Now, the xg[0] plays the role of resistance to ground 
//			so that we can monitor the vext[0] changes that shows 
//			the potential in the gap between two cells

begintemplate Ephaptic
objref gmat, cmat, bvec, e, xl, layer, sl, lm

proc init(){local i, ns
	// order is all a then all b, the e nodes are equivalent to a.vext
	sl = new SectionList()
	sl.append()
	insert extracellular
		// $2 = cleft to ground conductance in nanosiemens
		// if gmat -> inf ==> $2=(1/2*$2)*2
		xg[0] = $2
		xc[0] = 0
		// Since the section before the current section do not have an extracellular mechanism,  
		// the xraxial do not have any impact on the simulation
		xraxial[0] = 1e9
		xg[1] = $2
		xc[1] = 0
		xraxial[1] = 1e9
		e_extracellular = 0
	$o1.sec {
		insert extracellular  
			xg[0] = $2
			xc[0] = 0
			xraxial[0] = 1e9
			xg[1] = $2
			xc[1] = 0
			xraxial[1] = 1e9
			e_extracellular = 0
		sl.append()
	}
	gmat = new Matrix(2, 2, 2) // sparse
		gmat.x[0][0] =  $3
		gmat.x[0][1] = -$3
		gmat.x[1][1] =  $3
		gmat.x[1][0] = -$3
	
	cmat = new Matrix(2, 2, 2) // sparse and 0
	bvec = new Vector(2) // also 0
	e = new Vector(2) // extracellular potential
	
	// Since we assume that the nseg is one in the gap junction section x = .5
	xl = new Vector(2)
		xl.x[0] = .5
		xl.x[1] = .5
	layer = new Vector(2)
		layer.fill(1) // all connections are between a.vext[0] and b.vext[0]
	
	lm = new LinearMechanism(cmat, gmat, e, bvec, sl, xl, layer)
}
endtemplate Ephaptic

The test code:

Code: Select all

create a, b
objref secref, ephaptic
b secref = new SectionRef()
a ephaptic = new Ephaptic(secref,1e-3,1e9)

It seems, this mechanism works for me. Because I want to keep it stupidly simple (KISS), I still want to learn if it is possible to make a real shared node in the model? I also want to be sure about the value of the internal voltage node. Is it v+vext? If so, I think I found the problem of my previous nmodl mechanism.