calculating the "area under the curve"

[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]

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]

...
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]

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]

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]

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] ))
}