Axon metrics at zero for some files

[jonathanlurie]

Hello,

I am using NeuroM to compute some morph metrics on the soma, axon, basal dendrites and apical dendrite. In some cases, the axon metrics are all zeros even though there is an axon. From what we could gather, it seems that this issue is occuring in about 25% of cases among the ~200 morphs we processed. Note that things are working well for basal dendrite, apical dendrite (when present) and general morphology metrics.

Here is a short version of the code we use to compute the metrics:

METRIC_CONFIG = {
    'neurite': {
        'max_radial_distance': ['mean', 'std'],
        'number_of_sections': ['min', 'max', 'median', 'mean', 'std'],
        'number_of_bifurcations': ['min', 'max', 'median', 'mean', 'std'],
        'number_of_leaves': ['min', 'max', 'median', 'mean', 'std'],
        'total_length': ['min', 'max', 'median', 'mean', 'std'],
        'total_area': ['min', 'max', 'median', 'mean', 'std'],
        'total_volume': ['min', 'max', 'median', 'mean', 'std'],
        'section_lengths': ['min', 'max', 'median', 'mean', 'std'],
        'section_term_lengths': ['mean', 'std'],
        'section_bif_lengths': ['mean', 'std'],
        'section_branch_orders': ['min', 'max', 'mean', 'std'],
        'section_bif_branch_orders': ['mean', 'std'],
        'section_term_branch_orders': ['mean', 'std'],
        'section_path_distances': ['min', 'max', 'median', 'mean', 'std'],
        'section_taper_rates': ['min', 'max', 'median', 'mean', 'std'],
        'local_bifurcation_angles': ['min', 'max', 'median', 'mean', 'std'],
        'remote_bifurcation_angles': ['min', 'max', 'median', 'mean', 'std'],
        'partition_asymmetry': ['mean', 'std'],
        'partition_asymmetry_length': ['mean', 'std'],
        'sibling_ratios': ['min', 'max', 'median', 'mean', 'std'],
        'diameter_power_relations': ['min', 'max', 'median', 'mean', 'std'],
        'section_radial_distances': ['min', 'max', 'median', 'mean', 'std'],
        'section_term_radial_distances': ['mean', 'std'],
        'section_bif_radial_distances': ['mean', 'std'],
        'terminal_path_lengths': ['mean', 'std'],
        'section_volumes': ['min', 'max', 'median', 'mean', 'std'],
        'section_areas': ['min', 'max', 'median', 'mean', 'std'],
        'section_tortuosity': ['mean', 'std'],
        'section_strahler_orders': ['min', 'max', 'mean', 'std']
    },
    'morphology': {
        'soma_surface_area': ['min', 'max', 'median', 'mean', 'std'],
        'soma_radius': ['min', 'max', 'median', 'mean', 'std'],
        'max_radial_distance': ['mean', 'std'],
        'number_of_sections_per_neurite': ['min', 'max', 'median', 'mean', 'std'],
        'total_length_per_neurite': ['min', 'max', 'median', 'mean', 'std'],
        'total_area_per_neurite': ['min', 'max', 'median', 'mean', 'std'],
        'total_volume_per_neurite': ['min', 'max', 'median', 'mean', 'std'],
        'number_of_neurites': ['min', 'max', 'median', 'mean', 'std']
    },
    'neurite_type': ['AXON', 'BASAL_DENDRITE', 'APICAL_DENDRITE']
}

morph = nm.load_morphology(morphology_filepath)
metrics = morph_stats.extract_stats(morph, METRIC_CONFIG)

You can find below three examples of morphologies that are experiencing this issue. At first, this issue seemed related to the soma warning, though some morphologies to which the axon metrics are properly computed seem to also have a soma warning (so I just left it, just in case). I'll share the file provately so that you have a chance to rerun it.

• Morph 18454_00126.swc see it in the atlas

For this morph, NeuroM show the following message:

