OK with fixed step but not with cvode solvers

NMODL and the Channel Builder.
Post Reply
tardebut
Posts: 21
Joined: Thu Sep 18, 2014 8:23 am

OK with fixed step but not with cvode solvers

Post by tardebut »

Hi there,

I'm running a neuron model fetched from ModelDB (149100) to which I have added three new mechanism (I mean .mod files) to simulate the effects of neuromodulators. When I run the model using a fixed step solver (cvode.active(0)) the results are the expected: the spike train amplitude and frequency change according to the neuromodulation dynamics. As expected this simulation is very slow. Then, when I try an adaptative step solver (cvode.active(1)) the effects of neuromodulation are gone but still I got the spike train trace with no neuromodulatory effects.
Two of the three inserted neuromodulation mechanism don't have any ODE and their outputs feed into the third mechanism which has several ODEs (for the signaling cascade) and an output that is the one modifying the existing channel mechanisms (max conductances) of the model.

Why the neuromodulatory effects are not properly integrated using the adaptative step solver?

I look forward to hear from you.

Thanks in advance,

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

Re: OK with fixed step but not with cvode solvers

Post by ted »

Scrying bugs in unseen code is ordinarily very difficult, especially in a darkened LCD monitor, but Allhallows falls on Friday this week so perhaps I can hazard this guess: your NMODL code contains statements similar to this:

Code: Select all

if (t>tstart) {
  someparameter = value1
}
Guaranteed to fail with adaptive integration because the integrator doesn't know that there is anything special about tstart. Result: adaptive integration zooms past tstart and someparameter remains unperturbed.

The way to prevent that from happening is to use events to force parameter changes in the course of a simulation. This entry in ModelDB
https://senselab.med.yale.edu/modeldb/S ... odel=39949
shows how to do this with an event-driven mechanism specified with NMODL code.
This file
https://www.neuron.yale.edu/ftp/ted/neuron/ipulse.zip
shows an entirely hoc-based approach for forcing parameter changes that uses an FInitializeHandler and cvode.event.
tardebut
Posts: 21
Joined: Thu Sep 18, 2014 8:23 am

Re: OK with fixed step but not with cvode solvers

Post by tardebut »

Thanks Ted. You are right. I have that piece of code in the neuromodulator .mod.
I will try your suggestions, starting by the second.

Thanks again.

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

Re: OK with fixed step but not with cvode solvers

Post by ted »

The "pure hoc" approach, which uses cvode.event and FInitializeHandler, is best suited for simulation protocols that involve a few parameter changes at predetermined times.
tardebut
Posts: 21
Joined: Thu Sep 18, 2014 8:23 am

Re: OK with fixed step but not with cvode solvers

Post by tardebut »

Hi Ted here you have some partial results.
I simplified the problem to try the pure .hoc solution.
This solution worked. But in this simplified problem the old approach,

Code: Select all

if (t>tstart) {
  someparameter = value1
}
also worked.

This is the .hoc,

Code: Select all

load_file("nrngui.hoc")

objref cvode
cvode = new CVode()
cvode.active(1)
v_init = -75

create soma
access soma
soma {     L = 30 
	diam = 30 
	nseg = 1}

peakstart = 1000
peakdur = 2000

insert peaktraineventg 
startrans_peaktraineventg = peakstart
durtrans_peaktraineventg = peakdur
cvode.maxstep(100)

objref tvec, peak

peak = new Vector()
tvec = new Vector()

tstop = 10000

objref fih
fih = new FInitializeHandler("initi()")		
proc initi() {		
  		transient = 0
  		cvode.event(peakstart, "transientNM()")
	      }

proc transientNM() {
  		    if (transient==0) {
    		    			transient = 1
    		    			cvode.event(t + peakdur, "transientNM()")
					print "A transients was triggered at ", t, " ms" 
  		 			} else {
    						transient = 0
						print "Back to the basal level. The transient ended at ", t, " ms"
  		 			        }
		    
  		     soma {transient_peaktraineventg = transient }
		     cvode.re_init()    					
		    }


objref g 
g = new Graph() 
g.size(0,10000,0,1200)

cvode.record(&peak_peaktraineventg(0.5),peak,tvec)
g.addvar("soma.peak_peaktraineventg(0.5)") 

dt = 10

//init()

//proc go() { print cvode.active() 
//	    run() }



proc initialize() { finitialize(v_init) 
		    fcurrent()}

proc integrate() { g.begin()
		   while (t<tstop) { fadvance() 
				     g.plot(t)} 
		   		     g.flush()}

