f1/f0 demo

One can imagine several ways that a neuron might respond to a sinusoidal input, which is our normal visual input.

First, let's express the sinusoidal input signal S as a function of time. We will write out 20 seconds of time in 1ms steps. Before time 0, the signal will be 0, and it will start being sinusoidal at time 0.

start_time = -10;

stop_time = 20;

t = start_time:0.001:stop_time;

f = 4; % 4 Hz frequency, or 4 cycles per second

S = (t>=0).*sin(2*pi*f*t);

plot(t,S)

box off; % turn off the box around the plot

% zoom in a bit

axis([-1 1 -1.5 1.5])

ylabel('Signal')

xlabel('Time(s)')

There are many ways that our neuron could respond to this signal. Let's first think about its voltage response. The voltage could simply fluctuate sinusoidally as a function of time.

delay = 0.1; % the response is a bit delayed in a real system like a neuron

gain = 0.5; % 0.5 response units / input signal unit

Y_1 = (t>=0+delay).*gain.*sin(2*pi*f*(t-delay)); % .* means element-by-element multiplication

plot(t,S)

hold on;

plot(t,Y_1,'g');

box off; % turn off the box around the plot

% zoom in a bit

axis([-1 1 -1.5 1.5])

ylabel('Voltage (V)')

xlabel('Time(s)')


If we had the output signal Y_1, how can we measure the magnitude of the response? We'd ideally like to measure the amplitude of the sinusoidal response. The technique for doing this is called Fourier analysis (Lab 3.2). We have a Matlab function that performs Fourier analysis: **fouriercoeffs_tf2**. fouriercoeffs_tf2 takes the signal, the frequency to be examined, and the sample rate.

sample_rate = 1/0.001; % we had 1 ms steps, so the sample rate is 1000 Hz, or samples per second

A = fouriercoeffs_tf2( (Y_1((numel(Y_1)-1)/2:numel(Y_1)))', f, sample_rate)

A =

-0.2732 + 0.4187i

The fourier coefficient at 4 Hz is a complex number that has a magnitude and phase. Let's check out the magnitude by using the **abs** (absolute value) command:

abs(A)

ans =

0.5000

We got a value of 0.4940, which is very close to the magnitude of the acutal response (0.5). If there were no delay, then this would be exactly right.

Another cell might give a response that looks like this:

Y_0 = (t>=0+delay).*gain; % .* means element-by-element multiplication

plot(t,S)

hold on;

plot(t,Y_0,'g');

box off; % turn off the box around the plot

% zoom in a bit

axis([-1 1 -1.5 1.5])

ylabel('Voltage (V)')

xlabel('Time(s)')

Here, the voltage response to the input with a fixed step. Taking the mean of the response does a good job of capturing the response.

B = mean( (Y_0((numel(Y_0)-1)/2:numel(Y_0)))')

B =

0.5000

So, here again we have a value very close to the true function that was used to generate, with the only small difference arising from the delay in the response.

We can also imagine a mixed response, like the following:

Y = Y_0 + Y_1;

plot(t,S)

hold on;

plot(t,Y,'g');

box off; % turn off the box around the plot

% zoom in a bit

axis([-1 1 -1.5 1.5])

ylabel('Voltage (V)')

xlabel('Time(s)')

We can calculate the mean of this response, and the F1 response (the amplitude of the response at the fundamental stimulus):

B_ = mean( (Y((numel(Y)-1)/2:numel(Y)))')

B_ =

0.5000

Even though this captures the mean response very well, it does not describe the full response. It does not describe the sinusoidal modulation of the response. To capture this, we must also examine the F1 of the response:

A_ = fouriercoeffs_tf2( (Y_1((numel(Y_1)-1)/2:numel(Y_1)))', f, sample_rate)

abs(A_)

A_ =

-0.2732 + 0.4187i

ans =

0.5000

Now, together, these two numbers (the mean, also called F0, and F1, the amplitude of the signal at the stimulus frequency) capture the response of the cell pretty well.

We could also imagine that we do not have a voltage response, but, instead, we have a firing rate that is based on this voltage response. For the moment, we will let the **firing rate** of a neuron be equal to 10 spikes/(sec/volt) * Y.

FR = 10 * Y;

plot(t,FR);

box off; % turn off the box around the plot

% zoom in a bit

axis([-1 1 0 10])

ylabel('Firing rate (spikes/sec)')

xlabel('Time(s)')

We can calculate the probability that the cell fires in any given millisecond by calculating:

P = FR * 0.001;

plot(t,P);

box off; % turn off the box around the plot

% zoom in a bit

axis([-1 1 0 0.02])

ylabel('Probability of firing')

xlabel('Time(s)')

Now we can generate some spikes by drawing random numbers:

F = rand(size(P)) < P;

plot(t,S)

hold on;

plot(t,F,'r')

box off; % turn off the box around the plot

% zoom in a bit

axis([-1 3 0 1])

ylabel('Signal (blue) and spikes (red)')

xlabel('Time(s)')

We can calculate the mean response:

F_mean = sum(F)/(stop_time) % count the spikes, divide by the time

F_mean =

5.1500

The mean response has units of spikes/sec. But knowing the mean response doesn't fully describe the whole response. We need to consider the rhythmicity of the signal as well by calculating the F1. For calculating this for spikes, we will find the spike times and compute the fourier transform:

spiketimes = t(find(F));

F1_resp_manual = (2/stop_time)*sum(exp(-sqrt(-1)*2*pi*f*(spiketimes)))

abs(F1_resp_manual)

% or

F1_resp = fouriercoeffs_tf_spikes(spiketimes, f, stop_time)

abs(F1_resp)

F1_resp_manual =

-3.2691 + 3.7591i

ans =

4.9817

F1_resp =

-3.2691 + 3.7591i

ans =

4.9817

This indicates that the response has a high degree of rhythmicity; the magnitude is the same size as the mean response. It has the units of spikes/sec.