d_lambda setup

[nrnuser]

Hi, I have some issues with the d_lambda setup.

I was calculating the activating potential of a single axon. An extracellular potential was applied to the axon. When I was using the NEURON, I loaded the neuron .swc file by Import3D to the CellBuilder. In the CellBuilder, under Geometry, I used difference values of d_lambda from default 0.1 to 0.0001, and the resulting activating potential changed dramatically. I tried to see the convergence of the numerical error, but the numerical errors started to oscillate around 0.001.

My questions are --
1. Will the d_lambda affect the result's accuracy?
2. How much will it affect the results?
3. Could you give me any suggestion on my calculation of this problem?

Thank you very much!

[ted]

First, some important questions:
What is your operational definition of "activating potential" (by what experimental procedure did you determine its numerical value)?
Can you present a table of your results?

An extracellular potential was applied to the axon.

How? Specifically:
What was your conceptual model of extracellular stimulation (assumptions about geometry [size and shape] of the stimulating electrode(s), spatial relationship between the electrode(s) and the axon, geometry and electrical properties of the conductive medium)?
What was your conceptual model of the axon: length, diameter, presence or absence of myelination (and, if myelinated, node geometry)?
How did you implement extracellular stimulation in your program--by driving each segment's e_extracellular, or by injecting current into segments ("activating function" approach)?

Maybe it would be easier to just present your source code, if it's not too long, with some usage hints and a link to the swc file.

[nrnuser]

Dear Ted, Thank you very much for your quick response! The activating potential is the threshold extracellular voltage at which the axon is depolarized that we consider it as activated. For example, we have 1000 nodes on our axon in the .swc file, and for each node, we have a corresponding extracellular voltage value that we will use. We multiply the extracellular voltage by a scale, and see at which scale value, the axon is depolarized (activated).

In the NEURON simulation, we load the extracellular stimulation voltage that we have obtained outside NEURON, so we do not consider the electrode problem. The axon model is the simple one, without myelin, without branches, length around 150 microns, diameter varies from 0.15 to 3 microns. I used the code similar to your extracellular_stim_and_rec code, but modified it by importing the extracellular voltage values from a file.

My codes are not long. Sharing the codes is a good idea. How do I share the codes with me? Please let me know. Thank you very much!

[ted]

You could just zip up all the source files that are needed to reproduce your results (hoc, mod, ses, swc, txt, dat--but also see viewtopic.php?t=1048) and email them to ted dot carnevale at yale dot edu

[nrnuser]

Dear Ted, Thank you very much for offering help on my issue. I sent the zipped codes to your email address, and I removed the .o, .c, and .dll files as suggested in the link. Thank you very much!

[ted]

Thank you. I'll review it and reply accordingly.

In the meantime I decided to run a test of my own, using a simple model:
cylindrical axon
geometry: diam 1 um, length 1000 um, 0 end at 0,0,0, 1 end at 1000, 0, 0 (all units in um)
biophysics: hh membrane at standard operating temperature (6.3 deg C), cytoplasmic resistivity 100 ohm cm
extracellular medium: infinite, homogeneous, isotropic, with resistivity 35.4 ohm cm
stimulating electrode: monopolar electrode at -100,0,0 (i.e. 100 um "left" of 0 end of axon)
stimulus current: single rectangular pulse 0.1 ms duration starting at t == 1 ms
simulation parameters: baseline spatial discretization nseg = 37 (according to d_lambda rule with d_lambda = 0.1 * length constant at 100 Hz), temporal discretization (dt) 0.025 ms, simulation duration 15 ms

Protocol:

Code: Select all

Repeat
  note nseg
  find threshold stimulus amplitude (in mA) by binary search,
    stopping when (suprathreshold stim - subthreshold stim) is < 1e5 * suprathreshold stim
  increase nseg by a factor of 3
Until nseg > 37 * 3^5

Not that increasing nseg by the factor K reduces spatial error by the factor K^2. This means that tripling nseg reduces spatial error by almost an order of magnitude, so at the finest spatial resolution used in this experiment (nseg = 8991), spatial error is only 1.7e-5 of what it is with nseg = 37

Results:

Code: Select all

 nseg    threshold
37      -0.83311  
*3=111  -0.78491
*3=333  -0.77970
*3=999  -0.77913
*3=2997 -0.77907
*3=8991 -0.77906

I'll make the source code available online.

[ted]

Here's the code that I used. The main file is called test.hoc. It starts with

Code: Select all

// model of monopolar extracellular stimulation of an axon
// stimulus is a 100 us current pulse that drives extracellular potential
load_file("nrngui.hoc")

