Source code for bluepyopt.ephys.morphologies

"""Morphology classes"""

"""
Copyright (c) 2016-2020, EPFL/Blue Brain Project

 This file is part of BluePyOpt <https://github.com/BlueBrain/BluePyOpt>

 This library is free software; you can redistribute it and/or modify it under
 the terms of the GNU Lesser General Public License version 3.0 as published
 by the Free Software Foundation.

 This library is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
 details.

 You should have received a copy of the GNU Lesser General Public License
 along with this library; if not, write to the Free Software Foundation, Inc.,
 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
"""

# pylint: disable=W0511

import os
import platform
import logging
import pathlib
import bisect
import numpy
from bluepyopt.ephys.base import BaseEPhys
from bluepyopt.ephys.serializer import DictMixin
from bluepyopt.ephys.acc import arbor, ArbLabel

logger = logging.getLogger(__name__)

# TODO define an addressing scheme



[docs]
class Morphology(BaseEPhys):

    """Morphology class"""
    pass




[docs]
class NrnFileMorphology(Morphology, DictMixin):

    """Morphology loaded from a file"""
    SERIALIZED_FIELDS = ('morphology_path', 'do_replace_axon', 'do_set_nseg',
                         'replace_axon_hoc', 'nseg_frequency',
                         'morph_modifiers', 'morph_modifiers_hoc',
                         'morph_modifiers_kwargs')

    def __init__(
            self,
            morphology_path,
            do_replace_axon=False,
            do_set_nseg=True,
            comment='',
            replace_axon_hoc=None,
            axon_stub_length=60,
            axon_nseg_frequency=40,
            nseg_frequency=40,
            morph_modifiers=None,
            morph_modifiers_hoc=None,
            morph_modifiers_kwargs=None):
        """Constructor

        Args:
            morphology_path (str or Path): location of the file describing the
                morphology
            do_replace_axon (bool): Does the axon need to be replaced by an AIS
                stub with default function ?
            replace_axon_hoc (str): Translation in HOC language for the
                'replace_axon' method. This code will 'only' be used when
                calling create_hoc on a cell model. While the model is run in
                python, replace_axon is used instead. Must include
                'proc replace_axon(){ ... }
                If None,the default replace_axon is used
            axon_stub_length (float): Length of replacement axon
            axon_nseg_frequency (int): frequency of nseg, for axon
            nseg_frequency (float): frequency of nseg
            do_set_nseg (bool): if True, it will use nseg_frequency
            morph_modifiers (list): list of functions to modify the icell
                with (sim, icell) as arguments
            morph_modifiers_hoc (list): list of hoc strings corresponding
                to morph_modifiers
            morph_modifiers_kwargs (dict): kwargs for morph_modifiers functions
        """
        name = os.path.basename(morphology_path)
        super(NrnFileMorphology, self).__init__(name=name, comment=comment)
        # TODO speed up loading of morphologies from files
        # Path to morphology
        if isinstance(morphology_path, pathlib.Path):
            morphology_path = str(morphology_path)
        self.morphology_path = morphology_path
        self.do_replace_axon = do_replace_axon
        self.axon_stub_length = axon_stub_length
        self.axon_nseg_frequency = axon_nseg_frequency
        self.do_set_nseg = do_set_nseg
        self.nseg_frequency = nseg_frequency
        self.morph_modifiers = morph_modifiers
        self.morph_modifiers_hoc = morph_modifiers_hoc
        self.morph_modifiers_kwargs = morph_modifiers_kwargs
        if self.morph_modifiers_kwargs is None:
            self.morph_modifiers_kwargs = {}

        if replace_axon_hoc is None:
            self.replace_axon_hoc = self.default_replace_axon_hoc
        else:
            self.replace_axon_hoc = replace_axon_hoc

    def __str__(self):
        """Return string representation"""

        return self.morphology_path


[docs]
    def instantiate(self, sim=None, icell=None):
        """Load morphology"""

        logger.debug('Loading morphology %s', self.morphology_path)

        if not os.path.exists(self.morphology_path):
            raise IOError(
                'Morphology not found at \'%s\'' %
                self.morphology_path)

        sim.neuron.h.load_file('stdrun.hoc')
        sim.neuron.h.load_file('import3d.hoc')

        extension = self.morphology_path.split('.')[-1]

        if extension.lower() == 'swc':
            imorphology = sim.neuron.h.Import3d_SWC_read()
        elif extension.lower() == 'asc':
            imorphology = sim.neuron.h.Import3d_Neurolucida3()
        else:
            raise ValueError("Unknown filetype: %s" % extension)

        # TODO this is to get rid of stdout print of neuron
        # probably should be more intelligent here, and filter out the
        # lines we don't want

        imorphology.quiet = 1

        if platform.system() == 'Windows':
            sim.neuron.h.hoc_stdout('NUL')
        else:
            sim.neuron.h.hoc_stdout('/dev/null')

        imorphology.input(str(self.morphology_path))
        sim.neuron.h.hoc_stdout()

        morphology_importer = sim.neuron.h.Import3d_GUI(imorphology, 0)

        morphology_importer.instantiate(icell)

        # TODO Set nseg should be called after all the parameters have been
        # set
        # (in case e.g. Ra was changed)
        if self.do_set_nseg:
            self.set_nseg(icell)

        if self.do_replace_axon:
            self.replace_axon(sim=sim, icell=icell,
                              axon_stub_length=self.axon_stub_length,
                              axon_nseg_frequency=self.axon_nseg_frequency)

        if self.morph_modifiers is not None:
            for morph_modifier in self.morph_modifiers:
                morph_modifier(sim=sim, icell=icell,
                               **self.morph_modifiers_kwargs)



