15-883 Modeling Project
Computational Models of Neural Systems, Fall 2009
In this project you will replicate the short-term persistent spiking
buffer model of entorhinal cortex described by Randal A. Koene and
Michael E. Hasselmo in a series of recent papers.
Read these two key papers describing the model:
Additional papers related to the model you may find useful:
Write code to plot the conductance function gi(t) for any
of the pyramidal cell currents listed in Table 1 of the Brain Research
paper. (Plot the conductance values between 0 and 500 ms and assume a
spike occurs at 25 ms.) Also plot the calculated value of
tmax as a vertical line. Note: you will discover a problem
with the equations when you do this. Come up with a solution.
NOTE: Change to calculation of
anorm: The equation given in the paper for tmax is
incorrect. The correct equation is
tau_fall*tau_rise*log(tau_rise/tau_fall) / (tau_rise - tau_fall). To
solve for anorm directly, calculate the value of the bi-exponential
term in gi(t) for all t from 0 up to the value of tau_fall, in increments of 0.01 msec. Find the maximum of
these values. Set anorm equal to the reciprocal of the maximum. Once
you've done that, you will have guaranteed that gi(t) will peak at Gi.
The position of the maximum in the list of values will give you tmax,
which you only need for your conductance plots; it's not needed for
the simulation itself.
Create a data structure to describe the various currents that make up
a pyramidal cell, most of which are listed in Table 1. For example,
you might use a struct array whose first entry, for the AHP current,
looks like this:
In this example,
c.AHP = 1;
currents(c.AHP).tau_rise = 1e-4;
currents(c.AHP).tau_fall = 30;
currents(c.AHP).G = 23;
currents(c.AHP).Erev = -90;
currents(c.AHP).anorm = find_anorm(currents(c.AHP));
find_anorm is a function that computes
anorm, as described in the revised step 2 above. You will find it
advantageous to declare
currents to be
NOTE: Changes to parameter
values: For the ADP current, use tau_rise = 135 - 1e-4, and
tau_fall = 135. For the theta current, use tau_fall = 8. For the
gamma current, use G=150.
Create a struct array to describe the set of pyramidal cells in your
model. Two values you must keep for each cell are: its current
membrane voltage V, and the time it last emitted a spike. You may
wish to keep additional values around, such as a state indicator
(normal, spiking, or refractory), the cell's current total
conductance, and the most recent voltage adjustment, ΔV.
At the start of your simulation you will specify the number of
pyramidal cells P, and the length of the simulation run in
milliseconds Assume that time starts at t=0 and advances in increments
of Δt = 1 msec. From this you can calculate the number of
steps S in the simulation. Create an array
length S containing the time (in msec) at each step, and a history array
Vhist of size P × S that will hold the
membrane voltage of each pyramidal cell at each time step. You will use
Vhist for plotting the results of your simulation.
Write a function
updatePyramid(p,i) where p is the
pyramidal cell number (from 1 to P) and i is the current time step
(from 1 to S). Your function should compute the new membrane voltage
V(t+Δt) based on the previous voltage V(t) and the current
conductances gi(t). Start with a very simple cell that
just has a leak current.
Write the first version of the main loop of your simulation. It
should progress through all the time steps, recalculating the membrane
voltage of each pyramidal cell at each step. Set the cell's initial
membrane voltage to something other than the resting potential of -60
mV. When the main loop finishes, plot the membrane voltage history
and verify that the cell settled to its resting potential.
Theta modulation from the medial septum is simulated by a series of
theta "spikes"; each spike triggers an inhibitory conductance. Create
thetaSpikes of size S containing, for each
entry, the time of the most recent theta spike. (At a theta frequency
of 8 Hz, the last theta spike time should stay the same for 125 msec,
then jump ahead by 125 msec., etc.)
updatePyramid function to include the theta
current. Rerun the simulation and verify that theta ipsp's are
In the above plot, the black asterisks mark theta cell spikes.
Now we need to make the neuron spike. Create a P×S array called
inputSpikes that contains the last spike time of an input
to each of the P pyramidal cells. These spikes will control an
"input" current whose parameters are not listed in Table 1, but you
can find them on p. 3 of the "Reverse and Forward Buffering" paper;
use 0 mV as the reversal potential.
updatePyramid function to include the AHP,
ADP, and input currents, and add logic to make the cell spike if the
membrane voltage goes above threshold. When the cell enters the
"spiking" state it should hold its membrane voltage at 0 mV for 1
millisecond; then it should enter a refractory state for 1 millisecond
where it refuses to spike no matter what the membrane voltage is.
When the refractory period is over, the cell should return to its
normal state and can spike again if the membrane voltage goes above
Write a function
setInput(p,ts) that sets up an input
spike for pyramidal cell number p at time ts by writing appropriate
values into the
inputSpikes array. You call this
function at the beginning of your simulation for each input spike you
want to apply. Note that since
the time of the last spike up to and including step i of the
simulation, when you add a spike to the array you must modify all the
succeeding entries on that row.
Set up an input spike for the first pyramidal cell at time t=100, and
run your simulation with two cells. Verify that the first cell fires
repeatedly while the second one remains inactive, but both cells show
theta ipsps. Note: you must set the firing
rate threshold to -57 mV rather than the -50 mV value given in the
paper or the -54 mV value given previously.
Now set up an input to the second pyramidal cell at time t=225. Run
the simulation again. Notice that the relative spike times of the two
pyramidal cells are not stable. We will fix this in the next
The inhibitory gamma interneuron prevents pyramidal cells encoding
different items from firing at the same time. As soon as one
pyramidal cell fires, the interneuron supplies brief inhibition to all
the pyramidal cells, which will delay any cell that was just about to
To implement the gamma interneuron you will need to create a variable
gammaNeuron containing the same kind of structure as a
pyramidal cell. (There is only one of these interneurons in the
simulation so you don't need an array of them.) This neuron should
receive excitatory input from the pyramidal cells, which you can model
as a single input current using the last spike time of the most
recently-fired pyramidal cell.
updateGamma function to update the state of the
gamma interneuron, using the appropriate parameters for that neuron instead
of the ones for the pyramidal cell.
updatePyramid function to include inhibitory
input from the gamma interneuron. Modify your display code to display
the gamma interneuron's firing.
In the above plot, the black asterisks mark theta cell spikes, and the
membrane voltage of the gamma cell is shown as a black line. The
first input is at t=100 msec and the second at t=475 msec. The gamma
cell trace is scaled and translated downward to get it out of the way
of the pyramidal cell voltage traces.
Run your simulation with four pyramidal cells for 1510 time steps.
Input should come at t=100 for the first cell, t=225 for the second,
t=355 for the third, and t=605 for the fourth. Verify that the
relative spike times are stable.
Write a function
spikestats that is called after the simulation has run,
and prints out the following statistics separately for each pyramidal cell:
We're applying external inputs late in the theta cycle. If your
simulation is functioning properly you will see that with repeated
spiking the phase advances, but the phase of the earliest memory item
cannot advance too far due to theta inhibition, and the phases of the
remaining items cannot advance too far due to gamma inhibition
resulting from the firing of preceding items in the buffer.
- Time of each spike
- Inter-spike interval for each spike, i.e., the time between this spike and the previous one
- Phase of each spike relative to the theta cycle
What to Hand In
Email a zip file containing your code, plus a screen capture of the
plot your simulation produces, and the output of the spikestats
function. Your projects are due by 5:00 pm on December 14.
- See this suggestion about global variables.
- In order for gamma inhibition to work correctly you must have the
correct value of anorm; the instructions for
were modified to specify that you should search for the maximum
conductance using t values in increments of 0.01 msec. (If you use
increments of 1 msec you will get a poor result.) Also, for the gamma
current, the G parameter must be increased from 100 to 150.
- If you're still using the old Δt value of 0.1 msec, change it to
1 msec. Your simulation will run 10 times faster. Also, the
conductance parameters given here were tuned for Δt = 1 msec. To use a
finer-scale Δt value, which would reduce quantization error, the
conductances would have to be tweaked a bit. You're not being asked
to do that.
- The table of parameters in the paper shows the leak conductance
as 111, but the text says that gleak =
C/τleak, which might lead you to infer that the value
should be 1/9 = 0.111. They messed up on the units. The 111 value is
correct. In general, don't worry about getting the units right; just
use the numbers in Table 1.
Last modified: Mon Dec 14 04:32:32 EST 2009