///// create model cell
load_file("axon.hoc") // an axon with
  // geometry 1 um diam x 1000 um long, 0 end at 0,0,0, 1 end at 1000,0,0
  // biophysics hh, Ra = 100 ohm cm
  // spatial discretization uses d_lambda rule with d_lambda = 0.1

axon.hoc was exported from a CellBuilder.

Next in test.hoc is "instrumentation code" (specifies how the model cell is stimulated). xtra is a density mechanism specified in the NMODL file xtra.mod. xtra.mod and the four hoc files (anatscale.hoc, interpxyz.hoc, setpointers.hoc, and calcrxc.hoc) are all in https://www.neuron.yale.edu/ftp/ted/neu ... nd_rec.zip

Code: Select all

insert extracellular
insert xtra

load_file("anatscale.hoc") // show xyz scale bars

///// set up monopolar stimulation
load_file("interpxyz.hoc") // used by setpointers
  // calculates coords of segment centers
  // and stores in x_xtra, y_xtra, z_xtra

load_file("setpointers.hoc") // calls grindaway in interpxyz.hoc
  // and links extracellular and extra

load_file("calcrxc.hoc") // controls and shows where electrode is located
  // and calculates transfer resistances between electrode and segments 
  // default location is 50,0,0
  // change to -100, 0, 0
XE = -100 
setelec(XE, YE, ZE)

Next in test.hoc is a procedure that simplifies the task of changing nseg. To change nseg to some value foo, at the oc> prompt execute the command

changenseg(foo)

Code: Select all

// changing nseg to a new value will add and/or destroy some nodes
// so setpointers must be called to couple the new nodes to the stimulus
// and setrx must be called to make transfer resistances correct.
// proc changenseg takes care of all this
// argument is the new value of nseg

proc changenseg() {
  nseg = $1
  setpointers()
  setrx(XE, YE, ZE)
}

At this point I decided to make the user's life a little easier by providing a simple waveform generator, which produces a rectangular pulse, to drive is_xtra. For this I used an NMODL-specified class called Fsquare; this is specified by NMODL code in the file fsquare.mod, which is available with a small demo program https://www.neuron.yale.edu/ftp/ted/neuron/fsquare.zip

Code: Select all

///// stimulus

// Fsquare defined in fsquare.mod
// generates a train of one or more biphasic rectangular pulses
// we'll use just the first phase of a single pulse
DEL = 1   // start of first pulse
DP = 0.1  // duration of each phase in ms
NUM = 1   // number of pulses to generate
A1 = -1   // amplitude of phase 1 in mA
A2 = 0    // amplitude of phase 2 in mA

objref fsq
fsq = new Fsquare(0.5)
fsq.del = DEL
fsq.dp = DP
fsq.num = NUM
fsq.a1 = A1
fsq.a2 = A2

setpointer fsq.x, is_xtra

The last bit of code in test.hoc does two things. First, it loads a session file that recreates a RunControl panel with v_init == -65 mV, dt == 0.25 ms, and tstop == 15 ms, and a Graph that plots the axon's membrane potential at its 0, 0.5, and 1 nodes. It also provides code that uses a binary search to determine, within 1e-5 relative precision, the threshold current for triggering a spike. The returned result may be just above, or just below, spike threshold; if the latter, to get a spike simply increase the stimulus amplitude by a factor of 1e-5. Note that effective stimuli will be negative.

Code: Select all

///// simulation flow control

load_file("rig.ses") // RunControl panel and graph of v at 0, 0.5, and 1

// find threshold stimuls to within epsilon

objref apc
apc = new APCount(1) // put at distal end to detect propagating spikes 

// This function requires an existing APCount[0] at the user desired location
// returns 1 if the voltage passed the APCount[0].thresh
func thresh_excited() {
	run()
	return APCount[0].n > 0
}

// do a binary search for threshold
// for this particular model the call is
// threshold(&fsq.a1)
func threshold() { local sub, supra, epsilon
	// the "call by reference" arg $&1 is the independent variable

	// bounding the threshold first can save several iterations
	sub = 0 	supra = 1e5
	if ($&1 == 0) { $&1 = -1 }
	while (sub == 0 || supra == 1e5) {
		if (thresh_excited()) {
			supra = -1*$&1
			$&1 = -1*supra/2
		} else {
			sub = -1*$&1
			$&1 = -2*sub
		}
		if (stoprun) return $&1
		if (sub > supra) return -supra
	}

	// at this point supra is not more than twice the theshold and
	// sub is not less than half the threshold

	epsilon = 1e-9 + 1e-5 * supra	// five place precision
	/* now narrow the bounds */
	$&1 = -(supra + sub)/2
	while ( (supra - sub) > epsilon) {
		if (thresh_excited()) {
			supra = -1*$&1
		} else {
			sub = -1*$&1
		}
		$&1 = -(supra + sub)/2
		if (stoprun) break
	}
	return $&1
}