Warning: the soma does not conform the three point soma spec
The only valid neuro-morpho soma is:
1 1 x   y   z r -1
2 1 x (y-r) z r  1
3 1 x (y+r) z r  1

Got:
1 1 7397.3 3283.15 3618.82 6.13925 -1
2 1 7397.299805 3277.010742 3618.824951 6.139251 1
3 1 7397.299805 3289.289307 (exp. 3289.289062) 3618.824951 6.139251 1

• Morph 17302_00120.swc see it in the atlas

For this morph, NeuroM show the following message:

Warning: the soma does not conform the three point soma spec
The only valid neuro-morpho soma is:
1 1 x   y   z r -1
2 1 x (y-r) z r  1
3 1 x (y+r) z r  1

Got:
1 1 7397.3 3283.15 3618.82 6.13925 -1
2 1 7397.299805 3277.010742 3618.824951 6.139251 1
3 1 7397.299805 3289.289307 (exp. 3289.289062) 3618.824951 6.139251 1

• Morph 17302_00039.swc see it in the atlas

For this morph, NeuroM shows the following message:

Warning: found a disconnected neurite.
Neurites are not supposed to have parentId: -1
(although this is normal if this neuron has no soma)

[arnaudon]

Hello! You need to 'sanitize' your morphologies first with this: https://github.com/BlueBrain/NeuroR/blob/master/neuror/sanitize.py. BBP tools assume some minimal standard to deal with morphologies, and it seems (from the warning) that they don't meet that standard. This tool in NeuroR is fixing some issues and warning you about some. The point of our MCAR workflow will be to ensure all our morphologies are 'good enough. If sanitize does not do the trick, we can look closer.

My guess is: your neurite is disconnected, so it is not taken into account (likely to be the axon).

[arnaudon]

also, how do I dl your morphs? Can't find the button on the atlas UI

[jonathanlurie]

I can share in private

[mgeplf]

Please share w/ @eleftherioszisis too, thanks.

[jonathanlurie]

Yes, done.
After some investigation with @eleftherioszisis and @arnaudon , we found out that the morph causing the problem have one thing in common: their axon sprouts from a basal dendrite (close to the soma, like 3 nodes apart). From what I understand, NeuroM is not compatible with branching changing their type along the way, which leads to consider the whole branch as a basal dendrite (even though the cumulated length of this axon is more than 10x the length of the dendrite).

There are multiple ways to address this:

• make sure this forking is representing the truth (I'm suspected that this axon should instead be connected to the soma since it's only 2 or 3 nodes away)
• connect it straight to the soma instead of having a fork point (though this would need to be validated somehow)
• replicate the first part that happens between the soma and the axon-dendrite-fork to make one branch fully dendrite and the other fully axon

Not sure yet what is best but if a tool is developed to modify the structure, it could be part of NeuroR/sanitize.

Maybe @eleftherioszisis can give more info about the implications those workaround would have.

[jonathanlurie]

Another way to address this situation is to make NeuroM compatible with branching change of type, but this seems complex to introduce and there would probably be a lot of backward compatibility tests to perform.

[mgeplf]

make sure this forking is representing the truth (I'm suspected that this axon should instead be connected to the soma since it's only 2 or 3 nodes away)

For this axon sprout from basal-dendrite case, I wonder if we can make the error message better, so it's easier to understand the failure, and the feedback process of is this right? can be handled earlier.

NeuroM compatible with branching change of type

We'll have to look at this, but so far we've held off b/c (afaik) most examples of this have been an error - however, I'm not a neuroscientist, and I believe we know that heterogenous dendrite/axons exist, so it's something we should look into handling.

However, this would mean that these sorts of errors (if it is indeed an error) would not be uncovered, and I think they would cause problems down the toolchain.

[eleftherioszisis]

Reason & Example

When a morphology file is loaded, MorphIO does not know what a Neurite is; it only knows about root_sections.

NeuroM however supports neurites, the types of which are determined from the type of the morphio's root section:

