Tutorial in Systems Neurophysiology and Modeling by techmaster


									Tutorial in Systems Neurophysiology and Modeling
Ning Qian and Vincent Ferrera, Instructors
Fall 2003

Exercise #1: Non-oriented, Center-Surround Receptive Field Model.

Goal: To implement and test a simple RF model in Matlab.

1. Matlab Basics

When you start Matlab, it should open with a Command Window. This window allows you to
enter individual commands at the prompt (>>). Commands can define variables. For example:

       >> x = zeros(100);

creates a 2D array (100x100) and fills it with zeros. Variables are initialized whenever they
appear on the left hand side of an expression. You do not need to specify the storage type or
size, as in other languages. Variables must be initialized before they are used on the right hand
side of an expression. When variables are initialized they get stored in the Workspace. You can
view the variables in the Workspace by typing:

       >> who

You can see the contents of a single variable by typing the variable name at the prompt without a
semicolon at the end:

       >> x

Commands can also call Matlab functions, e.g.

       >> plot(x);

To find out more about the function “plot”, type:

       >> help plot

If you want to create a program of more than a few lines, it soon becomes very tedious to enter
and re-enter commands one at a time. Therefore, Matlab allows you to create M-files. M-files
are plain text files that contain Matlab expressions and function calls. M-file names end with the
extension “.m”. Usually, the Command Window has a button that will invoke an M-file editor.
You can type multiple commands into the editor window and save them. Usually, there will be a
button in the editor window that allows you to save and execute the file in one shot. Otherwise,
you can execute the M-file by typing its name, sans extension, at the command prompt. E.g. if
you created the file “myprog.m”, you would just type:

       >> myprog
Important note: Each time you run an M-file, the variables in the workspace are not cleared.
This can lead to a situation where you spend several hours working on a program and it seems to
run fine. Then you come back to it later, after starting a fresh Matlab session, and you get a
bunch of error messages. What happened? Well, your program was probably using variables
that were in the workspace but not actually initialized by the program itself. There is a simple
way to avoid this. Start every M-file with the following command:


BTW, lines in M-files that start with % are treated as comments.

2. A simple RF model.

Retinal ganglion cells and LGN relay cells have, to a first approximation, circularly symmetric
receptive fields. This makes them sensitive to light-dark contrast, but not to contour orientation.
Hence, they are called “non-oriented” or “isotropic” receptive fields. We will start with a one-
dimensional model based on a difference-of-Gaussians (DOG).

Q: Why is one spatial dimension sufficient to model these RFs?

The equation for the RF is the following:
                                                             k            
                                                                                            
                                                             2                               2
                                                   x  xe                          x  xi
                             RF(x)  k e * exp                       *exp             
                                                   2             i           2           
                                                     e                          i          

Q: Which terms in this equation correspond to the amplitude, peak position, and spread (width)
of the two Gaussians?

To test the response of this model, we will construct a stimulus that is described by a cosine light
intensity function:

                                           I(x)  cos(2 * *s * x  )

Q: Which terms in this equation correspond to the spatial frequency and phase of the cosine?

To get started, we need to define a vector that specifies the range and grain over which we will
compute RF(x) and I(x). So, open an M-file and enter the following commands:


       xmin            =   -1.0;
       xmax            =   1.0;
       dx              =   0.01;
       x               =   xmin:dx:xmax;

Now, we’ll define the excitatory and inhibitory parts of the RF separately
       ke              =   1.0;
       xe              =   0.0;
       sige            =   0.125;
       RFe             =   ke * exp(-(x-xe).^2/sige.^2);

Note: The “.^” is not a typo. Matlab uses the “.” to distinguish simple math functions from
matrix functions. Hence, x.^2 squares each element of x; whereas x^2 multiplies the vector x
with itself resulting in a 2D matrix. The “.” can precede any math function (*, /, …).

Now the inhibitory part:

       ki              =   -0.25;
       xi              =   0.0;
       sigi            =   0.5;
       RFi             =   ki * exp(-(x-xi).^2/sigi.^2);

The full RF is constructed simply by adding the two parts:

       RF              = RFe + RFi;

RFe and RFi are vectors of identical length and can therefore be added to produce another
vector, RF. Notice how we were able to do this without a “for” loop. That’s the beauty of
vectorized programming.

Q: What physiological variable might be represented by the vector “RF”? Firing rate?
Membrane potential? Sensitivity? Is it “physiologically plausible” for the RF model to use
negative numbers? How might we be able to construct and inhibitory surround without using
negative numbers?

At this point, it should be fairly obvious how to construct the stimulus, but I’ll walk you through
it anyway.

       ws              = 1.0;
       phi             = 0.0;
       I               = cos(2 * pi * ws * x + phi);

Notice how pi is already defined. Gotta love those Matlab guys.

To determine the response of our model RF to the stimulus, we need to compute the inner
product of the two vectors, RF and I. This sounds tricky, but it is really very simple:

       resp            = dx*sum(I.*RF);

Notice again the use of the “.” Why is this not needed between “dx” and “sum”? Why am I
multiplying by “dx” in the first place? (Hint: “resp” is the integral of the function resulting from
the multiplication of I and RF. In this formulation, we are implicitly implementing the
rectangular rule for integration. How would you implement the trapezoid or parabola rules?).
At this point, you might want to start plotting some things. Try this:

       figure(1); clf;

To finish the exercise, try the following on your own:

1. Compute the full frequency response by varying the spatial frequency, ws. A good way to do
this is to declare ws as a vector and then use a “for” loop:

       ws        = 0.125.*2.^(0:0.5:6);
       resp      = zeros(1,length(ws));
       for i = 1:length(ws)
            I          = cos(2 * pi * ws(i) * x + phi);
            resp(i)    = dx * sum(I.*RF);
       plot(ws, resp);

Note: if you’re clever, you can accomplish this without a “for” loop by vectorizing the code.

2. Compute the frequency response separately for the excitatory (RFe) and inhibitory (RFi) parts
of the RF. How does the frequency response of the full RF compare to those of its parts?

3. Vary some of the parameters of the RF model and see how this affects the frequency response.
What constraints on the RF parameters are required so that the frequency response doesn’t
become negative?

4. Find a way to test the phase response of the RF model. Design your own stimuli. The fun
never ends!

To top