Source code for hnn_core.cell

"""Establish class def for general cell features."""

# Authors: Mainak Jas <mjas@mgh.harvard.edu>
#          Sam Neymotin <samnemo@gmail.com>

from copy import deepcopy

import numpy as np
from numpy.linalg import norm

from neuron import h, nrn

from .viz import plot_cell_morphology
from .externals.mne import _validate_type, _check_option

# Units for e: mV
# Units for gbar: S/cm^2


def _get_cos_theta(sections, sec_name_apical):
    """Get cos(theta) to compute dipole along the apical dendrite."""
    a = (np.array(sections[sec_name_apical].end_pts[1]) -
         np.array(sections[sec_name_apical].end_pts[0]))
    cos_thetas = dict()
    for sec_name, section in sections.items():
        b = np.array(section.end_pts[1]) - np.array(section.end_pts[0])
        cos_thetas[sec_name] = np.dot(a, b) / (norm(a) * norm(b))
    return cos_thetas


def _calculate_gaussian(x_val, height, lamtha):
    """Return height of gaussian at x_val.

    Parameters
    ----------
    x_val : float
        Value on x-axis to query height of gaussian curve.
    height : float
        Height of the gaussian curve at zero.
    lamtha : float
        Space constant.

    Returns
    -------
    x_height : float
        Height of gaussian at x_val.

    Notes
    -----
    Gaussian curve is centered at zero and has a fixed peak height
    such the _calculate_gaussian(0, lamtha) returns 1 for all lamtha.
    """
    x_height = height * np.exp(-(x_val**2) / (lamtha**2))

    return x_height


def _get_gaussian_connection(src_pos, target_pos, nc_dict,
                             inplane_distance=1.):
    """Calculate distance dependent connection properties.

    Parameters
    ----------
    src_pos : float
        Position of source cell.
    target_pos : float
        Position of target cell.
    nc_dict : dict
        Dictionary with keys: pos_src, A_weight, A_delay, lamtha
        Defines the connection parameters
    inplane_distance : float
        The in plane-distance (in um) between pyramidal cell somas in the
        square grid. Default: 1.0 um.

    Returns
    -------
    weight : float
        Weight of the synaptic connection.
    delay : float
        Delay of synaptic connection.

    Notes
    -----
    Distance in xy plane is used for gaussian decay.
    """
    x_dist = target_pos[0] - src_pos[0]
    y_dist = target_pos[1] - src_pos[1]
    cell_dist = np.sqrt(x_dist**2 + y_dist**2)
    scaled_lamtha = nc_dict['lamtha'] * inplane_distance

    weight = _calculate_gaussian(
        cell_dist, nc_dict['A_weight'], scaled_lamtha)
    delay = nc_dict['A_delay'] / _calculate_gaussian(
        cell_dist, 1, scaled_lamtha)
    return weight, delay


class _ArtificialCell:
    """The ArtificialCell class for initializing a NEURON feed source.

    Parameters
    ----------
    event_times : list
        Spike times associated with a single feed source (i.e.,
        associated with a unique gid).
    threshold : float
        Membrane potential threshold that demarks a spike.
    gid : int or None (optional)
        Each cell in a network is uniquely identified by it's "global ID": GID.
        The GID is an integer from 0 to n_cells, or None if the cell is not
        yet attached to a network. Once the GID is set, it cannot be changed.

    Attributes
    ----------
    nrn_eventvec : instance of h.Vector()
        NEURON h.Vector() object of event times.
    nrn_vecstim : instance of h.VecStim()
        NEURON h.VecStim() object of spike sources created
        from nrn_eventvec.
    nrn_netcon : instance of h.NetCon()
        NEURON h.NetCon() object that creates the spike
        source-to-target references for nrn_vecstim.
    gid : int
        GID of the cell in a network (or None if not yet assigned)
    """
    def __init__(self, event_times, threshold, gid=None):
        # Convert event times into nrn vector
        self.nrn_eventvec = h.Vector()
        self.nrn_eventvec.from_python(event_times)

        # load eventvec into VecStim object
        self.nrn_vecstim = h.VecStim()
        self.nrn_vecstim.play(self.nrn_eventvec)

        # create the cell and artificial NetCon
        self.nrn_netcon = h.NetCon(self.nrn_vecstim, None)
        self.nrn_netcon.threshold = threshold

        self._gid = None
        if gid is not None:
            self.gid = gid  # use setter method to check input argument gid

    @property
    def gid(self):
        return self._gid

    @gid.setter
    def gid(self, gid):
        if not isinstance(gid, int):
            raise ValueError('gid must be an integer')
        if self._gid is None:
            self._gid = gid
        else:
            raise RuntimeError('Global ID for this cell already assigned!')


