20CS110001  Introduction to Computer Science  Fall 2010 


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.
Objective:
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. 
Background:
Neurons
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 V_{cc} is a constant voltage and V_{neuron} is the neuron's voltage which is initially 0. When the switch is closed at time t_{0}, V_{neuron} rises as shown on the right side of the figure until the discharge voltage of the neon lamp is reached at time t_{1}. 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 V_{neuron} 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 i_{C}=C(dV_{neuron}/dt), the current through the resistor is given by i_{R}=(V_{cc}V_{neuron})/R, and kirchoff's law says i_{R}=i_{C}, we can write: V_{cc}V_{neuron} = RC(dV_{neuron}/dt) → dV_{neuron}/dt + V_{neuron}/RC = V_{cc}/RC. This is a first order differential equation with the following solution for the appropriate initial conditions: V_{neuron} = (1e^{t/RC})V_{cc}.Therefore, the time profile of a spike has the neuron's voltage starting at 0, rising exponentially toward V_{cc}, 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 http://gauss.ececs.uc.edu/Courses/c110/labs/Video/Problemdoorway.wmv for a video showing a Parkinson's disease sufferer trying to get through a doorway and see http://gauss.ececs.uc.edu/Courses/c110/labs/Video/Solutiondoorway.wmv for a video showing the result of neuron spike stimulation. Finally see http://gauss.ececs.uc.edu/Courses/c112/PS/DBS.wmv for the result of deep brain simulation.
Figure 2. Schematic illustration of differences in spike patterns, normal (top) vs. epileptic (bottom).
The Geometric Distribution
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/2^{n}. Now, say the probability that the starter ignites the gas on a press is p. Then the probability that n1 consecutive failures are followed by ignition is p(1p)^{n1}. 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(1p)^{n1} = 1). Starting from the left we mark off distance p which is also 1(1p) 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(1p) which is placed a distance 1(1p)+p(1p) = 1(1p)^{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(1p)^{2} which is 1(1p)^{2}+p(1p)^{2} = 1(1p)^{3} from 0. Continue this indefinitely. Now, choose a random number from 0 to 1. Call it X. Find n so that 1(1p)^{n1} < X ≤ 1(1p)^{n}. The probability that X takes such a value is p(1p)^{n1} which is the probability that n1 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(1p)^{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(1X)/log(1p).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(1p)^{n}dn ≈ ∫_{0≤n≤∞} npe^{np}dn = ne^{pn}_{n=0,∞} + ∫_{0≤n≤∞} e^{pn}dn = (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(1rand(1))/log(1p);which has the same effect as n = log(rand(1))/log(1p);since X and 1X 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 hardtostart gas grills. 
The Problem:
Your assignment is to write a Matlab function that simulates neuron
firing. The longterm 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:
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:
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, neuron(3,10000,1,5);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 neuron(3,200,1,5);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*(ilast_start);where i is the current epoch, and you can make the test to determine whether the refractory period has finished as follows: if (ilast_peak)*timestep <= Refrac Vneuron(i) = 0; ... else ... endThe 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 = (1exp(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; ... endCreate a histogram vector ISI and fill it like this: for i=1:(length(peaks)1) ISI(i) = peaks(i+1)peaks(i); endPlot a histogram like this: figure(2) hist(ISI);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. 