20-CS-110-001 Introduction to Computer Science Fall 2010

Lab Number 6

Neuron Firing Simulation
Probability, Medical Research

Due: November 30, 2010
Submit one matlab program as two separate files via blackboard. See submit page for instructions.

    The objective is to practice the use of MATLAB random variables to better understand random variables, distributions, and their use in engineering and medical research problems.


According to the Wikipedia article entitled "Neuron":
    Neurons are responsive cells in the nervous system that process and transmit information by electrochemical signaling. They are the core components of the brain, the vertebrate spinal cord, the invertebrate ventral nerve cord, and the peripheral nerves. Motor neurons receive signals from the brain and spinal cord and cause muscle contractions and affect glands. Neurons respond to stimuli, and communicate the presence of stimuli to the central nervous system, which processes that information and sends responses to other parts of the body for action.

A neuron's response to a stimulus is analogous to that of computer circuitry. Imagine you are watching a little neuron in the brain (there are no big ones, I think - anyway they are ugly so no picture of them is shown). In response to an incoming electric current, the neuron's voltage increases. When the neuron reaches a specific "threshold" voltage (in this model), the neuron "fires" what is called an Action Potential (AP). Firing an AP causes changes in the conductances of ion channels on the neuron's cell membrane with the result that the neuron's voltage drops sharply back to 0. After a neuron fires an AP (this event is also known as a spike) there is a refractory period during which it is impossible to fire another AP due to the need for specific ion channels to "recover."

A simple electrical model of the spike is shown in Figure 1.

Figure 1. Electrical model of neuron firing. The resistance is R and capacitance is C.

In that model Vcc is a constant voltage and Vneuron is the neuron's voltage which is initially 0. When the switch is closed at time t0, Vneuron rises as shown on the right side of the figure until the discharge voltage of the neon lamp is reached at time t1. At that point, the capacitor discharges through the lamp, its voltage rapidly is reduced to 0, and the switch is then opened to keep it there. The rise of Vneuron may be determined easily from the resistor and capacitor. Let C be the capacitance of the capacitor and R the resistance of the resistor. Since the current through the capacitor is given by iC=C(dVneuron/dt), the current through the resistor is given by iR=(Vcc-Vneuron)/R, and kirchoff's law says iR=iC, we can write:

   Vcc-Vneuron = RC(dVneuron/dt) → dVneuron/dt + Vneuron/RC = Vcc/RC.

This is a first order differential equation with the following solution for the appropriate initial conditions:

   Vneuron = (1-e-t/RC)Vcc.
Therefore, the time profile of a spike has the neuron's voltage starting at 0, rising exponentially toward Vcc, and ending with a sharp drop to 0 beginning at the point in time where the neuron's voltage reaches the threshold voltage. The speed at which the voltage rises is determined by the product RC. Hence, RC is called the time constant and denoted by τ.

It is not known precisely what triggers a spike but it seems that normal nervous functions correlate with certain spike frequencies and patterns or codes whereas abnormal nervous functions correlate with different patterns as shown in Figure 2. Study of these patterns has resulted in significant progress in treating Parkinson's disease, for example, by introducing electrical stimuli that control neuron firing. See


for a video showing a Parkinson's disease sufferer trying to get through a doorway and see


for a video showing the result of neuron spike stimulation.

Finally see


for the result of deep brain simulation.

Figure 2. Schematic illustration of differences in spike patterns, normal (top) vs. epileptic (bottom).

The Geometric Distribution
Consider an ordinary, old outdoor gas grill that you might use once a year on July 4th to cook hamburgers. To light the grill you need to press a "starter" button which generates a spark that ignites the gas. Assume the probability that the spark ignites the gas on a single press is 1/2 (that is, half the time you have to try again). What is the probability that you have to press the starter n times before the grill lights? In other words, what is the probability that the starter fails n-1 times, then hits?

Let's assume that ignition is an event that is statistically independent of previous presses of the starter. If n is 1, the probability is obviously 1/2. If n is 2, the probability of a failure followed by ignition is (1/2)(1/2) = 1/4. If n is 3 the probability of two failures followed by ignition is 1/8. More generally, it's 1/2n.

Now, say the probability that the starter ignites the gas on a press is p. Then the probability that n-1 consecutive failures are followed by ignition is p(1-p)n-1. Let's lay out these probabilities consecutively on a line of length 1 (we can do this because the sum of all the probabilities must be 1 - that is, Σ1≤n≤∞p(1-p)n-1 = 1).

Starting from the left we mark off distance p which is also 1-(1-p) as shown by the first tick to the right of 0 in Figure 3.

Figure 3. Geometric distribution intervals on the line from 0 to 1.

Then we mark off an additional distance of p(1-p) which is placed a distance 1-(1-p)+p(1-p) = 1-(1-p)2 from 0 as shown by the second tick to the right of 0 in the figure. Next we mark off an additional distance of p(1-p)2 which is 1-(1-p)2+p(1-p)2 = 1-(1-p)3 from 0. Continue this indefinitely.

Now, choose a random number from 0 to 1. Call it X. Find n so that 1-(1-p)n-1 < X ≤ 1-(1-p)n. The probability that X takes such a value is p(1-p)n-1 which is the probability that n-1 failures occur before an ignition. Therefore, we can use a uniform random number generator to produce random values as though they originated from a geometric distribution. However, to more closely model reality, we are going to use the uniform random number generator to produce random values that are real numbers instead of integers according to X = 1-(1-p)n. If you negate both sides, add 1 to both sides, and then take the log of both sides you will be able to get an expression for n in terms of X. This is

   n = log(1-X)/log(1-p).