class Section:
    """Section class.

    Parameters
    ----------
    L : float
        length of a section in microns.
    diam : float
        diameter of a section in microns.
    cm : float
        membrane capacitance in micro-Farads.
    Ra : float
        axial resistivity in ohm-cm
    end_pts : list of [x, y, z]
        The start and stop points of the section.

    Attributes
    ----------
    mechs : dict
        Mechanisms to insert in this section. The keys
        are the names of the mechanisms and values
        are the properties. For e.g., {'ca': {'gbar_ca': 60}}
    syns : list of str
        The synaptic mechanisms to add in this section
    end_pts : list of [x, y, z]
        The start and stop points of the section. Cannot be changed.
    L : float
        length of a section in microns.
    diam : float
        diameter of a section in microns.
    cm : float
        membrane capacitance in micro-Farads.
    Ra : float
        axial resistivity in ohm-cm.
    """
    def __init__(self, L, diam, Ra, cm, end_pts=None):

        self._L = L
        self._diam = diam
        self._Ra = Ra
        self._cm = cm
        if end_pts is None:
            end_pts = list()
        self._end_pts = end_pts

        self.mechs = dict()
        self.syns = list()

    def __repr__(self):
        return f'L={self.L}, diam={self.diam}, cm={self.cm}, Ra={self.Ra}'

    @property
    def L(self):
        return self._L

    @property
    def diam(self):
        return self._diam

    @property
    def cm(self):
        return self._cm

    @property
    def Ra(self):
        return self._Ra

    @property
    def end_pts(self):
        return self._end_pts


