OK with fixed step but not with cvode solvers
OK with fixed step but not with cvode solvers
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
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

 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
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: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 eventdriven mechanism specified with NMODL code.
This file
https://www.neuron.yale.edu/ftp/ted/neuron/ipulse.zip
shows an entirely hocbased approach for forcing parameter changes that uses an FInitializeHandler and cvode.event.
Code: Select all
if (t>tstart) {
someparameter = value1
}
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 eventdriven mechanism specified with NMODL code.
This file
https://www.neuron.yale.edu/ftp/ted/neuron/ipulse.zip
shows an entirely hocbased approach for forcing parameter changes that uses an FInitializeHandler and cvode.event.
Re: OK with fixed step but not with cvode solvers
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
I will try your suggestions, starting by the second.
Thanks again.
Omar

 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
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.
Re: OK with fixed step but not with cvode solvers
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,
also worked.
This is the .hoc,
And this is the .mod
In the .mod you can comment out the eventbased 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 eventbased parameter change in the original model, the neuromodulatory effects are still absent.
I look forward to hear from you.
Thanks in advance,
Omar
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
}
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() }
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)/(k1k2)
K = 1/(exp(k1*tx)  exp(k2*tx))
tt = t  startrans
: The eventbased 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 }
}
This makes me think that this was not the source of the problem in the original model. Furthermore, when implementing the eventbased parameter change in the original model, the neuromodulatory effects are still absent.
I look forward to hear from you.
Thanks in advance,
Omar

 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
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.In the .mod you can comment out the eventbased or the 'old' approach and in both cases the results are the same using an adaptative step solver
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 advancejust once during initialization would be sufficient. So I suggest changing
Code: Select all
INITIAL {
peak = bsamp
}
BREAKPOINT {
. . . stuff . . .
}
Code: Select all
INITIAL {
peak = bsamp
tx = log(k1/k2)/(k1k2)
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
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 updatingsee
viewtopic.php?f=8&t=897&p=10447
Re: OK with fixed step but not with cvode solvers
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,
And this is the .hoc,
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
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 = 1e06 (snMminus)
krD1R_DA = 0.002 (kilohz)
kfD1R_Golf = 1e07 (snMminus)
krD1R_Golf = 0.002 (kilohz)
kfD1RDA_Golf = 1e05 (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)/(k1k2)
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)
}
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 = 1e06
krD1R_DA_D1RGolf1Cmp = 0.002
kfD1R_Golf_D1RGolf1Cmp = 1e07
krD1R_Golf_D1RGolf1Cmp = 0.002
kfD1RDA_Golf_D1RGolf1Cmp = 1e05
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()}
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
Re: OK with fixed step but not with cvode solvers
Hi there again!
Using the last example, I tried the approach with NET_RECEIVE and autoevents (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,
and the .hoc,
I look forward.
Regards,
Omar
Using the last example, I tried the approach with NET_RECEIVE and autoevents (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 = 1e06 (snMminus)
krD1R_DA = 0.002 (kilohz)
kfD1R_Golf = 1e07 (snMminus)
krD1R_Golf = 0.002 (kilohz)
kfD1RDA_Golf = 1e05 (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)/(k1k2)
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
}
}
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 = 1e06
injneurmod.krD1R_DA = 0.002
injneurmod.kfD1R_Golf = 1e07
injneurmod.krD1R_Golf = 0.002
injneurmod.kfD1RDA_Golf = 1e05
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()}
Regards,
Omar
Re: OK with fixed step but not with cvode solvers
Sorry, the BEFORE BREAKPOINT in the .mod below is as follows,
The rest, code and behavior, is the same.
Omar
Code: Select all
BEFORE BREAKPOINT {
tt = t  startrans
peak = bsamp + transient*mxamp*K*(exp(k1*tt)exp(k2*tt))
}
Omar
Re: OK with fixed step but not with cvode solvers
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
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