adding new template slows down

[pfortier]

I generated the following template (it works well):

Code: Select all

begintemplate tpTCN 
	external newseed, grDur, stimulatortime
	public soma,vcell,gnoise,electrode,stimulator
	create soma
	objref vcell,gnoise,electrode,stimulator
	
	proc init() { 
		access soma
		soma {
			
			nseg = 1
			diam = 76.58
			L = 100
			cm = 1
			
			insert hh2	
			ena = 50
			ek = -100
			vtraub_hh2 = -52
			gnabar_hh2 = 0.01	
			gkbar_hh2 = 0.01 	
			
			insert itGHK
			cai = 2.4e-4 
			cao = 2 
			eca = 120 
			shift_itGHK = -1	
			gcabar_itGHK = 0.0002
			qm_itGHK = 2.5
			qh_itGHK = 2.5
			pcabar_itGHK = 8e-5
			
			insert cad
			depth_cad = 0.1
			kt_cad = 0 
			kd_cad = 1e-4
			taur_cad = 5
			cainf_cad = 2.4e-4	
			
			insert pas
			g_pas = 3.79e-5
			e_pas = -76.5
			cm = 0.88
			Ra = 173
			
			gnoise = new Gfluct(0.5) 
			gnoise.g_e0 = 0.0711
			gnoise.std_e = 0.005 
			gnoise.g_i0 = 0.1
			gnoise.std_i = 0.02
			gnoise.tau_e = 3 
			gnoise.tau_i = 10
			gnoise.E_e = 0
			gnoise.E_i = -75
			gnoise.new_seed(newseed())
			
			electrode = new PointProcessManager(0) 
			electrode.pp=new IClamp(0)
			electrode.pp.del=0
			electrode.pp.dur=grDur
			electrode.pp.amp=0 
		}
		
		vcell = new Vector()
		vcell.record(&soma.v(0.5))
		
		stimulator=new VectorPlay(1)
		stimulator.vy = new Vector(grDur+1) 
		stimulator.vx = stimulatortime 
		stimulator.vy.play(&electrode.pp.amp, stimulator.vx)
	}
endtemplate tpTCN

and then I created many instances of this cell template using the following loop:

Code: Select all

rowsTCN=30
colsTCN=rowsTCN 
slabsTCN=2
objref TCN[colsTCN][rowsTCN][slabsTCN] 
for col=0, colsTCN-1 for row=0, rowsTCN-1 for slab=0, slabsTCN-1 {
	TCN[col][row][slab] = new tpTCN()
	TCN[col][row][slab].gnoise.g_e0 = .05
	print "created TCN[",col,"][",row,"][",slab,"]"
}

I didn't put the entire code of the program. The point is that the above loop for creation of instances of the cell starts executing rapidly (about 100 print statement lines per second) but then quickly slows down (1 print statement line per 20 seconds). The above 30x30x2 loop takes hours to complete.

Is this normal?
Is there a way to speed this up?

[ted]

First thing to do, before running any more tests: cut way back on the number of rows, columns, and slabs. You need to be able to find out quickly whether any revision improves the speed of model setup.

There are several possible explanations why this code executes slowly; the difficulty is to identify those that are most likely. My leading guess is that the problem is caused by using thousands of instances of the PointProcessManager and VectorPlay tools. Like most of the GUI's tools, these are meant to be used for interactive, exploratory simulations, in which there are only a very few instances of either of them.

What is the utility of creating any PointProcessManagers if you're (a) not going to display them on the screen, and (b) you're only going to use them as IClamps? Why not create IClamps in the first place? They're "lightweight" compared to PPMs.

Similar comments apply to the VectorPlay tool. I can't imagine a reason why anyone would need more than a couple of these things, and then only for a toy model involving at most two or three cells. Just use the Vector class's play method to drive IClamp.amp directly.

