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, 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 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.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) 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): """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 60, store diameter if sim.neuron.h.distance(1, 0.5, sec=section) > 60: 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.nseg = 1 section.L = 30 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('Replace axon with AIS')
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) 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)