recording the second derivative of an action potential
-
- Posts: 49
- Joined: Mon Mar 25, 2013 1:36 pm
- Location: france
- Contact:
recording the second derivative of an action potential
Hello,
I'd like to record from the second derivative (y") of an action potential and make graphs : y''=f(t) and a phase plot y''= f(x') = f(i_cap).
If I well understand y'' is the derivative of I-cap against t.
Should I create a mod file that read i or i_cap and the derive it?
Thanks a lot
Fabien
I'd like to record from the second derivative (y") of an action potential and make graphs : y''=f(t) and a phase plot y''= f(x') = f(i_cap).
If I well understand y'' is the derivative of I-cap against t.
Should I create a mod file that read i or i_cap and the derive it?
Thanks a lot
Fabien
-
- Site Admin
- Posts: 6299
- Joined: Wed May 18, 2005 4:50 pm
- Location: Yale University School of Medicine
- Contact:
Re: recording the second derivative of an action potential
No. Record the time course of membrane potential and i_cap to vectors during the simulation. After the simulation, use the Vector class's methods in your own algorithms to calculate the second derivative of the spike and generate the plots of interest. Your results will be most accurate if you use adaptive integration for the simulation. Yes, that produces results that are spaced unevenly in time, but you can use the Vector class's interpolate() method to resample the recorded values at evenly spaced intervals if you need to.fabien tell wrote:I'd like to record from the second derivative (y") of an action potential and make graphs : y''=f(t) and a phase plot y''= f(x') = f(i_cap).
If I well understand y'' is the derivative of I-cap against t.
Should I create a mod file that read i or i_cap and the derive it?
-
- Posts: 49
- Joined: Mon Mar 25, 2013 1:36 pm
- Location: france
- Contact:
Re: recording the second derivative of an action potential
Hello Ted,
Thanks. You mean that I can implement it in the hoc file that launches my simulation ?
I didn’t know we could use derivative functions here as in modl files.
Have a good day
Fabien
Thanks. You mean that I can implement it in the hoc file that launches my simulation ?
I didn’t know we could use derivative functions here as in modl files.
Have a good day
Fabien
-
- Site Admin
- Posts: 6299
- Joined: Wed May 18, 2005 4:50 pm
- Location: Yale University School of Medicine
- Contact:
Re: recording the second derivative of an action potential
Organize your program in this sequence.
1. Model specification code--anatomical and biophysical properties that are represented in the model
2. Instrumentation code--stimulus sources, graphical displays of variables during the simulation, recording of variables (e.g. v, spike times, calcium current) for post-simulation analysis
3. Simulation control code--custom initialization code, specification of simulation parameters (v_init, dt, tstop), simulation flow control e.g. "protocols" for changing parameters during a simulation (e.g. concentrations of ionic species or ligands)(implemented with FInitializeHandler, cvode.event).
Make sure that you can start a simulation by calling run(). When all is working properly, you're ready for the next step.
4. Post-simulation analysis code--create a procedure that is to be called after the simulation ends. This procedure, which you might want to name "proc postproc()" calls other procs and funcs that you wrote which operate on recorded variables to derive things such as derivatives, integrals, spike frequency etc.. It will look something like this:
Verify that postproc() works correctly.
Finally the last step:
5. Make a custom "run" procedure that launches a simulation and then performs the analysis.
And now you have automated launching a simulation and analyzing results afterwards.
1. Model specification code--anatomical and biophysical properties that are represented in the model
2. Instrumentation code--stimulus sources, graphical displays of variables during the simulation, recording of variables (e.g. v, spike times, calcium current) for post-simulation analysis
3. Simulation control code--custom initialization code, specification of simulation parameters (v_init, dt, tstop), simulation flow control e.g. "protocols" for changing parameters during a simulation (e.g. concentrations of ionic species or ligands)(implemented with FInitializeHandler, cvode.event).
Make sure that you can start a simulation by calling run(). When all is working properly, you're ready for the next step.
4. Post-simulation analysis code--create a procedure that is to be called after the simulation ends. This procedure, which you might want to name "proc postproc()" calls other procs and funcs that you wrote which operate on recorded variables to derive things such as derivatives, integrals, spike frequency etc.. It will look something like this:
Code: Select all
// assumes that
// tvec is a Vector to which time was recorded
// i_ca_soma is a Vector to which somatic i_ca was recorded, and that the soma has only 1 compartment
// netcasoma = 0 // will be total amount of calcium that entered the soma
proc postproc() {
netca = integr(i_ca_soma, tvec, t0, t1) // user-written func integr() returns integral of $o1 over interval t0,t1
// netca is now in ms*mA/cm2, but we want net charge
soma netca *= area(0.5)*(1e-8) // 1e-8 because area() returns area in um2
// netca is now in ms*mA = microcoulombs
. . . whatever else you need to do . . .
}
Finally the last step:
5. Make a custom "run" procedure that launches a simulation and then performs the analysis.
Code: Select all
proc myrun() {
run()
postproc()
}
I don't think I've seen a mod file that calculates a derivative. Wouldn't want to do that anyway. Better to record the variable of interest to a vector, then calculate the derivative after the simulation, because the code already exists--read about the Vector class's deriv() function--and it's easy to reuse that with just a little bit of user written code to do things like control the window over which integration is performed, resample unevenly spaced values at a fixed sampling rate (if the simulation used adaptive integration), or create a user-written hoc function that employs an integration algorithm designed to handle unevenly spaced data.I didn’t know we could use derivative functions here as in modl files.
-
- Posts: 49
- Joined: Mon Mar 25, 2013 1:36 pm
- Location: france
- Contact:
Re: recording the second derivative of an action potential
Thanks a lot Ted
Very helpful !!!
Very helpful !!!
-
- Posts: 49
- Joined: Mon Mar 25, 2013 1:36 pm
- Location: france
- Contact:
Re: recording the second derivative of an action potential
Hello,
I've been trying to obtain the derivative of v at soma uisng the vectors but something is wrong since the terminal "tells" me that it is not possible to obtain a derivative with less than two points which is true.
I get the following message :
I've been trying to obtain the derivative of v at soma uisng the vectors but something is wrong since the terminal "tells" me that it is not possible to obtain a derivative with less than two points which is true.
Code: Select all
create soma
soma {
nseg=1
L=35
diam=15
Ra = 150.0
cm=1
insert pasnts
insert kdrDA
gbar_kdrDA=100
insert Na16
gbar_Na16=300
insert ka
gbar_ka=2
ena=60
ek=-90
}
addgraph("soma.v(0.5)",-100,20)
addgraph("soma.i_cap(0.5)",-0.2,0.2)
/* if i'm right I below created 3 vectors, the "spike" contains the voltage time course in the soma and the "deriv_1" should contain the first derivative of it with a dt fixed at 0.025, and the deriv_2 should contain the 2nd deriv.
I 'm surely wrong somewhere */
objref spike, deriv_1, deriv_2
spike=new Vector()
deriv_1=new Vector()
deriv_2=new Vector()
spike.record(&soma.v(0.5))
deriv_1=deriv_1.deriv(spike, 0.025)
deriv_2= deriv_2.deriv(deriv_1,0.025)
run()
Thanks a lot/usr/local/nrn/x86_64/bin/nrniv: Can't take derivative of Vector with less than two points
in soma seul.hoc near line 79
deriv_1=deriv_1.deriv(spike, 0.025)
^
Vector[4].deriv(Vector[3], 0.025)
-
- Site Admin
- Posts: 6299
- Joined: Wed May 18, 2005 4:50 pm
- Location: Yale University School of Medicine
- Contact:
Re: recording the second derivative of an action potential
The vector that is to hold the time course of v is empty until after a simulation has been executed.
Defer execution of Vector.derivative until after run() has been called.
Defer execution of Vector.derivative until after run() has been called.
-
- Posts: 49
- Joined: Mon Mar 25, 2013 1:36 pm
- Location: france
- Contact:
Re: recording the second derivative of an action potential
Thanks again Ted,
I therefore wrote this :
It works perfectly but I have a last question. If I change the simulation parameters (e.g sodium conductance or soma size, ), the graph doesn’t change as if the vector had not been updated by a new run . If I erase it , I cannot plot it any more.
Thanks again
I therefore wrote this :
Code: Select all
load_file("nrngui.hoc")
tstop=50
dt=0.025
celsius=32
v_init=-55
finitialize(v_init)
create soma
soma {
nseg=1
L=35
diam=15
Ra = 150.0
cm=1
insert pasnts
insert kdrDA
gbar_kdrDA=100
insert Na12
Vmid_ac_Na12 = -35
gbar_Na12=250
insert ka
gbar_ka=2
ena=60
ek=-90
}
addgraph("soma.v(0.5)",-100,20)
addgraph("soma.i_cap(0.5)",-0.2,0.2)
run()
objref spike, deriv_1, deriv_2, g
spike=new Vector()
deriv_1=new Vector()
deriv_2=new Vector()
g=new Graph()
spike.record(&soma.v(0.5))
deriv_1.record(&soma.i_cap(0.5))
run()
deriv_2= deriv_2.deriv(deriv_1,0.025)
for i=0, 0.1 {deriv_2.line(g,deriv_1) deriv_2.rotate(1) }
Thanks again
-
- Site Admin
- Posts: 6299
- Joined: Wed May 18, 2005 4:50 pm
- Location: Yale University School of Medicine
- Contact:
Re: recording the second derivative of an action potential
It's looking good. Here are some suggestions.
1. The call to finitialize does nothing useful and can be eliminated. run() executes finitialize before it does anything else.
2. The first call to run() can be commented out--it isn't doing anything useful.
execute a simulation
calculate derivative(s)
plot results
I'm going to assume that you're capturing the time course of v and i_cap to a pair of Vectors called vvec and icvec, e.g. by doing thisin the "instrumentation" section of your program (that's after the model's anatomical and biophysical properties have been specified, and before run() is called).
It is convenient to call the procedure that automates simulation followed by data analysis and display something like myrun. myrun() is very simpleIt makes sense to put proc myrun() right after the instrumentation part of your program, and maybe precede it with the commentbecause that's what it does--it controls the sequence of events that happen during and after a simulation.
The role of proc postprocess() is to calculate the derivative of i_cap and then plot it vs. i_cap. Here's what I'd do--You can put this code right after the definition of proc myrun().
Now if you change parameters, just execute
myrun()
and g2 will automatically be updated.
You might be interested in plotting i_cap vs. v in the course of a simulation.
To do that, insert this at the end of your instrumentation code, i.e. just before the start of the simulation flow control codeIf you want this graph to rescale automatically, include the statement
icvsv.exec_menu("View = plot")
as the last thing that proc postprocess() does.
1. The call to finitialize does nothing useful and can be eliminated. run() executes finitialize before it does anything else.
2. The first call to run() can be commented out--it isn't doing anything useful.
Instead of calling run and calculating the derivatives at the top level of the interpreter, wrap that up in a procedure that automates the sequenceIf I change the simulation parameters (e.g sodium conductance or soma size, ), the graph doesn’t change as if the vector had not been updated by a new run
execute a simulation
calculate derivative(s)
plot results
I'm going to assume that you're capturing the time course of v and i_cap to a pair of Vectors called vvec and icvec, e.g. by doing this
Code: Select all
objref vvec, icvec
vvec = new Vector()
soma vvec.record(&v(0.5))
icvec = new Vector()
soma icvec.record(&i_cap(0.5))
It is convenient to call the procedure that automates simulation followed by data analysis and display something like myrun. myrun() is very simple
Code: Select all
proc myrun() {
run()
postprocess() // analyze and display results
}
Code: Select all
// simulation flow control
The role of proc postprocess() is to calculate the derivative of i_cap and then plot it vs. i_cap. Here's what I'd do--
Code: Select all
objref d2, g2 // d2 will contain derivative of i_cap, g2 will plot that vs. i_cap
d2 = new Vector()
objref nil
proc postprocess() {
d2.deriv(icvec, dt) // note use of dt, so you can use any dt>0
// 0.025 ms produces a rather crude phase plane plot
// instead of making dt smaller, make steps_per_ms larger
// and dt will automatically be made smaller
// 200 works nicely (dt will be 0.005 ms)
if (g2==nil) g2 = new Graph(0) // prevent discarding of previously "kept" traces
d2.plot(g2, icvec)
g2.size(-0.5,0.5,-15,15)
g2.view(-0.5, -15, 1, 30, 648, 296, 300.48, 200.32)
g2.exec_menu("View = plot") // just to make sure it's all visible
g2.label(0.05, 0.9, "x-axis: i_cap(0.5)", 2, 1, 0, 0, 1)
g2.label(0.05, 0.75, "y-axis: D_t(i_cap(0.5))", 2, 1, 0, 0, 1)
// maybe use g2.exec_menu("View = plot") to autoscale the graph
// g2.exec_menu("View = plot")
// to autoscale g2's axes
}
Now if you change parameters, just execute
myrun()
and g2 will automatically be updated.
You might be interested in plotting i_cap vs. v in the course of a simulation.
To do that, insert this at the end of your instrumentation code, i.e. just before the start of the simulation flow control code
Code: Select all
// plot i_cap vs. v during a simulation
objref icvsv
icvsv = new Graph(0)
icvsv.size(-100,20,-0.5,0.5)
icvsv.view(-100, -0.5, 120, 1, 648, 34, 300.48, 200.32)
graphList[3].append(icvsv) // to update it during a run
icvsv.xexpr("v(0.5)", 0)
icvsv.label(0.05, 0.9, "x-axis: v(0.5)", 2, 1, 0, 0, 1)
icvsv.addvar("soma.i_cap( 0.5 )", 2, 1, 0.05, 0.85, 2)
icvsv.exec_menu("View = plot")
as the last thing that proc postprocess() does.
-
- Posts: 49
- Joined: Mon Mar 25, 2013 1:36 pm
- Location: france
- Contact:
Re: recording the second derivative of an action potential
Thanks a lot Ted,
I 'll try it after the long week end and I will post the (commented) code for others.
Have a good Easter day. Go easy on chocolate !!!
Regards
fabien
I 'll try it after the long week end and I will post the (commented) code for others.
Have a good Easter day. Go easy on chocolate !!!
Regards
fabien
-
- Posts: 49
- Joined: Mon Mar 25, 2013 1:36 pm
- Location: france
- Contact:
Re: recording the second derivative of an action potential
Hello Ted,
I tried your solution. It works well. The only problem is that when I run "myrun ()" a second graph (the graph with the second derivative) is put on the screen with the same trace than the first one (which is blank before my run) then if I changed some parameters and do run plus my run the new results are superimposed on the old ones in both graphs plus a new third graph . If I redo it again myrun, I got a fourth graph with the three traces superimposed and so on.
Is it possible to clean graph and to have only one graph (unless I mark keep lines if I want to compare outcomes) ?
Thanks
I tried your solution. It works well. The only problem is that when I run "myrun ()" a second graph (the graph with the second derivative) is put on the screen with the same trace than the first one (which is blank before my run) then if I changed some parameters and do run plus my run the new results are superimposed on the old ones in both graphs plus a new third graph . If I redo it again myrun, I got a fourth graph with the three traces superimposed and so on.
Is it possible to clean graph and to have only one graph (unless I mark keep lines if I want to compare outcomes) ?
Thanks
-
- Site Admin
- Posts: 6299
- Joined: Wed May 18, 2005 4:50 pm
- Location: Yale University School of Medicine
- Contact:
Re: recording the second derivative of an action potential
Aha--a bug! The second run actually damages the first trace--the x coordinates of the first trace's plotted points are unaffected, but the y coordinates of those points are replaced by the y coordinates from the second run. Very bad.
The easy fix is to erase the trace from the previous run before plotting the results of the new run. In proc postprocess() change
toNow there is no persistence from one myrun() to the next.
Pick Vector
Next click on the displayed trace; it should turn red until you release the mouse button, but you might have to move the cursor a short distance along the graph to make this happen (the pointer must come close to a plotted data point).
Now Click on
NEURON Main Menu / Vector / Display
then in this new VecWrap window click on the menu box and select "Copy from Clipboard", then rescale the axes with either "View = plot" or "Set View".
You can paste as many new data traces as you like into a VecWrap, and you can have as many VecWraps as you like. It is not possible to eliminate individual traces in a VecWrap; Delete does remove individual displayed traces, but the undelying data seem to persist--copy a new trace from the clipboard into a VecWrap, and all of the old traces reappear, including any that were Erased or Deleted.
Using Lists to manage collections of Vectors is one way to have more control over what is plotted in a graph. For example, assuming that paramvec contains a sequence of parameter values that you want to try, then you could do something like this
The easy fix is to erase the trace from the previous run before plotting the results of the new run. In proc postprocess() change
Code: Select all
if (g2==nil) g2 = new Graph(0)
Code: Select all
// if (g2==nil) g2 = new Graph(0)
if (g2==nil) {
g2 = new Graph(0)
} else {
g2.exec_menu("Erase")
}
"Keep Lines" only affects traces that are managed by the standard run system. The traces in this particular graph are plotted by user-written code that lies outside of the standard run system. That's why previously generated traces persist (although damaged) unless explicitly erased by other user-written code. It is possible to write code that emulates the Keep Lines feature, but for initial exploratory simulations it is quicker and easier to just take advantage of NEURON's "Vector display" feature. In the phase plane graph click on the menu box and select(unless I mark keep lines if I want to compare outcomes) ?
Pick Vector
Next click on the displayed trace; it should turn red until you release the mouse button, but you might have to move the cursor a short distance along the graph to make this happen (the pointer must come close to a plotted data point).
Now Click on
NEURON Main Menu / Vector / Display
then in this new VecWrap window click on the menu box and select "Copy from Clipboard", then rescale the axes with either "View = plot" or "Set View".
You can paste as many new data traces as you like into a VecWrap, and you can have as many VecWraps as you like. It is not possible to eliminate individual traces in a VecWrap; Delete does remove individual displayed traces, but the undelying data seem to persist--copy a new trace from the clipboard into a VecWrap, and all of the old traces reappear, including any that were Erased or Deleted.
Using Lists to manage collections of Vectors is one way to have more control over what is plotted in a graph. For example, assuming that paramvec contains a sequence of parameter values that you want to try, then you could do something like this
Code: Select all
for i=0,paramvec.size()-1 {
paramname = paramvec.x[i] // replace paramname with the name of the parameter
myrun()
icveclist.append(icvec.c) // append a copy of icvec to icveclist
d2list.append(d2.c) // append a copy of d2 to d2list
}
plotem(icveclist, d2list, g2)
// $o1 list of x vectors
// $o2 list of y vectors
// $o3 graph in which these are to be plotted
proc plotem() { local i
for i=0,$o1.count()-1 $o2.o(i).plot($o3, $o1.o(i))
}
-
- Posts: 49
- Joined: Mon Mar 25, 2013 1:36 pm
- Location: france
- Contact:
Re: recording the second derivative of an action potential
Thanks again TED,
Actually, it's even easier if I simplify the code
For every sequence init and run followed by myrun(), I get a new fresh updated graph . Than if I want to keep traces I copy them using the tip given by Ted.
Thanks again Ted.
Regards
Actually, it's even easier if I simplify the code
Code: Select all
objref d2, g2 // d2 will contain derivative of i_cap, g2 will plot that vs. i_cap
d2 = new Vector()
// objref nil
proc postprocess() {
d2.deriv(icvec, dt)
//if (g2==nil) {
g2 = new Graph(0)
//} else {
//g2.exec_menu("Erase")
// }
d2.plot(g2, icvec) // plot
g2.size(-0.5,0.5,-15,15)
g2.view(-0.5, -15, 1, 30, 648, 296, 400.48, 300.32)
//g2.exec_menu("View = plot") // just to make sure it's all visible
g2.label(0.05, 0.9, "x: i_cap(0.5)", 2, 1, 0, 0, 2)
g2.label(0.05, 0.75, "y: D_t(i_cap(0.5))", 2, 1, 0, 0, 3)
g2.exec_menu("View = plot") // to autoscale the graph
icvsv.exec_menu("View = plot")
}
Thanks again Ted.
Regards
-
- Site Admin
- Posts: 6299
- Joined: Wed May 18, 2005 4:50 pm
- Location: Yale University School of Medicine
- Contact:
Re: recording the second derivative of an action potential
Good. Simpler is better. Easier to understand, easier to revise.