NeuronC Beginner's Course

Exercises

Section 1) How to run simulations

Section 2) Lessons

Section 3)

Section 4)





Section 1)

This beginner's manual is for the short course located in the Neuron-C distribution under "nc/ncourse".

Logging in

Turn on the computer, wait for the "login:" prompt, and login with your assigned login name. Once you have logged in, the computer prints a "%" at the beginning of every command line as a prompt to you that it has finished your previous command and is waiting for the next. Some simple (UNIX) commands are:
   %
   % ls                   ( directory listing )
   % cp file1 file2       ( copy file1 to make file2 ) 
   % logout               ( logs you off the computer )
The "" after these commands means that you enter the command into the computer by pressing the "return" or "enter" key.

Compiling nc

To compile "nc", run the "tar" command to retrieve the source code files:
    tar xvzf nc.tgz
This creates the folder "nc" and several subfolders. Go into the nc folder and run "make":
    cd nc
    make
This makes the simulator and several other useful programs. You must place them into a folder, such as ~/bin, your personal binary folder, or other bin folders that you have in your shell's path:
     cd bin
     cp nc stim plotmod vid ~/bin
 or e.g:
     su root
     cp nc stim plotmod vid neurc /usr/bin
Or else add the "nc/bin" folder to your path:
    vi ~/.cshrc
    add "~/nc/bin" to the path
Note that there are some useful programs in "nc/bin":
  bins.c      makes histogram
  column      cuts 1 column from an output file.r 
  ncd         runs nc to display morphology ("nc -d 1")
  ndv         runs nc to display morphology ("nc -d 1")
  ncv         runs nc with vid
  ncplot      displays file.r on vid ("plotmod $argv | xvid")
  replace     replaces string in files ("replace str1 str2 files")
  xvid        makes a larger window than vid
  etc.

Using the editor

There are many coding editors (kate, joe, vi) you can use to make a NeuronC program file. If you don't know how to use any of these editors, ask for help. Run the editor by typing its name (e.g. "kate") and a file name (the "%" is the computer's prompt):
For example, from the command-line prompt:

   % kate circuit1
Or, you can rely on a GUI editor like "soffice", or "openoffice". Run the editor by clicking on its icon, then click on "File" and then "New" to make a new file.

How to use editors

Vi