proc go() { initialize() 
	    integrate()
	    print cvode.active() }
And this is the .mod

Code: Select all

TITLE GlobalNeuromodulatorPeaks

: This is a neuromodulatory transient integrated with an adaptative step solver 

UNITS {
(molar) = (1/liter)
   (nM) = (nanomolar)
      }

NEURON {
  	SUFFIX peaktraineventg
	RANGE k1
        RANGE k2
	RANGE mxamp
	RANGE bsamp
	RANGE startrans
	RANGE durtrans
	RANGE peak
	RANGE transient
	}

PARAMETER {
	k1 = 0.0083 (kilohz)
	k2 = 0.00000333 (kilohz)
	bsamp = 10 (nM)
	mxamp = 1000 (nM)
	startrans = 4000 (ms)
	durtrans = 2000 (ms)
	transient = 0
	  }

ASSIGNED {
	  tt (ms) 
	  tx (ms)
	  K
	  peak (nM)
	  v (millivolt)
	 }

INITIAL { peak = bsamp }

BREAKPOINT {
	tx = log(k1/k2)/(k1-k2)
    	K = 1/(exp(-k1*tx) - exp(-k2*tx))
	tt = t - startrans

   : The event-based approach
	peak = bsamp + transient*mxamp*K*(exp(-k1*tt)-exp(-k2*tt))            

   : The 'old' approach
   :if (t < startrans+durtrans && t > startrans){ peak = bsamp + mxamp*K*(exp(-k1*tt)-exp(-k2*tt))}
   :				                                else { peak = bsamp }
	    }

In the .mod you can comment out the event-based or the 'old' approach and in both cases the results are the same using an adaptative step solver ... unless I'm doing something stupid.

This makes me think that this was not the source of the problem in the original model. Furthermore, when implementing the event-based parameter change in the original model, the neuromodulatory effects are still absent.

I look forward to hear from you.

Thanks in advance,

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

Re: OK with fixed step but not with cvode solvers

Post by ted »

In the .mod you can comment out the event-based or the 'old' approach and in both cases the results are the same using an adaptative step solver
That may have worked in a particular case, but I wouldn't bet a nickel that it will always work. Feelin' lucky? Save it for Vegas.

I'm glad you posted some demo code. It's just enough to demonstrate a couple of tricky problems that otherwise may have remained unrecognized, not to mention undiagnosed.

1. g is not properly initialized. The first plotted point in g is at t = 3.43226 ms. The first plotted point in any graph that has a time axis should be at the start of the simulation, i.e. t = 0 ms when using adaptive integration. It is best to use the NEURON's standard run system to handle model initialization and simulation execution (the latter includes graph initialization and updating). Attempts to replace it are likely to produce results that contain unexpected flaws, if not outright errors. So comment out procs initialize, integrate, and go, and simply call run() when you want to launch a simulation.

2. The plot of peak_peaktraineventg vs. t shows that peak_peak... misses the start of the transient. Nothing happens until t = 1003.44 ms, and the rising phase of the transient then lags 100 ms behind what it should be. At each new advance from t to t+dt, you want peak_peak... to get the value it should have at the new time, that is the value at t+dt. Instead, what's happening is that peak_peak... is being computed from the "old" time t. This happens because the peak... mechanism's model description is purely algebraic (no differential equations are involved); it wouldn't happen if the model description of this variable was specified with a pair of first order differential equations. Also, I see that the BREAKPOINT block contains other statements that really don't need to be evaulated at each advance--just once during initialization would be sufficient. So I suggest changing

Code: Select all

INITIAL {
  peak = bsamp
}
BREAKPOINT {
 . . . stuff . . .
}
to

Code: Select all

INITIAL {
  peak = bsamp
  tx = log(k1/k2)/(k1-k2)
  K = 1/(exp(-k1*tx) - exp(-k2*tx))
}
BEFORE BREAKPOINT {
  tt = t - startrans
  peak = bsamp + transient*mxamp*K*(exp(-k1*tt)-exp(-k2*tt))            
}
COMMENT
BREAKPOINT {
 . . . stuff . . .
}
ENDCOMMENT
and you'll find that peak_peak... now jumps up from the baseline right after the start of the transient, as it should.

Finally, I would caution against using cvode.active() to turn adaptive integration on and off. It's better to use cvode_active(), which is part of the standard run system and takes proper care of things like graph updating--see
viewtopic.php?f=8&t=897&p=10447
tardebut
Posts: 21
Joined: Thu Sep 18, 2014 8:23 am

