Problems Using FFT to Compute Impedance in a Model Neuron

Question

I'm a neuroscientist currently investigating the resonance properties of a single neuron model that a colleague and I have constructed. The language we code in is Julia, which I hope is similar enough to Python for most people to understand. I've run in to a very particular problem, so bear with me please.

In order to investigate the resonance properties of the model neuron, I am applying an exponentially swept sinusoid current to it and analyzing its voltage response. In my field, such a current is called a ZAP current, though a more familiar term is probably Chirp signal. Just so there is no confusion, the current is generated from the following function \begin{gather*} ZAP(t) = I_{app}sin(2\pi f(t)t)\\f(t) = (0.001)\frac{f_{min}}{Lt}(e^{Lt}-1)\\ L = \frac{\ln{\frac{f_{max}}{f_{min}}}}{t_{max}} \end{gather*} where $t_{max}$ is the time duration of the simulation, $f_{min}$ and $f_{max}$ are the starting and ending frequencies of the sweep, and $I_{app}$ is the magnitude of the sinusoid. The voltage response of the model neuron to this input looks like this:

Clearly, there is a resonance in the model at whatever the input frequency is at around 50,000 ms. This means there should be a peak in the impedance of the model at this frequency. The impedance is calculated as the ratio of the discrete Fourier transform of the output voltage to the input current $$ Z = \frac{FFT(V(t))}{FFT(ZAP(t))} $$ I then plot the magnitude of the impedance against the natural logarithm of the frequency. I expect to get something like the one on the left, but instead I obtain the one on the right:

As I vary the parameter $E_K$, the resonance frequency changes, and I expect to see this reflected in the impedance-frequency plot, but every one of the impedance-frequency plots looks almost identical to the one shown above despite the voltage responses looking quite different.

I have spent hours going through my code line by line and experimenting with different modifications to figure out where I'm going wrong, but have not been able to produce the desired output. Here is the stripped down code:

## Time Parameters
dt      = 0.01         # ms
tmax    = 100000        # ms
tspan   = convert(Int64,tmax/dt)
t       = linspace(dt,tmax,tspan)

## Simulation and Analysis
# magnitude of sinusoid
Iapp = 1.0
# minimum frequency
f0 = 0.1
# maximum frequency
ff = 15.0
# create an array of the frequency sweep values for plotting
Zfreq = ZAPfreq(f0, ff)
# optimize Fast Fourier Transform for an array of size and type ZAP()
P = plan_fft(ZAP(Iapp, f0, ff))
# calculate FFT for input ZAP()
inFFT = P*ZAP(Iapp, f0, ff)
output    = simulate(ZAP(Iapp, f0, ff))
# calculate the FFT of the membrane voltage response
outFFT = P*output
# calculate the impedance of the cell at each frequency ZAPfreq[i] from Z = V/I
Z = outFFT ./ inFFT
# calculate the magnitude of the complex impedance
MagZ = abs(Z)
# center the DC response (f = 0 at center of frequency axis)
shiftZ = fftshift(MagZ)
# plot the impedance at each frequency
plot(Zfreq, shiftZ, grid = false,
  legend = false,
  linecolor= :black,
  xscale= :ln,
  xlabel = "log(\$f\$)",
  ylabel = "Impedance (Ω)",
  title = "Kispersky-Caplan Impedance Response with \$E_K\$ = $EK",
  top_margin=5*mm,
  bottom_margin=5*mm,
  left_margin=5*mm,
  right_margin=10*mm,)   

Here, the $\texttt{simulate()}$ function is the one that actually carries out the simulation of the model. Any help in determining the problem is greatly appreciated.

Answer 1

The impedance diagram to the left looks like a Bode diagram. This means the ordinate should not be a linear scale, but rather a dB scale. Hence, you probably should be plotting 20log10(abs(Z)) and not abs(Z). Give this a try, see what happens.

P.S. For a proper Bode diagram, you should be plotting 20log10(abs(Z)) vs. log10(f), not ln(f).

Answer 2

I don't know how much background you have with FFT, so I apologize if I list a few things you already know.

In general, I think people have a tendency to throw an FFT at a problem and hope it works, but you may need to consider what the transform is really all about. First, when you take an FFT, it's implicit in the calculation that your signal is periodic. So, even though you're feeding it a finite number of samples with frequency content only between f0 and ff, the calculation makes the assumption that before time zero and after your samples, the waveform repeats forever. If your waveform's beginning and end don't match perfectly, there's going to be a lot of frequencies you wouldn't expect in your FFT. Basically, it's like there's a discontinuity. So, often, you'll window the signal to force it to zero at both ends so each periodic clone of your waveform meets the next smoothly. This obviously also has consequences for the frequency content of the FFT, so you need to do this carefully also.

Second, if your waveforms have any DC content (a constant offset, frequency = 0 Hz), it usually makes sense to remove it before your FFT unless DC is of interest to you. Often, low frequencies are huge and swamp the stuff you care about. Doing some high-pass filtering with a very low cutoff will help you see the interesting content of your plot more easily.

Third, when you view your plot on a linear frequency axis, you'll see that there are Fourier coefficients for both positive and negative frequencies, and the plot is symmetrical about 0 Hz. For qualitative assessments, you probably only need to look at the positive half.

Fourth, consider that you're taking the ratio of output over input, and anywhere your input FFT is near zero, that's going to blow up, and may make it hard to see what's going on.

So, I'd say, do some preprocessing of your signal to remove things you don't care about, and to prevent edge effects. Then post another plot and let's see how it looks. Since you're simulating the input and output here, there may not be a need, but in practice, I'll often use Welch's method to compute power spectral densities, and then compute an H1 or H2 estimation of the frequency response. Here, we can hopefully get a result with a simple analysis.

Answer 3

I have no idea about Neuron modelling so pardon me if my assumption is wrong but I can see that there appears to be a problem in in your electrical model.

Impedance is : $ I(\Omega) = V(t)/I(t).$ But ratio of frequency volatge (V(t)) and current I(t) is not impedance. This quantity (ratio of frequencies), might have some significance in your model but that is not impedance. I would suggest trying FFT of the impedance function instead, which is : $ Z = FFT(\ V(t)/ZAP(t)\ ) $

Also you seem to interested in input frequency at a particular time (where there is resonance), you can so this by STFT (Short Time Fourier Transform) or Wavelet transform.