Neural Tuning Curve

Computational Neuroscience
Chris Harris
31 May 2015

Problem Statement

Question 11

In the following three questions, we will explore Poisson neuron models and population coding.

This exercise is based on a set of artificial “experiments” that we’ve run on four simulated neurons that emulate the behavior found in the cercal organs of a cricket. Please note that all the supplied data is synthetic. Any resemblance to a real cricket is purely coincidental.

In the first set of experiments, we probed each neuron with a range of air velocity stimuli of uniform intensity and differing direction. We recorded the firing rate of each of the neurons in response to each of the stimulus values. Each of these recordings lasted 10 seconds and we repeated this process 100 times for each neuron-stimulus combination.

We’ve supplied you with a file containing data for each of the neurons that contains the recorded firing rates (in Hz). These are named neuron1, neuron2, neuron3, and neuron4. The stimulus, that is, the direction of the air velocity, is in the vector named stim.

The matrices contain the results of running a set of experiments in which we probed the synthetic neuron with the stimuli in stim. Each column of a neuron matrix contains the firing rate of that neuron (in Hz) in response to the corresponding stimulus value in stim. That is, the n^th column of neuron1 contains the 100 trials in which we applied the stimulus of value stim(n) to neuron1.

Plot the tuning curve– the mean firing rate of the neuron as a function of the stimulus– for each of the neurons.

Question 12

We have reason to suspect that one of the neurons is not like the others. Three of the neurons are Poisson neurons (they are accurately modeling using a Poisson process), but we believe that the remaining one might not be. Which of the neurons (if any) is NOT Poisson?

Hint: Think carefully about what it means for a neuron to be Poisson. You may find it useful to review the last lecture of week 2. Note that we give you the firing rate of each of the neurons, not the spike count. You may find it useful to convert the firing rates to spike counts in order to test for “Poisson-ness”, however this is not necessary.

In order to realize why this might be helpful, consider the fact that, for a constant a and a random variable X:

𝔼[aX] = a𝔼[X], however; Var(aX) = a²Var(X).

What might this imply about the Poisson statistics (like the Fano factor) when we convert the spike counts (the raw output of the Poisson spike generator) into a firing rate (what we gave you)?

Question 13

Finally, we ran an additional set of experiments in which we exposed each of the neurons to a single stimulus of unknown direction for 10 trials of 10 seconds each.

pop_coding contains four vectors named r1, r2, r3, and r4 that contain the responses (firing rate in Hz) of the four neurons to this mystery stimulus. It also contains four vectors named c1, c2, c3, and c4. These are the basis vectors corresponding to neuron 1, neuron 2, neuron 3, and neuron 4.

Decode the neural responses and recover the mystery stimulus vector by computing the population vector for these neurons. You should use the maximum average firing rate for a neuron as the value of r_max for that neuron. That is, r_max should be the maximum value in the tuning curve for that neuron.

What is the direction, in degrees, of the population vector? You should round your answer to the nearest degree. Your answer should contain the value only (no units!) and should be between 0^∘ and 360^∘. If your calculations give a negative number or a number greater than or equal to 360, convert it to a number in the proper range (you may use the mod function to do this).

You may need to convert your resulting vector from Cartesian coordinates to polar coordinates to find the angle. You may use the atan() function in MATLAB to do this. Note that the the convention we’re using defines 0^∘ to point in the direction of the positive y-axis and 90^∘ to point in the direction of the positive x-axis (i.e., 0^∘ = north, 90^∘ = east).

Python Code

tune.py:

import math
import pickle
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

class stat:

    def max(self,X):

        m = X[0]
        for e in X:
            if e > m:
                m = e
        return m

    def mean(self,X):

        s = 0.
        for e in X:
            s = s + e
        return s/len(X)

    def var(self,X,mx):

        s = 0.
        for e in X:
            s = s + (e-mx)**2
        return s/(len(X)-1)

    def cor(self,X,Y,mx,my,sx,sy):

        s = 0.
        for i,e in enumerate(X):
            s = s + (X[i]-mx)*(Y[i]-my)
        return s/sx/sy/(len(X)-1)

    def regression(self,X,Y):

        # mean
        mx = self.mean(X)
        my = self.mean(Y)

        # standard deviation
        sx = math.sqrt(self.var(X,mx))
        sy = math.sqrt(self.var(Y,my))

        # correlation
        rxy = self.cor(X,Y,mx,my,sx,sy)

        # slope, intercept, coefficient
        p    = [[] for n in range(3)]
        p[0] = rxy*sy/sx
        p[1] = my - p[0]*mx
        p[2] = rxy
        return p

