Different total area from python and hoc cell class

[opdalex]

Hi,

I have stumbled upon an issue with different measures of total area when trying to convert an old method to using a python class to represent a neuron.

The old method uses the NEURON gui; from there using import 3d for importing a Neurolucida reconstruction, then exporting to CellBuilder to edit geometry, biophysics and saving as a cell class file in hoc.

The alternative method avoids the gui through creating a cell class in python, which takes a Neurolucida reconstruction as a parameter to import the morphology and initialize it. The construction of a cell also call a helper function which sets the same geometrical and biophysical properties as when using CellBuilder in the older method.

I realize that there is already a difference between the methods here; the older uses CellBuilder, while the python class uses instantiate and edits the geometry and biophysics thereafter. I have attempted to use CellBuilder without a gui, but as I could not make this work, I attempted to stick closely to the python methodology used in the NEURON tutorial using the "BallAndStick" class as a framework. Therefore, I used instantiate on the imported morphology, and expected this to suffice.

However, the results were not the same when measuring the total area of morphology of the different methods.

This is the hoc cell class result:

Code: Select all

biomed0220678:~ alex$ python3
Python 3.9.17 (main, Jun 10 2023, 21:06:08) 
[Clang 9.1.0 (clang-902.0.39.2)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from neuron import h, gui
>>> h.load_file("/Users/alex/Documents/Master_Project_code_files_Alexander/m210309_2_class.hoc")
1.0
>>> h('objref cell')
1
>>> h('cell = new Cell()')
1
>>> cell = h.cell
>>> area = totalarea(h, cell)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
NameError: name 'totalarea' is not defined
>>> def totalarea(h,cell):
...     sum = 0
...     for sec in cell.all:
...          for seg in sec.allseg():
...              sum += seg.area()
...     return sum
... 
>>> area = totalarea(h, cell)
>>> area
2022.6127147520217
>>> 

This is the python cell class result:

Code: Select all

biomed0220678:~ alex$ python3
Python 3.9.17 (main, Jun 10 2023, 21:06:08) 
[Clang 9.1.0 (clang-902.0.39.2)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from neuron import h, gui
>>> import sys
>>> import os
>>> sys.path.append(os.path.abspath("/Users/alex/Documents/Master_Project_code_files_Alexander/new/scripts/import_scripts"))
>>> from AII_class import AmacrineAII as AII
>>> ASC_file = '/Users/alex/Documents/Master_Project_code_files_Alexander/fixed_ASC_files/m210309_2_v045fix.ASC'
>>> cell = AII(ASC_file)

4163 lines read
22 spines

/Users/alex/Documents/Master_Project_code_files_Alexander/fixed_ASC_files/m210309_2_v045fix.ASC problems


Main branch starting at line 274 is outside the soma bounding boxes
  Making a logical connection to center of nearest soma
>>> 	-65 
def totalarea(h,cell):
...     sum = 0
...     for sec in cell.all:
...          for seg in sec.allseg():
...              sum += seg.area()
...     return sum
... 
>>> area = totalarea(h, cell)
>>> area
2022.6125052233679
>>> 

AII class:

Code: Select all

import math
from neuron import h

h.load_file("stdlib.hoc")
h.load_file("stdrun.hoc") #for running simulation without gui
h.load_file("import3d.hoc") #for importing 3d morphology

class AmacrineAII:
    def __init__(self, morphology):#, gid, pos, theta):
        #use only name of specific file without .ASC suffix, not entire file path
        self.name = morphology[morphology.rfind('/') + 1 : -4]
        self.morphology = morphology
        self.load_morphology()
        self.discretize()

        self.add_channels()
        self.x = self.y = self.z = 0
        h.define_shape()
        #self.rotate_z(theta)
        #[x, y, z] = pos
        #self.set_position(x, y, z)
        #self._gid = gid

    def add_channels(self):
        #biophysical properties
        g_pas = 3e-5
        e_pas = -60
        Ra = 200
        Cm = 1

        #make cell with passive channels over the entire cell
        pas_locations = [sec for sec in self.all]
        h.pas.insert(pas_locations)
        for sec in pas_locations:
            sec.Ra = Ra
            sec.cm = Cm
            for seg in sec:
                seg.pas.g = g_pas
                seg.pas.e = e_pas

    def discretize(self):
        #geometrical finegrainedness
        d_lambda = 0.1

        #using d_lambda acquires good discretization often
        frequency = 100 #Hz
        for sec in self.all:
            #use hoc class standard nseg???
            nseg = math.ceil((sec.L / (d_lambda * h.lambda_f(frequency))) / 2) * 2 + 1
            if sec.nseg % 2 == 0:
                nseg += 1
            sec.nseg = int(nseg)


    def __str__(self): #string representation of object
        return self.name

    def load_morphology(self):
        cell = h.Import3d_Neurolucida3()
        cell.input(self.morphology)
        i3d = h.Import3d_GUI(cell, False) #import 3d
        i3d.instantiate(self) #construct object in neuron
        #use i3d.cellbuilder possible???

    def __repr__(self):
        return '{}[{}]'.format(self.name, self._gid)

    def set_position(self, x, y, z):
        for sec in self.all:
            for i in range(sec.n3d()):
                sec.pt3dchange(i,
                               x - self.x + sec.x3d(i),
                               y - self.y + sec.y3d(i),
                               z - self.z + sec.z3d(i),
                              sec.diam3d(i))
        self.x, self.y, self.z = x, y, z

    def rotate_z(self, theta):
        """Rotate the cell about the Z axis."""
        for sec in self.all:
            for i in range(sec.n3d()):
                x = sec.x3d(i)
                y = sec.y3d(i)
                c = h.cos(theta)
                s = h.sin(theta)
                xprime = x * c - y * s
                yprime = x * s + y * c
                sec.pt3dchange(i, xprime, yprime, sec.z3d(i), sec.diam3d(i))

I could also attach the Neurolucida file and the generated class file, but I am not sure where that would be appropriate to do (I am new here).

Hope to hear from someone soon!

[ted]

First, a question. What version of NEURON are you using?

Next, a suggestion. Use as simple a workflow as possible, that is, import the morphometric data, but _don't_ discretize the model cell, and don't add or alter any biophysical properties. What results do you get?

[opdalex]

1. I am using NEURON version 8.0.0 (compatibility with both python 2 and 3 is preferred to newest version)
2. I get the same results as before; 2022.6125052233679 is still the area using the python class.

[ramcdougal]

tl;dr: a difference of this magnitude is expected and should not meaningfully change your results. The Python version is ever so slightly more accurate.

There is a subtle difference in the implementation between HOC and Python that can create these tiny differences (here less than 0.001 µm^2):

In HOC, there are a couple of different code paths to instantiate a morphology, but they all involve taking perfectly good double-precision numbers, turning them into strings with only a relatively few digits of precision, then turning them back into double-precision before storing them in NEURON's internals as the 3D structure of the cell.

In Python, once the numbers are loaded from the file they stay as double-precision with no intermediate string conversion, before storing them in NEURON's internals.

To make matters messier still, NEURON's internals actually only store the coordinates at floating-point precision (not double).

[opdalex]

Thanks a lot, that cleared everything up nicely!

[ted]

Algorithmic details aside, there is one hard fact that should be kept in mind: most parameters of biological systems are known to only very limited precision. For a multitude of reasons, the precision of measurements of neuronal morphology obtained with light microscopy is at best about 0.2 um, but seldom as good as half a micron. This is a serious problem because neurites with diameter <= about 1.5 um account for most of the surface area of cells that have extensive axonal and dendritic fields. Diameter uncertainty of 33% means area uncertainty of ___ (you do the math--it ain't good). Also, equipment used to measure and record diameters typically has finite resolution; for some older systems, diameters could only be recorded as multiples of ~ 0.5 um. Then there is the fact that models based on light microscopy typically assume that neurites have a circular cross section; anybody who has ever seen an EM image of brain tissue from any animal will know that this assumption is way off.

Of course, key biophysical parameters aren't much better. Estimates of specific membrane capacitance have about 2 place precision. Cytoplasmic resistivity is much harder to measure and uncertainty is on the order of 30% or worse.

Viewed in this light, the difference between the hoc- and Python-calculated "total surface areas" seems vanishingly small.