#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Copyright (C) 2012 Computational Neuroscience Group, NMBU.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program 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 General Public License for more details.
"""
from __future__ import division
import os
import neuron
import numpy as np
import scipy.stats
import sys
import posixpath
from warnings import warn
import pickle
from .run_simulation import _run_simulation, _run_simulation_with_electrode
from .run_simulation import _collect_geometry_neuron
from .alias_method import alias_method
# check neuron version:
try:
try:
assert(neuron.version >= '7.6.4')
except:
warn('LFPy could not read NEURON version info. v7.6.4 or newer required')
except AssertionError:
warn('LFPy requires NEURON v7.6.4 or newer. Found v{}'.format(neuron.version))
[docs]class Cell(object):
"""
The main cell class used in LFPy.
Parameters
----------
morphology : str or neuron.h.SectionList
File path of morphology on format that NEURON can understand (w. file
ending .hoc, .asc, .swc or .xml), or neuron.h.SectionList instance
filled with references to neuron.h.Section instances.
v_init : float
Initial membrane potential. Defaults to -70 mV.
Ra : float
Axial resistance. Defaults to 35.4 Ohm*cm
cm : float
Membrane capacitance. Defaults to 1.0 uF/cm2.
passive : bool
Passive mechanisms are initialized if True. Defaults to False
passive_parameters : dict
parameter dictionary with values for the passive membrane mechanism in
NEURON ('pas'). The dictionary must contain keys 'g_pas' [S/cm^2] and
'e_pas' [mV], like the default:
passive_parameters=dict(g_pas=0.001, e_pas=-70)
extracellular : bool
Switch for NEURON's extracellular mechanism. Defaults to False
dt : float
simulation timestep. Defaults to 2^-4 ms
tstart : float
Initialization time for simulation <= 0 ms. Defaults to 0.
tstop : float
Stop time for simulation > 0 ms. Defaults to 100 ms.
nsegs_method : 'lambda100' or 'lambda_f' or 'fixed_length' or None
nseg rule, used by NEURON to determine number of compartments.
Defaults to 'lambda100'
max_nsegs_length : float or None
Maximum segment length for method 'fixed_length'. Defaults to None
lambda_f : int
AC frequency for method 'lambda_f'. Defaults to 100
d_lambda : float
Parameter for d_lambda rule. Defaults to 0.1
delete_sections : bool
Delete pre-existing section-references. Defaults to True
custom_code : list or None
List of model-specific code files ([.py/.hoc]). Defaults to None
custom_fun : list or None
List of model-specific functions with args. Defaults to None
custom_fun_args : list or None
List of args passed to custom_fun functions. Defaults to None
pt3d : bool
Use pt3d-info of the cell geometries switch. Defaults to False
celsius : float or None
Temperature in celsius. If nothing is specified here
or in custom code it is 6.3 celcius
verbose : bool
Verbose output switch. Defaults to False
Examples
--------
Simple example of how to use the Cell class with a passive-circuit
morphology (modify morphology path accordingly):
>>> import os
>>> import LFPy
>>> cellParameters = {
>>> 'morphology' : os.path.join('examples', 'morphologies', 'L5_Mainen96_LFPy.hoc'),
>>> 'v_init' : -65.,
>>> 'cm' : 1.0,
>>> 'Ra' : 150,
>>> 'passive' : True,
>>> 'passive_parameters' : {'g_pas' : 1./30000, 'e_pas' : -65},
>>> 'dt' : 2**-3,
>>> 'tstart' : 0,
>>> 'tstop' : 50,
>>> }
>>> cell = LFPy.Cell(**cellParameters)
>>> cell.simulate()
>>> print(cell.somav)
"""
def __init__(self, morphology,
v_init=-70.,
Ra=35.4,
cm=1.0,
passive=False,
passive_parameters = dict(
g_pas=0.001,
e_pas=-70.),
extracellular=False,
tstart=0.,
tstop=100.,
dt = 2**-4,
nsegs_method='lambda100',
lambda_f = 100,
d_lambda = 0.1,
max_nsegs_length=None,
delete_sections = True,
custom_code=None,
custom_fun=None,
custom_fun_args=None,
pt3d=False,
celsius=None,
verbose=False,
**kwargs):
"""
Initialization of the Cell object.
"""
self.verbose = verbose
self.pt3d = pt3d
# raise Exceptions on deprecated input arguments
for key in ['timeres_NEURON', 'timeres_python']:
if key in kwargs.keys():
raise DeprecationWarning('cell parameter {} is deprecated. Use dt=float instead'.format(key))
if 'tstartms' in kwargs.keys():
raise DeprecationWarning('cell parameter tstartms is deprecated. Use tstart=float instead')
if 'tstopms' in kwargs.keys():
raise DeprecationWarning('cell parameter tstopms is deprecated. Use tstop=float instead')
if 'rm' in kwargs.keys():
raise DeprecationWarning('Cell parameter rm is deprecated, set parameter passive_parameters=dict(g_pas=1/rm, e_pas=e_pas) instead')
if 'e_pas' in kwargs.keys():
raise DeprecationWarning('Cell parameter e_pas is deprecated, set parameter passive_parameters=dict(g_pas=1/rm, e_pas=e_pas) instead')
# check if there are un-used keyword arguments present in kwargs
for key, value in kwargs.items():
raise ValueError('The keyword and argument {}={} is not valid input to class LFPy.Cell'.format(key, value))
if passive:
try:
assert(type(passive_parameters) is dict)
except AssertionError:
raise AssertionError('passive_parameters must be a dictionary')
for key in ['g_pas', 'e_pas']:
try:
assert(key in passive_parameters.keys())
except AssertionError:
raise AssertionError('key {} not found in passive_parameters'.format(key))
if not hasattr(neuron.h, 'd_lambda'):
neuron.h.load_file('stdlib.hoc') #NEURON std. library
neuron.h.load_file('import3d.hoc') #import 3D morphology lib
if delete_sections:
numsec = 0
for numsec, sec in enumerate(neuron.h.allsec()):
pass
if numsec > 0 and self.verbose:
print('%s existing sections deleted from memory' % numsec)
neuron.h('forall delete_section()')
#print a warning if neuron have existing sections
numsec = 0
for numsec, sec in enumerate(neuron.h.allsec()):
pass
if numsec > 0 and self.verbose:
mssg = "%s sections detected! " % numsec + \
"Consider setting 'delete_sections=True'"
warn(mssg)
#load morphology
try:
assert(morphology is not None)
except AssertionError:
raise AssertionError('deprecated keyword argument morphology==None, value must be a file path or neuron.h.SectionList instance with neuron.h.Section instances')
if "win32" in sys.platform and type(morphology) is str:
# fix Path on windows
morphology = morphology.replace(os.sep, posixpath.sep)
self.morphology = morphology
if type(self.morphology) is str:
if os.path.isfile(self.morphology):
self._load_geometry()
else:
raise Exception('non-existent file %s' % self.morphology)
else:
try:
assert(type(self.morphology) is type(neuron.h.SectionList))
# #will try to import top level cell and create sectionlist,
# #in case there were no morphology file loaded
# neuron.h.define_shape()
# self._create_sectionlists()
except AssertionError:
raise Exception("Could not recognize Cell keyword argument morphology as neuron.h.SectionList instance")
# instantiate 3D geometry of all sections
neuron.h.define_shape()
# set some additional attributes
self._create_sectionlists_from_morphology_value()
#Some parameters and lists initialised
try:
assert(tstart <= 0)
except AssertionError:
raise AssertionError('tstart must be <= 0.')
try:
assert(dt in 2.**np.arange(-16, -1))
except AssertionError:
if tstart == 0.:
if self.verbose:
print('int(1./dt) not factorizable in base 2. '
'cell.tvec errors may occur, continuing initialization.')
elif tstart < 0:
raise AssertionError('int(1./dt) must be factorizable in base 2 if tstart < 0.')
self.dt = dt
self.tstart = tstart
self.tstop = tstop
self.synapses = []
self.synidx = []
self.pointprocesses = []
self.pointprocess_idx = []
self.v_init = v_init
self.default_rotation = self._get_rotation()
if passive:
#Set passive properties, insert passive on all segments
self.Ra = Ra
self.cm = cm
self.passive_parameters = passive_parameters
self._set_passive()
else:
self.Ra = Ra
self.cm = cm
self.passive_parameters = passive_parameters
if self.verbose:
print('No passive properties added')
#run user specified code and functions if argument given
if custom_code is not None or custom_fun is not None:
self._run_custom_codes(custom_code, custom_fun, custom_fun_args)
#Insert extracellular mech on all segments
self.extracellular = extracellular
if self.extracellular:
self._set_extracellular()
else:
if self.verbose:
print("no extracellular mechanism inserted")
#set number of segments accd to rule, and calculate the number
self._set_nsegs(nsegs_method, lambda_f, d_lambda, max_nsegs_length)
self.totnsegs = self._calc_totnsegs()
if self.verbose:
print("Total number of segments: %i" % self.totnsegs)
#extract pt3d info from NEURON, and set these with the same rotation
#and position in space as in our simulations, assuming RH rule, which
#NEURON do NOT use in shape plot
if self.pt3d:
self.x3d, self.y3d, self.z3d, self.diam3d = self._collect_pt3d()
#Gather geometry, set position and rotation of morphology
if self.pt3d:
self._update_pt3d()
else: # self._update_pt3d itself makes a call to self._collect_geometry()
self._collect_geometry()
if hasattr(self, 'somapos'):
self.set_pos()
else:
if self.verbose:
print('no soma, using the midpoint if initial segment.')
self.set_rotation(**self.default_rotation)
if celsius is not None:
if neuron.h.celsius != 6.3:
print("Overwriting custom temperature of %1.2f. New temperature is %1.2f"
% (neuron.h.celsius, celsius))
neuron.h.celsius = celsius
# initialize membrane voltage in all segments.
neuron.h.finitialize(self.v_init)
def _load_geometry(self):
"""Load the morphology-file in NEURON"""
try:
neuron.h.sec_counted = 0
except LookupError:
neuron.h('sec_counted = 0')
#import the morphology, try and determine format
fileEnding = self.morphology.split('.')[-1]
if fileEnding == 'hoc' or fileEnding == 'HOC':
neuron.h.load_file(1, self.morphology)
else:
neuron.h('objref this')
if fileEnding == 'asc' or fileEnding == 'ASC':
Import = neuron.h.Import3d_Neurolucida3()
if not self.verbose:
Import.quiet = 1
elif fileEnding == 'swc' or fileEnding == 'SWC':
Import = neuron.h.Import3d_SWC_read()
elif fileEnding == 'xml' or fileEnding == 'XML':
Import = neuron.h.Import3d_MorphML()
else:
raise ValueError('%s is not a recognised morphology file format!'
).with_traceback(
'Should be either .hoc, .asc, .swc, .xml!' % self.morphology)
#assuming now that morphologies file is the correct format
try:
Import.input(self.morphology)
except:
if not hasattr(neuron, 'neuroml'):
raise Exception('Can not import, try and copy the ' +
'nrn/share/lib/python/neuron/neuroml ' +
'folder into %s' % neuron.__path__[0])
else:
raise Exception('something wrong with file, see output')
try:
imprt = neuron.h.Import3d_GUI(Import, 0)
except:
raise Exception('See output, try to correct the file')
imprt.instantiate(neuron.h.this)
neuron.h.define_shape()
self._create_sectionlists()
def _run_custom_codes(self, custom_code, custom_fun, custom_fun_args):
"""Execute custom model code and functions with arguments"""
# load custom codes
if custom_code is not None:
for code in custom_code:
if "win32" in sys.platform:
code = code.replace(os.sep, posixpath.sep)
if code.split('.')[-1] == 'hoc':
try:
neuron.h.xopen(code)
except RuntimeError:
ERRMSG = '\n'.join(['',
'Could not load custom model code (%s)' %code,
'while creating a Cell object.',
'One possible cause is the NEURON mechanisms have',
'not been compiled, ',
'try running nrnivmodl or mknrndll (Windows) in ',
'the mod-file-containing folder. ',])
raise Exception(ERRMSG)
elif code.split('.')[-1] == 'py':
if sys.version >= "3.4":
exec(code, globals())
else:
exec(code)
else:
raise Exception('%s not a .hoc- nor .py-file' % code)
# run custom functions with arguments
i = 0
if custom_fun is not None:
for fun in custom_fun:
fun(**custom_fun_args[i])
i += 1
#recreate sectionlists in case something changed
neuron.h.define_shape()
self._create_sectionlists()
def _set_nsegs(self, nsegs_method, lambda_f, d_lambda, max_nsegs_length):
"""Set number of segments per section according to the lambda-rule,
or according to maximum length of segments"""
if nsegs_method == 'lambda100':
self._set_nsegs_lambda100(d_lambda)
elif nsegs_method == 'lambda_f':
self._set_nsegs_lambda_f(lambda_f, d_lambda)
elif nsegs_method == 'fixed_length':
self._set_nsegs_fixed_length(max_nsegs_length)
else:
if self.verbose:
print('No nsegs_method applied (%s)' % nsegs_method)
def _get_rotation(self):
"""Check if there exists a corresponding file
with rotation angles"""
if type(self.morphology) is str:
base = os.path.splitext(self.morphology)[0]
if os.path.isfile(base+'.rot'):
rotation_file = base+'.rot'
rotation_data = open(rotation_file)
rotation = {}
for line in rotation_data:
var, val = line.split('=')
val = val.strip()
val = float(str(val))
rotation[var] = val
else:
rotation = {}
else:
rotation = {}
return rotation
def _create_sectionlists(self):
"""Create section lists for different kinds of sections"""
#list with all sections
self.allsecnames = []
self.allseclist = neuron.h.SectionList()
for sec in neuron.h.allsec():
self.allsecnames.append(sec.name())
self.allseclist.append(sec=sec)
#list of soma sections, assuming it is named on the format "soma*"
self.nsomasec = 0
self.somalist = neuron.h.SectionList()
for sec in neuron.h.allsec():
if sec.name().find('soma') >= 0:
self.somalist.append(sec=sec)
self.nsomasec += 1
def _create_sectionlists_from_morphology_value(self):
"""Variant of Cell._create_sectionlists() used if keyword argument
morphology is a neuron.h.SectionList instance"""
#list with all sections
self.allsecnames = []
self.allseclist = self.morphology
for sec in self.allseclist:
self.allsecnames.append(sec.name())
#list of soma sections, assuming it is named on the format "soma*"
self.nsomasec = 0
self.somalist = neuron.h.SectionList()
for sec in self.allseclist:
if sec.name().find('soma') >= 0:
self.somalist.append(sec=sec)
self.nsomasec += 1
def _get_idx(self, seclist):
"""Return boolean vector which indexes where segments in seclist
matches segments in neuron.h.allsec(), rewritten from
LFPy.hoc function get_idx()"""
if neuron.h.allsec() == seclist:
return np.ones(self.totnsegs, dtype=bool)
else:
idxvec = np.zeros(self.totnsegs, dtype=bool)
#get sectionnames from seclist
seclistnames = []
for sec in seclist:
seclistnames.append(sec.name())
seclistnames = np.array(seclistnames, dtype='|S128')
segnames = np.empty(self.totnsegs, dtype='|S128')
i = 0
for sec in self.allseclist:
secname = sec.name()
for seg in sec:
segnames[i] = secname
i += 1
for name in seclistnames:
idxvec[segnames == name] = True
return idxvec
def _set_nsegs_lambda_f(self, frequency=100, d_lambda=0.1):
"""Set the number of segments for section according to the
d_lambda-rule for a given input frequency
Parameters
----------
frequency : float
frequency at which AC length constant is computed
d_lambda : float
"""
for sec in self.allseclist:
sec.nseg = int((sec.L / (d_lambda*neuron.h.lambda_f(frequency,
sec=sec)) + .9) / 2)*2 + 1
if self.verbose:
print("set nsegs using lambda-rule with frequency %i." % frequency)
def _set_nsegs_lambda100(self, d_lambda=0.1):
"""Set the numbers of segments using d_lambda(100)"""
self._set_nsegs_lambda_f(frequency=100, d_lambda=d_lambda)
def _set_nsegs_fixed_length(self, maxlength):
"""Set nseg for sections so that every segment L < maxlength"""
for sec in self.allseclist:
sec.nseg = int(sec.L / maxlength) + 1
def _calc_totnsegs(self):
"""Calculate the number of segments in the allseclist"""
i = 0
for sec in self.allseclist:
i += sec.nseg
return i
def _check_currents(self):
"""Check that the sum of all membrane and electrode currents over all
segments is sufficiently close to zero"""
raise NotImplementedError('this function need to be written')
def _set_passive(self):
"""Insert passive mechanism on all segments"""
for sec in self.allseclist:
sec.insert('pas')
sec.Ra = self.Ra
sec.cm = self.cm
sec.g_pas = self.passive_parameters['g_pas']
sec.e_pas = self.passive_parameters['e_pas']
def _set_extracellular(self):
"""Insert extracellular mechanism on all sections
to set an external potential V_ext as boundary condition.
"""
for sec in self.allseclist:
sec.insert('extracellular')
[docs] def set_synapse(self, idx, syntype,
record_current=False,
record_potential=False,
weight=None, **kwargs):
"""Insert synapse on cell segment
Parameters
----------
idx : int
Index of compartment where synapse is inserted
syntype : str
Type of synapse. Built-in types in NEURON: ExpSyn, Exp2Syn
record_current : bool
If True, record synapse current
record_potential : bool
If True, record postsynaptic potential seen by the synapse
weight : float
Strength of synapse
kwargs
arguments passed on from class Synapse
"""
if not hasattr(self, 'synlist'):
self.synlist = neuron.h.List()
if not hasattr(self, 'synireclist'):
self.synireclist = neuron.h.List()
if not hasattr(self, 'synvreclist'):
self.synvreclist = neuron.h.List()
if not hasattr(self, 'netstimlist'):
self.netstimlist = neuron.h.List()
if not hasattr(self, 'netconlist'):
self.netconlist = neuron.h.List()
if not hasattr(self, 'sptimeslist'):
self.sptimeslist = neuron.h.List()
i = 0
cmd = 'syn = neuron.h.{}(seg.x, sec=sec)'
for sec in self.allseclist:
for seg in sec:
if i == idx:
command = cmd.format(syntype)
if sys.version >= "3.4":
exec(command, locals(), globals())
else:
exec(command)
for param in list(kwargs.keys()):
try:
if sys.version >= "3.4":
exec('syn.' + param + '=' + str(kwargs[param]),
locals(), globals())
else:
exec('syn.' + param + '=' + str(kwargs[param]))
except:
pass
self.synlist.append(syn)
#create NetStim (generator) and NetCon (connection) objects
self.netstimlist.append(neuron.h.NetStim(0.5))
self.netstimlist[-1].number = 0
nc = neuron.h.NetCon(self.netstimlist[-1], syn)
nc.weight[0] = weight
self.netconlist.append(nc)
#record currents
if record_current:
synirec = neuron.h.Vector(int(self.tstop //
self.dt+1))
synirec.record(syn._ref_i, self.dt)
self.synireclist.append(synirec)
else:
synirec = neuron.h.Vector(0)
self.synireclist.append(synirec)
#record potential
if record_potential:
synvrec = neuron.h.Vector(int(self.tstop //
self.dt+1))
synvrec.record(seg._ref_v, self.dt)
self.synvreclist.append(synvrec)
else:
synvrec = neuron.h.Vector(0)
self.synvreclist.append(synvrec)
i += 1
return self.synlist.count() - 1
[docs] def set_point_process(self, idx, pptype, record_current=False,
record_potential=False, **kwargs):
"""Insert pptype-electrode type pointprocess on segment numbered
idx on cell object
Parameters
----------
idx : int
Index of compartment where point process is inserted
pptype : str
Type of pointprocess. Examples: SEClamp, VClamp,
IClamp, SinIClamp, ChirpIClamp
record_current : bool
Decides if current is stored
kwargs
Parameters passed on from class StimIntElectrode
"""
if not hasattr(self, 'stimlist'):
self.stimlist = neuron.h.List()
if not hasattr(self, 'stimireclist'):
self.stimireclist = neuron.h.List()
if not hasattr(self, 'stimvreclist'):
self.stimvreclist = neuron.h.List()
i = 0
cmd1 = 'stim = neuron.h.'
cmd2 = '(seg.x, sec=sec)'
for sec in self.allseclist:
for seg in sec:
if i == idx:
command = cmd1 + pptype + cmd2
if sys.version >= "3.4":
exec(command, locals(), globals())
else:
exec(command)
for param in list(kwargs.keys()):
try:
if sys.version >= "3.4":
exec('stim.' + param + '=' + str(kwargs[param]),
locals(), globals())
else:
exec('stim.' + param + '=' + str(kwargs[param]))
except SyntaxError:
ERRMSG = ''.join(['',
'Point process type "{0}" might not '.format(
pptype),
'recognize attribute "{0}". '.format(param),
'Check for misspellings'])
raise Exception(ERRMSG)
self.stimlist.append(stim)
#record current
if record_current:
stimirec = neuron.h.Vector(int(self.tstop /
self.dt+1))
stimirec.record(stim._ref_i, self.dt)
self.stimireclist.append(stimirec)
else:
stimirec = neuron.h.Vector(0)
self.stimireclist.append(stimirec)
# record potential
if record_potential:
stimvrec = neuron.h.Vector(int(self.tstop /
self.dt+1))
stimvrec.record(seg._ref_v, self.dt)
self.stimvreclist.append(stimvrec)
else:
stimvrec = neuron.h.Vector(0)
self.stimvreclist.append(stimvrec)
i += 1
return self.stimlist.count() - 1
def _collect_geometry(self):
"""Collects x, y, z-coordinates from NEURON"""
#None-type some attributes if they do not exis:
if not hasattr(self, 'xstart'):
self.xstart = None
self.ystart = None
self.zstart = None
self.xend = None
self.yend = None
self.zend = None
self.area = None
self.diam = None
self.length = None
_collect_geometry_neuron(self)
self._calc_midpoints()
self.somaidx = self.get_idx(section='soma')
if self.somaidx.size > 1:
xmids = self.xmid[self.somaidx]
ymids = self.ymid[self.somaidx]
zmids = self.zmid[self.somaidx]
self.somapos = np.zeros(3)
self.somapos[0] = xmids.mean()
self.somapos[1] = ymids.mean()
self.somapos[2] = zmids.mean()
elif self.somaidx.size == 1:
self.somapos = np.zeros(3)
self.somapos[0] = self.xmid[self.somaidx]
self.somapos[1] = self.ymid[self.somaidx]
self.somapos[2] = self.zmid[self.somaidx]
elif self.somaidx.size == 0:
if self.verbose:
print('There is no soma!')
print('using first segment as root point')
self.somaidx = np.array([0])
self.somapos = np.zeros(3)
self.somapos[0] = self.xmid[self.somaidx]
self.somapos[1] = self.ymid[self.somaidx]
self.somapos[2] = self.zmid[self.somaidx]
else:
raise Exception('Huh?!')
def _calc_midpoints(self):
"""Calculate midpoints of each segment"""
self.xmid = .5*(self.xstart+self.xend).flatten()
self.ymid = .5*(self.ystart+self.yend).flatten()
self.zmid = .5*(self.zstart+self.zend).flatten()
[docs] def get_idx(self, section='allsec', z_min=-np.inf, z_max=np.inf):
"""Returns compartment idx of segments from sections with names that match
the pattern defined in input section on interval [z_min, z_max].
Parameters
----------
section : str
Any entry in cell.allsecnames or just 'allsec'.
z_min : float
Depth filter. Specify minimum z-position
z_max : float
Depth filter. Specify maximum z-position
Examples
--------
>>> idx = cell.get_idx(section='allsec')
>>> print(idx)
>>> idx = cell.get_idx(section=['soma', 'dend', 'apic'])
>>> print(idx)
"""
if section == 'allsec':
seclist = neuron.h.allsec()
else:
seclist = neuron.h.SectionList()
if type(section) == str:
for sec in self.allseclist:
if sec.name().find(section) >= 0:
seclist.append(sec=sec)
elif type(section) == list:
for secname in section:
for sec in self.allseclist:
if sec.name().find(secname) >= 0:
seclist.append(sec=sec)
else:
if self.verbose:
print('%s did not match any section name' % str(section))
idx = self._get_idx(seclist)
sel_z_idx = (self.zmid[idx] > z_min) & (self.zmid[idx] < z_max)
return np.arange(self.totnsegs)[idx][sel_z_idx]
[docs] def get_closest_idx(self, x=0., y=0., z=0., section='allsec'):
"""Get the index number of a segment in specified section which
midpoint is closest to the coordinates defined by the user
Parameters
----------
x: float
x-coordinate
y: float
y-coordinate
z: float
z-coordinate
section: str
String matching a section-name. Defaults to 'allsec'.
"""
idx = self.get_idx(section)
dist = ((self.xmid[idx] - x)**2 +
(self.ymid[idx] - y)**2 +
(self.zmid[idx] - z)**2)
return idx[np.argmin(dist)]
[docs] def get_rand_idx_area_norm(self, section='allsec', nidx=1,
z_min=-1E6, z_max=1E6):
"""Return nidx segment indices in section with random probability
normalized to the membrane area of segment on
interval [z_min, z_max]
Parameters
----------
section : str
String matching a section-name
nidx : int
Number of random indices
z_min : float
Depth filter
z_max : float
Depth filter
"""
poss_idx = self.get_idx(section=section, z_min=z_min, z_max=z_max)
if nidx < 1:
print('nidx < 1, returning empty array')
return np.array([])
elif poss_idx.size == 0:
print('No possible segment idx match enquire! returning empty array')
return np.array([])
else:
area = self.area[poss_idx]
area /= area.sum()
return alias_method(poss_idx, area, nidx)
[docs] def get_rand_idx_area_and_distribution_norm(self, section='allsec', nidx=1,
z_min=-1E6, z_max=1E6,
fun=scipy.stats.norm,
funargs=dict(loc=0, scale=100),
funweights=None):
"""
Return nidx segment indices in section with random probability
normalized to the membrane area of each segment multiplied by
the value of the probability density function of "fun", a function
in the scipy.stats module with corresponding function arguments
in "funargs" on the interval [z_min, z_max]
Parameters
----------
section: str
string matching a section-name
nidx: int
number of random indices
z_min: float
depth filter
z_max: float
depth filter
fun : function or str, or iterable of function or str
if function a scipy.stats method, if str, must be method in
scipy.stats module with the same name (like 'norm'),
if iterable (list, tuple, numpy.array) of function or str some
probability distribution in scipy.stats module
funargs : dict or iterable
iterable (list, tuple, numpy.array) of dict, arguments to fun.pdf
method (e.g., w. keys 'loc' and 'scale')
funweights : None or iterable
iterable (list, tuple, numpy.array) of floats, scaling of each
individual fun (i.e., introduces layer specificity)
Examples
--------
>>> import LFPy
>>> import numpy as np
>>> import scipy.stats as ss
>>> import matplotlib.pyplot as plt
>>> from os.path import join
>>> cell = LFPy.Cell(morphology=join('cells', 'cells', 'j4a.hoc'))
>>> cell.set_rotation(x=4.99, y=-4.33, z=3.14)
>>> idx = cell.get_rand_idx_area_and_distribution_norm(nidx=10000,
fun=ss.norm,
funargs=dict(loc=0, scale=200))
>>> bins = np.arange(-30, 120)*10
>>> plt.hist(cell.zmid[idx], bins=bins, alpha=0.5)
>>> plt.show()
"""
poss_idx = self.get_idx(section=section, z_min=z_min, z_max=z_max)
if nidx < 1:
print('nidx < 1, returning empty array')
return np.array([])
elif poss_idx.size == 0:
print('No possible segment idx match enquire! returning empty array')
return np.array([])
else:
p = self.area[poss_idx]
# scale with density function
if type(fun) in [list, tuple, np.ndarray]:
assert(type(funargs) in [list, tuple, np.ndarray])
assert(type(funweights) in [list, tuple, np.ndarray])
assert((len(fun) == len(funargs)) & (len(fun) == len(funweights)))
mod = np.zeros(poss_idx.shape)
for f, args, scl in zip(fun, funargs, funweights):
if type(f) is str and f in dir(scipy.stats):
f = getattr(scipy.stats, f)
df = f(**args)
mod += df.pdf(x=self.zmid[poss_idx])*scl
p *= mod
else:
if type(fun) is str and fun in dir(scipy.stats):
fun = getattr(scipy.stats, fun)
df = fun(**funargs)
p *= df.pdf(x=self.zmid[poss_idx])
# normalize
p /= p.sum()
return alias_method(poss_idx, p, nidx)
[docs] def simulate(self, electrode=None, rec_imem=False, rec_vmem=False,
rec_ipas=False, rec_icap=False,
rec_current_dipole_moment=False,
rec_variables=[], variable_dt=False, atol=0.001,
to_memory=True, to_file=False, file_name=None,
dotprodcoeffs=None, **kwargs):
"""
This is the main function running the simulation of the NEURON model.
Start NEURON simulation and record variables specified by arguments.
Parameters
----------
electrode : :obj: or list, optional
Either an LFPy.RecExtElectrode object or a list of such.
If supplied, LFPs will be calculated at every time step
and accessible as `electrode.LFP`. If a list of objects
is given, accessible as `electrode[0].LFP` etc.
rec_imem : bool
If true, segment membrane currents will be recorded
If no electrode argument is given, it is necessary to
set rec_imem=True in order to calculate LFP later on.
Units of (nA).
rec_vmem : bool
Record segment membrane voltages (mV)
rec_ipas : bool
Record passive segment membrane currents (nA)
rec_icap : bool
Record capacitive segment membrane currents (nA)
rec_current_dipole_moment : bool
If True, compute and record current-dipole moment from
transmembrane currents as in Linden et al. (2010) J Comput Neurosci,
DOI: 10.1007/s10827-010-0245-4. Will set the `LFPy.Cell` attribute
`current_dipole_moment` as n_timesteps x 3 `np.ndarray` where the
last dimension contains the x,y,z components of the dipole moment.
rec_variables : list
List of variables to record, i.e arg=['cai', ]
variable_dt : bool
Use variable timestep in NEURON
atol : float
Absolute tolerance used with NEURON variable timestep
to_memory : bool
Only valid with electrode, store lfp in -> electrode.LFP
to_file : bool
Only valid with electrode, save LFPs in hdf5 file format
file_name : str
Name of hdf5 file, '.h5' is appended if it doesnt exist
dotprodcoeffs : list
List of N x Nseg ndarray. These arrays will at
every timestep be multiplied by the membrane currents.
Presumably useful for memory efficient csd or lfp calcs
"""
for key in kwargs.keys():
if key in ['rec_isyn', 'rec_vmemsyn', 'rec_istim', 'rec_vmemstim']:
raise DeprecationWarning('Cell.simulate parameter {} is deprecated.'.format(key))
self._set_soma_volt_recorder()
self._collect_tvec()
# set up integrator, use the CVode().fast_imem method by default
# as it doesn't hurt sim speeds much if at all.
cvode = neuron.h.CVode()
try:
cvode.use_fast_imem(1)
except AttributeError as ae:
raise Exception('neuron.h.CVode().use_fast_imem() method not found. Please update NEURON to v.7.4 or newer')
if rec_imem:
self._set_imem_recorders()
if rec_vmem:
self._set_voltage_recorders()
if rec_ipas:
self._set_ipas_recorders()
if rec_icap:
self._set_icap_recorders()
if rec_current_dipole_moment:
self._set_current_dipole_moment_array()
if len(rec_variables) > 0:
self._set_variable_recorders(rec_variables)
#run fadvance until t >= tstop, and calculate LFP if asked for
if electrode is None and dotprodcoeffs is None and not rec_current_dipole_moment:
if not rec_imem and self.verbose:
print("rec_imem = %s, membrane currents will not be recorded!"
% str(rec_imem))
_run_simulation(self, cvode, variable_dt, atol)
else:
#allow using both electrode and additional coefficients:
_run_simulation_with_electrode(self, cvode, electrode, variable_dt, atol,
to_memory, to_file, file_name,
dotprodcoeffs,
rec_current_dipole_moment)
# somatic trace
if self.nsomasec >= 1:
self.somav = np.array(self.somav)
if rec_imem:
self._calc_imem()
if rec_ipas:
self._calc_ipas()
if rec_icap:
self._calc_icap()
if rec_vmem:
self._collect_vmem()
self._collect_isyn()
self._collect_vsyn()
self._collect_istim()
self._collect_vstim()
if len(rec_variables) > 0:
self._collect_rec_variables(rec_variables)
if hasattr(self, 'netstimlist'):
del self.netstimlist
def _collect_tvec(self):
"""
Set the tvec to be a monotonically increasing numpy array after sim.
"""
self.tvec = np.arange(self.tstop / self.dt + 1) * self.dt
def _calc_imem(self):
"""
Fetch the vectors from the memireclist and calculate self.imem
containing all the membrane currents.
"""
self.imem = np.array(self.memireclist)
self.memireclist = None
del self.memireclist
def _calc_ipas(self):
"""
Get the passive currents
"""
self.ipas = np.array(self.memipasreclist)
for i in range(self.ipas.shape[0]):
self.ipas[i, ] *= self.area[i] * 1E-2
self.memipasreclist = None
del self.memipasreclist
def _calc_icap(self):
"""
Get the capacitive currents
"""
self.icap = np.array(self.memicapreclist)
for i in range(self.icap.shape[0]):
self.icap[i, ] *= self.area[i] * 1E-2
self.memicapreclist = None
del self.memicapreclist
def _collect_vmem(self):
"""
Get the membrane currents
"""
self.vmem = np.array(self.memvreclist)
self.memvreclist = None
del self.memvreclist
def _collect_isyn(self):
"""
Get the synaptic currents
"""
for syn in self.synapses:
if syn.record_current:
syn.collect_current(self)
self.synireclist = None
del self.synireclist
def _collect_vsyn(self):
"""
Collect the membrane voltage of segments with synapses
"""
for syn in self.synapses:
if syn.record_potential:
syn.collect_potential(self)
self.synvreclist = None
del self.synvreclist
def _collect_istim(self):
"""
Get the pointprocess currents
"""
for pp in self.pointprocesses:
if pp.record_current:
pp.collect_current(self)
self.stimireclist = None
del self.stimireclist
def _collect_vstim(self):
"""
Collect the membrane voltage of segments with stimulus
"""
for pp in self.pointprocesses:
if pp.record_potential:
pp.collect_potential(self)
self.stimvreclist = None
del self.stimvreclist
def _collect_rec_variables(self, rec_variables):
"""
Create dict of np.arrays from recorded variables, each dictionary
element named as the corresponding recorded variable name, i.e 'cai'
"""
self.rec_variables = {}
i = 0
for values in self.recvariablesreclist:
self.rec_variables.update({rec_variables[i] : np.array(values)})
if self.verbose:
print('collected recorded variable %s' % rec_variables[i])
i += 1
del self.recvariablesreclist
def _loadspikes(self):
"""
Initialize spiketimes from netcon if they exist
"""
if hasattr(self, 'synlist'):
if len(self.synlist) == len(self.sptimeslist):
for i in range(int(self.synlist.count())):
for ii in range(int(self.sptimeslist.o(i).size)):
self.netconlist.o(i).event(float(self.sptimeslist.o(i)[ii]))
def _set_soma_volt_recorder(self):
"""Record somatic membrane potential"""
if self.nsomasec == 0:
if self.verbose:
warn('Cell instance appears to have no somatic section. '
'No somav attribute will be set.')
elif self.nsomasec == 1:
self.somav = neuron.h.Vector(int(self.tstop / self.dt+1))
for sec in self.somalist:
self.somav.record(sec(0.5)._ref_v, self.dt)
elif self.nsomasec > 1:
self.somav = neuron.h.Vector(int(self.tstop / self.dt+1))
nseg = self.get_idx('soma').size
i, j = divmod(nseg, 2)
k = 1
for sec in self.somalist:
for seg in sec:
if nseg==2 and k == 1:
#if 2 segments, record from the first one:
self.somav.record(seg._ref_v, self.dt)
else:
if k == i*2:
#record from one of the middle segments:
self.somav.record(seg._ref_v, self.dt)
k += 1
def _set_imem_recorders(self):
"""
Record membrane currents for all segments
"""
self.memireclist = neuron.h.List()
for sec in self.allseclist:
for seg in sec:
memirec = neuron.h.Vector(int(self.tstop / self.dt+1))
memirec.record(seg._ref_i_membrane_, self.dt)
self.memireclist.append(memirec)
def _set_ipas_recorders(self):
"""
Record passive membrane currents for all segments
"""
self.memipasreclist = neuron.h.List()
for sec in self.allseclist:
for seg in sec:
memipasrec = neuron.h.Vector(int(self.tstop / self.dt+1))
memipasrec.record(seg._ref_i_pas, self.dt)
self.memipasreclist.append(memipasrec)
def _set_icap_recorders(self):
"""
Record capacitive membrane currents for all segments
"""
self.memicapreclist = neuron.h.List()
for sec in self.allseclist:
for seg in sec:
memicaprec = neuron.h.Vector(int(self.tstop / self.dt+1))
memicaprec.record(seg._ref_i_cap, self.dt)
self.memicapreclist.append(memicaprec)
def _set_voltage_recorders(self):
"""
Record membrane potentials for all segments
"""
self.memvreclist = neuron.h.List()
for sec in self.allseclist:
for seg in sec:
memvrec = neuron.h.Vector(int(self.tstop / self.dt+1))
memvrec.record(seg._ref_v, self.dt)
self.memvreclist.append(memvrec)
def _set_current_dipole_moment_array(self):
"""
Creates container for current dipole moment, an empty
n_timesteps x 3 `numpy.ndarray` that will be filled with values during
the course of each simulation
"""
self.current_dipole_moment = np.zeros((self.tvec.size, 3))
def _set_variable_recorders(self, rec_variables):
"""
Create a recorder for each variable name in list
rec_variables
Variables is stored in nested list self.recvariablesreclist
"""
self.recvariablesreclist = neuron.h.List()
for variable in rec_variables:
variablereclist = neuron.h.List()
self.recvariablesreclist.append(variablereclist)
for sec in self.allseclist:
for seg in sec:
recvector = neuron.h.Vector(int(self.tstop / self.dt + 1))
try:
recvector.record(getattr(seg, '_ref_%s' % variable),
self.dt)
except(NameError, AttributeError):
print('non-existing variable %s, section %s.%f' %
(variable, sec.name(), seg.x))
variablereclist.append(recvector)
[docs] def set_pos(self, x=0., y=0., z=0.):
"""Set the cell position.
Move the cell geometry so that midpoint of soma section is
in (x, y, z). If no soma pos, use the first segment
Parameters
----------
x : float
x position defaults to 0.0
y : float
y position defaults to 0.0
z : float
z position defaults to 0.0
"""
diffx = x-self.somapos[0]
diffy = y-self.somapos[1]
diffz = z-self.somapos[2]
#also update the pt3d_pos:
if self.pt3d and hasattr(self, 'x3d'):
self._set_pt3d_pos(diffx, diffy, diffz)
else:
self.somapos[0] = x
self.somapos[1] = y
self.somapos[2] = z
self.xstart += diffx
self.ystart += diffy
self.zstart += diffz
self.xend += diffx
self.yend += diffy
self.zend += diffz
self._calc_midpoints()
self._update_synapse_positions()
[docs] def strip_hoc_objects(self):
"""Destroy any NEURON hoc objects in the cell object"""
for varname in dir(self):
if type(getattr(self, varname)) == type(neuron.h.List()):
setattr(self, varname, None)
if self.verbose:
print('None-typed %s in cell instance' % varname)
[docs] def cellpickler(self, filename, pickler=pickle.dump):
"""Save data in cell to filename, using cPickle. It will however destroy
any neuron.h objects upon saving, as c-objects cannot be pickled
Parameters
----------
filename : str
Where to save cell
Examples
--------
To save a cell, use:
>>> cell.cellpickler('cell.cpickle')
To load this cell again in another session:
>>> import cPickle
>>> f = file('cell.cpickle', 'rb')
>>> cell = cPickle.load(f)
>>> f.close()
alternatively:
>>> import LFPy
>>> cell = LFPy.tools.load('cell.cpickle')
"""
self.strip_hoc_objects()
if pickler==pickle.dump:
filen = open(filename, 'wb')
pickle.dump(self, filen, protocol=2)
filen.close()
return None
elif pickler==pickle.dumps:
return pickle.dumps(self)
def _update_synapse_positions(self):
"""
Update synapse positions after rotation of morphology
"""
for i in range(len(self.synapses)):
self.synapses[i].update_pos(self)
[docs] def set_rotation(self, x=None, y=None, z=None, rotation_order='xyz'):
"""
Rotate geometry of cell object around the x-, y-, z-axis in the order
described by rotation_order parameter.
rotation_order should be a string with 3 elements containing x, y, and z
e.g. 'xyz', 'zyx'
Input should be angles in radians.
using rotation matrices, takes dict with rot. angles,
where x, y, z are the rotation angles around respective axes.
All rotation angles are optional.
Examples
--------
>>> cell = LFPy.Cell(**kwargs)
>>> rotation = {'x' : 1.233, 'y' : 0.236, 'z' : np.pi}
>>> cell.set_rotation(**rotation)
"""
if type(rotation_order) is not str:
raise AttributeError('rotation_order must be a string')
elif 'x' not in rotation_order or 'y' not in rotation_order or 'z' not in rotation_order:
raise AttributeError("'x', 'y', and 'z' must be in rotation_order")
elif len(rotation_order) != 3:
raise AttributeError("rotation_order should have 3 elements (e.g. 'zyx')")
for ax in rotation_order:
if ax == 'x' and x is not None:
theta = -x
rotation_x = np.array([[1, 0, 0],
[0, np.cos(theta), -np.sin(theta)],
[0, np.sin(theta), np.cos(theta)]])
rel_start, rel_end = self._rel_positions()
rel_start = np.dot(rel_start, rotation_x)
rel_end = np.dot(rel_end, rotation_x)
self._real_positions(rel_start, rel_end)
if self.verbose:
print('Rotated geometry %g radians around x-axis' % (-theta))
else:
if self.verbose:
print('Geometry not rotated around x-axis')
if ax == 'y' and y is not None:
phi = -y
rotation_y = np.array([[np.cos(phi), 0, np.sin(phi)],
[0, 1, 0],
[-np.sin(phi), 0, np.cos(phi)]])
rel_start, rel_end = self._rel_positions()
rel_start = np.dot(rel_start, rotation_y)
rel_end = np.dot(rel_end, rotation_y)
self._real_positions(rel_start, rel_end)
if self.verbose:
print('Rotated geometry %g radians around y-axis' % (-phi))
else:
if self.verbose:
print('Geometry not rotated around y-axis')
if ax == 'z' and z is not None:
gamma = -z
rotation_z = np.array([[np.cos(gamma), -np.sin(gamma), 0],
[np.sin(gamma), np.cos(gamma), 0],
[0, 0, 1]])
rel_start, rel_end = self._rel_positions()
rel_start = np.dot(rel_start, rotation_z)
rel_end = np.dot(rel_end, rotation_z)
self._real_positions(rel_start, rel_end)
if self.verbose:
print('Rotated geometry %g radians around z-axis' % (-gamma))
else:
if self.verbose:
print('Geometry not rotated around z-axis')
#rotate the pt3d geometry accordingly
if self.pt3d and hasattr(self, 'x3d'):
self._set_pt3d_rotation(x, y, z, rotation_order)
[docs] def chiral_morphology(self, axis='x'):
"""
Mirror the morphology around given axis, (default x-axis),
useful to introduce more heterogeneouties in morphology shapes
Parameters
----------
axis : str
'x' or 'y' or 'z'
"""
#morphology relative to soma-position
rel_start, rel_end = self._rel_positions()
if axis == 'x':
rel_start[:, 0] = -rel_start[:, 0]
rel_end[:, 0] = -rel_end[:, 0]
elif axis == 'y':
rel_start[:, 1] = -rel_start[:, 1]
rel_end[:, 1] = -rel_end[:, 1]
elif axis == 'z':
rel_start[:, 2] = -rel_start[:, 2]
rel_end[:, 2] = -rel_end[:, 2]
else:
raise Exception("axis must be either 'x', 'y' or 'z'")
if self.verbose:
print('morphology mirrored across %s-axis' % axis)
#set the proper 3D positions
self._real_positions(rel_start, rel_end)
def _squeeze_me_macaroni(self):
"""
Reducing the dimensions of the morphology matrices from 3D->1D
"""
self.xstart = np.array(self.xstart).flatten()
self.xend = np.array(self.xend).flatten()
self.ystart = np.array(self.ystart).flatten()
self.yend = np.array(self.yend).flatten()
self.zstart = np.array(self.zstart).flatten()
self.zend = np.array(self.zend).flatten()
def _rel_positions(self):
"""
Morphology relative to soma position
"""
rel_start = np.array([self.xstart-self.somapos[0],
self.ystart-self.somapos[1],
self.zstart-self.somapos[2]]).T
rel_end = np.array([self.xend-self.somapos[0],
self.yend-self.somapos[1],
self.zend-self.somapos[2]]).T
return rel_start, rel_end
def _real_positions(self, rel_start, rel_end):
"""
Morphology coordinates relative to Origo
"""
self.xstart = rel_start[:, 0] + self.somapos[0]
self.ystart = rel_start[:, 1] + self.somapos[1]
self.zstart = rel_start[:, 2] + self.somapos[2]
self.xend = rel_end[:, 0] + self.somapos[0]
self.yend = rel_end[:, 1] + self.somapos[1]
self.zend = rel_end[:, 2] + self.somapos[2]
self._squeeze_me_macaroni()
self._calc_midpoints()
self._update_synapse_positions()
[docs] def get_rand_prob_area_norm(self, section='allsec',
z_min=-10000, z_max=10000):
"""
Return the probability (0-1) for synaptic coupling on segments
in section sum(prob)=1 over all segments in section.
Probability normalized by area.
Parameters
----------
section : str
string matching a section-name. Defaults to 'allsec'
z_min : float
depth filter
z_max : float
depth filter
"""
idx = self.get_idx(section=section, z_min=z_min, z_max = z_max)
prob = self.area[idx] / sum(self.area[idx])
return prob
[docs] def get_rand_prob_area_norm_from_idx(self, idx=np.array([0])):
"""
Return the normalized probability (0-1) for synaptic coupling on
segments in idx-array.
Normalised probability determined by area of segments.
Parameters
----------
idx : ndarray, dtype=int.
array of segment indices
"""
prob = self.area[idx] / sum(self.area[idx])
return prob
[docs] def get_intersegment_vector(self, idx0=0, idx1=0):
"""Return the distance between midpoints of two segments with index
idx0 and idx1. The argument returned is a vector [x, y, z], where
x = self.xmid[idx1] - self.xmid[idx0] etc.
Parameters
----------
idx0 : int
idx1 : int
"""
vector = []
try:
if idx1 < 0 or idx0 < 0:
raise Exception('idx0 < 0 or idx1 < 0')
vector.append(self.xmid[idx1] - self.xmid[idx0])
vector.append(self.ymid[idx1] - self.ymid[idx0])
vector.append(self.zmid[idx1] - self.zmid[idx0])
return vector
except:
ERRMSG = 'idx0 and idx1 must be ints on [0, %i]' % self.totnsegs
raise ValueError(ERRMSG)
[docs] def get_intersegment_distance(self, idx0=0, idx1=0):
"""
Return the Euclidean distance between midpoints of two segments.
Parameters
----------
idx0 : int
idx1 : int
Returns
-------
float
Will return a float in unit of micrometers.
"""
try:
vector = np.array(self.get_intersegment_vector(idx0, idx1))
return np.sqrt((vector**2).sum())
except:
ERRMSG = 'idx0 and idx1 must be ints on [0, %i]' % self.totnsegs
raise ValueError(ERRMSG)
[docs] def get_idx_children(self, parent="soma[0]"):
"""Get the idx of parent's children sections, i.e. compartments ids
of sections connected to parent-argument
Parameters
----------
parent : str
name-pattern matching a sectionname. Defaults to "soma[0]"
"""
idxvec = np.zeros(self.totnsegs)
secnamelist = []
childseclist = []
#filling list of sectionnames for all sections, one entry per segment
for sec in self.allseclist:
for seg in sec:
secnamelist.append(sec.name())
#filling list of children section-names
sref = neuron.h.SectionRef(parent)
for sec in sref.child:
childseclist.append(sec.name())
#idxvec=1 where both coincide
i = 0
for sec in secnamelist:
for childsec in childseclist:
if sec == childsec:
idxvec[i] += 1
i += 1
[idx] = np.where(idxvec)
return idx
[docs] def get_idx_parent_children(self, parent="soma[0]"):
"""
Get all idx of segments of parent and children sections, i.e. segment
idx of sections connected to parent-argument, and also of the parent
segments
Parameters
----------
parent : str
name-pattern matching a sectionname. Defaults to "soma[0]"
"""
seclist = [parent]
sref = neuron.h.SectionRef(parent)
for sec in sref.child:
seclist.append(sec.name())
return self.get_idx(section=seclist)
[docs] def get_idx_name(self, idx=np.array([0], dtype=int)):
'''
Return NEURON convention name of segments with index idx.
The returned argument is a list of tuples with corresponding
segment idx, section name, and position along the section, like;
[(0, 'neuron.h.soma[0]', 0.5),]
kwargs:
::
idx : ndarray, dtype int
segment indices, must be between 0 and cell.totnsegs
'''
#ensure idx is array-like, or convert
if type(idx) == int or np.int64:
idx = np.array([idx])
elif len(idx) == 0:
return
else:
idx = np.array(idx).astype(int)
#ensure all idx are valid
if np.any(idx >= self.totnsegs):
wrongidx = idx[np.where(idx >= self.totnsegs)]
raise Exception('idx %s >= number of compartments' % str(wrongidx))
#create list of seg names:
allsegnames = []
segidx = 0
for sec in self.allseclist:
for seg in sec:
allsegnames.append((segidx, '%s' % sec.name(), seg.x))
segidx += 1
return np.array(allsegnames, dtype=object)[idx][0]
def _collect_pt3d(self):
"""collect the pt3d info, for each section"""
x = []
y = []
z = []
d = []
for sec in self.allseclist:
n3d = int(neuron.h.n3d())
x_i, y_i, z_i = np.zeros(n3d), np.zeros(n3d), np.zeros(n3d),
d_i = np.zeros(n3d)
for i in range(n3d):
x_i[i] = neuron.h.x3d(i)
y_i[i] = neuron.h.y3d(i)
z_i[i] = neuron.h.z3d(i)
d_i[i] = neuron.h.diam3d(i)
x.append(x_i)
y.append(y_i)
z.append(z_i)
d.append(d_i)
#remove offsets which may be present if soma is centred in Origo
if len(x) > 1:
xoff = x[0].mean()
yoff = y[0].mean()
zoff = z[0].mean()
for i in range(len(x)):
x[i] -= xoff
y[i] -= yoff
z[i] -= zoff
return x, y, z, d
def _update_pt3d(self):
"""
update the locations in neuron.hoc.space using neuron.h.pt3dchange()
"""
for i, sec in enumerate(self.allseclist):
n3d = int(neuron.h.n3d())
for n in range(n3d):
neuron.h.pt3dchange(n,
self.x3d[i][n],
self.y3d[i][n],
self.z3d[i][n],
self.diam3d[i][n])
#let NEURON know about the changes we just did:
neuron.h.define_shape()
#must recollect the geometry, otherwise we get roundoff errors!
self._collect_geometry()
def _set_pt3d_pos(self, diffx=0, diffy=0, diffz=0):
"""
Offset pt3d geometry with differential displacement
indicated in Cell.set_pos()
"""
for i in range(len(self.x3d)):
self.x3d[i] += diffx
self.y3d[i] += diffy
self.z3d[i] += diffz
self._update_pt3d()
def _set_pt3d_rotation(self, x=None, y=None, z=None, rotation_order='xyz'):
"""
Rotate pt3d geometry of cell object around the x-, y-, z-axis
in the order described by rotation_order parameter.
rotation_order should be a string with 3 elements containing x, y, and z
e.g. 'xyz', 'zyx'
Input should be angles in radians.
using rotation matrices, takes dict with rot. angles,
where x, y, z are the rotation angles around respective axes.
All rotation angles are optional.
Examples
--------
>>> cell = LFPy.Cell(**kwargs)
>>> rotation = {'x' : 1.233, 'y' : 0.236, 'z' : np.pi}
>>> cell.set_pt3d_rotation(**rotation)
"""
for ax in rotation_order:
if ax == 'x' and x is not None:
theta = -x
rotation_x = np.array([[1, 0, 0],
[0, np.cos(theta), -np.sin(theta)],
[0, np.sin(theta), np.cos(theta)]])
for i in range(len(self.x3d)):
rel_pos = self._rel_pt3d_positions(self.x3d[i],
self.y3d[i], self.z3d[i])
rel_pos = np.dot(rel_pos, rotation_x)
self.x3d[i], self.y3d[i], self.z3d[i] = \
self._real_pt3d_positions(rel_pos)
if self.verbose:
print(('Rotated geometry %g radians around x-axis' % (-theta)))
else:
if self.verbose:
print('Geometry not rotated around x-axis')
if ax == 'y' and y is not None:
phi = -y
rotation_y = np.array([[np.cos(phi), 0, np.sin(phi)],
[0, 1, 0],
[-np.sin(phi), 0, np.cos(phi)]])
for i in range(len(self.x3d)):
rel_pos = self._rel_pt3d_positions(self.x3d[i],
self.y3d[i], self.z3d[i])
rel_pos = np.dot(rel_pos, rotation_y)
self.x3d[i], self.y3d[i], self.z3d[i] = \
self._real_pt3d_positions(rel_pos)
if self.verbose:
print('Rotated geometry %g radians around y-axis' % (-phi))
else:
if self.verbose:
print('Geometry not rotated around y-axis')
if ax == 'z' and z is not None:
gamma = -z
rotation_z = np.array([[np.cos(gamma), -np.sin(gamma), 0],
[np.sin(gamma), np.cos(gamma), 0],
[0, 0, 1]])
for i in range(len(self.x3d)):
rel_pos = self._rel_pt3d_positions(self.x3d[i],
self.y3d[i], self.z3d[i])
rel_pos = np.dot(rel_pos, rotation_z)
self.x3d[i], self.y3d[i], self.z3d[i] = \
self._real_pt3d_positions(rel_pos)
if self.verbose:
print('Rotated geometry %g radians around z-axis' % (-gamma))
else:
if self.verbose:
print('Geometry not rotated around z-axis')
self._update_pt3d()
def _rel_pt3d_positions(self, x, y, z):
"""Morphology relative to soma position """
rel_pos = np.transpose(np.array([x - self.somapos[0],
y - self.somapos[1],
z - self.somapos[2]]))
return rel_pos
def _real_pt3d_positions(self, rel_pos):
"""Morphology coordinates relative to Origo """
x = rel_pos[:, 0] + self.somapos[0]
y = rel_pos[:, 1] + self.somapos[1]
z = rel_pos[:, 2] + self.somapos[2]
x = np.array(x).flatten()
y = np.array(y).flatten()
z = np.array(z).flatten()
return x, y, z
def _create_polygon(self, i, projection=('x', 'z')):
"""create a polygon to fill for each section"""
x = getattr(self, projection[0]+'3d')[i]
y = getattr(self, projection[1]+'3d')[i]
#x = self.x3d[i]
#z = self.z3d[i]
d = self.diam3d[i]
#calculate angles
dx = np.diff(x)
dy = np.diff(y)
theta = np.arctan2(dy, dx)
x = np.r_[x, x[::-1]]
y = np.r_[y, y[::-1]]
theta = np.r_[theta, theta[::-1]]
d = np.r_[d, d[::-1]]
#1st corner:
x[0] -= 0.5 * d[0] * np.sin(theta[0])
y[0] += 0.5 * d[0] * np.cos(theta[0])
##pt3d points between start and end of section, first side
x[1:dx.size] -= 0.25 * d[1:dx.size] * (
np.sin(theta[:dx.size-1]) + np.sin(theta[1:dx.size]))
y[1:dy.size] += 0.25 * d[1:dy.size] * (
np.cos(theta[:dy.size-1]) + np.cos(theta[1:dx.size]))
#end of section, first side
x[dx.size] -= 0.5 * d[dx.size] * np.sin(theta[dx.size])
y[dy.size] += 0.5 * d[dy.size] * np.cos(theta[dy.size])
#other side
#end of section, second side
x[dx.size+1] += 0.5 * d[dx.size+1] * np.sin(theta[dx.size])
y[dy.size+1] -= 0.5 * d[dy.size+1] * np.cos(theta[dy.size])
##pt3d points between start and end of section, second side
x[::-1][1:dx.size] += 0.25 * d[::-1][1:dx.size] * (
np.sin(theta[::-1][:dx.size-1]) + np.sin(theta[::-1][1:dx.size]))
y[::-1][1:dy.size] -= 0.25 * d[::-1][1:dy.size] * (
np.cos(theta[::-1][:dy.size-1]) + np.cos(theta[::-1][1:dx.size]))
#last corner:
x[-1] += 0.5 * d[-1] * np.sin(theta[-1])
y[-1] -= 0.5 * d[-1] * np.cos(theta[-1])
return x, y
[docs] def get_pt3d_polygons(self, projection=('x', 'z')):
"""For each section create a polygon in the plane determined by keyword
argument projection=('x', 'z'), that can be
visualized using e.g., plt.fill()
Returns
-------
list
list of (x, z) tuples giving the trajectory
of each section that can be plotted using PolyCollection
Examples
--------
>>> from matplotlib.collections import PolyCollection
>>> import matplotlib.pyplot as plt
>>> cell = LFPy.Cell(morphology='PATH/TO/MORPHOLOGY')
>>> zips = []
>>> for x, z in cell.get_pt3d_polygons(projection=('x', 'z')):
>>> zips.append(zip(x, z))
>>> polycol = PolyCollection(zips,
>>> edgecolors='none',
>>> facecolors='gray')
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.add_collection(polycol)
>>> ax.axis(ax.axis('equal'))
>>> plt.show()
"""
if len(projection) != 2:
raise ValueError("projection arg be a tuple like ('x', 'y')")
if 'x' in projection and 'y' in projection:
pass
elif 'x' in projection and 'z' in projection:
pass
elif 'y' in projection and 'z' in projection:
pass
else:
mssg = "projection must be a length 2 tuple of 'x', 'y' or 'z'!"
raise ValueError(mssg)
try:
assert(self.pt3d is True)
except AssertionError:
raise AssertionError('Cell keyword argument pt3d != True')
polygons = []
for i in range(len(self.x3d)):
polygons.append(self._create_polygon(i, projection))
return polygons
def _create_segment_polygon(self, i, projection=('x', 'z')):
"""Create a polygon to fill for segment i, in the plane
determined by kwarg projection"""
x = [getattr(self, projection[0]+'start')[i],
getattr(self, projection[0]+'end')[i]]
z = [getattr(self, projection[1]+'start')[i],
getattr(self, projection[1]+'end')[i]]
#x = [self.xstart[i], self.xend[i]]
#z = [self.zstart[i], self.zend[i]]
d = self.diam[i]
#calculate angles
dx = np.diff(x)
dz = np.diff(z)
theta = np.arctan2(dz, dx)
x = np.r_[x, x[::-1]]
z = np.r_[z, z[::-1]]
#1st corner:
x[0] -= 0.5 * d * np.sin(theta)
z[0] += 0.5 * d * np.cos(theta)
#end of section, first side
x[1] -= 0.5 * d * np.sin(theta)
z[1] += 0.5 * d * np.cos(theta)
#other side
#end of section, second side
x[2] += 0.5 * d * np.sin(theta)
z[2] -= 0.5 * d * np.cos(theta)
#last corner:
x[3] += 0.5 * d * np.sin(theta)
z[3] -= 0.5 * d * np.cos(theta)
return x, z
[docs] def get_idx_polygons(self, projection=('x', 'z')):
"""For each segment idx in cell create a polygon in the plane
determined by the projection kwarg (default ('x', 'z')),
that can be visualized using plt.fill() or
mpl.collections.PolyCollection
Parameters
----------
projection : tuple of strings
Determining projection. Defaults to ('x', 'z')
Returns
-------
polygons : list
list of (ndarray, ndarray) tuples
giving the trajectory of each section
Examples
--------
The most efficient way of using this would be something like
>>> from matplotlib.collections import PolyCollection
>>> import matplotlib.pyplot as plt
>>> cell = LFPy.Cell(morphology='PATH/TO/MORPHOLOGY')
>>> zips = []
>>> for x, z in cell.get_idx_polygons(projection=('x', 'z')):
>>> zips.append(zip(x, z))
>>> polycol = PolyCollection(zips,
>>> edgecolors='none',
>>> facecolors='gray')
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.add_collection(polycol)
>>> ax.axis(ax.axis('equal'))
>>> plt.show()
"""
if len(projection) != 2:
raise ValueError("projection arg be a tuple like ('x', 'y')")
if 'x' in projection and 'y' in projection:
pass
elif 'x' in projection and 'z' in projection:
pass
elif 'y' in projection and 'z' in projection:
pass
else:
mssg = "projection must be a length 2 tuple of 'x', 'y' or 'z'!"
raise ValueError(mssg)
polygons = []
for i in np.arange(self.totnsegs):
polygons.append(self._create_segment_polygon(i, projection))
return polygons
[docs] def insert_v_ext(self, v_ext, t_ext):
"""Set external extracellular potential around cell.
Playback of some extracellular potential v_ext on each cell.totnseg
compartments. Assumes that the "extracellular"-mechanism is inserted
on each compartment.
Can be used to study ephaptic effects and similar
The inputs will be copied and attached to the cell object as
cell.v_ext, cell.t_ext, and converted
to (list of) neuron.h.Vector types, to allow playback into each
compartment e_extracellular reference.
Can not be deleted prior to running cell.simulate()
Parameters
----------
v_ext : ndarray
Numpy array of size cell.totnsegs x t_ext.size, unit mV
t_ext : ndarray
Time vector of v_ext
Examples
--------
>>> import LFPy
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> #create cell
>>> cell = LFPy.Cell(morphology='morphologies/example_morphology.hoc',
>>> passive=True)
>>> #time vector and extracellular field for every segment:
>>> t_ext = np.arange(cell.tstop / cell.dt+ 1) * cell.dt
>>> v_ext = np.random.rand(cell.totnsegs, t_ext.size)-0.5
>>> #insert potentials and record response:
>>> cell.insert_v_ext(v_ext, t_ext)
>>> cell.simulate(rec_imem=True, rec_vmem=True)
>>> fig = plt.figure()
>>> ax1 = fig.add_subplot(311)
>>> ax2 = fig.add_subplot(312)
>>> ax3 = fig.add_subplot(313)
>>> eim = ax1.matshow(np.array(cell.v_ext), cmap='spectral')
>>> cb1 = fig.colorbar(eim, ax=ax1)
>>> cb1.set_label('v_ext')
>>> ax1.axis(ax1.axis('tight'))
>>> iim = ax2.matshow(cell.imem, cmap='spectral')
>>> cb2 = fig.colorbar(iim, ax=ax2)
>>> cb2.set_label('imem')
>>> ax2.axis(ax2.axis('tight'))
>>> vim = ax3.matshow(cell.vmem, cmap='spectral')
>>> ax3.axis(ax3.axis('tight'))
>>> cb3 = fig.colorbar(vim, ax=ax3)
>>> cb3.set_label('vmem')
>>> ax3.set_xlabel('tstep')
>>> plt.show()
"""
#test dimensions of input
try:
if v_ext.shape[0] != self.totnsegs:
raise ValueError("v_ext.shape[0] != cell.totnsegs")
if v_ext.shape[1] != t_ext.size:
raise ValueError('v_ext.shape[1] != t_ext.size')
except:
raise ValueError('v_ext, t_ext must both be np.array types')
if not self.extracellular:
raise Exception('LFPy.Cell arg extracellular != True')
#create list of extracellular potentials on each segment, time vector
self.t_ext = neuron.h.Vector(t_ext)
self.v_ext = []
for v in v_ext:
self.v_ext.append(neuron.h.Vector(v))
#play v_ext into e_extracellular reference
i = 0
for sec in self.allseclist:
for seg in sec:
self.v_ext[i].play(seg._ref_e_extracellular, self.t_ext)
i += 1
return
[docs] def get_axial_currents_from_vmem(self, timepoints=None):
"""Compute axial currents from cell sim: get current magnitude,
distance vectors and position vectors.
Parameters
----------
timepoints : ndarray, dtype=int
array of timepoints in simulation at which you want to compute
the axial currents. Defaults to False. If not given,
all simulation timesteps will be included.
Returns
-------
i_axial : ndarray, dtype=float
Shape ((cell.totnsegs-1)*2, len(timepoints)) array of axial current
magnitudes I in units of (nA) in cell at all timesteps in timepoints,
or at all timesteps of the simulation if timepoints=None.
Contains two current magnitudes per segment,
(except for the root segment): 1) the current from the mid point of
the segment to the segment start point, and 2) the current from
the segment start point to the mid point of the parent segment.
d_vectors : ndarray, dtype=float
Shape ((cell.totnsegs-1)*2, 3) array of distance vectors traveled by
each axial current in i_axial in units of (µm). The indices of the
first axis, correspond to the first axis of i_axial and pos_vectors.
pos_vectors : ndarray, dtype=float
Shape ((cell.totnsegs-1)*2, 3) array of position vectors pointing to
the mid point of each axial current in i_axial in units of (µm). The
indices of the first axis, correspond to the first axis
of i_axial and d_vectors.
Raises
------
AttributeError
Raises an exeption if the cell.vmem attribute cannot be found
"""
if not hasattr(self, 'vmem'):
raise AttributeError('no vmem, run cell.simulate(rec_vmem=True)')
i_axial = []
d_vectors = []
pos_vectors = []
dseg = np.c_[self.xmid - self.xstart,
self.ymid - self.ystart,
self.zmid - self.zstart]
dpar = np.c_[self.xend - self.xmid,
self.yend - self.ymid,
self.zend - self.zmid]
children_dict = self.get_dict_of_children_idx()
for sec in self.allseclist:
if not neuron.h.SectionRef(sec.name()).has_parent():
if sec.nseg == 1:
# skip soma, since soma is an orphan
continue
else:
# the first segment has more than one segment,
# need to compute axial currents within this section.
seg_idx = 1
parent_idx = 0
bottom_seg = False
first_sec = True
branch = False
parentsec = None
children_dict = None
connection_dict = None
conn_point = 1
else:
# section has parent section
first_sec = False
bottom_seg = True
secref = neuron.h.SectionRef(sec.name())
parentseg = secref.parent()
parentsec = parentseg.sec
children_dict = self.get_dict_of_children_idx()
branch = len(children_dict[parentsec.name()]) > 1
connection_dict = self.get_dict_parent_connections()
conn_point = connection_dict[sec.name()]
# find parent index
if conn_point == 1 or parentsec.nseg == 1:
internal_parent_idx = -1 # last seg in sec
elif conn_point == 0:
internal_parent_idx = 0 # first seg in sec
else:
# if parentseg is not first or last seg in parentsec
segment_xlist = np.array([segment.x for segment in parentsec])
internal_parent_idx = np.abs(segment_xlist - conn_point).argmin()
parent_idx = self.get_idx(section=parentsec.name())[internal_parent_idx]
# find segment index
seg_idx = self.get_idx(section=sec.name())[0]
for _ in sec:
if first_sec:
first_sec = False
continue
iseg, ipar = self._parent_and_segment_current(seg_idx,
parent_idx,
bottom_seg,
branch,
parentsec,
children_dict,
connection_dict,
conn_point,
timepoints
)
if bottom_seg:
# if a seg is connected to soma, it is
# connected to the middle of soma,
# and dpar needs to be altered.
par_dist = np.array([(self.xstart[seg_idx] -
self.xmid[parent_idx]),
(self.ystart[seg_idx] -
self.ymid[parent_idx]),
(self.zstart[seg_idx] -
self.zmid[parent_idx])])
else:
par_dist = dpar[parent_idx]
d_vectors.append(par_dist)
d_vectors.append(dseg[seg_idx])
i_axial.append(ipar)
i_axial.append(iseg)
pos_par = np.array([self.xstart[seg_idx],
self.ystart[seg_idx],
self.zstart[seg_idx]]) - 0.5*par_dist
pos_seg = np.array([self.xmid[seg_idx],
self.ymid[seg_idx],
self.zmid[seg_idx]]) - 0.5*dseg[seg_idx]
pos_vectors.append(pos_par)
pos_vectors.append(pos_seg)
parent_idx = seg_idx
seg_idx += 1
branch = False
bottom_seg = False
parent_ri = 0
return np.array(i_axial), np.array(d_vectors), np.array(pos_vectors)
[docs] def get_axial_resistance(self):
"""
Return NEURON axial resistance for all cell compartments.
Returns
-------
ri_list : ndarray, dtype=float
Shape (cell.totnsegs, ) array containing neuron.h.ri(seg.x) in units
of (MOhm) for all segments in cell calculated using the
neuron.h.ri(seg.x) method. neuron.h.ri(seg.x) returns the
axial resistance from the middle of the segment to the middle of
the parent segment. Note: If seg is the first segment in a section,
i.e. the parent segment belongs to a different section or there is
no parent section, then neuron.h.ri(seg.x) returns the axial
resistance from the middle of the segment to the node connecting the
segment to the parent section (or a ghost node if there is no
parent)
"""
ri_list = np.zeros(self.totnsegs)
comp = 0
for sec in self.allseclist:
for seg in sec:
ri_list[comp] = neuron.h.ri(seg.x)
comp += 1
return ri_list
[docs] def get_dict_of_children_idx(self):
"""
Return dictionary with children segment indices for all sections.
Returns
-------
children_dict : dictionary
Dictionary containing a list for each section,
with the segment index of all the section's children.
The dictionary is needed to find the
sibling of a segment.
"""
children_dict = {}
for sec in self.allseclist:
children_dict[sec.name()] = []
for child in neuron.h.SectionRef(sec.name()).child:
# add index of first segment of each child
children_dict[sec.name()].append(int(self.get_idx(
section=child.name())[0]))
return children_dict
[docs] def get_dict_parent_connections(self):
"""
Return dictionary with parent connection point for all sections.
Returns
-------
connection_dict : dictionary
Dictionary containing a float in range [0, 1] for each section
in cell. The float gives the location on the parent segment
to which the section is connected.
The dictionary is needed for computing axial currents.
"""
connection_dict = {}
for sec in self.allseclist:
connection_dict[sec.name()] = neuron.h.parent_connection()
return connection_dict
def _parent_and_segment_current(self, seg_idx, parent_idx, bottom_seg,
branch=False, parentsec=None,
children_dict=None, connection_dict=None,
conn_point=1, timepoints=None):
"""
Return axial current from segment (seg_idx) mid to segment start,
and current from parent segment (parent_idx) end to parent segment mid.
Parameters
----------
seg_idx : int
Segment index
parent_idx : int
Parent index
parent_ri : float
Axial resistance from parent end to mid in units of (MΩ)
bottom_seg : boolean
branch : boolean
parentsec : neuron.Section object
timepoints : ndarray, dtype=int
array of timepoints in simulation at which you want to compute
the axial currents. Defaults to None. If not given,
all simulation timesteps will be included.
Returns
-------
iseg : dtype=float
Axial current in units of (nA)
from segment mid point to segment start point.
ipar : dtype=float
Axial current in units of (nA)
from parent segment end point to parent segment mid point.
"""
# list of axial resistance between segments
ri_list = self.get_axial_resistance()
# axial resistance between segment mid and parent node
seg_ri = ri_list[seg_idx]
vmem = self.vmem
if timepoints is not None:
vmem = self.vmem[:,timepoints]
vpar = vmem[parent_idx]
vseg = vmem[seg_idx]
# if segment is the first in its section and it is connected to
# top or bottom of parent section, we need to find parent_ri explicitly
if bottom_seg and (conn_point == 0 or conn_point == 1):
if conn_point == 0:
parent_ri = ri_list[parent_idx]
else:
parent_ri = neuron.h.ri(0)
if not branch:
ri = parent_ri + seg_ri
iseg = (vpar - vseg) / ri
ipar = iseg
else:
# if branch, need to compute iseg and ipar separately
[sib_idcs] = np.take(children_dict[parentsec.name()],
np.where(children_dict[parentsec.name()]
!= seg_idx))
sibs = [self.get_idx_name(sib_idcs)[i][1] for i in range(len(sib_idcs))]
# compute potential in branch point between parent and siblings
v_branch_num = vpar/parent_ri + vseg/seg_ri
v_branch_denom = 1./parent_ri + 1./seg_ri
for sib_idx, sib in zip(sib_idcs, sibs):
sib_conn_point = connection_dict[sib]
if sib_conn_point == conn_point:
v_branch_num += vmem[sib_idx]/ri_list[sib_idx]
v_branch_denom += 1./ ri_list[sib_idx]
v_branch = v_branch_num/v_branch_denom
iseg = (v_branch - vseg)/seg_ri
# set ipar=iseg
# only fraction of total current into parent is added per sibling
ipar = iseg
else:
iseg = (vpar - vseg) / seg_ri
ipar = iseg
return iseg, ipar
[docs] def distort_geometry(self, factor=0., axis='z', nu=0.0):
"""
Distorts cellular morphology with a relative factor along a chosen axis
preserving Poisson's ratio. A ratio nu=0.5 assumes uncompressible and
isotropic media that embeds the cell. A ratio nu=0 will only affect
geometry along the chosen axis. A ratio nu=-1 will isometrically scale
the neuron geometry along each axis.
This method does not affect the underlying cable properties of the cell,
only predictions of extracellular measurements (by affecting the
relative locations of sources representing the compartments).
Parameters
----------
factor : float
relative compression/stretching factor of morphology. Default is 0
(no compression/stretching). Positive values implies a compression
along the chosen axis.
axis : str
which axis to apply compression/stretching. Default is "z".
nu : float
Poisson's ratio. Ratio between axial and transversal
compression/stretching. Default is 0.
"""
try:
assert(abs(factor) < 1.)
except AssertionError:
raise AssertionError('abs(factor) >= 1, factor must be in <-1, 1>')
try:
assert(axis in ['x', 'y', 'z'])
except AssertionError:
raise AssertionError('axis={} not "x", "y" or "z"'.format(axis))
for pos, dir_ in zip(self.somapos, 'xyz'):
geometry = np.c_[getattr(self, dir_+'start'),
getattr(self, dir_+'mid'),
getattr(self, dir_+'end')]
if dir_ == axis:
geometry -= pos
geometry *= (1. - factor)
geometry += pos
else:
geometry -= pos
geometry *= (1. + factor*nu)
geometry += pos
setattr(self, dir_+'start', geometry[:, 0])
setattr(self, dir_+'mid', geometry[:, 1])
setattr(self, dir_+'end', geometry[:, 2])
# recompute length of each segment
self.length = np.sqrt((self.xend - self.xstart)**2 +
(self.yend - self.ystart)**2 +
(self.zend - self.zstart)**2)
[docs] def get_multi_current_dipole_moments(self, timepoints=None):
'''
Return 3D current dipole moment vector and middle position vector
from each axial current in space.
Parameters
----------
timepoints : ndarray, dtype=int
array of timepoints at which you want to compute
the current dipole moments. Defaults to None. If not given,
all simulation timesteps will be included.
Returns
-------
multi_dipoles : ndarray, dtype = float
Shape (n_axial_currents, n_timepoints, 3) array
containing the x-,y-,z-components of the current dipole moment
from each axial current in cell, at all timepoints.
The number of axial currents, n_axial_currents = (cell.totnsegs-1)*2
and the number of timepoints, n_timepoints = cell.tvec.size.
The current dipole moments are given in units of (nA µm).
pos_axial : ndarray, dtype = float
Shape (n_axial_currents, 3) array containing the x-, y-, and
z-components giving the mid position in space of each multi_dipole
in units of (µm).
Examples
--------
Get all current dipole moments and positions from all axial currents in a
single neuron simulation.
>>> import LFPy
>>> import numpy as np
>>> cell = LFPy.Cell('PATH/TO/MORPHOLOGY', extracellular=False)
>>> syn = LFPy.Synapse(cell, idx=cell.get_closest_idx(0,0,1000),
>>> syntype='ExpSyn', e=0., tau=1., weight=0.001)
>>> syn.set_spike_times(np.mgrid[20:100:20])
>>> cell.simulate(rec_vmem=True, rec_imem=False)
>>> timepoints = np.array([1,2,3,4])
>>> multi_dipoles, dipole_locs = cell.get_multi_current_dipole_moments(timepoints=timepoints)
'''
i_axial, d_axial, pos_axial = self.get_axial_currents_from_vmem(timepoints=timepoints)
Ni, Nt = i_axial.shape
multi_dipoles = np.array([i_axial[i][:, np.newaxis]*d_axial[i] for i in range(Ni)])
return multi_dipoles, pos_axial