|Tutorial in Systems Neurophysiology and Modeling
Ning Qian and Vincent Ferrera, Instructors
Exercise #2: Linear and non-linear response properties of receptive fields.
For the purpose of this exercise, we will define a system as something that produces a measurable response, R, when presented with a stimulus, S. A system is said to be linear if it satisfies the following two properties:
1. The response scales proportionally to the input. If the amplitude of the input doubles, then the amplitude of the response also doubles. In general, if S -> R, then k*S -> k*R, where k is a constant.
2. The response to the sum of two (or more) stimuli is the sum of the responses to each stimulus independently. Thus, if S1 -> R1 and S2 -> R2, then S1 + S2 -> R1 + R2. This is also known as the Principle of Superposition.
Let’s start with the DOG model for RGCs from the last exercise and see if it is linear:
The stimulus we used last time was the following:
In reality, this stimulus is physically unrealizable, unless someone discovers “anti-photons” capable of creating negative light. So, let’s make it more realistic using the following equation:
Note: C ≤ 1. All we have done is to add a “mean luminance” around which the light intensity is modulated.
Q: Mean luminance is defined as the space-averaged value of I(x) and is usually denoted by Lmean. What is the mean luminance of the stimulus defined above? (Hint: replace the cos() term with its mean value and then do the math.)
Q: What are the maximum and minimum values of I(x) (i.e. Lmax, Lmin)? (Hint: replace the cos() term with it’s max or min value and do the math again.)
Q: If contrast is defined as the following:
Then substitute the expressions you derived for Lmax, Lmin and Lmean and derive an expression for contrast.
Note that we can re-write the expression for I(x) as:
Step 1 (scaling): Now, from the last exercise, you should have some code that computes the response of the RF model as a function of stimulus spatial frequency (the “frequency response”). Put this code inside a loop that varies Lmean. Plot the response of the RF as a function of spatial frequency for different values of Lmean. Next, make another loop that varies C. Plot the frequency response for different values of C.
To test the first condition for linearity, make an x-y plot of the response functions for two different Lmeans or contrasts. I.e. plot x vs. y where x = FR(ws, C1) and y = FR(ws, C2). (FR = frequency response; ws = spatial frequency). How does this plot demonstrate linear scaling? What is the slope of the line relating y to x?
Step 2 (superposition): To demonstrate superposition, set up two different stimuli I1 and I2. An easy way is just to use the expression for I(x) in both cases, but with different parameters, i.e. ws2 = 2*ws1. Now, calculated the response of the RF model to each stimulus independently:
R1 = dx * sum(RF .* I1);
R2 = dx * sum(RF .* I2);
Next, compute the response to the sum of the two stimuli:
R12 = dx * sum(RF .* (I1 + I2));
Confirm that R12 = R1 + R2. How does this demonstrate superposition? Does superposition hold for different values of ws? A nice way to test this is to make a loop that varies ws, so you end up with R1(ws), R2(ws) and R12(ws). Then you can plot R12(ws) with one symbol and R1(ws)+R2(ws) with another symbol (or line). The two curves should superimpose.
Step 3 (non-linearity): So far, so good, except we all know that the linear model is not entirely realistic. For one thing, as we increase Lmean, the response of the linear model keeps going up. As the stimulus becomes brighter, our model ganglion cells will be firing thousands, millions, or even billions of spikes per second. Those firing rates might be hard to sustain metabolically, not to mention the fact that we would need infinitessimal spike widths and a zero refractory period. Instead, let’s put a cap on the maximum firing rate by incorporating a saturating non-linearity after the linear filter. A good choice of non-linear function is the so-called Naka-Rushton (aka Michelis-Menten) function:
Rmax is the value of the function at y = infinity. k is the half-saturation constant. Convince yourself that R = Rmax/2 when y = k. The exponent, n, determines the steepness of the non-linearity. Note that I am using y as the argument because we already used x to represent space.
So, how do we incorporate this non-linearity into the RF model? The simplest way is just to use the response of the linear stage as the input argument for the non-linear stage. In other words:
Rlin = dx * sum(RF .* I);
Rnonlin = naka(Rlin, Rmax, k, n);
Notice that I have used a function called “naka” to represent the non-linearity. Where did this come from? Is it part of Matlab? Unfortunately, it is not. We have to create it ourselves. Here is the procedure:
1. Create a new M-file in the same directory as your current project. Name the new file “naka.m”
2. The first line of the new M-file should be:
function y = naka(x, Rmax, k, n);
3. The rest of the function can be written in one line, but let’s break it down for clarity:
temp1 = x.^n;
temp2 = k.^n;
y = Rmax .* (temp1 ./ (temp1 + temp2));
4. That’s it. Now when your main program calls naka with the appropriate arguments, the lines in naka.m will be executed and will return the value or vector y. Note: variables used in functions are not seen, nor are they available, to the main workspace, so it doesn’t matter if you use the same names in the function as you have in the M-file that calls the function.
Note: Matlab will search for any functions that are called by an M-file in the current PATH. The PATH is the current directory (folder) plus a list of other directories stored in a variable called “path”. To see the value of “path”, just type path at the command line.
You can now rewrite the M-file(s) you previously wrote to test RF linearity so as to incorporate the N-R non-linearity. Everywhere you see;
dx * sum(RF .* I);
Just replace it with:
naka(dx * sum(RF .* I), Rmax, k, n);
You can choose whatever values you like for Rmax, k, and n. Be careful so that the response stays positive for positive inputs. You might want to plot the function in a separate window. When you make the substitution, how does it affect the outcome of the linearity tests?