NeuroM/neurom/core/morphology.py

Lines 346 to 349 in 5a48859

	@property
	def type(self):
	"""The type of the root node."""
	return self.root_node.type

This works well for most neurites, but it breaks down for neurites in which an axon starts from a basal dendrite. I have created a test cell with one tree and two different subtree types

import neurom
n = neurom.load_morphology(
"""
1 1 7397.3 3283.1 3618.8 1.0 -1
2 1 7397.3 3282.1 3618.8 1.0 1
3 1 7397.3 3284.1 3618.8 1.0 1
4 3 7397.2 3283.0 3615.9 2.2 1
5 3 7397.3 3282.3 3614.3 2.2 4
6 3 7398.1 3281.4 3612.8 2.2 5
7 3 7398.6 3280.5 3611.4 2.2 6
8 3 7398.8 3280.0 3609.6 2.2 7
9 3 7399.4 3279.5 3608.9 2.2 8
10 3 7399.5 3279.5 3608.8 2.2 9
11 2 7400.3 3278.9 3607.2 1.0 7
12 2 7400.2 3277.6 3605.7 1.0 11
13 2 7400.3 3277.6 3605.7 1.0 12
""",
reader="swc")

print(n.neurites)
print(n.sections)

which outputs:

Neurites:
Neurite <type: NeuriteType.basal_dendrite>

Sections:
Section(id=0, type=3, n_points=4)<parent: Section(id=None), nchildren: 2>,
Section(id=1, type=3, n_points=4)<parent: Section(id=0), nchildren: 0>,
Section(id=2, type=2, n_points=4)<parent: Section(id=0), nchildren: 0>

When a feature, number_of_sections for instance, is calculated, it takes into account the neurite types for filtering:

get("number_of_sections", n, neurite_type=neurom.BASAL_DENDRITE)

outputs:

3

And of course:

get("number_of_sections", n, neurite_type=neurom.AXON)

outputs:

0

Therefore, @jonathanlurie this is the expected behavior with NeuroM given the current assumptions.

Solution Proposals

1. Adapt NeuroM to handle correctly the calculation of features on mixed type trees

Change how the tree/section representations and filtering are handled by the morph/neurite/section iterators. Certain features would require handling special cases so that they are calculated correctly if the starting point is not the soma.

Pros:

• Features would be calculated correctly on all kinds of trees

Cons:

• Many parts of NeuroM would require to be rewritten/revised
• Will possibly break tools dependent on NeuroM

2. Connect axon to the soma and remove it from the basal
A new sanitization function will be added to NeuroR that would remove the axonal part on the basal dendrite and connect it directly to the soma as if it was a standalone neurite.

Pros:

• Separates the mixed types
• No special logic is required in NeuroM

Cons

• Some axonal features will be wrongly estimated because of the sortest path connection
• Some features will be overestimated, such as the total length of the cell due to having two paths to the soma instead of a shared one

3. Clone shared branches and generate two different trees with an overlapping initial part

A new sanitization function will be added to NeuroR that would duplicate the shared neurite between the basal and the axon from the soma to the start of the latter. Then instead of one tree, two would be created for the basal and axon parts respectively.

Pros:

• Separates the mixed types
• No special logic is required in NeuroM
• Axon and basal features would be accurately calculated

Cons:

• Some features will be overestimated, such as the total length of the cell due to counting twice the shared part of the trees

IMHO the third solution looks like the most promising one as we would not need to change/break NeuroM for taking into account these not-so-frequent types of trees.

Thoughts? Suggestions? Concerns?
@mgeplf @lidakanari @arnaudon @wvangeit

[arnaudon]

I don't see much difference between 2 and 3, apart from adding a straight line to connect to soma, or a curly line. It will have the same cons.

Imo, from the science point of view:

