recording the second derivative of an action potential

NMODL and the Channel Builder.
Post Reply
fabien tell
Posts: 49
Joined: Mon Mar 25, 2013 1:36 pm
Location: france
Contact:

recording the second derivative of an action potential

Post by fabien tell »

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
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: recording the second derivative of an action potential

Post by ted »

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?
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
Posts: 49
Joined: Mon Mar 25, 2013 1:36 pm
Location: france
Contact:

Re: recording the second derivative of an action potential

Post by fabien tell »

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
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: recording the second derivative of an action potential

Post by ted »

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:

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 . . .
}
Verify that postproc() works correctly.

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()
}
And now you have automated launching a simulation and analyzing results afterwards.
I didn’t know we could use derivative functions here as in modl files.
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.
fabien tell
Posts: 49
Joined: Mon Mar 25, 2013 1:36 pm
Location: france
Contact:

Re: recording the second derivative of an action potential

Post by fabien tell »

Thanks a lot Ted

Very helpful !!!
fabien tell
Posts: 49
Joined: Mon Mar 25, 2013 1:36 pm
Location: france
Contact:

Re: recording the second derivative of an action potential

Post by fabien tell »

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.

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()
I get the following message :
/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)
Thanks a lot
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: recording the second derivative of an action potential

Post by ted »

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.
fabien tell
Posts: 49
Joined: Mon Mar 25, 2013 1:36 pm
Location: france
Contact:

Re: recording the second derivative of an action potential

Post by fabien tell »

Thanks again Ted,

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) }
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
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: recording the second derivative of an action potential

Post by ted »

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.
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
Instead of calling run and calculating the derivatives at the top level of the interpreter, wrap that up in a procedure that automates the sequence
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))
in 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 simple

Code: Select all

proc myrun() {
  run()
  postprocess() // analyze and display results
}
It makes sense to put proc myrun() right after the instrumentation part of your program, and maybe precede it with the comment

Code: Select all

// simulation flow control
because 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--

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
}
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 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)
If you want this graph to rescale automatically, include the statement
icvsv.exec_menu("View = plot")
as the last thing that proc postprocess() does.
fabien tell
Posts: 49
Joined: Mon Mar 25, 2013 1:36 pm
Location: france
Contact:

Re: recording the second derivative of an action potential

Post by fabien tell »

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
fabien tell
Posts: 49
Joined: Mon Mar 25, 2013 1:36 pm
Location: france
Contact:

Re: recording the second derivative of an action potential

Post by fabien tell »

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
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: recording the second derivative of an action potential

Post by ted »

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

Code: Select all

  if (g2==nil) g2 = new Graph(0)
to

Code: Select all

//  if (g2==nil) g2 = new Graph(0)
  if (g2==nil) {
    g2 = new Graph(0)
  } else {
    g2.exec_menu("Erase")
  }
Now there is no persistence from one myrun() to the next.
(unless I mark keep lines if I want to compare outcomes) ?
"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
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))
}
fabien tell
Posts: 49
Joined: Mon Mar 25, 2013 1:36 pm
Location: france
Contact:

Re: recording the second derivative of an action potential

Post by fabien tell »

Thanks again TED,

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")

}
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
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: recording the second derivative of an action potential

Post by ted »

Good. Simpler is better. Easier to understand, easier to revise.
Post Reply