[docs]
    def destroy(self, sim=None):
        """Destroy morphology instantiation"""
        pass



[docs]
    def set_nseg(self, icell):
        """Set the nseg of every section"""
        for section in icell.all:
            section.nseg = 1 + 2 * int(section.L / self.nseg_frequency)



[docs]
    @staticmethod
    def replace_axon(sim=None, icell=None,
                     axon_stub_length=60, axon_nseg_frequency=40):
        """Replace axon"""

        nsec = len([sec for sec in icell.axonal])

        if nsec == 0:
            ais_diams = [1, 1]
        elif nsec == 1:
            ais_diams = [icell.axon[0].diam, icell.axon[0].diam]
        else:
            ais_diams = [icell.axon[0].diam, icell.axon[0].diam]
            # Define origin of distance function
            sim.neuron.h.distance(0, 0.5, sec=icell.soma[0])

            for section in icell.axonal:
                # If distance to soma is larger than
                # axon_stub_length, store diameter
                if sim.neuron.h.distance(1, 0.5, sec=section) \
                        > axon_stub_length:
                    ais_diams[1] = section.diam
                    break

        for section in icell.axonal:
            sim.neuron.h.delete_section(sec=section)

        # Create new axon array
        sim.neuron.h.execute('create axon[2]', icell)

        for index, section in enumerate(icell.axon):
            section.L = axon_stub_length / 2
            section.nseg = 1 + 2 * int(section.L / axon_nseg_frequency)
            section.diam = ais_diams[index]
            icell.axonal.append(sec=section)
            icell.all.append(sec=section)

        icell.axon[0].connect(icell.soma[0], 1.0, 0.0)
        icell.axon[1].connect(icell.axon[0], 1.0, 0.0)

        logger.debug(f"Replace axon with AIS, {axon_stub_length=}")


    default_replace_axon_hoc = \
        '''
proc replace_axon(){ local nSec, D1, D2
  // preserve the number of original axonal sections
  nSec = sec_count(axonal)

  // Try to grab info from original axon
  if(nSec == 0) { //No axon section present
    D1 = D2 = 1
  } else if(nSec == 1) {
    axon[0] D1 = D2 = diam
  } else {
    axon[0] D1 = D2 = diam
    soma distance() //to calculate distance from soma
    forsec axonal{
      //if section is longer than 60um then store diam and exit from loop
      if(distance(0.5) > 60){
        D2 = diam
        break
      }
    }
  }

  // get rid of the old axon
  forsec axonal{
    delete_section()
  }

  create axon[2]

  axon[0] {
    L = 30
    diam = D1
    nseg = 1 + 2*int(L/40)
    all.append()
    axonal.append()
  }
  axon[1] {
    L = 30
    diam = D2
    nseg = 1 + 2*int(L/40)
    all.append()
    axonal.append()
  }
  nSecAxonal = 2
  soma[0] connect axon[0](0), 1
  axon[0] connect axon[1](0), 1
}
        '''




[docs]
class ArbFileMorphology(Morphology, DictMixin):
    """Arbor morphology utilities"""

    # Arbor morphology tags
    tags = dict(
        soma=1,
        axon=2,
        dend=3,
        apic=4,
        myelin=5
    )

    # Correspondence of BluePyOpt seclists to Arbor region labels
    # (renaming locations according to SWC convention: using
    # 'dend' for basal dendrite, 'apic' for apical dendrite)
    region_labels = dict(
        all=ArbLabel(
            type='region', name='all', s_expr='(all)'),
        somatic=ArbLabel(
            type='region', name='soma', s_expr='(tag %i)' % tags['soma']),
        axonal=ArbLabel(
            type='region', name='axon', s_expr='(tag %i)' % tags['axon']),
        basal=ArbLabel(
            type='region', name='dend', s_expr='(tag %i)' % tags['dend']),
        apical=ArbLabel(
            type='region', name='apic', s_expr='(tag %i)' % tags['apic']),
        myelinated=ArbLabel(
            type='region', name='myelin', s_expr='(tag %i)' % tags['myelin']),
    )