st = stat()

class tune:

    def __new__(self):

        key    = [['stim','neuron1','neuron2','neuron3','neuron4'],['r1','r2','r3','r4','c1','c2','c3','c4']]
        idim   = len(key[0])
        tarray = self.calc(key[0],idim)
        parray = self.pop(key[1])
        rmax   = self.peak(tarray,idim)
        self.angle(parray,rmax)
        self.plot(tarray,idim)

    def calc(key,idim):

        data   = open("../data/tuning.pickle",'rb')
        stream = pickle.load(data)
        data.close()

        array    = [[] for i in range(idim+4)]
        array[0] = stream[key[0]]
        jdim     = len(array[0])
        set      = []

        for i in range(1,idim):
            set = stream[key[i]]
            mx = []
            vx = []
            for j in range(jdim):
                row = []
                for k in range(jdim):
                    row.append(set.item((k,j)))
                mx.append(st.mean(row))
                vx.append(st.var(row,mx[-1]))
            array[i]   = mx
            array[i+4] = vx
        return array

    def pop(key):

        data   = open("../data/pop_coding.pickle",'rb')
        stream = pickle.load(data)
        data.close()

        dim   = len(key)
        array = [[] for i in range(dim)]
        for i in range(dim):
            array[i] = stream[key[i]]
        return array

    def peak(array,idim):

        rmax = []
        for i in range(1,idim):
            rmax.append(st.max(array[i]))
        print(f"\n\tRmax  = [{rmax[-4]:7.3f},{rmax[-3]:7.3f},{rmax[-2]:7.3f},{rmax[-1]:7.3f} ]")
        return rmax

    def angle(array,rmax):

        s = [0. for i in range(2)]
        for i in range(2):
            for e in array[i]:
                s[0] = s[0] + e/rmax[i]*array[i+4][0]
                s[1] = s[1] + e/rmax[i]*array[i+4][1]
        angle = math.degrees(math.atan(s[1]/s[0]))
        print(f"\tVpop  = [{s[0]:7.3f},{s[1]:7.3f} ]\n\tangle = {angle:10.1f} deg\n")

    def plot(tarray,idim):

        mpl.rc('text', color='#C8A078')
        mpl.rc('figure', facecolor='black', edgecolor='black')
        mpl.rc('axes', edgecolor='#C6BDBA', labelcolor='#C8A078', facecolor='black', linewidth=2)
        mpl.rc('xtick', color='#C8A078')
        mpl.rc('ytick', color='#C8A078')
        plt.tick_params(bottom=False, top=False, left=False, right=False)

        plt.figure(1)
        col = ['#640000','#004000','#004040','#640064']
        for i in range(1,idim):
            plt.plot(tarray[0], tarray[i], color=col[i-1], label="neuron %i" %i, linewidth=3)
        plt.legend(loc=(0.85,0.8), frameon=False, fontsize=20)
        plt.title('Tuning Curve', fontsize=25)
        plt.xlabel('Angle [deg]', fontsize=25)
        plt.ylabel('Neuron firing rate [Hz]', fontsize=25)
        plt.xticks(fontsize=20)
        plt.yticks(fontsize=20)

        plt.figure(2)
        plt.xlim(0,35)
        plt.ylim(0,5)
        mark = ['o','s','^','d']
        for i in range(1,idim):
            plt.plot(tarray[i], tarray[i+4], color=col[i-1], marker=mark[i-1], markersize=20,
                     markeredgewidth=3, fillstyle='none', label='neuron %i' %i, linewidth=0)
        for i in [1,3]:
            p = st.regression(tarray[i],tarray[i+4])
            Y = []
            for e in tarray[i]:
                Y.append(p[0]*e + p[1])
            plt.plot(tarray[i], Y, color='#F5A078', linewidth=3)
        plt.legend(loc=(0.85,0.05), frameon=False, numpoints=1, fontsize=20)
        plt.title('Poisson Check', fontsize=25)
        plt.xlabel('Neuron firing rate [Hz]', fontsize=25)
        plt.ylabel('Variance [Hz^2]', fontsize=25)
        plt.xticks(fontsize=20)
        plt.yticks(fontsize=20)
        plt.show()

def main():

    t = tune()

main()

Results

$ python tune.py

    Rmax  = [ 32.092, 21.625, 14.617, 26.700 ]
    Vpop  = [  9.095, -3.657 ]
    angle =      -21.9 deg

Figure 1: Tuning Curve

Figure 2: Poisson Check