Source code for LFPy.cell

#!/usr/bin/env python
'''
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.
'''

import os
import neuron
import numpy as np
import pickle
from LFPy import RecExtElectrode
from LFPy.run_simulation import _run_simulation, _run_simulation_with_electrode
from LFPy.run_simulation import _collect_geometry_neuron
from LFPy.alias_method import alias_method
import sys
from warnings import warn

[docs]class Cell(object): ''' The main cell class used in LFPy. Arguments: :: morphology : [str]: path/to/morphology/file v_init: [-65.]: initial potential passive: [True]/False: passive mechs are initialized if True Ra: [150.]: axial resistance rm: [30000]: membrane resistivity cm: [1.0]: membrane capacitance e_pas: [-65.]: passive mechanism reversal potential extracellular: [True]/False: switch for NEURON's extracellular mechanism timeres_NEURON: [0.1]: internal dt for NEURON simulation timeres_python: [0.1]: overall dt for python simulation tstartms: [0.]: initialization time for simulation <= 0 ms tstopms: [100.]: stop time for simulation > 0 ms nsegs_method: ['lambda100']/'lambda_f'/'fixed_length': nseg rule max_nsegs_length: [None]: max segment length for method 'fixed_length' lambda_f: [100]: AC frequency for method 'lambda_f' d_lambda: [0.1]: parameter for d_lambda rule delete_sections: [True]: delete pre-existing section-references custom_code: [None]: list of model-specific code files ([.py/.hoc]) custom_fun: [None]: list of model-specific functions with args custom_fun_args: [None]: list of args passed to custom_fun functions pt3d: True/[False]: use pt3d-info of the cell geometries switch celsius: [None]: Temperature in celsius. If nothing is specified here or in custom code it is 6.3 C verbose: True/[False]: verbose output switch Usage of cell class: :: import LFPy cellParameters = { 'morphology' : 'path/to/morphology', 'rm' : 30000, 'cm' : 1.0, 'Ra' : 150, 'timeres_NEURON' : 0.1, 'timeres_python' : 0.1, 'tstartms' : -50, 'tstopms' : 50, } cell = LFPy.Cell(**cellParameters) cell.simulate() ''' def __init__(self, morphology, v_init=-65., passive = True, Ra=150, rm=30000, cm=1.0, e_pas=-65., extracellular = True, timeres_NEURON=2**-3, timeres_python=2**-3, tstartms=0, tstopms=100, 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): ''' Initialization of the Cell object. ''' self.verbose = verbose self.pt3d = pt3d 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 self.morphology = morphology if self.morphology is not None: if os.path.isfile(self.morphology): self._load_geometry() else: raise Exception('non-existent file %s' % self.morphology) else: try: #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: raise Exception("Could not load existent top-level cell") #Some parameters and lists initialised if timeres_python not in 2.**np.arange(-16, -1) or timeres_NEURON \ not in 2.**np.arange(-16, -1): if self.verbose: print('timeres_python and or timeres_NEURON not a power of 2,') print('cell.tvec errors may occur.') print('Initialization will continue.') else: pass if timeres_python < timeres_NEURON: raise ValueError('timeres_python = %.3e < timeres_NEURON = %.3e' % (timeres_python, timeres_NEURON)) self.timeres_python = timeres_python self.timeres_NEURON = timeres_NEURON self.tstartms = tstartms self.tstopms = tstopms 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.rm = rm self.cm = cm self.e_pas = e_pas self._set_passive() else: 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, can't access imem!") #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 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 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. ',]) raise Exception(ERRMSG) elif code.split('.')[-1] == 'py': 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 self.morphology is not None: 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 _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 kwargs: :: frequency: float, frequency at whihc 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 = 1. / self.rm sec.e_pas = self.e_pas def _set_extracellular(self): '''Insert extracellular mechanism on all sections to access i_membrane''' 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 syntype (e.g. ExpSyn) synapse on segment with index idx, Arguments: :: idx : int syntype : str record_current : bool record_potential : bool weight : float 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 cmd1 = 'syn = neuron.h.' cmd2 = '(seg.x, sec=sec)' for sec in self.allseclist: for seg in sec: if i == idx: command = cmd1 + syntype + cmd2 exec(command) for param in list(kwargs.keys()): try: 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.tstopms / self.timeres_python+1)) synirec.record(syn._ref_i, self.timeres_python) self.synireclist.append(synirec) #record potential if record_potential: synvrec = neuron.h.Vector(int(self.tstopms / self.timeres_python+1)) synvrec.record(seg._ref_v, self.timeres_python) self.synvreclist.append(synvrec) i += 1 return self.synlist.count() - 1
[docs] def set_point_process(self, idx, pptype, record_current=False, **kwargs): ''' Insert pptype-electrode type pointprocess on segment numbered idx on cell object, with keyword arguments according to types SEClamp, VClamp, IClamp, SinIClamp, ChirpIClamp. Arguments: :: idx : int pptype : str record_current : bool kwargs : arguments 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() 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 exec(command) for param in list(kwargs.keys()): try: 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.tstopms / self.timeres_python+1)) stimirec.record(stim._ref_i, self.timeres_python) self.stimireclist.append(stimirec) 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=-10000, z_max=10000): ''' Returns neuron idx of segments from sections with names that match the pattern defined in input section on interval [z_min, z_max]. kwargs: :: section: str, any entry in cell.allsecnames or just 'allsec'. z_min: float, depth filter z_max: float depth filter Usage: :: 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 kwargs: :: x: float, coordinate y: float, coordinate z: float, coordinate section: str, string matching a section-name ''' idx = self.get_idx(section) dist = np.sqrt((self.xmid[idx] - x)**2 + (self.ymid[idx] - y)**2 + (self.zmid[idx] - z)**2) mindist = np.where(dist == np.min(dist)) return int(idx[mindist])
[docs] def get_rand_idx_area_norm(self, section='allsec', nidx=1, z_min=-10000, z_max=10000): '''Return nidx segment indices in section with random probability normalized to the membrane area of segment on interval [z_min, z_max] kwargs: :: 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() idx = alias_method(poss_idx, area, nidx) return idx
[docs] def simulate(self, electrode=None, rec_imem=False, rec_vmem=False, rec_ipas=False, rec_icap=False, rec_isyn=False, rec_vmemsyn=False, rec_istim=False, rec_variables=[], variable_dt=False, atol=0.001, to_memory=True, to_file=False, file_name=None, dotprodcoeffs=None): ''' This is the main function running the simulation of the NEURON model. Start NEURON simulation and record variables specified by arguments. Arguments: :: electrode: 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: 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: record segment membrane voltages (mV) rec_ipas: record passive segment membrane currents (nA) rec_icap: record capacitive segment membrane currents (nA) rec_isyn: record synaptic currents of from Synapse class (nA) rec_vmemsyn: record membrane voltage of segments with Synapse(mV) rec_istim: record currents of StimIntraElectrode (nA) rec_variables: list of variables to record, i.e arg=['cai', ] variable_dt: boolean, using variable timestep in NEURON atol: absolute tolerance used with NEURON variable timestep to_memory: only valid with electrode, store lfp in -> electrode.LFP to_file: only valid with electrode, save LFPs in hdf5 file format file_name: name of hdf5 file, '.h5' is appended if it doesnt exist dotprodcoeffs : list of N x Nseg np.ndarray. These arrays will at every timestep be multiplied by the membrane currents. Presumably useful for memory efficient csd or lfp calcs ''' self._set_soma_volt_recorder() self._collect_tvec() 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 len(rec_variables) > 0: self._set_variable_recorders(rec_variables) #run fadvance until t >= tstopms, and calculate LFP if asked for if electrode is None and dotprodcoeffs is None: if not rec_imem: print(("rec_imem = %s, membrane currents will not be recorded!" \ % str(rec_imem))) _run_simulation(self, variable_dt, atol) else: #allow using both electrode and additional coefficients: _run_simulation_with_electrode(self, electrode, variable_dt, atol, to_memory, to_file, file_name, dotprodcoeffs) #somatic trace 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() if rec_isyn: self._collect_isyn() if rec_vmemsyn: self._collect_vsyn() if rec_istim: self._collect_istim() 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.tstopms / self.timeres_python + 1) \ * self.timeres_python 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) for i in range(self.imem.shape[0]): self.imem[i, ] *= self.area[i] * 1E-2 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) else: raise Exception('must set record_current=True in Synapse class') self.synireclist = None del self.synireclist def _collect_vsyn(self): ''' Collect the membrane voltage of segments with synapses ''' for i in range(len(self.synapses)): self.synapses[i].collect_potential(self) self.synvreclist = None del self.synvreclist def _collect_istim(self): ''' Get the pointprocess currents ''' for i in range(len(self.pointprocesses)): if self.pointprocesses[i].record_current: self.pointprocesses[i].collect_current(self) else: raise Exception('must set record_current=True for pointp.') self.stimireclist = None del self.stimireclist 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])) # elif len(self.synlist) > 0 and len(self.sptimeslist) == 0: # errmsg = 'please run method "set_spike_times() for every' + \ # '\n' + 'instance of LFPy.pointprocess.Synapse' # raise Exception(errmsg) # else: # pass def _set_soma_volt_recorder(self): ''' Record somatic membrane potential ''' self.somav = neuron.h.Vector(int(self.tstopms / self.timeres_python+1)) if self.nsomasec == 0: pass elif self.nsomasec == 1: for sec in self.somalist: self.somav.record(sec(0.5)._ref_v, self.timeres_python) elif self.nsomasec > 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.timeres_python) else: if k == i*2: #record from one of the middle segments: self.somav.record(seg._ref_v, self.timeres_python) 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.tstopms / self.timeres_python+1)) memirec.record(seg._ref_i_membrane, self.timeres_python) 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.tstopms / self.timeres_python+1)) memipasrec.record(seg._ref_i_pas, self.timeres_python) 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.tstopms / self.timeres_python+1)) memicaprec.record(seg._ref_i_cap, self.timeres_python) 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.tstopms / self.timeres_python+1)) memvrec.record(seg._ref_v, self.timeres_python) self.memvreclist.append(memvrec) 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.tstopms / self.timeres_python + 1)) if hasattr(seg, variable): recvector.record(getattr(seg, '_ref_%s' % variable), self.timeres_python) else: print('non-existing variable %s, section %s.%f' % (variable, sec.name(), seg.x)) variablereclist.append(recvector)
[docs] def set_pos(self, xpos=0, ypos=0, zpos=0): ''' Move the cell geometry so that midpoint of soma section is in (xpos, ypos, zpos). If no soma pos, use the first segment ''' diffx = xpos-self.somapos[0] diffy = ypos-self.somapos[1] diffz = zpos-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] = xpos self.somapos[1] = ypos self.somapos[2] = zpos 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): ''' Save data in cell to filename, using cPickle. It will however destroy any neuron.h objects upon saving, as c-objects cannot be pickled Usage: :: 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() filen = open(filename, 'wb') pickle.dump(self, filen, protocol=2) filen.close()
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): ''' Rotate geometry of cell object around the x-, y-, z-axis in that order. 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. Usage: :: cell = LFPy.Cell(**kwargs) rotation = {'x' : 1.233, 'y' : 0.236, 'z' : np.pi} cell.set_rotation(**rotation) ''' if x is not None: theta = -x rotation_x = np.matrix([[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 = rel_start * rotation_x rel_end = 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 y is not None: phi = -y rotation_y = np.matrix([[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 = rel_start * rotation_y rel_end = 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 z is not None: gamma = -z rotation_z = np.matrix([[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 = rel_start * rotation_z rel_end = 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)
[docs] def chiral_morphology(self, axis='x'): ''' Mirror the morphology around given axis, (default x-axis), useful to introduce more heterogeneouties in morphology shapes kwargs: :: 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. Prob. determined by area. kwargs: :: section: str, string matching a section-name 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]), z_min=-10000, z_max=10000): ''' Return the normalized probability (0-1) for synaptic coupling on segments in idx-array. Normalised probability determined by area of segments. kwargs: :: idx : np.ndarray, dtype=int. array of segment indices z_min: float, depth filter z_max: float depth filter ''' 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. kwargs: :: idx0 : int idx1 : int ''' vector = [] try: if idx1 < 0 or idx0 < 0: raise Exception 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 with indices idx0 and idx1. 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 kwargs: :: parent: str name-pattern matching a sectionname ''' 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 kwargs: :: parent: str name-pattern matching a sectionname ''' idxvec = np.zeros(self.totnsegs) secnamelist = [] childseclist = [parent] #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 np.r_[self.get_idx(parent), idx]
[docs] def get_idx_name(self, idx=np.array([0])): ''' 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 : np.ndarray, dtype int segment indices, must be between 0 and cell.totnsegs ''' #ensure idx is array-like, or convert if type(idx) == int: 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 allsegnames[idx]
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() ''' i = 0 for sec in 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]) i += 1 #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): ''' Rotate pt3d geometry of cell object around the x-, y-, z-axis in that order. 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. Usage: :: cell = LFPy.Cell(**kwargs) rotation = {'x' : 1.233, 'y' : 0.236, 'z' : np.pi} cell.set_pt3d_rotation(**rotation) ''' if x is not None: theta = -x rotation_x = np.matrix([[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 = 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 y is not None: phi = -y rotation_y = np.matrix([[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 = 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 z is not None: gamma = -z rotation_z = np.matrix([[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 = 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() Returned argument is a list of (x, z) tuples giving the trajectory of each section that can be plotted using PolyCollection :: 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 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 celll 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 Returned argument is a list of (np.ndarray, np.ndarray) tuples giving the trajectory of each section kwargs: :: projection : ('x', 'z') tuple of two strings determining projection 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): ''' 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() Args: :: v_ext : cell.totnsegs x t_ext.size np.array, unit mV t_ext : np.array, time vector of v_ext Simple usage: :: import LFPy import numpy as np import matplotlib.pyplot as plt #create cell cell = LFPy.Cell(morphology='morphologies/example_morphology.hoc') #time vector and extracellular field for every segment: t_ext = np.arange(cell.tstopms / cell.timeres_python+ 1) * \ cell.timeres_python 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