Model Control: Range Variables

The task in the previous exercise was to characterize the effect of an independent variable (synaptic location) on a simple dependent variable (peak depolarization at the soma). Now we turn to a somewhat more complex problem: determining how the spatial distribution of biophysical properties affects the peak synaptic conductance needed to trigger a spike. Synaptic location is presumed to remain constant, although it is easy to imagine generalizing this question.

This exercise will show you how to

specify inhomogeneous range variables
change membrane properties under program control
use the APCount point process to detect spikes
use the Grapher to generate a family of simulations
write short (e.g. 3 line) procedures or functions that greatly enhance the power of the GUI

Physical System

The presence of voltage-gated channels distributed widely over the cell surface appears to be a general property of neurons. However, the density and even the biophysical properties of these channels may vary substantially with location. Although quantitative details generally remain uncertain, new experimental observations are beginning to provide an empirical basis for exploring how the spatial distribution of active currents may affect neuronal function.

Model

This exercise is based on the familiar ball and stick model. Imagine that HH channels are present on the soma and the proximal dendrite. Furthermore, channel density along the dendrite falls off linearly with distance from the soma, finally becoming An excitatory synapse is located at the distal end of the dendrite.

Simulation

Start NEURON with the init.hoc in the /course/range directory. This calls up the same init.hoc, ballstk.ses, and start.ses that we saw at the beginning of the "Simulation Families" exercise. Once again, be sure that the CellBuilder's Continuous Create button is checked. Otherwise the sections of the ball and stick model will not exist.

It's always a good idea to test things as you go along, so create an AlphaSynapse with tau = 1 ms, e = 0, and place it at the distal end of the dendrite. Try a few runs and manually adjust gmax until you trigger a spike. These tests are always reassuring when they work, but the dendrite has only one segment at this stage, so you shouldn't take the numerical value of gmax too seriously.

Be sure to include a "control" run with synaptic conductance = 0. Did soma.v drift away from the resting potential of -65mV? Why? (hint: do you need to change dendritic e_pas?) Strategy

You need to use the GUI and hoc code to build an interface that will automate the following tasks:

For selected locations along the dendrite
Set up channel distribution : create a linear gradient of HH channels in dend that falls to 0 at x
Find threshold gmax : the smallest synaptic conductance that triggers a somatic spike
Plot result : gmax vs. x

In the previous exercise, you used a for loop to generate a family of simulations, plotting the final results in a graph that you created with hoc code. For this exercise you will use the Grapher, which combines a for loop and a graph in a way that can be quite convenient for controlling a series of simulations.

You have already used the Grapher in a very direct way: to generate an arbitrary forcing function for a voltage clamp. The current exercise will show you how to exploit the Grapher in a way that is more subtle and far more powerful.

Preparing to set up channel distribution

The model needs some initial adjustments before you can use it to explore the effects of altered ion channel distributions.

The dendrite should have a reasonably fine spatial grid. Use the CellBuilder to change dendritic nseg to 25 (odd multiples of 5 guarantee that there will be nodes at 0.1, 0.3, 0.5, 0.7, and 0.9).
Before you can adjust the density of any ion channels, the conductance mechanism in question must first be inserted into the section. Any reference to a mechanism that has not been inserted will elicit an error message. For example, the first command below is OK but the second fails:
```
oc>print dend.g_pas(0.3)
0.0001 
oc>print dend.gnabar_hh(0.3)
hh mechanism not inserted in section dend
D:\NRN\BIN\NEURON.EXE:  near line 28
print dend.gnabar_hh(0.3)
```
Use the CellBuilder for this (in Biophysics / Specify Strategy make sure both pas and hh are checked for dend). Also change e_pas to -65 mV (why?). Finally, to avoid any confusion about the initial configuration of your model, it may be a good idea to reduce all of the hh conductances (including gl_hh) in dend to 0.

At this point you probably should save the CellBuilder by itself to a new session file (call it something other than ballstk.ses, maybe ballstktaper.ses).

Set up channel distribution

Next you will want to construct a procedure that establishes the tapering HH channel density in the dendrite, e.g.

proc hhtaper() {
	dend {
		gnabar_hh = 0
		gkbar_hh = 0
		gl_hh = 0
		if ($1 > 0) {
			if ($1 > 1) {$1 = 1}
			gnabar_hh(0:$1) = soma.gnabar_hh:0
			gkbar_hh(0:$1) = soma.gkbar_hh:0
			gl_hh(0:$1) = soma.gl_hh:0
		}
	}
}