Re: OK with fixed step but not with cvode solvers

Post by tardebut »

Thanks Ted.
I followed your suggestions. But still, the original problem of this post keeps being without solution at least with the approach based on cvode.event() that is the one I have tried so far.
Here you have an small size example suitable for posting that reproduce the problem.

This is the .mod,

Code: Select all

TITLE FromD1RGolf1Cmp
  
UNITS {
(molar) = (1/liter)
(nM) = (nanomolar)
(snMminus) = (kilohz/nM)
}
 
NEURON {
SUFFIX D1RGolf1Cmp
:EXTERNAL peak_peaktrainglobal
RANGE k1
RANGE k2
RANGE mxamp
RANGE bsamp
RANGE startrans
RANGE durtrans
RANGE peak
RANGE transient
GLOBAL GaolfGTPout
RANGE kGolfback
RANGE kGaolfGTPase
RANGE kfD1R_DA
RANGE krD1R_DA
RANGE kfD1R_Golf
RANGE krD1R_Golf
RANGE kfD1RDA_Golf
RANGE krD1RDA_Golf
RANGE kfD1RGolf_DA
RANGE krD1RGolf_DA
RANGE kactGolf
RANGE betaD1R
RANGE alphaD1R
RANGE KDD1R_DA
}
 
PARAMETER {
transient = 0
k1 = 0.002 (kilohz)
k2 = 0.003 (kilohz)
bsamp = 10 (nM)
mxamp = 1000 (nM)
startrans = 1000 (ms)
durtrans = 2000 (ms)
kGolfback = 0.1 (snMminus)
kGaolfGTPase = 0.001 (kilohz)
kfD1R_DA = 1e-06 (snMminus)
krD1R_DA = 0.002 (kilohz)
kfD1R_Golf = 1e-07 (snMminus)
krD1R_Golf = 0.002 (kilohz)
kfD1RDA_Golf = 1e-05 (snMminus)
krD1RDA_Golf = 0.002 (kilohz)
kfD1RGolf_DA = 0.0001 (snMminus)
krD1RGolf_DA = 0.002 (kilohz)
kactGolf = 0.003 (kilohz)
betaD1R = 0.1 (kilohz)
alphaD1R = 0.0001 (kilohz)
KDD1R_DA = 2 (kilohz)
}
 
ASSIGNED {
tt (ms) 
tx (ms)
K
peak (nM)
:peak_peaktrainglobal (nM)
DA (nM)
GaolfGTPout (nM)
}
 
STATE {
GaolfGTP (nM)
D1R (nM)
D1RDA (nM)
D1RGolf (nM)
Golf (nM)
D1RDAGolf (nM)
Gbgolf (nM)
GaolfGDP (nM)
}
 
INITIAL {
peak = bsamp
tx = log(k1/k2)/(k1-k2)
K = 1/(exp(-k1*tx) - exp(-k2*tx))
GaolfGTP = 30.8104
D1R = 658.9626
D1RDA = 5.6049
D1RGolf = 25.1604
Golf = 833.7484
D1RDAGolf = 10.2701
Gbgolf = 30.8198
GaolfGDP = 0.009997
}
 
BEFORE BREAKPOINT {
	:at_time(startrans)		:This was not required in any case
	:at_time(startrans+durtrans)	:This was not required in any case
	tt = t - startrans
	peak = bsamp + transient*mxamp*K*(exp(-k1*tt)-exp(-k2*tt))
	:
	: Comment the last line above and uncomment the two below to try the fixspet solver
	:
	:if (t < startrans+durtrans && t > startrans){ peak = bsamp + mxamp*K*(exp(-k1*tt)-exp(-k2*tt))}
	:			       			else { peak = bsamp }
		   }


BREAKPOINT {
SOLVE state METHOD derivimplicit
DA = peak
GaolfGTPout = GaolfGTP
}
 
