# 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.

### First Step

Read these two key papers describing the model:
Additional papers related to the model you may find useful:

### Second Step

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.

### Third Step

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:
```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));
```
In this example, `find_anorm` is a function that computes anorm, as described in the revised step 2 above. You will find it advantageous to declare `c` and `currents` to be global variables.

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 `timeline` of 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 an array `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.)

Extend your `updatePyramid` function to include the theta current. Rerun the simulation and verify that theta ipsp's are present. In the above plot, the black asterisks mark theta cell spikes.

### Fourth Step

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.

Modify the `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 threshold.

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 `inputSpikes(p,i)` encodes 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 step. ### Fifth Step

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 fire.

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.

Write an `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.

Modify your `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. ### Sixth Step

Write a function `spikestats` that is called after the simulation has run, and prints out the following statistics separately for each pyramidal cell:
• 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
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.

### 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.

### Last-Minute Notes

• 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 `find_anorm` 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.

Dave Touretzky
Last modified: Mon Dec 14 04:32:32 EST 2009