[docs]class Cell: """Create a cell object. Parameters ---------- name : str The name of the cell. pos : tuple The (x, y, z) coordinates. sections : dict of Section Dictionary with keys as section name. synapses : dict of dict Keys are name of synaptic mechanism. Each synaptic mechanism has keys for parameters of the mechanism, e.g., 'e', 'tau1', 'tau2'. topology : list of list The topology of cell sections. Each element is a list of 4 items in the format [parent_sec, parent_loc, child_sec, child_loc] where parent_sec and parent_loc are float between 0 and 1 specifying the location in the section to connect and parent_sec and child_sec are names of the connecting sections. sect_loc : dict of list Can have keys 'proximal' or 'distal' each containing names of section locations that are proximal or distal. gid : int or None (optional) Each cell in a network is uniquely identified by it's "global ID": GID. The GID is an integer from 0 to n_cells, or None if the cell is not yet attached to a network. Once the GID is set, it cannot be changed. Attributes ---------- pos : list of length 3 The position of the cell. sections : nested dict The section parameters. The key is the name of the section and the value is a dictionary parametrizing the morphology of the section and the mechanisms inserted. synapses : dict The synapses that the cell can use for connections. dipole_pp : list of h.Dipole() The Dipole objects (see dipole.mod). vsec : dict Recording of section specific voltage. Must be enabled by running simulate_dipole(net, record_vsec=True) or simulate_dipole(net, record_vsoma=True) isec : dict Contains recording of section specific currents indexed by synapse type (keys can be soma_gabaa, soma_gabab etc.). Must be enabled by running simulate_dipole(net, record_isec=True) or simulate_dipole(net, record_isoma=True) tonic_biases : list of h.IClamp The current clamps inserted at each section of the cell for tonic biasing inputs. gid : int GID of the cell in a network (or None if not yet assigned) sect_loc : dict of list Can have keys 'proximal' or 'distal' each containing names of section locations that are proximal or distal. Examples -------- >>> section_soma = Section( L=39, diam=20, cm=0.85, Ra=200., end_pts=[[0, 0, 0], [0, 39., 0]] ) """ def __init__(self, name, pos, sections, synapses, topology, sect_loc, gid=None): self.name = name self.pos = pos for section in sections.values(): if not isinstance(section, Section): raise ValueError(f'Items in section must be instances' f' of Section. Got {type(section)}') self.sections = sections self.synapses = synapses self.topology = topology self.sect_loc = sect_loc self._nrn_sections = dict() self._nrn_synapses = dict() self.dipole_pp = list() self.vsec = dict() self.isec = dict() # insert iclamp self.list_IClamp = list() self._gid = None self.tonic_biases = list() if gid is not None: self.gid = gid # use setter method to check input argument gid self._update_end_pts()
[docs] def __repr__(self): class_name = self.__class__.__name__ return f'<{class_name} | gid={self._gid}>'
@property def gid(self): return self._gid @gid.setter def gid(self, gid): if not isinstance(gid, int): raise ValueError('gid must be an integer') if self._gid is None: self._gid = gid else: raise RuntimeError('Global ID for this cell already assigned!') def _set_biophysics(self, sections): "Set the biophysics for the cell." # neuron syntax is used to set values for mechanisms # sec.gbar_mech = x sets value of gbar for mech to x for all segs # in a section. This method is significantly faster than using # a for loop to iterate over all segments to set mech values # If value depends on distance from the soma. Soma is set as # origin by passing cell.soma as a sec argument to h.distance() # Then iterate over segment nodes of dendritic sections # and set attribute depending on h.distance(seg.x), which returns # distance from the soma to this point on the CURRENTLY ACCESSED # SECTION!!! h.distance(sec=self._nrn_sections['soma']) for sec_name, section in sections.items(): sec = self._nrn_sections[sec_name] for mech_name, p_mech in section.mechs.items(): sec.insert(mech_name) for attr, val in p_mech.items(): if hasattr(val, '__call__'): sec.push() for seg in sec: setattr(seg, attr, val(h.distance(seg.x))) h.pop_section() else: setattr(sec, attr, val) def _create_synapses(self, sections, synapses): """Create synapses.""" for sec_name in sections: for receptor in sections[sec_name].syns: syn_key = f'{sec_name}_{receptor}' seg = self._nrn_sections[sec_name](0.5) self._nrn_synapses[syn_key] = self.syn_create( seg, **synapses[receptor]) def _create_sections(self, sections, topology): """Create soma and set geometry. Notes ----- By default neuron uses xy plane for height and xz plane for depth. This is opposite for model as a whole, but convention is followed in this function ease use of gui. """ if 'soma' not in self.sections: raise KeyError('soma must be defined for cell') # shift cell to self.pos and reorient apical dendrite # along z direction of self.pos dx = self.pos[0] - self.sections['soma'].end_pts[0][0] dy = self.pos[1] - self.sections['soma'].end_pts[0][1] dz = self.pos[2] - self.sections['soma'].end_pts[0][2] for sec_name in sections: sec = h.Section(name=f'{self.name}_{sec_name}') self._nrn_sections[sec_name] = sec h.pt3dclear(sec=sec) h.pt3dconst(0, sec=sec) # be explicit, see documentation for pt in sections[sec_name].end_pts: h.pt3dadd(pt[0] + dx, pt[1] + dy, pt[2] + dz, 1, sec=sec) # with pt3dconst==0, these will alter the 3d points defined above! sec.L = sections[sec_name].L sec.diam = sections[sec_name].diam sec.Ra = sections[sec_name].Ra sec.cm = sections[sec_name].cm if sec.L > 100.: # 100 um sec.nseg = int(sec.L / 50.) # make dend.nseg odd for all sections if not sec.nseg % 2: sec.nseg += 1 if topology is None: topology = list() # Connects sections of THIS cell together. for connection in topology: parent_sec = self._nrn_sections[connection[0]] child_sec = self._nrn_sections[connection[2]] parent_loc = connection[1] child_loc = connection[3] child_sec.connect(parent_sec, parent_loc, child_loc) # be explicit about letting sec.L dominate over the 3d points used by # h.pt3dadd(); see # https://nrn.readthedocs.io/en/latest/python/modelspec/programmatic/topology/geometry.html?highlight=pt3dadd#pt3dadd # noqa h.define_shape()
[docs] def build(self, sec_name_apical=None): """Build cell in Neuron and insert dipole if applicable. Parameters ---------- sec_name_apical : str | None If not None, a dipole will be inserted in this cell in alignment with this section. The section should belong to the apical dendrite of a pyramidal neuron. """ self._create_sections(self.sections, self.topology) self._create_synapses(self.sections, self.synapses) self._set_biophysics(self.sections) if sec_name_apical in self._nrn_sections: self._insert_dipole(sec_name_apical) elif sec_name_apical is not None: raise ValueError(f'sec_name_apical must be an existing ' f'section of the current cell or None. ' f'Got {sec_name_apical}.')
[docs] def copy(self): """Return copy of instance.""" return deepcopy(self)
# two things need to happen here for h: # 1. dipole needs to be inserted into each section # 2. a list needs to be created with a Dipole (Point Process) in each # section at position 1 # In Cell() and not Pyr() for future possibilities def _insert_dipole(self, sec_name_apical): """Insert dipole into each section of this cell. Parameters ---------- sec_name_apical : str The name of the section along which dipole moment is calculated. """ self.dpl_vec = h.Vector(1) self.dpl_ref = self.dpl_vec._ref_x[0] cos_thetas = _get_cos_theta(self.sections, 'apical_trunk') # setting pointers and ztan values for sect_name in self.sections: sect = self._nrn_sections[sect_name] sect.insert('dipole') dpp = h.Dipole(1, sec=sect) # defined in dipole_pp.mod self.dipole_pp.append(dpp) dpp.ri = h.ri(1, sec=sect) # assign internal resistance # sets pointers in dipole mod file to the correct locations dpp._ref_pv = sect(0.99)._ref_v dpp._ref_Qtotal = self.dpl_ref # gives INTERNAL segments of the section, non-endpoints # creating this because need multiple values simultaneously pos_all = np.array([seg.x for seg in sect.allseg()]) seg_lens = np.diff(pos_all) * sect.L seg_lens_z = seg_lens * cos_thetas[sect_name] # alternative procedure below with y_long(itudinal) # y_long = (h.y3d(1, sec=sect) - h.y3d(0, sec=sect)) * pos # y_diff = np.diff(y_long) # doing range to index multiple values of the same # np.array simultaneously for idx, pos in enumerate(pos_all[1:-1]): # assign the ri value to the dipole # ri not defined at 0 and L sect(pos).dipole.ri = h.ri(pos, sec=sect) # range variable 'dipole' # set pointers to previous segment's voltage, with # boundary condition sect(pos).dipole._ref_pv = sect(pos_all[idx])._ref_v # set aggregate pointers sect(pos).dipole._ref_Qsum = dpp._ref_Qsum sect(pos).dipole._ref_Qtotal = self.dpl_ref # add ztan values sect(pos).dipole.ztan = seg_lens_z[idx] # set the pp dipole's ztan value to the last value from seg_lens_z dpp.ztan = seg_lens_z[-1] self.dipole = h.Vector().record(self.dpl_ref)
[docs] def create_tonic_bias(self, amplitude, t0, tstop, loc=0.5): """Create tonic bias at the soma. Parameters ---------- amplitude : float The amplitude of the input. t0 : float The start time of tonic input (in ms). tstop : float The end time of tonic input (in ms). loc : float (0 to 1) The location of the input in the soma section. """ stim = h.IClamp(self._nrn_sections['soma'](loc)) stim.delay = t0 stim.dur = tstop - t0 stim.amp = amplitude self.tonic_biases.append(stim)
[docs] def record(self, record_vsec=False, record_isec=False): """ Record current and voltage from all sections Parameters ---------- record_vsec : 'all' | 'soma' | False Option to record voltages from all sections ('all'), or just the soma ('soma'). Default: False. record_isec : 'all' | 'soma' | False Option to record voltages from all sections ('all'), or just the soma ('soma'). Default: False. """ section_names = list(self.sections.keys()) # Logic checks if just recording soma, sections, or both if record_vsec == 'soma': self.vsec = dict.fromkeys(['soma']) elif record_vsec == 'all': self.vsec = dict.fromkeys(section_names) if record_vsec: for sec_name in self.vsec: self.vsec[sec_name] = h.Vector() self.vsec[sec_name].record( self._nrn_sections[sec_name](0.5)._ref_v) if record_isec == 'soma': self.isec = dict.fromkeys(['soma']) elif record_isec == 'all': self.isec = dict.fromkeys(section_names) if record_isec: for sec_name in self.isec: list_syn = [key for key in self._nrn_synapses.keys() if key.startswith(f'{sec_name}_')] self.isec[sec_name] = dict.fromkeys(list_syn) for syn_name in self.isec[sec_name]: self.isec[sec_name][syn_name] = h.Vector() self.isec[sec_name][syn_name].record( self._nrn_synapses[syn_name]._ref_i)
[docs] def syn_create(self, secloc, e, tau1, tau2): """Create an h.Exp2Syn synapse. Parameters ---------- secloc : instance of nrn.Segment The section location. E.g., soma(0.5). e: float Reverse potential (in mV) tau1: float Rise time (in ms) tau2: float Decay time (in ms) Returns ------- syn : instance of h.Exp2Syn A two state kinetic scheme synapse. """ if not isinstance(secloc, nrn.Segment): raise TypeError(f'secloc must be instance of' f'nrn.Segment. Got {type(secloc)}') syn = h.Exp2Syn(secloc) syn.e = e syn.tau1 = tau1 syn.tau2 = tau2 return syn
[docs] def setup_source_netcon(self, threshold): """Created for _PC.cell and specifies SOURCES of spikes. Parameters ---------- threshold : float The voltage threshold for action potential. """ nc = h.NetCon(self._nrn_sections['soma'](0.5)._ref_v, None, sec=self._nrn_sections['soma']) nc.threshold = threshold return nc
[docs] def parconnect_from_src(self, gid_presyn, nc_dict, postsyn, inplane_distance): """Parallel receptor-centric connect FROM presyn TO this cell, based on GID. Parameters ---------- gid_presyn : int The cell ID of the presynaptic neuron nc_dict : dict Dictionary with keys: pos_src, A_weight, A_delay, lamtha Defines the connection parameters postsyn : instance of h.Exp2Syn The postsynaptic cell object. inplane_distance : float The in plane-distance (in um) between pyramidal cell somas in the square grid. Returns ------- nc : instance of h.NetCon A network connection object. """ from .network_builder import _PC nc = _PC.gid_connect(gid_presyn, postsyn) # set props here. nc.threshold = nc_dict['threshold'] nc.weight[0], nc.delay = _get_gaussian_connection( nc_dict['pos_src'], self.pos, nc_dict, inplane_distance=inplane_distance) return nc
[docs] def plot_morphology(self, ax=None, cell_types=None, show=True): """Plot the cell morphology. Parameters ---------- ax : instance of Axes3D Matplotlib 3D axis show : bool If True, show the plot Returns ------- axes : instance of Axes3D The matplotlib 3D axis handle. """ return plot_cell_morphology(self, ax=ax, show=show)
def _update_end_pts(self): """"Create cell and copy coordinates to Section.end_pts""" self._create_sections(self.sections, self.topology) section_names = list(self.sections.keys()) for name in section_names: nrn_pts = self._nrn_sections[name].psection()['morphology'][ 'pts3d'] del self._nrn_sections[name] x0, y0, z0 = nrn_pts[0][0], nrn_pts[0][1], nrn_pts[0][2] x1, y1, z1 = nrn_pts[1][0], nrn_pts[1][1], nrn_pts[1][2] self.sections[name]._end_pts = [[x0, y0, z0], [x1, y1, z1]] self._nrn_sections = dict()
[docs] def modify_section(self, sec_name, L=None, diam=None, cm=None, Ra=None): """Change attributes of section specified by `sec_name` Parameters ---------- sec_name : str Name of section to be modified. Must be a key of Cell.sections L : float | int | None length of a section in microns. Default None. diam : float | int | None diameter of a section in microns. cm : float | int | None membrane capacitance in micro-Farads. Ra : float | int | None axial resistivity in ohm-cm. Notes ----- Leaving default of None produces no change. """ valid_sec_names = list(self.sections.keys()) _check_option('sec_name', sec_name, valid_sec_names) if L is not None: _validate_type(L, (float, int), 'L') self.sections[sec_name]._L = L if diam is not None: _validate_type(diam, (float, int), 'diam') self.sections[sec_name]._diam = diam if cm is not None: _validate_type(cm, (float, int), 'cm') self.sections[sec_name]._cm = cm if Ra is not None: _validate_type(Ra, (float, int), 'Ra') self.sections[sec_name]._Ra = Ra self._update_end_pts()