[ted]

Having read my previous post, someone might well ask
"OK, Mr. Wizard, I have all the source code. How do I proceed?"

First execute nrnivmodl to compile the mod files.

Next create your own rig.ses file.
1. Execute nrngui
2. In the NEURON Main Menu toolbar click on Build/single compartment
3. In the NMM toolbar click on Tools/RunControl. This pops up a RunControl panel. In that panel change the number in the field next to the Tstop button to 15.
4. In the NMM toolbar click on Graph/Voltage axis
5. On your desktop, drag the RunControl panel and the voltage axis graphs to convenient locations.
6. In the NMM toolbar click on File/save session. This pops up a panel that asks you to "Enter filename". Click in the adjacent edit field and enter the string
rig.ses
and then click on the Save button.
7. Exit NEURON.
8. Use a plain text file editor ("a programmer's editor") to open rig.ses. Find the line that starts with

//Begin SingleCompartment

and scan down to the line that starts with

//End SingleCompartment

Delete those two lines and all lines between them.
Save the edited rig.ses and exit the editing program.

Now execute
nrngui test.hoc

To launch a simulation, click on the Init & Run button in the RunControl panel.

To do a binary search for the threshold stimulus, at the oc> prompt execute the command
threshold(&fsq.a1)
NEURON will run some simulations, and then print
-0.83311081
which is the stimulus current that was used for the last run that it executed. You won't see a spike in the graph, so you know the stimulus was subthreshold, but the binary search returns a result that has 5 decimal places of accuracy, so you know that the threshold is between -0.83311 and -0.83312. To check this, at the oc> prompt execute the command
fsq.a1 = -0.83312
and then click on Init & Run. Got it!

To change the axon's spatial discretization, execute
changenseg(somenumber)
where somenumber is an odd natural number (1, 3, 5 . . . ) or the name of a variable whose value is an odd natural number. To see how this affects the estimate of threshold stimulus, execute
threshold(&fsq.a1)

Why use an odd natural number for nseg? Using an odd value for a section's nseg ensures that there will be an internal node at the middle (0.5 location) of that section. This is because NEURON uses the central difference approximation to the second spatial derivative, so each internal node is at the center of a (conceptual) compartment.

It also means that the error related to spatial discretization is inversely proportional to the square of compartment length (which is the same as the distance between adjacent internal nodes). Consequently, increasing a section's nseg by a factor of 3 reduces spatial error by almost an order of magnitude.

[nrnuser]

Dear Ted! Thank you very much for your detailed codes and explanations! They help me a lot! I like them very much!

I have a question on your codes. In the codes, you use CellBuilder to create a axon.hoc file, while I mostly save it as a .ses file. Are those two types of files saved in the same way by using SAVE in the GUI? What are the difference between those two in use? Which one is better or you recommend?

[ted]

The ses file recreates the CellBuilder tool. The hoc file creates the model cell.

If you have gone to the effort of using a CellBuilder to create a model specification, it's a good idea to save it to a ses file in case you want to change something later on (like disretization method, biophysical properties, the name of one or more sections, or maybe even the branched architecture of the model (e.g. create a new subset of sections, move or delete a section or a tree, add new sections or duplicate a tree) or export a cell class from the CellBuilder. So always save your configured CellBuilder to a ses file (in the FAQ list https://nrn.readthedocs.io/en/latest/guide/faq.html look for the answer to the question
I’ve used the NEURON Main Menu to construct and manage models. How can I save what I have done?)

The hoc file contains statements that are much more readable than the ses file, and has tutorial value because it is well organized and shows how many useful things are done. Also using the hoc file is "safer" than using the CellBuilder during exploratory modeling, because that avoids conflicts between the model specification in the CellBuilder and the model specifcation produced by your command line entries.

Confilicts? What conflicts?

During exploratory modeling you are likely to do some things with the GUI, and some things by executing commands at the oc> prompt. Unfortunately, communication between the CellBuilder and NEURON's computational engine is only one way--from the CellBuilder to the computational engine. The CellBuilder will remain ignorant of anything you did at the command line. So if you change a parameter of your model cell via the command line--whether it is a channel density, a reversal potential, or the value of nseg--and then you change some other parameter in the CellBuilder, the CellBuilder will send the entire model specification that it contains to NEURON's computational engine, and that will overwrite whatever parameter change you made via the command line.

You can prevent this by turning Continuous Create OFF, but then you won't be able to use the ChannelBuilder to change any parameter of the model that is in the computational engine without turning ContinuousCreate back ON . . .