1. if we have a fake axon bearing dendrite (reconstruction error), we can use 2 to sanitize it.
2. Instead, if it is a real one, it should be marked as non-valid, and the reconstructor should manually decide to use/discard the axon or move it to connect to soma somehow. This is only to comply with current bbp implementation: such cases are not handled at all in any tools ( I believe).
3. then, if somebody really wants to do science with these, in terms of simulations, etc... it's another matter, but I guess in this case, NEURON/morphio could handle it fine.

[lidakanari]

Here are some comments on this...

1. This is a real effect, it cannot be considered as an error. In some rodent brain regions this can be as frequent as 70% of the data! So imo this should not be removed / sanitized.
2. The main question is what does this axon starting from a dendrite affect? Electrophysiology-wise I think this does not have any impact on the process. So this leaves morphological analysis and connectivity. I am not well equipped to judge the effect on connectivity so I will focus on morphology.
3. Morphological analysis: In order to be 100% correct we would need to process the axon as starting from the soma. This should impact measurements such as total lengths etc. If I understand this correctly, it is not currently possible to have neurites starting from anywhere in neuroM. So solution # 1 might be difficult to implement @eleftherioszisis? If that's the case I would go with solution # 3. However, this will introduce a known error, it should be documented and the user has to be aware of this issue!

[mgeplf]

Thanks @eleftherioszisis for outlining possible solutions.

Thanks @lidakanari; that's valuable input. I agree 100% with So imo this should not be removed / sanitized.

I lean heavily to option 1; Adapt NeuroM to handle correctly the calculation of features on mixed type trees; since this is a real phenomenon, I think it would be best if we covered it without work-arounds that might result in surprise.

I looked through the features:

• neurom/features/morphology.py: all neurite related functions have an arguement: neurite_type=NeuriteType.all; fortunately, they are all (?) implemented w/ the is_type and iter_neurites helpers; so we just need to agree on how to modify this, and make it clear to the user what the behavior is.

  One option would be to add{axon,apical_dendrite,basal_dendrite}_heterogeneous; these types would only look at the soma connected portion, and ignore if the neurite was heterogeneous. The old axon,(apical|basal)dendrite would only apply themselves to to homogeneous neurites.

  Notes:

  • the behavior could be reversed, if that makes more sense
  • MorphIO could add a helper to check if a subtree is heterogeneous
    to make the computation faster
  • the same applies to neurom/features/population.py

• neurom/features/neurite.py: everything in is agnostic of type; if it's heterogeneous, it should be fine

• neurom/features/section.py: this can be left the same, as long as it seems 'natural' that things like section_path_length can cross heterogeneous boundaries

I'm sure I've missed other impacts of the heterogeneous change, so if you can think of them, please list them.

With respect to this being a breaking change; we could jump to version 4 to signal that people need to be aware of this.

[arnaudon]

+1 and it could help here as well: BlueBrain/morph-tool#70

[eleftherioszisis]

imo the axon, (apical|basal)_dendrite types should apply to both homogeneous and heterogeneous neurites. This way statistics for all axons, basals, and apicals would be possible without the user having to know the further subtype distinctions. That of course would break backward consistency of results if mixed neurites are included.

axon = axon_inhomogeneous | axon_homogeneous

Are there any other cases of mixed subtypes, or is the axon starting from a basal dendrite the only known exception?

[arnaudon]

I know these sub_neurite_types: apical_oblique, apical_tuft, apical_trunk, axon_ais, axon_collaterals, axon_main

[eleftherioszisis]

How are these annotated in the morphology file? Is there a specification with their identifiers?

[arnaudon]

they are not, the PR above tries to assign them (apical can be split if one knows apical point for example). So if a proper subtype option is available, handling these would be less hacky, which could be convenient! (but this is not a priority at all!)

[eleftherioszisis]

@arnaudon, this discussion is mainly about existing section types and neurites that may have mixed combinations of existing section types. Therefore, we should discuss that use-case in a different thread.

I would like to also share the following ideas:

