Source code for LFPy.recextelectrode

#!/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
GNU General Public License for more details.'''

import numpy as np
import warnings
from LFPy import lfpcalc, tools

[docs]class RecExtElectrodeSetup(object): ''' RecExtElectrode superclass. If (optional) cell argument is given then the it is imported, otherwise the cell argument has to be passed later on to calc_lfp. The argument cell can be an LFPy.cell.Cell or LFPy.templatecell.TemplateCell object Keyword arguments determine properties of later LFP-calculations Arguments: :: cell : object, LFPy.cell.Cell or LFPy.templatecell.TemplateCell sigma : float, extracellular conductivity x, y, z : np.ndarray, coordinates or arrays of coordinates. Must be same length N : np.ndarray, Normal vector [x, y, z] of contact surface, default None r : float, radius of contact surface, default None n : int, if N is not None and r > 0, the number of points to use for each contact point in order to calculate average 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, ['linesource']/'pointsource'/'som_as_point' switch color : str, color of electrode contact points in plots marker : str, marker of electrode contact points in plots from_file : Bool, if True, load cell object from file cellfile : str, path to cell pickle verbose : Bool, Flag for verbose output seedvalue : int, rand seed when finding random position on contact with r >0 ''' 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, shape=None, r_z=None, perCellLFP=False, method='linesource', color='g', marker='o', from_file=False, cellfile=None, verbose=False, seedvalue=None, **kwargs): '''Initialize class RecExtElectrodeSetup''' self.cell = cell self.sigma = sigma if type(x) is float or type(x) is int: self.x = np.array([x]) else: self.x = np.array(x).flatten() if type(y) is float or type(y) is int: self.y = np.array([y]) else: self.y = np.array(y).flatten() if type(z) is float or type(z) is 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 as ae: raise ae, "The number of elements in [x, y, z] must be identical" self.color = color self.marker = marker 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 shape is None: self.shape = 'circle' elif shape in ['circle', 'square']: self.shape = shape else: raise ValueError('The shape argument must be either: None, \'circle\', \'square\'') self.r_z = r_z 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 = 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 self.cell is not None: self._test_imem_sum() 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') if self.cell is not None: sum_imem = self.cell.imem.sum(axis=0) #check if eye matrix is supplied: if np.any(sum_imem == np.ones(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 else: pass
[docs]class RecExtElectrode(RecExtElectrodeSetup): ''' RecExtElectrode class with inheritance from LFPy.RecExtElectrodeSetup able to actually calculate local field potentials from LFPy.Cell objects. Usage: :: import numpy as np import import matplotlib.pyplot as plt import LFPy 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, } cellParameters = { 'morphology' : 'L5_Mainen96_LFPy.hoc', # morphology file 'rm' : 30000, # membrane resistivity 'cm' : 1.0, # membrane capacitance 'Ra' : 150, # axial resistivity 'timeres_NEURON' : 2**-4, # dt for NEURON sim. 'timeres_python' : 2**-4, # dt for python output 'tstartms' : -50, # start t of simulation 'tstopms' : 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.PointProcessSynapse(cell, **synapseParameters) synapse.set_spike_times(cell, np.array([10, 15, 20, 25])) cell.simulate() electrode = LFPy.RecExtElectrode(cell, **electrodeParameters) electrode.calc_lfp() plt.matshow(electrode.LFP) ''' 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=0, shape=None, r_z=None, perCellLFP=False, method='linesource', color='g', marker='o', from_file=False, cellfile=None, verbose=False, seedvalue=None, **kwargs): '''This is the regular implementation of the RecExtElectrode class that calculates the LFP serially using a single core Arguments: :: cell : LFPy.Cell like object sigma : ''' RecExtElectrodeSetup.__init__(self, cell, sigma, x, y, z, N, r, n, shape, r_z, perCellLFP, method, color, marker, from_file, cellfile, verbose, seedvalue, **kwargs)
[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' ''' if cell is not None: self.cell = cell self._test_imem_sum() if not hasattr(self, 'LFP'): if t_indices is not None: self.LFP = np.zeros((self.x.size, t_indices.size)) else: self.LFP = np.zeros((self.x.size, self.cell.imem.shape[1])) 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.circle_circum, self.offsets, LFP_temp] = \ self._lfp_el_pos_calc_dist(t_indices=t_indices, r_limit=self.cell.diam/2) if self.verbose: print('calculations finished, %s, %s' % (str(self), str(self.cell))) else: LFP_temp = self._loop_over_contacts(t_indices=t_indices, r_limit=self.cell.diam/2) if self.verbose: print('calculations finished, %s, %s' % (str(self), str(self.cell))) #dump results: self.LFP = LFP_temp
def _loop_over_contacts(self, r_limit=None, timestep=None, t_indices=None): '''Loop over electrode contacts, and will return LFPs across channels''' if t_indices is not None: LFP_temp = np.zeros((self.x.size, t_indices.size)) else: LFP_temp = np.zeros((self.x.size, self.cell.imem.shape[1])) for i in range(self.x.size): LFP_temp[i, :] = LFP_temp[i, :] + \ lfpcalc.calc_lfp_choose(self.cell, x = self.x[i], y = self.y[i], z = self.z[i], sigma = self.sigma, r_limit = r_limit, timestep = timestep, t_indices = t_indices, method = self.method, **self.kwargs) return LFP_temp def _lfp_el_pos_calc_dist(self, r_limit=None, m=50, t_indices=None, ): ''' 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) offsets = {} circle_circ = {} def create_crcl(m, i): '''make circumsize of contact point''' crcl = np.zeros((m, 3)) for j in range(m): 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(m, i): '''make circle in which square contact is circumscribed''' sqr = np.zeros((m, 3)) for j in range(m): 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.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.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 = lfpcalc.calc_lfp_choose(self.cell, x = x_n[j], y = y_n[j], z = z_n[j], r_limit = r_limit, sigma = self.sigma, t_indices = t_indices, method = self.method, **self.kwargs) if j == 0: lfp_e = tmp else: lfp_e = np.r_['0,2', lfp_e, tmp] #no longer needed del tmp return lfp_e.mean(axis=0) #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 lfp_el_pos[i] = loop_over_points(x_n, y_n, z_n) #lfp_e.mean(axis=0) ##no longer needed #del lfp_e else: lfp_el_pos[i] = lfpcalc.calc_lfp_choose(self.cell, x=self.x[i], y=self.y[i], z=self.z[i], r_limit = r_limit, sigma=self.sigma, t_indices=t_indices, **self.kwargs) offsets[i] = { 'x_n' : x_n, 'y_n' : y_n, 'z_n' : z_n } #fetch circumscribed circle around contact if self.shape is 'circle': crcl = create_crcl(m, i) circle_circ[i] = { 'x' : crcl[0], 'y' : crcl[1], 'z' : crcl[2], } elif self.shape is 'square': sqr = create_sqr(m, i) circle_circ[i] = { 'x': sqr[0], 'y': sqr[1], 'z': sqr[2], } return circle_circ, offsets, lfp_el_pos