There are two modes, visual mode and command mode. In visual mode (where you initally get to) typing "i" gets you into insert mode. Any keys you type at this point will go into the file. To get out of insert mode, hit the key; this gets you back into visual mode. You can use the arrow keys on the keypad to move around. A colon ":" gets you into command mode.

   arrow keys      move around
   :wq             save and exit.
   i               insert text
   a               insert after char
   <Esc>     stop inserting text
  • more vi commands

    Emacs

    Enter text by typing normally. You can move using the arrow keys. Other "Ctrl" key sequences allow you to delete and cut and paste. Mouse button Function left Set point. middle Paste text. right Cut text into X cut buffer. Shift-middle Cut text into X cut buffer. Shift-right Paste text. Ctrl-middle Cut text into X cut buffer and delete it. Ctrl-right Select window, then split into two.

    Running Simulations

    After you type in your NeuronC program, you run the program by typing one of the following commands:

      neurc file.n                      Displays graph on screen.
      nc -t file.n                      Displays only numbers on output    .
                                        Used for "debugging" a model.
      nc file.n > file.r                Saves numbers in file. 
                                           Used for running long simulations.
      nc -v file.n | vid                (equivalent to "neurc file.n")
      plotmod file.r | vid              Displays graph on screen from file.
      plotmod file.r | vid -c >file.ps  Creates PostScript file 
      ps2pdf file.ps                    Creates .pdf file for emailing
    

    How the simulator works

    NeuronC is a language and simulator program for simulating complete physiology experiments on neural circuits. With the NeuronC language, you can describe and simulate a neural circuit as well as a complete stimulus and recording paradigm. Stimuli in NeuronC can be voltage or current clamps at single locations in the neural circuit or may be 2-dimensional patterns of light which cause simulated photoreceptors to transduce visual signals. Thus, NeuronC is especially designed for simulating visual neural circuits.

    Construction mode

    NeuronC has 2 distinct modes of operation, "construction" and "run". In the "construction" mode, NeuronC builds a neural circuit from programmed instructions, along with stimuli and recordings.

    In the "construction" mode, you can perform simple calculations, program loops, procedures and subroutines to define neural elements. You may organize these procedures to create a neuron or circuits of many neurons. You may also include stimuli and records in loops and procedures. "Construction" mode is familiar to computer programmers because it is similar to a language like "C" in which procedures are used to accomplish tasks.

    A simple calculation:
    
         print 54/7;                  // (prints simple calculations)
    
    A print statement within a loop:
    
         for (i=0; i<10; i++)
                print i, sqrt(i);     // (prints square roots)
    

    Run mode

    NeuronC's "run" mode stops "construction" mode and translates the neural circuit already defined into difference equations that are numerically integrated at high speed. During "run" mode, NeuronC numerically integrates voltage, runs optical and elecrode stimuli, and plots parameters. Further "construction" must wait until the run is completed.

         statement 1;   (these statements construct neural circuit)
         statement 2;
         statement 3;
         .
         .
         .
         step .05;      (stops circuit construction, runs 50 msec)
    

    Nodes, neural elements

    NeuronC describes circuits as a set of neural elements ("cable, "sphere", "synapse", etc.) connected to each other at nodes. A node is a common point at which neural elements connect electrically to each other. A node is not a neural element; it only defines connections between neural elements.

    NeuronC has several types of neural elements. These connect to either 1 or 2 nodes, and the nodes must be given when the element is defined. The types are:

    
       neural element      part of a neuron (see list below)
       node                point to place neural element, can have (x,y,z) location
    
       Element    Nodes    Description
    
       cable       2       defines multiple compartments.
       sphere      1       defines one compartment.
       synapse     2       defines connection between 2 nodes.
       gap junc    2       defines resistive connection between 2 nodes.
       photorecep  1       light transducer connected to node.
       load        1       resistor to ground.
       resistor    2       defines resistive conn. between 2 nodes (same as g.j.)
       cap         2       series capacitor between 2 nodes.
       gndcap      1       capacitor to ground; adds capacitance to node.
       batt        2       series battery between 2 nodes.
       gndbatt     1       battery to ground.
       chan        1       active set of channels in membrane.
    
    For neural elements that have 1 node connection, use the word "at" followed by the node number:
         at 20 sphere dia 5
    
         (at node 20, construct a sphere with diameter 5 microns)
    
    For neural elements that have 2 node connections, use the word "conn" followed by the first node number, followed by "to" and the second node number:
         conn 20 to 30 cable length 100 dia 2;
         (between nodes 20 and 30, connect a cable ...)
    
    Each neural element has a set of parameters which may be specified when the individual neural element is declared. If a parameter is not defined when the neural element is defined, it has a default value specified by a predefined variable. For instance the default value of membrane resistivity ("Rm") for passive dendrites is set by the "drm" variable. This variable may be set with an assignment statment (e.g. "drm = 3000") like any other variable.

    Predefined variables

    These variables are set to the following "reasonable" values, but you can change them with an assignment statement, e.g. endexp = 0.1;
    name          default   unit    meaning of variable
    ----------------------------------------------------------------
    timinc           1e-4    sec      time step for numeric integration (sec)
    endexp           .05     sec      time for end of experiment
    ploti            .001    sec      time step for plot increment (sec)
    
    vk              -.08     V        potassium battery voltage  
    vna              .04     V        sodium battery voltage  
    vcl             -.07     V        chloride leakage battery voltage  
    
    dcm                1    ufd/cm2   default membrane capacitance
    dri              200    ohm cm        default axial resistance 
    drg              5e6    ohm um2   default gap junction resistance 
    drm            40000    ohm cm2   default membrane resistance ohm cm2
    
    

    Units

    NeuronC assumes certain units for biophysical parameters, and cannot understand units if you type them in.
        diameters and lengths          microns
        membrane parameters            see units in above table.  
        time                           seconds, exception synaptic filters = msec
        voltages                       volts 
        current                        amperes
        resistance                     Ohms
        conductance                    Siemens
    

    Constructing neural circuits

    NeuronC defines a neural circuit with a sequence of statements much like any programming language. Each statement is either 1) a neural element statement, or 2) a standard programming statement (as in "C") which controls how many neural elements will be created. The analogy to a programming language is that NeuronC elements are like the printed output from a program; if a print statement is included inside a "for" loop, many printed lines will appear in the output. To create many neural elements, include their definition in a "for" loop in NeuronC.

    To build a neural circuit, you need to include at least one neural element statement and some stimulus and plot or print statements:

          conn 1 to 2 cable length 100 dia 1 rm 5000;   // constructs cable
          stim node 1 vclamp -.01 start .05 dur .01;    // adds stimulus
          plot V[1];                                    // adds plot
          plot V[2];
          step .1;                                      // runs experiment
    
    This simple NeuronC program builds a cable and defines a stimulus of a voltage clamp to minus 10 mv, starting at 50 msec into the experiment and stopping after 60 msec. The voltage at each end of the cable is plotted. The plotting stops at 0.1 seconds.


    Section 2)

    Lesson 1) Synaptic input to neuron.

    Background:

    A neuron receives signals through synaptic connections from other neurons. The synapse evokes a post-synaptic potential (PSP) by opening or closing channels in the post-synaptic membrane.

    Model 1:

    Construct a presynaptic terminal, a synapse, and a segment of a neural dendrite. Stimulate the presynaptic terminal with a voltage clamp, and record the post-synaptic potential. Use the editor to make a file (or use existing "mod1.n"):

    (Enter the following lines in the editor:)
    
    (c o d e:)                               ( c o m m e n t s: )
    
    ploti=1e-4;
    
    pre = 100;                              /* convenience definition */
    at pre (0,10) sphere dia 1;             /* pre = node for presynaptic term. */
    
                                           /* make synapse */
    conn pre to 1 loc (0,0) synapse expon 5 thresh -.05 maxcond 5e-9 vrev 0;
    
    conn 1 to 2 loc (120,0) cable length 120 dia 1 rm 5000;   /* make cable */
    
    stim node pre vclamp -.045 start .005 dur .01;       /* voltage clamp */
    
    plot V[1] max 0 min -.07;               /* plot the voltage at node 1 */
    plot V[pre] max 0 min -.07;
    
    if (disp) {
      display center (0,0,0);
      display size (250);
      display matching [-1][-1][-1];
      display node matching [-1][-1][-1] color 7;
      exit;
    };
    
    endexp = .05;                           /* x axis length on output plot */
    step .05;                               /* length of simulation in sec. */
    

    Next, run the program "mod1.n" by typing:

    
         % neurc mod1
    

    NeuronC runs a voltage-clamp stimulus to open synaptic channels at 5 msec, and closes the channels 10 msec later. NeuronC plots a graph showing voltage and current vs. time on the screen. Each "plot" statement makes a different color line on the graph.

    Notes:

    1) The green line represents the presynaptic potential, or -45 mv in this model. After the voltage clamp is turned off at 15 msec, the green line decays to resting potential (-70 mV).

    2) The synapse is depolarizing. The blue line represents the time-course of the post-synaptic potential (PSP), which peaks at about 20 mV above resting potential.

    3) The post-synaptic membrane potential does not reach the synaptic reversal potential (0 mv) because the membrane conductance (1/resistance) shunts current from the synaptic conductance. The post-synaptic cable segment is quite long and has a large surface area, so its membrane conductance is relatively high.

    4) The voltage decay of the post-synaptic membrane reflects membrane properties (tau), as well as synaptic properties (default synaptic delay dfta = dftb = 0.2 msec).

         tau, the time constant = Rm   * Cm
                                = 5000 * 1e-6
                                = 5 msec
    

    5) It may be helpful and fun to try changing some of the parameters. What happens if you change synaptic "thresh", or the "vclamp" value?



    Lesson 2) Electrotonic decay of voltage in cable

    Background: Neurons in the retina have special shapes (morphology) tuned for their function. The morphology of some neurons imparts a "weighting function" to synaptic inputs that originates in "electrotonic decay".

    A long neural dendrite approximates an infinite cable which causes a voltage at one point to decay exponentially along the dendrite. An important parameter for a neural cable is lambda, or "length constant", which is the distance along the cable that is required for the voltage to decay to 1/e or 37%.

      Definitions:
    
      Rm = membrane resistivity (ohm*cm2)
      rm = membrane resistance  (ohm*cm)
      rm = Rm / 2 * PI * radius
    
      Ri = axial resistivity (ohm*cm)
      ri = axial resistance  (ohm/cm)
      ri = Ri / PI * radius2
    
      lambda = sqrt (rm/ri)
      tau = Rm * Cm  =  rm * cm
    

    Problem: How far along a dendrite must a neural signal travel for its voltage to decay to 37% (1/e)?

    Model 2: Construct a cable long enough to observe exponential decay and stimulate it by injecting current with the synapse created in "mod1" at one end. Record voltages in the cable at various locations. Use the editor to make a file containing the following lines (or use "mod2t.n"): % vi mod2t

    ( c o d e )                               ( c o m m e n t s)
    
    ploti = 1e-4;
    drm = 5000;                               /* default Rm = 5000 Ohm-cm2 */
    pre = 100;                                /* pre-synaptic node = 100 */
    
    at pre loc (0,10) sphere dia 1;
    conn pre to 1 loc (0,0) synapse expon 5 thresh -.05 maxcond 5e-9 vrev 0;
                                       /* thresh 50 mv, max con 5e-9 S, etc. */
    
    conn 1 to 2 loc (120,0) cable length 120 dia 1;       /* construct cable */ 
    conn 2 to 3 loc (240,0) cable length 120 dia 1;
    conn 3 to 4 loc (360,0) cable length 120 dia 1;
    conn 4 to 5 loc (480,0) cable length 120 dia 1;
    conn 5 to 6 loc (600,0) cable length 120 dia 1;
    
    stim node pre vclamp -.01 start 0 dur .01;     /* voltage clamp */
    
    
    plot V[1] max 0 min -.07;                 /* max and min define  */
    plot V[2] max 0 min -.07;                 /*   the Y-axis of graph */
    plot V[3] max 0 min -.07;
    plot V[4] max 0 min -.07;
    plot V[5] max 0 min -.07;
    plot V[6] max 0 min -.07;
    
    plot V[pre] max 0 min -.07;
    
    if (disp) {
         display center (0,0,0);
         display size (1000);
         display matching [-1][-1][-1];
         exit;
       };
    
    endexp = .015;                            /* defines length of X-axis */
    step .015;                                /* run model for 15 msec */
    
    
    Run the program mod2t by typing:
         % neurc mod2t.n 
    
    The model consists of a synapse and 5 cable segments, each 120 microns long, connected together between 5 nodes. This makes a continuous cable 600 microns long. It is constructed from 5 segments to allow recording the voltage at several nodes.

    Notes for mod2:

    1) Verify that the time constant (time to fall to 63% of maximum voltage) of the presynaptic node is about 5 msec. This can be calculated as Rm= 5000 Ohm-cm2 * Cm = uF/cm2.

    2) Verify that the voltage on the red line (node 4), after it reaches an equilibrium value at 10 msec is near 37% of maximum on the blue line (node 1). Thus, lambda is about 360 microns (3 * 120) for the cable.

    3) Note that the stimulus takes a few msec. to reach the distal end of the cable. The conduction velocity may be calculated as:

         velocity = tau/lambda
    
    Next, run the same model with a "space plot" instead of a "time plot". The space plot is "mod2s.n":

    
       drm = 5000;                               /* default rm = 5000 */
       pre = 100;                                /* pre-synaptic node = 100 */
    
       at pre loc (10,0) sphere dia 1;
       conn pre to 1 (0,0) synapse expon 5 thresh -.05 maxcond 5e-9 vrev 0;
    				   /* thresh 50 mv, max con 5e-9 S, etc. */
    
       conn 1 to 2 loc (120,0) cable length 120 dia 1;       /* construct cable */
       conn 2 to 3 loc (240,0) cable length 120 dia 1;
       conn 3 to 4 loc (360,0) cable length 120 dia 1;
       conn 4 to 5 loc (480,0) cable length 120 dia 1;
       conn 5 to 6 loc (600,0) cable length 120 dia 1;
    
       if (disp) {
         display center (0,0,0);
         display size (1000);
         display matching [-1][-1][-1];
         exit;
       };
    
       stim node pre vclamp -.01 start 0 dur .01;
    
       step .01;                            /* run model for 10 msec */
    
       graph X max 600 min 0;               /* set X-axis scale */
       graph Y max 0 min -.07;              /* set Y-axis scale */ 
       graph init;                          /* draw X- and Y-axes */ 
    
       graph (0,  V[1]);                    /* graph volts vs. distance */
       graph (120,V[2]);
       graph (240,V[3]);
       graph (360,V[4]);
       graph (480,V[5]);
       graph (600,V[6]);
    
       gpen (7);                            /* make labels for graph */
       gmove (0.02,0.85);
       gtext ("Volts");
       gmove (0.45,0.01);
       gtext ("microns");
    
    Run the program mod2s by typing:
         % neurc mod2s.n
    
    Notes for mod2s:

    1) Remember that the signal is measured by subtracting the resting potential (-70 mV) from the voltage measured in the cable. Verify from the plot that lambda (37% of the maximum signal) is about 350 microns. In this model, Rm=5000 ohm-cm2, Ri=100 Ohm- cm, so lambda is calculated to be 354 microns for a 1 micron diameter cable.

              lambda  =  sqrt (rm/ri)  =  sqrt ((Rm / Ri) * (dia / 4)) 
    
    2) In an infinite cable, electric current flows away from the source of current, attenuated by leakage through membrane conductances. This produces exponential decay, which causes the voltage to decrease by a constant ratio over each segment of cable.

    3) Since the cable is not infinite, its decay is not exactly a simple exponential. Verify that the voltage decay in "mod2s" is not exponential near the distal end (nodes 4-6). At the distal end (away from the synapse), no current leaks through the end of the cable. Therefore the voltage does not decay there, and the voltage vs. distance plot approaches a level slope.

    Effect of cable diameter

    Next, shrink the cable diameter to diameter 0.1 (micron) instead of 1. To run a different model, you should copy "mod2s" into a new file and modify the new file:
         % cp mod2s mod2ss 
    
    Run the simulation with the smaller diameter, and compare it with the original:
       % neurc mod2ss mod2s
    
    Notes for mod2ss:

    1) Note that the post-synaptic potential is larger. This is because the cable's resistance (both rm and ri) is higher with a smaller diameter. Now the voltage at node 1 approaches the reversal potential for the synapse (near 0 mV).

    2) Verify that decreasing diameter by a factor of 10 decreases lambda by a factor of about 3 (square root of 10). Lambda (distance for decay to 37% of max) is now about 100 microns.

    3) Verify that the voltage decay is exponential between segments (voltage decreases by roughly a constant factor at each segment) for the smaller 0.1 micron diameter cable. Voltage decay is more exponential because the cable is several "lambdas" long, so the distal end does not affect the proximal end much.

    Next, change "drm" (default Rm) from 5000 to 50000 and observe the corresponding changes. Note that this increases lambda to near 350 um, but has little effect on the PSP amplitude. This is because the synaptic output is already near its "driving" potential, and the post-synaptic potential is said to be "saturated".

    Note: A change in either Rm or diameter causes an equivalent change in lambda.




    Lesson 3) Electrotonic decay in dendrite with spines

    Background: Many neurons receive input at "spines" which are short dendritic branches with extremely fine necks (about 0.1 micron). These spines restrict the flow of electric current because they are so narrow. The function of spines for retinal horizontal cells may be to shape their spatial responses.

    Problem: How does a spine affect a synaptic input signal?

    Model 3: Construct a shorter cable than in the previous model, and attach spines at the 2 ends. Use the editor to create a new file (or use existing "mod3.n"):

       % vi mod3
    
      drm = 5000;
    
      pre = 100;                           /* convenient node definitions */
      spine1 = 10;
      spine2 = 20;
    
      numsegs = 4;                         /* number of segments */
      seglen = 60;                         /* length of segments */ 
      spinelen = 5;
      totlen = numsegs * seglen + 2 * spinelen;
    
      at pre sphere dia 1;
      conn pre to spine1 synapse expon 5 thresh -.05 maxcond 2e-9;
    
      conn spine1 to 1 cable length spinelen dia 0.1; /* spine at beg. of cable */
    
                                            /* "for" loop saves typing */
      for (i=1; i<=numsegs; i++)           /* "for segments 1 to numsegs */
            conn i to i+1 cable length seglen dia 1;
    
      conn spine2 to 5 cable length spinelen dia 0.1; /* spine at end of cable */
    
       stim node pre vclamp -.01 start 0 dur .01;
    
       step .01;                                 /* run model for 10 msec */
    
       graph X max totlen min 0;                 /* commands to scale graph */
       graph Y max 0 min -.07;
       graph init;
    
       graph pen (12);
       graph (0, V[spine1]);                     /* post-synaptic potential */
       for (i=0; i<=numsegs; i++) {
          if (i==1) graph pen (2);                 /* i==1 -> (double equals) */
          graph (i*seglen+spinelen, V[i+1]);     /* graph volts vs. distance */
       };
       graph pen (1);
       graph (totlen, V[spine2]);                /* output signal */    
    
       gpen (7);                                 /* labels for graph */
       gmove (0.02,0.85);
       gtext ("Volts");
       gmove (0.45,0.01);
       gtext ("microns");
    
    
    Run the program mod3 by typing:
         % neurc mod3.n
    
    The model consists of a synapse connecting to the first spine that in turn connects to a cable of the same diameter as in "mod2". In a retinal horizontal cell, each spine is also a site for the neuron's output. In this model, the neuron's output is located at the tip of the second spine.

    Notes for mod3:

    1) Verify that the voltage difference across the length of spine 1 (the red line) is about 20 mV, and that the decay along the larger dendrite (the green line) is about 5 mV.

    2) The axial resistance "ri" in the spine neck is much greater than the axial resistance in the main dendrite. Therfore, most of the voltage drop occurs in the spine neck.

    3) To understand this reasoning, imagine injecting water into a leaky garden hose with a hypodermic needle.

    4) The value of the axial resistance "ri" is more affected by reduced cable diameter than rm, because ri is inversely proportional to the square of diameter:

              rm = Rm / 2 * PI * radius 
              ri = Ri / PI * radius2
    
    When the synaptic input current passes through a spine to enter a neuron, the loss of current through the spine's membrane resistance is insignificant compared to the voltage drop across ri.

    5) The output spine (spine 2 in "mod3", the blue line) does not cause a large voltage drop. This is because its resistance is high and therefore only a small current enters it. Thus, the voltage drop across ri is also small. Imagine the minor loss of pressure when filling a hypodermic needle with a garden hose.



    Lesson 4) Model of Horizontal cell

    Background: The retinal horizontal cell is roughly circularly symmetric, having 5 to 8 branched dendrites emanating from a large cell body (soma). Input to horizontal cells comes from photoreceptor synapses distributed throughout the dendritic tree, and their output synapses are also distributed throughout the dendritic tree. Morphology varies with type of horizontal cell (dendrite length/diameter/ branching pattern) but it is likely that all horizontal cells function in a similar way.

    The effect of a synaptic input on a neuron's signal can vary markedly, depending on the location of the cell's input and output. A synapse has a relatively high resistance compared with the total resistance of the neuron. If Rin (input resistance of the neuron) seen by the input synapse is high, the synapse drives the postsynaptic membrane easily but the membrane voltage may saturate, causing a non-linear signal transformation. However, if Rin seen by an input synapse is low, the synapse cannot drive its post- synaptic membrane well. For this reason, the electric signal from a single synapse within a neuron does not transfer well from dendrites to the low-resistance soma, because the low resistance shunts the signal.

    If a neuron's output is derived from this low resistance soma (e.g. as in a ganglion cell), then the output signal collected from an individual synaptic input is reduced by this low resistance. However, the electric signal collected by the soma from many input synapses does reflect the average synaptic input (e.g. ganglion cell, horizontal cell). This average signal transfers well to distal dendrites, if they are thick and have low axial resistance (e.g. horizontal cell), because dendrites usually have a higher input resistance than the soma.

    If an output synapse is located in a tapered distal dendrite with relatively low axial resistance (e.g. horizontal cell), the soma's signal can effectively drive the output synapse. If, however, the output synapse is located near input synapses in a fine distal dendrite with high axial resistance (e.g. wide-field amacrine cells), its signal may reflect nearby input synapses rather than the soma's signal.

    Thus, the effectiveness of synaptic transfer from input to output reflects both dendritic morphology and relative position of the synaptic inputs and outputs.

    Problem: How does anatomy and membrane biophysics of horizontal cell affect voltage response at different parts of cell?


    Model 4) Add a soma to model 3 and test the effect of varying soma location with respect to the input and output spines. Copy the previous modeling file ("mod3.n") to make another file:

       % cp mod3.n mod4.n
    
    Edit the file ("mod4.n"):
       % vi mod4.n
    
    (Modify "mod4" in the editor to look like this:)
    (Soma added!)
    ( c o d e )                             ( c o m m e n t s ) 
    
       drm = 5000;
    
       pre = 100;                           /* convenient node definitions */
       spine1 = 10;
       spine2 = 20;
    
       numsegs = 4;                         /* number of segments */
       seglen = 60;                         /* lengths of segments */ 
       spinelen = 5;
       totlen = numsegs * seglen + 2 * spinelen;
       soma = numsegs+1;                    /* soma is connected to node 5 */ 
    
       at pre sphere dia 2;
       conn pre to spine1 synapse expon 5 thresh -.05 maxcond 2e-9;
    
       conn spine1 to 1         cable length spinelen dia 0.1;
       conn spine2 to numsegs+2 cable length spinelen dia 0.1;
    
       for (i=1; i<=numsegs; i++)
            conn i to i+1 cable length seglen dia 1;
       conn i to i+1 cable length seglen dia 1;
    
       at soma sphere dia 20;
    
       stim node pre vclamp -0.01 start 0 dur 0.01;
    
       step 0.01;                                /* run model for 10 msec */
    
       graph X max totlen min 0;                 /* commands to scale graph */
       graph Y max 0 min -.07;
       graph init;
    
       graph pen (12);
       graph (0, V[spine1]);                     /* post-synaptic potential */
       for (i=0; i<=numsegs; i++) {
          if (i==1) graph pen (2);
          graph (i*seglen+spinelen, V[i+1]);     /* graph volts vs. distance */
       };
       graph pen (1);
       graph (totlen, V[spine2]);                /* output signal at spine 2*/    
    
       gpen (7);                                 /* labels for graph */
       gmove (0.02,0.85);
       gtext ("Volts");
       gmove (0.45,0.01);
       gtext ("microns");
    
    Run the program mod4 by typing:
         % neurc mod4
    
    Run the simulation:
         % neurc mod4
    
    Notes for mod4:

    1) Verify that the soma increases voltage drop along spine 1 (red line) and the main dendrite, but does not affect voltage drop in spine 2 (blue line).

    2) Verify that voltage decay is now greater than exponential, that is, voltage decay is now relatively linear in each segment, not constant in ratio.

    Move soma to proximal location:

    Next, copy "mod4" and make a new version "mod4a" in which the soma is located at the branchpoint for the input spine (node 1). This time, run the 2 versions together so you can compare them:

              neurc mod4 mod4a
    

    Note that the effect of locating the soma directly at the input spine is to increase voltage drop across the spine neck. However, now there is a slightly larger signal at the soma, because the signal is not reduced by voltage drop due to the axial resistance (ri) along the main dendrite. But for the signal at spine 2, the soma location doesn't matter much.

    Move soma halfway:

    Make another version of "mod4", say "mod4b". Modify "mod4b" to locate the soma halfway between spine 1 and spine 2 (node 3), and run the model again. Run the 3 versions of mod4 together:

         neurc mod4 mod4a mod4b
    
    Note that the effect of placing the soma halfway along the main dendrite is to reduce the length of the dendrite and slightly increase signal at the soma. The signal arriving at spine 2 is only slightly reduced by traversing the distal part of the main dendrite (between soma and spine 2).

    The effect of additional horizontal cell branches on the signal from one synapse can be simulated by the soma. The soma's effect is to shunt electric current and therefore to increase voltage drop along the axial resistance of the dendrites and spines.

    To see a display of the neuron in mod4.n, you can add display statements. However to correctly display the neuron, all of its nodes must be given locations with the "loc" parameter. The mod4d.n model is the same as mod4.n but has the node locations defined and includes the display statements. Run mod4d.n like this:

         neurc -d 1 mod4d.n
    
    To see a display of the node numbers, you can run mod4d.n like this:
         neurc -d 9 mod4d.n
    



    Lesson 5) Model of photoreceptor network

    Background: Photoreceptors in all vertebrates are coupled by gap junctions, which spread electrical signals laterally. Gap junctions therefore couple photoreceptors into a large 2-dimensional syncytium (coupled neural array). In effect, this coupling "blurs" visual signals after they are transduced by photoreceptor outer segments. This blur is an essential part of retinal signal processing. Its function is to filter out spatial high-frequency signals before they traverse the first chemical synapse.

    Problem: What is the effect of electrical coupling in a small array of cones on the signal from one cone?

    Model 5) Create a small array of cones (modeled here as single compartments) and couple them into a 2- dimensional array with gap junctions. Also create an isolated cone for comparison. Stimulate the array with a small spot over the central cone, and record the responses of a line of cones. With the editor, create a new file (or use existing "mod5t.n"):

       % vi mod5t.n
    
    Enter the following lines:
    (c o d e)                                 /* c o m m e n t s */
    
    endexp = .15;                             /* x axis length for plot */
    drm = 5000;
    
    for (row=0; row<35; row+=5)               /* create cone array */
          for (col=0; col<35; col+=5) {       /* nested "for" loops */
              at [row][col] cone (row,col) maxcond 1e-9; 
              at [row][col] sphere dia 10;
          };
    
    
    drg = 1;                                /* calibration for gsiz in Siemens */ 
    gsiz = 1e-8;                           /* 1 connexon = 100 pS */
    
    for (row=0; row<30; row+=5)                      /* connect gj's between rows */
          for (col=0; col<35;  col+=5) {
              conn [row][col] to [row+5][col] gj gsiz;
          };
    
    for (row=0; row<35; row+=5)                      /* connect gj's between columns*/
          for (col=0; col<30; col+=5) {
              conn [row][col] to [row][col+5] gj gsiz;
          };
    
    cent = 15;
    
    stim spot 1 loc (cent,cent) inten 1e7 start .01 dur 1e-3;
      /* spot 1 um diam;        inten 1e7 phot/um2/sec; dur 1 msec */
    
    plot V[0][0]     max -.031 min -.035;         /* plot corner response */
    for (dist=0; dist<35; dist+= 5) {            /* plot line though center */
       plot V[cent][dist] max -.031 min -.035;
    };
    
    plot L[cent][cent]   max 5e7  min 0;         /* plot light flash */
    
    step .15;
    
    Run the simulation:
         % neurc mod5t
    
    Notes: 1) Note that mod5t gives a plot of the cones' response voltage versus time. The cones take about 10 msec. to initially reach an equilibrium potential. Compare the response of the central cone with the others in the array. The central cone's response is considerably reduced because of lateral current flow through the gap junctions.

    2) Note that the response of cones near the edge of the array is lower in amplitude and delayed in time, compared with the central cone's response.

    Next, modify "mod5t.n" for a spatial plot. Copy mod5t.n into a new file "mod5s.n" and modify it to look like this:

    (mod5s should look like this:)

    drm = 5000;
    
    for (row=0; row<35; row+=5)                  /* create cone array */
          for (col=0; col<35; col+=5) {
              at [row][col] cone (row,col) maxcond 1e-9; 
              at [row][col] sphere dia 10;
          };
    
    drg = 1;                               /* calibration for gsiz in Siemens */ 
    gsiz = 1e-8;
    
    for (row=0; row<30; row+=5)                      /* connect gj's between rows */
          for (col=0; col<35;  col+=5) {
              conn [row][col] to [row+5][col] gj gsiz;
          };
    
    for (row=0; row<35; row+=5)                      /* connect gj's between columns*/
          for (col=0; col<30; col+=5) {
              conn [row][col] to [row][col+5] gj gsiz;
          };
    
    cent = 15;
    stim spot 1 loc (cent,cent) inten 1e7 start .01 dur 1e-3;
    
    step .07;
    
    graph X max 30 min 0;                   /* commands to scale graph */
    graph Y max -.031 min -.035;
    graph init;
    
    graph pen (12);
    graph (0, V[0][0]);                     /* corner cone */
    for (dist=0; dist<35; dist += 5) {      /* plot line though center */
          graph (dist, V[cent][dist]);      /* graph volts vs. distance */
    };
    
    gpen (7);                               /* labels for graph */
    gmove (0.02,0.85);
    gtext ("Volts");
    gmove (0.45,0.01);
    gtext ("microns");
    
    
    Remember that mod5s takes about 30 seconds to run. Be patient while the cone responses are recorded at their maximum.

    By recording from a line of cones, "mod5s" gives a spatial plot similar to a receptive field map of a single cone. Note that the spatial response of the cone array looks somewhat exponential (actually a modified Bessel function). The space constant for the array's response is approximately 2 cones or 10 um distance.

    The space constant can be varied by changing the gap junction conductance. Try increasing and decreasing the conductance by factors of 2, and then try increasing by a factor of 1000 to simulate closed gap junctions.

    In real experiments, a stimulus always has some optical blur. To see the effect of optical blur, try increasing the size of the stimulus light spot to 15 microns in diameter. This causes a larger response by stimulating several cone photoreceptors, and changes the shape of the cone receptive field.

    To see the effect of totally isolating the cones, you can try copying "mod5t" to make "mod5ti", and add the following lines to "mod5ti":

    row = 100;
    col = 100;
    cent = 10;
    
    at [row][col] cone (cent,cent);              /* make single isolated cone */
    at [row][col] sphere dia 5;
    
    plot V[row][col] max -.02 min -.07;          /* plot the isolated cone */
    
    
    But perhaps an easier way to look at "isolated" cones in this type of model is to stimulate all of them with a spot large enough to "cover" all of them. Try increasing the size of the spot to 100 microns. In this case, there is no lateral coupling because each cone has an identical response. This approach requires more computation, though!




    Lesson 6) Model of wide-field amacrine cell.

    Background: Many amacrine cells have a distinctive morphology consisting of extremely long, fine dendritic processes emanating from the soma. The dendritic processes enlarge to form small beads or "varicosities" every 5 to 10 um. Each varicosity typically contains one input and one output synapse. It is likely that the space constant of the fine dendritic processes is relatively short (i.e. less than 50 um); therefore the function of this type of amacrine cell is probably local processing.

    Problem: What is the space constant of a fine dendritic process that contains varicosities and synaptic inputs?

    Model: Construct a fine (0.1 um diameter) dendritic cable with varicosities (4 um diameter) spaced every 10 um. Locate a synaptic input at one varicosity. Stimulate the synapse, and record voltage in the local varicosity as well as in other distal varicosities. With the editor, create a new file ("mod6t.n"):

       % vi mod6t.n
    
    Enter the following lines:
    /* amacrine dendrite with varicosities */
    /*  temporal plot */
    
    drm = 5000;
    endexp = 0.02;
    
    nsegs = 10;
    seglen = 10;
    totlen = seglen * nsegs;
    varicos = 4;
    synap = 100;
    
    for (i=1; i<=nsegs; i++) {
      at   [i][synap] sphere dia 3;                 /* presynaptic terminal */
      conn [i][synap] to [i] synapse expon 5 maxcond 5e-9 thresh -0.04 vrev -0.01;
      at i sphere dia varicos;                      /* the varicosity */
      conn i to i+1 cable length seglen dia 0.1;    /* the fine cable */
      if (i==5) stim node [i][synap] vclamp -.025 start 0 dur 0.01
    /*else      stim node [i][synap] vclamp -.030 start 0 dur 0.05; */
    };
    at i sphere dia 5;
    
    for (i=1; i<=nsegs; i++) {
      plot V[i] max -0.01 min -0.07;
    };
    step .02;
    
    
    Now, exit the editor. Run the simulation:
         % neurc mod6t
    
    Notes

    This model consists of 10 varicosities and fine cable segments connected as beads on a string, connected on one end to a small soma. A synapse is located at varicosity 5, which is near the center of the cable. Verify that the voltage at this node is greater than other nodes, and that it reaches an equilibrium potential near 10 msec, when the stimulus gets turned off.

    Now, try the space plot. Make another file "mod6s" by copying "mod6t". Erase the "plot" lines and put in the "graph" lines:

    /* amacrine dendrite with varicosities */
    /*  spatial plot */
    
    drm = 5000;
    
    nsegs = 10;
    seglen = 10;
    totlen = seglen * nsegs;
    varicos = 4;
    synap = 100;
    
    for (i=1; i<=nsegs; i++) {
      at   [i][synap] sphere dia 3;                 /* presynaptic terminal */
      conn [i][synap] to [i] synapse expon 5 maxcond 5e-9 thresh -.04 vrev -.01;
      at i sphere dia varicos;                      /* the varicosity */
      conn i to i+1 cable length seglen dia 0.1;    /* the fine cable */
      if (i==5) stim node [i][synap] vclamp -0.025 start 0 dur 0.01;
    /*  else    stim node [i][synap] vclamp -0.030 start 0 dur 0.05; */
    };
    at i sphere dia 5;
    
    graph X max totlen min 0;
    graph Y max -0.01 min -0.07;
    graph init;
    gpen (7);                                 /* labels for graph */
    gmove (0.02,0.85);
    gtext ("Volts");
    gmove (0.45,0.01);
    gtext ("microns");
    
    stim node [seglen][synap] vclamp -.045 start 0.01 dur 0.02;
    
    graph pen (2);
    stepsiz=0.001;
    for (t=0; t<0.02; t+=stepsiz) {
       step stepsiz;
       graph restart;
       if (t>=.01) graph pen (4);
       if (t==0.01) stim node [i][synap] vclamp -0.045 start 0.01 dur 0.02;
       for (i=0; i<=nsegs; i++) {
          graph (i*seglen, V[i+1]);           /* graph volts vs. distance */
       };
    };
    
    Now, exit the editor. Run the simulation:
         % neurc mod6s
    
    Notes:

    1) Note that the dynamic range of the signal is about 55 mV. This results because the synapse is turned completely on and then off.

    2) Verify that lambda is about 3 cable segments long (30 microns). The green graph lines represent the voltage along the dendrite during rising stimulus, and the red lines represent the falling stimulus. This is longer than lambda for a 0.1 micron dendrite (see Model 2) which is near 100 microns. The effect of the varicosities is to act as "radiators" for the electric current. The effect of the varicosities is to decrease membrane resistance for the fine dendrite.

    To investigate the effect of synaptic conductance at all the varicosities, locate synaptic inputs at all varicosities. This time, stimulate the original synapse as before, but stimulate the others as well.

    Construct additional extra synapses by removing the comments around the second "stim" statment in "mod6s". Also remove the semicolon at the end of the first stim statement. The two stim statments should look like this:

         if (i==5) stim node [i][synap] vclamp -.025 start 0 dur .01
         else      stim node [i][synap] vclamp -.030 start 0 dur .05;
    
    Also, modify the Y min value for the graph:
         graph Y max -.01 min -.04;
    
    Then run the model:
         neurc mod6s
    
    Notes:

    1) The extra synapses are activated throughout the duration of the model and therefore do not directly add to the signal. However, the effect of the extra synapses is much like a varicosity: an activated synapse can reduce membrane resistance and therefore passively reduce signal spread.

    2) Verify that the dynamic range (from activated on to inactivated) of the signal is only about 8 mV now. This occurs because the decreased membrane resistance shunts electric current laterally.

    3) Note that after the stimulus is turned off, the red line shows an "inverted V" which is the result of turning off the signal synapse at varicosity 5 while the other synapses remain activated.

    4) Verify that lambda is now about 10 microns. This much reduced lambda implies that under certain conditions a) synaptic connections can control lambda in fine dendritic trees, and b) the signal in each varicosity can be effectively controlled by its input synapse, irrespective of signals in neighboring varicosities.


    Section 3)

    To learn about spiking neurons, you can run the "scomp" series of models in nc/tcomp. They are single-compartment models with several membrane channels that together form a spike generator.
    scomp_fm.n       single compartment model based on Fohlmeister & Miller (1997a)
    scomp_fm2.n      single compartmen model based on Fohlmeister & Miller (1997a), modified
    scomp_a.n        shows spike adaptation.
    scomp_b.n        similar to scompa.n, tuned to 37 deg
    scomp_c.n        derived from scompa.n @ 22 deg C subst "ndensity" for "density"
    scomp_d.n        similar to scompc.n, tuned to 37 deg
    scomp_persis.n   sets up persistent Na current, controlled by ifactor
    scomp_ad*.n      series of models tuned to 37 deg, used to develop Mark van Rossum's (2003) model
    scomp_ad6.n      contains several experiments, many params, can be run at diff temperatures
    

    Scomp models: single compartment spike generator

    The scomp models look very simple. They consist of just one compartment with several membrane channels, and they generate spikes according to the current injected into the compartment. However, the simplicity of the model can be misleading at first. When more than 1 or 2 membrane channel types are included in a compartment, the complexity factor rises exponentially, so that it can be very difficult to predict with intuition what will happen, or even to understand the role of each channel.

    The Fohmmeister & Miller (1997) model

    The first scomp model is "nc/tcomp/scomp_fm.n". This model is a close match to the model presented in Fohlmeister & Miller (1997a), in which a single compartment model of a spike generater from salamander ganglion cells is described.

    The model starts with default values (before line 61, "x = setvar()"), which sets parameters from the command lines, possibly overriding the default values. Then some of the kinetic rate functions for activation, deactivation, inactivation, and recovery from inactivation are specified. Although these all have default definitions in the simulator (which vary by channel type), they are specified explicitly in scomp_fm.n so you can see the original definitions that Fohlmeister & Miller (1997) used. This also allows you to change the rate functions so you can see the effect.

    Next, the single compartment cell is defined:

    at soma sphere dia 25 vrest grest vrev -.06;
    
    This specified a cell with 25 um diameter soma, with Rm = 50000 (set above by drm=50000), and a battery voltage for Rm of -60 mV.

    Next, the channels are added to the soma, each with its type specified, its density (in S/cm2), possibly with noise properties (chnoise), and a name for the channel so it can be plotted as a trace on the output graph.

    The "onplot.m" file contains functions to run at plot time. The "onplot()" procedure in "onplot.m" runs at the plot interval but runs before any of the other plots are generated. In this case, it calls the spikplot() procedure which calculates the instantaneous spike frequency from the spike interval, and it also calls the calc_totcur() function in the scomp_fm.n model, which calculates the total channel current in the model.

    Next the experiments are defined. Each one is called when its name is found in the "expt" parameter, set by default to "spike". You can change the experiment that runs by changing the value of "expt" on the command line: neurc --expt spike scomp_fm.n or neurc --expt isi scomp_fm.n or neurc --expt f/i scomp_fm.n

    The spike experiment

    The "scomp_fm.n" model file contains several experiments that are run on the same spike generator. The "spike" experiment is run by default. It generates 1 spike, showing the currents during the spike:

    neurc scomp_fm.n | vid
    or
    neurc --expt spike scomp_fm.n | vid
    
    This experiment is:
    if (expt=="spike") {                    /* look at spike shape */
    
      stim node soma cclamp 250e-12 start 0 dur .0005;
      stim node soma cclamp 150e-12  start .0005 dur .0015;
    
      endexp = .006;
    
      plgain = 5e-9;
    
      naoff = 0;
      nagain = plgain;
      caoff = 0;
      cagain = plgain;
      kdroff = 0;
      kdrgain = plgain;
      kaoff = 0;
      kagain = plgain;
      skoff = 0;
      skgain = plgain;
      bkoff = 0;
      bkgain = plgain;
    
      plot G(I) naxx  pen 1  max (naoff+1) * nagain   min (naoff-1) * nagain;    // blue trace
      plot G(I) kdr   pen 2  max (kdroff+1)* kdrgain  min (kdroff-1)* kdrgain;   // green trace
      plot G(I) ka    pen 14 max (kaoff+1) * kagain   min (kaoff-1) * kagain;    // yellow trace
      plot G(I) ca1   pen 3  max (caoff+1) * cagain   min (caoff-1) * cagain;    // cyan trace
    
    if (!notinit(skca1))
        plot G(I) skca1 pen 5 max (skoff+1)*skgain   min (skoff-1) * skgain;     // magenta trace
    if (!notinit(skca2))
        plot G(I) skca2 pen 8 max (skoff+1)*skgain   min (skoff-1) * skgain;     // gray trace (if uncommented)
    if (!notinit(bkca))
      plot G(I) bkca    pen 7 max (bkoff+1)*bkgain   min (bkoff-1) * bkgain;     // white trace (if uncommented)
    
      plot V[soma]    pen 4  max .16 min -.08;                                   // red trace
      plot S totcur pen 12 max plgain min -plgain;  /* total current from onplot(), dark red trace  */ 
    }
    

    Several of the plots show the separate currents from each channel type. There is also a plot of the compartment voltate showing the spike, and there is also a plot of the total current from all the channels (the dark red trace). The skca5 and bkca channels are disabled by default, but you can enable them by uncommenting them.

      plot G(I) naxx  pen 1  max (naoff+1) * nagain   min (naoff-1) * nagain;    // blue trace
    
    The "naoff" and "nagain" variables in the plot statement are suppoed to allow the user to offset the traces and change the relative scaling, to make the overlapping traces easier to see. For example, the Na current is much larger initially than the other currents, so sometimes it is helpful to increase nagain which then gives a smaller scaling to the Na current, allowing the complete waveshape to be easily seen. For most purposes, though, one wants to see all the currents at the same scale, so all the gains are set to the same value: +/- 5e-9.

    The lesson in the model is to try to understand the role of each current in generating the fast-rising and fast-falling spike. The model starts out with a 250 pA current to activate the Na current, ant at 0.5 msec the current is reduced to 150 pA. From there, the Na current slowly builds up until about 2 msec, after which it rapidly becomes regenerative. This meas that when the Na current activates, it depolarizes the cell which causes it to be more activated. This effect is known as "positive feedback" in control systems theory, and is responsible for the very fast rise of the spike's leading edge.

    After the spike has depolarized to its peak, the Na current inactivates. This allows the spike to repolarize (hyperpolarize) but it would not do this very fast unless K currents are present. So at the same time, the Kdr ("delayed rectifier") current gets activated by the depolarization of the spike, and starts to hyperpolarize the cell, once the Na current has decreased by inactivation. There is another effect on the Na current -- it decreases just when the spike is most depolarized because the driving force (Vna - Vm) for the Na current is less.

    Once the Kdr current turns on, the spike repolarizes in a msec or so -- and it actually hyperpolarizes below the original voltage. This hyperpolarization deactivates the Kdr channel -- it does not inactivate automatically as the Na channel does, but deactivates when the depolarization that originally activated it is not present. In addition, the KA channel is activated by the leading edge of the spike. It, however, inactivates like the Na channel soon after it is activated. Together the Kdr and KA currents terminate the spike. During their activation, and immediately following it, there is said to be a "refractory period" during which it is not possible to generate another spike, even with a large injected current. One function of the refractory period is to allow all the voltage-gated channels to recover from their activation during the spike, to allow another spike to occur with the same waveshape as the first one.

    The inter-spike interval

    In some spike generators that have only Na, Kdr, and KA channels, the KA current plays a special role in which it activates a little, soon after the refractory period, at a relatively hyperpolarized voltage. This slows the ramping up depolarization to the next spike, which may be important for some spike generators that require large currents to start spiking. By slowing the ramping depolarization, this allows the spike generator to initiate spikes farther apart in time (longer inter-spike interval), i.e. to make a slower spike train. But then KA will inactivate, depending on how much depolarization it sees, also at a relatively hyperpolarized voltage, which then allows faster ramping depolarization, and a spike to occur more immediately when the Na channel gets activated.

    In this scomp_fm.n spike generator, the role of KA appears to be more important in terminating the spike than in regulating the inter-spike interval. The reason is that other channels take the role of regulating the inter-spike interval. The Ca channel is a little slower than the Na channel, and but it turns on during each spike and allows some calcium to flow into the cell. This calcium is then a signal for how many and how fast the spikes occur. There is also a calcium pump that removes calcium from the cell and pumps it to extracellular space. The SK and BK Ca channels sense the calcium and hyperpolarize the cell when the calcium level goes up. The BK channel (commented out by default) is sensitive to both Vm (membrane voltage) and [Ca]i (the internal calcium level), so it is activated only during a spike when [Ca]i is high. It is very fast, which allows it to shorten the spike duration when the [Ca]i is too high, which will then reduce the amount of calcium that flows into the cell. BK channels are often involved in generating a strong after-hyperpolarization, but like Kdr they are also deactivated by hyperpolarization so they soon turn off after a few msec.

    The SK channels are activated only by internal calcium [Ca]i, not by voltage, so they act only on the accumlulated calcium in the cell. This makes their role very straightforward -- they hyperpolarize the cell after a burst of spikes when the [Ca]i gets high. This allows the cell to have spike frequency adaptation. While there are several other mechanisms that can accomplish spike frequency adaptation, such as slow inactivation of the Na channels (sometimes by delayed recovery from inactivation), the SK channels are a well-known regulator of spike frequency. The SKCa1 channel is much slower than Kdr, and has a time constant of ~50 msec, and the SKCa2 channel is even slower with a time constant of ~400-500 msec. This gives them full control over very long inter-spike intervals.

    To run the "isi" experment, enter:

    neurc --expt isi scomp_fm.n | vid
    
    else if (expt=="isi") {         /* look at inter-spike interval */
    
      waittime = 0.0004;
      stimdur  = 1;
      afterstim = 0.10;
      stimampl = 45e-12;
    
      stim node soma cclamp 200e-12 start 0 dur waittime;
      stim node soma cclamp stimampl start waittime dur stimdur;
    
      endexp = .025;
      //endexp = .5;
    
      plgain = 2e-10;
    
      naoff = 0;
      nagain = plgain;
      caoff = 0;
      cagain = plgain;
      kdroff = 0;
      kdrgain = plgain;
      kaoff = 0;
      kagain = plgain;
      kacgain = plgain*10;
      skoff = 0;
      skgain = plgain;
      bkoff = 0;
      bkgain = plgain;
      moff = 1;
      mgain = plgain;
    
      plot G(I) naxx  pen 1  max (naoff+1) * nagain   min (naoff-1) * nagain;
      plot G(I) kdr   pen 2  max (kdroff+1)* kdrgain  min (kdroff-1)* kdrgain;
      plot G(I) ka    pen 14 max (kaoff+1) * kagain   min (kaoff-1) * kagain;
    //plot G(0) ka    pen 2 max (kaoff+1) * kacgain   min (kaoff-1) * kacgain;
      plot G(I) ca1   pen 3  max (caoff+1) * cagain   min (caoff-1) * cagain;
    
    if (!notinit(skca1))
        plot G(I) skca1 pen 5 max (skoff+1)*skgain   min (skoff-1) * skgain;
    if (!notinit(skca2))
        plot G(I) skca2 pen 8 max (skoff+1)*skgain   min (skoff-1) * skgain;
    if (!notinit(bkca))
      plot G(I) bkca    pen 7 max (bkoff+1)*bkgain   min (bkoff-1) * bkgain;
    
      plot V [soma]    pen 4  max .16 min -.08;
      plot Im[soma]    pen 12  max (moff+1)*mgain min (moff-1) * mgain;
      plot Ca(1)[soma] pen 6 max 20e-6 min 0;      /* [Ca]i */
    
    }
    

    Role of currents in inter-spike interval

    In this "isi" (inter-spike interval) experiment, the emphasis is on what happens between the spikes. During the first spike, calcium in the cell [Ca]i rises, and then slowly drops due to being pumped out of the cell. But this calcium activates the SKCa1 channel, which can then regulate the ramping up depoarization to the next spike.

    Note that the largest current in the beginning of the inter-spike interval is the SKCa1 current. This means that it is the current that controls the overall spike rate (given a specific injected stimulus current), not the Na or KDR, KA channels. Typically in the axon of a neuron, there are only Na and Kdr channels, which take on the role of transmitting spikes quickly along the length of the axon. In the axon the Na and Kdr channels have overlapping voltage activations, i.e. both currents are typically activated at rest. This prevents any small currents from affecting when the next spike will occur -- a spike can only occur if a spike traveling along the axon depolarizes the Na channels in the presence of the continuous Kdr channel activation.

    However, in a spike generator, it is important that the Na and Kdr, KA currents are small, essentially off, at the resting membrane potential, i.e. that other currents control the resting membrane potential so they can regulate the inter-spike interval and thus the overall spike rate. The BK and SK channel currents can control spike rate because they are larger than the Na, KDR, KA currents. But the BK and SK currents are small. They have to be, to allow the synaptic input to the cell to be the main control for the spike rate. This is the main function of a spike generator -- to generate spikes according to very small depolarizations (EPSPs) from synaptic currents.

    If you uncomment the "skca2" and "bkca" channels, they add to the hyperpolarizing influence in the cell and increase the inter-spike interval. This makes the plot of all the currents a little more complex, and to see the next few spikes, you need to uncomment the second "endexp = 0.5" statement. This shows what happens in a typicalburst of spikes. The first spike lets in some calcium, but not enough to activate the BK and SK channels much. But after the second spike, the BK channels turn on at a more hyperpolarized voltage, and the SKCa1 chanel gets activated, which slow down the ramping up depolarization after the intial after-hyperpolarization. The spike rate eventually slows to a plateau frequency that is set by the BK and SK currents in response to the spikes generated by the injected current.

    If you then replace the comments around the bk and skca2 channel definintions, and run "expt f/i", you can see the spike rate at different currents. The little circles show the instantaneous spike rate from one spike to the next. If you then uncomment bk and skca2 and change their density with the "skca1_mult", "skca2_mult", and "bkca_mult" parameters, you can create different amounts of spike frequency adaptation.

    neurc --expt f/i --bkca_mult 0 --skca2_mult 1 scomp_fm.n | vid
    
    To see the influence of the random gating of the channels, run the model with "channoise 1" on the command line:
    neurc --expt isi --bkca_mult 0 --skca2_mult 1 --channoise 1 scomp_fm.n | vid
    
    This gives some intuition about what the noise that typical neurons have to deal with. The greater density of channels, the less the noise problemm becomes.

    Other scomp models

    The other scomp models are built upon scomp_fm.n, but are designed to describe a mammalian cell, which lives at warmer temperatures than a salamander. Their "tempcel" (temperature of the cell in Celsius) is set to 35 degrees which greatly increases the kinetic rate of the activation and inactivation of the channel currents. The Na current rises much faster, and the depolarization of the spike rises much faster, typically in 0.1 - 0.2 msec. In this case, the capacitance of the cell is much more of a factor, i.e. since the rise time is much faster, for the same capacitance (set by the cell dimensions) the voltage-gated currents must be much larger, typically 3-fold greater (140e-3 vs.50e-3 S/cm2). In scomp[a-e].n the channel densities are set using "ndensity" which is the number of channels per um2. This is calculated from the standard for channel densities, S/cm2, by dividing by the unitary channel current amplitude, given by qcond(unitarycurr), where qcond() is adjusts the unitary current from cm2 to um2 also by the model temperature. In the "scomp_ad?.n" models, the ndensity values are given as n/um2 which is equivalent to the other way of specifying, either using qcond() or using "density". This is supposed to help the user get more intuition about the actual density of channels in a typical cell.

    Scomp summary

    Try modifying any of these models, by changing the channel densities and looking at the different experiments. See if you can understand how the currents develop at the beginning of the spike, how the spike ends, and how the spike interval is controlled. You can do this by adding parameter values on the command line.

    Overall, you can see how complex a single compartment model can be. All the channels are nonlinear systems in themselves, and when 2-5 channel types are combined, the outcome can be difficult to predict. In a real neuron, of course, there are voltage-gated, ligand-gated, and calcium-gated channels throughout the cell, so a multi-compartment model is necessary to capture realistic behavior. Interactions, i.e. lateral current flows, between compartments make a multi-compartment much more difficult to model realistically.

    Multi compartment spiking neuron

    The second paper in the Fohlmeister & Miller (1997b) series shows how adding a dendritic arbor to a neuron greatly affects how the spike generator works. Before this classic paper, it was thought that ganglion cell dendrites were essentially passive collectors of the synaptic currents that define the spike rate. However, Fohlmeister discovered that adding a dendritic arbor to a single compartment spike generator that had been already calibrated caused the spike generator to fail.

    Each spike that was initiated in the soma back-propagated out into the dendritic arbor, and if the dendritic arbor was passive, i.e. it contained no voltage-gated channels, the spike decayed as it propagated out to the dendritic tips. This charged up the capacitance of the dendritic arbor, which is typically much greater than the capacitance of the soma, and the charge then traveled back to the soma, where it initiated the next spike too soon. It was a classic case of "impedance mismatch". A passive dendritic arbor makes the cell spike too fast! By adding Na and Kdr channels at a slightly lower density than at the soma, spikes can then back-propagate actively, and the Kdr channels would discharge the capacitance as the spike propagate out to the dendritic tips, so it cannot propagate back to the soma. This arrangement allows the spike generator to fire spikes slowly, as did the original single comparment model.

    Overall, Fohlmeister found that different dendritic morphologies and their channel densities affect the spike rate generated by a given input current. It was later verified by Velte & Masland (1999) that indeed ganglion cells do actively back-propagate spikes in their dendrites. Although large currents injected into dendrites can initiate spikes due to the active voltage-gated channels present, generally the Na channel density is low enough that typical light-driven synapic inputs don't initiate spikes in most ganglion cell dendrites. The synaptic currents normally propagate to the soma and initial segment of the axon to initiate a spike there, and then actively back-propagate into the dendrites. The one exception is the direction-selective ganglion cell (DSGC) that was found to initiate spikes in its dendrites evoked by light input by Oesch et al (2005).

    To see a multi-compartment model of a spike generator, you can go to nc/models/spikegen and run "nc -v spikegen.n | vid". This model was developed for the van Rossum, O'Brien & Smith (2003) paper which investigated the influence of noise properties in multi-compartment mammalian model. The "spikegen" model also contains several experiments that can be specified by setting the "--expt" parameter on the command line.

    When noise properties are added, it is even more difficult to achieve realism -- but noise properties are essential for understanding the problems that neurons face.




    Section 4)

    Design philosophy

    Design philosophy: "keep it simple!"

    Often when setting up a neural circuit model it is helpful to remember that a simple model is easier to understand and possibly more believable than a more complex one. This is generally true because the utility of a computational model is often that it quantitatively assists (stretches) our imagination. If a model is so complex that we don't understand it, then it does not help us. If a model uses many parameters for which numerical values are not known, then it may be a focus for disbelief.

    "Do it in stages"

    It is best to first make a simple model to verify the operation of the model's components. Then enlarge the model after it gives a threshold of confidence. For instance, to model a large circuit with many neurons involved in negative synaptic feedback, it is best to start with one cell and get its anatomy and membrane properties worked out. Next set up 2 cells and their synaptic feedback and verify that the feedback is working properly using appropriate stimulus and recording sites.

    After the anatomy and basic synaptic model is all working properly, you should expand the model to include more cells and complexity. The advantage of this technique is that any problems that occur can usually be localized to the addition of a few neural circuit elements, because everything was developed and verified in successive stages.

    If multiple experiments on a specific model are put into one model file, this is generally OK if the file isn't too long and the model is the same for each experiment. But if different experiments need variations on the construction of the model, it is best to make a new model file. However, this can cause another problem, of having too many files, all closely related -- this can cause confusion and slow development of new models. It is better to share common parts of multiple related files -- to reduce the number of files, for clarity and simplicity.

    Retsim

    The "spikegen.n" model described above was the precursor for the "retsim" retinal simulator. It called procedures from several other files in the "nc" interpreter language -- these gradually developed into an Applications Programmer Interface (API) that could construct different neurons and connect them synaptically. The original "retsim1.n" grew out of this project. The rationale behind "retsim" is to construct the model with a common set of algorithms using standard neuron morphologies, biophysical properties, and synaptic connection types. The experiment file requests the model to be constructed in this standard manner, and then need only describe the experimental protocol that is run on the model, defining stimuli, recording points and plots.

    Retsim was next revised using C/C++ as the script language (i.e. not directly using the nc interpreter) using the original nc API which is written in C/C++, but revising retsim1.n into retsim.cc. For retsim.cc, each experiment file is a separate module, dynamically linked at runtime, and generally contains one experiment on one specific type of model. This keeps each model simple and easy to test. Once the model is set up correctly in the experiment file, then mulitiple experiments can be run on it, either in the same experiment, or in different but related experiment files.




                   Copyright (C) 2024
    
                     Robert G. Smith
                 Department of Neuroscience
              Univ. of PA Medical School
                Phila., PA   19104-6058
                          USA
    
                  All Rights Reserved