diff --git a/PySONIC/core/pneuron.py b/PySONIC/core/pneuron.py index ab311c5..e37e51b 100644 --- a/PySONIC/core/pneuron.py +++ b/PySONIC/core/pneuron.py @@ -1,475 +1,474 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-03 11:53:04 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-11-22 17:41:59 +# @Last Modified time: 2018-11-29 21:03:25 import os import time import pickle import abc import numpy as np from scipy.integrate import odeint import pandas as pd from ..postpro import findPeaks from ..constants import * from ..utils import si_format, logger, ESTIM_filecode from ..batches import xlslog class PointNeuron(metaclass=abc.ABCMeta): ''' Abstract class defining the common API (i.e. mandatory attributes and methods) of all subclasses implementing the channels mechanisms of specific point neurons. The mandatory attributes are: - **name**: a string defining the name of the mechanism. - **Cm0**: a float defining the membrane resting capacitance (in F/m2) - **Vm0**: a float defining the membrane resting potential (in mV) - **states_names**: a list of strings defining the names of the different state probabilities governing the channels behaviour (i.e. the differential HH variables). - **states0**: a 1D array of floats (NOT integers !!!) defining the initial values of the different state probabilities. - **coeff_names**: a list of strings defining the names of the different coefficients to be used in effective simulations. The mandatory methods are: - **currNet**: compute the net ionic current density (in mA/m2) across the membrane, given a specific membrane potential (in mV) and channel states. - **steadyStates**: compute the channels steady-state values for a specific membrane potential value (in mV). - **derStates**: compute the derivatives of channel states, given a specific membrane potential (in mV) and channel states. This method must return a list of derivatives ordered identically as in the states0 attribute. - **getEffRates**: get the effective rate constants of ion channels to be used in effective simulations. This method must return an array of effective rates ordered identically as in the coeff_names attribute. - **derStatesEff**: compute the effective derivatives of channel states, based on 1-dimensional linear interpolators of "effective" coefficients. This method must return a list of derivatives ordered identically as in the states0 attribute. ''' def __repr__(self): return self.__class__.__name__ def pprint(self): return '{} neuron'.format(self.__class__.__name__) @property @abc.abstractmethod def name(self): return 'Should never reach here' @property @abc.abstractmethod def Cm0(self): return 'Should never reach here' @property @abc.abstractmethod def Vm0(self): return 'Should never reach here' @abc.abstractmethod def currNet(self, Vm, states): ''' Compute the net ionic current per unit area. :param Vm: membrane potential (mV) :states: state probabilities of the ion channels :return: current per unit area (mA/m2) ''' @abc.abstractmethod def steadyStates(self, Vm): ''' Compute the channels steady-state values for a specific membrane potential value. :param Vm: membrane potential (mV) :return: array of steady-states ''' @abc.abstractmethod def derStates(self, Vm, states): ''' Compute the derivatives of channel states. :param Vm: membrane potential (mV) :states: state probabilities of the ion channels :return: current per unit area (mA/m2) ''' @abc.abstractmethod def getEffRates(self, Vm): ''' Get the effective rate constants of ion channels, averaged along an acoustic cycle, for future use in effective simulations. :param Vm: array of membrane potential values for an acoustic cycle (mV) :return: an array of rate average constants (s-1) ''' @abc.abstractmethod def derStatesEff(self, Qm, states, interp_data): ''' Compute the effective derivatives of channel states, based on 1-dimensional linear interpolation of "effective" coefficients that summarize the system's behaviour over an acoustic cycle. :param Qm: membrane charge density (C/m2) :states: state probabilities of the ion channels :param interp_data: dictionary of 1D vectors of "effective" coefficients over the charge domain, for specific frequency and amplitude values. ''' def getGates(self): ''' Retrieve the names of the neuron's states that match an ion channel gating. ''' gates = [] for x in self.states_names: if 'alpha{}'.format(x.lower()) in self.coeff_names: gates.append(x) return gates def getRates(self, Vm): ''' Compute the ion channels rate constants for a given membrane potential. :param Vm: membrane potential (mV) :return: a dictionary of rate constants and their values at the given potential. ''' rates = {} for x in self.getGates(): x = x.lower() alpha_str, beta_str = ['{}{}'.format(s, x.lower()) for s in ['alpha', 'beta']] inf_str, tau_str = ['{}inf'.format(x.lower()), 'tau{}'.format(x.lower())] if hasattr(self, 'alpha{}'.format(x)): alphax = getattr(self, alpha_str)(Vm) betax = getattr(self, beta_str)(Vm) elif hasattr(self, '{}inf'.format(x)): xinf = getattr(self, inf_str)(Vm) taux = getattr(self, tau_str)(Vm) alphax = xinf / taux betax = 1 / taux - alphax rates[alpha_str] = alphax rates[beta_str] = betax return rates def vtrap(self, x, y): ''' Generic function used to compute rate constants. ''' return x / (np.exp(x / y) - 1) def Vderivatives(self, y, t, Iinj): ''' Compute the derivatives of a V-cast HH system for a specific value of injected current. :param y: vector of HH system variables at time t :param t: time value (s, unused) :param Iinj: injected current (mA/m2) :return: vector of HH system derivatives at time t ''' - Vm, *states = y Iionic = self.currNet(Vm, states) # mA/m2 dVmdt = (- Iionic + Iinj) / self.Cm0 # mV/s dstates = self.derStates(Vm, states) return [dVmdt, *dstates] def Qderivatives(self, y, t, Cm=None): ''' Compute the derivatives of the n-ODE HH system variables, based on a value of membrane capacitance. :param y: vector of HH system variables at time t :param t: specific instant in time (s) :param Cm: membrane capacitance (F/m2) :return: vector of HH system derivatives at time t ''' if Cm is None: Cm = self.Cm0 Qm, *states = y Vm = Qm / Cm * 1e3 # mV dQm = - self.currNet(Vm, states) * 1e-3 # A/m2 dstates = self.derStates(Vm, states) return [dQm, *dstates] def checkInputs(self, Astim, tstim, toffset, PRF, DC): ''' Check validity of electrical stimulation parameters. :param Astim: pulse amplitude (mA/m2) :param tstim: pulse duration (s) :param toffset: offset duration (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) ''' # Check validity of stimulation parameters if not all(isinstance(param, float) for param in [Astim, tstim, toffset, DC]): raise TypeError('Invalid stimulation parameters (must be float typed)') if tstim <= 0: raise ValueError('Invalid stimulus duration: {} ms (must be strictly positive)' .format(tstim * 1e3)) if toffset < 0: raise ValueError('Invalid stimulus offset: {} ms (must be positive or null)' .format(toffset * 1e3)) if DC <= 0.0 or DC > 1.0: raise ValueError('Invalid duty cycle: {} (must be within ]0; 1])'.format(DC)) if DC < 1.0: if not isinstance(PRF, float): raise TypeError('Invalid PRF value (must be float typed)') if PRF is None: raise AttributeError('Missing PRF value (must be provided when DC < 1)') if PRF < 1 / tstim: raise ValueError('Invalid PRF: {} Hz (PR interval exceeds stimulus duration)' .format(PRF)) def simulate(self, Astim, tstim, toffset, PRF=None, DC=1.0): ''' Compute solutions of a neuron's HH system for a specific set of electrical stimulation parameters, using a classic integration scheme. :param Astim: pulse amplitude (mA/m2) :param tstim: pulse duration (s) :param toffset: offset duration (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) :return: 3-tuple with the time profile and solution matrix and a state vector ''' # Check validity of stimulation parameters self.checkInputs(Astim, tstim, toffset, PRF, DC) # Determine system time step dt = DT_ESTIM # if CW stimulus: divide integration during stimulus into single interval if DC == 1.0: PRF = 1 / tstim # Compute vector sizes npulses = int(np.round(PRF * tstim)) Tpulse_on = DC / PRF Tpulse_off = (1 - DC) / PRF # For high-PRF pulsed protocols: adapt time step to ensure minimal # number of samples during TON or TOFF dt_warning_msg = 'high-PRF protocol: lowering time step to %.2e s to properly integrate %s' for key, Tpulse in {'TON': Tpulse_on, 'TOFF': Tpulse_off}.items(): if Tpulse > 0 and Tpulse / dt < MIN_SAMPLES_PER_PULSE_INT: dt = Tpulse / MIN_SAMPLES_PER_PULSE_INT logger.warning(dt_warning_msg, dt, key) n_pulse_on = int(np.round(Tpulse_on / dt)) n_pulse_off = int(np.round(Tpulse_off / dt)) # Compute offset size n_off = int(np.round(toffset / dt)) # Set initial conditions y0 = [self.Vm0, *self.states0] nvar = len(y0) # Initialize global arrays t = np.array([0.]) states = np.array([1]) y = np.array([y0]).T # Initialize pulse time and states vectors t_pulse0 = np.linspace(0, Tpulse_on + Tpulse_off, n_pulse_on + n_pulse_off) states_pulse = np.concatenate((np.ones(n_pulse_on), np.zeros(n_pulse_off))) # Loop through all pulse (ON and OFF) intervals for i in range(npulses): # Construct and initialize arrays t_pulse = t_pulse0 + t[-1] y_pulse = np.empty((nvar, n_pulse_on + n_pulse_off)) # Integrate ON system y_pulse[:, :n_pulse_on] = odeint( self.Vderivatives, y[:, -1], t_pulse[:n_pulse_on], args=(Astim,)).T # Integrate OFF system if n_pulse_off > 0: y_pulse[:, n_pulse_on:] = odeint( self.Vderivatives, y_pulse[:, n_pulse_on - 1], t_pulse[n_pulse_on:], args=(0.0,)).T # Append pulse arrays to global arrays states = np.concatenate([states, states_pulse[1:]]) t = np.concatenate([t, t_pulse[1:]]) y = np.concatenate([y, y_pulse[:, 1:]], axis=1) # Integrate offset interval if n_off > 0: t_off = np.linspace(0, toffset, n_off) + t[-1] states_off = np.zeros(n_off) y_off = odeint(self.Vderivatives, y[:, -1], t_off, args=(0.0, )).T # Concatenate offset arrays to global arrays states = np.concatenate([states, states_off[1:]]) t = np.concatenate([t, t_off[1:]]) y = np.concatenate([y, y_off[:, 1:]], axis=1) # Return output variables return (t, y, states) def titrate(self, tstim, toffset, PRF=None, DC=1.0, Arange=(0., 2 * TITRATION_ESTIM_A_MAX)): ''' Use a dichotomic recursive search to determine the threshold amplitude needed to obtain neural excitation for a given duration, PRF and duty cycle. :param tstim: duration of US stimulation (s) :param toffset: duration of the offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) :param Arange: search interval for Astim, iteratively refined :return: 5-tuple with the determined threshold, time profile, solution matrix, state vector and response latency ''' Astim = (Arange[0] + Arange[1]) / 2 # Run simulation and detect spikes t0 = time.time() (t, y, states) = self.simulate(Astim, tstim, toffset, PRF, DC) tcomp = time.time() - t0 dt = t[1] - t[0] ipeaks, *_ = findPeaks(y[0, :], SPIKE_MIN_VAMP, int(np.ceil(SPIKE_MIN_DT / dt)), SPIKE_MIN_VPROM) nspikes = ipeaks.size latency = t[ipeaks[0]] if nspikes > 0 else None logger.debug('A = %sA/m2 ---> %s spike%s detected', si_format(Astim * 1e-3, 2, space=' '), nspikes, "s" if nspikes > 1 else "") # If accurate threshold is found, return simulation results if (Arange[1] - Arange[0]) <= TITRATION_ESTIM_DA_MAX and nspikes == 1: return (Astim, t, y, states, latency, tcomp) # Otherwise, refine titration interval and iterate recursively else: if nspikes == 0: # if Astim too close to max then stop if (TITRATION_ESTIM_A_MAX - Astim) <= TITRATION_ESTIM_DA_MAX: return (np.nan, t, y, states, latency, tcomp) Arange = (Astim, Arange[1]) else: Arange = (Arange[0], Astim) return self.titrate(tstim, toffset, PRF, DC, Arange=Arange) def runAndSave(self, outdir, tstim, toffset, PRF=None, DC=1.0, Astim=None): ''' Run a simulation of the point-neuron Hodgkin-Huxley system with specific parameters, and save the results in a PKL file. :param outdir: full path to output directory :param tstim: stimulus duration (s) :param toffset: stimulus offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: stimulus duty cycle (-) :param Astim: stimulus amplitude (mA/m2) ''' # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") if Astim is not None: logger.info('%s: simulation @ A = %sA/m2, t = %ss (%ss offset)%s', self, si_format(Astim * 1e-3, 2, space=' '), *si_format([tstim, toffset], 1, space=' '), (', PRF = {}Hz, DC = {:.2f}%'.format(si_format(PRF, 2, space=' '), DC * 1e2) if DC < 1.0 else '')) # Run simulation tstart = time.time() t, y, states = self.simulate(Astim, tstim, toffset, PRF, DC) Vm, *channels = y tcomp = time.time() - tstart # Detect spikes on Vm signal dt = t[1] - t[0] ipeaks, *_ = findPeaks(Vm, SPIKE_MIN_VAMP, int(np.ceil(SPIKE_MIN_DT / dt)), SPIKE_MIN_VPROM) nspikes = ipeaks.size lat = t[ipeaks[0]] if nspikes > 0 else 'N/A' outstr = '{} spike{} detected'.format(nspikes, 's' if nspikes > 1 else '') else: logger.info('%s: titration @ t = %ss%s', self, si_format(tstim, 1, space=' '), (', PRF = {}Hz, DC = {:.2f}%'.format(si_format(PRF, 2, space=' '), DC * 1e2) if DC < 1.0 else '')) # Run titration Astim, t, y, states, lat, tcomp = self.titrate(tstim, toffset, PRF, DC) Vm, *channels = y nspikes = 1 if Astim is np.nan: outstr = 'no spikes detected within titration interval' nspikes = 0 else: nspikes = 1 outstr = 'Athr = {}A/m2'.format(si_format(Astim * 1e-3, 2, space=' ')) logger.debug('completed in %s, %s', si_format(tcomp, 1), outstr) sr = np.mean(1 / np.diff(t[ipeaks])) if nspikes > 1 else None # Store dataframe and metadata df = pd.DataFrame({ 't': t, 'states': states, 'Vm': Vm, 'Qm': Vm * self.Cm0 * 1e-3 }) for j in range(len(self.states_names)): df[self.states_names[j]] = channels[j] meta = { 'neuron': self.name, 'Astim': Astim, 'tstim': tstim, 'toffset': toffset, 'PRF': PRF, 'DC': DC, 'tcomp': tcomp } # Export into to PKL file simcode = ESTIM_filecode(self.name, Astim, tstim, PRF, DC) outpath = '{}/{}.pkl'.format(outdir, simcode) with open(outpath, 'wb') as fh: pickle.dump({'meta': meta, 'data': df}, fh) logger.debug('simulation data exported to "%s"', outpath) # Export key metrics to log file logpath = os.path.join(outdir, 'log_ESTIM.xlsx') logentry = { 'Date': date_str, 'Time': daytime_str, 'Neuron Type': self.name, 'Astim (mA/m2)': Astim, 'Tstim (ms)': tstim * 1e3, 'PRF (kHz)': PRF * 1e-3 if DC < 1 else 'N/A', 'Duty factor': DC, '# samples': t.size, 'Comp. time (s)': round(tcomp, 2), '# spikes': nspikes, 'Latency (ms)': lat * 1e3 if isinstance(lat, float) else 'N/A', 'Spike rate (sp/ms)': sr * 1e-3 if isinstance(sr, float) else 'N/A' } if xlslog(logpath, logentry) == 1: logger.debug('log exported to "%s"', logpath) else: logger.error('log export to "%s" aborted', self.logpath) return outpath def findRheobaseAmps(self, DCs, Vthr, curr='net'): ''' Find the rheobase amplitudes (i.e. threshold amplitudes of infinite duration that would result in excitation) of a specific neuron for various stimulation duty cycles. :param DCs: duty cycles vector (-) :param Vthr: threshold membrane potential above which the neuron necessarily fires (mV) :return: rheobase amplitudes vector (mA/m2) ''' # Compute the pulse average net (or leakage) current along the amplitude space if curr == 'net': iNet = self.currNet(Vthr, self.steadyStates(Vthr)) elif curr == 'leak': iNet = self.currL(Vthr) # Compute rheobase amplitudes return iNet / np.array(DCs) diff --git a/PySONIC/neurons/__init__.py b/PySONIC/neurons/__init__.py index 4db7a38..4b75c8d 100644 --- a/PySONIC/neurons/__init__.py +++ b/PySONIC/neurons/__init__.py @@ -1,23 +1,24 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-06-06 13:36:00 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-09-25 14:54:59 +# @Last Modified time: 2018-11-29 20:41:39 import inspect import sys from .cortical import CorticalRS, CorticalFS, CorticalLTS, CorticalIB from .thalamic import ThalamicRE, ThalamoCortical from .leech import LeechTouch, LeechPressure, LeechRetzius +from .stn import OtsukaSTN def getNeuronsDict(): current_module = sys.modules[__name__] neurons_dict = {} for _, obj in inspect.getmembers(current_module): if inspect.isclass(obj) and isinstance(obj.name, str): neurons_dict[obj.name] = obj return neurons_dict diff --git a/PySONIC/neurons/cortical.py b/PySONIC/neurons/cortical.py index cbda424..8754472 100644 --- a/PySONIC/neurons/cortical.py +++ b/PySONIC/neurons/cortical.py @@ -1,811 +1,811 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-07-31 15:19:51 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-09-28 14:05:49 +# @Last Modified time: 2018-11-30 10:33:49 import numpy as np from ..core import PointNeuron class Cortical(PointNeuron): ''' Class defining the generic membrane channel dynamics of a cortical neuron with 4 different current types: - Inward Sodium current - Outward, delayed-rectifier Potassium current - Outward, slow non.inactivating Potassium current - Non-specific leakage current This generic class cannot be used directly as it does not contain any specific parameters. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* ''' # Generic biophysical parameters of cortical cells Cm0 = 1e-2 # Cell membrane resting capacitance (F/m2) Vm0 = 0.0 # Dummy value for membrane potential (mV) VNa = 50.0 # Sodium Nernst potential (mV) VK = -90.0 # Potassium Nernst potential (mV) VCa = 120.0 # # Calcium Nernst potential (mV) def __init__(self): ''' Constructor of the class ''' # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 'p'] self.states0 = np.array([]) # Names of the different coefficients to be averaged in a lookup table. self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphap', 'betap'] # Charge interval bounds for lookup creation self.Qbounds = (np.round(self.Vm0 - 25.0) * 1e-5, 50.0e-5) def alpham(self, Vm): ''' Compute the alpha rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = 0.32 * self.vtrap(13 - Vdiff, 4) # ms-1 return alpha * 1e3 # s-1 def betam(self, Vm): ''' Compute the beta rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = 0.28 * self.vtrap(Vdiff - 40, 5) # ms-1 return beta * 1e3 # s-1 def alphah(self, Vm): ''' Compute the alpha rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = (0.128 * np.exp(-(Vdiff - 17) / 18)) # ms-1 return alpha * 1e3 # s-1 def betah(self, Vm): ''' Compute the beta rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = (4 / (1 + np.exp(-(Vdiff - 40) / 5))) # ms-1 return beta * 1e3 # s-1 def alphan(self, Vm): ''' Compute the alpha rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = 0.032 * self.vtrap(15 - Vdiff, 5) # ms-1 return alpha * 1e3 # s-1 def betan(self, Vm): ''' Compute the beta rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = (0.5 * np.exp(-(Vdiff - 10) / 40)) # ms-1 return beta * 1e3 # s-1 def pinf(self, Vm): ''' Compute the asymptotic value of the open-probability of slow non-inactivating Potassium channels. :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1 + np.exp(-(Vm + 35) / 10)) # prob def taup(self, Vm): ''' Compute the decay time constant for adaptation of slow non-inactivating Potassium channels. :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' return self.TauMax / (3.3 * np.exp((Vm + 35) / 20) + np.exp(-(Vm + 35) / 20)) # s def derM(self, Vm, m): ''' Compute the evolution of the open-probability of Sodium channels. :param Vm: membrane potential (mV) :param m: open-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alpham(Vm) * (1 - m) - self.betam(Vm) * m def derH(self, Vm, h): ''' Compute the evolution of the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :param h: inactivation-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphah(Vm) * (1 - h) - self.betah(Vm) * h def derN(self, Vm, n): ''' Compute the evolution of the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :param n: open-probability of delayed-rectifier Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphan(Vm) * (1 - n) - self.betan(Vm) * n def derP(self, Vm, p): ''' Compute the evolution of the open-probability of slow non-inactivating Potassium channels. :param Vm: membrane potential (mV) :param p: open-probability of slow non-inactivating Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.pinf(Vm) - p) / self.taup(Vm) def currNa(self, m, h, Vm): ''' Compute the inward Sodium current per unit area. :param m: open-probability of Sodium channels :param h: inactivation-probability of Sodium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GNa = self.GNaMax * m**3 * h return GNa * (Vm - self.VNa) def currK(self, n, Vm): ''' Compute the outward, delayed-rectifier Potassium current per unit area. :param n: open-probability of delayed-rectifier Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GKMax * n**4 return GK * (Vm - self.VK) def currM(self, p, Vm): ''' Compute the outward, slow non-inactivating Potassium current per unit area. :param p: open-probability of the slow non-inactivating Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GM = self.GMMax * p return GM * (Vm - self.VK) - def currL(self, Vm): + def currLeak(self, Vm): ''' Compute the non-specific leakage current per unit area. :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' - return self.GL * (Vm - self.VL) + return self.GLeak * (Vm - self.VLeak) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, p = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + - self.currM(p, Vm) + self.currL(Vm)) # mA/m2 + self.currM(p, Vm) + self.currLeak(Vm)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) peq = self.pinf(Vm) return np.array([meq, heq, neq, peq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, p = states dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dpdt = self.derP(Vm, p) return [dmdt, dhdt, dndt, dpdt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) Tp = self.taup(Vm) pinf = self.pinf(Vm) ap_avg = np.mean(pinf / Tp) bp_avg = np.mean(1 / Tp) - ap_avg # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, ap_avg, bp_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) m, h, n, p = states dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dpdt = rates[6] * (1 - p) - rates[7] * p return [dmdt, dhdt, dndt, dpdt] class CorticalRS(Cortical): ''' Specific membrane channel dynamics of a cortical regular spiking, excitatory pyramidal neuron. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* ''' # Name of channel mechanism name = 'RS' # Cell-specific biophysical parameters Vm0 = -71.9 # Cell membrane resting potential (mV) GNaMax = 560.0 # Max. conductance of Sodium current (S/m^2) GKMax = 60.0 # Max. conductance of delayed Potassium current (S/m^2) GMMax = 0.75 # Max. conductance of slow non-inactivating Potassium current (S/m^2) - GL = 0.205 # Conductance of non-specific leakage current (S/m^2) - VL = -70.3 # Non-specific leakage Nernst potential (mV) + GLeak = 0.205 # Conductance of non-specific leakage current (S/m^2) + VLeak = -70.3 # Non-specific leakage Nernst potential (mV) VT = -56.2 # Spike threshold adjustment parameter (mV) TauMax = 0.608 # Max. adaptation decay of slow non-inactivating Potassium current (s) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_K\ kin.': ['n'], 'i_M\ kin.': ['p'], - 'I': ['iNa', 'iK', 'iM', 'iL', 'iNet'] + 'I': ['iNa', 'iK', 'iM', 'iLeak', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' # Instantiate parent class super().__init__() # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) class CorticalFS(Cortical): ''' Specific membrane channel dynamics of a cortical fast-spiking, inhibitory neuron. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* ''' # Name of channel mechanism name = 'FS' # Cell-specific biophysical parameters Vm0 = -71.4 # Cell membrane resting potential (mV) GNaMax = 580.0 # Max. conductance of Sodium current (S/m^2) GKMax = 39.0 # Max. conductance of delayed Potassium current (S/m^2) GMMax = 0.787 # Max. conductance of slow non-inactivating Potassium current (S/m^2) - GL = 0.38 # Conductance of non-specific leakage current (S/m^2) - VL = -70.4 # Non-specific leakage Nernst potential (mV) + GLeak = 0.38 # Conductance of non-specific leakage current (S/m^2) + VLeak = -70.4 # Non-specific leakage Nernst potential (mV) VT = -57.9 # Spike threshold adjustment parameter (mV) TauMax = 0.502 # Max. adaptation decay of slow non-inactivating Potassium current (s) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_K\ kin.': ['n'], 'i_M\ kin.': ['p'], - 'I': ['iNa', 'iK', 'iM', 'iL', 'iNet'] + 'I': ['iNa', 'iK', 'iM', 'iLeak', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' # Instantiate parent class super().__init__() # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) class CorticalLTS(Cortical): ''' Specific membrane channel dynamics of a cortical low-threshold spiking, inhibitory neuron with an additional inward Calcium current due to the presence of a T-type channel. References: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* *Huguenard, J.R., and McCormick, D.A. (1992). Simulation of the currents involved in rhythmic oscillations in thalamic relay neurons. J. Neurophysiol. 68, 1373–1383.* ''' # Name of channel mechanism name = 'LTS' # Cell-specific biophysical parameters Vm0 = -54.0 # Cell membrane resting potential (mV) GNaMax = 500.0 # Max. conductance of Sodium current (S/m^2) GKMax = 40.0 # Max. conductance of delayed Potassium current (S/m^2) GMMax = 0.28 # Max. conductance of slow non-inactivating Potassium current (S/m^2) GTMax = 4.0 # Max. conductance of low-threshold Calcium current (S/m^2) - GL = 0.19 # Conductance of non-specific leakage current (S/m^2) - VL = -50.0 # Non-specific leakage Nernst potential (mV) + GLeak = 0.19 # Conductance of non-specific leakage current (S/m^2) + VLeak = -50.0 # Non-specific leakage Nernst potential (mV) VT = -50.0 # Spike threshold adjustment parameter (mV) TauMax = 4.0 # Max. adaptation decay of slow non-inactivating Potassium current (s) Vx = -7.0 # Voltage-dependence uniform shift factor at 36°C (mV) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_K\ kin.': ['n'], 'i_M\ kin.': ['p'], 'i_T\ kin.': ['s', 'u'], - 'I': ['iNa', 'iK', 'iM', 'iT', 'iL', 'iNet'] + 'I': ['iNa', 'iK', 'iM', 'iT', 'iLeak', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' # Instantiate parent class super().__init__() # Add names of cell-specific Calcium channel probabilities self.states_names += ['s', 'u'] # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) # Define the names of the different coefficients to be averaged in a lookup table. self.coeff_names += ['alphas', 'betas', 'alphau', 'betau'] def sinf(self, Vm): ''' Compute the asymptotic value of the open-probability of the S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1.0 + np.exp(-(Vm + self.Vx + 57.0) / 6.2)) # prob def taus(self, Vm): ''' Compute the decay time constant for adaptation of S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' tmp = np.exp(-(Vm + self.Vx + 132.0) / 16.7) + np.exp((Vm + self.Vx + 16.8) / 18.2) return 1.0 / 3.7 * (0.612 + 1.0 / tmp) * 1e-3 # s def uinf(self, Vm): ''' Compute the asymptotic value of the open-probability of the U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1.0 + np.exp((Vm + self.Vx + 81.0) / 4.0)) # prob def tauu(self, Vm): ''' Compute the decay time constant for adaptation of U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' if Vm + self.Vx < -80.0: return 1.0 / 3.7 * np.exp((Vm + self.Vx + 467.0) / 66.6) * 1e-3 # s else: return 1.0 / 3.7 * (np.exp(-(Vm + self.Vx + 22) / 10.5) + 28.0) * 1e-3 # s def derS(self, Vm, s): ''' Compute the evolution of the open-probability of the S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :param s: open-probability of S-type Calcium activation gates (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.sinf(Vm) - s) / self.taus(Vm) def derU(self, Vm, u): ''' Compute the evolution of the open-probability of the U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :param u: open-probability of U-type Calcium inactivation gates (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.uinf(Vm) - u) / self.tauu(Vm) def currCa(self, s, u, Vm): ''' Compute the inward, low-threshold Calcium current per unit area. :param s: open-probability of the S-type activation gate of Calcium channels :param u: open-probability of the U-type inactivation gate of Calcium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GT = self.GTMax * s**2 * u return GT * (Vm - self.VCa) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, p, s, u = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currM(p, Vm) + - self.currCa(s, u, Vm) + self.currL(Vm)) # mA/m2 + self.currCa(s, u, Vm) + self.currLeak(Vm)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Call parent method to compute Sodium and Potassium channels gates steady-states NaK_eqstates = super().steadyStates(Vm) # Compute Calcium channel gates steady-states seq = self.sinf(Vm) ueq = self.uinf(Vm) Ca_eqstates = np.array([seq, ueq]) # Merge all steady-states and return return np.concatenate((NaK_eqstates, Ca_eqstates)) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack input states *NaK_states, s, u = states # Call parent method to compute Sodium and Potassium channels states derivatives NaK_derstates = super().derStates(Vm, NaK_states) # Compute Calcium channels states derivatives dsdt = self.derS(Vm, s) dudt = self.derU(Vm, u) # Merge all states derivatives and return return NaK_derstates + [dsdt, dudt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Call parent method to compute Sodium and Potassium effective rate constants NaK_rates = super().getEffRates(Vm) # Compute Calcium effective rate constants Ts = self.taus(Vm) as_avg = np.mean(self.sinf(Vm) / Ts) bs_avg = np.mean(1 / Ts) - as_avg Tu = np.array([self.tauu(v) for v in Vm]) au_avg = np.mean(self.uinf(Vm) / Tu) bu_avg = np.mean(1 / Tu) - au_avg Ca_rates = np.array([as_avg, bs_avg, au_avg, bu_avg]) # Merge all rates and return return np.concatenate((NaK_rates, Ca_rates)) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' # Unpack input states *NaK_states, s, u = states # Call parent method to compute Sodium and Potassium channels states derivatives NaK_dstates = super().derStatesEff(Qm, NaK_states, interp_data) # Compute Calcium channels states derivatives Ca_rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names[8:]]) dsdt = Ca_rates[0] * (1 - s) - Ca_rates[1] * s dudt = Ca_rates[2] * (1 - u) - Ca_rates[3] * u # Merge all states derivatives and return return NaK_dstates + [dsdt, dudt] class CorticalIB(Cortical): ''' Specific membrane channel dynamics of a cortical intrinsically bursting neuron with an additional inward Calcium current due to the presence of a L-type channel. References: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* *Reuveni, I., Friedman, A., Amitai, Y., and Gutnick, M.J. (1993). Stepwise repolarization from Ca2+ plateaus in neocortical pyramidal cells: evidence for nonhomogeneous distribution of HVA Ca2+ channels in dendrites. J. Neurosci. 13, 4609–4621.* ''' # Name of channel mechanism name = 'IB' # Cell-specific biophysical parameters Vm0 = -71.4 # Cell membrane resting potential (mV) GNaMax = 500 # Max. conductance of Sodium current (S/m^2) GKMax = 50 # Max. conductance of delayed Potassium current (S/m^2) GMMax = 0.3 # Max. conductance of slow non-inactivating Potassium current (S/m^2) GCaLMax = 1.0 # Max. conductance of L-type Calcium current (S/m^2) - GL = 0.1 # Conductance of non-specific leakage current (S/m^2) - VL = -70 # Non-specific leakage Nernst potential (mV) + GLeak = 0.1 # Conductance of non-specific leakage current (S/m^2) + VLeak = -70 # Non-specific leakage Nernst potential (mV) VT = -56.2 # Spike threshold adjustment parameter (mV) TauMax = 0.608 # Max. adaptation decay of slow non-inactivating Potassium current (s) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_K\ kin.': ['n'], 'i_M\ kin.': ['p'], 'i_{CaL}\ kin.': ['q', 'r', 'q2r'], - 'I': ['iNa', 'iK', 'iM', 'iCaL', 'iL', 'iNet'] + 'I': ['iNa', 'iK', 'iM', 'iCaL', 'iLeak', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' # Instantiate parent class super().__init__() # Add names of cell-specific Calcium channel probabilities self.states_names += ['q', 'r'] # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) # Define the names of the different coefficients to be averaged in a lookup table. self.coeff_names += ['alphaq', 'betaq', 'alphar', 'betar'] def alphaq(self, Vm): ''' Compute the alpha rate for the open-probability of L-type Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = 0.055 * self.vtrap(-(Vm + 27), 3.8) # ms-1 return alpha * 1e3 # s-1 def betaq(self, Vm): ''' Compute the beta rate for the open-probability of L-type Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 0.94 * np.exp(-(Vm + 75) / 17) # ms-1 return beta * 1e3 # s-1 def alphar(self, Vm): ''' Compute the alpha rate for the inactivation-probability of L-type Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = 0.000457 * np.exp(-(Vm + 13) / 50) # ms-1 return alpha * 1e3 # s-1 def betar(self, Vm): ''' Compute the beta rate for the inactivation-probability of L-type Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 0.0065 / (np.exp(-(Vm + 15) / 28) + 1) # ms-1 return beta * 1e3 # s-1 def derQ(self, Vm, q): ''' Compute the evolution of the open-probability of the Q (activation) gate of L-type Calcium channels. :param Vm: membrane potential (mV) :param q: open-probability of Q gate (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphaq(Vm) * (1 - q) - self.betaq(Vm) * q def derR(self, Vm, r): ''' Compute the evolution of the open-probability of the R (inactivation) gate of L-type Calcium channels. :param Vm: membrane potential (mV) :param r: open-probability of R gate (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphar(Vm) * (1 - r) - self.betar(Vm) * r def currCaL(self, q, r, Vm): ''' Compute the inward L-type Calcium current per unit area. :param q: open-probability of Q gate (prob) :param r: open-probability of R gate (prob) :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GCaL = self.GCaLMax * q**2 * r return GCaL * (Vm - self.VCa) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, p, q, r = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currM(p, Vm) + - self.currCaL(q, r, Vm) + self.currL(Vm)) # mA/m2 + self.currCaL(q, r, Vm) + self.currLeak(Vm)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Call parent method to compute Sodium and Potassium channels gates steady-states NaK_eqstates = super().steadyStates(Vm) # Compute L-type Calcium channel gates steady-states qeq = self.alphaq(Vm) / (self.alphaq(Vm) + self.betaq(Vm)) req = self.alphar(Vm) / (self.alphar(Vm) + self.betar(Vm)) CaL_eqstates = np.array([qeq, req]) # Merge all steady-states and return return np.concatenate((NaK_eqstates, CaL_eqstates)) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack input states *NaK_states, q, r = states # Call parent method to compute Sodium and Potassium channels states derivatives NaK_derstates = super().derStates(Vm, NaK_states) # Compute L-type Calcium channels states derivatives dqdt = self.derQ(Vm, q) drdt = self.derR(Vm, r) # Merge all states derivatives and return return NaK_derstates + [dqdt, drdt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Call parent method to compute Sodium and Potassium effective rate constants NaK_rates = super().getEffRates(Vm) # Compute Calcium effective rate constants aq_avg = np.mean(self.alphaq(Vm)) bq_avg = np.mean(self.betaq(Vm)) ar_avg = np.mean(self.alphar(Vm)) br_avg = np.mean(self.betar(Vm)) CaL_rates = np.array([aq_avg, bq_avg, ar_avg, br_avg]) # Merge all rates and return return np.concatenate((NaK_rates, CaL_rates)) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' # Unpack input states *NaK_states, q, r = states # Call parent method to compute Sodium and Potassium channels states derivatives NaK_dstates = super().derStatesEff(Qm, NaK_states, interp_data) # Compute Calcium channels states derivatives CaL_rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names[8:]]) dqdt = CaL_rates[0] * (1 - q) - CaL_rates[1] * q drdt = CaL_rates[2] * (1 - r) - CaL_rates[3] * r # Merge all states derivatives and return return NaK_dstates + [dqdt, drdt] diff --git a/PySONIC/neurons/leech.py b/PySONIC/neurons/leech.py index 23d08ef..a25e425 100644 --- a/PySONIC/neurons/leech.py +++ b/PySONIC/neurons/leech.py @@ -1,1076 +1,1076 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-07-31 15:20:54 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-09-28 14:06:06 +# @Last Modified time: 2018-11-30 10:34:22 from functools import partialmethod import numpy as np from ..core import PointNeuron class LeechTouch(PointNeuron): ''' Class defining the membrane channel dynamics of a leech touch sensory neuron. with 4 different current types: - Inward Sodium current - Outward Potassium current - Inward Calcium current - Non-specific leakage current - Calcium-dependent, outward Potassium current - Outward, Sodium pumping current Reference: *Cataldo, E., Brunelli, M., Byrne, J.H., Av-Ron, E., Cai, Y., and Baxter, D.A. (2005). Computational model of touch sensory cells (T Cells) of the leech: role of the afterhyperpolarization (AHP) in activity-dependent conduction failure. J Comput Neurosci 18, 5–24.* ''' # Name of channel mechanism name = 'LeechT' # Cell-specific biophysical parameters Cm0 = 1e-2 # Cell membrane resting capacitance (F/m2) Vm0 = -53.58 # Cell membrane resting potential (mV) VNa = 45.0 # Sodium Nernst potential (mV) VK = -62.0 # Potassium Nernst potential (mV) VCa = 60.0 # Calcium Nernst potential (mV) - VL = -48.0 # Non-specific leakage Nernst potential (mV) + VLeak = -48.0 # Non-specific leakage Nernst potential (mV) VPumpNa = -300.0 # Sodium pump current reversal potential (mV) GNaMax = 3500.0 # Max. conductance of Sodium current (S/m^2) GKMax = 900.0 # Max. conductance of Potassium current (S/m^2) GCaMax = 20.0 # Max. conductance of Calcium current (S/m^2) GKCaMax = 236.0 # Max. conductance of Calcium-dependent Potassium current (S/m^2) - GL = 1.0 # Conductance of non-specific leakage current (S/m^2) + GLeak = 1.0 # Conductance of non-specific leakage current (S/m^2) GPumpNa = 20.0 # Max. conductance of Sodium pump current (S/m^2) taum = 0.1e-3 # Sodium activation time constant (s) taus = 0.6e-3 # Calcium activation time constant (s) # Original conversion constants from inward ion current (nA) to build-up of # intracellular ion concentration (arb.) K_Na_original = 0.016 # iNa to intracellular [Na+] K_Ca_original = 0.1 # iCa to intracellular [Ca2+] # Constants needed to convert K from original model (soma compartment) # to current model (point-neuron) surface = 6434.0e-12 # surface of cell assumed as a single soma (m2) curr_factor = 1e6 # mA to nA # Time constants for the removal of ions from intracellular pools (s) tau_Na_removal = 16.0 # Na+ removal tau_Ca_removal = 1.25 # Ca2+ removal # Time constants for the iPumpNa and iKCa currents activation # from specific intracellular ions (s) tau_PumpNa_act = 0.1 # iPumpNa activation from intracellular Na+ tau_KCa_act = 0.01 # iKCa activation from intracellular Ca2+ # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h', 'm3h'], 'i_K\ kin.': ['n'], 'i_{Ca}\ kin.': ['s'], 'pools': ['C_Na_arb', 'C_Na_arb_activation', 'C_Ca_arb', 'C_Ca_arb_activation'], - 'I': ['iNa', 'iK', 'iCa', 'iKCa', 'iPumpNa', 'iL', 'iNet'] + 'I': ['iNa', 'iK', 'iCa', 'iKCa', 'iPumpNa', 'iLeak', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'C_Na', 'A_Na', 'C_Ca', 'A_Ca'] self.states0 = np.array([]) # Names of the channels effective coefficients self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas'] self.K_Na = self.K_Na_original * self.surface * self.curr_factor self.K_Ca = self.K_Ca_original * self.surface * self.curr_factor # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) # Charge interval bounds for lookup creation self.Qbounds = (np.round(self.Vm0 - 10.0) * 1e-5, 50.0e-5) # ----------------- Generic ----------------- def _xinf(self, Vm, halfmax, slope, power): ''' Generic function computing the steady-state activation/inactivation of a particular ion channel at a given voltage. :param Vm: membrane potential (mV) :param halfmax: half-(in)activation voltage (mV) :param slope: slope parameter of (in)activation function (mV) :param power: power exponent multiplying the exponential expression (integer) :return: steady-state (in)activation (-) ''' return 1 / (1 + np.exp((Vm - halfmax) / slope))**power def _taux(self, Vm, halfmax, slope, tauMax, tauMin): ''' Generic function computing the voltage-dependent, activation/inactivation time constant of a particular ion channel at a given voltage. :param Vm: membrane potential (mV) :param halfmax: voltage at which (in)activation time constant is half-maximal (mV) :param slope: slope parameter of (in)activation time constant function (mV) :return: steady-state (in)activation (-) ''' return (tauMax - tauMin) / (1 + np.exp((Vm - halfmax) / slope)) + tauMin def _derC_ion(self, Cion, Iion, Kion, tau): ''' Generic function computing the time derivative of the concentration of a specific ion in its intracellular pool. :param Cion: ion concentration in the pool (arbitrary unit) :param Iion: ionic current (mA/m2) :param Kion: scaling factor for current contribution to pool (arb. unit / nA???) :param tau: time constant for removal of ions from the pool (s) :return: variation of ionic concentration in the pool (arbitrary unit /s) ''' return (Kion * (-Iion) - Cion) / tau def _derA_ion(self, Aion, Cion, tau): ''' Generic function computing the time derivative of the concentration and time dependent activation function, for a specific pool-dependent ionic current. :param Aion: concentration and time dependent activation function (arbitrary unit) :param Cion: ion concentration in the pool (arbitrary unit) :param tau: time constant for activation function variation (s) :return: variation of activation function (arbitrary unit / s) ''' return (Cion - Aion) / tau # ------------------ Na ------------------- minf = partialmethod(_xinf, halfmax=-35.0, slope=-5.0, power=1) hinf = partialmethod(_xinf, halfmax=-50.0, slope=9.0, power=2) tauh = partialmethod(_taux, halfmax=-36.0, slope=3.5, tauMax=14.0e-3, tauMin=0.2e-3) def derM(self, Vm, m): ''' Instantaneous derivative of Sodium activation. ''' return (self.minf(Vm) - m) / self.taum # s-1 def derH(self, Vm, h): ''' Instantaneous derivative of Sodium inactivation. ''' return (self.hinf(Vm) - h) / self.tauh(Vm) # s-1 # ------------------ K ------------------- ninf = partialmethod(_xinf, halfmax=-22.0, slope=-9.0, power=1) taun = partialmethod(_taux, halfmax=-10.0, slope=10.0, tauMax=6.0e-3, tauMin=1.0e-3) def derN(self, Vm, n): ''' Instantaneous derivative of Potassium activation. ''' return (self.ninf(Vm) - n) / self.taun(Vm) # s-1 # ------------------ Ca ------------------- sinf = partialmethod(_xinf, halfmax=-10.0, slope=-2.8, power=1) def derS(self, Vm, s): ''' Instantaneous derivative of Calcium activation. ''' return (self.sinf(Vm) - s) / self.taus # s-1 # ------------------ Pools ------------------- def derC_Na(self, C_Na, I_Na): ''' Derivative of Sodium concentration in intracellular pool. ''' return self._derC_ion(C_Na, I_Na, self.K_Na, self.tau_Na_removal) def derA_Na(self, A_Na, C_Na): ''' Derivative of Sodium pool-dependent activation function for iPumpNa. ''' return self._derA_ion(A_Na, C_Na, self.tau_PumpNa_act) def derC_Ca(self, C_Ca, I_Ca): ''' Derivative of Calcium concentration in intracellular pool. ''' return self._derC_ion(C_Ca, I_Ca, self.K_Ca, self.tau_Ca_removal) def derA_Ca(self, A_Ca, C_Ca): ''' Derivative of Calcium pool-dependent activation function for iKCa. ''' return self._derA_ion(A_Ca, C_Ca, self.tau_KCa_act) # ------------------ Currents ------------------- def currNa(self, m, h, Vm): ''' Sodium inward current. ''' return self.GNaMax * m**3 * h * (Vm - self.VNa) def currK(self, n, Vm): ''' Potassium outward current. ''' return self.GKMax * n**2 * (Vm - self.VK) def currCa(self, s, Vm): ''' Calcium inward current. ''' return self.GCaMax * s * (Vm - self.VCa) def currKCa(self, A_Ca, Vm): ''' Calcium-activated Potassium outward current. ''' return self.GKCaMax * A_Ca * (Vm - self.VK) def currPumpNa(self, A_Na, Vm): ''' Outward current mimicking the activity of the NaK-ATPase pump. ''' return self.GPumpNa * A_Na * (Vm - self.VPumpNa) - def currL(self, Vm): + def currLeak(self, Vm): ''' Leakage current. ''' - return self.GL * (Vm - self.VL) + return self.GLeak * (Vm - self.VLeak) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, _, A_Na, _, A_Ca = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currCa(s, Vm) + - self.currL(Vm) + self.currPumpNa(A_Na, Vm) + self.currKCa(A_Ca, Vm)) # mA/m2 + self.currLeak(Vm) + self.currPumpNa(A_Na, Vm) + self.currKCa(A_Ca, Vm)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Standard gating dynamics: Solve the equation dx/dt = 0 at Vm for each x-state meq = self.minf(Vm) heq = self.hinf(Vm) neq = self.ninf(Vm) seq = self.sinf(Vm) # PumpNa pool concentration and activation steady-state INa_eq = self.currNa(meq, heq, Vm) CNa_eq = self.K_Na * (-INa_eq) ANa_eq = CNa_eq # KCa current pool concentration and activation steady-state ICa_eq = self.currCa(seq, Vm) CCa_eq = self.K_Ca * (-ICa_eq) ACa_eq = CCa_eq return np.array([meq, heq, neq, seq, CNa_eq, ANa_eq, CCa_eq, ACa_eq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack states m, h, n, s, C_Na, A_Na, C_Ca, A_Ca = states # Standard gating states derivatives dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) # PumpNa current pool concentration and activation state I_Na = self.currNa(m, h, Vm) dCNa_dt = self.derC_Na(C_Na, I_Na) dANa_dt = self.derA_Na(A_Na, C_Na) # KCa current pool concentration and activation state I_Ca = self.currCa(s, Vm) dCCa_dt = self.derC_Ca(C_Ca, I_Ca) dACa_dt = self.derA_Ca(A_Ca, C_Ca) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dCNa_dt, dANa_dt, dCCa_dt, dACa_dt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants Tm = self.taum minf = self.minf(Vm) am_avg = np.mean(minf / Tm) bm_avg = np.mean(1 / Tm) - am_avg Th = self.tauh(Vm) hinf = self.hinf(Vm) ah_avg = np.mean(hinf / Th) bh_avg = np.mean(1 / Th) - ah_avg Tn = self.taun(Vm) ninf = self.ninf(Vm) an_avg = np.mean(ninf / Tn) bn_avg = np.mean(1 / Tn) - an_avg Ts = self.taus sinf = self.sinf(Vm) as_avg = np.mean(sinf / Ts) bs_avg = np.mean(1 / Ts) - as_avg # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) Vmeff = np.interp(Qm, interp_data['Q'], interp_data['V']) # Unpack states m, h, n, s, C_Na, A_Na, C_Ca, A_Ca = states # Standard gating states derivatives dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s # PumpNa current pool concentration and activation state I_Na = self.currNa(m, h, Vmeff) dCNa_dt = self.derC_Na(C_Na, I_Na) dANa_dt = self.derA_Na(A_Na, C_Na) # KCa current pool concentration and activation state I_Ca_eff = self.currCa(s, Vmeff) dCCa_dt = self.derC_Ca(C_Ca, I_Ca_eff) dACa_dt = self.derA_Ca(A_Ca, C_Ca) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dCNa_dt, dANa_dt, dCCa_dt, dACa_dt] class LeechMech(PointNeuron): ''' Class defining the basic dynamics of Sodium, Potassium and Calcium channels for several neurons of the leech. Reference: *Baccus, S.A. (1998). Synaptic facilitation by reflected action potentials: enhancement of transmission when nerve impulses reverse direction at axon branch points. Proc. Natl. Acad. Sci. U.S.A. 95, 8345–8350.* ''' alphaC_sf = 1e-5 # Calcium activation rate constant scaling factor (M) betaC = 0.1e3 # beta rate for the open-probability of Ca2+-dependent Potassium channels (s-1) T = 293.15 # Room temperature (K) Rg = 8.314 # Universal gas constant (J.mol^-1.K^-1) Faraday = 9.6485e4 # Faraday constant (C/mol) def alpham(self, Vm): ''' Compute the alpha rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = -0.03 * (Vm + 28) / (np.exp(- (Vm + 28) / 15) - 1) # ms-1 return alpha * 1e3 # s-1 def betam(self, Vm): ''' Compute the beta rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 2.7 * np.exp(-(Vm + 53) / 18) # ms-1 return beta * 1e3 # s-1 def alphah(self, Vm): ''' Compute the alpha rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = 0.045 * np.exp(-(Vm + 58) / 18) # ms-1 return alpha * 1e3 # s-1 def betah(self, Vm): ''' Compute the beta rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) .. warning:: the original paper contains an error (multiplication) in the expression of this rate constant, corrected in the mod file on ModelDB (division). ''' beta = 0.72 / (np.exp(-(Vm + 23) / 14) + 1) # ms-1 return beta * 1e3 # s-1 def alphan(self, Vm): ''' Compute the alpha rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = -0.024 * (Vm - 17) / (np.exp(-(Vm - 17) / 8) - 1) # ms-1 return alpha * 1e3 # s-1 def betan(self, Vm): ''' Compute the beta rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 0.2 * np.exp(-(Vm + 48) / 35) # ms-1 return beta * 1e3 # s-1 def alphas(self, Vm): ''' Compute the alpha rate for the open-probability of Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = -1.5 * (Vm - 20) / (np.exp(-(Vm - 20) / 5) - 1) # ms-1 return alpha * 1e3 # s-1 def betas(self, Vm): ''' Compute the beta rate for the open-probability of Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 1.5 * np.exp(-(Vm + 25) / 10) # ms-1 return beta * 1e3 # s-1 def alphaC(self, C_Ca_in): ''' Compute the alpha rate for the open-probability of Calcium-dependent Potassium channels. :param C_Ca_in: intracellular Calcium concentration (M) :return: rate constant (s-1) ''' alpha = 0.1 * C_Ca_in / self.alphaC_sf # ms-1 return alpha * 1e3 # s-1 def derM(self, Vm, m): ''' Compute the evolution of the open-probability of Sodium channels. :param Vm: membrane potential (mV) :param m: open-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alpham(Vm) * (1 - m) - self.betam(Vm) * m def derH(self, Vm, h): ''' Compute the evolution of the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :param h: inactivation-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphah(Vm) * (1 - h) - self.betah(Vm) * h def derN(self, Vm, n): ''' Compute the evolution of the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :param n: open-probability of delayed-rectifier Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphan(Vm) * (1 - n) - self.betan(Vm) * n def derS(self, Vm, s): ''' Compute the evolution of the open-probability of Calcium channels. :param Vm: membrane potential (mV) :param s: open-probability of Calcium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphas(Vm) * (1 - s) - self.betas(Vm) * s def derC(self, c, C_Ca_in): ''' Compute the evolution of the open-probability of Calcium-dependent Potassium channels. :param c: open-probability of Calcium-dependent Potassium channels (prob) :param C_Ca_in: intracellular Calcium concentration (M) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphaC(C_Ca_in) * (1 - c) - self.betaC * c def currNa(self, m, h, Vm, C_Na_in): ''' Compute the inward Sodium current per unit area. :param m: open-probability of Sodium channels :param h: inactivation-probability of Sodium channels :param Vm: membrane potential (mV) :param C_Na_in: intracellular Sodium concentration (M) :return: current per unit area (mA/m2) ''' GNa = self.GNaMax * m**4 * h VNa = self.nernst(self.Z_Na, C_Na_in, self.C_Na_out) # Sodium Nernst potential return GNa * (Vm - VNa) def currK(self, n, Vm): ''' Compute the outward, delayed-rectifier Potassium current per unit area. :param n: open-probability of delayed-rectifier Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GKMax * n**2 return GK * (Vm - self.VK) def currCa(self, s, Vm, C_Ca_in): ''' Compute the inward Calcium current per unit area. :param s: open-probability of Calcium channels :param Vm: membrane potential (mV) :param C_Ca_in: intracellular Calcium concentration (M) :return: current per unit area (mA/m2) ''' GCa = self.GCaMax * s VCa = self.nernst(self.Z_Ca, C_Ca_in, self.C_Ca_out) # Calcium Nernst potential return GCa * (Vm - VCa) def currKCa(self, c, Vm): ''' Compute the outward Calcium-dependent Potassium current per unit area. :param c: open-probability of Calcium-dependent Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GKCa = self.GKCaMax * c return GKCa * (Vm - self.VK) - def currL(self, Vm): + def currLeak(self, Vm): ''' Compute the non-specific leakage current per unit area. :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' - return self.GL * (Vm - self.VL) + return self.GLeak * (Vm - self.VLeak) class LeechPressure(LeechMech): ''' Class defining the membrane channel dynamics of a leech pressure sensory neuron. with 7 different current types: - Inward Sodium current - Outward Potassium current - Inward high-voltage-activated Calcium current - Non-specific leakage current - Calcium-dependent, outward Potassium current - Sodium pump current - Calcium pump current Reference: *Baccus, S.A. (1998). Synaptic facilitation by reflected action potentials: enhancement of transmission when nerve impulses reverse direction at axon branch points. Proc. Natl. Acad. Sci. U.S.A. 95, 8345–8350.* ''' # Name of channel mechanism name = 'LeechP' # Cell-specific biophysical parameters Cm0 = 1e-2 # Cell membrane resting capacitance (F/m2) Vm0 = -48.865 # Cell membrane resting potential (mV) C_Na_out = 0.11 # Sodium extracellular concentration (M) C_Ca_out = 1.8e-3 # Calcium extracellular concentration (M) C_Na_in0 = 0.01 # Initial Sodium intracellular concentration (M) C_Ca_in0 = 1e-7 # Initial Calcium intracellular concentration (M) # VNa = 60 # Sodium Nernst potential, from MOD file on ModelDB (mV) # VCa = 125 # Calcium Nernst potential, from MOD file on ModelDB (mV) VK = -68.0 # Potassium Nernst potential (mV) - VL = -49.0 # Non-specific leakage Nernst potential (mV) + VLeak = -49.0 # Non-specific leakage Nernst potential (mV) INaPmax = 70.0 # Maximum pump rate of the NaK-ATPase (mA/m2) khalf_Na = 0.012 # Sodium concentration at which NaK-ATPase is at half its maximum rate (M) ksteep_Na = 1e-3 # Sensitivity of NaK-ATPase to varying Sodium concentrations (M) iCaS = 0.1 # Calcium pump current parameter (mA/m2) GNaMax = 3500.0 # Max. conductance of Sodium current (S/m^2) GKMax = 60.0 # Max. conductance of Potassium current (S/m^2) GCaMax = 0.02 # Max. conductance of Calcium current (S/m^2) GKCaMax = 8.0 # Max. conductance of Calcium-dependent Potassium current (S/m^2) - GL = 5.0 # Conductance of non-specific leakage current (S/m^2) + GLeak = 5.0 # Conductance of non-specific leakage current (S/m^2) diam = 50e-6 # Cell soma diameter (m) Z_Na = 1 # Sodium valence Z_Ca = 2 # Calcium valence # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h', 'm4h'], 'i_K\ kin.': ['n'], 'i_{Ca}\ kin.': ['s'], 'i_{KCa}\ kin.': ['c'], 'pools': ['C_Na', 'C_Ca'], - 'I': ['iNa2', 'iK', 'iCa2', 'iKCa2', 'iPumpNa2', 'iPumpCa2', 'iL', 'iNet'] + 'I': ['iNa2', 'iK', 'iCa2', 'iKCa2', 'iPumpNa2', 'iPumpCa2', 'iLeak', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' SV_ratio = 6 / self.diam # surface to volume ratio of the (spherical) cell soma # Conversion constant from membrane ionic currents into # change rate of intracellular ionic concentrations self.K_Na = SV_ratio / (self.Z_Na * self.Faraday) * 1e-6 # Sodium (M/s) self.K_Ca = SV_ratio / (self.Z_Ca * self.Faraday) * 1e-6 # Calcium (M/s) # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'c', 'C_Na', 'C_Ca'] self.states0 = np.array([]) # Names of the channels effective coefficients self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas'] # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) # Charge interval bounds for lookup creation self.Qbounds = (np.round(self.Vm0 - 10.0) * 1e-5, 60.0e-5) def nernst(self, z_ion, C_ion_in, C_ion_out): ''' Return the Nernst potential of a specific ion given its intra and extracellular concentrations. :param z_ion: ion valence :param C_ion_in: intracellular ion concentration (M) :param C_ion_out: extracellular ion concentration (M) :return: ion Nernst potential (mV) ''' return (self.Rg * self.T) / (z_ion * self.Faraday) * np.log(C_ion_out / C_ion_in) * 1e3 def currPumpNa(self, C_Na_in): ''' Outward current mimicking the activity of the NaK-ATPase pump. :param C_Na_in: intracellular Sodium concentration (M) :return: current per unit area (mA/m2) ''' return self.INaPmax / (1 + np.exp((self.khalf_Na - C_Na_in) / self.ksteep_Na)) def currPumpCa(self, C_Ca_in): ''' Outward current representing the activity of a Calcium pump. :param C_Ca_in: intracellular Calcium concentration (M) :return: current per unit area (mA/m2) ''' return self.iCaS * (C_Ca_in - self.C_Ca_in0) / 1.5 def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, c, C_Na_in, C_Ca_in = states return (self.currNa(m, h, Vm, C_Na_in) + self.currK(n, Vm) + self.currCa(s, Vm, C_Ca_in) + - self.currKCa(c, Vm) + self.currL(Vm) + + self.currKCa(c, Vm) + self.currLeak(Vm) + (self.currPumpNa(C_Na_in) / 3.) + self.currPumpCa(C_Ca_in)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Intracellular concentrations C_Na_eq = self.C_Na_in0 C_Ca_eq = self.C_Ca_in0 # Standard gating dynamics: Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) seq = self.alphas(Vm) / (self.alphas(Vm) + self.betas(Vm)) ceq = self.alphaC(C_Ca_eq) / (self.alphaC(C_Ca_eq) + self.betaC) return np.array([meq, heq, neq, seq, ceq, C_Na_eq, C_Ca_eq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack states m, h, n, s, c, C_Na_in, C_Ca_in = states # Standard gating states derivatives dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) dcdt = self.derC(c, C_Ca_in) # Intracellular concentrations dCNa_dt = - (self.currNa(m, h, Vm, C_Na_in) + self.currPumpNa(C_Na_in)) * self.K_Na # M/s dCCa_dt = -(self.currCa(s, Vm, C_Ca_in) + self.currPumpCa(C_Ca_in)) * self.K_Ca # M/s # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dCNa_dt, dCCa_dt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) as_avg = np.mean(self.alphas(Vm)) bs_avg = np.mean(self.betas(Vm)) # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) Vmeff = np.interp(Qm, interp_data['Q'], interp_data['V']) # Unpack states m, h, n, s, c, C_Na_in, C_Ca_in = states # Standard gating states derivatives dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s # KCa current gating state derivative dcdt = self.derC(c, C_Ca_in) # Intracellular concentrations dCNa_dt = - (self.currNa(m, h, Vmeff, C_Na_in) + self.currPumpNa(C_Na_in)) * self.K_Na # M/s dCCa_dt = -(self.currCa(s, Vmeff, C_Ca_in) + self.currPumpCa(C_Ca_in)) * self.K_Ca # M/s # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dCNa_dt, dCCa_dt] class LeechRetzius(LeechMech): ''' Class defining the membrane channel dynamics of a leech Retzius neuron. with 5 different current types: - Inward Sodium current - Outward Potassium current - Inward high-voltage-activated Calcium current - Non-specific leakage current - Calcium-dependent, outward Potassium current References: *Vazquez, Y., Mendez, B., Trueta, C., and De-Miguel, F.F. (2009). Summation of excitatory postsynaptic potentials in electrically-coupled neurones. Neuroscience 163, 202–212.* *ModelDB link: https://senselab.med.yale.edu/modeldb/ShowModel.cshtml?model=120910* ''' # Name of channel mechanism # name = 'LeechR' # Cell-specific biophysical parameters Cm0 = 5e-2 # Cell membrane resting capacitance (F/m2) Vm0 = -44.45 # Cell membrane resting potential (mV) VNa = 50.0 # Sodium Nernst potential, from retztemp.ses file on ModelDB (mV) VCa = 125.0 # Calcium Nernst potential, from cachdend.mod file on ModelDB (mV) VK = -79.0 # Potassium Nernst potential, from retztemp.ses file on ModelDB (mV) - VL = -30.0 # Non-specific leakage Nernst potential, from leakdend.mod file on ModelDB (mV) + VLeak = -30.0 # Non-specific leakage Nernst potential, from leakdend.mod file on ModelDB (mV) GNaMax = 1250.0 # Max. conductance of Sodium current (S/m^2) GKMax = 10.0 # Max. conductance of Potassium current (S/m^2) GAMax = 100.0 # Max. conductance of transient Potassium current (S/m^2) GCaMax = 4.0 # Max. conductance of Calcium current (S/m^2) GKCaMax = 130.0 # Max. conductance of Calcium-dependent Potassium current (S/m^2) - GL = 1.25 # Conductance of non-specific leakage current (S/m^2) + GLeak = 1.25 # Conductance of non-specific leakage current (S/m^2) Vhalf = -73.1 # mV C_Ca_in = 5e-8 # Calcium intracellular concentration, from retztemp.ses file (M) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h', 'm4h'], 'i_K\ kin.': ['n'], 'i_A\ kin.': ['a', 'b', 'ab'], 'i_{Ca}\ kin.': ['s'], 'i_{KCa}\ kin.': ['c'], - 'I': ['iNa', 'iK', 'iCa', 'iKCa2', 'iA', 'iL', 'iNet'] + 'I': ['iNa', 'iK', 'iCa', 'iKCa2', 'iA', 'iLeak', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'c', 'a', 'b'] self.states0 = np.array([]) # Names of the channels effective coefficients self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas', 'alphac', 'betac', 'alphaa', 'betaa' 'alphab', 'betab'] # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) def ainf(self, Vm): ''' Steady-state activation probability of transient Potassium channels. Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' Vth = -55.0 # mV return 0 if Vm <= Vth else min(1, 2 * (Vm - Vth)**3 / ((11 - Vth)**3 + (Vm - Vth)**3)) def taua(self, Vm): ''' Activation time constant of transient Potassium channels. (assuming T = 20°C). Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' x = -1.5 * (Vm - self.Vhalf) * 1e-3 * self.Faraday / (self.Rg * self.T) # [-] alpha = np.exp(x) # ms-1 beta = np.exp(0.7 * x) # ms-1 return max(0.5, beta / (0.3 * (1 + alpha))) * 1e-3 # s def binf(self, Vm): ''' Steady-state inactivation probability of transient Potassium channels. Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' return 1. / (1 + np.exp((self.Vhalf - Vm) / -6.3)) def taub(self, Vm): ''' Inactivation time constant of transient Potassium channels. (assuming T = 20°C). Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' x = 2 * (Vm - self.Vhalf) * 1e-3 * self.Faraday / (self.Rg * self.T) alpha = np.exp(x) beta = np.exp(0.65 * x) return max(7.5, beta / (0.02 * (1 + alpha))) * 1e-3 # s def derA(self, Vm, a): ''' Compute the evolution of the activation-probability of transient Potassium channels. :param Vm: membrane potential (mV) :param a: activation-probability of transient Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.ainf(Vm) - a) / self.taua(Vm) def derB(self, Vm, b): ''' Compute the evolution of the inactivation-probability of transient Potassium channels. :param Vm: membrane potential (mV) :param b: inactivation-probability of transient Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.binf(Vm) - b) / self.taub(Vm) def currA(self, a, b, Vm): ''' Compute the outward, transient Potassium current per unit area. :param a: open-probability of transient Potassium channels :param b: inactivation-probability of transient Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GAMax * a * b return GK * (Vm - self.VK) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, c, a, b = states - return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currCa(s, Vm) + self.currL(Vm) + + return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currCa(s, Vm) + self.currLeak(Vm) + self.currKCa(c, Vm) + self.currA(a, b, Vm)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Standard gating dynamics: Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) seq = self.alphas(Vm) / (self.alphas(Vm) + self.betas(Vm)) ceq = self.alphaC(self.C_Ca_in) / (self.alphaC(self.C_Ca_in) + self.betaC) aeq = self.ainf(Vm) beq = self.binf(Vm) return np.array([meq, heq, neq, seq, ceq, aeq, beq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack states m, h, n, s, c, a, b = states # Standard gating states derivatives dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) dcdt = self.derC(c, self.C_Ca_in) dadt = self.derA(Vm, a) dbdt = self.derB(Vm, b) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dadt, dbdt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) as_avg = np.mean(self.alphas(Vm)) bs_avg = np.mean(self.betas(Vm)) Ta = self.taua(Vm) ainf = self.ainf(Vm) aa_avg = np.mean(ainf / Ta) ba_avg = np.mean(1 / Ta) - aa_avg Tb = self.taub(Vm) binf = self.binf(Vm) ab_avg = np.mean(binf / Tb) bb_avg = np.mean(1 / Tb) - ab_avg # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg, aa_avg, ba_avg, ab_avg, bb_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) # Unpack states m, h, n, s, c, a, b = states # Standard gating states derivatives dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s dadt = rates[8] * (1 - a) - rates[9] * a dbdt = rates[10] * (1 - b) - rates[11] * b # KCa current gating state derivative dcdt = self.derC(c, self.C_Ca_in) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dadt, dbdt] diff --git a/PySONIC/neurons/stn.py b/PySONIC/neurons/stn.py new file mode 100644 index 0000000..85f706e --- /dev/null +++ b/PySONIC/neurons/stn.py @@ -0,0 +1,556 @@ + +import numpy as np +from ..core import PointNeuron +from ..utils import si_format + + +class OtsukaSTN(PointNeuron): + ''' Class defining the Otsuka model of sub-thalamic nucleus neuron + with 5 different current types: + - Inward Sodium current (iNa) + - Outward, delayed-rectifer Potassium current (iK) + - Inward, A-type Potassium current (iA) + - Inward, low-threshold Calcium current (iT) + - Inward, high-threshold Calcium current (iL) + - Outward, Calcium-dependent Potassium current (iCaK) + - Non-specific leakage current (iLeak) + + Reference: + *Tarnaud, T., Joseph, W., Martens, L., and Tanghe, E. (2018). Computational Modeling + of Ultrasonic Subthalamic Nucleus Stimulation. IEEE Trans Biomed Eng.* + ''' + + name = 'STN' + + # Generic biophysical parameters of STN cells + Cm0 = 1e-2 # Cell membrane resting capacitance (F/m2) + Vm0 = -58.0 # Resting membrane potential (mV) + + VNa = 60.0 # Sodium Nernst potential (mV) + VK = -90.0 # Potassium Nernst potential (mV) + VLeak = -60.0 # Leakage reversal potential (mV) + GNaMax = 490.0 # Max. conductance of Sodium current (S/m^2) + GKMax = 570.0 # Max. conductance of delayed-rectifier Potassium current (S/m^2) + GAMax = 50.0 # Max. conductance of A-type Potassium current (S/m^2) + GTMax = 50.0 # Max. conductance of low-threshold Calcium current (S/m^2) + GLMax = 150.0 # Max. conductance of high-threshold Calcium current (S/m^2) + GCaKMax = 10.0 # Max. conductance of Calcium-dependent Potassium current (S/m^2) + GLeak = 3.5 # Conductance of non-specific leakage current (S/m^2) + + Faraday = 9.64853415e4 # Faraday constant for (Coulomb / mole) + Rg = 8.314 # Universal gas constant (Pa.m^3.mol^-1.K^-1) + T = 306.15 # K (33°C) + + nCa = 2 # Calcium ion valence + CCa_out = 2e3 # uM (2 mM) + KCa = 2e3 # s-1 + CCa_in0 = 5e-3 # uM + + # Steady-state and time constants parameters + thetax_a = -45 # mV + thetax_b = -90 # mV + thetax_c = -30.6 # mV + thetax_d1 = -60 # mV + thetax_d2 = 0.1 # uM + thetax_m = -40 # mV + thetax_h = -45.5 # mV + thetax_n = -41 # mV + thetax_p = -56 # mV + thetax_q = -85 # mV + thetax_r = 0.17 # uM + + kx_a = -14.7 # mV + kx_b = 7.5 # mV + kx_c = -5 # mV + kx_d1 = 7.5 # mV + kx_d2 = 0.02 # uM + kx_m = -8 # mV + kx_h = 6.4 # mV + kx_n = -14 # mV + kx_p = -6.7 # mV + kx_q = 5.8 # mV + kx_r = -0.08 # uM + + tau0_a = 1 * 1e-3 # s + tau1_a = 1 * 1e-3 # s + tau0_b = 0 * 1e-3 # s + tau1_b = 200 * 1e-3 # s + tau0_c = 45 * 1e-3 # s + tau1_c = 10 * 1e-3 # s + tau0_d1 = 400 * 1e-3 # s + tau1_d1 = 500 * 1e-3 # s + tau_d2 = 130 * 1e-3 # s + tau0_m = 0.2 * 1e-3 # s + tau1_m = 3 * 1e-3 # s + tau0_h = 0 * 1e-3 # s + tau1_h = 24.5 * 1e-3 # s + tau0_n = 0 * 1e-3 # s + tau1_n = 11 * 1e-3 # s + tau0_p = 5 * 1e-3 # s + tau1_p = 0.33 * 1e-3 # s + tau0_q = 0 * 1e-3 # s + tau1_q = 400 * 1e-3 # s + tau0_r = 2 * 1e-3 # s + tau_r = 2 * 1e-3 # s + + thetaT_a = -40 # mV + thetaT1_b = -60 # mV + thetaT2_b = -40 # mV + thetaT1_c = -27 # mV + thetaT2_c = -50 # mV + thetaT1_d1 = -40 # mV + thetaT2_d1 = -20 # mV + thetaT_m = -53 # mV + thetaT1_h = -50 # mV + thetaT2_h = -50 # mV + thetaT1_n = -40 # mV + thetaT2_n = -40 # mV + thetaT1_p = -27 # mV + thetaT2_p = -102 # mV + thetaT1_q = -50 # mV + thetaT2_q = -50 # mV + + sigmaT_a = -0.5 # mV + sigmaT1_b = -30 # mV + sigmaT2_b = -10 # mV + sigmaT1_c = -20 # mV + sigmaT2_c = 15 # mV + sigmaT1_d1 = -15 # mV + sigmaT2_d1 = 20 # mV + sigmaT_m = -0.7 # mV + sigmaT1_h = -15 # mV + sigmaT2_h = 16 # mV + sigmaT1_n = -40 # mV + sigmaT2_n = 50 # mV + sigmaT1_p = -10 # mV + sigmaT2_p = 15 # mV + sigmaT1_q = -15 # mV + sigmaT2_q = 16 # mV + + # Default plotting scheme + pltvars_scheme = { + 'i_{Na}\ kin.': ['m', 'h', 'm3h'], + 'i_K\ kin.': ['n', 'n4'], + 'i_A\ kin.': ['a', 'b', 'a2b'], + 'i_T\ kin.': ['p', 'q', 'p2q'], + 'i_L\ kin.': ['c', 'd1', 'd2', 'c2d1d2'], + 'i_{CaK}\ kin.': ['r', 'r2'], + 'I': ['iNa', 'iK', 'iA', 'iT2', 'iL', 'iCaK', 'iLeak', 'iNet'] + } + + + def __init__(self): + ''' Constructor of the class ''' + + # Names and initial states of the channels state probabilities + self.states_names = ['a', 'b', 'c', 'd1', 'd2', 'm', 'h', 'n', 'p', 'q', 'r', 'C_Ca'] + + # Names of the different coefficients to be averaged in a lookup table. + self.coeff_names = [ + 'alphaa', 'betaa', + 'alphab', 'betab', + 'alphac', 'betac', + 'alphad1', 'betad1', + 'alpham', 'betam', + 'alphah', 'betah', + 'alphan', 'betan', + 'alphap', 'betap', + 'alphaq', 'betaq', + ] + + # Compute deff for Cai = 5 nM + self.VCa = self.nernst(self.CCa_out, self.CCa_in0) + iT = self.currT(self.pinf(self.Vm0), self.qinf(self.Vm0), self.Vm0) + iL = self.currL(self.cinf(self.Vm0), self.d1inf(self.Vm0), self.d2inf(self.CCa_in0), self.Vm0) + self.deff = -(iT + iL) / (self.nCa * self.Faraday * self.KCa * self.CCa_in0) # m + print('deff = {}m'.format(si_format(self.deff, 2))) + + # Compute conversion factor from electrical current (mA/m2) to Calcium concentration (uM) + self.i2CCa = 1 / (self.nCa * self.deff * self.Faraday) + + # Initial states + self.states0 = self.steadyStates(self.Vm0) + + # Charge interval bounds for lookup creation + self.Qbounds = (np.round(self.Vm0 - 25.0) * 1e-5, 50.0e-5) + + + def nernst(self, xout, xin): + return self.Rg * self.T / (2 * self.Faraday) * np.log(xout / xin) + + + def _xinf(self, var, theta, k): + ''' Generic function computing the steady-state activation/inactivation of a + particular ion channel at a given voltage or ion concentration. + + :param var: membrane potential (mV) or ion concentration (mM) + :param theta: half-(in)activation voltage or concentration (mV or mM) + :param k: slope parameter of (in)activation function (mV or mM) + :return: steady-state (in)activation (-) + ''' + return 1 / (1 + np.exp((var - theta) / k)) + + def ainf(self, Vm): + return self._xinf(Vm, self.thetax_a, self.kx_a) + + def binf(self, Vm): + return self._xinf(Vm, self.thetax_b, self.kx_b) + + def cinf(self, Vm): + return self._xinf(Vm, self.thetax_c, self.kx_c) + + def d1inf(self, Vm): + return self._xinf(Vm, self.thetax_d1, self.kx_d1) + + def d2inf(self, Cai): + return self._xinf(Cai, self.thetax_d2, self.kx_d2) + + def minf(self, Vm): + return self._xinf(Vm, self.thetax_m, self.kx_m) + + def hinf(self, Vm): + return self._xinf(Vm, self.thetax_h, self.kx_h) + + def ninf(self, Vm): + return self._xinf(Vm, self.thetax_n, self.kx_n) + + def pinf(self, Vm): + return self._xinf(Vm, self.thetax_p, self.kx_p) + + def qinf(self, Vm): + return self._xinf(Vm, self.thetax_q, self.kx_q) + + def rinf(self, Cai): + return self._xinf(Cai, self.thetax_r, self.kx_r) + + + def _taux1(self, Vm, theta, sigma, tau0, tau1): + ''' Generic function computing the voltage-dependent, activation/inactivation time constant + of a particular ion channel at a given voltage (first variant). + + :param Vm: membrane potential (mV) + :param theta: voltage at which (in)activation time constant is half-maximal (mV) + :param sigma: slope parameter of (in)activation time constant function (mV) + :param tau0: minimal time constant (s) + :param tau1: modulated time constant (s) + :return: (in)activation time constant (s) + ''' + return tau0 + tau1 / (1 + np.exp(-(Vm - theta) / sigma)) + + def taua(self, Vm): + return self._taux1(Vm, self.thetaT_a, self.sigmaT_a, self.tau0_a, self.tau1_a) + + def taum(self, Vm): + return self._taux1(Vm, self.thetaT_m, self.sigmaT_m, self.tau0_m, self.tau1_m) + + + def _taux2(self, Vm, theta1, theta2, sigma1, sigma2, tau1, tau0): + ''' Generic function computing the voltage-dependent, activation/inactivation time constant + of a particular ion channel at a given voltage (second variant). + + :param Vm: membrane potential (mV) + :param theta: voltage at which (in)activation time constant is half-maximal (mV) + :param sigma: slope parameter of (in)activation time constant function (mV) + :param tau0: minimal time constant (s) + :param tau1: modulated time constant (s) + :return: (in)activation time constant (s) + ''' + return tau0 + tau1 / (np.exp(-(Vm - theta1) / sigma1) + np.exp(-(Vm - theta2) / sigma2)) + + def taub(self, Vm): + return self._taux2(Vm, self.thetaT1_b, self.thetaT2_b, self.sigmaT1_b, self.sigmaT2_b, + self.tau0_b, self.tau1_b) + + def tauc(self, Vm): + return self._taux2(Vm, self.thetaT1_c, self.thetaT2_c, self.sigmaT1_c, self.sigmaT2_c, + self.tau0_c, self.tau1_c) + + def taud1(self, Vm): + return self._taux2(Vm, self.thetaT1_d1, self.thetaT2_d1, self.sigmaT1_d1, self.sigmaT2_d1, + self.tau0_d1, self.tau1_d1) + + def tauh(self, Vm): + return self._taux2(Vm, self.thetaT1_h, self.thetaT2_h, self.sigmaT1_h, self.sigmaT2_h, + self.tau0_h, self.tau1_h) + + def taun(self, Vm): + return self._taux2(Vm, self.thetaT1_n, self.thetaT2_n, self.sigmaT1_n, self.sigmaT2_n, + self.tau0_n, self.tau1_n) + + def taup(self, Vm): + return self._taux2(Vm, self.thetaT1_p, self.thetaT2_p, self.sigmaT1_p, self.sigmaT2_p, + self.tau0_p, self.tau1_p) + + def tauq(self, Vm): + return self._taux2(Vm, self.thetaT1_q, self.thetaT2_q, self.sigmaT1_q, self.sigmaT2_q, + self.tau0_q, self.tau1_q) + + + def derA(self, Vm, a): + return (self.ainf(Vm) - a) / self.taua(Vm) + + def derB(self, Vm, b): + return (self.binf(Vm) - b) / self.taub(Vm) + + def derC(self, Vm, c): + return (self.cinf(Vm) - c) / self.tauc(Vm) + + def derD1(self, Vm, d1): + return (self.d1inf(Vm) - d1) / self.taud1(Vm) + + def derD2(self, Cai, d2): + return (self.d2inf(Cai) - d2) / self.tau_d2 + + def derM(self, Vm, m): + return (self.minf(Vm) - m) / self.taum(Vm) + + def derH(self, Vm, h): + return (self.hinf(Vm) - h) / self.tauh(Vm) + + def derN(self, Vm, n): + return (self.ninf(Vm) - n) / self.taun(Vm) + + def derP(self, Vm, p): + return (self.pinf(Vm) - p) / self.taup(Vm) + + def derQ(self, Vm, q): + return (self.qinf(Vm) - q) / self.tauq(Vm) + + def derR(self, Cai, r): + return (self.rinf(Cai) - r) / self.tau_r + + + def derC_Ca(self, C_Ca, iT, iL): + ''' Compute the evolution of the Calcium concentration in submembranal space. + + :param Vm: membrane potential (mV) + :param C_Ca: Calcium concentration in submembranal space (uM) + :param iT: inward, low-threshold Calcium current (mA/m2) + :param iL: inward, high-threshold Calcium current (mA/m2) + :return: derivative of Calcium concentration in submembranal space w.r.t. time (uM/s) + ''' + return - self.i2CCa * (iT + iL) - C_Ca * self.KCa + + + def currNa(self, m, h, Vm): + ''' Compute the inward Sodium current per unit area. + + :param m: open-probability of m-gate + :param h: inactivation-probability of h-gate + :param Vm: membrane potential (mV) + :return: current per unit area (mA/m2) + ''' + return self.GNaMax * m**3 * h * (Vm - self.VNa) + + + def currK(self, n, Vm): + ''' Compute the outward delayed-rectifier Potassium current per unit area. + + :param n: open-probability of n-gate + :param Vm: membrane potential (mV) + :return: current per unit area (mA/m2) + ''' + return self.GKMax * n**4 * (Vm - self.VK) + + + def currA(self, a, b, Vm): + ''' Compute the outward A-type Potassium current per unit area. + + :param a: open-probability of a-gate + :param b: open-probability of b-gate + :param Vm: membrane potential (mV) + :return: current per unit area (mA/m2) + ''' + return self.GAMax * a**2 * b * (Vm - self.VK) + + + def currT(self, p, q, Vm): + ''' Compute the inward low-threshold Calcium current per unit area. + + :param p: open-probability of p-gate + :param q: open-probability of q-gate + :param Vm: membrane potential (mV) + :return: current per unit area (mA/m2) + ''' + return self.GTMax * p**2 * q * (Vm - self.VCa) + + + def currL(self, c, d1, d2, Vm): + ''' Compute the inward high-threshold Calcium current per unit area. + + :param c: open-probability of c-gate + :param d1: open-probability of d1-gate + :param d2: open-probability of d2-gate + :param Vm: membrane potential (mV) + :return: current per unit area (mA/m2) + ''' + return self.GLMax * c**2 * d1 * d2 * (Vm - self.VCa) + + + def currCaK(self, r, Vm): + ''' Compute the outward, Calcium activated Potassium current per unit area. + + :param r: open-probability of r-gate + :param Vm: membrane potential (mV) + :return: current per unit area (mA/m2) + ''' + return self.GCaKMax * r**2 * (Vm - self.VK) + + + def currLeak(self, Vm): + ''' Compute the non-specific leakage current per unit area. + + :param Vm: membrane potential (mV) + :return: current per unit area (mA/m2) + ''' + return self.GLeak * (Vm - self.VLeak) + + + def currNet(self, Vm, states): + ''' Compute net membrane current per unit area. ''' + + a, b, c, d1, d2, m, h, n, p, q, r, CCa_in = states + + return ( + self.currNa(m, h, Vm) + + self.currK(n, Vm) + + self.currA(a, b, Vm) + + self.currT(p, q, Vm) + + self.currL(c, d1, d2, Vm) + + self.currCaK(r, Vm) + + self.currLeak(Vm) + ) # mA/m2 + + + def steadyStates(self, Vm): + ''' Concrete implementation of the abstract API method. ''' + + # Solve the equation dx/dt = 0 at Vm for each x-state + aeq = self.ainf(Vm) + beq = self.binf(Vm) + ceq = self.cinf(Vm) + d1eq = self.d1inf(Vm) + meq = self.minf(Vm) + heq = self.hinf(Vm) + neq = self.ninf(Vm) + peq = self.pinf(Vm) + qeq = self.qinf(Vm) + + d2eq = self.d2inf(self.CCa_in0) + req = self.rinf(self.CCa_in0) + + return np.array([aeq, beq, ceq, d1eq, d2eq, meq, heq, neq, peq, qeq, req, self.CCa_in0]) + + + def derStates(self, Vm, states): + ''' Concrete implementation of the abstract API method. ''' + + a, b, c, d1, d2, m, h, n, p, q, r, CCa_in = states + dadt = self.derA(Vm, a) + dbdt = self.derB(Vm, b) + dcdt = self.derC(Vm, c) + dd1dt = self.derD1(Vm, d1) + dd2dt = self.derD2(CCa_in, d2) + dmdt = self.derM(Vm, m) + dhdt = self.derH(Vm, h) + dndt = self.derN(Vm, n) + dpdt = self.derP(Vm, p) + dqdt = self.derQ(Vm, q) + drdt = self.derR(CCa_in, r) + + iT = self.currT(p, q, Vm) + iL = self.currL(c, d1, d2, Vm) + dCCaindt = self.derC_Ca(CCa_in, iT, iL) + + return [dadt, dbdt, dcdt, dd1dt, dd2dt, dmdt, dhdt, dndt, dpdt, dqdt, drdt, dCCaindt] + + + def getEffRates(self, Vm): + ''' Concrete implementation of the abstract API method. ''' + + # Compute average cycle value for rate constants + Ta = self.taua(Vm) + alphaa_avg = np.mean(self.ainf(Vm) / Ta) + betaa_avg = np.mean(1 / Ta) - alphaa_avg + + Tb = self.taub(Vm) + alphab_avg = np.mean(self.binf(Vm) / Tb) + betab_avg = np.mean(1 / Tb) - alphab_avg + + Tc = self.tauc(Vm) + alphac_avg = np.mean(self.cinf(Vm) / Tc) + betac_avg = np.mean(1 / Tc) - alphac_avg + + Td1 = self.taud1(Vm) + alphad1_avg = np.mean(self.ainf(Vm) / Td1) + betad1_avg = np.mean(1 / Td1) - alphad1_avg + + Td2 = self.tau_d2 + alphad2_avg = np.mean(self.d2inf(Vm)) / Td2 + betad2_avg = 1 / Td2 - alphad2_avg + + Tm = self.taum(Vm) + alpham_avg = np.mean(self.minf(Vm) / Tm) + betam_avg = np.mean(1 / Tm) - alpham_avg + + Th = self.tauh(Vm) + alphah_avg = np.mean(self.hinf(Vm) / Th) + betah_avg = np.mean(1 / Th) - alphah_avg + + Tn = self.taun(Vm) + alphan_avg = np.mean(self.ninf(Vm) / Tn) + betan_avg = np.mean(1 / Tn) - alphan_avg + + Tp = self.taup(Vm) + alphap_avg = np.mean(self.pinf(Vm) / Tp) + betap_avg = np.mean(1 / Tp) - alphap_avg + + Tq = self.tauq(Vm) + alphaq_avg = np.mean(self.qinf(Vm) / Tq) + betaq_avg = np.mean(1 / Tq) - alphaq_avg + + Tr = self.tau_r + alphar_avg = np.mean(self.rinf(Vm)) / Tr + betar_avg = 1 / Tr - alphar_avg + + # Return array of coefficients + return np.array([ + alphaa_avg, betaa_avg, + alphab_avg, betab_avg, + alphac_avg, betac_avg, + alphad1_avg, betad1_avg, + alphad2_avg, betad2_avg, + alpham_avg, betam_avg, + alphah_avg, betah_avg, + alphan_avg, betan_avg, + alphap_avg, betap_avg, + alphaq_avg, betaq_avg, + alphar_avg, betar_avg, + ]) + + + def derStatesEff(self, Qm, states, interp_data): + ''' Concrete implementation of the abstract API method. ''' + + rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) + for rn in self.coeff_names]) + Vmeff = np.interp(Qm, interp_data['Q'], interp_data['V']) + + a, b, c, d1, d2, m, h, n, p, q, r, CCa_in = states + dadt = rates[0] - (1 - a) - rates[1] * a + dbdt = rates[2] - (1 - b) - rates[3] * b + dcdt = rates[4] - (1 - c) - rates[5] * c + dd1dt = rates[6] - (1 - d1) - rates[7] * d1 + dd2dt = rates[8] - (1 - d2) - rates[9] * d2 + dmdt = rates[10] - (1 - m) - rates[11] * m + dhdt = rates[12] - (1 - h) - rates[13] * h + dndt = rates[14] - (1 - n) - rates[15] * n + dpdt = rates[16] - (1 - p) - rates[17] * p + dqdt = rates[18] - (1 - q) - rates[19] * q + drdt = rates[20] - (1 - r) - rates[21] * r + + iT = self.currT(p, q, Vmeff) + iL = self.currL(self.c, d1, d2, Vmeff) + dCCaindt = self.derC_Ca(CCa_in, iT, iL) + + return [dadt, dbdt, dcdt, dd1dt, dd2dt, dmdt, dhdt, dndt, dpdt, dqdt, drdt, dCCaindt] diff --git a/PySONIC/neurons/thalamic.py b/PySONIC/neurons/thalamic.py index 37a653d..0a2f649 100644 --- a/PySONIC/neurons/thalamic.py +++ b/PySONIC/neurons/thalamic.py @@ -1,786 +1,786 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-07-31 15:20:54 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-09-28 14:06:20 +# @Last Modified time: 2018-11-30 10:34:03 import numpy as np from ..core import PointNeuron class Thalamic(PointNeuron): ''' Class defining the generic membrane channel dynamics of a thalamic neuron with 4 different current types: - Inward Sodium current - Outward Potassium current - Inward Calcium current - Non-specific leakage current This generic class cannot be used directly as it does not contain any specific parameters. Reference: *Plaksin, M., Kimmel, E., and Shoham, S. (2016). Cell-Type-Selective Effects of Intramembrane Cavitation as a Unifying Theoretical Framework for Ultrasonic Neuromodulation. eNeuro 3.* ''' # Generic biophysical parameters of thalamic cells Cm0 = 1e-2 # Cell membrane resting capacitance (F/m2) Vm0 = 0.0 # Dummy value for membrane potential (mV) VNa = 50.0 # Sodium Nernst potential (mV) VK = -90.0 # Potassium Nernst potential (mV) VCa = 120.0 # Calcium Nernst potential (mV) def __init__(self): ''' Constructor of the class ''' # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'u'] self.states0 = np.array([]) # Names of the different coefficients to be averaged in a lookup table. self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas', 'alphau', 'betau'] # Charge interval bounds for lookup creation self.Qbounds = (np.round(self.Vm0 - 25.0) * 1e-5, 50.0e-5) def alpham(self, Vm): ''' Compute the alpha rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = 0.32 * self.vtrap(13 - Vdiff, 4) # ms-1 return alpha * 1e3 # s-1 def betam(self, Vm): ''' Compute the beta rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = 0.28 * self.vtrap(Vdiff - 40, 5) # ms-1 return beta * 1e3 # s-1 def alphah(self, Vm): ''' Compute the alpha rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = (0.128 * np.exp(-(Vdiff - 17) / 18)) # ms-1 return alpha * 1e3 # s-1 def betah(self, Vm): ''' Compute the beta rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = (4 / (1 + np.exp(-(Vdiff - 40) / 5))) # ms-1 return beta * 1e3 # s-1 def alphan(self, Vm): ''' Compute the alpha rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = 0.032 * self.vtrap(15 - Vdiff, 5) # ms-1 return alpha * 1e3 # s-1 def betan(self, Vm): ''' Compute the beta rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = (0.5 * np.exp(-(Vdiff - 10) / 40)) # ms-1 return beta * 1e3 # s-1 def derM(self, Vm, m): ''' Compute the evolution of the open-probability of Sodium channels. :param Vm: membrane potential (mV) :param m: open-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alpham(Vm) * (1 - m) - self.betam(Vm) * m def derH(self, Vm, h): ''' Compute the evolution of the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :param h: inactivation-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphah(Vm) * (1 - h) - self.betah(Vm) * h def derN(self, Vm, n): ''' Compute the evolution of the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :param n: open-probability of delayed-rectifier Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphan(Vm) * (1 - n) - self.betan(Vm) * n def derS(self, Vm, s): ''' Compute the evolution of the open-probability of the S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :param s: open-probability of S-type Calcium activation gates (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.sinf(Vm) - s) / self.taus(Vm) def derU(self, Vm, u): ''' Compute the evolution of the open-probability of the U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :param u: open-probability of U-type Calcium inactivation gates (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.uinf(Vm) - u) / self.tauu(Vm) def currNa(self, m, h, Vm): ''' Compute the inward Sodium current per unit area. :param m: open-probability of Sodium channels :param h: inactivation-probability of Sodium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GNa = self.GNaMax * m**3 * h return GNa * (Vm - self.VNa) def currK(self, n, Vm): ''' Compute the outward delayed-rectifier Potassium current per unit area. :param n: open-probability of delayed-rectifier Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GKMax * n**4 return GK * (Vm - self.VK) def currCa(self, s, u, Vm): ''' Compute the inward Calcium current per unit area. :param s: open-probability of the S-type activation gate of Calcium channels :param u: open-probability of the U-type inactivation gate of Calcium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GT = self.GTMax * s**2 * u return GT * (Vm - self.VCa) - def currL(self, Vm): + def currLeak(self, Vm): ''' Compute the non-specific leakage current per unit area. :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' - return self.GL * (Vm - self.VL) + return self.GLeak * (Vm - self.VLeak) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, u = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + - self.currCa(s, u, Vm) + self.currL(Vm)) # mA/m2 + self.currCa(s, u, Vm) + self.currLeak(Vm)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) seq = self.sinf(Vm) ueq = self.uinf(Vm) return np.array([meq, heq, neq, seq, ueq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, u = states dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) dudt = self.derU(Vm, u) return [dmdt, dhdt, dndt, dsdt, dudt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) Ts = self.taus(Vm) as_avg = np.mean(self.sinf(Vm) / Ts) bs_avg = np.mean(1 / Ts) - as_avg Tu = np.array([self.tauu(v) for v in Vm]) au_avg = np.mean(self.uinf(Vm) / Tu) bu_avg = np.mean(1 / Tu) - au_avg # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg, au_avg, bu_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) m, h, n, s, u = states dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s dudt = rates[8] * (1 - u) - rates[9] * u return [dmdt, dhdt, dndt, dsdt, dudt] class ThalamicRE(Thalamic): ''' Specific membrane channel dynamics of a thalamic reticular neuron. References: *Destexhe, A., Contreras, D., Steriade, M., Sejnowski, T.J., and Huguenard, J.R. (1996). In vivo, in vitro, and computational analysis of dendritic calcium currents in thalamic reticular neurons. J. Neurosci. 16, 169–185.* *Huguenard, J.R., and Prince, D.A. (1992). A novel T-type current underlies prolonged Ca(2+)-dependent burst firing in GABAergic neurons of rat thalamic reticular nucleus. J. Neurosci. 12, 3804–3817.* ''' # Name of channel mechanism name = 'RE' # Cell-specific biophysical parameters Vm0 = -89.5 # Cell membrane resting potential (mV) GNaMax = 2000.0 # Max. conductance of Sodium current (S/m^2) GKMax = 200.0 # Max. conductance of Potassium current (S/m^2) GTMax = 30.0 # Max. conductance of low-threshold Calcium current (S/m^2) - GL = 0.5 # Conductance of non-specific leakage current (S/m^2) - VL = -90.0 # Non-specific leakage Nernst potential (mV) + GLeak = 0.5 # Conductance of non-specific leakage current (S/m^2) + VLeak = -90.0 # Non-specific leakage Nernst potential (mV) VT = -67.0 # Spike threshold adjustment parameter (mV) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h', 'm3h'], 'i_K\ kin.': ['n'], 'i_{TS}\ kin.': ['s', 'u', 's2u'], - 'I': ['iNa', 'iK', 'iTs', 'iL', 'iNet'] + 'I': ['iNa', 'iK', 'iTs', 'iLeak', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' # Instantiate parent class super().__init__() # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) def sinf(self, Vm): ''' Compute the asymptotic value of the open-probability of the S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1.0 + np.exp(-(Vm + 52.0) / 7.4)) # prob def taus(self, Vm): ''' Compute the decay time constant for adaptation of S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' return (1 + 0.33 / (np.exp((Vm + 27.0) / 10.0) + np.exp(-(Vm + 102.0) / 15.0))) * 1e-3 # s def uinf(self, Vm): ''' Compute the asymptotic value of the open-probability of the U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1.0 + np.exp((Vm + 80.0) / 5.0)) # prob def tauu(self, Vm): ''' Compute the decay time constant for adaptation of U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' return (28.3 + 0.33 / (np.exp((Vm + 48.0) / 4.0) + np.exp(-(Vm + 407.0) / 50.0))) * 1e-3 # s class ThalamoCortical(Thalamic): ''' Specific membrane channel dynamics of a thalamo-cortical neuron, with a specific hyperpolarization-activated, mixed cationic current and a leakage Potassium current. References: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* *Destexhe, A., Bal, T., McCormick, D.A., and Sejnowski, T.J. (1996). Ionic mechanisms underlying synchronized oscillations and propagating waves in a model of ferret thalamic slices. J. Neurophysiol. 76, 2049–2070.* *McCormick, D.A., and Huguenard, J.R. (1992). A model of the electrophysiological properties of thalamocortical relay neurons. J. Neurophysiol. 68, 1384–1400.* ''' # Name of channel mechanism name = 'TC' # Cell-specific biophysical parameters # Vm0 = -63.4 # Cell membrane resting potential (mV) Vm0 = -61.93 # Cell membrane resting potential (mV) GNaMax = 900.0 # Max. conductance of Sodium current (S/m^2) GKMax = 100.0 # Max. conductance of Potassium current (S/m^2) GTMax = 20.0 # Max. conductance of low-threshold Calcium current (S/m^2) GKL = 0.138 # Conductance of leakage Potassium current (S/m^2) GhMax = 0.175 # Max. conductance of mixed cationic current (S/m^2) - GL = 0.1 # Conductance of non-specific leakage current (S/m^2) + GLeak = 0.1 # Conductance of non-specific leakage current (S/m^2) Vh = -40.0 # Mixed cationic current reversal potential (mV) - VL = -70.0 # Non-specific leakage Nernst potential (mV) + VLeak = -70.0 # Non-specific leakage Nernst potential (mV) VT = -52.0 # Spike threshold adjustment parameter (mV) Vx = 0.0 # Voltage-dependence uniform shift factor at 36°C (mV) tau_Ca_removal = 5e-3 # decay time constant for intracellular Ca2+ dissolution (s) CCa_min = 50e-9 # minimal intracellular Calcium concentration (M) deff = 100e-9 # effective depth beneath membrane for intracellular [Ca2+] calculation F_Ca = 1.92988e5 # Faraday constant for bivalent ion (Coulomb / mole) nCa = 4 # number of Calcium binding sites on regulating factor k1 = 2.5e22 # intracellular Ca2+ regulation factor (M-4 s-1) k2 = 0.4 # intracellular Ca2+ regulation factor (s-1) k3 = 100.0 # intracellular Ca2+ regulation factor (s-1) k4 = 1.0 # intracellular Ca2+ regulation factor (s-1) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_K\ kin.': ['n'], 'i_{T}\ kin.': ['s', 'u'], 'i_{H}\ kin.': ['O', 'OL', 'O + 2OL'], - 'I': ['iNa', 'iK', 'iT', 'iH', 'iKL', 'iL', 'iNet'] + 'I': ['iNa', 'iK', 'iT', 'iH', 'iKL', 'iLeak', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' # Instantiate parent class super().__init__() # Compute current to concentration conversion constant self.iT_2_CCa = 1e-6 / (self.deff * self.F_Ca) # Define names of the channels state probabilities self.states_names += ['O', 'C', 'P0', 'C_Ca'] # Define the names of the different coefficients to be averaged in a lookup table. self.coeff_names += ['alphao', 'betao'] # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) def sinf(self, Vm): ''' Compute the asymptotic value of the open-probability of the S-type, activation gate of Calcium channels. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1.0 + np.exp(-(Vm + self.Vx + 57.0) / 6.2)) # prob def taus(self, Vm): ''' Compute the decay time constant for adaptation of S-type, activation gate of Calcium channels. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' tmp = np.exp(-(Vm + self.Vx + 132.0) / 16.7) + np.exp((Vm + self.Vx + 16.8) / 18.2) return 1.0 / 3.7 * (0.612 + 1.0 / tmp) * 1e-3 # s def uinf(self, Vm): ''' Compute the asymptotic value of the open-probability of the U-type, inactivation gate of Calcium channels. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1.0 + np.exp((Vm + self.Vx + 81.0) / 4.0)) # prob def tauu(self, Vm): ''' Compute the decay time constant for adaptation of U-type, inactivation gate of Calcium channels. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' if Vm + self.Vx < -80.0: return 1.0 / 3.7 * np.exp((Vm + self.Vx + 467.0) / 66.6) * 1e-3 # s else: return 1 / 3.7 * (np.exp(-(Vm + self.Vx + 22) / 10.5) + 28.0) * 1e-3 # s def derS(self, Vm, s): ''' Compute the evolution of the open-probability of the S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :param s: open-probability of S-type Calcium activation gates (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.sinf(Vm) - s) / self.taus(Vm) def derU(self, Vm, u): ''' Compute the evolution of the open-probability of the U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :param u: open-probability of U-type Calcium inactivation gates (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.uinf(Vm) - u) / self.tauu(Vm) def oinf(self, Vm): ''' Voltage-dependent steady-state activation of hyperpolarization-activated cation current channels. Reference: *Huguenard, J.R., and McCormick, D.A. (1992). Simulation of the currents involved in rhythmic oscillations in thalamic relay neurons. J. Neurophysiol. 68, 1373–1383.* :param Vm: membrane potential (mV) :return: steady-state activation (-) ''' return 1.0 / (1.0 + np.exp((Vm + 75.0) / 5.5)) def tauo(self, Vm): ''' Time constant for activation of hyperpolarization-activated cation current channels. Reference: *Huguenard, J.R., and McCormick, D.A. (1992). Simulation of the currents involved in rhythmic oscillations in thalamic relay neurons. J. Neurophysiol. 68, 1373–1383.* :param Vm: membrane potential (mV) :return: time constant (s) ''' return 1 / (np.exp(-14.59 - 0.086 * Vm) + np.exp(-1.87 + 0.0701 * Vm)) * 1e-3 def alphao(self, Vm): ''' Transition rate between closed and open form of hyperpolarization-activated cation current channels. :param Vm: membrane potential (mV) :return: transition rate (s-1) ''' return self.oinf(Vm) / self.tauo(Vm) def betao(self, Vm): ''' Transition rate between open and closed form of hyperpolarization-activated cation current channels. :param Vm: membrane potential (mV) :return: transition rate (s-1) ''' return (1 - self.oinf(Vm)) / self.tauo(Vm) def derC(self, C, O, Vm): ''' Compute the evolution of the proportion of hyperpolarization-activated cation current channels in closed state. Kinetics scheme of Calcium dependent activation derived from: *Destexhe, A., Bal, T., McCormick, D.A., and Sejnowski, T.J. (1996). Ionic mechanisms underlying synchronized oscillations and propagating waves in a model of ferret thalamic slices. J. Neurophysiol. 76, 2049–2070.* :param Vm: membrane potential (mV) :param C: proportion of Ih channels in closed state (-) :param O: proportion of Ih channels in open state (-) :return: derivative of proportion w.r.t. time (s-1) ''' return self.betao(Vm) * O - self.alphao(Vm) * C def derO(self, C, O, P0, Vm): ''' Compute the evolution of the proportion of hyperpolarization-activated cation current channels in open state. Kinetics scheme of Calcium dependent activation derived from: *Destexhe, A., Bal, T., McCormick, D.A., and Sejnowski, T.J. (1996). Ionic mechanisms underlying synchronized oscillations and propagating waves in a model of ferret thalamic slices. J. Neurophysiol. 76, 2049–2070.* :param Vm: membrane potential (mV) :param C: proportion of Ih channels in closed state (-) :param O: proportion of Ih channels in open state (-) :param P0: proportion of Ih channels regulating factor in unbound state (-) :return: derivative of proportion w.r.t. time (s-1) ''' return - self.derC(C, O, Vm) - self.k3 * O * (1 - P0) + self.k4 * (1 - O - C) def derP0(self, P0, C_Ca): ''' Compute the evolution of the proportion of Ih channels regulating factor in unbound state. Kinetics scheme of Calcium dependent activation derived from: *Destexhe, A., Bal, T., McCormick, D.A., and Sejnowski, T.J. (1996). Ionic mechanisms underlying synchronized oscillations and propagating waves in a model of ferret thalamic slices. J. Neurophysiol. 76, 2049–2070.* :param Vm: membrane potential (mV) :param P0: proportion of Ih channels regulating factor in unbound state (-) :param C_Ca: Calcium concentration in effective submembranal space (M) :return: derivative of proportion w.r.t. time (s-1) ''' return self.k2 * (1 - P0) - self.k1 * P0 * C_Ca**self.nCa def derC_Ca(self, C_Ca, ICa): ''' Compute the evolution of the Calcium concentration in submembranal space. Model of Ca2+ buffering and contribution from iCa derived from: *McCormick, D.A., and Huguenard, J.R. (1992). A model of the electrophysiological properties of thalamocortical relay neurons. J. Neurophysiol. 68, 1384–1400.* :param Vm: membrane potential (mV) :param C_Ca: Calcium concentration in submembranal space (M) :param ICa: inward Calcium current filling up the submembranal space with Ca2+ (mA/m2) :return: derivative of Calcium concentration in submembranal space w.r.t. time (s-1) ''' return (self.CCa_min - C_Ca) / self.tau_Ca_removal - self.iT_2_CCa * ICa def currKL(self, Vm): ''' Compute the voltage-dependent leak Potassium current per unit area. :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GKL * (Vm - self.VK) def currH(self, O, C, Vm): ''' Compute the outward mixed cationic current per unit area. :param O: proportion of the channels in open form :param OL: proportion of the channels in locked-open form :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' OL = 1 - O - C return self.GhMax * (O + 2 * OL) * (Vm - self.Vh) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, u, O, C, _, _ = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currCa(s, u, Vm) + - self.currKL(Vm) + self.currH(O, C, Vm) + self.currL(Vm)) # mA/m2 + self.currKL(Vm) + self.currH(O, C, Vm) + self.currLeak(Vm)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Call parent method to compute Sodium, Potassium and Calcium channels gates steady-states NaKCa_eqstates = super().steadyStates(Vm) # Compute steady-state Calcium current seq = NaKCa_eqstates[3] ueq = NaKCa_eqstates[4] iTeq = self.currCa(seq, ueq, Vm) # Compute steady-state variables for the kinetics system of Ih CCa_eq = self.CCa_min - self.tau_Ca_removal * self.iT_2_CCa * iTeq P0_eq = self.k2 / (self.k2 + self.k1 * CCa_eq**self.nCa) BA = self.betao(Vm) / self.alphao(Vm) O_eq = self.k4 / (self.k3 * (1 - P0_eq) + self.k4 * (1 + BA)) C_eq = BA * O_eq kin_eqstates = np.array([O_eq, C_eq, P0_eq, CCa_eq]) # Merge all steady-states and return return np.concatenate((NaKCa_eqstates, kin_eqstates)) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, u, O, C, P0, C_Ca = states NaKCa_states = [m, h, n, s, u] NaKCa_derstates = super().derStates(Vm, NaKCa_states) dO_dt = self.derO(C, O, P0, Vm) dC_dt = self.derC(C, O, Vm) dP0_dt = self.derP0(P0, C_Ca) ICa = self.currCa(s, u, Vm) dCCa_dt = self.derC_Ca(C_Ca, ICa) return NaKCa_derstates + [dO_dt, dC_dt, dP0_dt, dCCa_dt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute effective coefficients for Sodium, Potassium and Calcium conductances NaKCa_effrates = super().getEffRates(Vm) # Compute effective coefficients for Ih conductance ao_avg = np.mean(self.alphao(Vm)) bo_avg = np.mean(self.betao(Vm)) iH_effrates = np.array([ao_avg, bo_avg]) # Return array of coefficients return np.concatenate((NaKCa_effrates, iH_effrates)) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) Vmeff = np.interp(Qm, interp_data['Q'], interp_data['V']) # Unpack states m, h, n, s, u, O, C, P0, C_Ca = states # INa, IK, ICa effective states derivatives dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s dudt = rates[8] * (1 - u) - rates[9] * u # Ih effective states derivatives dC_dt = rates[11] * O - rates[10] * C dO_dt = - dC_dt - self.k3 * O * (1 - P0) + self.k4 * (1 - O - C) dP0_dt = self.derP0(P0, C_Ca) ICa_eff = self.currCa(s, u, Vmeff) dCCa_dt = self.derC_Ca(C_Ca, ICa_eff) # Merge derivatives and return return [dmdt, dhdt, dndt, dsdt, dudt, dO_dt, dC_dt, dP0_dt, dCCa_dt] diff --git a/PySONIC/plt/pltvars.py b/PySONIC/plt/pltvars.py index bbd7107..f8071bb 100644 --- a/PySONIC/plt/pltvars.py +++ b/PySONIC/plt/pltvars.py @@ -1,614 +1,732 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-21 14:33:36 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-05-25 22:01:35 +# @Last Modified time: 2018-11-30 10:32:06 ''' Dictionary of plotting settings for output variables of the model. ''' pltvars = { 't_ms': { 'desc': 'time', 'label': 'time', 'unit': 'ms', 'factor': 1e3, 'onset': 3e-3 }, 't_us': { 'desc': 'time', 'label': 'time', 'unit': 'us', 'factor': 1e6, 'onset': 1e-6 }, 'Z': { 'desc': 'leaflets deflection', 'label': 'Z', 'unit': 'nm', 'factor': 1e9, 'min': -1.0, 'max': 10.0 }, 'ng': { 'desc': 'gas content', 'label': 'n_g', 'unit': '10^{-22}\ mol', 'factor': 1e22, 'min': 1.0, 'max': 15.0 # 'unit': 'ymol', # 'factor': 1e24, # 'min': 100.0, # 'max': 1500.0 }, 'Pac': { 'desc': 'acoustic pressure', 'label': 'P_{AC}', 'unit': 'kPa', 'factor': 1e-3, 'alias': 'bls.Pacoustic(t, meta["Adrive"] * states, Fdrive)' }, 'Pmavg': { 'desc': 'average intermolecular pressure', 'label': 'P_M', 'unit': 'kPa', 'factor': 1e-3, 'alias': 'bls.PMavgpred(df["Z"].values)' }, 'Telastic': { 'desc': 'leaflet elastic tension', 'label': 'T_E', 'unit': 'mN/m', 'factor': 1e3, 'alias': 'bls.TEleaflet(df["Z"].values)' }, 'Qm': { 'desc': 'charge density', # 'label': 'Q_m', 'label': 'charge', 'unit': 'nC/cm^2', 'factor': 1e5, 'min': -100, 'max': 50 }, 'Cm': { 'desc': 'membrane capacitance', 'label': 'C_m', 'unit': 'uF/cm^2', 'factor': 1e2, 'min': 0.0, 'max': 1.5, 'alias': 'bls.v_Capct(df["Z"].values)' }, 'Vm': { 'desc': 'membrane potential', 'label': 'V_m', 'unit': 'mV', 'factor': 1, }, 'm': { 'desc': 'iNa activation gate opening', 'label': 'm-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'h': { 'desc': 'iNa inactivation gate opening', 'label': 'h-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'm2h': { 'desc': 'iNa relative conductance', 'label': 'm^2h', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': 'df["m"].values**2 * df["h"].values' }, 'm3h': { 'desc': 'iNa relative conductance', 'label': 'm^3h', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': 'df["m"].values**3 * df["h"].values' }, 'm4h': { 'desc': 'iNa relative conductance', 'label': 'm^4h', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': 'df["m"].values**4 * df["h"].values' }, 'n': { 'desc': 'iK activation gate opening', 'label': 'n-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, + 'n4': { + 'desc': 'iK relative conductance', + 'label': 'n^4', + 'unit': None, + 'factor': 1, + 'min': -0.1, + 'max': 1.1, + 'alias': 'df["n"].values**4' + }, + 'p': { 'desc': 'iM activation gate opening', 'label': 'p-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 's': { 'desc': 'iCa activation gates opening', 'label': 's-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'u': { 'desc': 'iCa inactivation gates opening', 'label': 'u-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 's2u': { 'desc': 'iT relative conductance', 'label': 's^2u', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': 'df["s"].values**2 * df["u"].values' }, + 'p': { + 'desc': 'iT activation gates opening', + 'label': 'p-gate', + 'unit': None, + 'factor': 1, + 'min': -0.1, + 'max': 1.1 + }, 'q': { - 'desc': 'iCaL activation gates opening', + 'desc': 'iT inactivation gates opening', 'label': 'q-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, + 'p2q': { + 'desc': 'iT relative conductance', + 'label': 'p^2q', + 'unit': None, + 'factor': 1, + 'min': -0.1, + 'max': 1.1, + 'alias': 'df["p"].values**2 * df["q"].values' + }, + 'r': { - 'desc': 'iCaL inactivation gates opening', + 'desc': 'iCaK Ca2+-gated activation gates opening', 'label': 'r-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, - 'q2r': { - 'desc': 'iCaL relative conductance', - 'label': 'q^2r', + 'r2': { + 'desc': 'iCaK relative conductance', + 'label': 'r^2', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, - 'alias': 'df["q"].values**2 * df["r"].values' + 'alias': 'df["r"].values**2' + }, + + # 'q': { + # 'desc': 'iCaL activation gates opening', + # 'label': 'q-gate', + # 'unit': None, + # 'factor': 1, + # 'min': -0.1, + # 'max': 1.1 + # }, + + # 'r': { + # 'desc': 'iCaL inactivation gates opening', + # 'label': 'r-gate', + # 'unit': None, + # 'factor': 1, + # 'min': -0.1, + # 'max': 1.1 + # }, + + # 'q2r': { + # 'desc': 'iCaL relative conductance', + # 'label': 'q^2r', + # 'unit': None, + # 'factor': 1, + # 'min': -0.1, + # 'max': 1.1, + # 'alias': 'df["q"].values**2 * df["r"].values' + # }, + + # 'c': { + # 'desc': 'iKCa activation gates opening', + # 'label': 'c-gate', + # 'unit': None, + # 'factor': 1, + # 'min': -0.1, + # 'max': 1.1 + # }, + + 'a': { + 'desc': 'iA activation gates opening', + 'label': 'a-gate', + 'unit': None, + 'factor': 1, + 'min': -0.1, + 'max': 1.1 }, + 'b': { + 'desc': 'iA inactivation gates opening', + 'label': 'b-gate', + 'unit': None, + 'factor': 1, + 'min': -0.1, + 'max': 1.1 + }, + + 'a2b': { + 'desc': 'iA relative conductance', + 'label': 'ab', + 'unit': None, + 'factor': 1, + 'min': -0.1, + 'max': 1.1, + 'alias': 'df["a"].values**2 * df["b"].values' + }, + + 'c': { - 'desc': 'iKCa activation gates opening', + 'desc': 'iL activation gates opening', 'label': 'c-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, - 'a': { - 'desc': 'iKA activation gates opening', - 'label': 'a-gate', + 'd1': { + 'desc': 'iL inactivation gates opening', + 'label': 'd1-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, - 'b': { - 'desc': 'iKA inactivation gates opening', - 'label': 'b-gate', + 'd2': { + 'desc': 'iL Ca2+-gated inactivation gates opening', + 'label': 'd2-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, - 'ab': { - 'desc': 'iKA relative conductance', - 'label': 'ab', + 'c2d1d2': { + 'desc': 'iL relative conductance', + 'label': 'c^2d_1d_2', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, - 'alias': 'df["a"].values * df["b"].values' + 'alias': 'df["c"].values**2 * df["d1"].values * df["d2"].values' }, 'O': { 'desc': 'iH activation gate opening', 'label': 'O', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'OL': { 'desc': 'iH activation gate locked-opening', 'label': 'O_L', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': '1 - df["O"].values - df["C"].values' }, 'O + 2OL': { 'desc': 'iH activation gate relative conductance', 'label': 'O\ +\ 2O_L', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 2.1, 'alias': 'df["O"].values + 2 * (1 - df["O"].values - df["C"].values)' }, 'P0': { 'desc': 'iH regulating factor activation', 'label': 'P_0', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'C_Ca': { 'desc': 'sumbmembrane Ca2+ concentration', 'label': '[Ca^{2+}]_i', 'unit': 'uM', 'factor': 1e6 # 'min': 0, # 'max': 150.0 }, 'C_Na': { 'desc': 'sumbmembrane Na+ concentration', 'label': '[Na^+]_i', 'unit': 'uM', 'factor': 1e6 # 'min': 0, # 'max': 150.0 }, 'C_Na_arb': { 'key': 'C_Na', 'desc': 'submembrane Na+ concentration', 'label': '[Na^+]', 'unit': 'arb.', 'factor': 1 }, 'C_Na_arb_activation': { 'key': 'A_Na', 'desc': 'Na+ dependent PumpNa current activation', 'label': 'A_{Na^+}', 'unit': 'arb', 'factor': 1 }, 'C_Ca_arb': { 'key': 'C_Ca', 'desc': 'submembrane Ca2+ concentration', 'label': '[Ca^{2+}]', 'unit': 'arb.', 'factor': 1 }, 'C_Ca_arb_activation': { 'key': 'A_Ca', 'desc': 'Ca2+ dependent Potassium current activation', 'label': 'A_{Ca^{2+}}', 'unit': 'arb', 'factor': 1 }, - 'VL': { - 'constant': 'neuron.VL', + 'VLeak': { + 'constant': 'neuron.VLeak', 'desc': 'non-specific leakage current resting potential', - 'label': 'V_L', + 'label': 'V_{leak}', 'unit': 'mV', 'factor': 1e0 }, - 'iL': { + 'iLeak': { 'desc': 'leakage current', - 'label': 'I_L', + 'label': 'I_{Leak}', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currL(df["Vm"].values)' + 'alias': 'neuron.currLeak(df["Vm"].values)' }, 'iNa': { 'desc': 'Sodium current', 'label': 'I_{Na}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currNa(df["m"].values, df["h"].values, df["Vm"].values)' }, 'iNa2': { 'desc': 'Sodium current', 'label': 'I_{Na}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currNa(df["m"].values, df["h"].values, df["Vm"].values, df["C_Na"].values)' }, 'iK': { 'desc': 'delayed-rectifier Potassium current', 'label': 'I_K', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currK(df["n"].values, df["Vm"].values)' }, 'iM': { 'desc': 'slow non-inactivating Potassium current', 'label': 'I_M', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currM(df["p"].values, df["Vm"].values)' }, 'iA': { 'desc': 'transient Potassium current', 'label': 'I_A', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currA(df["a"].values, df["b"].values, df["Vm"].values)' }, 'iT': { 'desc': 'low-threshold Calcium current', 'label': 'I_T', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currCa(df["s"].values, df["u"].values, df["Vm"].values)' }, + 'iT2': { + 'desc': 'low-threshold Calcium current', + 'label': 'I_T', + 'unit': 'A/m^2', + 'factor': 1e-3, + 'alias': 'neuron.currT(df["p"].values, df["q"].values, df["Vm"].values)' + }, + + 'iL': { + 'desc': 'high-threshold Calcium current', + 'label': 'I_L', + 'unit': 'A/m^2', + 'factor': 1e-3, + 'alias': 'neuron.currL(df["c"].values, df["d1"].values, df["d2"].values, df["Vm"].values)' + }, + + 'iCaK': { + 'desc': 'Calcium-activated Potassium current', + 'label': 'I_{CaK}', + 'unit': 'A/m^2', + 'factor': 1e-3, + 'alias': 'neuron.currCaK(df["r"].values, df["Vm"].values)' + }, + 'iCaL': { 'desc': 'L-type Calcium current', 'label': 'I_{CaL}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currCaL(df["q"].values, df["r"].values, df["Vm"].values)' }, 'iTs': { 'desc': 'low-threshold Calcium current', 'label': 'I_{TS}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currCa(df["s"].values, df["u"].values, df["Vm"].values)' }, 'iCa': { 'desc': 'leech Calcium current', 'label': 'I_{Ca}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currCa(df["s"].values, df["Vm"].values)' }, 'iCa2': { 'desc': 'leech Calcium current', 'label': 'I_{Ca}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currCa(df["s"].values, df["Vm"].values, df["C_Ca"].values)' }, 'iH': { 'desc': 'hyperpolarization-activated cationic current', 'label': 'I_h', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currH(df["O"].values, df["C"].values, df["Vm"].values)' }, 'iKL': { 'desc': 'leakage Potassium current', 'label': 'I_{KL}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currKL(df["Vm"].values)' }, 'iKCa': { 'desc': 'Calcium-activated Potassium current', 'label': 'I_{KCa}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currKCa(df["A_Ca"].values, df["Vm"].values)' }, 'iKCa2': { 'desc': 'Calcium-activated Potassium current', 'label': 'I_{KCa}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currKCa(df["c"].values, df["Vm"].values)' }, 'iPumpNa': { 'desc': 'Outward current mimicking the activity of the NaK-ATPase pump', 'label': 'I_{PumpNa}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currPumpNa(df["A_Na"].values, df["Vm"].values)' }, 'iPumpNa2': { 'desc': 'Outward current mimicking the activity of the NaK-ATPase pump', 'label': 'I_{PumpNa}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currPumpNa(df["C_Na"].values)' }, 'iPumpCa2': { 'desc': 'Outward current describing the removal of Ca2+ from the intracellular space', 'label': 'I_{PumpCa}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currPumpCa(df["C_Ca"].values)' }, 'iNet': { 'desc': 'net current', 'label': 'I_{net}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.currNet(df["Vm"].values, neuron_states)' }, 'Veff': { 'key': 'V', 'desc': 'effective membrane potential', 'label': 'V_{m, eff}', 'unit': 'mV', 'factor': 1e0 }, 'alpham': { 'desc': 'iNa m-gate activation rate', 'label': '\\alpha_m', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betam': { 'desc': 'iNa m-gate inactivation rate', 'label': '\\beta_m', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphah': { 'desc': 'iNa h-gate activation rate', 'label': '\\alpha_h', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betah': { 'desc': 'iNa h-gate inactivation rate', 'label': '\\beta_h', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphan': { 'desc': 'iK n-gate activation rate', 'label': '\\alpha_n', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betan': { 'desc': 'iK n-gate inactivation rate', 'label': '\\beta_n', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphap': { 'desc': 'iM p-gate activation rate', 'label': '\\alpha_p', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betap': { 'desc': 'iM p-gate inactivation rate', 'label': '\\beta_p', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphas': { 'desc': 'iT s-gate activation rate', 'label': '\\alpha_s', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betas': { 'desc': 'iT s-gate inactivation rate', 'label': '\\beta_s', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphau': { 'desc': 'iT u-gate activation rate', 'label': '\\alpha_u', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betau': { 'desc': 'iT u-gate inactivation rate', 'label': '\\beta_u', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphao': { 'desc': 'iH channels activation rate (between closed and open forms)', 'label': '\\alpha_O', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betao': { 'desc': 'iH channels inactivation rate (between closed and open forms)', 'label': '\\beta_O', 'unit': 'ms^{-1}', 'factor': 1e-3 } } diff --git a/scripts/run_estim.py b/scripts/run_estim.py index 8412ef2..9cd46e9 100644 --- a/scripts/run_estim.py +++ b/scripts/run_estim.py @@ -1,114 +1,114 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-24 11:55:07 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-11-22 17:52:42 +# @Last Modified time: 2018-11-30 10:04:25 ''' Run E-STIM simulations of a specific point-neuron. ''' import os import logging import numpy as np import matplotlib.pyplot as plt from argparse import ArgumentParser from PySONIC.utils import logger, selectDirDialog from PySONIC.neurons import * from PySONIC.batches import createEStimQueue, runBatch from PySONIC.plt import plotBatch # Default parameters defaults = dict( neuron='RS', amp=[10.0], # mA/m2 duration=[100.0], # ms PRF=[100.0], # Hz DC=[100.0], # % offset=[50.], # ms method='sonic' ) def runEStimBatch(outdir, neuron, stim_params, mpi=False): ''' Run batch E-STIM simulations of the system for various neuron types and stimulation parameters. :param outdir: full path to output directory :param stim_params: dictionary containing sweeps for all stimulation parameters :param mpi: boolean statting wether or not to use multiprocessing :return: list of full paths to the output files ''' mandatory_params = ['durations', 'offsets', 'PRFs', 'DCs'] for mparam in mandatory_params: if mparam not in stim_params: raise ValueError('Missing stimulation parameter field: "{}"'.format(mparam)) logger.info("Starting E-STIM simulation batch") # Generate simulations queue queue = createEStimQueue( stim_params.get('amps', [None]), stim_params['durations'], stim_params['offsets'], stim_params['PRFs'], stim_params['DCs'] ) # Run batch return runBatch(neuron, 'runAndSave', queue, extra_params=[outdir], mpi=mpi) def main(): ap = ArgumentParser() # Runtime options ap.add_argument('--mpi', default=False, action='store_true', help='Use multiprocessing') ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity') ap.add_argument('-p', '--plot', default=False, action='store_true', help='Plot results') ap.add_argument('-o', '--outputdir', type=str, default=None, help='Output directory') ap.add_argument('-t', '--titrate', default=False, action='store_true', help='Perform titration') # Stimulation parameters ap.add_argument('-n', '--neuron', type=str, default=defaults['neuron'], help='Neuron name (string)') ap.add_argument('-A', '--amp', nargs='+', type=float, help='Injected current density (mA/m2)') ap.add_argument('-d', '--duration', nargs='+', type=float, help='Stimulus duration (ms)') ap.add_argument('--offset', nargs='+', type=float, help='Offset duration (ms)') ap.add_argument('--PRF', nargs='+', type=float, help='PRF (Hz)') ap.add_argument('--DC', nargs='+', type=float, help='Duty cycle (%%)') # Parse arguments args = {key: value for key, value in vars(ap.parse_args()).items() if value is not None} loglevel = logging.DEBUG if args['verbose'] is True else logging.INFO logger.setLevel(loglevel) outdir = args['outputdir'] if 'outputdir' in args else selectDirDialog() mpi = args['mpi'] plot = args['plot'] titrate = args['titrate'] neuron_str = args['neuron'] stim_params = dict( amps=np.array(args.get('amp', defaults['amp'])), # mA/m2 durations=np.array(args.get('duration', defaults['duration'])) * 1e-3, # s PRFs=np.array(args.get('PRF', defaults['PRF'])), # Hz DCs=np.array(args.get('DC', defaults['DC'])) * 1e-2, # (-) offsets=np.array(args.get('offset', defaults['offset'])) * 1e-3 # s ) if titrate: stim_params['amps'] = [None] # Run E-STIM batch if neuron_str not in getNeuronsDict(): logger.error('Unknown neuron type: "%s"', neuron_str) return neuron = getNeuronsDict()[neuron_str]() pkl_filepaths = runEStimBatch(outdir, neuron, stim_params, mpi=mpi) pkl_dir, _ = os.path.split(pkl_filepaths[0]) # Plot resulting profiles if plot: - plotBatch(pkl_filepaths, vars_dict={'V_m': ['Vm']}) + plotBatch(pkl_filepaths) plt.show() if __name__ == '__main__': main()