A    Introduction

The linear-nonlinear model of a neuron is a way to predict the neural response to a stimulus in terms of an estimated firing rate for every time point in the stimulus sequence. Not every neuron behaves according to this simple model, but it can be a starting point for a more complete characterization.

The model makes the assumption that the neuron scans the input stimulus for the presence of a specific feature.  The more closely the stimulus resembles this trigger feature, the greater the probability that the neuron will fire an action potential. What this feature is can be extracted by computing the spike-triggered average (see below).

The simple assumption that twice as much trigger feature will cause twice as much response is often too simple, as neurons exhibit nonlinear behavior, e.g. saturation at high stimulus intensities. The same is true for low trigger feature intensities as these are usually below the threshold for firing. A good description is a sigmoidal curve connecting trigger feature intensity with spiking rate.

So what is the linear part and what is the nonlinearity? The linear part of the model is the assumption that the neuron runs a linear filter (the spike-triggered average) over the data to find the intensity of the trigger feature.  The filter output at each point in time is then transformed into a spiking rate through a nonlinear, often sigmoidal, function. In addition, we assume that the nonlinearity is static: it will not change over time. One reason that this assumption can be wrong is that the neuron may adapt to stimulus contrast. This happens, for example, in the retina.

B    The Spike-Triggered Average

B.1    Introduction

For certain stimuli, an unbiased estimate of the linear filter is the spike-triggered average.  In this method, the neuron is probed with a white noise stimulus sequence. From time to time, the neuron will fire an action potential because, by chance, the noise resembled part of the “trigger feature” for the cell. Because a long sequence of white noise probes all aspects of a possible trigger feature with equal probability, we will be able to reconstruct this optimal stimulus piece by piece. When we build the average of all white noise sequences that precede a spike, unimportant features in the noise average out and important features add up, giving in theory the stimulus to which the cell’s response is tuned.

B.2    Adaptation of Retinal Ganglion Cells to Different Contrasts Levels

Read pages 19 - 23 in Dayan and Abbott to understand the details of the spike-triggered average analysis using white noise stimuli for neurons of the early visual system. Download the file contrast response.mat from Canvas containing the stimulus and the response of a retinal ganglion cell to different contrasts of a random full-field flicker stimulus, i.e.  randomly changing the brightness of the entire field of vision.  (The data is from experiments by Bart Borghuis and Kristy Simmons, Balasubramanian lab.) The data contains nine trials each 100 sec long at 2 kHz sampling rate. Each trial has a different contrast (standard deviation divided by mean) of random full-field flicker. The contrast steps are given in the variable called contrasts. The variable stimulus contains the random stimulus sequence for  one  trial.  The stimulus at each time is given as the deviation of the monitor intensity from its mean value.  All nine trials use this identical stimulus sequence, scaled by a factor contrasts(i).  The variable spikeTimes{i} contains all spike times (in samples) for the ith  trial relative to the beginning of that trial.

Problem 1: (a) Display a plot of the spikes in the first 20 seconds of the trial with the highest contrast.

(b) For the trial with the highest contrast, plot the time-binned firing rate for the first 20 seconds.  To do that you can use the output of MATLAB’s histogram function. Use 100ms time bins. Be sure to scale the output so that the firing rate is measured in spikes/s.

(c) The firing rate of a neuron is often estimated by assuming it changes smoothly across time. Accordingly, a common step before further analyses are done is to smooth the spike vector (i.e.  a binary vector with 0 indicating there is no spike recorded during a particular sample, and 1 indicating there was a spike recorded at said sample) rather than using set bins as we did in problem problem (1.b). This smoothing is done by convolving the spike vector with a Gaussian with a set standard deviation; for our purposes 35 ms will work well for the latter. Also note, the normalization is important here: ensure that the values in each Gaussian sum up to the sampling rate. MATLAB and Python both have their own built-in convolution functions that you can use.  Pay attention to the input arguments for these functions, in particular, their order; you will need them again for later questions.

(d) Plot the overall average firing rate as a function of contrast.