Do the IClamps have to deliver a waveform that varies at every time step, or are they just used to deliver something very simple like a pulse or a series of pulses? If the latter, your t and i Vectors only need to contain the values at the "break points" (e.g. a 1 nA x 1 ms pulse that starts at 3.1 ms is described completely by these 6 pairs of t, i values:

Code: Select all

t i
0 0
3.1 0
3.1 1
4.1 1
4.1 0
5 0

(the final pair is needed by adaptive integration when Vector play with interpolation is used (vsrc.play(&var or stmt, Dt or Tvec, continuous) where continuous == 1)).

Progress indicators are nice, but sometimes they waste time. You can probably gain some speedup by not printing anything after each cell has been constructed.

A comment about code organization:
It is a violation of an important software design principle (modularity) and a source of needless difficulty to include instrumentation and control code in the template that defines a cell class. It is best to keep "specification of biology" separate from "specification of instrumentation" and "specification of simulation control." The template should contain only the statements that create the section and specify its anatomical and biophysical properties. IClamps are instrumentation, and should be created by procedure(s) that are not part of the template. Vector play can be regarded as either instrumentation or control code (because it defines an "experimental protocol") and should likewise be handled by procedure(s) that are not mixed into the template.

If you reorganize your program in a modular fashion

Code: Select all

template that defines cell class
create cell instances
attach IClamps
set up Vector play

and test in an incremental manner

Code: Select all

1.  create a program that does this
  template that defines cell class
  create cell instances
and test it

2.  add code that attaches IClamps and run a new test

3.  add code that sets up Vector play and run a new test

you will discover for yourself the source of any bottlenecks that you encounter, at a stage in program development when they are fixable with minimal effort.

[pfortier]

Thanks for the design tips.

The exponential slowing down of new instances was due to the multiple use of PointProcessManager. I used this because I was reproducing an example I found of VectorPlay. Replacing the PointProcessManager with IClamp eliminated the exponential slowing down. I didn't know that I could use the Vector class's play method to drive IClamp.amp directly. Replacing VectorPlay with the Vector class's play method increased execution speed.

I need to stimulate and save the voltage response of all 30x30x2 neurons (all graphs can be done after the simulation). Each cell receives the same sine wave but the onset of stimulation is at a different delay for each cell. To reduce memory requirements by having a single sine wave (rather than 30x30x2 copies with different delays), I tried stimulator.play(&electrode.amp, 1) and set electrode.del=delay but the delay is ignored (seems to make sense that stimulator.play would override or ignore the electrode.del variable setting). I'm now testing the following algorithm that would read from a single copy of the stimulation sine wave (stored as double array) at the appropriate delay and save the voltage responses to a multidimensional array (stored as double array).

Code: Select all

proc run() {
	stdinit()
	while (t<grDur) {
		record voltage for all cells
		set amp for all cells
		fadvance()
	}
	record voltage for all cells
}

Do you see problems with this approach or a better approach?

[ted]

I didn't know that I could use the Vector class's play method to drive IClamp.amp directly.

Strong suggestion: read the Programmer's Reference. Browse around in it now and then, even if you don't have a specific problem in mind. Contains lots of useful stuff, likely to spark new ideas.

I need to stimulate and save the voltage response of all 30x30x2 neurons (all graphs can be done after the simulation). Each cell receives the same sine wave but the onset of stimulation is at a different delay for each cell. To reduce memory requirements by having a single sine wave (rather than 30x30x2 copies with different delays), I tried stimulator.play(&electrode.amp, 1) and set electrode.del=delay but the delay is ignored

delay specifies the start time of the current pulse generated by the IClamp (described in the Programmer's Reference).

I'm now testing the following algorithm . . . Do you see problems with this approach or a better approach?

Way too complex. Better to write a mod file that defines a point process or a density mechanism that delivers a sinusoidal current. Easier to use (each cell gets one of its own), and imposes much less storage and computational overhead.

This is not something you'll find in the Programmer's Reference, so tell me exactly what you need and I'll create it for you.
Do you want to specify frequency in Hz or radians/s? Phase in degrees or radians?
Do you want the current waveform to start at t = 0, an arbitrary time >=0 , or at the time of the first positive-going zero crossing?
Do you want it to continue forever, or only for a finite time? If the latter, do you want to specify a particular duration of current delivery, a particular end time, or a particular number of cycles?

[pfortier]

A wonderful feature of NEURON is that it allows simulations to be coded in many different ways. Unfortunately, when I find one way that works, I often forget to look for a better way. I'll try and develop the reflex to look in the manuals for alternate ways.

Thank you for offering to write a mod file. I think that it should offer the following parameters:

sineobj.freq = // frequency in radians/s
sineobj.phase = // phase in radians
sineobj.del = // delay in ms (if phase=0 and del=100 then values start from 0 at 100 ms)
sineobj.dur = // duration of stimulation in ms
sineobj.amp = // amplitude of stimulation in nA (I'm assuming that the sine function would give values between 1 and -1 and that this would be multiplied by amp; for amp=0.2 then the sine wave would vary between 0.2 and -0.2 and for amp=-0.3 then the sine wave would vary between -0.3 and 0.3)

[ted]

At the risk of making it look too easy (but then, Vladimir Horowitz always made it look easy), here's the code.

Code: Select all

COMMENT
isin.mod
Delivers sinusoidal current i = amp*sin(f*(t - del) + phi)
starting at t==del and ending at t==del+dur
where f is frequency in radians/s and phi is phase in radians
3/13/2009 NTC
ENDCOMMENT

NEURON {
  POINT_PROCESS Isin
  RANGE del, dur, f, phi, amp, i
  ELECTRODE_CURRENT i
}

UNITS {
  (nA) = (nanoamp)
  PI = (pi) (1)
  (Radian) = (1)
}

PARAMETER {
  del (ms)           : time at which current starts
  dur (ms) <0, 1e9>  : duration
  f (Radian/s) <0, 1e9> : frequency in radians/s
  phi (Radian)       : phase shift in radians
  amp (nA)
}

ASSIGNED {
  i (nA)
  on (1)
}

INITIAL {
  i = 0
  on = 0
  net_send(del, 1)
}

FUNCTION fsin(t (ms)) {
  fsin = sin(f*(t - del)/(1000) + phi)
}

BREAKPOINT {
  i = on*amp*fsin(t)
}

NET_RECEIVE (w) {
  : ignore any but self-events with flag == 1
  if (flag == 1) {
    if (on == 0) {
      : turn it on
      on = 1
      : prepare to turn it off
      net_send(dur, 1)
    } else {
      : turn it off
      on = 0
    }
  }
}

My bill is in the mail.

[pfortier]

Simplicity can be very attractive. I would not have ended up with this code had I based it on IClamp (stim.mod) because it uses at_time(). I had to do some reading to figure out what net_send(), NET_RECEIVE() and flag were. Thank you very much, it works great.

To summarize, I started with a working simulation that took 10 hours to run because I used 1800 instances of the PointProcessManager and VectorPlay objects to provide different stimulus patterns to 1800 neurons. The time was shortened to minutes by replacing the PointProcessManager with IClamp (creating multiple instances of PointProcessManager was the source of the EXPONENTIAL SLOWING of the simulation) and the VectorPlay with the play method of Vector objects. The memory requirements were dropped enourmously by replacing the 1800 stimulus vectors with a mod file. Wow, this was a great learning experience. Thanks again!

[ted]

The event delivery system makes it easy to implement "state machines" with NMODL and/or hoc that are useful for all sorts of nifty purposes.

Glad to be of help. Thank you for using NEURON in your research! Please be sure to let me know of any papers you publish that report work done with NEURON, so we can include them in the on-line bibliography.