"""
Neuronal morphology module.
This module defines classes to load and build neuronal morphologies.
"""
import abc
import numbers
import os
from abc import abstractmethod
from collections import OrderedDict, defaultdict, namedtuple
from brian2 import numpy as np
from brian2.units.allunits import meter
from brian2.units.fundamentalunits import Quantity, check_units, have_same_dimensions
from brian2.units.stdunits import um
from brian2.utils.logger import get_logger
logger = get_logger(__name__)
__all__ = ["Morphology", "Section", "Cylinder", "Soma"]
Node = namedtuple("Node", field_names="index,comp_name,x,y,z,diameter,parent,children")
def _to_meters(value):
"""
Helper function to convert a floating point value (or array) to a `Quantity`
in units of "meter", but also allow for ``None`` and return it as it is.
"""
if value is None:
return None
else:
return Quantity(value, dim=meter.dim)
def _from_morphology(variable, i, j):
"""
Helper function to return coordinates from a main morphology (used by
`SubMorphology`), dealing with ``None``.
"""
if variable is None:
return None
return variable[i:j]
[docs]
class MorphologyIndexWrapper:
"""
A simpler version of `~brian2.groups.group.IndexWrapper`, not allowing for
string indexing (`Morphology` is not a `Group`). It allows to use
``morphology.indices[...]`` instead of ``morphology[...]._indices()``.
"""
def __init__(self, morphology):
self.morphology = morphology
def __getitem__(self, item):
if isinstance(item, str):
raise NotImplementedError("Morphologies do not support string indexing")
assert isinstance(self.morphology, (SubMorphology, Morphology))
return self.morphology._indices(item)
def _calc_start_idx(section):
"""
Calculate the absolute start index that will be used by a flattened
representation.
"""
# calculate the absolute start index of this section
# 1. find the root of the tree
root = section
while root._parent is not None:
root = root._parent
# 2. go down from the root and advance the indices until we find
# the current section
start_idx, found = _find_start_index(root, section)
assert found
return start_idx
def _find_start_index(current, target_section, index=0):
if current == target_section:
return index, True
index += current.n
for child in current.children:
if child == target_section:
return index, True
else:
index, found = _find_start_index(child, target_section, index)
if found:
return index, True
return index, False
[docs]
class Topology:
"""
A representation of the topology of a `Morphology`. Has a useful string
representation, inspired by NEURON's ``topology`` function.
"""
def __init__(self, morphology):
self.morphology = morphology
def __str__(self):
# TODO: Make sure that the shown compartments do not get out of hand
divisor = 1
return Topology._str_topology(self.morphology, compartments_divisor=divisor)
@staticmethod
def _str_topology(
morphology, indent=0, named_path="", compartments_divisor=1, parent=None
):
"""
A simple string-based representation of a morphology. Inspired by
NEURON's ``topology`` function.
"""
description = " " * indent
length = max([1, morphology.n // compartments_divisor])
if parent is not None:
description += "`"
if isinstance(morphology, Soma):
description += "( )"
else:
description += "-" * length
description += "|"
if len(named_path) == 0:
description += " [root] \n"
else:
description += f" {named_path}\n"
for child in morphology.children:
name = morphology.children.name(child)
description += Topology._str_topology(
child,
indent=indent + 2 + length,
named_path=f"{named_path}.{name}",
compartments_divisor=compartments_divisor,
parent=morphology,
)
return description
__repr__ = __str__
def _rotate(vec, axis, angle):
"""
Rotate a vector around an arbitrary axis.
Parameters
----------
vec : `ndarray`
The vector to rotate.
axis : `ndarray`
The axis around which the vector should be rotated.
angle : float
The rotation angle (in radians).
Returns
-------
rotated : `ndarray`
The rotated vector.
"""
return (
vec * np.cos(angle)
- np.cross(axis, vec) * np.sin(angle)
+ axis * np.dot(axis, vec) * (1 - np.cos(angle))
)
def _perturb(vec, sigma):
if sigma == 0:
return vec
# Get an arbitrary orthogonal vector
if vec[1] != 0 or vec[0] != 0:
orthogonal = np.hstack([vec[1], vec[0], 0])
else: # special case for the [0, 0, 1] vector
orthogonal = np.array([1, 0, 0])
# Rotate the orthogonal vector
orthogonal = _rotate(orthogonal, vec, np.random.rand() * np.pi * 2)
# Use an exponentially distributed angle for the perturbation
perturbation = np.random.exponential(sigma, 1)
return _rotate(vec, orthogonal, perturbation)
def _add_coordinates(
orig_morphology,
root=None,
parent=None,
name=None,
section_randomness=0.0,
compartment_randomness=0.0,
n_th_child=0,
total_children=0,
overwrite_existing=False,
):
# Note that in the following, all values are without physical units
# The new direction is based on the direction of the parent section
if parent is None:
section_dir = np.array([0, 0, 0])
else:
section_dir = np.hstack(
[
np.asarray(parent.end_x[-1] - parent.start_x[0]),
np.asarray(parent.end_y[-1] - parent.start_y[0]),
np.asarray(parent.end_z[-1] - parent.start_z[0]),
]
)
parent_dir_norm = np.sqrt(np.sum(section_dir**2))
if parent_dir_norm != 0:
section_dir /= parent_dir_norm
else:
section_dir = np.array([0, 0, 0])
if not overwrite_existing and orig_morphology.x is not None:
section = orig_morphology.copy_section()
elif isinstance(orig_morphology, Soma):
# No perturbation for the soma
section = Soma(
diameter=orig_morphology.diameter,
x=section_dir[0] * meter,
y=section_dir[1] * meter,
z=section_dir[2] * meter,
)
else:
if np.sum(section_dir**2) == 0:
# We don't have any direction to base this section on (most common
# case is that the root section is a soma)
# We stay in the x-y plane and distribute all children in a 360 degree
# circle around (0, 0, 0)
section_dir = np.array([1, 0, 0])
rotation_axis = np.array([0, 0, 1])
angle_increment = 2 * np.pi / total_children
rotation_angle = np.pi / 2 + angle_increment * n_th_child
section_dir = _rotate(section_dir, rotation_axis, rotation_angle)
else:
if (
section_randomness == 0 and section_dir[2] == 0
): # If we are in the x-y plane, stay there
rotation_axis = np.array([0, 0, 1])
else:
rotation_axis = np.array([-section_dir[1], section_dir[2], 0])
if section_randomness == 0:
angle_increment = np.pi / (total_children + 1)
rotation_angle = -np.pi / 2 + angle_increment * (n_th_child + 1)
section_dir = _rotate(section_dir, rotation_axis, rotation_angle)
if section_randomness > 0:
# Rotate randomly
section_dir = _perturb(section_dir, section_randomness)
section_dir_norm = np.sqrt(np.sum(section_dir**2))
section_dir /= section_dir_norm
# For a soma, we let child sections begin at the surface of the sphere
if isinstance(parent, Soma):
origin = parent.diameter / 2 * section_dir
else:
origin = (0, 0, 0) * um
coordinates = np.zeros((orig_morphology.n + 1, 3)) * meter
start_coords = origin
coordinates[0, :] = origin
# Perturb individual compartments as well
for idx, length in enumerate(orig_morphology.length):
compartment_dir = _perturb(section_dir, compartment_randomness)
compartment_dir_norm = np.sqrt(np.sum(compartment_dir**2))
compartment_dir /= compartment_dir_norm
current_coords = start_coords + length * compartment_dir
coordinates[idx + 1, :] = current_coords
start_coords = current_coords
if isinstance(orig_morphology, Cylinder) and compartment_randomness == 0:
section = Cylinder(
n=orig_morphology.n,
diameter=orig_morphology.diameter[0],
x=coordinates[[0, -1], 0],
y=coordinates[[0, -1], 1],
z=coordinates[[0, -1], 2],
type=orig_morphology.type,
)
elif isinstance(orig_morphology, Section):
section = Section(
n=orig_morphology.n,
diameter=np.hstack(
[orig_morphology.start_diameter[0], orig_morphology.end_diameter]
)
* meter,
x=coordinates[:, 0],
y=coordinates[:, 1],
z=coordinates[:, 2],
type=orig_morphology.type,
)
else:
raise NotImplementedError(
f"Do not know how to deal with section of type {type(orig_morphology)}."
)
if parent is None:
root = section
else:
parent.children.add(name, section)
for idx, child in enumerate(orig_morphology.children):
_add_coordinates(
child,
root=root,
parent=section,
name=orig_morphology.children.name(child),
n_th_child=idx,
total_children=len(orig_morphology.children),
section_randomness=section_randomness,
compartment_randomness=compartment_randomness,
overwrite_existing=overwrite_existing,
)
return section
[docs]
class Children:
"""
Helper class to represent the children (sub trees) of a section. Can be
used like a dictionary (mapping names to `Morphology` objects), but iterates
over the values (sub trees) instead of over the keys (names).
"""
def __init__(self, owner):
self._owner = owner
self._counter = 0
self._children = []
self._named_children = {}
self._given_name = defaultdict(lambda: None)
def __iter__(self):
return iter(self._children)
def __len__(self):
return len(self._children)
def __contains__(self, item):
return item in self._named_children
[docs]
def name(self, child):
"""
Return the given name (i.e. not the automatic name such as ``1``) for a
child subtree.
Parameters
----------
child : `Morphology`
Returns
-------
name : str
The given name for the ``child``.
"""
return self._given_name[child]
def __getitem__(self, item):
if isinstance(item, str):
return self._named_children[item]
else:
raise TypeError("Index has to be an integer or a string.")
[docs]
def add(self, name, subtree, automatic_name=False):
"""
Add a new child to the morphology.
Parameters
----------
name : str
The name (e.g. ``"axon"``, ``"soma"``) to use for this sub tree.
subtree : `Morphology`
The subtree to link as a child.
automatic_name : bool, optional
Whether to chose a new name automatically, if a subtree of the same
name already exists (uses e.g. ``"dend2"`` instead ``"dend"``).
Defaults to ``False`` and will raise an error instead.
"""
if name in self._named_children and self._named_children[name] is not subtree:
if automatic_name:
basename = name
counter = 1
while name in self._named_children:
counter += 1
name = basename + str(counter)
else:
raise AttributeError(f"The name {name} is already used for a subtree.")
if subtree not in self._children:
self._counter += 1
self._children.append(subtree)
self._named_children[str(self._counter)] = subtree
self._given_name[subtree] = name
if name is not None:
self._named_children[name] = subtree
subtree._parent = self._owner
[docs]
def remove(self, name):
"""
Remove a subtree from this morphology.
Parameters
----------
name : str
The name of the sub tree to remove.
"""
if name not in self:
raise AttributeError(f"The subtree {name} does not exist")
subtree = self._named_children[name]
del self._named_children[name]
self._children.remove(subtree)
subtree._parent = None
def __repr__(self):
n = len(self._children)
s = f"<{int(n)} children"
if n > 0:
name_dict = {self.name(sec): sec for sec in self._children}
s += f": {name_dict!r}"
return f"{s}>"
[docs]
class Morphology(metaclass=abc.ABCMeta):
"""
Neuronal morphology (tree structure).
The data structure is a tree where each node is an un-branched section
consisting of a number of connected compartments, each one defined by its
geometrical properties (length, area, diameter, position).
Notes
-----
You cannot create objects of this class, create a `Soma`, a `Section`, or
a `Cylinder` instead.
"""
@check_units(n=1)
def __init__(self, n, type=None):
if isinstance(n, str):
raise TypeError(
"Need the number of compartments, not a string. "
"If you want to load a morphology from a file, "
"use 'Morphology.from_file' instead."
)
self._n = int(n)
if self._n != n:
raise TypeError("The number of compartments n has to be an integer value.")
if n <= 0:
raise ValueError("The number of compartments n has to be at least 1.")
self.type = type
self._children = Children(self)
self._parent = None
self.indices = MorphologyIndexWrapper(self)
def __getitem__(self, item):
"""
Return the subtree with the given name/index.
Ex.: ```neuron['axon']``` or ```neuron['11213']```
```neuron[10*um:20*um]``` returns the subbranch from 10 um to 20 um.
```neuron[10*um]``` returns one compartment.
```neuron[5]``` returns compartment number 5.
"""
if isinstance(item, slice): # neuron[10*um:20*um] or neuron[1:3]
using_lengths = all(
[
arg is None or have_same_dimensions(arg, meter)
for arg in [item.start, item.stop]
]
)
using_ints = all(
[
arg is None or int(arg) == float(arg)
for arg in [item.start, item.stop]
]
)
if not (using_lengths or using_ints):
raise TypeError("Index slice has to use lengths or integers")
if using_lengths:
if item.step is not None:
raise TypeError(
"Cannot provide a step argument when slicing with lengths"
)
pos = np.cumsum(np.asarray(self.length)) # coordinate on the section
# We use a special handling for values very close to the points
# between the compartments to avoid non-intuitive rounding
# effects: a point closer than 1e-12*length of section will be
# considered to be within the following section (for a start
# index), respectively within the previous section (for an end
# index)
if item.start is None:
i = 0
else:
diff = np.abs(float(item.start) - pos)
if min(diff) < 1e-12 * pos[-1]:
i = np.argmin(diff) + 1
else:
i = np.searchsorted(pos, item.start)
if item.stop is None:
j = len(pos)
else:
diff = np.abs(float(item.stop) - pos)
if min(diff) < 1e-12 * pos[-1]:
j = np.argmin(diff) + 1
else:
j = np.searchsorted(pos, item.stop) + 1
else: # integers
i, j, step = item.indices(self.n)
if step != 1:
raise TypeError("Can only slice a contiguous segment")
elif isinstance(item, Quantity) and have_same_dimensions(item, meter):
pos = np.hstack(
[0, np.cumsum(np.asarray(self.length))]
) # coordinate on the section
if float(item) < 0 or float(item) > (1 + 1e-12) * pos[-1]:
raise IndexError(
f"Invalid index {item}, has to be in the interval "
f"[{0 * meter!s}, {pos[-1] * meter!s}]."
)
diff = np.abs(float(item) - pos)
if min(diff) < 1e-12 * pos[-1]:
i = np.argmin(diff)
else:
i = np.searchsorted(pos, item) - 1
j = i + 1
elif isinstance(item, numbers.Integral): # int: returns one compartment
if item < 0: # allows e.g. to use -1 to get the last compartment
item += self.n
if item >= self.n:
raise IndexError(f"Invalid index {item} for {self.n} compartments")
i = item
j = i + 1
elif isinstance(item, str):
item = str(item) # convert int to string
if (len(item) > 1) and all(
[c in "LR123456789" for c in item]
): # binary string of the form LLLRLR or 1213 (or mixed)
return self._children[item[0]][item[1:]]
elif item in self._children:
return self._children[item]
else:
raise AttributeError(f"The subtree {item} does not exist")
else:
raise TypeError(f"Index of type {type(item)} not understood")
return SubMorphology(self, i, j)
def __setitem__(self, item, child):
"""
Inserts the subtree and name it ``item``.
Ex.: ``neuron['axon']`` or ``neuron['11213']``
"""
item = str(item) # convert int to string
if (len(item) > 1) and all([c in "LR123456789" for c in item]):
# binary string of the form LLLRLR or 1213 (or mixed)
self.children[item[0]][item[1:]] = child
else:
self.children.add(item, child)
def __delitem__(self, item):
"""
Remove the subtree ``item``.
"""
item = str(item) # convert int to string
if (len(item) > 1) and all([c in "LR123456789" for c in item]):
# binary string of the form LLLRLR or 1213 (or mixed)
del self._children[item[0]][item[1:]]
else:
self._children.remove(item)
def __getattr__(self, item):
"""
Return the subtree named ``item``.
Ex.: ``axon = neuron.axon``
"""
if item.startswith("_"):
return super(object, self).__getattr__(item)
else:
return self[item]
def __setattr__(self, item, child):
"""
Attach a subtree and name it ``item``.
Ex.: ``neuron.axon = Soma(diameter=10*um)``
"""
if isinstance(child, Morphology) and not item.startswith("_"):
self[item] = child
else: # If it is not a subtree, then it's a normal class attribute
object.__setattr__(self, item, child)
def __delattr__(self, item):
"""
Remove the subtree ``item``.
"""
del self[item]
def _indices(self, item=None, index_var="_idx"):
"""
Return compartment indices for the main section, relative to the
original morphology.
"""
if index_var != "_idx":
raise AssertionError(f"Unexpected index {index_var}")
if not (item is None or item == slice(None)):
if isinstance(item, slice):
# So that this always returns an array of values, even if it is
# just a single value
return self[item]._indices(slice(None))
else:
return self[item]._indices(None)
else:
start_idx = _calc_start_idx(self)
if self.n == 1 and item is None:
return start_idx
else:
return np.arange(start_idx, start_idx + self.n)
[docs]
def topology(self):
"""
Return a representation of the topology
Returns
-------
topology : `Topology`
An object representing the topology (can be converted to a string
by using ``str(...)`` or simply by printing it with `print`.)
"""
return Topology(self)
[docs]
def generate_coordinates(
self,
section_randomness=0.0,
compartment_randomness=0.0,
overwrite_existing=False,
):
r"""
Create a new `Morphology`, with coordinates filled in place where the
previous morphology did not have any. This is mostly useful for
plotting a morphology, it does not affect its electrical properties.
Parameters
----------
section_randomness : float, optional
The randomness when deciding the direction vector for each new
section. The given number is the :math:`\beta` parameter of an
exponential distribution (in degrees) which will be used to
determine the deviation from the direction of the parent section.
If the given value equals 0 (the default), then a deterministic
algorithm will be used instead.
compartment_randomness : float, optional
The randomness when deciding the direction vector for each
compartment within a section. The given number is the :math:`\beta`
parameter of an exponential distribution (in degrees) which will be
used to determine the deviation from the main direction of the
current section. If the given value equals 0 (the default), then all
compartments will be along a straight line.
overwrite_existing : bool, optional
Whether to overwrite existing coordinates in the morphology. This
is by default set to ``False``, meaning that only sections that do
not currently have any coordinates set will get new coordinates.
This allows to conveniently generate a morphology that can be
plotted for a morphology that is based on points but also has
artificially added sections (the most common case: an axon added
to a reconstructed morphology). If set to ``True``, all sections
will get new coordinates. This can be useful to either get a
schematic representation of the morphology (with
``section_randomness`` and ``compartment_randomness`` both 0) or to
simply generate a new random variation of a morphology (which will
still be electrically equivalent, of course).
Returns
-------
morpho_with_coordinates : `Morphology`
The same morphology, but with coordinates
"""
# Convert to radians
section_randomness *= np.pi / 180
compartment_randomness *= np.pi / 180
return _add_coordinates(
self,
section_randomness=section_randomness,
compartment_randomness=compartment_randomness,
overwrite_existing=overwrite_existing,
)
[docs]
@abstractmethod
def copy_section(self):
"""
Create a copy of the current section (attributes of this section only,
not re-creating the parent/children relation)
Returns
-------
copy : `Morphology`
A copy of this section (without the links to the parent/children)
"""
raise NotImplementedError()
@property
def n(self):
"""
The number of compartments in this section.
"""
return self._n
def __len__(self):
"""
This is not well-defined, use `Morphology.n` or
`Morphology.total_compartments` instead.
"""
raise TypeError(
"The 'length' of a Morphology is ambiguous, use its "
"'n' attribute for the number of compartments in this "
"section or the 'total_compartments' attribute for the "
"total number of compartments in the whole sub-tree."
)
@property
def total_compartments(self):
"""
The total number of compartments in this subtree (i.e. the number of
compartments in this section plus all the compartments in the sections
deeper in the tree).
"""
return self.n + sum(c.total_compartments for c in self.children)
@property
def total_sections(self):
"""
The total number of sections in this subtree.
"""
return 1 + sum(c.total_sections for c in self.children)
@property
def parent(self):
"""
The parent section of this section.
"""
return self._parent
@property
def children(self):
"""
The children (as a `Children` object) of this section.
"""
return self._children
@property
@abc.abstractmethod
def end_distance(self):
"""
The distance to the root of the morphology at the end of this section.
"""
raise NotImplementedError()
# Per-compartment attributes
@property
@abc.abstractmethod
def area(self):
"""
The membrane surface area of each compartment in this section.
"""
raise NotImplementedError()
@property
@abc.abstractmethod
def volume(self):
"""
The volume of each compartment in this section.
"""
raise NotImplementedError()
@property
@abc.abstractmethod
def length(self):
"""
The length of each compartment in this section.
"""
raise NotImplementedError()
@property
@abc.abstractmethod
def r_length_1(self):
"""
The geometry-dependent term to calculate the conductance between the
start and the midpoint of each compartment. Dividing this value by the
Intracellular resistivity gives the conductance.
"""
raise NotImplementedError()
@property
@abc.abstractmethod
def r_length_2(self):
"""
The geometry-dependent term to calculate the conductance between the
midpoint and the end of each compartment. Dividing this value by the
Intracellular resistivity gives the conductance.
"""
raise NotImplementedError()
# At-midpoint attributes
@property
@abc.abstractmethod
def diameter(self):
"""
The diameter at the middle of each compartment in this section.
"""
raise NotImplementedError()
@property
@abc.abstractmethod
def distance(self):
"""
The total distance between the midpoint of each compartment and the root
of the morphology.
"""
raise NotImplementedError()
@property
def start_x(self):
"""
The x coordinate at the beginning of each compartment. Returns ``None``
for morphologies without coordinates.
"""
return _to_meters(self.start_x_)
@property
def start_y(self):
"""
The y coordinate at the beginning of each compartment. Returns ``None``
for morphologies without coordinates.
"""
return _to_meters(self.start_y_)
@property
def start_z(self):
"""
The z coordinate at the beginning of each compartment. Returns ``None``
for morphologies without coordinates.
"""
return _to_meters(self.start_z_)
@property
@abc.abstractmethod
def start_x_(self):
"""
The x coordinate (as a unitless floating point number) at the beginning
of each compartment. Returns ``None`` for morphologies without
coordinates.
"""
raise NotImplementedError()
@property
@abc.abstractmethod
def start_y_(self):
"""
The y coordinate (as a unitless floating point number) at the beginning
of each compartment. Returns ``None`` for morphologies without
coordinates.
"""
raise NotImplementedError()
@property
@abc.abstractmethod
def start_z_(self):
"""
The z coordinate (as a unitless floating point number) at the beginning
of each compartment. Returns ``None`` for morphologies without
coordinates.
"""
raise NotImplementedError()
@property
def x(self):
"""
The x coordinate at the midpoint of each compartment. Returns ``None``
for morphologies without coordinates.
"""
return _to_meters(self.x_)
@property
def y(self):
"""
The y coordinate at the midpoint of each compartment. Returns ``None``
for morphologies without coordinates.
"""
return _to_meters(self.y_)
@property
def z(self):
"""
The y coordinate at the midpoint of each compartment. Returns ``None``
for morphologies without coordinates.
"""
return _to_meters(self.z_)
@property
@abc.abstractmethod
def x_(self):
"""
The x coordinate (as a unitless floating point number) at the midpoint
of each compartment. Returns ``None`` for morphologies without
coordinates.
"""
raise NotImplementedError()
@property
@abc.abstractmethod
def y_(self):
"""
The y coordinate (as a unitless floating point number) at the midpoint
of each compartment. Returns ``None`` for morphologies without
coordinates.
"""
raise NotImplementedError()
@property
@abc.abstractmethod
def z_(self):
"""
The z coordinate (as a unitless floating point number) at the midpoint
of each compartment. Returns ``None`` for morphologies without
coordinates.
"""
raise NotImplementedError()
@property
def end_x(self):
"""
The x coordinate at the end of each compartment. Returns ``None``
for morphologies without coordinates.
"""
return _to_meters(self.end_x_)
@property
def end_y(self):
"""
The y coordinate at the end of each compartment. Returns ``None``
for morphologies without coordinates.
"""
return _to_meters(self.end_y_)
@property
def end_z(self):
"""
The z coordinate at the end of each compartment. Returns ``None``
for morphologies without coordinates.
"""
return _to_meters(self.end_z_)
@property
@abc.abstractmethod
def end_x_(self):
"""
The x coordinate (as a unitless floating point number) at the end of
each compartment. Returns ``None`` for morphologies without coordinates.
"""
raise NotImplementedError()
@property
@abc.abstractmethod
def end_y_(self):
"""
The y coordinate (as a unitless floating point number) at the end of
each compartment. Returns ``None`` for morphologies without coordinates.
"""
raise NotImplementedError()
@property
@abc.abstractmethod
def end_z_(self):
"""
The z coordinate (as a unitless floating point number) at the end of
each compartment. Returns ``None`` for morphologies without coordinates.
"""
raise NotImplementedError()
@property
def coordinates(self):
r"""
Array with all coordinates at the start- and end-points of each
compartment in this section. The array has size :math:`(n+1) \times 3`,
where :math:`n` is the number of compartments in this section. Each
row is one point (start point of first compartment, end point of first
compartment, end point of second compartment, ...), with the columns
being the x, y, and z coordinates. Returns ``None`` for morphologies
without coordinates.
"""
if self.x_ is None:
return None
else:
return Quantity(self.coordinates_, dim=meter.dim)
@property
def coordinates_(self):
r"""
Array with all coordinates (as unitless floating point numbers) at the
start- and end-points of each compartment in this section. The array has
size :math:`(n+1) \times 3`, where :math:`n` is the number of
compartments in this section. Each row is one point (start point of
first compartment, end point of first compartment, end point of second
compartment, ...), with the columns being the x, y, and z coordinates.
Returns ``None`` for morphologies without coordinates.
"""
if self.x_ is None:
return None
else:
return np.vstack(
[
np.hstack([self.start_x_[0], self.end_x_[:]]),
np.hstack([self.start_y_[0], self.end_y_[:]]),
np.hstack([self.start_z_[0], self.end_z_[:]]),
]
).T
@staticmethod
def _create_section(compartments, name, parent, sections, spherical_soma):
if (
spherical_soma
and len(compartments) == 1
and compartments[0].comp_name == "soma"
):
soma = compartments[0]
section = Soma(
diameter=soma.diameter * um, x=soma.x * um, y=soma.y * um, z=soma.z * um
)
else:
sec_x, sec_y, sec_z, sec_diameter = zip(
*[(c.x, c.y, c.z, c.diameter) for c in compartments]
)
# Add a point for the end of the parent_idx compartment
if parent is not None:
n = len(compartments)
if parent.comp_name is not None and parent.comp_name.lower() == "soma":
# For a Soma, we don't use its diameter
start_diameter = sec_diameter[0]
else:
start_diameter = parent.diameter
# Use relative coordinates
sec_x = np.array(sec_x) - parent.x
sec_y = np.array(sec_y) - parent.y
sec_z = np.array(sec_z) - parent.z
start_x = start_y = start_z = 0.0
else:
n = len(compartments) - 1
start_diameter = sec_diameter[0]
sec_diameter = sec_diameter[1:]
start_x = sec_x[0]
start_y = sec_y[0]
start_z = sec_z[0]
sec_x = sec_x[1:]
sec_y = sec_y[1:]
sec_z = sec_z[1:]
diameter = np.hstack([start_diameter, sec_diameter]) * um
x = np.hstack([start_x, sec_x]) * um
y = np.hstack([start_y, sec_y]) * um
z = np.hstack([start_z, sec_z]) * um
section = Section(n=n, diameter=diameter, x=x, y=y, z=z, type=name)
# Add the section as a child to its parent
if parent is not None:
parent_sec = sections[parent.index]
parent_sec.children.add(name, section, automatic_name=True)
return section
@staticmethod
def _compartments_to_sections(
compartment, spherical_soma, current_compartments=None, sections=None
):
# Merge all unbranched compartments of the same type into a single
# section
if sections is None:
sections = OrderedDict()
if current_compartments is None:
current_compartments = []
# We have to create a new section, if we are either
# 1. at a leaf of the tree or at a branching point, or
# 2. if the compartment type changes
while True:
current_compartments.append(compartment)
if (
len(compartment.children) != 1
or compartment.children[0].comp_name != compartment.comp_name
):
break
compartment = compartment.children[0]
parent = current_compartments[0].parent
section = Morphology._create_section(
current_compartments,
compartment.comp_name,
parent=parent,
sections=sections,
spherical_soma=spherical_soma,
)
sections[current_compartments[-1].index] = section
# If we are at a branching point, recurse into all subtrees
for child in compartment.children:
Morphology._compartments_to_sections(
child,
spherical_soma=spherical_soma,
current_compartments=None,
sections=sections,
)
return sections
@staticmethod
def _replace_three_point_soma(compartment, all_compartments):
# Replace a three-point/two-cylinder soma by a single spherical soma
# if possible (see http://neuromorpho.org/SomaFormat.html for some
# details)
# We are looking for a node with two children of the soma type (and
# other childen of other types), where the two children don't have any
# children of their own
# To avoid recursion errors we are not searching any further
soma_children = [c for c in compartment.children if c.comp_name == "soma"]
if (
compartment.comp_name == "soma"
and len(soma_children) == 2
and all(len(c.children) == 0 for c in soma_children)
):
# This is a 3-point soma to replace
soma_c = [compartment] + soma_children
if not all(abs(c.diameter - soma_c[0].diameter) < 1e-15 for c in soma_c):
indices = ", ".join(str(c.index) for c in soma_c)
raise ValueError(
f"Found a '3-point-soma' (lines: {indices}), but not "
"all the diameters are "
"identical."
)
diameter = soma_c[0].diameter
point_0 = np.array([soma_c[0].x, soma_c[0].y, soma_c[0].z])
point_1 = np.array([soma_c[1].x, soma_c[1].y, soma_c[1].z])
point_2 = np.array([soma_c[2].x, soma_c[2].y, soma_c[2].z])
length_1 = np.sqrt(np.sum((point_1 - point_0) ** 2))
length_2 = np.sqrt(np.sum((point_2 - point_0) ** 2))
if (
np.abs(length_1 - diameter / 2) > 0.1
or np.abs(length_2 - diameter / 2) > 0.1
):
raise ValueError(
"Cannot replace '3-point-soma' by a single "
"point, the second and third points should "
"be positioned one radius away from the "
f"first point. Distances are {length_1:.3f}um and "
f"{length_2:.3f}um, respectively, while the "
f"radius is {diameter / 2:.3f}um."
)
children = [c for c in compartment.children if c not in soma_c]
compartment = Node(
index=compartment.index,
comp_name="soma",
x=point_0[0],
y=point_0[1],
z=point_0[2],
diameter=diameter,
parent=compartment.parent,
children=children,
)
all_compartments[compartment.index] = compartment
del all_compartments[soma_children[0].index]
del all_compartments[soma_children[1].index]
else:
# Raise a warning if there is more than one soma compartment
somata = [c for c in all_compartments.values() if c.comp_name == "soma"]
if len(somata) > 1:
logger.warning(
f"Found {len(somata)} soma compartments. If you have a 3-point soma, "
"make sure it is at the beginning.",
name_suffix="soma_compartments",
)
[docs]
@staticmethod
def from_points(points, spherical_soma=True):
"""
Create a morphology from a sequence of points (similar to the ``SWC``
format, see `Morphology.from_swc_file`). Each point has to be
a 7-tuple: ``(index, name, x, y, z, diameter, parent)``
Note that the values should not use units, but are instead all taken
to be in micrometers.
Parameters
----------
points : sequence of 7-tuples
The points of the morphology.
spherical_soma : bool, optional
Whether to model a soma as a sphere.
Returns
-------
morphology : `Morphology`
Notes
-----
This format closely follows the SWC format (see
`Morphology.from_swc_file`) with two differences: the ``type`` should
be a string (e.g. ``'soma'``) instead of an integer and the 6-th element
should be the diameter and not the radius.
"""
# First pass through all points to get the dependency structure
compartments = OrderedDict()
for counter, point in enumerate(points):
if len(point) != 7:
raise ValueError(
"Each point needs to be described by 7 values, got"
f" {len(point)} instead."
)
index, name, x, y, z, diameter, parent_idx = point
if index in compartments:
raise ValueError(f"Two compartments with index {int(index)}")
if parent_idx == index:
raise ValueError(
f"Compartment {int(index)} lists itself as the parent compartment."
)
if counter == 0 and parent_idx == -1:
parent = None # The first compartment does not have a parent
elif parent_idx not in compartments:
raise ValueError(
f"Did not find the compartment {parent_idx} (parent "
f"compartment of compartment {index}). Make sure "
"that parent compartments are listed before "
"their children."
)
else:
parent = compartments[parent_idx]
children = []
node = Node(index, name, x, y, z, diameter, parent, children)
compartments[index] = node
if parent is not None:
parent.children.append(node)
if spherical_soma:
Morphology._replace_three_point_soma(
list(compartments.values())[0], compartments
)
sections = Morphology._compartments_to_sections(
list(compartments.values())[0], spherical_soma
)
# Go through all the sections again and add standard names for all
# sections (potentially in addition to the name they already have):
# "L" + "R" for one or two children, "1", "2", "3", etc. otherwise
children_counter = defaultdict(int)
for section in sections.values():
parent = section.parent
if parent is not None:
children_counter[parent] += 1
children = parent.children
nth_child = children_counter[parent]
if len(children) <= 2:
name = "L" if nth_child == 1 else "R"
else:
name = f"{int(nth_child)}"
children.add(name, section)
# There should only be one section without parents
root = [sec for sec in sections.values() if sec.parent is None]
assert len(root) == 1
return root[0]
[docs]
@staticmethod
def from_swc_file(filename, spherical_soma=True):
"""
Load a morphology from a ``SWC`` file. A large database of morphologies
in this format can be found at http://neuromorpho.org
The format consists of an optional header of lines starting with ``#``
(ignored), followed by a sequence of points, each described in a line
following the format::
index type x y z radius parent
``index`` is an integer label (starting at 1) that identifies the
current point and increases by one each line. ``type`` is an integer
representing the type of the neural segment. The only type that changes
the interpretation by Brian is the type ``1`` which signals a soma.
Types ``2`` (axon), ``3`` (dendrite), and ``4`` (apical dendrite) are
used to give corresponding names to the respective sections. All other
types are ignored. ``x``, ``y``, and ``z`` are the cartesian coordinates
at each point and ``r`` is its radius. ``parent`` refers to the index
of the parent point or is ``-1`` for the root point.
Parameters
----------
filename : str
The name of the ``SWC`` file.
spherical_soma : bool, optional
Whether to model the soma as a sphere.
Returns
-------
morpho : `Morphology`
The morphology stored in the given file.
"""
swc_types = defaultdict(lambda: None)
# The following names will be translated into names, all other will be
# ignored
swc_types.update({"1": "soma", "2": "axon", "3": "dend", "4": "apic"})
with open(filename) as f:
points = []
for line_no, line in enumerate(f):
line = line.strip()
if line.startswith("#") or len(line) == 0:
# Ignore comments or empty lines
continue
splitted = line.split()
if len(splitted) != 7:
raise ValueError(
"Each line of an SWC file has to contain "
"7 space-separated entries, but line "
f"{line_no + 1} contains {len(splitted)}."
)
index, comp_type, x, y, z, radius, parent = splitted
points.append(
(
int(index),
swc_types[comp_type],
float(x),
float(y),
float(z),
2 * float(radius),
int(parent),
)
)
return Morphology.from_points(points, spherical_soma=spherical_soma)
[docs]
@staticmethod
def from_file(filename, spherical_soma=True):
"""
Convencience method to load a morphology from a given file. At the
moment, only ``SWC`` files are supported, calling this function is
therefore equivalent to calling `Morphology.from_swc_file` directly.
Parameters
----------
filename : str
The name of a file storing a morphology.
spherical_soma : bool, optional
Whether to model the soma as a sphere.
Returns
-------
morphology : `Morphology`
The morphology stored in the given file.
"""
_, ext = os.path.splitext(filename)
if ext.lower() == ".swc":
return Morphology.from_swc_file(filename, spherical_soma=spherical_soma)
else:
raise NotImplementedError(
"Currently, SWC is the only supported file format."
)
[docs]
class SubMorphology:
"""
A view on a subset of a section in a morphology.
"""
def __init__(self, morphology, i, j):
self._morphology = morphology
self.indices = MorphologyIndexWrapper(self)
self._i = i
self._j = j
def _indices(self, item=None):
if not (item is None or item == slice(None)):
raise IndexError("Cannot index a view on a subset of a section further")
# Start index of the main section
start_idx = _calc_start_idx(self._morphology)
if item is None and self.n == 1:
return start_idx + self._i
else:
return np.arange(start_idx + self._i, start_idx + self._j)
@property
def n(self):
"""
The number of compartments in this sub-section.
"""
return self._j - self._i
def __len__(self):
return self.n
@property
def n_sections(self):
"""
The number of sections in this sub-section (always 1).
"""
return 1
# Per-compartment attributes
@property
def area(self):
"""
The membrane surface area of each compartment in this sub-section.
"""
return self._morphology.area[self._i : self._j]
@property
def volume(self):
"""
The volume of each compartment in this sub-section.
"""
return self._morphology.volume[self._i : self._j]
@property
def length(self):
"""
The length of each compartment in this sub-section.
"""
return self._morphology.length[self._i : self._j]
@property
def r_length_1(self):
"""
The geometry-dependent term to calculate the conductance between the
start and the midpoint of each compartment in this sub-section.
Dividing this value by the Intracellular resistivity gives the
conductance.
"""
return self._morphology.r_length_1[self._i : self._j]
@property
def r_length_2(self):
"""
The geometry-dependent term to calculate the conductance between the
midpoint and the end of each compartment in this sub-section. Dividing
this value by the Intracellular resistivity gives the conductance.
"""
return self._morphology.r_length_2[self._i : self._j]
# At-midpoint attributes
@property
def diameter(self):
"""
The diameter at the middle of each compartment in this sub-section.
"""
return self._morphology.diameter[self._i : self._j]
@property
def distance(self):
"""
The total distance between the midpoint of each compartment in this
sub-section and the root of the morphology.
"""
return self._morphology.distance[self._i : self._j]
@property
def start_x(self):
"""
The x coordinate at the beginning of each compartment in this
sub-section. Returns ``None`` for morphologies without coordinates.
"""
return _to_meters(self.start_x_)
@property
def start_y(self):
"""
The y coordinate at the beginning of each compartment in this
sub-section. Returns ``None`` for morphologies without coordinates.
"""
return _to_meters(self.start_y_)
@property
def start_z(self):
"""
The x coordinate at the beginning of each compartment in this
sub-section. Returns ``None`` for morphologies without coordinates.
"""
return _to_meters(self.start_z_)
@property
def start_x_(self):
"""
The x coordinate (as a unitless floating point number) at the beginning
of each compartment in this sub-section. Returns ``None`` for
morphologies without coordinates.
"""
return _from_morphology(self._morphology.start_x_, self._i, self._j)
@property
def start_y_(self):
"""
The y coordinate (as a unitless floating point number) at the beginning
of each compartment in this sub-section. Returns ``None`` for
morphologies without coordinates.
"""
return _from_morphology(self._morphology.start_y_, self._i, self._j)
@property
def start_z_(self):
"""
The z coordinate (as a unitless floating point number) at the beginning
of each compartment in this sub-section. Returns ``None`` for
morphologies without coordinates.
"""
return _from_morphology(self._morphology.start_z_, self._i, self._j)
@property
def x(self):
"""
The x coordinate at the midpoint of each compartment in this
sub-section. Returns ``None`` for morphologies without coordinates.
"""
return _to_meters(self.x_)
@property
def y(self):
"""
The y coordinate at the midpoint of each compartment in this
sub-section. Returns ``None`` for morphologies without coordinates.
"""
return _to_meters(self.y_)
@property
def z(self):
"""
The z coordinate at the midpoint of each compartment in this
sub-section. Returns ``None`` for morphologies without coordinates.
"""
return _to_meters(self.z_)
@property
def x_(self):
"""
The x coordinate (as a unitless floating point number) at the midpoint
of each compartment in this sub-section. Returns ``None`` for
morphologies without coordinates.
"""
return _from_morphology(self._morphology.x_, self._i, self._j)
@property
def y_(self):
"""
The y coordinate (as a unitless floating point number) at the midpoint
of each compartment in this sub-section. Returns ``None`` for
morphologies without coordinates.
"""
return _from_morphology(self._morphology.y_, self._i, self._j)
@property
def z_(self):
"""
The z coordinate (as a unitless floating point number) at the midpoint
of each compartment in this sub-section. Returns ``None`` for
morphologies without coordinates.
"""
return _from_morphology(self._morphology.z_, self._i, self._j)
@property
def end_x(self):
"""
The x coordinate at the end of each compartment in this sub-section.
Returns ``None`` for morphologies without coordinates.
"""
return _to_meters(self.end_x_)
@property
def end_y(self):
"""
The y coordinate at the end of each compartment in this sub-section.
Returns ``None`` for morphologies without coordinates.
"""
return _to_meters(self.end_y_)
@property
def end_z(self):
"""
The z coordinate at the end of each compartment in this sub-section.
Returns ``None`` for morphologies without coordinates.
"""
return _to_meters(self.end_z_)
@property
def end_x_(self):
"""
The x coordinate (as a unitless floating point number) at the end of
each compartment in this sub-section. Returns ``None`` for morphologies
without coordinates.
"""
return _from_morphology(self._morphology.end_x_, self._i, self._j)
@property
def end_y_(self):
"""
The y coordinate (as a unitless floating point number) at the end of
each compartment in this sub-section. Returns ``None`` for morphologies
without coordinates.
"""
return _from_morphology(self._morphology.end_y_, self._i, self._j)
@property
def end_z_(self):
"""
The z coordinate (as a unitless floating point number) at the end of
each compartment in this sub-section. Returns ``None`` for morphologies
without coordinates.
"""
return _from_morphology(self._morphology.end_z_, self._i, self._j)
[docs]
class Soma(Morphology):
"""
A spherical, iso-potential soma.
Parameters
----------
diameter : `Quantity`
Diameter of the sphere.
x : `Quantity`, optional
The x coordinate of the position of the soma.
y : `Quantity`, optional
The y coordinate of the position of the soma.
z : `Quantity`, optional
The z coordinate of the position of the soma.
type : str, optional
The ``type`` of this section, defaults to ``'soma'``.
"""
@check_units(diameter=meter, x=meter, y=meter, z=meter)
def __init__(self, diameter, x=None, y=None, z=None, type="soma"):
Morphology.__init__(self, n=1, type=type)
if diameter.shape != () and len(diameter) != 1:
raise TypeError("Diameter has to be a scalar value.")
for coord in [x, y, z]:
if coord is not None and coord.shape != () and len(coord) != 1:
raise TypeError("Coordinates have to be scalar values.")
self._diameter = np.ones(1) * diameter
if any(coord is not None for coord in (x, y, z)):
default_value = np.array([0.0])
else:
default_value = None
self._x = np.atleast_1d(np.asarray(x)) if x is not None else default_value
self._y = np.atleast_1d(np.asarray(y)) if y is not None else default_value
self._z = np.atleast_1d(np.asarray(z)) if z is not None else default_value
def __repr__(self):
s = f"{self.__class__.__name__}(diameter={self.diameter[0]!r}"
if self._x is not None:
s += f", x={self.x[0]!r}, y={self.y[0]!r}, z={self.z[0]!r}"
if self.type != "soma":
s += f", type={self.type!r}"
return f"{s})"
[docs]
def copy_section(self):
return Soma(self.diameter, x=self.x, y=self.y, z=self.z, type=self.type)
# Note that the per-compartment properties should always return 1D arrays,
# i.e. for the soma arrays of length 1 instead of scalar values
@property
def area(self):
"""
The membrane surface area of this section (as an array of length 1).
"""
return np.pi * self.diameter**2
@property
def volume(self):
"""
The volume of this section (as an array of length 1).
"""
return (np.pi * self.diameter**3) / 6
@property
def length(self):
"""
The "length" (equal to `diameter`) of this section (as an array of
length 1).
"""
return self.diameter
@property
def r_length_1(self):
"""
The geometry-dependent term to calculate the conductance between the
start and the midpoint of each compartment. Returns a fixed (high)
value for a `Soma`, corresponding to a section with very low
intracellular resistance.
"""
return [1] * meter
@property
def r_length_2(self):
"""
The geometry-dependent term to calculate the conductance between the
midpoint and the end of each compartment. Returns a fixed (high)
value for a `Soma`, corresponding to a section with very low
intracellular resistance.
"""
return [1] * meter
@property
def diameter(self):
"""
The diameter of this section (as an array of length 1).
"""
return self._diameter
@property
def distance(self):
"""
The total distance between the midpoint of this section and the root
of the morphology. The `Soma` is most likely the root of the
morphology, and therefore the `distance` is 0.
"""
dist = self._parent.distance[-1:] if self._parent is not None else [0] * um
return dist
@property
def start_x_(self):
"""
The x-coordinate of the current section (as an array of length 1). Note
that a `Soma` is modelled as a "point" with finite surface/volume,
equivalent to that of a sphere with the given `diameter`. It's start-,
midpoint-, and end-coordinates are therefore identical.
"""
return self._x
@property
def start_y_(self):
"""
The y-coordinate of the current section (as an array of length 1). Note
that a `Soma` is modelled as a "point" with finite surface/volume,
equivalent to that of a sphere with the given `diameter`. It's start-,
midpoint-, and end-coordinates are therefore identical.
"""
return self._y
@property
def start_z_(self):
"""
The z-coordinate of the current section (as an array of length 1). Note
that a `Soma` is modelled as a "point" with finite surface/volume,
equivalent to that of a sphere with the given `diameter`. It's start-,
midpoint-, and end-coordinates are therefore identical.
"""
return self._z
@property
def x_(self):
"""
The x-coordinate of the current section (as an array of length 1). Note
that a `Soma` is modelled as a "point" with finite surface/volume,
equivalent to that of a sphere with the given `diameter`. It's start-,
midpoint-, and end-coordinates are therefore identical.
"""
return self._x
@property
def y_(self):
"""
The y-coordinate of the current section (as an array of length 1). Note
that a `Soma` is modelled as a "point" with finite surface/volume,
equivalent to that of a sphere with the given `diameter`. It's start-,
midpoint-, and end-coordinates are therefore identical.
"""
return self._y
@property
def z_(self):
"""
The z-coordinate of the current section (as an array of length 1). Note
that a `Soma` is modelled as a "point" with finite surface/volume,
equivalent to that of a sphere with the given `diameter`. It's start-,
midpoint-, and end-coordinates are therefore identical.
"""
return self._z
@property
def end_x_(self):
"""
The x-coordinate of the current section (as an array of length 1). Note
that a `Soma` is modelled as a "point" with finite surface/volume,
equivalent to that of a sphere with the given `diameter`. It's start-,
midpoint-, and end-coordinates are therefore identical.
"""
return self._x
@property
def end_y_(self):
"""
The y-coordinate of the current section (as an array of length 1). Note
that a `Soma` is modelled as a "point" with finite surface/volume,
equivalent to that of a sphere with the given `diameter`. It's start-,
midpoint-, and end-coordinates are therefore identical.
"""
return self._y
@property
def end_z_(self):
"""
The z-coordinate of the current section (as an array of length 1). Note
that a `Soma` is modelled as a "point" with finite surface/volume,
equivalent to that of a sphere with the given `diameter`. It's start-,
midpoint-, and end-coordinates are therefore identical.
"""
return self._z
@property
def end_distance(self):
"""
The distance to the root of the morphology at the end of this section.
Note that since a `Soma` is modeled as a point (see docs of `x`, etc.),
it does not add anything to the total distance, e.g. a section
connecting to a `Soma` has a `distance` of 0 um at its start.
"""
dist = self._parent.end_distance if self._parent is not None else 0 * um
return dist
[docs]
class Section(Morphology):
"""
A section (unbranched structure), described as a sequence of truncated
cones with potentially varying diameters and lengths per compartment.
Parameters
----------
diameter : `Quantity`
Either a single value (the constant diameter along the whole section),
or a value of length ``n+1``. When ``n+1`` values are given, they
will be interpreted as the diameters at the start of the first
compartment and the diameters at the end of each compartment (which is
equivalent to: the diameter at the start of each compartment and the
diameter at the end of the last compartment.
n : int, optional
The number of compartments in this section. Defaults to 1.
length : `Quantity`, optional
Either a single value (the total length of the section), or a value of
length ``n``, the length of each individual compartment. Cannot be
combined with the specification of coordinates.
x : `Quantity`, optional
``n+1`` values, specifying the x coordinates of the start point of the
first compartment and the end-points of all compartments (which is
equivalent to: the start point of all compartments and the end point of
the last compartment). The coordinates are interpreted as relative to
the end point of the parent compartment (if any), so in most cases the
start point should be ``0*um``. The common exception is a cylinder
connecting to a `Soma`, here the start point can be used to make the
cylinder start at the surface of the sphere instead of at its center.
You can specify all of ``x``, ``y``, or ``z`` to specify
a morphology in 3D, or only one or two out of them to specify a
morphology in 1D or 2D.
y : `Quantity`, optional
See ``x``
z : `Quantity`, optional
See ``x``
type : str, optional
The type (e.g. ``"axon"``) of this `Section`.
"""
@check_units(
n=1,
length=meter,
diameter=meter,
start_diameter=meter,
x=meter,
y=meter,
z=meter,
)
def __init__(
self,
diameter,
n=1,
length=None,
x=None,
y=None,
z=None,
start_diameter=None,
origin=None,
type=None,
):
n = int(n)
Morphology.__init__(self, n=n, type=type)
if diameter.ndim != 1 or len(diameter) != n + 1:
raise TypeError(
"The diameter argument has to be a one-dimensional array of length "
f"{int(n + 1)}"
)
self._diameter = Quantity(diameter, copy=True).reshape((n + 1,))
if (x is not None or y is not None or z is not None) and length is not None:
raise TypeError("Cannot specify coordinates and length at the same time.")
if length is not None:
# Length
if length.ndim != 1 or len(length) != n:
raise TypeError(
"The length argument has to be a one-dimensional array of length "
f"{int(n)}"
)
self._length = Quantity(length, copy=True).reshape((n,))
self._x = self._y = self._z = None
else:
# Coordinates
if x is None and y is None and z is None:
raise TypeError(
"No length specified, need to specify at least "
"one out of x, y, or z."
)
for name, value in [("x", x), ("y", y), ("z", z)]:
if value is not None and (value.ndim != 1 or len(value) != n + 1):
raise TypeError(
f"'{name}' needs to be a 1-dimensional array of length {n + 1}."
)
self._x = (
np.asarray(x).reshape((n + 1,)) if x is not None else np.zeros(n + 1)
)
self._y = (
np.asarray(y).reshape((n + 1,)) if y is not None else np.zeros(n + 1)
)
self._z = (
np.asarray(z).reshape((n + 1,)) if z is not None else np.zeros(n + 1)
)
length = np.sqrt(
(self.end_x - self.start_x) ** 2
+ (self.end_y - self.start_y) ** 2
+ (self.end_z - self.start_z) ** 2
)
self._length = length
def __repr__(self):
if all(
np.abs(self.end_diameter - self.end_diameter[0])
< self.end_diameter[0] * 1e-12
):
# Constant diameter
diam = self.end_diameter[0]
else:
diam = (
np.hstack(
[np.asarray(self.start_diameter[0]), np.asarray(self.end_diameter)]
)
* meter
)
s = f"{self.__class__.__name__}(diameter={diam!r}"
if self.n != 1:
s += f", n={self.n}"
if self._x is not None:
s += f", x={self._x!r}, y={self._y!r}, z={self._z!r}"
else:
s += f", length={sum(self._length)!r}"
if self.type is not None:
s += f", type={self.type!r}"
return f"{s})"
[docs]
def copy_section(self):
if self.x is None:
x, y, z = None, None, None
length = self.length
else:
x, y, z = self._x * meter, self._y * meter, self._z * meter
length = None
return Section(
diameter=self._diameter,
n=self.n,
x=x,
y=y,
z=z,
length=length,
type=self.type,
)
@property
def area(self):
r"""
The membrane surface area of each compartment in this section. The
surface area of each compartment is calculated as
:math:`\frac{\pi}{2}(d_1 + d_2)\sqrt{\frac{(d_1 - d_2)^2}{4} + l^2)}`,
where :math:`l` is the length of the compartment, and :math:`d_1` and
:math:`d_2` are the diameter at the start and end of the compartment,
respectively. Note that this surface area does not contain the area of
the two disks at the two sides of the truncated cone.
"""
d_1 = self.start_diameter
d_2 = self.end_diameter
return (
np.pi / 2 * (d_1 + d_2) * np.sqrt(((d_1 - d_2) ** 2) / 4 + self._length**2)
)
@property
def volume(self):
r"""
The volume of each compartment in this section. The volume of each
compartment is calculated as
:math:`\frac{\pi}{12} l (d_1^2 + d_1 d_2 + d_2^2)`,
where :math:`l` is the length of the compartment, and :math:`d_1` and
:math:`d_2` are the diameter at the start and end of the compartment,
respectively.
"""
d_1 = self.start_diameter
d_2 = self.end_diameter
return np.pi * self._length * (d_1**2 + d_1 * d_2 + d_2**2) / 12
@property
def length(self):
"""
The length of each compartment in this section.
"""
return self._length
@property
def start_diameter(self):
"""
The diameter at the start of each compartment in this section.
"""
return Quantity(self._diameter[:-1], copy=True)
@property
def end_diameter(self):
"""
The diameter at the end of each compartment in this section.
"""
return Quantity(self._diameter[1:], copy=True)
@property
def diameter(self):
"""
The diameter at the middle of each compartment in this section.
"""
d_1 = self.start_diameter
d_2 = self.end_diameter
# Diameter at the center
return 0.5 * (d_1 + d_2)
@property
def distance(self):
"""
The total distance between the midpoint of each compartment and the root
of the morphology.
"""
dist = self._parent.end_distance if self._parent is not None else 0 * um
return dist + np.cumsum(self.length) - 0.5 * self.length
@property
def end_distance(self):
"""
The distance to the root of the morphology at the end of this section.
"""
return self.distance[-1] + 0.5 * self.length[-1]
@property
def r_length_1(self):
"""
The geometry-dependent term to calculate the conductance between the
start and the midpoint of each compartment. Dividing this value by the
Intracellular resistivity gives the conductance.
"""
d_1 = self.start_diameter
d_2 = (self.start_diameter + self.end_diameter) * 0.5
return np.pi / 2 * (d_1 * d_2) / self._length
@property
def r_length_2(self):
"""
The geometry-dependent term to calculate the conductance between the
midpoint and the end of each compartment. Dividing this value by the
Intracellular resistivity gives the conductance.
"""
d_1 = (self.start_diameter + self.end_diameter) * 0.5
d_2 = self.end_diameter
return np.pi / 2 * (d_1 * d_2) / self._length
@property
def start_x_(self):
"""
The x coordinate (as a unitless floating point number) at the beginning
of each compartment. Returns ``None`` for morphologies without
coordinates.
"""
if self._x is None:
return None
parent_end_x = self.parent.end_x_ if self.parent is not None else None
if parent_end_x is not None:
return parent_end_x[-1] + self._x[:-1]
else:
return self._x[:-1]
@property
def start_y_(self):
"""
The y coordinate (as a unitless floating point number) at the beginning
of each compartment. Returns ``None`` for morphologies without
coordinates.
"""
if self._y is None:
return None
parent_end_y = self.parent.end_y_ if self.parent is not None else None
if parent_end_y is not None:
return parent_end_y[-1] + self._y[:-1]
else:
return self._y[:-1]
@property
def start_z_(self):
"""
The z coordinate (as a unitless floating point number) at the beginning
of each compartment. Returns ``None`` for morphologies without
coordinates.
"""
if self._z is None:
return None
parent_end_z = self.parent.end_z_ if self.parent is not None else None
if parent_end_z is not None:
return parent_end_z[-1] + self._z[:-1]
else:
return self._z[:-1]
@property
def x_(self):
"""
The x coordinate (as a unitless floating point number) at the midpoint
of each compartment. Returns ``None`` for morphologies without
coordinates.
"""
if self._x is None:
return None
start_x = self.start_x_
diff_x = self.end_x_ - start_x
return start_x + 0.5 * diff_x
@property
def y_(self):
"""
The y coordinate (as a unitless floating point number) at the midpoint
of each compartment. Returns ``None`` for morphologies without
coordinates.
"""
if self._y is None:
return None
start_y = self.start_y_
diff_y = self.end_y_ - start_y
return start_y + 0.5 * diff_y
@property
def z_(self):
"""
The z coordinate (as a unitless floating point number) at the midpoint
of each compartment. Returns ``None`` for morphologies without
coordinates.
"""
if self._z is None:
return None
start_z = self.start_z_
diff_z = self.end_z_ - start_z
return start_z + 0.5 * diff_z
@property
def end_x_(self):
"""
The x coordinate (as a unitless floating point number) at the end of
each compartment. Returns ``None`` for morphologies without coordinates.
"""
if self._x is None:
return None
parent_end_x = self.parent.end_x_ if self.parent is not None else None
if parent_end_x is not None:
return parent_end_x[-1] + self._x[1:]
else:
return self._x[1:]
@property
def end_y_(self):
"""
The y coordinate (as a unitless floating point number) at the end of
each compartment. Returns ``None`` for morphologies without coordinates.
"""
if self._y is None:
return None
parent_end_y = self.parent.end_y_ if self.parent is not None else None
if parent_end_y is not None:
return parent_end_y[-1] + self._y[1:]
else:
return self._y[1:]
@property
def end_z_(self):
"""
The z coordinate (as a unitless floating point number) at the end of
each compartment. Returns ``None`` for morphologies without coordinates.
"""
if self._z is None:
return None
parent_end_z = self.parent.end_z_ if self.parent is not None else None
if parent_end_z is not None:
return parent_end_z[-1] + self._z[1:]
else:
return self._z[1:]
[docs]
class Cylinder(Section):
"""
A cylindrical section. For sections with more complex geometry (varying
length and/or diameter of each compartment), use the `Section` class.
Parameters
----------
diameter : `Quantity`
The diameter of the cylinder.
n : int, optional
The number of compartments in this section. Defaults to 1.
length : `Quantity`, optional
The length of the cylinder. Cannot be combined with the specification
of coordinates.
x : `Quantity`, optional
A sequence of two values, the start and the end point of the cylinder.
The coordinates are interpreted as relative to the end point of the
parent compartment (if any), so in most cases the start point should
be ``0*um``. The common exception is a cylinder connecting to a `Soma`,
here the start point can be used to make the cylinder start at the
surface of the sphere instead of at its center.
You can specify all of ``x``, ``y``, or ``z`` to specify
a morphology in 3D, or only one or two out of them to specify a
morphology in 1D or 2D.
y : `Quantity`, optional
See ``x``
z : `Quantity`, optional
See ``x``
type : str, optional
The type (e.g. ``"axon"``) of this `Cylinder`.
"""
@check_units(n=1, length=meter, diameter=meter, x=meter, y=meter, z=meter)
def __init__(self, diameter, n=1, length=None, x=None, y=None, z=None, type=None):
n = int(n)
Morphology.__init__(self, n=n, type=type)
# Diameter
if diameter.shape != () and (diameter.ndim > 1 or len(diameter) != 1):
raise TypeError("The diameter argument has to be a single value.")
diameter = np.ones(n) * diameter
self._diameter = diameter
if (x is not None or y is not None or z is not None) and length is not None:
raise TypeError("Cannot specify coordinates and length at the same time.")
if length is not None:
# Length
if length.shape != () and (length.ndim > 1 or len(length) != 1):
raise TypeError("The length argument has to be a single value.")
self._length = np.ones(n) * (length / n) # length was total length
self._x = self._y = self._z = None
else:
# Coordinates
if x is None and y is None and z is None:
raise TypeError(
"No length specified, need to specify at least "
"one out of x, y, or z."
)
for name, value in [("x", x), ("y", y), ("z", z)]:
if value is not None and (value.ndim != 1 or len(value) != 2):
raise TypeError(
f"{name} needs to be a 1-dimensional array of length 2 (start "
"and end point)"
)
self._x = (
np.asarray(np.linspace(x[0], x[1], n + 1))
if x is not None
else np.zeros(n + 1)
)
self._y = (
np.asarray(np.linspace(y[0], y[1], n + 1))
if y is not None
else np.zeros(n + 1)
)
self._z = (
np.asarray(np.linspace(z[0], z[1], n + 1))
if z is not None
else np.zeros(n + 1)
)
length = np.sqrt(
(self.end_x - self.start_x) ** 2
+ (self.end_y - self.start_y) ** 2
+ (self.end_z - self.start_z) ** 2
)
self._length = length
def __repr__(self):
s = f"{self.__class__.__name__}(diameter={self.diameter[0]!r}"
if self.n != 1:
s += f", n={self.n}"
if self._x is not None:
s += (
f", x={self._x[[0, -1]]!r}, y={self._y[[0, -1]]!r},"
f" z={self._z[[0, -1]]!r}"
)
else:
s += f", length={sum(self._length)!r}"
if self.type is not None:
s += f", type={self.type!r}"
return f"{s})"
[docs]
def copy_section(self):
if self.x is None:
return Cylinder(
self.diameter[0], n=self.n, length=self.length, type=self.type
)
else:
return Cylinder(
self.diameter[0],
n=self.n,
x=self._x[[0, -1]],
y=self._y[[0, -1]],
z=self._z[[0, -1]],
type=self.type,
)
# Overwrite the properties that differ from `Section`
@property
def area(self):
r"""
The membrane surface area of each compartment in this section. The
surface area of each compartment is calculated as
:math:`\pi d l`,
where :math:`l` is the length of the compartment, and :math:`d` is its
diameter. Note that this surface area does not contain the area of
the two disks at the two sides of the cylinder.
"""
return np.pi * self._diameter * self.length
@property
def start_diameter(self):
"""
The diameter at the start of each compartment in this section.
"""
return self._diameter
@property
def diameter(self):
"""
The diameter at the middle of each compartment in this section.
"""
return self._diameter
@property
def end_diameter(self):
"""
The diameter at the end of each compartment in this section.
"""
return self._diameter
@property
def volume(self):
r"""
The volume of each compartment in this section. The volume of each
compartment is calculated as
:math:`\pi \frac{d}{2}^2 l` ,
where :math:`l` is the length of the compartment, and :math:`d` is its
diameter.
"""
return np.pi * (self._diameter / 2) ** 2 * self.length
@property
def r_length_1(self):
"""
The geometry-dependent term to calculate the conductance between the
start and the midpoint of each compartment. Dividing this value by the
Intracellular resistivity gives the conductance.
"""
return np.pi / 2 * (self._diameter**2) / self.length
@property
def r_length_2(self):
"""
The geometry-dependent term to calculate the conductance between the
midpoint and the end of each compartment. Dividing this value by the
Intracellular resistivity gives the conductance.
"""
return np.pi / 2 * (self._diameter**2) / self.length