(e) Generate and plot the spike-triggered averages for all nine different contrast levels.  Normalize each of them so that the sum of the squared values is one. You will have to find out how far back in the past from the spikes you have to go to catch the complete spike-triggered average.

(f) Describe how the spike-triggered average changes with contrast. What does that tell us about the behavior of this retinal ganglion cell?

C    The Static Nonlinearity

Before you start this question, you should make sure to complete the previous one, and pay attention to parts (c) and (e) in which the normalization of the spike triggered average is important.  Read pages 48 - 51 of the textbook to understand the concept of static nonlinearities.

Problem 2: (a) For every trial, use the corresponding STA from problem 1 as a linear filter (as described in the lecture) to extract the trigger feature intensity as a function of time. Be sure to scale the stimulus by the appropriate contrast value, in order to match the actual stimulus that was displayed during the experiment.

Hint:  If using the conv or convn functions in MATLAB, be aware that the kernel you are convolving with is flip from left to right before the convolution takes place. You can get around this by using the fliplr function on the kernel before convolving.

(b) Compare your smooth spike rate curves from problem 1c with the trigger feature intensities from problem 2a.  Tune the width of your Gaussian kernel until you get the best agreement between the two curves (in terms of how smeared out the bumps in both graphs look). A good way to do that is to temporally set all values of the trigger feature curve that are smaller than zero equal to zero so you get only positive bumps as in the firing rate plot. Also have a look at figure 2.26. This step is essential to get good results! It should be helpful to start with the highest intensity trial. Plot your result showing both the trigger feature intensity graph and the smoothed spiking rates for some of the contrasts. Be sure not to introduce an artificial offset between the two curves.

Note:   The tip to set values of the trigger feature curve that are less than zero equal to zero is  only for help comparing the graphs visually.  Do not have these values set to zero for the later questions.  This introduces a rectifying non-linearity before you calculate the non-linearity and will negatively impact your answer.

Now we want to generate the static nonlinearities for all nine trials.

(c) Plot a scatter plot of the filter outputs (x-values) against the firing rates (y-values) by using the two curves you just constructed. Using every 10th  time point in both plots should give you enough data points. This gives you something like figure 2.2 in Dayan and Abbott (for now without the fitted curve).  Do this for all nine trials.

(d) Average the points in your scatter plot over the x-axis until it is possible to draw a nice graph of the nonlinearity for all nine trials. You want each value of x in your graph to correspond to a single value of y. (You should aim for 20-40 points along the x-axis) What do you observe when you compare the nonlinearities for the different intensities?

(e) Apply your estimated nonlinearity to the filter output to generate the model’s predicted firing rate. This can be done with the interp1 command.  interp1(x,y,x0) interpolates between a list of (x, y) points to give the function values at the points x0.  (x,y) will be the points on the nonlinearity that you found in problem 2d, and x0 will be the full filter output. Plot your result together with the actual firing rate. Is the linear-nonlinear model a fair description of the response of the neuron? Check different contrast intensities.

D    Simulating V1 Simple Cells

The primary visual cortex, also called V1, is an early place (and functional unit) in the brain where the visual image coming from the retina is analyzed.  One purpose of V1 is to build a ”sparse” representation of the visual scene in terms of Gabor functions. Gabor functions will represent the image as a collection of moving edges which then are further analyzed in associated visual cortex areas (edge detection).  You can find a typical sparse representation in terms of Gabor functions in figure 10.7 of our book. Building a sparse (which means efficient and compressed) representation of an image using Gabor functions is also used in image compression.  The so-called simple cells in V1 are the key elements in building that representation. There are also cells in V1 doing ”higher order” processing called complex cells.  They are the first step in getting a translationally invariant representation of objects. You want to identify an object independently of its exact position in space! It is strongly recommended to read the introduction to V1 cells in Dayan and Abbott pages 58 - 76 for a better understanding of this homework set.