In the above equation n is a random value and p is related to the mean of n, denoted μn. This relationship is derived as follows:
  μn = 0≤n≤∞ np(1-p)ndn ≈ 0≤n≤∞ npe-npdn = -ne-pn|n=0,∞ + 0≤n≤∞ e-pndn
     = -(1/p)e-pn|n=0,∞ = 1/p.  
Since X translates to rand(1) in Matlab (it's a random number between 0 and 1) this results in an expression that is written in Matlab as follows:
   n = log(1-rand(1))/log(1-p);
which has the same effect as
   n = log(rand(1))/log(1-p);
since X and 1-X have the same probability of being generated by rand(1). In this lab the mean of the distribution of n will be given as Mean instead of p. From the above, μn=1/p, so the Matlab expression becomes
   n = log(rand(1))/log(1-(1/Mean));

The Matlab statement above will be used below to generate random effects on neuron spikes, particularly in modeling neuron firings as though they are hard-to-start gas grills.

The Problem:
    Your assignment is to write a Matlab function that simulates neuron firing. The long-term objective of such simulations is to produce a mathematical model for accurate prediction of neuron firings in normal and abnormal circumstances. The function will do the following:
  1. Generate a plot of a neuron's spikes over time
  2. Record the times at which the neuron fires
  3. Generate a plot of the Inter-Spike-Intervals (ISI) as a histogram
  4. Calculate the coefficient of variation of ISI
where the coefficient of variation of a random variable X is defined as the ratio of the standard deviation of X to its mean and the Inter-Spike-Interval between two consecutive spikes is the time that has elapsed from the moment the earlier spike reaches its peak to the moment the later spike reaches its peak. The simulation will span a specified time period called tot_time, measured in seconds. Total time will be partitioned into a number of points in time, called epochs, which are spaced timestep seconds apart. All changes to voltage values will be made at epochs only. Epochs are numbered in increasing order of occurrence in time from the start of the simulation, beginning from 1. Spike voltage (Vneuron) will be calculated at epoch i using information obtained at epoch i-1 and at the most recent spike peak. We use Matlab notation Vneuron(i) to denote the value of Vneuron at epoch i, Vcc to denote Vcc, and Vthr to denote threshold voltage. Values for Vneuron(i) are computed as follows:
  1. If the time elapsed since the last peak is no greater than 5 milliseconds then Vneuron(i)=0;
  2. otherwise, let Dt denote the time since the start of the current spike (a spike starts at some epoch that occurs at least 5 milliseconds after the previous peak and Vneuron is 0 at the start epoch of a spike) and compute Vtmp=(1-exp(-(Dt/RC)))*Vcc. If Vtmp is greater than Vthr then a peak has been reached, it must be recorded, and Vneuron(i)=0, otherwise Vneuron(i)=Vtmp.

Your function, called neuron, will take four input arguments, all scalars. The prototype for the function should be the following:

          function neuron(dist, tot_time, timestep, Mean)
The argument dist will take one of the values in the set {0,1,2,3}. If the value of dist is 0, the value of Vcc will be fixed at 5 and the threshold voltage Vthr will be fixed at 1. If the value of dist is 1, the threshold voltage Vthr will be 1 plus a number determined randomly (recalculated for each spike) from a uniform distribution with mean Mean no greater than 1.5, and Vcc will be fixed at 5. That will simulate a leaky neuron. If the value of dist is 2, Vcc will be 5 plus a number determined randomly for each spike from an exponential distribution of mean Mean, and Vthr will be fixed at 1. If dist is 3, then Vthr will be fixed at 1, Vcc will be fixed at 5, and the start time of a spike will be 5 plus a number determined randomly from the continuous version of the geometric distribution of mean Mean.

The constants used by the function are:

  1. Refractory period: 5 milliseconds
  2. The product τ=RC: 70 milliseconds

Important Notes: Be careful not to choose variable names that are thr same as Matlab functions - for example, I used Mean instead of mean above as the name for the variable that holds the mean of the distributions considered in this lab.

The units of tot_time and timestep are always milliseconds. The units of Mean are milliseconds when dist takes value 3, and volts otherwise. Therefore,

results in a simulation of 10 seconds, using timesteps of 1 millisecond, with a mean firing delay of 5 milliseconds after the refractory period. The histogram of such a simulation is shown in Figure 4.

Figure 4. ISI histogram from neuron(3,10000,1,5).

A plot of

is shown in Figure 5.

Figure 5. ISI plot from neuron(3,200,1,5).

Analysis and Coding:
    You probably want to define variables that save the epoch on which the most recent peak occurred and the epoch at which the current, increasing spike began. Call the first one last_peak and the second one last_start. Then you can compute Dt as follows:
     Dt = timestep*(i-last_start);
where i is the current epoch, and you can make the test to determine whether the refractory period has finished as follows:
     if (i-last_peak)*timestep <= Refrac
        Vneuron(i) = 0;
The number of epochs through which the simulation should run is determined by tot_time/timestep. You will have to define a vector, say called peaks, that will record the epochs at which peaks occur. A variable, say k, initialized to 1, can be used to index into peak to record values. At every epoch during which time a spike is rising you may want to compute
     Vtmp = (1-exp(-Dt/RC))*Vcc;
and then determine whether the peak has been reached using
     if Vthr >= Vtmp
        last_peak = i;
        peaks(k) = i*timestep;
        k = k+1;
Create a histogram vector ISI and fill it like this:
     for i=1:(length(peaks)-1)
        ISI(i) = peaks(i+1)-peaks(i);
Plot a histogram like this:
Use the Matlab functions std and mean to compute the standard deviation and mean of the inter spike intervals so you can compute the coefficient of variation of ISI.