DERIVATIVE state {
GaolfGTP' =  - (kGaolfGTPase*GaolfGTP) + (kactGolf*D1RDAGolf)
D1R' =  - (kfD1R_DA*D1R*DA - krD1R_DA*D1RDA) - (kfD1R_Golf*D1R*Golf - krD1R_Golf*D1RGolf)
D1RDA' = (kfD1R_DA*D1R*DA - krD1R_DA*D1RDA) - (kfD1RDA_Golf*D1RDA*Golf - krD1RDA_Golf*D1RDAGolf) + (kactGolf*D1RDAGolf)
D1RGolf' = (kfD1R_Golf*D1R*Golf - krD1R_Golf*D1RGolf) - (kfD1RGolf_DA*D1RGolf*DA - krD1RGolf_DA*D1RDAGolf)
Golf' = (kGolfback*GaolfGDP*Gbgolf) - (kfD1R_Golf*D1R*Golf - krD1R_Golf*D1RGolf) - (kfD1RDA_Golf*D1RDA*Golf - krD1RDA_Golf*D1RDAGolf)
D1RDAGolf' = (kfD1RDA_Golf*D1RDA*Golf - krD1RDA_Golf*D1RDAGolf) + (kfD1RGolf_DA*D1RGolf*DA - krD1RGolf_DA*D1RDAGolf) - (kactGolf*D1RDAGolf)
Gbgolf' =  - (kGolfback*GaolfGDP*Gbgolf) + (kactGolf*D1RDAGolf)
GaolfGDP' =  - (kGolfback*GaolfGDP*Gbgolf) + (kGaolfGTPase*GaolfGTP)
}
And this is the .hoc,

Code: Select all

load_file("nrngui.hoc")
create soma
access soma
soma {     L = 30 
	diam = 30 
	nseg = 1}

insert D1RGolf1Cmp

startrans = 3000
durtrans = 2000

neuromodulator = 1

transient = 0

k1_D1RGolf1Cmp = 0.002
k2_D1RGolf1Cmp = 0.003
bsamp_D1RGolf1Cmp = 10
mxamp_D1RGolf1Cmp = 1000
startrans_D1RGolf1Cmp = startrans
durtrans_D1RGolf1Cmp = durtrans

kGolfback_D1RGolf1Cmp = 0.1
kGaolfGTPase_D1RGolf1Cmp = 0.001
kfD1R_DA_D1RGolf1Cmp = 1e-06
krD1R_DA_D1RGolf1Cmp = 0.002
kfD1R_Golf_D1RGolf1Cmp = 1e-07
krD1R_Golf_D1RGolf1Cmp = 0.002
kfD1RDA_Golf_D1RGolf1Cmp = 1e-05
krD1RDA_Golf_D1RGolf1Cmp = 0.002
kfD1RGolf_DA_D1RGolf1Cmp = 0.0001
krD1RGolf_DA_D1RGolf1Cmp = 0.002
kactGolf_D1RGolf1Cmp = 0.003
betaD1R_D1RGolf1Cmp = 0.1
alphaD1R_D1RGolf1Cmp = 0.0001
KDD1R_DA_D1RGolf1Cmp = 2

//*

// This is for the adaptative step 

objref cvode
cvode = new CVode()
cvode_active(1)
cvode.maxstep(100)

objref fih
fih = new FInitializeHandler("initi()")		
proc initi() {		
  		transient = 0
  		cvode.event(startrans, "transientNM()")
	     }

proc transientNM() {
  		    if ((transient==0)&&neuromodulator) {
    					 		transient = 1
    					 		cvode.event(t + durtrans, "transientNM()")
					 		print "A neuromodulatory transient was triggered at ", t, " ms"
  		 		       			} else {
    								transient = 0
								print "Back to the basal level. The neuromodulatory transient ended at ", t, " ms"
  		 			       			}

			transient_D1RGolf1Cmp = transient

			if (cvode_active()) { cvode.re_init() }    					
			} 

//*/


///////////////////////////// Plotting Results /////////////////////////////

objref g 
g = new Graph() 
g.size(0,10000,0,1200)
graphList[0].append(g)
g.addvar("GaolfGTPout_D1RGolf1Cmp")
g.addvar("peak_D1RGolf1Cmp(0.5)")

///////////////////////////* Simulation Control *///////////////////////////

dt = 100 
tstop = 10000 
v_init = -65

init()

proc go() {print "Neuromodulation is ", neuromodulator
	      print "Adaptative solver is ", cvode_active()
	      run()}

Either with the fixed or the adaptative step methods the algebraic equation "peak = ..." is computed.
However, just with the fixed step method the ODEs are solved.

I look forward to hear from you or from any other who can give me a hint on how to solve this.

Thanks in advance,

Omar
tardebut
Posts: 21
Joined: Thu Sep 18, 2014 8:23 am

Re: OK with fixed step but not with cvode solvers

Post by tardebut »

Hi there again!