The simple cells in V1 build up their receptive field from a clever combination of inputs from retinal ganglion cells.  How this is done is still not understood completely.  Instead, we can measure the receptive field of simple cells directly via STA analysis without caring about how the cell is actually getting its receptive field. A good approximation of a receptive field of a V1 simple cell is the so-called Gabor function. A Gabor function is a cosine wave with a Gaussian envelope: a ”wave packet.” The functional description of a Gabor function is:

Ds(x,y) =  exp  −  −  cos(kx − ϕ)                                  (2.27)

where σx  and σy  determine the extent of the receptive field in the x and y directions. ϕ is the preferred spatial phase and k is the preferred spatial frequency.

At the same time the cell responds to temporal light variation. This can be described with the following functional form:

Dt(τ) = αexp( −ατ)  −                                              (2.29)

where α sets the timescale. At time τ = 0 the neuron is most likely to fire an action potential. Bigger τ means going backwards in time. It is important to note that the function Dt(τ) is only different from zero for τ > 0.  This will become important in part 2.  A full description of the spatiotemporal receptive field (which we could get with an STA analysis) can be given as a product of the time-dependent part Dt(τ) and the space-dependent part Ds(x,y):

D(x,y,τ) = Ds(x,y)Dt(τ)                                                        (2.a)

Because V1 simple cells are a kind of spatial frequency detector (because the Gabor function looks like a localized cosine wave), the cell can be probed with a time varying counterphase sinusoidal grating.  The functional form of that grating is:

s(x,y,t) = Acos(KxcosΘ + Ky sinΘ − Φ)cos(ωt)                                   (2.18)

K is the spatial frequency, ω describes the temporal frequency, Θ is the orientation of the grating, Φ is the spatial phase and A is the overall amplitude.  We will set A = 1.  The parameters of the counterphase sinusoidal grating are visualized in figure 2.8.  The counterphase sinusoidal grating is not moving towards a direction.  Instead it is oscillating with frequency ω/2π.  By changing the parameters of the grating we can understand what angles or spatial frequencies a specific V1 cell is tuned to.  By doing that, one finds that there are all kinds of simple cells in V1, tuned to all possible angles, spatial frequencies, and temporal frequencies.

Because of the special structure 2.a of the receptive field, the spatial and the temporal parts of the receptive field response separate:

L(t) = LsLt(t)                                                                 (2.b)

The spatial linear response Ls  is just the grating filtered with the spatial part of the receptive field Ds(x,y):

Ls = Zdxdy Ds(x,y)Acos(Kxcos(Θ) + Ky sin(Θ) − Φ)                               (2.31)

The temporal linear response function Lt(t) is the temporal oscillation filtered with the temporal part of the receptive field Dt(τ):

Lt(t) = Z0 ∞ dτ Dt(τ)cos(ω(t − τ))                                                (2.32)

Problem 3:

First we want to construct the receptive field of a V1 cell according to equations 2.27, 2.29 and 2.a.

(a) Construct the spatial part of a vertically oriented V1 simple cell using function 2.27. Use the values from figure 2.15, as we will need them later. In particular, pick the special case k = 1 and σx = σy  = σ = 2. Take the phase ϕ = 0 as given in figure 2.15. Make plots of your spatial receptive field equivalent to figure 2.10 A. Instead of doing a contour plot, you can also make a color image (MATLAB function image) where different colors (or brightnesses) correspond to the bright or dark sensitive areas of the receptive field. The image function of MATLAB operates in a certain number range, so you probably have to scale your values to make the image visible.  You can use help image to see alternative functions for plotting as well (e.g. imshow, imagesc, etc; in Python, matplotlib has a similar imshow function) Use whichever one works best for you.

(b) Construct the temporal part of the corresponding receptive field as given in formula 2.29.   Use α = 1/(15ms). Convince yourself that you get a graph qualitatively equivalent to figure 2.14.

(c) Combine the spatial and the temporal parts of your receptive field (according to 2.a) to get a spa- tiotemporal field. Make a series of plots similar to figure 2.13 to see the temporal dynamics of your receptive field.

Now we want to construct the stimulus (a counterphase sinusoidal grating) according to equation 2.18.