Source code for LFPy.recextelectrode

#!/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 sys
import warnings
import numpy as np
from . import lfpcalc, tools

[docs]class RecExtElectrode(object): """class RecExtElectrode Main class that represents an extracellular electric recording devices such as a laminar probe. Parameters ---------- cell : None or object If not None, instantiation of LFPy.Cell, LFPy.TemplateCell or similar. sigma : float or list/ndarray of floats extracellular conductivity in units of (S/m). A scalar value implies an isotropic extracellular conductivity. If a length 3 list or array of floats is provided, these values corresponds to an anisotropic conductor with conductivities [sigma_x, sigma_y, sigma_z] accordingly. x, y, z : np.ndarray coordinates or arrays of coordinates in units of (um). Must be same length N : None or list of lists Normal vectors [x, y, z] of each circular electrode contact surface, default None r : float radius of each contact surface, default None n : int if N is not None and r > 0, the number of discrete points used to compute the n-point average potential on each circular contact point. contact_shape : str 'circle'/'square' (default 'circle') defines the contact point shape If 'circle' r is the radius, if 'square' r is the side length method : str switch between the assumption of 'linesource', 'pointsource', 'soma_as_point' to represent each compartment when computing extracellular potentials from_file : bool if True, load cell object from file cellfile : str path to cell pickle verbose : bool Flag for verbose output, i.e., print more information seedvalue : int random seed when finding random position on contact with r > 0 Examples -------- Compute extracellular potentials after simulating and storage of transmembrane currents with the LFPy.Cell class: >>> import numpy as np >>> import matplotlib.pyplot as plt >>> import LFPy >>> >>> cellParameters = { >>> 'morphology' : 'examples/morphologies/L5_Mainen96_LFPy.hoc', # morphology file >>> 'v_init' : -65, # initial voltage >>> 'cm' : 1.0, # membrane capacitance >>> 'Ra' : 150, # axial resistivity >>> 'passive' : True, # insert passive channels >>> 'passive_parameters' : {"g_pas":1./3E4, "e_pas":-65}, # passive params >>> 'dt' : 2**-4, # simulation time res >>> 'tstart' : 0., # start t of simulation >>> 'tstop' : 50., # end t of simulation >>> } >>> cell = LFPy.Cell(**cellParameters) >>> >>> synapseParameters = { >>> 'idx' : cell.get_closest_idx(x=0, y=0, z=800), # compartment >>> 'e' : 0, # reversal potential >>> 'syntype' : 'ExpSyn', # synapse type >>> 'tau' : 2, # syn. time constant >>> 'weight' : 0.01, # syn. weight >>> 'record_current' : True # syn. current record >>> } >>> synapse = LFPy.Synapse(cell, **synapseParameters) >>> synapse.set_spike_times(np.array([10., 15., 20., 25.])) >>> >>> cell.simulate(rec_imem=True) >>> >>> N = np.empty((16, 3)) >>> for i in xrange(N.shape[0]): N[i,] = [1, 0, 0] #normal vec. of contacts >>> electrodeParameters = { #parameters for RecExtElectrode class >>> 'sigma' : 0.3, #Extracellular potential >>> 'x' : np.zeros(16)+25, #Coordinates of electrode contacts >>> 'y' : np.zeros(16), >>> 'z' : np.linspace(-500,1000,16), >>> 'n' : 20, >>> 'r' : 10, >>> 'N' : N, >>> } >>> electrode = LFPy.RecExtElectrode(cell, **electrodeParameters) >>> electrode.calc_lfp() >>> plt.matshow(electrode.LFP) >>> plt.colorbar() >>> plt.axis('tight') >>> plt.show() Compute extracellular potentials during simulation (recommended): >>> import numpy as np >>> import matplotlib.pyplot as plt >>> import LFPy >>> >>> cellParameters = { >>> 'morphology' : 'examples/morphologies/L5_Mainen96_LFPy.hoc', # morphology file >>> 'v_init' : -65, # initial voltage >>> 'cm' : 1.0, # membrane capacitance >>> 'Ra' : 150, # axial resistivity >>> 'passive' : True, # insert passive channels >>> 'passive_parameters' : {"g_pas":1./3E4, "e_pas":-65}, # passive params >>> 'dt' : 2**-4, # simulation time res >>> 'tstart' : 0., # start t of simulation >>> 'tstop' : 50., # end t of simulation >>> } >>> cell = LFPy.Cell(**cellParameters) >>> >>> synapseParameters = { >>> 'idx' : cell.get_closest_idx(x=0, y=0, z=800), # compartment >>> 'e' : 0, # reversal potential >>> 'syntype' : 'ExpSyn', # synapse type >>> 'tau' : 2, # syn. time constant >>> 'weight' : 0.01, # syn. weight >>> 'record_current' : True # syn. current record >>> } >>> synapse = LFPy.Synapse(cell, **synapseParameters) >>> synapse.set_spike_times(np.array([10., 15., 20., 25.])) >>> >>> N = np.empty((16, 3)) >>> for i in xrange(N.shape[0]): N[i,] = [1, 0, 0] #normal vec. of contacts >>> electrodeParameters = { #parameters for RecExtElectrode class >>> 'sigma' : 0.3, #Extracellular potential >>> 'x' : np.zeros(16)+25, #Coordinates of electrode contacts >>> 'y' : np.zeros(16), >>> 'z' : np.linspace(-500,1000,16), >>> 'n' : 20, >>> 'r' : 10, >>> 'N' : N, >>> } >>> electrode = LFPy.RecExtElectrode(**electrodeParameters) >>> >>> cell.simulate(electrode=electrode) >>> >>> plt.matshow(electrode.LFP) >>> plt.colorbar() >>> plt.axis('tight') >>> plt.show() """ def __init__(self, cell=None, sigma=0.3, x=np.array([0]), y=np.array([0]), z=np.array([0]), N=None, r=None, n=None, contact_shape='circle', perCellLFP=False, method='linesource', from_file=False, cellfile=None, verbose=False, seedvalue=None, **kwargs): """Initialize RecExtElectrode class""" self.sigma = sigma if type(sigma) in [list, np.ndarray]: self.sigma = np.array(sigma) if not self.sigma.shape == (3,): raise ValueError("Conductivity, sigma, should be float " "or array of length 3: " "[sigma_x, sigma_y, sigma_z]") self.anisotropic = True else: self.sigma = sigma self.anisotropic = False if type(x) in [float, int]: self.x = np.array([x]) else: self.x = np.array(x).flatten() if type(y) in [float, int]: self.y = np.array([y]) else: self.y = np.array(y).flatten() if type(z) in [float, int]: self.z = np.array([z]) else: self.z = np.array(z).flatten() try: assert((self.x.size==self.y.size) and (self.x.size==self.z.size)) except AssertionError: raise AssertionError("The number of elements in [x, y, z] must be identical") if N is not None: if type(N) != np.array: try: N = np.array(N) except: print('Keyword argument N could not be converted to a ' 'numpy.ndarray of shape (n_contacts, 3)') print(sys.exc_info()[0]) raise if N.shape[-1] == 3: self.N = N else: self.N = N.T if N.shape[-1] != 3: raise Exception('N.shape must be (n_contacts, 1, 3)!') else: self.N = N self.r = r self.n = n if contact_shape is None: self.contact_shape = 'circle' elif contact_shape in ['circle', 'square']: self.contact_shape = contact_shape else: raise ValueError('The contact_shape argument must be either: ' 'None, \'circle\', \'square\'') self.perCellLFP = perCellLFP self.method = method self.verbose = verbose self.seedvalue = seedvalue self.kwargs = kwargs #None-type some attributes created by the Cell class self.electrodecoeff = None self.circle = None self.offsets = None if from_file: if type(cellfile) == type(str()): cell = tools.load(cellfile) elif type(cellfile) == type([]): cell = [] for fil in cellfile: cell.append(tools.load(fil)) else: raise ValueError('cell either string or list of strings') if cell is not None: self.set_cell(cell) if method == 'soma_as_point': if self.anisotropic: self.lfp_method = lfpcalc.calc_lfp_soma_as_point_anisotropic else: self.lfp_method = lfpcalc.calc_lfp_soma_as_point elif method == 'som_as_point': raise RuntimeError('The method "som_as_point" is deprecated.' 'Use "soma_as_point" instead') elif method == 'linesource': if self.anisotropic: self.lfp_method = lfpcalc.calc_lfp_linesource_anisotropic else: self.lfp_method = lfpcalc.calc_lfp_linesource elif method == 'pointsource': if self.anisotropic: self.lfp_method = lfpcalc.calc_lfp_pointsource_anisotropic else: self.lfp_method = lfpcalc.calc_lfp_pointsource else: raise ValueError("LFP method not recognized. " "Should be 'soma_as_point', 'linesource' " "or 'pointsource'")
[docs] def set_cell(self, cell): """Set the supplied cell object as attribute "cell" of the RecExtElectrode object Parameters ---------- cell : obj `LFPy.Cell` or `LFPy.TemplateCell` instance. Returns ------- None """ self.cell = cell if self.cell is not None: self.r_limit = self.cell.diam/2 self.mapping = np.zeros((self.x.size, len(cell.xmid)))
def _test_imem_sum(self, tolerance=1E-8): """Test that the membrane currents sum to zero""" if type(self.cell) == dict or type(self.cell) == list: raise DeprecationWarning('no support for more than one cell-object') sum_imem = self.cell.imem.sum(axis=0) #check if eye matrix is supplied: if ((self.cell.imem.shape == (self.cell.totnsegs, self.cell.totnsegs)) and (np.all(self.cell.imem == np.eye(self.cell.totnsegs)))): pass else: if abs(sum_imem).max() >= tolerance: warnings.warn('Membrane currents do not sum to zero') [inds] = np.where((abs(sum_imem) >= tolerance)) if self.cell.verbose: for i in inds: print('membrane current sum of celltimestep %i: %.3e' % (i, sum_imem[i])) else: pass
[docs] def calc_mapping(self, cell): """Creates a linear mapping of transmembrane currents of each segment of the supplied cell object to contribution to extracellular potential at each electrode contact point of the RexExtElectrode object. Sets the class attribute "mapping", which is a shape (n_contact, n_segs) ndarray, such that the extracellular potential at the contacts phi = np.dot(mapping, I_mem) where I_mem is a shape (n_segs, n_tsteps) ndarray with transmembrane currents for each time step of the simulation. Parameters ---------- cell : obj `LFPy.Cell` or `LFPy.TemplateCell` instance. Returns ------- mapping : ndarray The attribute RecExtElectrode.mapping is returned (optional) """ if cell is not None: self.set_cell(cell) if self.n is not None and self.N is not None and self.r is not None: if self.n <= 1: raise ValueError("n = %i must be larger that 1" % self.n) else: pass self._lfp_el_pos_calc_dist() if self.verbose: print('calculations finished, %s, %s' % (str(self), str(self.cell))) else: self._loop_over_contacts() if self.verbose: print('calculations finished, %s, %s' % (str(self), str(self.cell))) # return mapping return self.mapping
[docs] def calc_lfp(self, t_indices=None, cell=None): """Calculate LFP on electrode geometry from all cell instances. Will chose distributed calculated if electrode contain 'n', 'N', and 'r' Parameters ---------- cell : obj, optional `LFPy.Cell` or `LFPy.TemplateCell` instance. Must be specified here if it was not specified at the initiation of the `RecExtElectrode` class t_indices : np.ndarray Array of timestep indexes where extracellular potential should be calculated. """ self.calc_mapping(cell) if t_indices is not None: currmem = self.cell.imem[:, t_indices] else: currmem = self.cell.imem self._test_imem_sum() self.LFP = np.dot(self.mapping, currmem)
# del self.mapping def _loop_over_contacts(self, **kwargs): """Loop over electrode contacts, and return LFPs across channels""" for i in range(self.x.size): self.mapping[i, :] = self.lfp_method(self.cell, x = self.x[i], y = self.y[i], z = self.z[i], sigma = self.sigma, r_limit = self.r_limit, **kwargs) def _lfp_el_pos_calc_dist(self, **kwargs): """ Calc. of LFP over an n-point integral approximation over flat electrode surface: circle of radius r or square of side r. The locations of these n points on the electrode surface are random, within the given surface. """ # lfp_el_pos = np.zeros(self.LFP.shape) self.offsets = {} self.circle_circ = {} def create_crcl(i): """make circumsize of contact point""" crcl = np.zeros((self.n, 3)) for j in range(self.n): B = [(np.random.rand()-0.5), (np.random.rand()-0.5), (np.random.rand()-0.5)] crcl[j, ] = np.cross(self.N[i, ], B) crcl[j, ] = crcl[j, ]/np.sqrt(crcl[j, 0]**2 + crcl[j, 1]**2 + crcl[j, 2]**2)*self.r crclx = crcl[:, 0] + self.x[i] crcly = crcl[:, 1] + self.y[i] crclz = crcl[:, 2] + self.z[i] return crclx, crcly, crclz def create_sqr(i): """make circle in which square contact is circumscribed""" sqr = np.zeros((self.n, 3)) for j in range(self.n): B = [(np.random.rand() - 0.5), (np.random.rand() - 0.5), (np.random.rand() - 0.5)] sqr[j,] = np.cross(self.N[i,], B)/np.linalg.norm(np.cross(self.N[i,], B)) * self.r * np.sqrt(2)/2 sqrx = sqr[:, 0] + self.x[i] sqry = sqr[:, 1] + self.y[i] sqrz = sqr[:, 2] + self.z[i] return sqrx, sqry, sqrz def calc_xyz_n(i): """calculate some offsets""" #offsets and radii init offs = np.zeros((self.n, 3)) r2 = np.zeros(self.n) #assert the same random numbers are drawn every time if self.seedvalue is not None: np.random.seed(self.seedvalue) if self.contact_shape is 'circle': for j in range(self.n): A = [(np.random.rand()-0.5)*self.r*2, (np.random.rand()-0.5)*self.r*2, (np.random.rand()-0.5)*self.r*2] offs[j, ] = np.cross(self.N[i, ], A) r2[j] = offs[j, 0]**2 + offs[j, 1]**2 + offs[j, 2]**2 while r2[j] > self.r**2: A = [(np.random.rand()-0.5)*self.r*2, (np.random.rand()-0.5)*self.r*2, (np.random.rand()-0.5)*self.r*2] offs[j, ] = np.cross(self.N[i, ], A) r2[j] = offs[j, 0]**2 + offs[j, 1]**2 + offs[j, 2]**2 elif self.contact_shape is 'square': for j in range(self.n): A = [(np.random.rand()-0.5), (np.random.rand()-0.5), (np.random.rand()-0.5)] offs[j, ] = np.cross(self.N[i, ], A)*self.r r2[j] = offs[j, 0]**2 + offs[j, 1]**2 + offs[j, 2]**2 x_n = offs[:, 0] + self.x[i] y_n = offs[:, 1] + self.y[i] z_n = offs[:, 2] + self.z[i] return x_n, y_n, z_n def loop_over_points(x_n, y_n, z_n): #loop over points on contact for j in range(self.n): tmp = self.lfp_method(self.cell, x = x_n[j], y = y_n[j], z = z_n[j], r_limit = self.r_limit, sigma = self.sigma, **kwargs ) if j == 0: lfp_e = tmp else: lfp_e += tmp #no longer needed del tmp return lfp_e / self.n #loop over contacts for i in range(len(self.x)): if self.n > 1: #fetch offsets: x_n, y_n, z_n = calc_xyz_n(i) #fill in with contact average self.mapping[i] = loop_over_points(x_n, y_n, z_n) #lfp_e.mean(axis=0) else: self.mapping[i] = self.lfp_method(self.cell, x=self.x[i], y=self.y[i], z=self.z[i], r_limit = self.r_limit, sigma=self.sigma, **kwargs) self.offsets[i] = {'x_n' : x_n, 'y_n' : y_n, 'z_n' : z_n} #fetch circumscribed circle around contact if self.contact_shape is 'circle': crcl = create_crcl(i) self.circle_circ[i] = { 'x' : crcl[0], 'y' : crcl[1], 'z' : crcl[2], } elif self.contact_shape is 'square': sqr = create_sqr(i) self.circle_circ[i] = { 'x': sqr[0], 'y': sqr[1], 'z': sqr[2], }
class RecMEAElectrode(RecExtElectrode): """class RecMEAElectrode Electrode class that represents an extracellular in vitro slice recording as a Microelectrode Array (MEA). Inherits RecExtElectrode class Set-up: Above neural tissue (Saline) -> sigma_S <----------------------------------------------------> z = z_shift + h Neural Tissue -> sigma_T o -> source_pos = [x',y',z'] <-----------*----------------------------------------> z = z_shift + 0 \-> elec_pos = [x,y,z] Below neural tissue (MEA Glass plate) -> sigma_G Parameters ---------- cell : None or object If not None, instantiation of LFPy.Cell, LFPy.TemplateCell or similar. sigma_T : float extracellular conductivity of neural tissue in unit (S/m) sigma_S : float conductivity of saline bath that the neural slice is immersed in [1.5] (S/m) sigma_G : float conductivity of MEA glass electrode plate. Most commonly assumed non-conducting [0.0] (S/m) h : float, int Thickness in um of neural tissue layer containing current the current sources (i.e., in vitro slice or cortex) z_shift : float, int Height in um of neural tissue layer bottom. If e.g., top of neural tissue layer should be z=0, use z_shift=-h. Defaults to z_shift = 0, so that the neural tissue layer extends from z=0 to z=h. squeeze_cell_factor : float or None Factor to squeeze the cell in the z-direction. This is needed for large cells that are thicker than the slice, since no part of the cell is allowed to be outside the slice. The squeeze is done after the neural simulation, and therefore does not affect neuronal simulation, only calculation of extracellular potentials. x, y, z : np.ndarray coordinates or arrays of coordinates in units of (um). Must be same length N : None or list of lists Normal vectors [x, y, z] of each circular electrode contact surface, default None r : float radius of each contact surface, default None n : int if N is not None and r > 0, the number of discrete points used to compute the n-point average potential on each circular contact point. contact_shape : str 'circle'/'square' (default 'circle') defines the contact point shape If 'circle' r is the radius, if 'square' r is the side length method : str switch between the assumption of 'linesource', 'pointsource', 'soma_as_point' to represent each compartment when computing extracellular potentials from_file : bool if True, load cell object from file cellfile : str path to cell pickle verbose : bool Flag for verbose output, i.e., print more information seedvalue : int random seed when finding random position on contact with r > 0 Examples See also examples/example_MEA.py >>> import numpy as np >>> import matplotlib.pyplot as plt >>> import LFPy >>> >>> cellParameters = { >>> 'morphology' : 'examples/morphologies/L5_Mainen96_LFPy.hoc', # morphology file >>> 'v_init' : -65, # initial voltage >>> 'cm' : 1.0, # membrane capacitance >>> 'Ra' : 150, # axial resistivity >>> 'passive' : True, # insert passive channels >>> 'passive_parameters' : {"g_pas":1./3E4, "e_pas":-65}, # passive params >>> 'dt' : 2**-4, # simulation time res >>> 'tstart' : 0., # start t of simulation >>> 'tstop' : 50., # end t of simulation >>> } >>> cell = LFPy.Cell(**cellParameters) >>> cell.set_rotation(x=np.pi/2, z=np.pi/2) >>> cell.set_pos(z=100) >>> synapseParameters = { >>> 'idx' : cell.get_closest_idx(x=800, y=0, z=100), # compartment >>> 'e' : 0, # reversal potential >>> 'syntype' : 'ExpSyn', # synapse type >>> 'tau' : 2, # syn. time constant >>> 'weight' : 0.01, # syn. weight >>> 'record_current' : True # syn. current record >>> } >>> synapse = LFPy.Synapse(cell, **synapseParameters) >>> synapse.set_spike_times(np.array([10., 15., 20., 25.])) >>> >>> MEA_electrode_parameters = { >>> 'sigma_T' : 0.3, # extracellular conductivity >>> 'sigma_G' : 0.0, # MEA glass electrode plate conductivity >>> 'sigma_S' : 1.5, # Saline bath conductivity >>> 'x' : np.linspace(0, 1200, 16), # electrode requires 1d vector of positions >>> 'y' : np.zeros(16), >>> 'z' : np.zeros(16), >>> "method": "pointsource", >>> "h": 300, >>> "squeeze_cell_factor": 0.3, >>> } >>> MEA = LFPy.RecMEAElectrode(cell, **MEA_electrode_parameters) >>> >>> cell.simulate(electrode=MEA) >>> >>> plt.matshow(MEA.LFP) >>> plt.colorbar() >>> plt.axis('tight') >>> plt.show() """ def __init__(self, cell=None, sigma_T=0.3, sigma_S=1.5, sigma_G=0.0, h=300., z_shift=0., steps=20, x=np.array([0]), y=np.array([0]), z=np.array([0]), N=None, r=None, n=None, perCellLFP=False, method='linesource', from_file=False, cellfile=None, verbose=False, seedvalue=None, squeeze_cell_factor=None, **kwargs): RecExtElectrode.__init__(self, cell=cell, x=x, y=y, z=z, N=N, r=r, n=n, perCellLFP=perCellLFP, method=method, from_file=from_file, cellfile=cellfile, verbose=verbose, seedvalue=seedvalue, **kwargs) self.sigma_G = sigma_G self.sigma_T = sigma_T self.sigma_S = sigma_S self.sigma = None self.h = h self.z_shift = z_shift self.steps = steps self.squeeze_cell_factor = squeeze_cell_factor self.moi_param_kwargs = {"h": self.h, "steps": self.steps, "sigma_G": self.sigma_G, "sigma_T": self.sigma_T, "sigma_S": self.sigma_S, } if cell is not None: self.set_cell(cell) if method == 'pointsource': self.lfp_method = lfpcalc.calc_lfp_pointsource_moi elif method == "linesource": if (np.abs(z - self.z_shift) > 1e-9).any(): raise NotImplementedError("The method 'linesource' is only " "supported for electrodes at the " "z=0 plane. Use z=0 or method " "'pointsource'.") if np.abs(self.sigma_G) > 1e-9: raise NotImplementedError("The method 'linesource' is only " "supported for sigma_G=0. Use " "sigma_G=0 or method " "'pointsource'.") self.lfp_method = lfpcalc.calc_lfp_linesource_moi elif method == "soma_as_point": if (np.abs(z - self.z_shift) > 1e-9).any(): raise NotImplementedError("The method 'soma_as_point' is only " "supported for electrodes at the " "z=0 plane. Use z=0 or method " "'pointsource'.") if np.abs(self.sigma_G) > 1e-9: raise NotImplementedError("The method 'soma_as_point' is only " "supported for sigma_G=0. Use " "sigma_G=0 or method " "'pointsource'.") self.lfp_method = lfpcalc.calc_lfp_soma_as_point_moi else: raise ValueError("LFP method not recognized. " "Should be 'soma_as_point', 'linesource' " "or 'pointsource'") def _squeeze_cell_in_depth_direction(self): """Will squeeze self.cell centered around the soma by a scaling factor, so that it fits inside the slice. If scaling factor is not big enough, a RuntimeError is raised. """ self.cell.distort_geometry(factor=self.squeeze_cell_factor) if (np.max([self.cell.zstart, self.cell.zend]) > self.h + self.z_shift or np.min([self.cell.zstart, self.cell.zend]) < self.z_shift): bad_comps, reason = self._return_comp_outside_slice() msg = ("Compartments {} of cell ({}) has cell.{} slice. " "Increase squeeze_cell_factor, move or rotate cell." ).format(bad_comps, self.cell.morphology, reason) raise RuntimeError(msg) def _return_comp_outside_slice(self): """ Assuming part of the cell is outside the valid region, i.e, not in the slice (self.z_shift < z < self.z_shift + self.h) this function check what array (cell.zstart or cell.zend) that is outside, and if it is above or below the valid region. Raises: RuntimeError If no compartment is outside valid region. Returns: array, str Numpy array with the compartments that are outside the slice, and a string with additional information on the problem. """ zstart_above = np.where(self.cell.zstart > self.z_shift + self.h)[0] zend_above = np.where(self.cell.zend > self.z_shift + self.h)[0] zend_below = np.where(self.cell.zend < self.z_shift)[0] zstart_below = np.where(self.cell.zstart < self.z_shift)[0] if len(zstart_above) > 0: return zstart_above, "zstart above" if len(zstart_below) > 0: return zstart_below, "zstart below" if len(zend_above) > 0: return zend_above, "zend above" if len(zend_below) > 0: return zend_below, "zend below" raise RuntimeError("This function should only be called if cell" "extends outside slice") def test_cell_extent(self): """ Test if the cell is confined within the slice. If class argument "squeeze_cell" is True, cell is squeezed to to fit inside slice. """ if self.cell is None: raise RuntimeError("Does not have cell instance.") if (np.max([self.cell.zstart, self.cell.zend]) > self.z_shift + self.h or np.min([self.cell.zstart, self.cell.zend]) < self.z_shift): if self.verbose: print("Cell extends outside slice.") if self.squeeze_cell_factor is not None: if not self.z_shift < self.cell.zmid[0] < self.z_shift + self.h: raise RuntimeError("Soma position is not in slice.") self._squeeze_cell_in_depth_direction() else: bad_comps, reason = self._return_comp_outside_slice() msg = ("Compartments {} of cell ({}) has cell.{} slice " "and argument squeeze_cell_factor is None." ).format(bad_comps, self.cell.morphology, reason) raise RuntimeError(msg) else: if self.verbose: print("Cell position is good.") if self.squeeze_cell_factor is not None: if self.verbose: print("Squeezing cell anyway.") self._squeeze_cell_in_depth_direction() def calc_mapping(self, cell): """Creates a linear mapping of transmembrane currents of each segment of the supplied cell object to contribution to extracellular potential at each electrode contact point of the RexExtElectrode object. Sets the class attribute "mapping", which is a shape (n_contact, n_segs) ndarray, such that the extracellular potential at the contacts phi = np.dot(mapping, I_mem) where I_mem is a shape (n_segs, n_tsteps) ndarray with transmembrane currents for each time step of the simulation. Parameters ---------- cell : obj `LFPy.Cell` or `LFPy.TemplateCell` instance. Returns ------- None """ if cell is not None: self.set_cell(cell) self.test_cell_extent() # Temporarily shift coordinate system so middle layer extends # from z=0 to z=h self.z = self.z - self.z_shift self.cell.zstart = self.cell.zstart - self.z_shift self.cell.zmid = self.cell.zmid - self.z_shift self.cell.zend = self.cell.zend - self.z_shift if self.n is not None and self.N is not None and self.r is not None: if self.n <= 1: raise ValueError("n = %i must be larger that 1" % self.n) else: pass self._lfp_el_pos_calc_dist(**self.moi_param_kwargs) if self.verbose: print('calculations finished, %s, %s' % (str(self), str(self.cell))) else: self._loop_over_contacts(**self.moi_param_kwargs) if self.verbose: print('calculations finished, %s, %s' % (str(self), str(self.cell))) # Shift coordinate system back so middle layer extends # from z=z_shift to z=z_shift + h self.z = self.z + self.z_shift self.cell.zstart = self.cell.zstart + self.z_shift self.cell.zmid = self.cell.zmid + self.z_shift self.cell.zend = self.cell.zend + self.z_shift def calc_lfp(self, t_indices=None, cell=None): """Calculate LFP on electrode geometry from all cell instances. Will chose distributed calculated if electrode contain 'n', 'N', and 'r' Parameters ---------- cell : obj, optional `LFPy.Cell` or `LFPy.TemplateCell` instance. Must be specified here if it was not specified at the initiation of the `RecExtElectrode` class t_indices : np.ndarray Array of timestep indexes where extracellular potential should be calculated. """ self.calc_mapping(cell) if t_indices is not None: currmem = self.cell.imem[:, t_indices] else: currmem = self.cell.imem self._test_imem_sum() self.LFP = np.dot(self.mapping, currmem) # del self.mapping