Gap Junction Linear Mechanism cvode

Anything that doesn't fit elsewhere.
Post Reply
oren
Posts: 55
Joined: Fri Mar 22, 2013 1:03 am

Gap Junction Linear Mechanism cvode

Post by oren »

Hello,
I have a network simulation with gap junction, I am using pc.source pc.target for the Gap Junctions.
I want to test if my results stay the same when I use linear mechaisem instead of pc.source and pc.target, but it is extremely slow, so I tried to cvode ( is it possible to use cvode and linear mechanism? )

after I do cvod.active(1) (With the linear mechanism) and then run() I get the following error:

Code: Select all

NEURON: NrnDAEs only work with secondorder==0 or daspk
 near line 0
 {LinearMGap[832].set_gm()}
                           ^
        finitialize(-58.1897)
      init()
    stdinit()
  run()
I don't know where to start with this error massage,



Another question, I read once that when the Gap Junction conductance is not to high there will be no difference between linear mechanism and pc.transfer , if this is true why?

Thank You,
Oren.

EDIT

this is the code that I use for the Gap Junction

Code: Select all

begintemplate LinearMGap
public src, target, g, valid, pr, set_gm, gm, cm
objref srcsec, targetsec, cm, gm, y, b, xvec, sl, lm, this, fih
strdef tstr

proc init(){
	g_ = 10
	valid_ = 0
	cm = new Matrix(2,2,2)
	gm = new Matrix(2,2)
	y = new Vector(2)
	b = new Vector(2)
	xvec = new Vector(2)
}

func src() {
	srcsec = new SectionRef()
	xvec.x[0] = $1
	valid_ = 0
	return valid()
}

func target() {
	targetsec = new SectionRef()
	xvec.x[1] = $1
	valid_ = 0
	return valid()
}

func g() {
	if (numarg() > 0) {
		g_ = $1
		if (valid_) {
			set_gm()
		}
	}
	return g_
}

proc set_gm() { local us, a// conductance in nanosiemens
	if (valid_ == 0) { return }
	us = .001*g_
	srcsec.sec { a = 100/area(xvec.x[0]) }
	gm.x[0][0] = us*a
	gm.x[0][1] = -us*a
	targetsec.sec { a = 100/area(xvec.x[1]) }
	gm.x[1][1] = us*a
	gm.x[1][0] = -us*a
}

func valid() {
	if (valid_ == 0) {
		if (object_id(srcsec) && object_id(targetsec)) {
			mkgap()		
		}
	}
	return valid_
}

proc mkgap() {
	sl = new SectionList()
	srcsec.sec sl.append()
	targetsec.sec sl.append()
	valid_ = 1
	set_gm()
	lm = new LinearMechanism(cm, gm, y, b, sl, xvec)
	// only necessary because we anticpate changes in diameter
	// fih = new FInitializeHandler(0, "set_gm()")
	// unfortunately up through the  5.6 2004/02/09 Main (44)
	// version there is an error in parsing the third arg, so
	sprint(tstr, "%s.set_gm()", this)
	fih = new FInitializeHandler(0, tstr)
}

proc pr() {
	if (valid_) {
srcsec.sec printf("%s %s(%g)", this, secname(), xvec.x[0])
targetsec.sec printf("---%s(%g) \tg = %g (ns)\n", secname(), xvec.x[1], g())
	}else{
		printf("%s not used\n", this)
	}
}
endtemplate LinearMGap
ted
Site Admin
Posts: 6300
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Gap Junction Linear Mechanism cvode

Post by ted »

oren wrote:is it possible to use cvode and linear mechanism?
Yes. I do it all the time.
after I do cvod.active(1) (With the linear mechanism) and then run() I get the following error:
. . .
NEURON: NrnDAEs only work with secondorder==0 or daspk
In your program, change
cvode.active(1)
to
cvode_active(1)
and see if that fixes the problem.
I read once that when the Gap Junction conductance is not to high there will be no difference between linear mechanism and pc.transfer , if this is true why?
You probably refer to something related to this statement
"Gap junctions are assumed to couple cells relatively weakly so that the modified euler method is acceptable for accuracy and stability."
in the documentation of Parallel Transfer methods
http://www.neuron.yale.edu/neuron/stati ... eltransfer
If you use linear mechanism to represent a gap junction, the model's Jacobian will include the appropriate off-diagonal terms that represent the conductance between the coupled compartments, and advancing from one time step to the next can use a fully implicit integration method. If instead you use parallel transfer to implement gap junctions, the Jacobian will not include these terms, and when the solution is advanced step t to t+dt, the current between the coupled compartments will have the value that was computed at time t. In other words, the solution will be generated by a modified Euler method. The difference between the two solutions depends on the degree of coupling; weak coupling means the conductances are small, so the presence or absence of the off-diagonal terms will have little effect on the stability or accuracy of the solution.
Post Reply