1. Introduce a general heterogenous type for neurites that have more than one section types.
2. If a filter encounters a heterogeneous type, instead of filtering it as a whole, it could delegate the filtering on a per section basis. Of course, this relaxation of type checking may have side effects if users provide trees without contiguous subtrees.
3. Features that require performing tree traversals (e.g. path_length) need to be extended so that they can stop on different nodes than root/terminations.

An alternative route, which is still vague in my mind would be to introduce a NeuriteType dataclass with subtypes.

@dataclass
class NeuriteType:
    subtypes: Optional[set] = None

heterogeneous = NeuriteType({SectionType.axon, SectionType.basal_dendrite})

[arnaudon]

Sure, but it does not make sense to have a super-neurite which contains an axon and a basal imo. An axon is an axon, regardless if it is attached to the soma, or any other sections of basal. I would rather relax the root_section to be something else than the soma, if possible.

[eleftherioszisis]

NeuroM is not designed for introducing mutations. Therefore, breaking trees and making intermediate containers that store subtrees is not the way to go imho, as we should stay as close as possible to the biological data.

Essentially, we are talking about the same thing with the only difference being that the distinction/processing is made while traversing the data as opposed to changing the data so that we can do the classical processing. If the latter was to be done I would suggest (as I did in the first set of suggestions), to do it before NeuroM and pass the curated data to NeuroM for classical analysis.

[arnaudon]

Funny enough, I just bumped into that issue now. So I looked at neurom a bit more, and found that in core/morphology.py

443     @property                                                                                       
444     def neurites(self):                                                                                                     
445         """The list of neurites."""                                                                                                              
446         return [Neurite(root_section) for root_section in self.root_sections]                                                             

where self.root_sections are the sections emerging from the soma red by morphio. I believe that if I simply append the root section of the axon manually (I can do a search for non-contiguous sections), neurom will behave fine (maybe one needs to stop iteration over sections if type changes or something).

[eleftherioszisis]

Creating a neurite with a different section than root_section does not mean that the section becomes a root_section. Features that traverse the tree upstream (path_length) will still stop at the initial root_section.

[arnaudon]

Use a self.is_root here to stop when you want?

 57     @property                                                                                       
 58     def parent(self):                                                                               
 59         """Returns the parent section if non root section else None."""                             
 60         if self.morphio_section.is_root:                                                            
 61             return None                                                                             
 62         return Section(self.morphio_section.parent)                                                 
 63              

[eleftherioszisis]

Not sure I follow the "Use a self.is_root here to stop when you want?".

[arnaudon]

instead of asking morphio if a section is root (which would be True only for root_sections), set an attribute is_root for neurom sections (which is == morphio_section.is_root for most of the cases), which we can set to the 'custom' root_section to stop the ipustream before you jump into another dendrite. This is the possible issue you were referring to, no?

[eleftherioszisis]

Ah, I see. I am afraid we cannot bypass the graph representation just to serve a specific use case. That's a hack and I am sure there is a better way to do this without breaking the neurom/morphio integration.

[arnaudon]

How would be break morphio/neurom? Morphio does not know about neurites, and this is only about neurite?

[eleftherioszisis]

The integration of morphio into neurom, not morphio itself. On the neurom side we would have a state that is not present on the morphio side. is_root is a section method, not a neurite one.

[arnaudon]

yes, I mean neuron Section, it has an is_root which is not doing much, we can add a .is_initial or something. Why is it an issue to add a new property in neurom section?

Anyway, I'm curious about this stuff now, I'll play with the code a little next week. Have a good week end!

[arnaudon]

Attempt here:#977

[eleftherioszisis]

A check for the heterogeneity of trees has been added to MorphIO (BlueBrain/MorphIO#360). This will allow NeuroM to easily check if a neurite is inhomogeneous so that it can process it in a special way.

Before going deeper into the implementation details, however, we need to first pin down the behavior of this feature from the user side:

1. Should the user have access to the inhomogeneous processing, using the get feature function?
2. Should it be optional because it will not be useful/preferable to everyone?