Filtering a voltage trace before fitting

[gouwens@mac.com]

I'm trying to fit a compartmental model to current-clamp data using the Multiple Run Fitter. My data are run through a low-pass analog filter before being digitized, and so I wanted to do the same thing to the voltages calculated by NEURON before they are fit by the MRF. However, I'm not exactly sure how best to implement the low-pass filter in NEURON. I see that there is a filter() function under the "Fourier" section of the Vector class, but I'm not sure how to use that. The documentation reads:

SYNTAX

Code: Select all

obj = vdest.filter(src,filter)
obj = vsrcdest.filter(filter)

DESCRIPTION
Digital filter implemented by taking the inverse fft of filter and convolving it with vec1. vec and vec1 are in the time domain and filter is in the frequency domain.

It says that the filter vector is in the frequency domain, but how does NEURON know the scaling of that vector (i.e., which values correspond to which frequencies)? Is there an easier way of doing what I'm trying to accomplish than using the Vector filter() function?

Thanks in advance.

[ted]

There are lots of ways to do this. Which is "best" depends your particular circumstances. Could you be more specific about the properties of your analog filter?

[gouwens@mac.com]

The filter that's on my rig is a 4-pole low-pass Bessel filter, with the -3dB frequency set to 5 kHz.

One idea I had was to model the actual circuitry of the filter in the Linear Circuit builder, but I thought it might be more complicated to do it that way (and so easier to do incorrectly).

The other way I had thought of was to use a Gaussian filter, because the chapter on analysis by Colquhoun and Sigworth in the "Single-Channel Recording" book mentions that a Bessel filter is well approximated by a Gaussian digital filter. When I was trying to implement that, I ran into my issues with the filter() function of the Vector class.

[ted]

I don't know of anyone who has used filter() successfully. The filter() method actually belongs to the Fourier class, the implementation of which is based on related stuff in the "Numerical Methods in Fortran | Pascal | C | C++" books. Unlike the Vector class, the Fourier class hasn't seen much use. I suspect that's because of scaling issues similar to what you encountered. A few years ago I had to write a program that used convolution and deconvolution, and most of my effort went into writing and running test code just to discover how to get consistent results when transforming data between the time and frequency domains. When that task was complete, the momentum of success and the spur of curiosity drove me to try to make sense out of filter(). But scaling was again a problem, and what I had learned from wrestling with fft and inverse fft didn't seem to help. Then I had to give up and work on a grant proposal, and I never bothered with filter() again.

You could just work entirely in the time domain, as Colquhoun and Sigworth describe in that chapter. With NEURON's Vector class you should be able to implement a hoc version of the algorithm that is the basis of the FORTRAN and Modula-2 code they present in appendix A.3.1. (well, it's A.3.1. in the 2nd edition of Single Channel Recording). This requires you to capture the time course of membrane potential to a Vector during a run, then apply the filter afterwards, e.g. with a custom run() procedure that automates the "simulate -> apply filter" sequence.

Code: Select all

proc run() {
  stdinit()
  continuerun(tstop) // time course of some variable is captured to a "raw data" Vector()
  apply_filter() // the filtered time course now exists in a "filtered data" Vector()
}

The Linear Circuit Builder approach has the advantage of producing immediate results that require no postprocessing. You'd only need a couple of op amps and a few resistors and capacitors. Circuits and parameter tables are available in many books and online sources.

[gouwens@mac.com]

Thanks for the information. I think I've got it basically working using the convlv() function, but it's a little hack-y, so I may end up going with the LinearCircuit approach anyway.

Convolving the voltage with a kernel of the form described by Colquhoun and Sigworth seems to work. There seems to be a little strangeness at the edges, though I don't need to fit at the edges. Another issue I ran into is that it looks like the MRF gets its data by using record(), which, if I understand correctly, is set in fadvance(). Since my filter is applied after continuerun() ends, I wasn't sure how to get it into the MRF without setting the data vector used by the MRF directly. It doesn't seem like an ideal way of doing it, but it apparently works. I haven't tested it very extensively yet, though.

Anyway, here's the code I've come up with (I haven't done a lot of hoc programming, so I'm sure it's not very idiomatic):

Code: Select all

objref filt_v, filt_kernel

proc run() { 
	running_ = 1
	stdinit()
	continuerun(tstop)
	apply_filter()
}

proc apply_filter() { local n, i, s, fc
  n = FitnessGenerator[0].yveclist.object(0).size()
  
  // Gaussian filter kernel from Colquhoun and Sigworth (1995), Section 2.2.1, Eq. 6 
  filt_kernel = new Vector(n)
  fc = 5000 / 1000 // cutoff frequency of 5 kHz
  for i = 0, filt_kernel.size()-1 {
    filt_kernel.x[i] = 3.011 * 5 * exp(-(5.336*fc*i*dt))
  }
  
  // Put the kernel in "wrap-around" order (see Numerical Recipes in C 13.1)
  for i = int(filt_kernel.size()/2), filt_kernel.size()-1 {
    filt_kernel.x[i] = filt_kernel.x[filt_kernel.size() - i]
  }
  
  // Normalize the kernel
  s = filt_kernel.sum()
  filt_kernel.div(s)
  
  filt_v = new Vector()
  filt_v = filt_v.convlv(FitnessGenerator[0].yveclist.object(0), filt_kernel)

  for i = 0, FitnessGenerator[0].yveclist.object(0).size()-1 {
    FitnessGenerator[0].yveclist.object(0).x[i] = filt_v.x[i]
  }
}