Using the last example, I tried the approach with NET_RECEIVE and auto-events (net_send). First, the previous distributed mechanism was converted into a point process which is the one compatible with NET_RECEIVE. The post http://www.neuron.yale.edu/phpBB/viewto ... alue#p8937 was of some help here.

The results are the same than with cvode.event() in the distributed mechanism. The algebraic equation is computed with both fixed and adaptative methods but the adaptative does not integrates the ODEs while the fixed does. However, If the algebraic equations are placed inside the BREAKPOINT block no computation is made with the adaptative step solver, while the fixed step solver keeps computing the algebraic equation and integrating the ODEs just fine. The BREAKPOINT block seems plain dead for the adaptative solver with this problem. By the way, all of this keeps the same disregarding whether I use METHOD derivimplicit or cnexp.

Am I missing something obvious here?

Here you have the .mod,

Code: Select all

TITLE FromD1RGolf1Cmp
  
UNITS {
(molar) = (1/liter)
(nM) = (nanomolar)
(snMminus) = (kilohz/nM)
}
 
NEURON {
POINT_PROCESS D1RGolf1Cmp
RANGE k1
RANGE k2
RANGE mxamp
RANGE bsamp
RANGE startrans
RANGE durtrans
RANGE peak
RANGE transient
RANGE GaolfGTPout		:GLOBAL
RANGE kGolfback
RANGE kGaolfGTPase
RANGE kfD1R_DA
RANGE krD1R_DA
RANGE kfD1R_Golf
RANGE krD1R_Golf
RANGE kfD1RDA_Golf
RANGE krD1RDA_Golf
RANGE kfD1RGolf_DA
RANGE krD1RGolf_DA
RANGE kactGolf
RANGE betaD1R
RANGE alphaD1R
RANGE KDD1R_DA
}
 
PARAMETER {
transient = 0
w = 0
k1 = 0.002 (kilohz)
k2 = 0.003 (kilohz)
bsamp = 10 (nM)
mxamp = 1000 (nM)
startrans = 1000 (ms)
durtrans = 2000 (ms)
kGolfback = 0.1 (snMminus)
kGaolfGTPase = 0.001 (kilohz)
kfD1R_DA = 1e-06 (snMminus)
krD1R_DA = 0.002 (kilohz)
kfD1R_Golf = 1e-07 (snMminus)
krD1R_Golf = 0.002 (kilohz)
kfD1RDA_Golf = 1e-05 (snMminus)
krD1RDA_Golf = 0.002 (kilohz)
kfD1RGolf_DA = 0.0001 (snMminus)
krD1RGolf_DA = 0.002 (kilohz)
kactGolf = 0.003 (kilohz)
betaD1R = 0.1 (kilohz)
alphaD1R = 0.0001 (kilohz)
KDD1R_DA = 2 (kilohz)
}
 
ASSIGNED {
tt (ms) 
tx (ms)
K
peak (nM)
DA (nM)
GaolfGTPout (nM)
}
 
STATE {
GaolfGTP (nM)
D1R (nM)
D1RDA (nM)
D1RGolf (nM)
Golf (nM)
D1RDAGolf (nM)
Gbgolf (nM)
GaolfGDP (nM)
}
 
INITIAL {
net_send(startrans,1)
peak = bsamp
tx = log(k1/k2)/(k1-k2)
K = 1/(exp(-k1*tx) - exp(-k2*tx))

GaolfGTP = 30.8104
D1R = 658.9626
D1RDA = 5.6049
D1RGolf = 25.1604
Golf = 833.7484
D1RDAGolf = 10.2701
Gbgolf = 30.8198
GaolfGDP = 0.009997
:GaolfGTPout = 30.8104
}


BEFORE BREAKPOINT {
	 tt = t - startrans
	 peak = bsamp + transient*mxamp*K*(exp(-k1*tt)-exp(-k2*tt))	
	
	:Comment the last line above and uncomment the two below to try the fix step solver
	
	:if (t < startrans+durtrans && t > startrans){ peak = bsamp + mxamp*K*(exp(-k1*tt)-exp(-k2*tt))}
	:			       			else { peak = bsamp }
		   }


BREAKPOINT {
SOLVE state METHOD derivimplicit
DA = peak
GaolfGTPout = GaolfGTP
}


