Getting peaks of two voltage plots

The basics of how to develop, test, and use models.
Post Reply
babaji
Posts: 3
Joined: Sat Feb 23, 2008 1:47 pm

Getting peaks of two voltage plots

Post by babaji »

Hello there, I'm a real newbie, and I was wondering if there's a way to accomplish the following thing.
I have two cells connected by a cable. I stimulate one cell and i get two plots: one is the membrane voltage (function of time) on the first cell, the second is the membrane voltage in the second cell. I get two spikes, one for each cell: the 2nd cell's spike has some delay.
I would like to plot this delay as a function of the diameter of the cable that links them. At the moment I'm simply picking the values at the peaks by hand, changing the diameter, subtracting the values, and plotting these data in Matlab.
Is there a way to do all this straight into a .hoc file? I've seen some Maximum() function, but I've jet to figure out how to use it.

Thank you very much.
B.
ted
Site Admin
Posts: 6305
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Post by ted »

First make sure you have properly conceptualized the problem. Start with the definition of
"interspike latency." If you were doing a wet lab experiment, you would discover that
measurements of peak times are harder to do, and less accurate, than measurements of
when a variable crosses a threshold.

It would be better to pick a threshold level for spike detection, say halfway between
resting potential and the peak spike amplitude, then define the interspike latency as t1 - t0
where t0 and t1 are the times when the two cells cross this level in a depolarizing
direction. With this definition, you can use the NetCon class's record() method to capture
t0 and t1. You could then automate the process by using a procedure to iterate over the
range of cable diameters, running a simulation and recording the ISI for each diameter.
This will do it:

Code: Select all

// after specification of anatomical and biophysical properties of the model
objref nc[2], nil
soma0 nc[0] = new NetCon(&v(0.5), nil)
t0 = 0
nc[0].record("t0 = t")
soma1 nc[1] = new NetCon(&v(0.5), nil)
t1 = 0
nc[1].record("t1 = t")
for ii=0,1 nc[ii].threshold = a value halfway between rest and spike peak

objref dvec, isivec
// args are min diam, diam increment, and # of runs
proc batchrun() { local ii
  dvec = new Vector($3)
  isivec = new Vector($3)
  for ii=0,$3-1 {
    cable.diam = $1 + ii*$2
    run()
    dvec.x[ii] = cable.diam
    isivec.x[ii] = t1 - t0
  }
}
Of course this assumes that changing the diameter of the cable does not significantly
affect the interval between spike detection threshold and spike peak.
babaji
Posts: 3
Joined: Sat Feb 23, 2008 1:47 pm

Post by babaji »

Ted thank you for your invaluable suggestions. My resulting plots are coherent with those obtained with the previous method. I'm even exporting some data to Matlab to get more manageable plots :)
Still, please consider the following: suppose, for instance, one of the two cells was passive, and what I wanted to grab was just the delay between the spike-peak in the active, stimulated one, and the faint (diameter-dependent) voltage peak in the passive. In that case I wouldn't be able t use the threshold method, right? I tried setting the threshold to a very low value (ie a few units bigger than the membrane resting potential) but I cannot really predict the passive cell's voltage amplitude (as a function of the diameter) and so I'm sure I'm missing some peaks. Have you got any suggestions?

Thank you again
ted
Site Admin
Posts: 6305
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Post by ted »

If you absoultely must have peak amplitudes, then record the time course of the variable
in question to a Vector, and use the Vector class's max() method to discover the value of
the maximum; record t as well, if you need to know the time at which the maximum occurs.
See the entries on the Vector class's record and max methods in the Programmer's
Reference
http://www.neuron.yale.edu/neuron/stati ... tml#Vector

You might also find this exercise from the NEURON Summer Course helpful:
6. Model Control : Simulation Families
--see http://www.neuron.yale.edu/course/handson.html
babaji
Posts: 3
Joined: Sat Feb 23, 2008 1:47 pm

Post by babaji »

Thank you once again Ted!
rosebud

Post by rosebud »

Hi Ted,
ted wrote:If you absoultely must have peak amplitudes, then record the time course of the variable
in question to a Vector, and use the Vector class's max() method to discover the value of
the maximum; record t as well, if you need to know the time at which the maximum occurs.
See the entries on the Vector class's record and max methods in the Programmer's
Reference
http://www.neuron.yale.edu/neuron/stati ... tml#Vector

You might also find this exercise from the NEURON Summer Course helpful:
6. Model Control : Simulation Families
--see http://www.neuron.yale.edu/course/handson.html
I am also trying to get the voltage values of a soma during simulation time. I've read the documentation you suggest and some other posts and ended up with the code:

objref vvec, tvec
vvec = new Vector()
tvec = new Vector()
soma vvec.record(&soma.v(.5),tvec)
run()

However when I use vvec.printf() I only get 0.
When I'm plotting the graph of the voltage axis for the soma I get a spike!
Can you tell me what is wrong?

Thanks in advance.
ted
Site Admin
Posts: 6305
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Post by ted »

Sure, and right away, if you let me look over your shoulder while you do it. Have you already
set up secure* remote desktop access, and do you have Skype working? If so, contact
me by email
ted dot carnevale at yale dot edu
to set up a time that we can work on this over the Internet.

If not, then here's how to start debugging the problem yourself.

Start with a very simple model--just one compartment, with pas or hh, and one IClamp
attached to it that delivers a brief current pulse starting at 1 ms. Make the run time just 5 ms.
Test it to make sure it works. Then set up Vector recording, and run a simulation to see if
the contents of vvec and tvec are what you'd expect from the "Voltage axis" graph of v(.5).
Tell me what you discover.

*--via SSH tunnel or a VPN
rosebud

Post by rosebud »

Hi Ted,

the bug was trivial (and rather silly!!) afterall.
I didn't record time...
I added:

tvec.record(&t)

and everything worked out just fine!!

Thanks for your interest and for being so devoted :)
ted
Site Admin
Posts: 6305
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Post by ted »

rosebud wrote:the bug was trivial (and rather silly!!)
Many are, and everybody makes them. An interesting note: it often happens that a trivial
bug in one's own code is very hard to find, but the same bug in somebody else's code
seems to jump right out of the page.
Post Reply