[docs]
    @staticmethod
    def load(morpho_filename, replace_axon):
        '''Load morphology and optionally perform axon replacement

        Args:
            morpho_filename (str): Path to file with original morphology.
            replace_axon (): Path to/ACC string for morphology to replace
            axon with (if not None).
        '''

        morpho_suffix = pathlib.Path(morpho_filename).suffix

        if morpho_suffix == '.acc':
            morpho = arbor.load_component(morpho_filename).component
        elif morpho_suffix == '.swc':
            morpho = arbor.load_swc_arbor(morpho_filename)
            # turn loaded_morphology into morphology type
            morpho = morpho.morphology
        elif morpho_suffix == '.asc':
            morpho = arbor.load_asc(morpho_filename).morphology
        else:
            raise RuntimeError(
                'Unsupported morphology %s' % morpho_filename +
                ' (only .swc and .asc supported)')

        if replace_axon is not None:
            replacement = arbor.load_component(replace_axon).component
            morpho = ArbFileMorphology.replace_axon(morpho, replacement)

        return morpho



[docs]
    @staticmethod
    def extract_nrn_seclists(icell, seclists):
        '''Extract section lists from an instantiated cell (axon replacement)

        Args:
            icell (): Instantiated cell model in the NEURON simulator.
            seclists (): List of section lists to extract
            (typically ['axon'] or ['axon', 'myelin']).
        '''
        replace_axon = arbor.segment_tree()
        nrn_seg_to_dist = dict()
        nrn_seg_to_arb_seg = dict()
        for sec in seclists:
            for section in getattr(icell, sec):

                if replace_axon.size == 0:  # root
                    arb_parent_seg = arbor.mnpos
                else:
                    parent_seg = section.parentseg()
                    parent_sec = parent_seg.sec.name()
                    parent_x = parent_seg.x
                    parent_seg_id = bisect.bisect_left(
                        nrn_seg_to_dist[parent_sec],
                        parent_x)
                    arb_parent_seg = \
                        nrn_seg_to_arb_seg[parent_sec][parent_seg_id]

                pts3d = section.psection()['morphology']['pts3d']
                if len(pts3d) == 0:
                    # stylized geometry, must use sim.neuron.h.define_shape()
                    raise ValueError('Before exporting to ACC, embed'
                                     ' stylized geometry in 3d'
                                     ' by instantiating morphology with'
                                     ' cell_model.instantiate_morphology_3d.')

                pts3d = numpy.array(pts3d)
                dist_x = numpy.cumsum(
                    numpy.linalg.norm(
                        pts3d[1:, :3] - pts3d[:-1, :3], axis=1)) /\
                    section.psection()['morphology']['L']
                relative_length_err = abs(1. - dist_x[-1])
                if relative_length_err > 1e-4:
                    logger.warn('pts3d length does not add up to'
                                ' section length, relative error = %s' %
                                relative_length_err)
                    if relative_length_err > 1e-2:
                        raise ValueError('pts3d length inconsistent'
                                         ' with section length, relative'
                                         ' error = %s' %
                                         relative_length_err)

                dist_x[-1] = 1.
                nrn_seg_to_dist[section.name()] = dist_x

                arb_seg_ids = []
                for i in range(1, len(pts3d)):
                    prox = pts3d[i - 1]
                    dist = pts3d[i]
                    arb_parent_seg = replace_axon.append(
                        arb_parent_seg,
                        arbor.mpoint(*prox[:3], 0.5 * prox[3]),
                        arbor.mpoint(*dist[:3], 0.5 * dist[3]),
                        ArbFileMorphology.tags[sec])
                    arb_seg_ids.append(arb_parent_seg)
                # dist, arb_seg_id pairs
                nrn_seg_to_arb_seg[section.name()] = arb_seg_ids

        replace_axon = arbor.morphology(replace_axon)

        return replace_axon



[docs]
    @staticmethod
    def replace_axon(morphology, replacement=None):
        '''return a morphology with the axon replaced by another morphology

        Args:
            morphology (arbor.morphology): The original Arbor morphology
            replacement (): An Arbor morphology to replace the axon with
        '''

        # Check if tag_roots is available
        if not hasattr(arbor.segment_tree, 'tag_roots'):
            raise NotImplementedError(
                "Need a newer version of Arbor for axon replacement.")

        # Arbor tags
        axon_tag = ArbFileMorphology.tags['axon']

        # prune morphology at axon root
        st = morphology.to_segment_tree()
        axon_roots = st.tag_roots(axon_tag)
        if len(axon_roots) > 1:
            raise ValueError("Axon replacement is only supported for "
                             "morphologies with a single axon root.")
        elif len(axon_roots) == 1:
            axon_root = axon_roots[0]
            logger.debug('Axon replacement: splitting segment tree'
                         ' at segment %d.', axon_root)
            pruned_st, axon_st = st.split_at(axon_root)
            axon_parent = st.parents[axon_root]
        else:
            pruned_st = st
            axon_parent = arbor.mnpos

        # join pruned segment tree and replacement at axon parent
        axon_replacement_st = replacement.to_segment_tree()

        logger.debug('Axon replacement: joining replacement onto'
                     'pruned tree at parent segment %d.', axon_parent)
        joined_st = pruned_st.join_at(
            axon_parent, axon_replacement_st)

        return arbor.morphology(joined_st)