DERIVATIVE state {
GaolfGTP' =  - (kGaolfGTPase*GaolfGTP) + (kactGolf*D1RDAGolf)
D1R' =  - (kfD1R_DA*D1R*DA - krD1R_DA*D1RDA) - (kfD1R_Golf*D1R*Golf - krD1R_Golf*D1RGolf)
D1RDA' = (kfD1R_DA*D1R*DA - krD1R_DA*D1RDA) - (kfD1RDA_Golf*D1RDA*Golf - krD1RDA_Golf*D1RDAGolf) + (kactGolf*D1RDAGolf)
D1RGolf' = (kfD1R_Golf*D1R*Golf - krD1R_Golf*D1RGolf) - (kfD1RGolf_DA*D1RGolf*DA - krD1RGolf_DA*D1RDAGolf)
Golf' = (kGolfback*GaolfGDP*Gbgolf) - (kfD1R_Golf*D1R*Golf - krD1R_Golf*D1RGolf) - (kfD1RDA_Golf*D1RDA*Golf - krD1RDA_Golf*D1RDAGolf)
D1RDAGolf' = (kfD1RDA_Golf*D1RDA*Golf - krD1RDA_Golf*D1RDAGolf) + (kfD1RGolf_DA*D1RGolf*DA - krD1RGolf_DA*D1RDAGolf) - (kactGolf*D1RDAGolf)
Gbgolf' =  - (kGolfback*GaolfGDP*Gbgolf) + (kactGolf*D1RDAGolf)
GaolfGDP' =  - (kGolfback*GaolfGDP*Gbgolf) + (kGaolfGTPase*GaolfGTP)
}


NET_RECEIVE (w) {
  	     	if (flag==1) {
    			      transient = 1
			      net_send(durtrans,2)
  			      }

	     	if (flag==2) { 
    			      transient = 0
  			      }

	       }
and the .hoc,

Code: Select all

load_file("nrngui.hoc")
create soma
access soma
soma {     L = 30 
	diam = 30 
	nseg = 1}

objref injneurmod
soma injneurmod = new D1RGolf1Cmp(0.5)

startrans = 1000
durtrans = 2000

neuromodulator = 1

injneurmod.k1 = 0.002
injneurmod.k2 = 0.003
injneurmod.bsamp = 10
injneurmod.mxamp = 1000
injneurmod.startrans = startrans
injneurmod.durtrans = durtrans

injneurmod.kGolfback = 0.1
injneurmod.kGaolfGTPase = 0.001
injneurmod.kfD1R_DA = 1e-06
injneurmod.krD1R_DA = 0.002
injneurmod.kfD1R_Golf = 1e-07
injneurmod.krD1R_Golf = 0.002
injneurmod.kfD1RDA_Golf = 1e-05
injneurmod.krD1RDA_Golf = 0.002
injneurmod.kfD1RGolf_DA = 0.0001
injneurmod.krD1RGolf_DA = 0.002
injneurmod.kactGolf = 0.003
injneurmod.betaD1R = 0.1
injneurmod.alphaD1R = 0.0001
injneurmod.KDD1R_DA = 2

// This is for the adaptative step

/*

objref cvode
cvode = new CVode()
cvode_active(1)
cvode.maxstep(10)

*/


///////////////////////////// Plotting Results /////////////////////////////

objref g 
g = new Graph() 
g.size(0,10000,0,1200)
graphList[1].append(g)
g.addvar("injneurmod.GaolfGTPout")
g.addvar("injneurmod.peak")

///////////////////////////* Simulation Control *///////////////////////////

dt = 100 
tstop = 10000 
v_init = -65

init()

proc go() {print "Neuromodulation is ", neuromodulator
	   print "Adaptative solver is ", cvode_active()
	   run()}
I look forward.

Regards,

Omar
tardebut
Posts: 21
Joined: Thu Sep 18, 2014 8:23 am

Re: OK with fixed step but not with cvode solvers

Post by tardebut »

Sorry, the BEFORE BREAKPOINT in the .mod below is as follows,

Code: Select all

BEFORE BREAKPOINT {
	 tt = t - startrans
	 peak = bsamp + transient*mxamp*K*(exp(-k1*tt)-exp(-k2*tt))	
		   }
The rest, code and behavior, is the same.

Omar
tardebut
Posts: 21
Joined: Thu Sep 18, 2014 8:23 am

Re: OK with fixed step but not with cvode solvers

Post by tardebut »

After some critical assistance, it become clear that the solution to this problem is to move DA = peak to the beginning of the DERIVATIVE block and GaolfGTPout = GaolfGTP to the end of the DERIVATIVE block.
This is because these are not currents and the extension of BREAKPOINT to work with variables that are not currents is limited to fixed step methods.

Thanks,

Omar
Post Reply