Source code for allensdk.model.glif.glif_neuron_methods

# Allen Institute Software License - This software license is the 2-clause BSD
# license plus a third clause that prohibits redistribution for commercial
# purposes without further permission.
#
# Copyright 2017. Allen Institute. All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimer in the documentation
# and/or other materials provided with the distribution.
#
# 3. Redistributions for commercial purposes are not permitted without the
# Allen Institute's written permission.
# For purposes of this license, commercial purposes is the incorporation of the
# Allen Institute's software into anything for which you will charge fees or
# other compensation. Contact terms@alleninstitute.org for commercial licensing
# opportunities.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
#
""" The methods in this module are used for configuring dynamics and reset rules for the GlifNeuron.  
For more details on how to use these methods, see :doc:`glif_models`.
"""
import functools
import numpy as np


[docs]class GlifNeuronMethod( object ): """ A simple class to keep track of the name and parameters associated with a neuron method. This class is initialized with a name, function, and parameters to pass to the function. The function then has those passed parameters fixed to a partial function using functools.partial. This class then mimics a function itself using the __call__ convention. Parameters that are not fixed in this way are assumed to be passed into the method when it is called. If the passed parameters contain an argument that is not part of the function signature, an exception will be raised. Parameters ---------- method_name : string A shorthand name that will be used to reference this method in the `GlifNeuron`. method : function A python function to be called when this instance is called. method_params : dict A dictionary mapping function arguments to values for values that should be fixed. """ def __init__(self, method_name, method, method_params): self.name = method_name self.params = method_params self.method = functools.partial(method, **method_params) def __call__(self, *args, **kwargs): """ Defining this method allows an instance to be called like a function """ return self.method(*args, **kwargs)
[docs] def to_dict(self): return { 'name': self.name, 'params': self.params }
[docs] def modify_parameter(self, param, operator): """ Modify a function parameter needs to be modified after initialization. Parameters ---------- param : string the name of the parameter to modify operator : callable a function or lambda that returns the desired modified value Returns ------- type the new value of the variable that was just modified. """ value = operator(self.method.keywords[param]) self.method.keywords[param] = value return value
[docs]def max_of_line_and_const(x,b,c,d): #TODO: move to other library """ Find the maximum of a value and a position on a line Parameters ---------- x: float x position on line 1 c: float slope of line 1 d: float y-intercept of line 1 b: float y-intercept of line 2 Returns ------- float the max of a line value and a constant """ one = b two = c*x+d return np.maximum(one,two)
[docs]def min_of_line_and_zero(x,c,d): #TODO: move to other library """ Find the minimum of a value and a position on a line Parameters ---------- x: float x position on line 1 c: float slope of line 1 d: float y-intercept of line 1 b: float y-intercept of line 2 Returns ------- float the max of a line value and a constant """ one = 0 two = c*x+d return np.minimum(one,two)
[docs]def dynamics_AScurrent_exp(neuron, AScurrents_t0, time_step, spike_time_steps): """ Exponential afterspike current dynamics method takes a current at t0 and returns the current at a time step later. """ return AScurrents_t0*np.exp(-neuron.k*neuron.dt)
[docs]def dynamics_AScurrent_none(neuron, AScurrents_t0, time_step, spike_time_steps): """ This method always returns zeros for the afterspike currents, regardless of input. """ return np.zeros(len(AScurrents_t0))
[docs]def dynamics_voltage_linear_forward_euler(neuron, voltage_t0, AScurrents_t0, inj): """ (TODO) Linear voltage dynamics. """ return voltage_t0 + (inj + np.sum(AScurrents_t0) - neuron.G * neuron.coeffs['G'] * (voltage_t0 - neuron.El)) * neuron.dt / (neuron.C * neuron.coeffs['C'])
[docs]def dynamics_voltage_linear_exact(neuron, voltage_t0, AScurrents_t0, inj): """ (TODO) Linear voltage dynamics. """ C = (neuron.C * neuron.coeffs['C']) I = inj + np.sum(AScurrents_t0) g = neuron.G * neuron.coeffs['G'] tau = g/C N = (I+ g*neuron.El)/C return voltage_t0*np.exp(-neuron.dt*tau) + N*(1-np.exp(-tau*neuron.dt))/tau
[docs]def spike_component_of_threshold_forward_euler(th_t0, b_spike, dt): '''Spike component of threshold modeled as an exponential decay. Implemented here for forward Euler Parameters ---------- th_t0 : float threshold input to function b_spike : float decay constant of exponential dt : float time step ''' b_spike=-b_spike #TODO: this is here because b_spike is always input as positive although it is negative return th_t0 + th_t0*b_spike * dt
[docs]def spike_component_of_threshold_exact(th0, b_spike, t): '''Spike component of threshold modeled as an exponential decay. Implemented here as exact analytical solution. Parameters ---------- th0 : float threshold input to function b_spike : float decay constant of exponential t : float or array time step if used in an Euler setup time if used analytically ''' b_spike=-b_spike return th0*np.exp(b_spike * t)
[docs]def voltage_component_of_threshold_forward_euler(th_t0, v_t0, dt, a_voltage, b_voltage, El): '''Equation 2.1 of Mihalas and Nieber, 2009 implemented for use in forward Euler. Note here all variables are in reference to threshold infinity. Therefore thr_inf is zero here (replaced threshold_inf with 0 in the equation to be verbose). This is done so that th_inf can be optimized without affecting this function. Parameters ---------- th_t0 : float threshold input to function v_t0 : float voltage input to function dt : float time step a_voltage : float constant a b_voltage : float constant b El : float reversal potential ''' return th_t0 + (a_voltage*(v_t0-El)-b_voltage*(th_t0-0))*dt
[docs]def voltage_component_of_threshold_exact(th0, v0, I, t, a_voltage, b_voltage, C, g, El): '''Note this function is the exact formulation; however, dt is used because t0 is the initial time and dt is the time the function is exactly evaluated at. Note: that here, this equation is in reference to th_inf. Therefore th0 is the total threshold-thr_inf (threshold_inf replaced with 0 in the equation to be verbose). This is done so that th_inf can be optimized without affecting this function. Parameters ---------- th0 : float threshold input to function v0 : float voltage input to function I : float total current entering neuron (note if there are after spike currents these must be included in this value) t : float or array time step if used in an Euler setup time if used analytically a_voltage : float constant a b_voltage : float constant b C : float capacitance g : float conductance (1/resistance) El : float reversal potential ''' beta=(I+g*El)/g phi=a_voltage/(b_voltage-g/C) return phi*(v0-beta)*np.exp(-g*t/C)+1/(np.exp(b_voltage*t))*(th0-phi*(v0-beta)- (a_voltage/b_voltage)*(beta-El)-0) +(a_voltage/b_voltage)*(beta-El) +0
[docs]def dynamics_threshold_three_components_exact(neuron, threshold_t0, voltage_t0, AScurrents_t0, inj, a_spike, b_spike, a_voltage, b_voltage): """Analytical solution for threshold dynamics. The threshold will adapt via two mechanisms: 1. a voltage dependent adaptation. 2. a component initiated by a spike which decays as an exponential. These two component are in reference to threshold infinity and are recorded in the neuron's threshold components. The third component refers to th_inf which is added separately as opposed to being included in the voltage component of the threshold as is done in equation 2.1 of Mihalas and Nieber 2009. Threshold infinity is removed for simple optimization. Parameters ---------- neuron : class threshold_t0 : float threshold input to function voltage_t0 : float voltage input to function AScurrents_t0 : vector values of after spike currents inj : float current injected into the neuron """ #TODO: just having the get_threshold_components added an erroneous zero to the beginning of the list if neuron.threshold_components is None: neuron.threshold_components = { 'spike': [], 'voltage': [] } th_spike = 0 th_voltage = 0 else: tcs = neuron.threshold_components th_spike = tcs['spike'][-1] th_voltage = tcs['voltage'][-1] a_voltage = a_voltage * neuron.coeffs['a'] b_voltage = b_voltage * neuron.coeffs['b'] I = inj + np.sum(AScurrents_t0) C = neuron.C * neuron.coeffs['C'] g = neuron.G * neuron.coeffs['G'] voltage_component=voltage_component_of_threshold_exact(th_voltage, voltage_t0, I, neuron.dt, a_voltage, b_voltage, C, g, neuron.El) spike_component = spike_component_of_threshold_exact(th_spike, b_spike, neuron.dt) #------update the voltage and spiking values of the the neuron.append_threshold_components(spike_component, voltage_component) return voltage_component+spike_component+neuron.th_inf * neuron.coeffs['th_inf']
[docs]def dynamics_threshold_spike_component(neuron, threshold_t0, voltage_t0, AScurrents_t0, inj, a_spike, b_spike, a_voltage, b_voltage): """Analytical solution for spike component of threshold. The threshold will adapt via a component initiated by a spike which decays as an exponential. The component is in reference to threshold infinity and are recorded in the neuron's threshold components. The voltage component of the threshold is set to zero in the threshold components because it is zero here The third component refers to th_inf which is added separately as opposed to being included in the voltage component of the threshold as is done in equation 2.1 of Mihalas and Nieber 2009. Threshold infinity is removed for simple optimization. Parameters ---------- neuron : class threshold_t0 : float threshold input to function voltage_t0 : float voltage input to function AScurrents_t0 : vector values of after spike currents inj : float current injected into the neuron """ #TODO: just having the get_threshold_components added an erroneous zero to the beginning of the list if neuron.threshold_components is None: neuron.threshold_components = { 'spike': [], 'voltage': [] } th_spike = 0 th_voltage = 0 else: tcs = neuron.threshold_components th_spike = tcs['spike'][-1] th_voltage = tcs['voltage'][-1] spike_component = spike_component_of_threshold_exact(th_spike, b_spike, neuron.dt) #------update the voltage and spiking values of the the neuron.append_threshold_components(spike_component, 0.0) return spike_component+neuron.th_inf * neuron.coeffs['th_inf']
[docs]def dynamics_threshold_inf(neuron, threshold_t0, voltage_t0, AScurrents_t0, inj): """ Set threshold to the neuron's instantaneous threshold. Parameters ---------- neuron : class threshold_t0 : not used here voltage_t0 : not used here AScurrents_t0 : not used here inj : not used here AScurrents_t0 : not used here inj : not used here """ return neuron.coeffs['th_inf'] * neuron.th_inf
[docs]def reset_AScurrent_sum(neuron, AScurrents_t0, r): """ Reset afterspike currents by adding summed exponentials. Left over currents from last spikes as well as newly initiated currents from current spike. Currents amplitudes in neuron.asc_amp_array need to be the amplitudes advanced though the spike cutting. I.e. In the preprocessor if the after spike currents are calculated via the GLM from spike initiation the amplitude at the time after the spike cutting needs to be calculated and neuron.asc_amp_array needs to be set to this value. Parameters ---------- r : np.ndarray a coefficient vector applied to the afterspike currents """ new_currents=neuron.asc_amp_array * neuron.coeffs['asc_amp_array'] #neuron.asc_amp_array are amplitudes initiating after the spike is cut left_over_currents=AScurrents_t0 * r * np.exp(-(neuron.k * neuron.dt * neuron.spike_cut_length)) #advancing cut currents though the spike return new_currents+left_over_currents
[docs]def reset_AScurrent_none(neuron, AScurrents_t0): """ Reset afterspike currents to zero. """ if np.sum(AScurrents_t0)!=0: raise Exception('You are running a LIF but the AScurrents are not zero!') return np.zeros(len(AScurrents_t0))
[docs]def reset_voltage_v_before(neuron, voltage_t0, a, b): """ Reset voltage to the previous value with a scale and offset applied. Parameters ---------- a : float voltage scale constant b : float voltage offset constant """ return a*(voltage_t0)+b
[docs]def reset_voltage_zero(neuron, voltage_t0): """ Reset voltage to zero. """ return 0.0
[docs]def reset_threshold_inf(neuron, threshold_t0, voltage_v1): """ Reset the threshold to instantaneous threshold. """ return neuron.coeffs['th_inf'] * neuron.th_inf
[docs]def reset_threshold_three_components(neuron, threshold_t0, voltage_v1, a_spike, b_spike): '''This method calculates the two components of the threshold: a spike (fast) component and a voltage (slow) component. The threshold_components vectors are then updated so that the traces match the voltage, current, and total threshold traces. The spike component of the threshold decays via an exponential fit specified by the amplitude a_spike and the time constant b_spike fit via the multiblip data. The voltage component does not change during the duration of the spike. The spike component are threshold component are summed along with threshold infinity to return the total threshold. Note that in the current implementation a_spike is added to the last value of the threshold_components which means that a_spike is the amplitude after spike cutting (if there is any). Inputs: neuron: class contains attributes of the neuron threshold_t0, voltage_t0: float are not used but are here for consistency with other methods a_spike: float amplitude of the exponential decay of spike component of threshold after spike cutting has been implemented. b_spike: float amplitude of the exponential decay of spike component of threshold Outputs: Returns: float the total threshold which is the sum of the spike component of threshold, the voltage component of threshold and threshold infinity (with it's corresponding coefficient) neuron.threshold_components: dictionary containing a spike: list vector of spiking component of threshold that corresponds to the voltage, current, and total threshold traces b_spike: list vector of voltage component of threshold that corresponds to the voltage, current, and total threshold traces. Note that this function can be changed to use a_spike at the time of the spike and then have the the spike component plus the residual decay thought the spike. There are benefits and drawbacks to this. This potential change would be beneficial as it perhaps makes more biological sense for the threshold to go up at the time of spike if the traces are ever used. Also this would mean that a_spike would not have to be adjusted thought the spike cutting after the multiblip fit. However the current implementation makes sense in that it is similar to how afterspike currents are implemented. ''' if neuron.threshold_components is None: raise Exception('reset should never happen at the beginning of a trace') tcs = neuron.threshold_components #for ease of updating # note that these values are at the indicie of the time of the spike which is the index right after the voltage crosses # threshold since the neuron.threshold_components are updated by the dynamics method which is called before the reset. th_spike=tcs['spike'][-1] #this needs to decay through the spike must be very particular about how many indicies to decay th_voltage= tcs['voltage'][-1] # calculate spike component decay though spike from time =1 (not zero because zero is already in neuron.threshold_components # via the dynamics method) though the end of the spike cutting spike_comp_decay=spike_component_of_threshold_exact(th_spike, b_spike, np.arange(1,neuron.spike_cut_length+1)*neuron.dt) #Note that the plus one is that one needs to know the decay and the inital condition for next starting point #update neuron.threshold_components via pass by reference. [tcs['voltage'].append(value) for value in np.ones(neuron.spike_cut_length)*th_voltage] #note that here I don't need the plus one because I am starting from zero [tcs['spike'].append(value) for value in spike_comp_decay] # add the amplitude of the spike component decay to last value of vector (reseting) tcs['spike'][-1]=tcs['spike'][-1]+a_spike return tcs['spike'][-1] + tcs['voltage'][-1] + neuron.th_inf * neuron.coeffs['th_inf']
#: The METHOD_LIBRARY constant groups dynamics and reset methods by group name (e.g. 'voltage_dynamics_method'). #Those groups assign each method in this file a string name. This is used by the GlifNeuron when initializing #its dynamics and reset methods. METHOD_LIBRARY = { 'AScurrent_dynamics_method': { 'exp': dynamics_AScurrent_exp, 'none': dynamics_AScurrent_none }, 'voltage_dynamics_method': { 'linear_forward_euler': dynamics_voltage_linear_forward_euler }, 'threshold_dynamics_method': { 'spike_component': dynamics_threshold_spike_component, 'inf': dynamics_threshold_inf, 'three_components_exact': dynamics_threshold_three_components_exact }, 'AScurrent_reset_method': { 'sum': reset_AScurrent_sum, 'none': reset_AScurrent_none }, 'voltage_reset_method': { 'v_before': reset_voltage_v_before, 'zero': reset_voltage_zero }, 'threshold_reset_method': { 'inf': reset_threshold_inf, 'three_components': reset_threshold_three_components } }