calculating the "area under the curve"

The basics of how to develop, test, and use models.
Post Reply
udi
Posts: 13
Joined: Wed Aug 10, 2005 3:25 am
Location: Hebrew University, School of Medicine

calculating the "area under the curve"

Post by udi »

Hi,
I want to compare the size of different ADPs. How do I calculate the area under the curve (e.g - start point is the end of the fAHP / repolarization, end point is RP)?

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

integrating simulation results

Post by ted »

This is duck soup for the Vector class, whose methods include
many powerful tools for postprocessing in the time and frequency
domains.

Here are the basic steps to accomplish your task. You should read
the documentation of these various methods in the Programmers'
Reference (AKA "Documentation" in the NEURON program group,
if you're using MSWin).

1. Record the waveform to a Vector. Better, record both the y and t
data to a pair of vectors.
objref vvec, tvec
vvec = new Vector()
tvec = new Vector()
// replace secname and x with the name of the section and the "range"
// at which you want to record v
secname vvec.record(&v(x), tvec)

2. Isolate the portion of the y data that contains the transient you want
to integrate. For the sake of this example, suppose you are interested
in the values of v for t in the interval t1...t2
i1 = tvec.indvwhere(">=", t1)
i2 = tvec.indvwhere(">=", t2)
objref vx // "extract" of vvec
vx = vvec.c(i1+1, i2)
A question for the reader: why start at i1+1?

3. Subtract the baseline. Let's say you have decided that the baseline
is the average v over the time interval tb1...tb2. Use indvwhere() to
find the indices ib1 and ib2 that correspond to tb1 and tb2 ("this is left
as an exercise for the reader"). Then
vbaseline = vvec.mean(ib1, ib2)
vx.sub(vbaseline)
Now vx contains the v data minus the baseline value.

4. Integrate the baseline-subtracted v data. You just want the total area,
so
varea = dt*vx.sum() // assuming fixed dt simulation

Two more questions for the reader:
1. How would you change these steps to implement Simpson's rule
integration?
2. How would you compute the integral if the simulation used
variable dt?
eacheon
Posts: 97
Joined: Wed Jan 18, 2006 2:20 pm

Re: integrating simulation results

Post by eacheon »

ted wrote: ...
1. Record the waveform to a Vector. Better, record both the y and t
data to a pair of vectors.
objref vvec, tvec
vvec = new Vector()
tvec = new Vector()
// replace secname and x with the name of the section and the "range"
// at which you want to record v
secname vvec.record(&v(x), tvec)

2. Isolate the portion of the y data that contains the transient you want
to integrate. For the sake of this example, suppose you are interested
in the values of v for t in the interval t1...t2
i1 = tvec.indvwhere(">=", t1)
i2 = tvec.indvwhere(">=", t2)
objref vx // "extract" of vvec
vx = vvec.c(i1+1, i2)
A question for the reader: why start at i1+1?

3. Subtract the baseline. Let's say you have decided that the baseline
is the average v over the time interval tb1...tb2. Use indvwhere() to
find the indices ib1 and ib2 that correspond to tb1 and tb2 ("this is left
as an exercise for the reader"). Then
vbaseline = vvec.mean(ib1, ib2)
vx.sub(vbaseline)
Now vx contains the v data minus the baseline value.

4. Integrate the baseline-subtracted v data. You just want the total area,
so
varea = dt*vx.sum() // assuming fixed dt simulation

Two more questions for the reader:
1. How would you change these steps to implement Simpson's rule
integration?
2. How would you compute the integral if the simulation used
variable dt?
Ted. How can I compute the integral if the simulation is using variable dt? I cannot figure out a good answer for myself. Can you detail on the answer a little bit?

Would summing a series of integral by Simpson's rule will get me 2 times the answer I wanted, if I divide the recording into (t1,t2,t3), (t2,t3,t4), (t3,t4,t5)... and integrate them individually? or simply use Vector.dot to do a dot product and stop worrying about the accuracy?

How should I process the first and the last interval then? Interval(t1,t2) and (t(n-1), t(n)) are only counted once, so the sum is not exactly twice of the "accurate integral", but some how a little less than that. How can I deal with this?
ted
Site Admin
Posts: 5892
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Post by ted »

How can I compute the integral if the simulation is using variable dt?
EITHER
use a method designed to handle a continuous variable sampled at irregular intervals (there
must be at least one such method in the vast literature on numerical computation)
OR
sample the data stream at regular intervals, then do fixed interval integration.
Sampling at regular intervals can be done in the course of simulation by using the Vector
class's record() method with the vdest.record(&var, Dt) syntax. Or you could apply the
Vector class's interpolate() method to a data stream sampled at irregular intervals.
eacheon
Posts: 97
Joined: Wed Jan 18, 2006 2:20 pm

Post by eacheon »

ted wrote:
How can I compute the integral if the simulation is using variable dt?
EITHER
use a method designed to handle a continuous variable sampled at irregular intervals (there
must be at least one such method in the vast literature on numerical computation)
OR
sample the data stream at regular intervals, then do fixed interval integration.
Sampling at regular intervals can be done in the course of simulation by using the Vector
class's record() method with the vdest.record(&var, Dt) syntax. Or you could apply the
Vector class's interpolate() method to a data stream sampled at irregular intervals.
Thanks, I am doing the latter. But kind of interested in the first choice. not well educated in numerical computation, I only learned from your previous post that Simpson's rule is one way that does irregular intervals integration. Am I on the right track to the solution?
ted
Site Admin
Posts: 5892
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Post by ted »

eacheon wrote:I only learned from your previous post that Simpson's rule is one way that does irregular intervals integration.
I didn't say that, as you will discover if you re-read my prior messages. There is
something called "adaptive Simpson's rule" which can be used for irregular intervals,
but what I had in mind was the plain old composite Simpson's rule, which requires fixed
intervals. However, that has certain prerequisites (requires an even number of intervals,
and the first 4 derivatives of the function to be integrated must exist and be continuous)
that may limit its general application, and its 4th order accuracy is probably higher than is
justified by the data that are to be integrated.

So it's probably best to stick with the composite trapezoidal rule for fixed intervals, which
has second order accuracy. This function should do the job:

Code: Select all

// does trapezoidal integration of data sampled at fixed intervals
// expects 2 args: vector to be integrated, and data sampling interval

func trapint() {
  return $2*($o1.sum() - 0.5 * ( $o1.x[0] + $o1.x[$o1.size()-1] ))
}
Post Reply