diff --git a/PySONIC/core/pneuron.py b/PySONIC/core/pneuron.py index 34e83e4..c098db3 100644 --- a/PySONIC/core/pneuron.py +++ b/PySONIC/core/pneuron.py @@ -1,567 +1,567 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2017-08-03 11:53:04 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2020-07-03 10:51:15 +# @Last Modified time: 2020-07-15 19:59:11 import abc import inspect import numpy as np from .protocols import * from .model import Model from .lookups import EffectiveVariablesLookup from .solvers import EventDrivenSolver from .drives import Drive, ElectricDrive from ..postpro import detectSpikes, computeFRProfile from ..constants import * from ..utils import * class PointNeuron(Model): ''' Generic point-neuron model interface. ''' tscale = 'ms' # relevant temporal scale of the model simkey = 'ESTIM' # keyword used to characterize simulations made with this model def __repr__(self): return self.__class__.__name__ def copy(self): return self.__class__() def __eq__(self, other): if not isinstance(other, PointNeuron): return False return self.name == other.name @property @classmethod @abc.abstractmethod def name(cls): ''' Neuron name. ''' raise NotImplementedError @property @classmethod @abc.abstractmethod def Cm0(cls): ''' Neuron's resting capacitance (F/m2). ''' raise NotImplementedError @property @classmethod @abc.abstractmethod def Vm0(cls): ''' Neuron's resting membrane potential(mV). ''' raise NotImplementedError @property def Qm0(self): return self.Cm0 * self.Vm0 * 1e-3 # C/m2 @property def tau_pas(self): ''' Passive membrane time constant (s). ''' return self.Cm0 / self.gLeak @property def meta(self): return {'neuron': self.name} @staticmethod def inputs(): return ElectricDrive.inputs() @classmethod def filecodes(cls, drive, pp): return { 'simkey': cls.simkey, 'neuron': cls.name, 'nature': pp.nature, **drive.filecodes, **pp.filecodes } @classmethod def normalizedQm(cls, Qm): ''' Compute membrane charge density normalized by resting capacitance. :param Qm: membrane charge density (Q/m2) :return: normalized charge density (mV) ''' return Qm / cls.Cm0 * 1e3 @classmethod def getPltVars(cls, wrapleft='df["', wrapright='"]'): pltvars = { 'Qm': { 'desc': 'membrane charge density', 'label': 'Q_m', 'unit': 'nC/cm^2', 'factor': 1e5, 'bounds': ((cls.Vm0 - 20.0) * cls.Cm0 * 1e2, 60) }, 'Qm/Cm0': { 'desc': 'membrane charge density over resting capacitance', 'label': 'Q_m / C_{m0}', 'unit': 'mV', 'bounds': (-150, 70), 'func': f"normalizedQm({wrapleft}Qm{wrapright})" }, 'Vm': { 'desc': 'membrane potential', 'label': 'V_m', 'unit': 'mV', 'bounds': (-150, 70) }, 'ELeak': { 'constant': 'obj.ELeak', 'desc': 'non-specific leakage current resting potential', 'label': 'V_{leak}', 'unit': 'mV', 'ls': '--', 'color': 'k' } } for cname in cls.getCurrentsNames(): cfunc = getattr(cls, cname) cargs = inspect.getargspec(cfunc)[0][1:] pltvars[cname] = { 'desc': inspect.getdoc(cfunc).splitlines()[0], 'label': f'I_{{{cname[1:]}}}', 'unit': 'A/m^2', 'factor': 1e-3, 'func': f"{cname}({', '.join([f'{wrapleft}{a}{wrapright}' for a in cargs])})" } for var in cargs: if var != 'Vm': pltvars[var] = { 'desc': cls.states[var], 'label': var, 'bounds': (-0.1, 1.1) } pltvars['iNet'] = { 'desc': inspect.getdoc(getattr(cls, 'iNet')).splitlines()[0], 'label': 'I_{net}', 'unit': 'A/m^2', 'factor': 1e-3, 'func': f'iNet({wrapleft}Vm{wrapright}, {wrapleft[:-1]}{cls.statesNames()}{wrapright[1:]})', 'ls': '--', 'color': 'black' } pltvars['dQdt'] = { 'desc': inspect.getdoc(getattr(cls, 'dQdt')).splitlines()[0], 'label': 'dQ_m/dt', 'unit': 'A/m^2', 'factor': 1e-3, 'func': f'dQdt({wrapleft}Vm{wrapright}, {wrapleft[:-1]}{cls.statesNames()}{wrapright[1:]})', 'ls': '--', 'color': 'black' } pltvars['iCap'] = { 'desc': inspect.getdoc(getattr(cls, 'iCap')).splitlines()[0], 'label': 'I_{cap}', 'unit': 'A/m^2', 'factor': 1e-3, 'func': f'iCap({wrapleft}t{wrapright}, {wrapleft}Vm{wrapright})' } for rate in cls.rates: if 'alpha' in rate: prefix, suffix = 'alpha', rate[5:] else: prefix, suffix = 'beta', rate[4:] pltvars[rate] = { 'label': '\\{}_{{{}}}'.format(prefix, suffix), 'unit': 'ms^{-1}', 'factor': 1e-3 } pltvars['FR'] = { 'desc': 'riring rate', 'label': 'FR', 'unit': 'Hz', 'factor': 1e0, # 'bounds': (0, 1e3), 'func': f'firingRateProfile({wrapleft[:-2]})' } return pltvars @classmethod def iCap(cls, t, Vm): ''' Capacitive current. ''' dVdt = np.insert(np.diff(Vm) / np.diff(t), 0, 0.) return cls.Cm0 * dVdt @property def pltScheme(self): pltscheme = { 'Q_m': ['Qm'], 'V_m': ['Vm'] } pltscheme['I'] = self.getCurrentsNames() + ['iNet'] for cname in self.getCurrentsNames(): if 'Leak' not in cname: key = f'i_{{{cname[1:]}}}\ kin.' cargs = inspect.getargspec(getattr(self, cname))[0][1:] pltscheme[key] = [var for var in cargs if var not in ['Vm', 'Cai']] return pltscheme @classmethod def statesNames(cls): ''' Return a list of names of all state variables of the model. ''' return list(cls.states.keys()) @classmethod @abc.abstractmethod def derStates(cls): ''' Dictionary of states derivatives functions ''' raise NotImplementedError @classmethod def getDerStates(cls, Vm, states): ''' Compute states derivatives array given a membrane potential and states dictionary ''' return np.array([cls.derStates()[k](Vm, states) for k in cls.statesNames()]) @classmethod @abc.abstractmethod def steadyStates(cls): ''' Return a dictionary of steady-states functions ''' raise NotImplementedError @classmethod def getSteadyStates(cls, Vm): ''' Compute array of steady-states for a given membrane potential ''' return np.array([cls.steadyStates()[k](Vm) for k in cls.statesNames()]) @classmethod def getDerEffStates(cls, lkp, states): ''' Compute effective states derivatives array given lookups and states dictionaries. ''' return np.array([cls.derEffStates()[k](lkp, states) for k in cls.statesNames()]) @classmethod def getEffRates(cls, Vm): ''' Compute array of effective rate constants for a given membrane potential vector. ''' return {k: np.mean(np.vectorize(v)(Vm)) for k, v in cls.effRates().items()} def getLookup(self): ''' Get lookup of membrane potential rate constants interpolated along the neuron's charge physiological range. ''' Qmin, Qmax = expandRange(*self.Qbounds, exp_factor=10.) Qref = np.arange(Qmin, Qmax, 1e-5) # C/m2 Vref = Qref / self.Cm0 * 1e3 # mV tables = {k: np.vectorize(v)(Vref) for k, v in self.effRates().items()} return EffectiveVariablesLookup({'Q': Qref}, {'V': Vref, **tables}) @classmethod @abc.abstractmethod def currents(cls): ''' Dictionary of ionic currents functions (returning current densities in mA/m2) ''' @classmethod def iNet(cls, Vm, states): ''' net membrane current :param Vm: membrane potential (mV) :states: states of ion channels gating and related variables :return: current per unit area (mA/m2) ''' return sum([cfunc(Vm, states) for cfunc in cls.currents().values()]) @classmethod def dQdt(cls, Vm, states): ''' membrane charge density variation rate :param Vm: membrane potential (mV) :states: states of ion channels gating and related variables :return: variation rate (mA/m2) ''' return -cls.iNet(Vm, states) @classmethod def titrationFunc(cls, *args, **kwargs): ''' Default titration function. ''' return cls.isExcited(*args, **kwargs) @staticmethod def currentToConcentrationRate(z_ion, depth): ''' Compute the conversion factor from a specific ionic current (in mA/m2) into a variation rate of submembrane ion concentration (in M/s). :param: z_ion: ion valence :param depth: submembrane depth (m) :return: conversion factor (Mmol.m-1.C-1) ''' return 1e-6 / (z_ion * depth * FARADAY) @staticmethod def nernst(z_ion, Cion_in, Cion_out, T): ''' Nernst potential of a specific ion given its intra and extracellular concentrations. :param z_ion: ion valence :param Cion_in: intracellular ion concentration :param Cion_out: extracellular ion concentration :param T: temperature (K) :return: ion Nernst potential (mV) ''' return (Rg * T) / (z_ion * FARADAY) * np.log(Cion_out / Cion_in) * 1e3 @staticmethod def vtrap(x, y): ''' Generic function used to compute rate constants. ''' return x / (np.exp(x / y) - 1) @staticmethod def efun(x): ''' Generic function used to compute rate constants. ''' return x / (np.exp(x) - 1) @classmethod def ghkDrive(cls, Vm, Z_ion, Cion_in, Cion_out, T): ''' Use the Goldman-Hodgkin-Katz equation to compute the electrochemical driving force of a specific ion species for a given membrane potential. :param Vm: membrane potential (mV) :param Cin: intracellular ion concentration (M) :param Cout: extracellular ion concentration (M) :param T: temperature (K) :return: electrochemical driving force of a single ion particle (mC.m-3) ''' x = Z_ion * FARADAY * Vm / (Rg * T) * 1e-3 # [-] eCin = Cion_in * cls.efun(-x) # M eCout = Cion_out * cls.efun(x) # M return FARADAY * (eCin - eCout) * 1e6 # mC/m3 @classmethod def xBG(cls, Vref, Vm): ''' Compute dimensionless Borg-Graham ratio for a given voltage. :param Vref: reference voltage membrane (mV) :param Vm: membrane potential (mV) :return: dimensionless ratio ''' return (Vm - Vref) * FARADAY / (Rg * cls.T) * 1e-3 # [-] @classmethod def alphaBG(cls, alpha0, zeta, gamma, Vref, Vm): ''' Compute the activation rate constant for a given voltage and temperature, using a Borg-Graham formalism. :param alpha0: pre-exponential multiplying factor :param zeta: effective valence of the gating particle :param gamma: normalized position of the transition state within the membrane :param Vref: membrane voltage at which alpha = alpha0 (mV) :param Vm: membrane potential (mV) :return: rate constant (in alpha0 units) ''' return alpha0 * np.exp(-zeta * gamma * cls.xBG(Vref, Vm)) @classmethod def betaBG(cls, beta0, zeta, gamma, Vref, Vm): ''' Compute the inactivation rate constant for a given voltage and temperature, using a Borg-Graham formalism. :param beta0: pre-exponential multiplying factor :param zeta: effective valence of the gating particle :param gamma: normalized position of the transition state within the membrane :param Vref: membrane voltage at which beta = beta0 (mV) :param Vm: membrane potential (mV) :return: rate constant (in beta0 units) ''' return beta0 * np.exp(zeta * (1 - gamma) * cls.xBG(Vref, Vm)) @classmethod def getCurrentsNames(cls): return list(cls.currents().keys()) @staticmethod def firingRateProfile(*args, **kwargs): return computeFRProfile(*args, **kwargs) @property def Qbounds(self): ''' Determine bounds of membrane charge physiological range for a given neuron. ''' return np.array([np.round(self.Vm0 - 25.0), 50.0]) * self.Cm0 * 1e-3 # C/m2 @classmethod def isVoltageGated(cls, state): ''' Determine whether a given state is purely voltage-gated or not.''' return f'alpha{state.lower()}' in cls.rates @classmethod @Model.checkOutputDir def simQueue(cls, amps, durations, offsets, PRFs, DCs, **kwargs): ''' Create a serialized 2D array of all parameter combinations for a series of individual parameter sweeps. :param amps: list (or 1D-array) of acoustic amplitudes :param durations: list (or 1D-array) of stimulus durations :param offsets: list (or 1D-array) of stimulus offsets (paired with durations array) :param PRFs: list (or 1D-array) of pulse-repetition frequencies :param DCs: list (or 1D-array) of duty cycle values :return: list of parameters (list) for each simulation ''' if amps is None: amps = [None] drives = ElectricDrive.createQueue(amps) protocols = PulsedProtocol.createQueue(durations, offsets, PRFs, DCs) queue = [] for drive in drives: for pp in protocols: queue.append([drive, pp]) return queue @classmethod @Model.checkOutputDir def simQueueBurst(cls, amps, durations, PRFs, DCs, BRFs, nbursts, **kwargs): if amps is None: amps = [None] drives = ElectricDrive.createQueue(amps) protocols = BurstProtocol.createQueue(durations, PRFs, DCs, BRFs, nbursts) queue = [] for drive in drives: for pp in protocols: queue.append([drive, pp]) return queue @staticmethod def checkInputs(drive, pp): ''' Check validity of electrical stimulation parameters. :param drive: electric drive object :param pp: pulse protocol object ''' if not isinstance(drive, Drive): raise TypeError(f'Invalid "drive" parameter (must be an "Drive" object)') if not isinstance(pp, TimeProtocol): raise TypeError('Invalid time protocol (must be "TimeProtocol" instance)') def chooseTimeStep(self): ''' Determine integration time step based on intrinsic temporal properties. ''' return DT_EFFECTIVE @classmethod def derivatives(cls, t, y, Cm=None, drive=None): ''' Compute system derivatives for a given membrane capacitance and injected current. :param t: specific instant in time (s) :param y: vector of HH system variables at time t :param Cm: membrane capacitance (F/m2) :param Iinj: injected current (mA/m2) :return: vector of system derivatives at time t ''' if Cm is None: Cm = cls.Cm0 Qm, *states = y Vm = Qm / Cm * 1e3 # mV states_dict = dict(zip(cls.statesNames(), states)) dQmdt = - cls.iNet(Vm, states_dict) # mA/m2 if drive is not None: dQmdt += drive.compute(t) dQmdt *= 1e-3 # A/m2 # dQmdt = (Iinj - cls.iNet(Vm, states_dict)) * 1e-3 # A/m2 return [dQmdt, *cls.getDerStates(Vm, states_dict)] @Model.logNSpikes @Model.checkTitrate @Model.addMeta @Model.logDesc @Model.checkSimParams def simulate(self, drive, pp): ''' Simulate a specific neuron model for a set of simulation parameters, and return output data in a dataframe. :param drive: electric drive object :param pp: pulse protocol object :return: output DataFrame ''' # Set initial conditions y0 = { 'Qm': self.Qm0, **{k: self.steadyStates()[k](self.Vm0) for k in self.statesNames()} } # Initialize solver and compute solution solver = EventDrivenSolver( lambda x: setattr(solver.drive, 'xvar', drive.xvar * x), # eventfunc y0.keys(), # variables lambda t, y: self.derivatives(t, y, drive=solver.drive), # dfunc event_params={'drive': drive.copy().updatedX(0.)}, # event parameters dt=self.chooseTimeStep()) # time step data = solver(y0, pp.stimEvents(), pp.tstop) # Add Vm timeries to solution data = addColumn(data, 'Vm', data['Qm'].values / self.Cm0 * 1e3, preceding_key='Qm') # Return solution dataframe return data def desc(self, meta): return f'{self}: simulation @ {meta["drive"].desc}, {meta["pp"].desc}' @staticmethod def getNSpikes(data): ''' Compute number of spikes in charge profile of simulation output. :param data: dataframe containing output time series :return: number of detected spikes ''' return detectSpikes(data)[0].size @staticmethod def getStabilizationValue(data): ''' Determine stabilization value from the charge profile of a simulation output. :param data: dataframe containing output time series :return: charge stabilization value (or np.nan if no stabilization detected) ''' # Extract charge signal posterior to observation window t, Qm = [data[key].values for key in ['t', 'Qm']] if t.max() <= TMIN_STABILIZATION: raise ValueError('solution length is too short to assess stabilization') Qm = Qm[t > TMIN_STABILIZATION] # Compute variation range Qm_range = np.ptp(Qm) logger.debug('%.2f nC/cm2 variation range over the last %.0f ms, Qmf = %.2f nC/cm2', Qm_range * 1e5, TMIN_STABILIZATION * 1e3, Qm[-1] * 1e5) # Return final value only if stabilization is detected if np.ptp(Qm) < QSS_Q_DIV_THR: return Qm[-1] else: return np.nan @classmethod def isExcited(cls, data): ''' Determine if neuron is excited from simulation output. :param data: dataframe containing output time series :return: boolean stating whether neuron is excited or not ''' return cls.getNSpikes(data) > 0 @classmethod def isSilenced(cls, data): ''' Determine if neuron is silenced from simulation output. :param data: dataframe containing output time series :return: boolean stating whether neuron is silenced or not ''' return not np.isnan(cls.getStabilizationValue(data)) def getArange(self, drive): return drive.xvar_range diff --git a/PySONIC/neurons/mrg.py b/PySONIC/neurons/mrg.py index 16a05fb..8dd444b 100644 --- a/PySONIC/neurons/mrg.py +++ b/PySONIC/neurons/mrg.py @@ -1,207 +1,172 @@ # -*- coding: utf-8 -*- # @Author: Mariia Popova # @Email: theo.lemaire@epfl.ch # @Date: 2020-02-27 21:24:05 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2020-07-15 10:36:43 +# @Last Modified time: 2020-07-15 21:23:41 import numpy as np from ..core import PointNeuron, addSonicFeatures @addSonicFeatures class MRGNode(PointNeuron): ''' Mammalian myelinated fiber node. Reference: *McIntyre, C.C., Richardson, A.G., and Grill, W.M. (2002). Modeling the excitability of mammalian nerve fibers: influence of afterpotentials on the recovery cycle. J. Neurophysiol. 87, 995–1006.* ''' # Mechanism name name = 'MRGnode' # ------------------------------ Biophysical parameters ------------------------------ # Resting parameters Cm0 = 2e-2 # Membrane capacitance (F/m2) Vm0 = -80. # Membrane potential (mV) # Reversal potentials (mV) ENa = 50. # Sodium EK = -90. # Potassium ELeak = -90. # Non-specific leakage # Maximal channel conductances (S/m2) gNafbar = 3e4 # Fast Sodium gNapbar = 100. # Persistent Sodium gKsbar = 800. # Slow Potassium gLeak = 70. # Non-specific leakage # Additional parameters celsius = 36.0 # Temperature (Celsius) celsius_Schwarz = 20.0 # Temperature in Schwarz 1995 (Celsius) celsius_Ks = 36.0 # Temperature iused for Ks channels (unknown ref.) mhshift = 3. # m and h gates voltage shift (mV) vtraub = -80. # Reference voltage for the definition of the s rate constants (mV) # ------------------------------ States names & descriptions ------------------------------ states = { 'm': 'iNaf activation gate', 'h': 'iNaf inactivation gate', 'p': 'iNap activation gate', 's': 'iKs activation gate', } # ------------------------------ Gating states kinetics ------------------------------ def __new__(cls): cls.q10_mp = 2.2**((cls.celsius - cls.celsius_Schwarz) / 10) # from Schwarz 1987 cls.q10_h = 2.9**((cls.celsius - cls.celsius_Schwarz) / 10) # from Schwarz 1987 cls.q10_s = 3.0**((cls.celsius - cls.celsius_Ks) / 10) # from ??? return super(MRGNode, cls).__new__(cls) # iNaf kinetics: adapted from Schwarz 1995, with notable changes: # - Q10 correction to account for temperature adaptation from 20 degrees reference # - 3 mV offset to m and h gates to shift voltage dependence in the hyperpolarizing direction, # to increase the conduction velocity and strength-duration time constant (from McIntyre 2002) # - increase in tau_h to slow down inactivation (achieved through alphah?) @classmethod def alpham(cls, Vm): Vm += cls.mhshift return cls.q10_mp * 1.86 * cls.vtrap(-(Vm + 18.4), 10.3) * 1e3 # s-1 @classmethod def betam(cls, Vm): Vm += cls.mhshift return cls.q10_mp * 0.086 * cls.vtrap(Vm + 22.7, 9.16) * 1e3 # s-1 @classmethod def alphah(cls, Vm): Vm += cls.mhshift return cls.q10_h * 0.062 * cls.vtrap(Vm + 111.0, 11.0) * 1e3 # s-1 @classmethod def betah(cls, Vm): Vm += cls.mhshift return cls.q10_h * 2.3 / (1 + np.exp(-(Vm + 28.8) / 13.4)) * 1e3 # s-1 # iNap kinetics: adapted from ???, with notable changes: # - Q10 correction to account for temperature adaptation from 20 degrees reference # - increase in tau_p in order to extend the duration and amplitude of the DAP. @classmethod def alphap(cls, Vm): return cls.q10_mp * 0.01 * cls.vtrap(-(Vm + 27.), 10.2) * 1e3 # s-1 @classmethod def betap(cls, Vm): return cls.q10_mp * 0.00025 * cls.vtrap(Vm + 34., 10.) * 1e3 # s-1 # iKs kinetics: adapted from ???, with notable changes: # - Q10 correction to account for temperature adaptation from 36 degrees reference # - increase in tau_s in order to to create an AHP that matches experimental records. @classmethod def alphas(cls, Vm): Vm -= cls.vtraub return cls.q10_s * 0.3 / (1 + np.exp(-(Vm - 27.) / 5.)) * 1e3 # s-1 @classmethod def betas(cls, Vm): Vm -= cls.vtraub return cls.q10_s * 0.03 / (1 + np.exp(-(Vm + 10.) / 1.)) * 1e3 # s-1 # ------------------------------ States derivatives ------------------------------ @classmethod def derStates(cls): return { 'm': lambda Vm, x: cls.alpham(Vm) * (1 - x['m']) - cls.betam(Vm) * x['m'], 'h': lambda Vm, x: cls.alphah(Vm) * (1 - x['h']) - cls.betah(Vm) * x['h'], 'p': lambda Vm, x: cls.alphap(Vm) * (1 - x['p']) - cls.betap(Vm) * x['p'], 's': lambda Vm, x: cls.alphas(Vm) * (1 - x['s']) - cls.betas(Vm) * x['s'], } # ------------------------------ Steady states ------------------------------ @classmethod def steadyStates(cls): return { 'm': lambda Vm: cls.alpham(Vm) / (cls.alpham(Vm) + cls.betam(Vm)), 'h': lambda Vm: cls.alphah(Vm) / (cls.alphah(Vm) + cls.betah(Vm)), 'p': lambda Vm: cls.alphap(Vm) / (cls.alphap(Vm) + cls.betap(Vm)), 's': lambda Vm: cls.alphas(Vm) / (cls.alphas(Vm) + cls.betas(Vm)), } # ------------------------------ Membrane currents ------------------------------ @classmethod def iNaf(cls, m, h, Vm): ''' fast Sodium current. ''' return cls.gNafbar * m**3 * h * (Vm - cls.ENa) # mA/m2 @classmethod def iNap(cls, p, Vm): ''' persistent Sodium current. ''' return cls.gNapbar * p**3 * (Vm - cls.ENa) # mA/m2 @classmethod def iKs(cls, s, Vm): ''' slow Potassium current ''' return cls.gKsbar * s * (Vm - cls.EK) # mA/m2 @classmethod def iLeak(cls, Vm): ''' non-specific leakage current ''' return cls.gLeak * (Vm - cls.ELeak) # mA/m2 @classmethod def currents(cls): return { 'iNaf': lambda Vm, x: cls.iNaf(x['m'], x['h'], Vm), 'iNap': lambda Vm, x: cls.iNap(x['p'], Vm), 'iKs': lambda Vm, x: cls.iKs(x['s'], Vm), 'iLeak': lambda Vm, _: cls.iLeak(Vm) } def chooseTimeStep(self): ''' neuron-specific time step for fast dynamics. ''' return super().chooseTimeStep() * 1e-2 - - -@addSonicFeatures -class SilentMRGNode(MRGNode): - ''' MRG node with an additional axial leakage current to mimick the influence of - the neighboring internodal sections in preventing spontaneous activity. - ''' - name = 'silent_MRGnode' - - # Axial leakage conductance yielding identical threshold for intracellular stimulation - # at the central node of a 10 um myelinated fiber with a 100 us pulse - gAxialLeak = 5.37e3 # S/m2 - - # Corresponding code snippet: - ''' - from PySONIC.core import PulsedProtocol - from ExSONIC.core import MRGFiber, IntracellularCurrent - fiber = MRGFiber(10e-6, 11) - Ithr = fiber.titrate( - IntracellularCurrent(fiber.central_ID, None), - PulsedProtocol(tstim=1e-4, toffset=3e-3)) # A - Am = fiber.nodes[fiber.central_ID].Am # cm2 - ithr = (Ithr * 1e3) / (Am * 1e-4) # mA/m2 - print(f'Ithr = {Ithr * 1e9:.2f} nA, Am = {Am:.2e} cm2 -> ithr = {ithr:.2e} mA/m2') - ''' - - @classmethod - def iAxialLeak(cls, Vm): - ''' axial leakage current ''' - return cls.gAxialLeak * (Vm - cls.Vm0) # mA/m2 - - @classmethod - def currents(cls): - ''' updated currents dictionary ''' - return {**MRGNode.currents(), 'iAxialLeak': lambda Vm, _: cls.iAxialLeak(Vm)}