This procedure could be placed in the ballstk.hoc file; an argument could be made for putting it in a separate file that is xopen'ed by init.hoc. An appropriate name for such a file might be taper.hoc.

Test your code with a few different arguments (e.g. 0, 0.5, and 1). First, verify that the conductances have the correct profile and change appropriately. Then see if increasing the dendritic extent of HH channels has the expected effect on threshold gmax.
Hint: use the GUI to make a space plot that shows gnabar_hh throughout the cell; after invoking hhtaper() with a new argument, update this graph by pressing RunControl's Init button. Alternatively, try the command
dend for (x) print x, " ", x*L, " ", gnabar_hh(x)

Find threshold

The basic idea is to repeat these actions

Run a simulation
If a spike occurred
reduce gmax
otherwise increase gmax

You have already seen that the hoc statement run() is equivalent to pressing the Init & Run button.

Whether or not an action potential occurs can be most easily determined by an APCount point process located at the soma. If APCount[0].n > 1 after a run, then voltage passed APCount[0].thresh (default -20 mV) at least once during the run.

Attach an APCount to the soma (hint: use a new Point Process Manager) and verify that APCount[0].n depends on whether a spike occurs.

The only remaining trick is to reduce or increase gmax by an amount that is governed by a clever rule that progressively homes in on the threshold.

A binary search is a convenient and reasonably efficient way to find a threshold. The kernel of the algorithm is:

while ( high - low > resolution) {
	run()
	if (APCount[0].n > 0) {
		high = independent_variable	// upper bound
	}else{
		low = independent_variable	// lower bound
	}
	independent_variable = (high + low)/2	//average
}

where the interval gets smaller by half on each iteration of the loop. After 10 iterations, the interval would be 1/1024 of its original size -- much better than testing each of the 1024 subintervals.

Actually, some elaboration of this algorithm is required to make a sensible choice of initial values of high and low. Also, the "independent_variable" has to be the one YOU choose. For this reason I am supplying a thresh.hoc file which contains a binary search implementation of

	final_val = threshold(&var_name)

You will need to edit your init.hoc file to load thresh.hoc by adding the line

	load_file("thresh.hoc")

You can test the function with

	threshold(&AlphaSynapse[0].gmax)

While testing, it is a good idea to display the soma voltage and the parameters of AlphaSynapse[0] and APCount[0].

Hints to keep run time short:
1. Set Tstop to 15 ms (we happen to know that no spike will start later than 15 ms)
2. Use variable dt (gains a large speedup at the cost of slightly lower accuracy)

Plot result

Bring up a Grapher tool (NEURON Main Menu / Graph / Grapher). You have already read about it but you might want to try this simple example just to brush up:
PlotWhat : sin(t*t)
Plot

To use the Grapher to generate a family of simulations and plot the results, we want to plot

	AlphaSynapse[0].gmax vs taper

If we define "taper" as our independent variable, we need a generator such as

	hhtaper(taper) threshold(&AlphaSynapse[0].gmax)

This is just a "one-liner" (old BASIC hackers will remember one-liners) that first executes the procedure hhtaper(taper) and then applies the threshold() function, adjusting gmax until the threshold value is found.

Alternatively we could leave the Generator field empty and just plot a function like

func g_thresh_taper() {
	hhtaper($1)
	return threshold(&AlphaSynapse[0].gmax)
}

Don't forget that taper should be allowed to range only from 0 to 1.

To save time (you'll be executing a lot of runs!), set the Grapher's Steps parameter to 5 or 10 at most.

Footnotes

Speed

Notice that the threshold function executes a family of runs. Since we are using the Grapher to plot the threshold value as a function of the spatial extent of HH channels, we effectively have a family of families (or, if we changed synaptic location as well, a family of families of families . . . ), which may involve hundreds of runs. Since plotting is often the rate limiting step for a simulation, the entire process can be speeded up considerably by not plotting the soma voltage (in the NEURON Main Menu, toggle "Quiet" on).

Using the Grapher to plot threshold gmax as a function of ion channel distribution

Exercises

Try putting the AlphaSynapse at different points in the model. Turn on Keep Lines in the Grapher so you can compare results for different synaptic locations.