Import3D in batch-process mode?

[phendric]

Guys,

As I outlined in another post earlier today (http://www.neuron.yale.edu/phpBB/viewto ... 2161#p8770), I'm using L-NEURON (http://krasnow.gmu.edu/cn3/L-Neuron/index.htm) to generate unique morphologies of cells within a certain morphological class. L-NEURON outputs these virtual cells in SWC format, but I need them in NEURON format to use them.

I know that Import3D has a GUI interface to make importing single SWC cell morphologies into NEURON easy, but does it have a command-line, batch processing mode to facilitate automatically importing lots of SWC files simultaneously? I'm thinking about the time in the near future when I'm going to automatically generate 500 or 5,000 unique cell morphologies in SWC format, and I'm not looking forward to having to import each one into NEURON via the Import3D GUI by hand.

Note: I'm using NEURON under Linux. Someone pointed me to a tool on neuroland.org that does what I'm asking for, but it doesn't work very well, and it's Windows only.

Phillip

[hines]

The following code will do mostly what you want. I got the swc file from downloading
http://neuromorpho.org/neuroMorpho/dabl ... 61.CNG.swc
and removing all the type 10 lines at the end of the file.

Code: Select all

{load_file("nrngui.hoc")}
{load_file("import3d.hoc")}

begintemplate Cell
public soma, axon, dend, apic         
create soma[1], axon[1], dend[1], apic[1]
public all, somatic, axonal, basal, apical 
objref all, somatic, axonal, basal, apical
proc init() {
    all = new SectionList()
    somatic = new SectionList()
    axonal = new SectionList()
    basal = new SectionList()
    apical = new SectionList()
}
endtemplate Cell

// $s1 swc morphology file name
// on exit the return object is a Cell instance with
// the morpology specified by the $s1 file
obfunc mkcell() { localobj import, morph, cell
    cell = new Cell()
    morph = new Import3d_SWC_read()
    morph.input($s1)
    import = new Import3d_GUI(morph, 0)
    execute("forall delete_section()", cell)
    import.instantiate(cell)
    return cell
}

objref cell, shape
cell = mkcell("c10261.CNG.swc")
shape = new Shape(cell.all)

This makes use of an unusual feature of NEURON which allows sections to be deleted within a cell object and new ones created. The limitation is that you cannot create new names (that is wny you must remove the type 10 lines from the swc file --- those would have created dendrite sections with a special name unknown to Cell.)

Clearly the Cell class can be endlessly elaborated to specify membrane properties and proper compartmentalization (nseg) using the section lists. The Blue Brain project uses an analogous style which creates cells from a cell type file, a membrane file, and a morphology file. They need to create 10k cells using dozens of cell types and membrane templates and hundreds if not thousands of 3-d neurolucida reconstructions. That is why I added a "quiet" variable to Import3d_Neurolucida3

[R_Jarvis]

Hi I was just wondering.

After you have succeeded in batch importing a lot of morphologies and assuming that you had managed to check them all for artifact errors somehow. How would you then go on to add subsets and biophysics?

I think you are alluding to the Blue Brain Project batch importing biophysics files as well as morphology files.
are the import statements like those below which import biophysics, rather than morphology?

morph = new Import3d_SWC_read()
morph.input($s1)
import = new Import3d_GUI(morph, 0)

I think I could use load_file(""), where the argument of load_file specificies a .hoc file that defines biophysics somehow. But I am unsure how to get this information incorporated into the new cell object in the statement:

cell = new Cell()
etc.

Thanks again for any help or suggestions.

[hines]

Your question is posed at quite a high design level and so I would answer by pointing you to the hoc file constructed by the CellBuild tool.
The cell specification has a procedure that specifies the channel insertion and channel parameters in terms of the cell's subsets, e.g. all, apical
or whatever you created and specific sections such as the soma. As long as every cell type you create has at least those same section lists and
section names, you can factor out that procedure and add an objref arg. Then every statement which iterates over a section list or pushes a named section
would just be prefixed with a $o1. ie.
forsec all insert hh
would become
forsec $o1.all insert hh
and
soma { gnabar = ... }
would become
$o1.soma { gnabar = ... }

You could then have a bunch of such procedures and use them at will for any cell.

[R_Jarvis]

Hi,

Thanks for your reply.

I wouldn't be able to write the body of the procedure that you mention above, as I will explain below.

I imported a 3d morphology, and then designed subsets in the GUI, then exported to .hoc. But the code that was exported was not that helpful, because I want to write code that can generalise to any array length for Cell[0].axon and Cell[0].dend. The exported code automatically knows how long the axon and dend arrays are but I don't, and it is dynamically accessing these array lengths.

I am still trying to understand how I can iterate through all the cells and their sections. Is there a function that can give me the length of a dend array? That way I can iterate over various lengths of the dendrite in the for loop and assign different biophysics according to where I am in the length of dendrite.

I would like to be able to assign biophysics like the following:

array_length=Cell[0].axon.length()
for i=0, array_length Cell[0].axon {insert hh}

And to use the above on different lengths of dendrites as well, maybe in a nested for loop that iterates cells, axons, dends, everything.

Are the statements forsec and forall designed to avoid manual iteration through arrays?

>forall Cell[14].axon {insert hh}

Also I would like to know how come the following yields a syntax error

oc>Cell[0].dend[0]
/usr/local/nrn/i686/bin/nrniv: stack empty
near line 93
Cell[0].dend[0]

And the following is permitted?

oc>access Cell[0].dend[0]

Also to go back to the original response you give to batch import of morphology question, I wonder if there is a reason why you chose to make an array of cell objects rather than appending each new cell objects to a list of cells? Maybe you where anticipating that an array of cells would have advantages in this context?

I have been searching a lot for common conceptual errors regarding hoc syntax and errors. Are there any documents that address misconceptions?

Will there ever be a NEURON course in Australia?

Thanks again for your help.

[hines]

One of the hoc procedures emitted by the CellBuild tool is the subsets definition. If you have a bunch of axons then I would recommend constructing a subset within the cell builder
called axons, and define it as all the sections containing the name axon (note: the cell builder would emit code like
axons = new SectionList()
forsec "axon" { axons.append() }
If the subset definition procedure is executed after you fill in the new morphology, the axons subset will have all the axon sections.
However such statements must be executed within the context of an object so that they refer only to the sections within that object.
After a SectionList is specified it can be used anywhere.

Yes, forsec and forall provide flexible iteration over sections with similar names.

I kind of agree that an array should have size property accessible from hoc. However, sometimes section arrays have holes (from specific delete_section() calls.
If you subsets are inhomogeneous, I would consider using the SubsetDomainIterator concept available in the CellBuild tool.
http://www.neuron.yale.edu/neuron/stati ... tline.html

The point I am trying to make is that if you get a hoc file from a CellBuild tool that has all the subsets needed to define the biophysics, You can make an instance of the cell,
throw away the existing sections, import a 3-d morpology into that instance (as long as you do not need new base names), and then call the subsets procedure to create the subsets
and then the biophysical procedures either for constant parameters over subsets or inhomogeneous functions over the subsets.

Cell[0].dend[0]
is a syntax error because the proper syntax is
Cell[0].dend[0] stmt
where stmt may be a multiline statement enclosed in {}
The
access section
is a complete statement

Generally I prefer lists of cells instead of arrays of cells. But nothing wrong with arrays as long as you explicitly have the size tag along with the array when passing as an arg to a function.

I wish there was a document that addressed misconceptions. Let this be an open call for help in assembling a list.

We'd be delighted to teach a 3 or 5 day course in Australia. However I'm not sure I could survive an economy flight.

[R_Jarvis]

Hi, your suggestions are giving me useful ideas to test but I am still thoroughly confused on this issue.

Below (1) is the relevant code that was exported out of the cellbuild tool after I imported a .swc, and then exported to cell builder made subsets and saved as .hoc

Based on this code, I have redesigned a template further below (2) that also includes subsets and biophysics, then further below again I made a procedure (3) which tries to change the biophysics of sections in already instantiated cells. None of the approaches work.

(Note in the template below I have not designed the biophysics to be realistic I only inserted different mechanisms for each subsection so that they where distinuishable for debugging subsets and biophysics togethor).

Code: Select all

//(1)

proc subsets() { local i
  objref all, somatic, axonal, basal, apical
  all = new SectionList()
    soma all.append()
/*   for i=0, 44 dend[i] all.append()
    for i=0, 321 axon[i] all.append()
    for i=0, 16 apic[i] all.append() */

  somatic = new SectionList()
    soma somatic.append()

  axonal = new SectionList()
//    for i=0, 321 axon[i] axonal.append()

  basal = new SectionList()
  //  for i=0, 44 dend[i] basal.append()

  apical = new SectionList()
    //  for i=0, 16 apic[i] apical.append()

}

proc biophys() { 
  forsec somatic {
    insert Nasoma
    insert somacar
  }
  forsec axonal {
    insert Kaxon
  }
  forsec basal {
    insert nmda
  }
  forsec apical {
    insert pas
      g_pas = 0.001
      e_pas = -70
    insert Nadend
  }
}
access soma

Below I used the general idea from above and tried to make some subsets that might be generalisable in the cell template below.
I think what you are saying is. Since the morphology and cell template are separated, you can make an abstract subsets in the template and no indexs are required to be known, since the morphology doesn't exist at this point.

And then when the morphology is imported at about the same point the cell is instantiated, the sections from the morphology can be made to belong to abstract subsets from the cell template.

I wonder why the line: execute("forall delete_section()",cell) at (2) is necessary since no sections have been made yet?

Code: Select all

obfunc mkcell() { localobj import, morph, cell
	
		    cell = new Cell() 
	 	
		    morph = new Import3d_SWC_read()
		    morph.input($s1)
		    import = new Import3d_GUI(morph, 0)
         		 execute("forall delete_section()", cell) //(2)
		    import.instantiate(cell)
 
        	    return cell
		}

objref cell1, cell2, CELLSgr[45], CELLSp[73], CELLSbask[55], CELLSgl[1],  CELLSch[3], shape



/* Granule cell 
Turner 
 Rat 
 Hippocampus */

CELLSgr[0]=mkcell("n220.CNG.swc") 
//////////////////////

Code: Select all

//The code above is instantiating a template like the following:


//(2) 
begintemplate Cell
public soma, axon, apic, nclist, dend, syn         
create soma[1], axon[1], dend[1], apic[1] 
objectvar f, nclist
public all, somatic, axonal, apical, dends, basal 
objref all, somatic, axonal, apical, dends, basal

proc biophys() {

  forall { insert hha2}
  
  somatic {
    insert Ih
  }
  axonal {
    insert Kaxon
  }

  dends{
    insert IA
  }
  
   basalar {
   insert Nadend
    
  }
  apical {
    insert hha2 
  }
 
}


geom_nseg() {
   forsec somatic { nseg = 5  }
   forsec axonal { nseg = 5  }
   forsec basal { nseg = 5  }
   forsec apical { nseg = 5  }
//   soma { nseg = 10  }
}

proc subsets() {objref all, somatic, axonal, basalar, apical, dends
   all = new SectionList()
 
  somatic = new SectionList()
   soma somatic.append()
  soma all.append()
 
  axonal = new SectionList()
  axon axonal.append()
   axon all.append()

   dends = new SectionList()
   dend dends.append()
   
  apical = new SectionList()
  apic apical.append()
  apic dends.append() 
  apic all.append()

  basalar = new SectionList()
  basal basalar.append()
  dend dends.append()
  basal all.append()
}





proc init() {
  somatic = new SectionList()
  axonal = new SectionList()
  basalar = new SectionList()
  dends = new SectionList()
  apical = new SectionList()
  nclist = new List()
  geom_nseg()	    
  subsets() 
  biophys() 
	  	
}

After I have run the batch import script I try the tests:

Code: Select all

oc>Cell[175].apical psection()
oc>Cell[175].dends psection()
 

both return information not about cell[175] but instead about Cell[0].soma which somehow got all of the mechanisms inserted into
it from all of the subsets above, not just the somatic one which is puzzeling.



An example from a similar file produced:

Code: Select all

oc>Cell[175].axonal psection()
yields:
Cell[0].axon[0] { nseg=100  L=100  Ra=35.4
	/*location 0 attached to cell 0*/
	/* First segment only */
	insert capacitance { cm=1}
	insert morphology { diam=500}
	insert Kaxon {}
	insert Nadend {}
	insert IA {}
	insert Ih {}
	insert hha2 {}
	insert hh { gnabar_hh=0.12 gkbar_hh=0.036 gl_hh=0.0003 el_hh=-54.3}
	insert na_ion { ena=50}
	insert k_ion { ek=-77}
}

I also made a procedure outside the cell template that is supposed to insert mechanisms into already instantiated cell objects however it doesn't really work.

Code: Select all

//(3)
proc dendrites2() {
	        forsec $o1.dends { insert hh } 
		forsec $o1.dends { insert Kaxon }
		    $o1.dends { insert Naaxon }
		    $o1.dends { insert hNa }
		    $o1.dends { insert hh }
		    $o1.dends { insert Ih }
		    $o1.dends { insert Kdend }
}

//function call to test
for i=170,175 { dendrites2(Cell[i]) }

oc>Cell[170].dends psection()

If I used an approach like this on already instantiated objects I don't understant how I could do it without knowing all of the indexs. As I don't know any statements that are flexible enough to accomodate these changing size section arrays without using indexs.

I have noticed that the statement :

Code: Select all

		forall $o1.dends { insert Kaxon }

reliably changes every section listed in psection(), but it is a bit useless to do that since you lose the specifity of subsets which is obviously the whole point.


Also I can't find any .hoc code for the implementation of the topology() function. I thought this would be helpful since it somehow is able to iterate the length of section arrays without exceeding the size and causing an error as a corollary. Am I wrong to assume that topology() would be in the directory:

/usr/local/nrn/share/nrn/lib/hoc/

?

I tried $grep -r 'topology' /usr/local/nrn/share/nrn/lib/hoc/

But I couldn't find anything relevant.

I suspected that the file might be similar:

celbild/celtopol.hoc

?

Thanks so much for anyhelp.
R

[hines]

I wonder why the line: execute("forall delete_section()",cell) at (2) is necessary since no sections have been made yet?

When Cell is instatiated, the full cell is created acording to the template and that indeed does create sections. See the
create statements in the template (2nd line after begintemplate Cell). So you have one soma, axon, dend, and apic.

oc>Cell[175].apical psection()

That is syntax for a section, not a SectionList. So you are getting the default currently accessed section. Try:
forsec Cell[175].apical psection()

topology() is implemented in c in src/nrnoc/solve.c
Subsets from 3-d reconstructions are often indicated in the 3-d format files. Is there no indication in the swc file which sections
are apical, basal, etc. If there is, then import3d.instantiate would create the sections, specifiy the 3-d points, connect the sections,
and create the subsets for sections indicated as somatic, axonal, basal, apical, dendritic, and, for unknown types, "minus_<-type>set"

[R_Jarvis]

Some of the reconstructions in .CNG.swc from neuromorpho.org are missing axons.

If I have an array of cell objects of say ca1 pyramidal cells can I replicate the axon from one cell object, and connect it to an instantiated reconstructed ca1 pyr neuron with no axon? I would try to edit the axon lines of the .CNG.swc file, but not all the reconstructions are translated to the same (0,0,0) origin.

I also tried importing a cell that did have an axon into cell builder and then deleting everything except the axon, and then saving that as a hoc file and adding it to various reconstructions that don't come with axons with a statement like:

connect Cellsx[0].soma[0](1), axon[0](0)

But it doesn't seem like a very good solution, I wonder if anyone has any better ideas for trying to keep as many different cell morphologies as possible by repairing incomplete neurons?

Thanks for any suggestions.

[hines]

can I replicate the axon from one cell object, and connect it to an instantiated reconstructed ca1 pyr neuron with no axon?

Yes. But I recommend keeping all sections of a cell "owned" by the cell. I.e created in the contest of the
cell object. That is how the "proc instantiate" of import3d_gui.hoc does it. ie. it bullds statement strings
and executes them in the context of the cell using execute(stmt, obj) .
It is generally important for the sections in a cell to be appended to the proper cell section lists.

[R_Jarvis]

I found the procedure in /opt/nrn/share/nrn/hoc etc.

The file I am using is a hoc file called: "just_axon.hoc" I am having trouble relating the imported contents of this file (visible as sections from topology) to tstr, below.

Should I make an object template for the geometry/subsets/topology in just_axon.hoc and then do something with morph.execute() as suggested below?

proc instantiate() {local i, j, min localobj sec, xx, yy, zz, dd
chk_valid()
min = set_nameindex(swc.sections)
// create
for i = 0, types.size-1 {
type2name(i+min, tstr)
if (types.x == 1) {
sprint(tstr1, "~create %s[1]\n", tstr)
execute(tstr1, $o1)
}else if (types.x > 1) {
sprint(tstr1, "~create %s[%d]\n", tstr, types.x)
execute(tstr1, $o1)
}
if ($o1 != nil) { mksubset($o1, i+min, tstr) }
}
if ($o1 != nil) {execute("forall all.append", $o1) }
// connect
for i = 0, swc.sections.count-1 {
sec = swc.sections.object(i)
if (sec.is_subsidiary) { continue }
name(sec, tstr)
if (i == 0) {
sprint(tstr1, "access %s", tstr)
if ($o1 == nil) {
execute(tstr1, $o1)
}
}
if (sec.parentsec != nil) {
name(sec.parentsec, tstr1)
sprint(tstr1, "%s connect %s(0), %g", tstr1, tstr, sec.parentx)
execute(tstr1, $o1)
}
// 3-d point info
if (sec.first == 1) {
sprint(tstr1, "%s { pt3dstyle(1, %g, %g, %g) }", tstr, sec.raw.x[0][0], sec.raw.x[1][0], sec.raw.x[2][0])
execute(tstr1, $o1)
}
j = sec.first
xx = sec.raw.getrow(0).c(j)
yy = sec.raw.getrow(1).c(j)
zz = sec.raw.getrow(2).c(j)
dd = sec.d.c(j)
if (sec.iscontour_) {
contour2centroid(xx, yy, zz, dd, sec)
}
for j = 0, dd.size-1 {
sprint(tstr1, "%s { pt3dadd(%g, %g, %g, %g) }",\
tstr,xx.x[j], yy.x[j], zz.x[j], dd.x[j])
execute(tstr1, $o1)
}
}
}

So in the obfunc mkcell. Instead of finishing with import.instantiate(cell)



obfunc mkcell() { localobj import, morph, cell
cell = new Cell()
morph = new Import3d_SWC_read()
morph.input($s1)
import = new Import3d_GUI(morph, 0)
execute("forall delete_section()", cell)
import.instantiate(cell)
//
//proposed new lines:
objref obj
obj = new just_axon()
morph.execute(stmt, obj)
//
return cell
}

can I call execute(stmt, obj) with an object made out of a template from the axon geometry/topology .hoc file?

Thanks.
Russell.

[hines]

In looking at your last fragment, i'm clearly misunderstanding your intention. I believe I understand up through the import.instantiate(cell)
statement. I.e. you replaced the original morphology in your cell with the morphology in the $s1 swc file. Then I get confused. You build a new
just_axon instance and then you call a non-existent execute function on the Import3d_SWC_read instance. Did you want something like
obj.copy2cell(cell) in which the just_axon instance would create an axon in cell and copy the just_cell's axon properties into the cell's axon?
In that case, The just_axon type would not even have to have its only copy but just a recipe for building the axon in cell. In fact just_axon could
be a proc with parameters.