Automatic segmentation in multicompartmental models

[romainbrette]

When creating a multicompartmental model, a basic question is how many compartments we need to have. This issue also arises when loading files coming from reconstructions, which often have a huge number of compartments, and merging is desired.

One method is to split a branch evenly in units of space length. But the problem is the space length is calculated for DC inputs, and it depends on channel conductances. At high frequencies (transients), the space length is smaller, and if frequency is high enough, it doesn't depend on channel conductances anymore:

lambda(f) = 1/2 * (d/(pi*f*Ri*Cm))^(1/2)

This is what Neuron uses, with f = 100 Hz (which is of course arbitrary) and standard compartment length 0.1*lambda(f). What we notice is that in both cases (lambda(0) and lambda(f)), the space length scales as sqrt(d), ie, like axial resistance. Cm is almost always a constant (except for myelin), of around 1 µF/cm^2; Ri is normally in the order of 100 ohm*cm. Therefore, we can simply consider that, for the purpose of resegmentation, lambda(f) = k*sqrt(d), where k is a constant.

In summary, I propose that automatic resegmentation is based on splitting a branch into compartments of equal length, in units of sqrt(d). One could specify a total number of compartments, or a precision (dimensionless number). A cylinder could also be initialized without specifying the number of compartments, or by specifying a precision value.

There is one issue, which is to deal with truncated cones. Initially, it may be fine to simply consider the average diameter (it won't make much of a difference anyway). Otherwise, we can do a more specific calculation.

[rcaze]

We had a long discussion on segmenting with @mstimberg in issue #569 and the segmeter function is segmenting a morphological object

import unittest
import numpy as np
from itertools import izip
import brian2 as br2

def dist2p(pt1,  pt2):
    """Compute the Euclidean distance between two points"""
    return np.sqrt(np.sum([(a-b)**2 for a, b in izip(pt1, pt2)]))

def segmenter(sec, nseg=2):
    """Segment the coordinate of all segments given the coordinates of
    the section and the number of segments
    """
    # Give the distance between all the points of the section and onvert in m
    lengths = br2.asarray(sec.length)*1e6
    # Compute the total distance after each points
    distances = [np.sum(lengths[:i+1]) for i in range(len(lengths))]
    # Give the total length
    tot_length = distances[-1]
    # Give the distances after each segments
    seg_dists = [0.] + [(i*tot_length)/float(nseg)
                        for i in range(1, nseg+1)]

    # We use the distance to set the new coordinates in Brian
    # Set the new coordinates of the section
    seg_length = tot_length/float(nseg)

    # Add the points to interpolate
    total_dists = distances + seg_dists

    # Turn into arrays and uniquify or sort
    total_dists = np.unique(np.array(total_dists, np.float))

    # Use the interpolation to obtain the new coordinates
    new_x = np.interp(total_dists, distances, br2.asarray(sec.x)*1e6)
    new_y = np.interp(total_dists, distances, br2.asarray(sec.y)*1e6)
    new_z = np.interp(total_dists, distances, br2.asarray(sec.z)*1e6)

    segs = []
    for i in range(0, nseg):
        mask1 = total_dists >= seg_dists[i]
        mask2 = total_dists <= seg_dists[i+1]
        mask = mask1 * mask2
        segs.append(np.array([new_x[mask], new_y[mask], new_z[mask]]))

    return segs, seg_length

Now the function outputs n sets of coordinates defining (possibly non-straight) segments. The bit that I do not know is to turn these coordinates into meaningful new Morphology objects and put them together into one "section" object. Note that in this first draft of function there is no care about the diameters. It seems the first step toward segmenting.