diff --git a/PySONIC/core/nbls.py b/PySONIC/core/nbls.py index 8b2b1a8..5a0b26c 100644 --- a/PySONIC/core/nbls.py +++ b/PySONIC/core/nbls.py @@ -1,698 +1,693 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2016-09-29 16:16:19 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-25 17:53:56 +# @Last Modified time: 2019-06-25 18:00:35 from copy import deepcopy import logging import numpy as np import pandas as pd from scipy.integrate import solve_ivp from .simulators import PWSimulator, HybridSimulator, PeriodicSimulator from .bls import BilayerSonophore from .pneuron import PointNeuron from .model import Model from .batches import createQueue from ..neurons import getNeuronLookup from ..utils import * from ..constants import * from ..postpro import getFixedPoints class NeuronalBilayerSonophore(BilayerSonophore): ''' This class inherits from the BilayerSonophore class and receives an PointNeuron instance at initialization, to define the electro-mechanical NICE model and its SONIC variant. ''' tscale = 'ms' # relevant temporal scale of the model simkey = 'ASTIM' # keyword used to characterize simulations made with this model def __init__(self, a, pneuron, Fdrive=None, embedding_depth=0.0): ''' Constructor of the class. :param a: in-plane radius of the sonophore structure within the membrane (m) :param pneuron: point-neuron model :param Fdrive: frequency of acoustic perturbation (Hz) :param embedding_depth: depth of the embedding tissue around the membrane (m) ''' # Check validity of input parameters if not isinstance(pneuron, PointNeuron): raise ValueError('Invalid neuron type: "{}" (must inherit from PointNeuron class)' .format(pneuron.name)) self.pneuron = pneuron - # Parse pneuron object to create derEffStates method (along with corresponding expression) - self.pneuron.parse() - # self.pneuron.printDerEffStates() - # self.pneuron.printEffRates() - # Initialize BilayerSonophore parent object BilayerSonophore.__init__(self, a, pneuron.Cm0, pneuron.Qm0, embedding_depth) def __repr__(self): s = '{}({:.1f} nm, {}'.format(self.__class__.__name__, self.a * 1e9, self.pneuron) if self.d > 0.: s += ', d={}m'.format(si_format(self.d, precision=1, space=' ')) return s + ')' def params(self): return {**super().params(), **self.pneuron.params()} def getPltVars(self, wrapleft='df["', wrapright='"]'): return {**super().getPltVars(wrapleft, wrapright), **self.pneuron.getPltVars(wrapleft, wrapright)} def getPltScheme(self): return self.pneuron.getPltScheme() def filecode(self, *args): return Model.filecode(self, *args) def inputVars(self): # Get parent input vars and supress irrelevant entries bls_vars = super().inputVars() pneuron_vars = self.pneuron.inputVars() del bls_vars['Qm'] del pneuron_vars['Astim'] # Fill in current input vars in appropriate order inputvars = bls_vars inputvars.update(pneuron_vars) inputvars['method'] = None return inputvars def filecodes(self, Fdrive, Adrive, tstim, toffset, PRF, DC, method='sonic'): # Get parent codes and supress irrelevant entries bls_codes = super().filecodes(Fdrive, Adrive, 0.0) pneuron_codes = self.pneuron.filecodes(0.0, tstim, toffset, PRF, DC) for x in [bls_codes, pneuron_codes]: del x['simkey'] del bls_codes['Qm'] del pneuron_codes['Astim'] # Fill in current codes in appropriate order codes = { 'simkey': self.simkey, 'neuron': pneuron_codes.pop('neuron'), 'nature': pneuron_codes.pop('nature') } codes.update(bls_codes) codes.update(pneuron_codes) codes['method'] = method return codes def fullDerivatives(self, t, y, Adrive, Fdrive, phi): ''' Compute the derivatives of the (n+3) ODE full NBLS system variables. :param t: specific instant in time (s) :param y: vector of state variables :param Adrive: acoustic drive amplitude (Pa) :param Fdrive: acoustic drive frequency (Hz) :param phi: acoustic drive phase (rad) :return: vector of derivatives ''' dydt_mech = BilayerSonophore.derivatives(self, t, y[:3], Adrive, Fdrive, y[3], phi) dydt_elec = self.pneuron.Qderivatives(t, y[3:], self.Capct(y[1])) return dydt_mech + dydt_elec def effDerivatives(self, t, y, lkp1D): ''' Compute the derivatives of the n-ODE effective HH system variables, based on 1-dimensional linear interpolation of "effective" coefficients that summarize the system's behaviour over an acoustic cycle. :param t: specific instant in time (s) :param y: vector of HH system variables at time t :param lkp: dictionary of 1D data points of "effective" coefficients over the charge domain, for specific frequency and amplitude values. :return: vector of effective system derivatives at time t ''' Qm, *states = y states_dict = dict(zip(self.pneuron.statesNames(), states)) lkp0D = lkp1D.interpolate1D('Q', Qm) dQmdt = - self.pneuron.iNet(lkp0D['V'], states_dict) * 1e-3 return [dQmdt, *self.pneuron.getDerEffStates(lkp0D, states_dict)] def interpOnOffVariable(self, key, Qm, stim, lkps): ''' Interpolate Q-dependent effective variable along ON and OFF periods of a solution. :param key: lookup variable key :param Qm: charge density solution vector :param stim: stimulation state solution vector :param lkp: dictionary of lookups for ON and OFF states :return: interpolated effective variable vector ''' x = np.zeros(stim.size) x[stim == 0] = lkps['OFF'].interpVar(Qm[stim == 0], 'Q', key) x[stim == 1] = lkps['ON'].interpVar(Qm[stim == 1], 'Q', key) return x def runFull(self, Fdrive, Adrive, tstim, toffset, PRF, DC, phi=np.pi): ''' Compute solutions of the full electro-mechanical system for a specific set of US stimulation parameters, using a classic integration scheme. The first iteration uses the quasi-steady simplification to compute the initiation of motion from a flat leaflet configuration. Afterwards, the ODE system is solved iteratively until completion. :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :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 phi: acoustic drive phase (rad) :return: 2-tuple with the output dataframe and computation time. ''' # Determine time step dt = 1 / (NPC_DENSE * Fdrive) # Compute non-zero deflection value for a small perturbation (solving quasi-steady equation) Pac = self.Pacoustic(dt, Adrive, Fdrive, phi) Z0 = self.balancedefQS(self.ng0, self.Qm0, Pac) # Set initial conditions y0 = np.concatenate(( [0., Z0, self.ng0, self.Qm0], self.pneuron.getSteadyStates(self.pneuron.Vm0))) # Initialize simulator and compute solution logger.debug('Computing detailed solution') simulator = PWSimulator( lambda t, y: self.fullDerivatives(t, y, Adrive, Fdrive, phi), lambda t, y: self.fullDerivatives(t, y, 0., 0., 0.)) (t, y, stim), tcomp = simulator( y0, dt, tstim, toffset, PRF, DC, print_progress=logger.getEffectiveLevel() <= logging.INFO, target_dt=CLASSIC_TARGET_DT, monitor_time=True) logger.debug('completed in %ss', si_format(tcomp, 1)) # Store output in dataframe data = pd.DataFrame({ 't': t, 'stimstate': stim, 'Z': y[:, 1], 'ng': y[:, 2], 'Qm': y[:, 3] }) data['Vm'] = data['Qm'].values / self.v_Capct(data['Z'].values) * 1e3 # mV for i in range(len(self.pneuron.states)): data[self.pneuron.statesNames()[i]] = y[:, i + 4] # Return dataframe and computation time return data, tcomp def runHybrid(self, Fdrive, Adrive, tstim, toffset, PRF, DC, phi=np.pi): ''' Compute solutions of the system for a specific set of US stimulation parameters, using a hybrid integration scheme. :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param tstim: duration of US stimulation (s) :param toffset: duration of the offset (s) :param phi: acoustic drive phase (rad) :return: 3-tuple with the time profile, the solution matrix and a state vector ''' # Determine time steps dt_dense, dt_sparse = [1. / (n * Fdrive) for n in [NPC_DENSE, NPC_SPARSE]] # Compute non-zero deflection value for a small perturbation (solving quasi-steady equation) Pac = self.Pacoustic(dt_dense, Adrive, Fdrive, phi) Z0 = self.balancedefQS(self.ng0, self.Qm0, Pac) # Set initial conditions y0 = np.concatenate(( [0., Z0, self.ng0, self.Qm0], self.pneuron.getSteadyStates(self.pneuron.Vm0))) is_dense_var = np.array([True] * 3 + [False] * (len(self.pneuron.states) + 1)) # Initialize simulator and compute solution logger.debug('Computing hybrid solution') simulator = HybridSimulator( lambda t, y: self.fullDerivatives(t, y, Adrive, Fdrive, phi), lambda t, y: self.fullDerivatives(t, y, 0., 0., 0.), lambda t, y, Cm: self.pneuron.Qderivatives(t, y, Cm), lambda yref: self.Capct(yref[1]), is_dense_var, ivars_to_check=[1, 2]) (t, y, stim), tcomp = simulator( y0, dt_dense, dt_sparse, 1. / Fdrive, tstim, toffset, PRF, DC, monitor_time=True) logger.debug('completed in %ss', si_format(tcomp, 1)) # Store output in dataframe data = pd.DataFrame({ 't': t, 'stimstate': stim, 'Z': y[:, 1], 'ng': y[:, 2], 'Qm': y[:, 3] }) data['Vm'] = data['Qm'].values / self.v_Capct(data['Z'].values) * 1e3 # mV for i in range(len(self.pneuron.states)): data[self.pneuron.statesNames()[i]] = y[:, i + 4] # Return dataframe and computation time return data, tcomp def computeEffVars(self, Fdrive, Adrive, Qm, fs): ''' Compute "effective" coefficients of the HH system for a specific combination of stimulus frequency, stimulus amplitude and charge density. A short mechanical simulation is run while imposing the specific charge density, until periodic stabilization. The HH coefficients are then averaged over the last acoustic cycle to yield "effective" coefficients. :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param Qm: imposed charge density (C/m2) :param fs: list of sonophore membrane coverage fractions :return: list with computation time and a list of dictionaries of effective variables ''' # Run simulation and retrieve deflection and gas content vectors from last cycle data, tcomp = BilayerSonophore.simulate(self, Fdrive, Adrive, Qm) Z_last = data.loc[-NPC_DENSE:, 'Z'].values # m Cm_last = self.v_Capct(Z_last) # F/m2 # For each coverage fraction effvars = [] for x in fs: # Compute membrane capacitance and membrane potential vectors Cm = x * Cm_last + (1 - x) * self.Cm0 # F/m2 Vm = Qm / Cm * 1e3 # mV # Compute average cycle value for membrane potential and rate constants effvars.append({**{'V': np.mean(Vm)}, **self.pneuron.getEffRates(Vm)}) # Log process log = '{}: lookups @ {}Hz, {}Pa, {:.2f} nC/cm2'.format( self, *si_format([Fdrive, Adrive], precision=1, space=' '), Qm * 1e5) if len(fs) > 1: log += ', fs = {:.0f} - {:.0f}%'.format(fs.min() * 1e2, fs.max() * 1e2) log += ', tcomp = {:.3f} s'.format(tcomp) logger.info(log) # Return effective coefficients return [tcomp, effvars] def runSONIC(self, Fdrive, Adrive, tstim, toffset, PRF, DC): ''' Compute solutions of the system for a specific set of US stimulation parameters, using charge-predicted "effective" coefficients to solve the HH equations at each step. :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :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 (-) :return: 3-tuple with the time profile, the effective solution matrix and a state vector ''' # Load appropriate 2D lookups lkp2D = getNeuronLookup(self.pneuron.name).projectN({'a': self.a, 'f': Fdrive}) # Interpolate 2D lookups at zero and US amplitude logger.debug('Interpolating lookups at A = %.2f kPa and A = 0', Adrive * 1e-3) lkps1D = {'ON': lkp2D.project('A', Adrive), 'OFF': lkp2D.project('A', 0.)} # Set initial conditions y0 = np.concatenate(([self.Qm0], self.pneuron.getSteadyStates(self.pneuron.Vm0))) # Initialize simulator and compute solution logger.debug('Computing effective solution') simulator = PWSimulator( lambda t, y: self.effDerivatives(t, y, lkps1D['ON']), lambda t, y: self.effDerivatives(t, y, lkps1D['OFF'])) (t, y, stim), tcomp = simulator( y0, DT_EFFECTIVE, tstim, toffset, PRF, DC, monitor_time=True) logger.debug('completed in %ss', si_format(tcomp, 1)) # Store output in dataframe data = pd.DataFrame({ 't': t, 'stimstate': stim, 'Qm': y[:, 0] }) data['Vm'] = self.interpOnOffVariable('V', data['Qm'].values, stim, lkps1D) for key in ['Z', 'ng']: data[key] = np.full(t.size, np.nan) for i in range(len(self.pneuron.states)): data[self.pneuron.statesNames()[i]] = y[:, i + 1] # Return dataframe and computation time return data, tcomp def meta(self, Fdrive, Adrive, tstim, toffset, PRF, DC, method): ''' Return information about object and simulation parameters. :param Fdrive: US frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param tstim: stimulus duration (s) :param toffset: stimulus offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: stimulus duty cycle (-) :param method: integration method :return: meta-data dictionary ''' return { 'simkey': self.simkey, 'neuron': self.pneuron.name, 'a': self.a, 'd': self.d, 'Fdrive': Fdrive, 'Adrive': Adrive, 'tstim': tstim, 'toffset': toffset, 'PRF': PRF, 'DC': DC, 'method': method } def simulate(self, Fdrive, Adrive, tstim, toffset, PRF=100., DC=1.0, method='sonic'): ''' Simulate the electro-mechanical model for a specific set of US stimulation parameters, and return output data in a dataframe. :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :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 method: selected integration method :return: 2-tuple with the output dataframe and computation time. ''' logger.info( '%s: simulation @ f = %sHz, A = %sPa, t = %ss (%ss offset)%s', self, si_format(Fdrive, 0, space=' '), si_format(Adrive, 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 '')) # Check validity of stimulation parameters BilayerSonophore.checkInputs(self, Fdrive, Adrive, 0.0, 0.0) self.pneuron.checkInputs(Adrive, tstim, toffset, PRF, DC) # Call appropriate simulation function try: simfunc = { 'full': self.runFull, 'hybrid': self.runHybrid, 'sonic': self.runSONIC }[method] except KeyError: raise ValueError('Invalid integration method: "{}"'.format(method)) data, tcomp = simfunc(Fdrive, Adrive, tstim, toffset, PRF, DC) # Log number of detected spikes nspikes = self.pneuron.getNSpikes(data) logger.debug('{} spike{} detected'.format(nspikes, plural(nspikes))) # Return dataframe and computation time return data, tcomp @logCache(os.path.join(os.path.split(__file__)[0], 'astim_titrations.log')) def titrate(self, Fdrive, tstim, toffset, PRF=100., DC=1., method='sonic', xfunc=None, Arange=None): ''' Use a binary search to determine the threshold amplitude needed to obtain neural excitation for a given frequency, duration, PRF and duty cycle. :param Fdrive: US frequency (Hz) :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 method: integration method :param xfunc: function determining whether condition is reached from simulation output :param Arange: search interval for Adrive, iteratively refined :return: determined threshold amplitude (Pa) ''' # Default output function if xfunc is None: xfunc = self.pneuron.titrationFunc # Default amplitude interval if Arange is None: Arange = [0., getNeuronLookup(self.pneuron.name).refs['A'].max()] return binarySearch( lambda x: xfunc(self.simulate(*x)[0]), [Fdrive, tstim, toffset, PRF, DC, method], 1, Arange, THRESHOLD_CONV_RANGE_ASTIM ) def simQueue(self, freqs, amps, durations, offsets, PRFs, DCs, methods, outputdir=None): ''' Create a serialized 2D array of all parameter combinations for a series of individual parameter sweeps, while avoiding repetition of CW protocols for a given PRF sweep. :param freqs: list (or 1D-array) of US frequencies :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 :params methods: integration methods :return: list of parameters (list) for each simulation ''' method_ids = list(range(len(methods))) if amps is None: amps = [np.nan] DCs = np.array(DCs) queue = [] if 1.0 in DCs: queue += createQueue(freqs, amps, durations, offsets, min(PRFs), 1.0, method_ids) if np.any(DCs != 1.0): queue += createQueue(freqs, amps, durations, offsets, PRFs, DCs[DCs != 1.0], method_ids) for item in queue: if np.isnan(item[1]): item[1] = None item[-1] = methods[int(item[-1])] if outputdir is not None: for item in queue: item.insert(0, outputdir) return queue def quasiSteadyStates(self, Fdrive, amps=None, charges=None, DCs=1.0, squeeze_output=False): ''' Compute the quasi-steady state values of the neuron's gating variables for a combination of US amplitudes, charge densities and duty cycles, at a specific US frequency. :param Fdrive: US frequency (Hz) :param amps: US amplitudes (Pa) :param charges: membrane charge densities (C/m2) :param DCs: duty cycle value(s) :return: 4-tuple with reference values of US amplitude and charge density, as well as interpolated Vmeff and QSS gating variables ''' # Get DC-averaged lookups interpolated at the appropriate amplitudes and charges amps, charges, lookups = getLookupsDCavg( self.pneuron.name, self.a, Fdrive, amps, charges, DCs) # Compute QSS states using these lookups nA, nQ, nDC = lookups['V'].shape QSS = {k: np.empty((nA, nQ, nDC)) for k in self.pneuron.statesNames()} for iA in range(nA): for iDC in range(nDC): QSS_1D = self.pneuron.quasiSteadyStates( {k: v[iA, :, iDC] for k, v in lookups.items()}) for k in QSS.keys(): QSS[k][iA, :, iDC] = QSS_1D[k] # Compress outputs if needed if squeeze_output: QSS = {k: v.squeeze() for k, v in QSS.items()} lookups = {k: v.squeeze() for k, v in lookups.items()} # Return reference inputs and outputs return amps, charges, lookups, QSS def iNetQSS(self, Qm, Fdrive, Adrive, DC): ''' Compute quasi-steady state net membrane current for a given combination of US parameters and a given membrane charge density. :param Qm: membrane charge density (C/m2) :param Fdrive: US frequency (Hz) :param Adrive: US amplitude (Pa) :param DC: duty cycle (-) :return: net membrane current (mA/m2) ''' _, _, lookups, QSS = self.quasiSteadyStates( Fdrive, amps=Adrive, charges=Qm, DCs=DC, squeeze_output=True) return self.pneuron.iNet(lookups['V'], np.array(list(QSS.values()))) # mA/m2 def evaluateStability(self, Qm0, states0, lkp): ''' Integrate the effective differential system from a given starting point, until clear convergence or clear divergence is found. :param Qm0: initial membrane charge density (C/m2) :param states0: dictionary of initial states values :param lkp: dictionary of 1D data points of "effective" coefficients over the charge domain, for specific frequency and amplitude values. :return: boolean indicating convergence state ''' # Initialize y0 vector t0 = 0. y0 = np.array([Qm0] + list(states0.values())) # Initialize simulator and compute solution # simulator = PeriodicSimulator( # lambda t, y: self.effDerivatives(t, y, lkp), # ivars_to_check=[0]) # simulator.stopfunc = simulator.stopFuncTmp # simulator.refs = [Qm0] # simulator.conv_thr = [QSS_Q_CONV_THR] # simulator.div_thr = [QSS_Q_DIV_THR] # simulator.t_history = QSS_HISTORY_INTERVAL # t, y, stim = simulator.compute( # y0, # QSS_INTEGRATION_INTERVAL, # QSS_INTEGRATION_INTERVAL, # nmax=int(QSS_MAX_INTEGRATION_DURATION // QSS_INTEGRATION_INTERVAL) # ) # conv = simulator.isAsymptoticallyStable(t, y, QSS_INTEGRATION_INTERVAL) != -1 # tf = t[-1] # dQ = [y[-1, 0]] # Initializing empty list to record evolution of charge deviation n = int(QSS_HISTORY_INTERVAL // QSS_INTEGRATION_INTERVAL) # size of history dQ = [] # As long as there is no clear charge convergence or divergence conv, div = False, False tf, yf = t0, y0 while not conv and not div: # Integrate system for small interval and retrieve final charge deviation t0, y0 = tf, yf sol = solve_ivp( lambda t, y: self.effDerivatives(t, y, lkp), [t0, t0 + QSS_INTEGRATION_INTERVAL], y0, method='LSODA' ) tf, yf = sol.t[-1], sol.y[:, -1] dQ.append(yf[0] - Qm0) # logger.debug('{:.0f} ms: dQ = {:.5f} nC/cm2, avg dQ = {:.5f} nC/cm2'.format( # tf * 1e3, dQ[-1] * 1e5, np.mean(dQ[-n:]) * 1e5)) # If last charge deviation is too large -> divergence if np.abs(dQ[-1]) > QSS_Q_DIV_THR: div = True # If last charge deviation or average deviation in recent history # is small enough -> convergence for x in [dQ[-1], np.mean(dQ[-n:])]: if np.abs(x) < QSS_Q_CONV_THR: conv = True # If max integration duration is been reached -> error if tf > QSS_MAX_INTEGRATION_DURATION: logger.warning('too many iterations (dQ = {:.5f} nC/cm2)'.format(dQ[-1] * 1e5)) conv = True logger.debug('{}vergence after {:.0f} ms: dQ = {:.5f} nC/cm2'.format( {True: 'con', False: 'di'}[conv], tf * 1e3, dQ[-1] * 1e5)) return conv def fixedPointsQSS(self, Fdrive, Adrive, DC, lkp, dQdt): ''' Compute QSS fixed points along the charge dimension for a given combination of US parameters, and determine their stability. :param Fdrive: US frequency (Hz) :param Adrive: US amplitude (Pa) :param DC: duty cycle (-) :param lkp: lookup dictionary for effective variables along charge dimension :param dQdt: charge derivative profile along charge dimension :return: 2-tuple with values of stable and unstable fixed points ''' logger.debug('A = {:.2f} kPa, DC = {:.0f}%'.format(Adrive * 1e-3, DC * 1e2)) # Extract stable and unstable fixed points from QSS charge variation profile def dfunc(Qm): return - self.iNetQSS(Qm, Fdrive, Adrive, DC) SFP_candidates = getFixedPoints(lkp['Q'], dQdt, filter='stable', der_func=dfunc).tolist() UFPs = getFixedPoints(lkp['Q'], dQdt, filter='unstable', der_func=dfunc).tolist() SFPs = [] pltvars = self.getPltVars() # For each candidate SFP for i, Qm in enumerate(SFP_candidates): logger.debug('Q-SFP = {:.2f} nC/cm2'.format(Qm * 1e5)) # Re-compute QSS *_, QSS_FP = self.quasiSteadyStates(Fdrive, amps=Adrive, charges=Qm, DCs=DC, squeeze_output=True) # Simulate from unperturbed QSS and evaluate stability if not self.evaluateStability(Qm, QSS_FP, lkp): logger.warning('diverging system at ({:.2f} kPa, {:.2f} nC/cm2)'.format( Adrive * 1e-3, Qm * 1e5)) UFPs.append(Qm) else: # For each state unstable_states = [] for x in self.pneuron.statesNames(): pltvar = pltvars[x] unit_str = pltvar.get('unit', '') factor = pltvar.get('factor', 1) is_stable_direction = [] for sign in [-1, +1]: # Perturb state with small offset QSS_perturbed = deepcopy(QSS_FP) QSS_perturbed[x] *= (1 + sign * QSS_REL_OFFSET) # If gating state, bound within [0., 1.] if self.pneuron.isVoltageGated(x): QSS_perturbed[x] = np.clip(QSS_perturbed[x], 0., 1.) logger.debug('{}: {:.5f} -> {:.5f} {}'.format( x, QSS_FP[x] * factor, QSS_perturbed[x] * factor, unit_str)) # Simulate from perturbed QSS and evaluate stability is_stable_direction.append( self.evaluateStability(Qm, QSS_perturbed, lkp)) # Check if system shows stability upon x-state perturbation # in both directions if not np.all(is_stable_direction): unstable_states.append(x) # Classify fixed point as stable only if all states show stability is_stable_FP = len(unstable_states) == 0 {True: SFPs, False: UFPs}[is_stable_FP].append(Qm) logger.info('{}stable fixed-point at ({:.2f} kPa, {:.2f} nC/cm2){}'.format( '' if is_stable_FP else 'un', Adrive * 1e-3, Qm * 1e5, '' if is_stable_FP else ', caused by {} states'.format(unstable_states))) return SFPs, UFPs def isStableQSS(self, Fdrive, Adrive, DC): _, Qref, lookups, QSS = self.quasiSteadyStates( Fdrive, amps=Adrive, DCs=DC, squeeze_output=True) lookups['Q'] = Qref dQdt = -self.pneuron.iNet( lookups['V'], np.array([QSS[k] for k in self.pneuron.statesNames()])) # mA/m2 SFPs, _ = self.fixedPointsQSS(Fdrive, Adrive, DC, lookups, dQdt) return len(SFPs) > 0 def titrateQSS(self, Fdrive, DC=1., Arange=None): # Default amplitude interval if Arange is None: Arange = [0., getNeuronLookup(self.pneuron.name).refs['A'].max()] # Titration function def xfunc(x): if self.pneuron.name == 'STN': return self.isStableQSS(*x) else: return not self.isStableQSS(*x) return binarySearch( xfunc, [Fdrive, DC], 1, Arange, THRESHOLD_CONV_RANGE_ASTIM) diff --git a/PySONIC/core/pneuron.py b/PySONIC/core/pneuron.py index 9747bb2..44fa8fa 100644 --- a/PySONIC/core/pneuron.py +++ b/PySONIC/core/pneuron.py @@ -1,726 +1,700 @@ # -*- 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: 2019-06-25 17:55:35 +# @Last Modified time: 2019-06-25 18:15:45 import abc import re import pprint import inspect import numpy as np import pandas as pd from .batches import createQueue from .model import Model from .simulators import PWSimulator from ..postpro import findPeaks, computeFRProfile from ..constants import * from ..utils import si_format, logger, plural, binarySearch 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 __init__(self): + # Determine neuron's resting charge density self.Qm0 = self.Cm0 * self.Vm0 * 1e-3 # C/cm2 - if hasattr(self, 'states'): - self.rates = self.getRatesNames(self.statesNames()) + + # Parse pneuron object to create derEffStates and effRates methods + self.parse() + # self.printDerEffStates() + # self.printEffRates() def __repr__(self): return self.__class__.__name__ def inputVars(self): return { 'Astim': { 'desc': 'current density amplitude', 'label': 'A', 'unit': 'mA/m2', 'factor': 1e0, 'precision': 1 }, 'tstim': { 'desc': 'stimulus duration', 'label': 't_{stim}', 'unit': 'ms', 'factor': 1e3, 'precision': 0 }, 'toffset': { 'desc': 'offset duration', 'label': 't_{offset}', 'unit': 'ms', 'factor': 1e3, 'precision': 0 }, 'PRF': { 'desc': 'pulse repetition frequency', 'label': 'PRF', 'unit': 'Hz', 'factor': 1e0, 'precision': 0 }, 'DC': { 'desc': 'duty cycle', 'label': 'DC', 'unit': '%', 'factor': 1e2, 'precision': 2 } } def filecodes(self, Astim, tstim, toffset, PRF, DC): is_CW = DC == 1. return { 'simkey': self.simkey, 'neuron': self.name, 'nature': 'CW' if is_CW else 'PW', 'Astim': '{:.1f}mAm2'.format(Astim), 'tstim': '{:.0f}ms'.format(tstim * 1e3), 'toffset': None, 'PRF': 'PRF{:.2f}Hz'.format(PRF) if not is_CW else None, 'DC': 'DC{:.2f}%'.format(DC * 1e2) if not is_CW else None } @property @abc.abstractmethod def name(self): raise NotImplementedError @property @abc.abstractmethod def Cm0(self): raise NotImplementedError @property @abc.abstractmethod def Vm0(self): raise NotImplementedError def statesNames(self): return list(self.states.keys()) @abc.abstractmethod def steadyStates(self): ''' Dictionary of steady-states functions ''' def getSteadyStates(self, Vm): ''' Compute array of steady-states for a given membrane potential ''' return np.array([self.steadyStates()[k](Vm) for k in self.statesNames()]) @abc.abstractmethod def derStates(self): ''' Dictionary of states derivatives functions ''' def getDerStates(self, Vm, states): ''' Compute states derivatives array given a membrane potential and states dictionary ''' return np.array([self.derStates()[k](Vm, states) for k in self.statesNames()]) def parse(self): ''' Construct a neuron-specific dictionary of effective derivative functions. ''' def getLambdaSource(dict_entry): clean_dict_entry = re.sub( ' +', ' ', inspect.getsource(dict_entry).replace('\n', ' ')).strip() if clean_dict_entry[-1] == ',': clean_dict_entry = clean_dict_entry[:-1] fsource = clean_dict_entry.split(':', 1)[-1].strip() lambda_pattern = 'lambda ([a-z_A-Z,0-9\s]*): (.+)' return re.match(lambda_pattern, fsource).groups() def getFuncCalls(s): func_pattern = '([a-z_A-Z]*).([a-z_A-Z][a-z_A-Z0-9]*)\(([^\)]*)\)' return re.finditer(func_pattern, s) def createStateEffectiveDerivative(expr): lambda_expr = 'lambda lkp, x: {}'.format(expr) lambda_expr_self = 'lambda self, lkp, x: {}'.format(expr) f = eval(lambda_expr_self) return lambda_expr, lambda *args: f(self, *args) def defineConstLambda(const): return lambda _: const def addEffRates(expr, d, dstr): # Define patterns suffix_pattern = '[A-Za-z0-9_]+' xinf_pattern = re.compile('^({})inf$'.format(suffix_pattern)) taux_pattern = re.compile('^tau({})$'.format(suffix_pattern)) alphax_pattern = re.compile('^alpha({})$'.format(suffix_pattern)) betax_pattern = re.compile('^beta({})$'.format(suffix_pattern)) err_str = 'gating states must be defined via the alphaX-betaX or Xinf-tauX paradigm' # If expression matches alpha or beta rate -> return corresponding # effective rate function if alphax_pattern.match(expr) or betax_pattern.match(expr): try: ratex = getattr(self, expr) d[expr] = np.vectorize(lambda Vm: np.mean(ratex(Vm))) except AttributeError: raise ValueError(err_str) dstr[expr] = 'lambda Vm: np.mean(self.{}(Vm))'.format(expr) # If expression matches xinf or taux -> add corresponding alpha and beta # effective rates functions else: for pattern in [taux_pattern, xinf_pattern]: m = pattern.match(expr) if m: k = m.group(1) alphax_str, betax_str = ['{}{}'.format(p, k) for p in ['alpha', 'beta']] xinf_str, taux_str = ['{}inf'.format(k), 'tau{}'.format(k)] try: xinf, taux = [getattr(self, s) for s in [xinf_str, taux_str]] # If taux is a constant, define a lambda function that returns it if not callable(taux): taux = defineConstLambda(taux) d[alphax_str] = np.vectorize( lambda Vm: np.mean(xinf(Vm) / taux(Vm))) d[betax_str] = np.vectorize( lambda Vm: np.mean((1 - xinf(Vm)) / taux(Vm))) except AttributeError: raise ValueError(err_str) dstr.update({ alphax_str: 'lambda Vm: np.mean(self.{}(Vm) / self.{}(Vm))'.format( xinf_str, taux_str), betax_str: 'lambda Vm: np.mean((1 - self.{}(Vm)) / self.{}(Vm))'.format( xinf_str, taux_str) }) # Initialize empty dictionaries to gather effective derivatives functions # and effective rates functions eff_dstates, eff_dstates_str = {}, {} eff_rates, eff_rates_str = {}, {} # For each state derivative for k, dfunc in self.derStates().items(): # Get derivative function source code dfunc_args, dfunc_exp = getLambdaSource(dfunc) # For each internal function call in the derivative function expression matches = getFuncCalls(dfunc_exp) for m in matches: # Determine function arguments fprefix, fname, fargs = m.groups() fcall = '{}({})'.format(fname, fargs) if fprefix: fcall = '{}.{}'.format(fprefix, fcall) args_list = fargs.split(',') # If sole argument is Vm if len(args_list) == 1 and args_list[0] == 'Vm': # Replace function call by lookup retrieval in expression dfunc_exp = dfunc_exp.replace(fcall, "lkp['{}']".format(fname)) # Add the corresponding effective rate function(s) to the dictionnary addEffRates(fname, eff_rates, eff_rates_str) # Replace Vm by lkp['V'] in expression dfunc_exp = dfunc_exp.replace('Vm', "lkp['V']") # Create the modified lambda expression and evaluate it eff_dstates_str[k], eff_dstates[k] = createStateEffectiveDerivative(dfunc_exp) + self.rates = list(eff_rates.keys()) + # Define methods that return dictionaries of effective states derivatives # and effective rates functions, along with their corresponding string equivalents if not hasattr(self, 'derEffStates'): self.derEffStates = lambda: eff_dstates self.printDerEffStates = lambda: pprint.PrettyPrinter(indent=4).pprint(eff_dstates_str) if not hasattr(self, 'effRates'): self.effRates = lambda: eff_rates self.printEffRates = lambda: pprint.PrettyPrinter(indent=4).pprint(eff_rates_str) def getDerEffStates(self, lkp, states): ''' Compute states effective derivatives array given a dictionary of lookup vectors ''' return np.array([ self.derEffStates()[k](lkp, states) for k in self.statesNames()]) def getEffRates(self, Vm): return {k: v(Vm) for k, v in self.effRates().items()} @abc.abstractmethod def currents(self): ''' Dictionary of ionic currents functions (returning current densities in mA/m2) ''' def iNet(self, 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 self.currents().values()]) def dQdt(self, 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 -self.iNet(Vm, states) def titrationFunc(self, *args, **kwargs): ''' Default titration function. ''' return self.isExcited(*args, **kwargs) def currentToConcentrationRate(self, 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) def nernst(self, 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 def vtrap(self, x, y): ''' Generic function used to compute rate constants. ''' return x / (np.exp(x / y) - 1) def efun(self, x): ''' Generic function used to compute rate constants. ''' return x / (np.exp(x) - 1) def ghkDrive(self, 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 * self.efun(-x) # M eCout = Cion_out * self.efun(x) # M return FARADAY * (eCin - eCout) * 1e6 # mC/m3 def getCurrentsNames(self): return list(self.currents().keys()) def getPltScheme(self): pltscheme = { 'Q_m': ['Qm'], 'V_m': ['Vm'] } pltscheme['I'] = self.getCurrentsNames() + ['iNet'] for cname in self.getCurrentsNames(): if 'Leak' not in cname: key = 'i_{{{}}}\ kin.'.format(cname[1:]) cargs = inspect.getargspec(getattr(self, cname))[0][1:] pltscheme[key] = [var for var in cargs if var not in ['Vm', 'Cai']] return pltscheme def getPltVars(self, wrapleft='df["', wrapright='"]'): ''' Return a dictionary with information about all plot variables related to the neuron. ''' pltvars = { 'Qm': { 'desc': 'membrane charge density', 'label': 'Q_m', 'unit': 'nC/cm^2', 'factor': 1e5, 'bounds': (-100, 50) }, 'Vm': { 'desc': 'membrane potential', 'label': 'V_m', 'unit': 'mV', 'y0': self.Vm0, '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 self.getCurrentsNames(): cfunc = getattr(self, cname) cargs = inspect.getargspec(cfunc)[0][1:] pltvars[cname] = { 'desc': inspect.getdoc(cfunc).splitlines()[0], 'label': 'I_{{{}}}'.format(cname[1:]), 'unit': 'A/m^2', 'factor': 1e-3, 'func': '{}({})'.format(cname, ', '.join(['{}{}{}'.format(wrapleft, a, wrapright) for a in cargs])) } for var in cargs: if var != 'Vm': pltvars[var] = { 'desc': self.states[var], 'label': var, 'bounds': (-0.1, 1.1) } pltvars['iNet'] = { 'desc': inspect.getdoc(getattr(self, 'iNet')).splitlines()[0], 'label': 'I_{net}', 'unit': 'A/m^2', 'factor': 1e-3, 'func': 'iNet({0}Vm{1}, {2}{3}{4})'.format( wrapleft, wrapright, wrapleft[:-1], self.statesNames(), wrapright[1:]), 'ls': '--', 'color': 'black' } pltvars['dQdt'] = { 'desc': inspect.getdoc(getattr(self, 'dQdt')).splitlines()[0], 'label': 'dQ_m/dt', 'unit': 'A/m^2', 'factor': 1e-3, 'func': 'dQdt({0}Vm{1}, {2}{3}{4}.values.T)'.format( wrapleft, wrapright, wrapleft[:-1], self.statesNames(), wrapright[1:]), 'ls': '--', 'color': 'black' } - for x in self.getGates(): - for rate in ['alpha', 'beta']: - pltvars['{}{}'.format(rate, x)] = { - 'label': '\\{}_{{{}}}'.format(rate, x), + for rate in self.rates: + if 'alpha' in rate: + prefix, suffix = 'alpha', rate[5:] + else: + prefix, suffix = 'beta', rate[4:] + pltvars['{}'.format(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': 'firingRateProfile({0}t{1}.values, {0}Qm{1}.values)'.format(wrapleft, wrapright) } return pltvars def firingRateProfile(*args, **kwargs): return computeFRProfile(*args, **kwargs) - def getRatesNames(self, states): - ''' Return a list of names of the alpha and beta rates of the neuron. ''' - return list(self.effRates.keys()) - 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 def isVoltageGated(self, state): ''' Determine whether a given state is purely voltage-gated or not.''' return 'alpha{}'.format(state.lower()) in self.rates - def getGates(self): - ''' Retrieve the names of the neuron's states that match an ion channel gating. ''' - gates = [] - for x in self.statesNames(): - if self.isVoltageGated(x): - gates.append(x) - return gates - def qsStates(self, lkp, states): ''' Compute a collection of quasi steady states using the standard xinf = ax / (ax + Bx) equation. :param lkp: dictionary of 1D vectors of "effective" coefficients over the charge domain, for specific frequency and amplitude values. :return: dictionary of quasi-steady states ''' return { x: lkp['alpha{}'.format(x)] / (lkp['alpha{}'.format(x)] + lkp['beta{}'.format(x)]) for x in states } def quasiSteadyStates(self, lkp): ''' Compute the quasi-steady states of a neuron for a range of membrane charge densities, based on 1-dimensional lookups interpolated at a given sonophore diameter, US frequency, US amplitude and duty cycle. :param lkp: dictionary of 1D vectors of "effective" coefficients over the charge domain, for specific frequency and amplitude values. :return: dictionary of quasi-steady states ''' return self.qsStates(lkp, self.statesNames()) # return {k: func(lkp['Vm']) for k, func in self.steadyStates().items()} - 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 Vderivatives(self, t, y, Iinj): ''' Compute the derivatives of a V-cast HH system for a specific value of injected current. :param t: time value (s, unused) :param y: vector of HH system variables at time t :param Iinj: injected current (mA/m2) :return: vector of HH system derivatives at time t ''' Vm, *states = y states_dict = dict(zip(self.statesNames(), states)) Iionic = self.iNet(Vm, states_dict) # mA/m2 dVmdt = (- Iionic + Iinj) / self.Cm0 # mV/s return [dVmdt, *self.getDerStates(Vm, states_dict)] def Qderivatives(self, t, y, Cm=None): ''' Compute the derivatives of the n-ODE HH system variables, based on a value of membrane capacitance. :param t: specific instant in time (s) :param y: vector of HH system variables at time t :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 states_dict = dict(zip(self.statesNames(), states)) dQmdt = - self.iNet(Vm, states_dict) * 1e-3 # A/m2 return [dQmdt, *self.getDerStates(Vm, states_dict)] 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 meta(self, Astim, tstim, toffset, PRF, DC): ''' Return information about object and simulation parameters. :param Astim: stimulus amplitude (mA/m2) :param tstim: stimulus duration (s) :param toffset: stimulus offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: stimulus duty cycle (-) :return: meta-data dictionary ''' return { 'simkey': self.simkey, 'neuron': self.name, 'Astim': Astim, 'tstim': tstim, 'toffset': toffset, 'PRF': PRF, 'DC': DC } def simulate(self, Astim, tstim, toffset, PRF=100., DC=1.0): ''' Simulate a specific neuron model for a specific set of electrical parameters, and return output data in a dataframe. :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: 2-tuple with the output dataframe and computation time. ''' logger.info( '%s: simulation @ A = %sA/m2, t = %ss (%ss offset)%s', self, si_format(Astim, 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 '')) # Check validity of stimulation parameters self.checkInputs(Astim, tstim, toffset, PRF, DC) # Set initial conditions y0 = np.array((self.Vm0, *self.getSteadyStates(self.Vm0))) # Initialize simulator and compute solution logger.debug('Computing solution') simulator = PWSimulator( lambda t, y: self.Vderivatives(t, y, Astim), lambda t, y: self.Vderivatives(t, y, 0.)) (t, y, stim), tcomp = simulator( y0, DT_EFFECTIVE, tstim, toffset, PRF, DC, monitor_time=True) logger.debug('completed in %ss', si_format(tcomp, 1)) # Store output in dataframe data = pd.DataFrame({ 't': t, 'stimstate': stim, 'Vm': y[:, 0], 'Qm': y[:, 0] * self.Cm0 * 1e-3 }) data['Qm'] = data['Vm'].values * self.Cm0 * 1e-3 for i in range(len(self.states)): data[self.statesNames()[i]] = y[:, i + 1] # Log number of detected spikes nspikes = self.getNSpikes(data) logger.debug('{} spike{} detected'.format(nspikes, plural(nspikes))) # Return dataframe and computation time return data, tcomp def simQueue(self, amps, durations, offsets, PRFs, DCs, outputdir=None): ''' Create a serialized 2D array of all parameter combinations for a series of individual parameter sweeps, while avoiding repetition of CW protocols for a given PRF sweep. :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 = [np.nan] DCs = np.array(DCs) queue = [] if 1.0 in DCs: queue += createQueue(amps, durations, offsets, min(PRFs), 1.0) if np.any(DCs != 1.0): queue += createQueue(amps, durations, offsets, PRFs, DCs[DCs != 1.0]) for item in queue: if np.isnan(item[0]): item[0] = None if outputdir is not None: for item in queue: item.insert(0, outputdir) return queue def getNSpikes(self, data): ''' Compute number of spikes in charge profile of simulation output. :param data: dataframe containing output time series :return: number of detected spikes ''' dt = np.diff(data.ix[:1, 't'].values)[0] ipeaks, *_ = findPeaks( data['Qm'].values, SPIKE_MIN_QAMP, int(np.ceil(SPIKE_MIN_DT / dt)), SPIKE_MIN_QPROM ) return ipeaks.size def getStabilizationValue(self, 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 def isExcited(self, 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 self.getNSpikes(data) > 0 def isSilenced(self, 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(self.getStabilizationValue(data)) def titrate(self, tstim, toffset, PRF, DC, xfunc=None, Arange=(0., 2 * AMP_UPPER_BOUND_ESTIM)): ''' Use a binary 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 xfunc: function determining whether condition is reached from simulation output :param Arange: search interval for Astim, iteratively refined :return: excitation threshold amplitude (mA/m2) ''' # Default output function if xfunc is None: xfunc = self.titrationFunc return binarySearch( lambda x: xfunc(self.simulate(*x)[0]), [tstim, toffset, PRF, DC], 0, Arange, THRESHOLD_CONV_RANGE_ESTIM ) diff --git a/PySONIC/neurons/leech.py b/PySONIC/neurons/leech.py index 9200a45..3c2d295 100644 --- a/PySONIC/neurons/leech.py +++ b/PySONIC/neurons/leech.py @@ -1,570 +1,568 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2017-07-31 15:20:54 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-25 17:53:09 +# @Last Modified time: 2019-06-25 18:09:08 from functools import partialmethod import numpy as np from ..core import PointNeuron from ..constants import FARADAY, Rg, Z_Na, Z_Ca class LeechTouch(PointNeuron): ''' Leech touch sensory neuron 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.* ''' # Neuron name name = 'LeechT' # ------------------------------ Biophysical parameters ------------------------------ # Resting parameters Cm0 = 1e-2 # Membrane capacitance (F/m2) Vm0 = -53.58 # Membrane potential (mV) # Reversal potentials (mV) ENa = 45.0 # Sodium EK = -62.0 # Potassium ECa = 60.0 # Calcium ELeak = -48.0 # Non-specific leakage EPumpNa = -300.0 # Sodium pump # Maximal channel conductances (S/m2) gNabar = 3500.0 # Sodium gKdbar = 900.0 # Delayed-rectifier Potassium gCabar = 20.0 # Calcium gKCabar = 236.0 # Calcium-dependent Potassium gLeak = 1.0 # Non-specific leakage gPumpNa = 20.0 # Sodium pump # Activation time constants (s) taum = 0.1e-3 # Sodium taus = 0.6e-3 # Calcium # Original conversion constants from inward ionic 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) taur_Na = 16.0 # Sodium taur_Ca = 1.25 # Calcium # Time constants for the PumpNa and KCa currents activation # from specific intracellular ions (s) taua_PumpNa = 0.1 # PumpNa current activation from intracellular Na+ taua_KCa = 0.01 # KCa current activation from intracellular Ca2+ # ------------------------------ States names & descriptions ------------------------------ states = { 'm': 'iNa activation gate', 'h': 'iNa inactivation gate', 'n': 'iKd gate', 's': 'iCa gate', 'Nai': 'submembrane Na+ concentration (arbitrary unit)', 'ANa': 'Na+ dependent iPumpNa gate', 'Cai': 'submembrane Ca2+ concentration (arbitrary unit)', 'ACa': 'Ca2+ dependent iKCa gate' } def __init__(self): super().__init__() - self.rates = self.getRatesNames(['m', 'h', 'n', 's']) self.K_Na = self.K_Na_original * self.surface * self.curr_factor self.K_Ca = self.K_Ca_original * self.surface * self.curr_factor # ------------------------------ Gating states kinetics ------------------------------ def _xinf(self, Vm, halfmax, slope, power): ''' Generic function computing the steady-state open-probability of a particular ion channel gate at a given voltage. :param Vm: membrane potential (mV) :param halfmax: half-activation voltage (mV) :param slope: slope parameter of activation function (mV) :param power: power exponent multiplying the exponential expression (integer) :return: steady-state open-probability (-) ''' return 1 / (1 + np.exp((Vm - halfmax) / slope))**power def _taux(self, Vm, halfmax, slope, tauMax, tauMin): ''' Generic function computing the voltage-dependent, adaptation time constant of a particular ion channel gate at a given voltage. :param Vm: membrane potential (mV) :param halfmax: voltage at which adaptation time constant is half-maximal (mV) :param slope: slope parameter of adaptation time constant function (mV) :return: adptation time constant (s) ''' return (tauMax - tauMin) / (1 + np.exp((Vm - halfmax) / slope)) + tauMin def _derCion(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 _derAion(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 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) 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) sinf = partialmethod(_xinf, halfmax=-10.0, slope=-2.8, power=1) # ------------------------------ States derivatives ------------------------------ def derNai(self, Nai, m, h, Vm): ''' Evolution of submembrane Sodium concentration ''' return self._derCion(Nai, self.iNa(m, h, Vm), self.K_Na, self.taur_Na) # M/s def derCai(self, Cai, s, Vm): ''' Evolution of submembrane Calcium concentration ''' return self._derCion(Cai, self.iCa(s, Vm), self.K_Ca, self.taur_Ca) # M/s def derANa(self, ANa, Nai): ''' Evolution of Na+ dependent iPumpNa gate ''' return self._derAion(ANa, Nai, self.taua_PumpNa) def derACa(self, ACa, Cai): ''' Evolution of Ca2+ dependent iKCa gate ''' return self._derAion(ACa, Cai, self.taua_KCa) def derStates(self): return { 'm': lambda Vm, x: (self.minf(Vm) - x['m']) / self.taum, 'h': lambda Vm, x: (self.hinf(Vm) - x['h']) / self.tauh(Vm), 'n': lambda Vm, x: (self.ninf(Vm) - x['n']) / self.taun(Vm), 's': lambda Vm, x: (self.sinf(Vm) - x['s']) / self.taus, 'Nai': lambda Vm, x: self.derNai(x['Nai'], x['m'], x['h'], Vm), 'ANa': lambda Vm, x: self.derANa(x['ANa'], x['Nai']), 'Cai': lambda Vm, x: self.derCai(x['Cai'], x['s'], Vm), 'ACa': lambda Vm, x: self.derACa(x['ACa'], x['Cai']) } # ------------------------------ Steady states ------------------------------ def Naiinf(self, Vm): ''' Steady.state Sodium intracellular concentration. ''' return -self.K_Na * self.iNa(self.minf(Vm), self.hinf(Vm), Vm) def Caiinf(self, Vm): ''' Steady.state Calcium intracellular concentration. ''' return -self.K_Ca * self.iCa(self.sinf(Vm), Vm) def steadyStates(self): return { 'm': lambda Vm: self.minf(Vm), 'h': lambda Vm: self.hinf(Vm), 'n': lambda Vm: self.ninf(Vm), 's': lambda Vm: self.sinf(Vm), 'Nai': lambda Vm: self.Naiinf(Vm), 'ANa': lambda Vm: self.Naiinf(Vm), 'Cai': lambda Vm: self.Caiinf(Vm), 'ACa': lambda Vm: self.Caiinf(Vm) } # def quasiSteadyStates(self, lkp): # qsstates = self.qsStates(lkp, ['m', 'h', 'n', 's']) # qsstates['Nai'] = - self.K_Na * self.iNa(qsstates['m'], qsstates['h'], lkp['V']) # qsstates['ANa'] = qsstates['Nai'] # qsstates['Cai'] = - self.K_Ca * self.iCa(qsstates['s'], lkp['V']) # qsstates['ACa'] = qsstates['Cai'] # return qsstates # ------------------------------ Membrane currents ------------------------------ def iNa(self, m, h, Vm): ''' Sodium current ''' return self.gNabar * m**3 * h * (Vm - self.ENa) # mA/m2 def iKd(self, n, Vm): ''' Delayed-rectifier Potassium current ''' return self.gKdbar * n**2 * (Vm - self.EK) # mA/m2 def iCa(self, s, Vm): ''' Calcium current ''' return self.gCabar * s * (Vm - self.ECa) # mA/m2 def iKCa(self, ACa, Vm): ''' Calcium-activated Potassium current ''' return self.gKCabar * ACa * (Vm - self.EK) # mA/m2 def iPumpNa(self, ANa, Vm): ''' NaK-ATPase pump current ''' return self.gPumpNa * ANa * (Vm - self.EPumpNa) # mA/m2 def iLeak(self, Vm): ''' Non-specific leakage current ''' return self.gLeak * (Vm - self.ELeak) # mA/m2 def currents(self): return { 'iNa': lambda Vm, x: self.iNa(x['m'], x['h'], Vm), 'iKd': lambda Vm, x: self.iKd(x['n'], Vm), 'iCa': lambda Vm, x: self.iCa(x['s'], Vm), 'iPumpNa': lambda Vm, x: self.iPumpNa(x['ANa'], Vm), 'iKCa': lambda Vm, x: self.iKCa(x['ACa'], Vm), 'iLeak': lambda Vm, _: self.iLeak(Vm) } class LeechMech(PointNeuron): ''' Generic leech neuron 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.* ''' # ------------------------------ Biophysical parameters ------------------------------ alphaC_sf = 1e-5 # Calcium activation rate constant scaling factor (M) betaC = 0.1e3 # beta rate for the open-probability of iKCa channels (s-1) T = 293.15 # Room temperature (K) # ------------------------------ Gating states kinetics ------------------------------ def alpham(self, Vm): return -0.03 * (Vm + 28) / (np.exp(- (Vm + 28) / 15) - 1) * 1e3 # s-1 def betam(self, Vm): return 2.7 * np.exp(-(Vm + 53) / 18) * 1e3 # s-1 def alphah(self, Vm): return 0.045 * np.exp(-(Vm + 58) / 18) * 1e3 # s-1 def betah(self, Vm): ''' .. warning:: the original paper contains an error (multiplication) in the expression of this rate constant, corrected in the mod file on ModelDB (division). ''' return 0.72 / (np.exp(-(Vm + 23) / 14) + 1) * 1e3 # s-1 def alphan(self, Vm): return -0.024 * (Vm - 17) / (np.exp(-(Vm - 17) / 8) - 1) * 1e3 # s-1 def betan(self, Vm): return 0.2 * np.exp(-(Vm + 48) / 35) * 1e3 # s-1 def alphas(self, Vm): return -1.5 * (Vm - 20) / (np.exp(-(Vm - 20) / 5) - 1) * 1e3 # s-1 def betas(self, Vm): return 1.5 * np.exp(-(Vm + 25) / 10) * 1e3 # s-1 def alphaC(self, Cai): return 0.1 * Cai / self.alphaC_sf * 1e3 # s-1 # ------------------------------ States derivatives ------------------------------ def derC(self, c, Cai): ''' Evolution of the c-gate open-probability ''' return self.alphaC(Cai) * (1 - c) - self.betaC * c # s-1 def derStates(self): return { 'm': lambda Vm, x: self.alpham(Vm) * (1 - x['m']) - self.betam(Vm) * x['m'], 'h': lambda Vm, x: self.alphah(Vm) * (1 - x['h']) - self.betah(Vm) * x['h'], 'n': lambda Vm, x: self.alphan(Vm) * (1 - x['n']) - self.betan(Vm) * x['n'], 's': lambda Vm, x: self.alphas(Vm) * (1 - x['s']) - self.betas(Vm) * x['s'], 'c': lambda Vm, x: self.derC(x['c'], x['Cai']) } # ------------------------------ Steady states ------------------------------ def steadyStates(self): return { 'm': lambda Vm: self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)), 'h': lambda Vm: self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)), 'n': lambda Vm: self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)), 's': lambda Vm: self.alphas(Vm) / (self.alphas(Vm) + self.betas(Vm)), } # ------------------------------ Membrane currents ------------------------------ def iNa(self, m, h, Vm, Nai): ''' Sodium current ''' ENa = self.nernst(Z_Na, Nai, self.Nao, self.T) # mV return self.gNabar * m**4 * h * (Vm - ENa) # mA/m2 def iKd(self, n, Vm): ''' Delayed-rectifier Potassium current ''' return self.gKdbar * n**2 * (Vm - self.EK) # mA/m2 def iCa(self, s, Vm, Cai): ''' Calcium current ''' ECa = self.nernst(Z_Ca, Cai, self.Cao, self.T) # mV return self.gCabar * s * (Vm - ECa) # mA/m2 def iKCa(self, c, Vm): ''' Calcium-activated Potassium current ''' return self.gKCabar * c * (Vm - self.EK) # mA/m2 def iLeak(self, Vm): ''' Non-specific leakage current ''' return self.gLeak * (Vm - self.ELeak) # mA/m2 def currents(self): return { 'iNa': lambda Vm, x: self.iNa(x['m'], x['h'], Vm, x['Nai']), 'iKd': lambda Vm, x: self.iKd(x['n'], Vm), 'iCa': lambda Vm, x: self.iCa(x['s'], Vm, x['Cai']), 'iKCa': lambda Vm, x: self.iKCa(x['c'], Vm), 'iLeak': lambda Vm, _: self.iLeak(Vm) } class LeechPressure(LeechMech): ''' Leech pressure sensory neuron 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.* ''' # Neuron name name = 'LeechP' # ------------------------------ Biophysical parameters ------------------------------ # Resting parameters Cm0 = 1e-2 # Membrane capacitance (F/m2) Vm0 = -48.865 # Membrane potential (mV) Nai0 = 0.01 # Intracellular Sodium concentration (M) Cai0 = 1e-7 # Intracellular Calcium concentration (M) # Reversal potentials (mV) # ENa = 60 # Sodium (from MOD file on ModelDB) # ECa = 125 # Calcium (from MOD file on ModelDB) EK = -68.0 # Potassium ELeak = -49.0 # Non-specific leakage # Maximal channel conductances (S/m2) gNabar = 3500.0 # Sodium gKdbar = 60.0 # Delayed-rectifier Potassium gCabar = 0.02 # Calcium gKCabar = 8.0 # Calcium-dependent Potassium gLeak = 5.0 # Non-specific leakage # Ionic concentrations (M) Nao = 0.11 # Extracellular Sodium Cao = 1.8e-3 # Extracellular Calcium # Additional parameters 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) diam = 50e-6 # Cell soma diameter (m) # ------------------------------ States names & descriptions ------------------------------ states = { 'm': 'iNa activation gate', 'h': 'iNa inactivation gate', 'n': 'iKd gate', 's': 'iCa gate', 'c': 'iKCa gate', 'Nai': 'submembrane Na+ concentration (M)', 'Cai': 'submembrane Ca2+ concentration (M)' } def __init__(self): ''' Constructor of the class. ''' super().__init__() - self.rates = self.getRatesNames(['m', 'h', 'n', 's']) # Surface to volume ratio of the (spherical) cell soma (m-1) SV_ratio = 6 / self.diam # Conversion constants from membrane ionic currents into # change rate of intracellular ionic concentrations (M/s) self.K_Na = SV_ratio / (Z_Na * FARADAY) * 1e-6 # Sodium self.K_Ca = SV_ratio / (Z_Ca * FARADAY) * 1e-6 # Calcium # ------------------------------ States derivatives ------------------------------ def derStates(self): return {**super().derStates(), **{ 'Nai': lambda Vm, x: -(self.iNa(x['m'], x['h'], Vm, x['Nai']) + self.iPumpNa(x['Nai'])) * self.K_Na, 'Cai': lambda Vm, x: -(self.iCa(x['s'], Vm, x['Cai']) + self.iPumpCa(x['Cai'])) * self.K_Ca }} # ------------------------------ Steady states ------------------------------ def cinf(self, Cai): return self.alphaC(Cai) / (self.alphaC(Cai) + self.betaC) def steadyStates(self): return {**super().steadyStates(), **{ 'Nai': lambda _: self.Nai0, 'Cai': lambda _: self.Cai0, 'c': lambda _: self.cinf(self.Cai0) }} # def quasiSteadyStates(self, lkp): # qsstates = self.qsStates(lkp, ['m', 'h', 'n', 's']) # qsstates.update({ # 'Nai': self.Nai0, # 'Cai': self.Cai0 # }) # qsstates['c'] = self.alphaC(qsstates['Cai']) / (self.alphaC(qsstates['Cai']) + self.betaC) # return qsstates # ------------------------------ Membrane currents ------------------------------ def iPumpNa(self, Nai): ''' NaK-ATPase pump current ''' return self.INaPmax / (1 + np.exp((self.khalf_Na - Nai) / self.ksteep_Na)) # mA/m2 def iPumpCa(self, Cai): ''' Calcium pump current ''' return self.iCaS * (Cai - self.Cai0) / 1.5 # mA/m2 def currents(self): return {**super().currents(), **{ 'iPumpNa': lambda Vm, x: self.iPumpNa(x['Nai']) / 3., 'iPumpCa': lambda Vm, x: self.iPumpCa(x['Cai']) }} class LeechRetzius(LeechMech): ''' Leech Retzius neuron 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* iA current reference: *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.* ''' # Neuron name # name = 'LeechR' # ------------------------------ Biophysical parameters ------------------------------ # Resting parameters Cm0 = 5e-2 # Membrane capacitance (F/m2) Vm0 = -44.45 # Membrane resting potential (mV) # Reversal potentials (mV) ENa = 50.0 # Sodium (from retztemp.ses file on ModelDB) EK = -79.0 # Potassium (from retztemp.ses file on ModelDB) ECa = 125.0 # Calcium (from cachdend.mod file on ModelDB) ELeak = -30.0 # Non-specific leakage (from leakdend.mod file on ModelDB) # Maximal channel conductances (S/m2) gNabar = 1250.0 # Sodium current gKdbar = 10.0 # Delayed-rectifier Potassium GAMax = 100.0 # Transient Potassium gCabar = 4.0 # Calcium current gKCabar = 130.0 # Calcium-dependent Potassium gLeak = 1.25 # Non-specific leakage # Ionic concentrations (M) Cai = 5e-8 # Intracellular Calcium (from retztemp.ses file) # Additional parameters Vhalf = -73.1 # half-activation voltage (mV) # ------------------------------ States names & descriptions ------------------------------ states = { 'm': 'iNa activation gate', 'h': 'iNa inactivation gate', 'n': 'iKd gate', 's': 'iCa gate', 'c': 'iKCa gate', 'a': 'iA activation gate', 'b': 'iA inactivation gate', } # ------------------------------ Gating states kinetics ------------------------------ def ainf(self, Vm): 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): x = -1.5 * (Vm - self.Vhalf) * 1e-3 * FARADAY / (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): return 1. / (1 + np.exp((self.Vhalf - Vm) / -6.3)) def taub(self, Vm): x = 2 * (Vm - self.Vhalf) * 1e-3 * FARADAY / (Rg * self.T) # [-] alpha = np.exp(x) # ms-1 beta = np.exp(0.65 * x) # ms-1 return max(7.5, beta / (0.02 * (1 + alpha))) * 1e-3 # s # ------------------------------ States derivatives ------------------------------ def derStates(self, Vm, states): return {**super().derStates(Vm, states), **{ 'a': lambda Vm, x: (self.ainf(Vm) - x['a']) / self.taua(Vm), 'b': lambda Vm, x: (self.binf(Vm) - x['b']) / self.taub(Vm) }} # ------------------------------ Steady states ------------------------------ def steadyStates(self): return {**super().steadyStates(), **{ 'a': lambda Vm: self.ainf(Vm), 'b': lambda Vm: self.binf(Vm) }} # def quasiSteadyStates(self, lkp): # qsstates = self.qsStates(lkp, ['m', 'h', 'n', 's', 'a', 'b']) # qsstates['c'] = self.alphaC(self.Cai) / (self.alphaC(self.Cai) + self.betaC), # return qsstates # ------------------------------ Membrane currents ------------------------------ def iA(self, a, b, Vm): ''' Transient Potassium current ''' return self.GAMax * a * b * (Vm - self.EK) # mA/m2 def currents(self): return {**super().currents(), **{ 'iA': lambda Vm, x: self.iA(x['a'], x['b'], Vm) }} diff --git a/PySONIC/neurons/stn.py b/PySONIC/neurons/stn.py index 38dc510..8dc6f77 100644 --- a/PySONIC/neurons/stn.py +++ b/PySONIC/neurons/stn.py @@ -1,424 +1,423 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2018-11-29 16:56:45 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-25 17:27:07 +# @Last Modified time: 2019-06-25 18:09:48 import numpy as np from scipy.optimize import brentq from ..core import PointNeuron from ..constants import FARADAY, Z_Ca class OtsukaSTN(PointNeuron): ''' Sub-thalamic nucleus neuron References: *Otsuka, T., Abe, T., Tsukagawa, T., and Song, W.-J. (2004). Conductance-Based Model of the Voltage-Dependent Generation of a Plateau Potential in Subthalamic Neurons. Journal of Neurophysiology 92, 255–264.* *Tarnaud, T., Joseph, W., Martens, L., and Tanghe, E. (2018). Computational Modeling of Ultrasonic Subthalamic Nucleus Stimulation. IEEE Trans Biomed Eng.* ''' # Neuron name name = 'STN' # ------------------------------ Biophysical parameters ------------------------------ # Resting parameters Cm0 = 1e-2 # Membrane capacitance (F/m2) Vm0 = -58.0 # Membrane potential (mV) Cai0 = 5e-9 # Intracellular Calcium concentration (M) # Reversal potentials (mV) ENa = 60.0 # Sodium EK = -90.0 # Potassium ELeak = -60.0 # Non-specific leakage # Maximal channel conductances (S/m2) gNabar = 490.0 # Sodium gLeak = 3.5 # Non-specific leakage gKdbar = 570.0 # Delayed-rectifier Potassium gCaTbar = 50.0 # Low-threshold Calcium gCaLbar = 150.0 # High-threshold Calcium gAbar = 50.0 # A-type Potassium gKCabar = 10.0 # Calcium-dependent Potassium # Physical constants T = 306.15 # K (33°C) # Calcium dynamics Cao = 2e-3 # extracellular Calcium concentration (M) taur_Cai = 0.5e-3 # decay time constant for intracellular Ca2+ dissolution (s) # Fast Na current m-gate thetax_m = -40 # mV kx_m = -8 # mV tau0_m = 0.2e-3 # s tau1_m = 3e-3 # s thetaT_m = -53 # mV sigmaT_m = -0.7 # mV # Fast Na current h-gate thetax_h = -45.5 # mV kx_h = 6.4 # mV tau0_h = 0e-3 # s tau1_h = 24.5e-3 # s thetaT1_h = -50 # mV thetaT2_h = -50 # mV sigmaT1_h = -15 # mV sigmaT2_h = 16 # mV # Delayed rectifier K+ current n-gate thetax_n = -41 # mV kx_n = -14 # mV tau0_n = 0e-3 # s tau1_n = 11e-3 # s thetaT1_n = -40 # mV thetaT2_n = -40 # mV sigmaT1_n = -40 # mV sigmaT2_n = 50 # mV # T-type Ca2+ current p-gate thetax_p = -56 # mV kx_p = -6.7 # mV tau0_p = 5e-3 # s tau1_p = 0.33e-3 # s thetaT1_p = -27 # mV thetaT2_p = -102 # mV sigmaT1_p = -10 # mV sigmaT2_p = 15 # mV # T-type Ca2+ current q-gate thetax_q = -85 # mV kx_q = 5.8 # mV tau0_q = 0e-3 # s tau1_q = 400e-3 # s thetaT1_q = -50 # mV thetaT2_q = -50 # mV sigmaT1_q = -15 # mV sigmaT2_q = 16 # mV # L-type Ca2+ current c-gate thetax_c = -30.6 # mV kx_c = -5 # mV tau0_c = 45e-3 # s tau1_c = 10e-3 # s thetaT1_c = -27 # mV thetaT2_c = -50 # mV sigmaT1_c = -20 # mV sigmaT2_c = 15 # mV # L-type Ca2+ current d1-gate thetax_d1 = -60 # mV kx_d1 = 7.5 # mV tau0_d1 = 400e-3 # s tau1_d1 = 500e-3 # s thetaT1_d1 = -40 # mV thetaT2_d1 = -20 # mV sigmaT1_d1 = -15 # mV sigmaT2_d1 = 20 # mV # L-type Ca2+ current d2-gate thetax_d2 = 0.1e-6 # M kx_d2 = 0.02e-6 # M tau_d2 = 130e-3 # s # A-type K+ current a-gate thetax_a = -45 # mV kx_a = -14.7 # mV tau0_a = 1e-3 # s tau1_a = 1e-3 # s thetaT_a = -40 # mV sigmaT_a = -0.5 # mV # A-type K+ current b-gate thetax_b = -90 # mV kx_b = 7.5 # mV tau0_b = 0e-3 # s tau1_b = 200e-3 # s thetaT1_b = -60 # mV thetaT2_b = -40 # mV sigmaT1_b = -30 # mV sigmaT2_b = 10 # mV # Ca2+-activated K+ current r-gate thetax_r = 0.17e-6 # M kx_r = -0.08e-6 # M tau_r = 2e-3 # s # ------------------------------ States names & descriptions ------------------------------ states = { 'm': 'iNa activation gate', 'h': 'iNa inactivation gate', 'n': 'iKd gate', 'a': 'iA activation gate', 'b': 'iA inactivation gate', 'p': 'iCaT activation gate', 'q': 'iCaT inactivation gate', 'c': 'iCaL activation gate', 'd1': 'iCaL inactivation gate 1', 'd2': 'iCaL inactivation gate 2', 'r': 'iCaK gate', 'Cai': 'submembrane Calcium concentration (M)' } def __init__(self): super().__init__() - self.rates = self.getRatesNames(['a', 'b', 'c', 'd1', 'm', 'h', 'n', 'p', 'q']) self.deff = self.getEffectiveDepth(self.Cai0, self.Vm0) # m self.iCa_to_Cai_rate = self.currentToConcentrationRate(Z_Ca, self.deff) # Mmol.m-1.C-1 def getPltScheme(self): pltscheme = super().getPltScheme() pltscheme['[Ca^{2+}]_i'] = ['Cai'] return pltscheme def getPltVars(self, wrapleft='df["', wrapright='"]'): pltvars = super().getPltVars(wrapleft, wrapright) pltvars['Cai'] = { 'desc': 'submembrane Ca2+ concentration', 'label': '[Ca^{2+}]_i', 'unit': 'uM', 'factor': 1e6 } return pltvars def titrationFunc(self, *args, **kwargs): return self.isSilenced(*args, **kwargs) def getEffectiveDepth(self, Cai, Vm): ''' Compute effective depth that matches a given membrane potential and intracellular Calcium concentration. :return: effective depth (m) ''' iCaT = self.iCaT(self.pinf(Vm), self.qinf(Vm), Vm, Cai) # mA/m2 iCaL = self.iCaL(self.cinf(Vm), self.d1inf(Vm), self.d2inf(Cai), Vm, Cai) # mA/m2 return -(iCaT + iCaL) / (Z_Ca * FARADAY * Cai / self.taur_Cai) * 1e-6 # m # ------------------------------ Gating states kinetics ------------------------------ def _xinf(self, var, theta, k): ''' Generic function computing the steady-state opening of a particular channel gate 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 opening (-) ''' 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, tau0, tau1): ''' 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) # ------------------------------ States derivatives ------------------------------ def derCai(self, p, q, c, d1, d2, Cai, Vm): iCaT = self.iCaT(p, q, Vm, Cai) iCaL = self.iCaL(c, d1, d2, Vm, Cai) return - self.iCa_to_Cai_rate * (iCaT + iCaL) - Cai / self.taur_Cai # M/s def derStates(self): return { 'a': lambda Vm, x: (self.ainf(Vm) - x['a']) / self.taua(Vm), 'b': lambda Vm, x: (self.binf(Vm) - x['b']) / self.taub(Vm), 'c': lambda Vm, x: (self.cinf(Vm) - x['c']) / self.tauc(Vm), 'd1': lambda Vm, x: (self.d1inf(Vm) - x['d1']) / self.taud1(Vm), 'd2': lambda Vm, x: (self.d2inf(x['Cai']) - x['d2']) / self.tau_d2, 'm': lambda Vm, x: (self.minf(Vm) - x['m']) / self.taum(Vm), 'h': lambda Vm, x: (self.hinf(Vm) - x['h']) / self.tauh(Vm), 'n': lambda Vm, x: (self.ninf(Vm) - x['n']) / self.taun(Vm), 'p': lambda Vm, x: (self.pinf(Vm) - x['p']) / self.taup(Vm), 'q': lambda Vm, x: (self.qinf(Vm) - x['q']) / self.tauq(Vm), 'r': lambda Vm, x: (self.rinf(x['Cai']) - x['r']) / self.tau_r, 'Cai': lambda Vm, x: self.derCai(x['p'], x['q'], x['c'], x['d1'], x['d2'], x['Cai'], Vm) } # ------------------------------ Steady states ------------------------------ def Caiinf(self, Vm): ''' Steady-state intracellular Calcium concentration ''' if isinstance(Vm, np.ndarray): return np.array([self.Caiinf(v) for v in Vm]) else: return brentq( lambda x: self.derCai(self.pinf(Vm), self.qinf(Vm), self.cinf(Vm), self.d1inf(Vm), self.d2inf(x), x, Vm), self.Cai0 * 1e-4, self.Cai0 * 1e3, xtol=1e-16 ) def steadyStates(self): return { 'a': lambda Vm: self.ainf(Vm), 'b': lambda Vm: self.binf(Vm), 'c': lambda Vm: self.cinf(Vm), 'd1': lambda Vm: self.d1inf(Vm), 'd2': lambda Vm: self.d2inf(self.Caiinf(Vm)), 'm': lambda Vm: self.minf(Vm), 'h': lambda Vm: self.hinf(Vm), 'n': lambda Vm: self.ninf(Vm), 'p': lambda Vm: self.pinf(Vm), 'q': lambda Vm: self.qinf(Vm), 'r': lambda Vm: self.rinf(self.Caiinf(Vm)), 'Cai': lambda Vm: self.Caiinf(Vm) } # def quasiSteadyStates(self, lkp): # qsstates = self.qsStates(lkp, ['a', 'b', 'c', 'd1', 'm', 'h', 'n', 'p', 'q']) # qsstates['Cai'] = self.Caiinf(lkp['V'], qsstates['p'], qsstates['q'], qsstates['c'], # qsstates['d1']) # qsstates['d2'] = self.d2inf(qsstates['Cai']) # qsstates['r'] = self.rinf(qsstates['Cai']) # return qsstates # ------------------------------ Membrane currents ------------------------------ def iNa(self, m, h, Vm): ''' Sodium current ''' return self.gNabar * m**3 * h * (Vm - self.ENa) # mA/m2 def iKd(self, n, Vm): ''' delayed-rectifier Potassium current ''' return self.gKdbar * n**4 * (Vm - self.EK) # mA/m2 def iA(self, a, b, Vm): ''' A-type Potassium current ''' return self.gAbar * a**2 * b * (Vm - self.EK) # mA/m2 def iCaT(self, p, q, Vm, Cai): ''' low-threshold (T-type) Calcium current ''' return self.gCaTbar * p**2 * q * (Vm - self.nernst(Z_Ca, Cai, self.Cao, self.T)) # mA/m2 def iCaL(self, c, d1, d2, Vm, Cai): ''' high-threshold (L-type) Calcium current ''' return self.gCaLbar * c**2 * d1 * d2 * ( Vm - self.nernst(Z_Ca, Cai, self.Cao, self.T)) # mA/m2 def iKCa(self, r, Vm): ''' Calcium-activated Potassium current ''' return self.gKCabar * r**2 * (Vm - self.EK) # mA/m2 def iLeak(self, Vm): ''' non-specific leakage current ''' return self.gLeak * (Vm - self.ELeak) # mA/m2 def currents(self): return { 'iNa': lambda Vm, x: self.iNa(x['m'], x['h'], Vm), 'iKd': lambda Vm, x: self.iKd(x['n'], Vm), 'iA': lambda Vm, x: self.iA(x['a'], x['b'], Vm), 'iCaT': lambda Vm, x: self.iCaT(x['p'], x['q'], Vm, x['Cai']), 'iCaL': lambda Vm, x: self.iCaL( x['c'], x['d1'], x['d2'], Vm, x['Cai']), 'iKCa': lambda Vm, x: self.iKCa(x['r'], Vm), 'iLeak': lambda Vm, _: self.iLeak(Vm) } # ------------------------------ Other methods ------------------------------ def getLowIntensities(self): ''' Return an array of acoustic intensities (W/m2) used to study the STN neuron in Tarnaud, T., Joseph, W., Martens, L., and Tanghe, E. (2018). Computational Modeling of Ultrasonic Subthalamic Nucleus Stimulation. IEEE Trans Biomed Eng. ''' return np.hstack(( np.arange(10, 101, 10), np.arange(101, 131, 1), np.array([140]) )) # W/m2 diff --git a/PySONIC/neurons/thalamic.py b/PySONIC/neurons/thalamic.py index f1961ea..3b51caf 100644 --- a/PySONIC/neurons/thalamic.py +++ b/PySONIC/neurons/thalamic.py @@ -1,374 +1,372 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2017-07-31 15:20:54 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-25 17:36:35 +# @Last Modified time: 2019-06-25 18:07:25 import numpy as np from ..core import PointNeuron from ..constants import Z_Ca class Thalamic(PointNeuron): ''' Generic thalamic neuron 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.* ''' # ------------------------------ Biophysical parameters ------------------------------ # Resting parameters Cm0 = 1e-2 # Membrane capacitance (F/m2) # Reversal potentials (mV) ENa = 50.0 # Sodium EK = -90.0 # Potassium ECa = 120.0 # Calcium # ------------------------------ Gating states kinetics ------------------------------ def alpham(self, Vm): return 0.32 * self.vtrap(13 - (Vm - self.VT), 4) * 1e3 # s-1 def betam(self, Vm): return 0.28 * self.vtrap((Vm - self.VT) - 40, 5) * 1e3 # s-1 def alphah(self, Vm): return 0.128 * np.exp(-((Vm - self.VT) - 17) / 18) * 1e3 # s-1 def betah(self, Vm): return 4 / (1 + np.exp(-((Vm - self.VT) - 40) / 5)) * 1e3 # s-1 def alphan(self, Vm): return 0.032 * self.vtrap(15 - (Vm - self.VT), 5) * 1e3 # s-1 def betan(self, Vm): return 0.5 * np.exp(-((Vm - self.VT) - 10) / 40) * 1e3 # s-1 # ------------------------------ States derivatives ------------------------------ def derStates(self): return { 'm': lambda Vm, x: self.alpham(Vm) * (1 - x['m']) - self.betam(Vm) * x['m'], 'h': lambda Vm, x: self.alphah(Vm) * (1 - x['h']) - self.betah(Vm) * x['h'], 'n': lambda Vm, x: self.alphan(Vm) * (1 - x['n']) - self.betan(Vm) * x['n'], 's': lambda Vm, x: (self.sinf(Vm) - x['s']) / self.taus(Vm), 'u': lambda Vm, x: (self.uinf(Vm) - x['u']) / self.tauu(Vm) } # ------------------------------ Steady states ------------------------------ def steadyStates(self): return { 'm': lambda Vm: self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)), 'h': lambda Vm: self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)), 'n': lambda Vm: self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)), 's': lambda Vm: self.sinf(Vm), 'u': lambda Vm: self.uinf(Vm) } # ------------------------------ Membrane currents ------------------------------ def iNa(self, m, h, Vm): ''' Sodium current ''' return self.gNabar * m**3 * h * (Vm - self.ENa) # mA/m2 def iKd(self, n, Vm): ''' delayed-rectifier Potassium current ''' return self.gKdbar * n**4 * (Vm - self.EK) def iCaT(self, s, u, Vm): ''' low-threshold (Ts-type) Calcium current ''' return self.gCaTbar * s**2 * u * (Vm - self.ECa) # mA/m2 def iLeak(self, Vm): ''' non-specific leakage current ''' return self.gLeak * (Vm - self.ELeak) # mA/m2 def currents(self): return { 'iNa': lambda Vm, states: self.iNa(states['m'], states['h'], Vm), 'iKd': lambda Vm, states: self.iKd(states['n'], Vm), 'iCaT': lambda Vm, states: self.iCaT(states['s'], states['u'], Vm), 'iLeak': lambda Vm, _: self.iLeak(Vm) } # ------------------------------ Other methods ------------------------------ def computeEffRates(self, Vm): return { 'alpham': np.mean(self.alpham(Vm)), 'betam': np.mean(self.betam(Vm)), 'alphah': np.mean(self.alphah(Vm)), 'betah': np.mean(self.betah(Vm)), 'alphan': np.mean(self.alphan(Vm)), 'betan': np.mean(self.betan(Vm)), 'alphas': np.mean(self.sinf(Vm) / self.taus(Vm)), 'betas': np.mean((1 - self.sinf(Vm)) / self.taus(Vm)), 'alphau': np.mean(self.uinf(Vm) / self.tauu(Vm)), 'betau': np.mean((1 - self.uinf(Vm)) / self.tauu(Vm)) } class ThalamicRE(Thalamic): ''' 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.* ''' # Neuron name name = 'RE' # ------------------------------ Biophysical parameters ------------------------------ # Resting parameters Vm0 = -89.5 # Membrane potential (mV) # Reversal potentials (mV) ELeak = -90.0 # Non-specific leakage # Maximal channel conductances (S/m2) gNabar = 2000.0 # Sodium gKdbar = 200.0 # Delayed-rectifier Potassium gCaTbar = 30.0 # Low-threshold Calcium gLeak = 0.5 # Non-specific leakage # Additional parameters VT = -67.0 # Spike threshold adjustment parameter (mV) # ------------------------------ States names & descriptions ------------------------------ states = { 'm': 'iNa activation gate', 'h': 'iNa inactivation gate', 'n': 'iKd gate', 's': 'iCaT activation gate', 'u': 'iCaT inactivation gate' } # ------------------------------ Gating states kinetics ------------------------------ def sinf(self, Vm): return 1.0 / (1.0 + np.exp(-(Vm + 52.0) / 7.4)) def taus(self, Vm): 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): return 1.0 / (1.0 + np.exp((Vm + 80.0) / 5.0)) def tauu(self, Vm): 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): ''' Thalamo-cortical neuron 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.* Model of Ca2+ buffering and contribution from iCaT 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.* ''' # Neuron name name = 'TC' # ------------------------------ Biophysical parameters ------------------------------ # Resting parameters # Vm0 = -63.4 # Membrane potential (mV) Vm0 = -61.93 # Membrane potential (mV) # Reversal potentials (mV) EH = -40.0 # Mixed cationic current ELeak = -70.0 # Non-specific leakage # Maximal channel conductances (S/m2) gNabar = 900.0 # Sodium gKdbar = 100.0 # Delayed-rectifier Potassium gCaTbar = 20.0 # Low-threshold Calcium gKLeak = 0.138 # Leakage Potassium gHbar = 0.175 # Mixed cationic current gLeak = 0.1 # Non-specific leakage # Additional parameters VT = -52.0 # Spike threshold adjustment parameter (mV) Vx = 0.0 # Voltage-dependence uniform shift factor at 36°C (mV) taur_Cai = 5e-3 # decay time constant for intracellular Ca2+ dissolution (s) Cai_min = 50e-9 # minimal intracellular Calcium concentration (M) deff = 100e-9 # effective depth beneath membrane for intracellular [Ca2+] calculation 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) # ------------------------------ States names & descriptions ------------------------------ states = { 'm': 'iNa activation gate', 'h': 'iNa inactivation gate', 'n': 'iKd gate', 's': 'iCaT activation gate', 'u': 'iCaT inactivation gate', 'C': 'iH gate closed state', 'O': 'iH gate open state', 'Cai': 'submembrane Ca2+ concentration (M)', 'P0': 'proportion of unbound iH regulating factor', } def __init__(self): super().__init__() - self.rates = self.getRatesNames(['m', 'h', 'n', 's', 'u', 'O']) self.iCa_to_Cai_rate = self.currentToConcentrationRate(Z_Ca, self.deff) - # self.states += ['O', 'C', 'P0', 'Cai'] def OL(self, O, C): ''' O-gate locked-open probability ''' return 1 - O - C def getPltScheme(self): pltscheme = super().getPltScheme() pltscheme['i_{H}\\ kin.'] = ['O', 'OL', 'P0'] pltscheme['[Ca^{2+}]_i'] = ['Cai'] return pltscheme def getPltVars(self, wrapleft='df["', wrapright='"]'): return {**super().getPltVars(wrapleft, wrapright), **{ 'Cai': { 'desc': 'sumbmembrane Ca2+ concentration', 'label': '[Ca^{2+}]_i', 'unit': 'uM', 'factor': 1e6 }, 'OL': { 'desc': 'iH O-gate locked-opening', 'label': 'O_L', 'bounds': (-0.1, 1.1), 'func': 'OL({0}O{1}, {0}C{1})'.format(wrapleft, wrapright) }, 'P0': { 'desc': 'iH regulating factor activation', 'label': 'P_0', 'bounds': (-0.1, 1.1) } }} # ------------------------------ Gating states kinetics ------------------------------ def sinf(self, Vm): return 1.0 / (1.0 + np.exp(-(Vm + self.Vx + 57.0) / 6.2)) def taus(self, Vm): x = 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 / x) * 1e-3 # s def uinf(self, Vm): return 1.0 / (1.0 + np.exp((Vm + self.Vx + 81.0) / 4.0)) def tauu(self, Vm): 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 oinf(self, Vm): return 1.0 / (1.0 + np.exp((Vm + 75.0) / 5.5)) def tauo(self, Vm): return 1 / (np.exp(-14.59 - 0.086 * Vm) + np.exp(-1.87 + 0.0701 * Vm)) * 1e-3 # s def alphao(self, Vm): return self.oinf(Vm) / self.tauo(Vm) # s-1 def betao(self, Vm): return (1 - self.oinf(Vm)) / self.tauo(Vm) # s-1 # ------------------------------ States derivatives ------------------------------ def derCai(self, Cai, s, u, Vm): return (self.Cai_min - Cai) / self.taur_Cai -\ self.iCa_to_Cai_rate * self.iCaT(s, u, Vm) # M/s def derStates(self): return {**super().derStates(), **{ 'C': lambda Vm, x: self.betao(Vm) * x['O'] - self.alphao(Vm) * x['C'], 'O': lambda Vm, x: (self.alphao(Vm) * x['C'] - self.betao(Vm) * x['O'] - self.k3 * x['O'] * (1 - x['P0']) + self.k4 * (1 - x['O'] - x['C'])), 'P0': lambda _, x: self.k2 * (1 - x['P0']) - self.k1 * x['P0'] * x['Cai']**self.nCa, 'Cai': lambda Vm, x: ((self.Cai_min - x['Cai']) / self.taur_Cai - self.iCa_to_Cai_rate * self.iCaT(x['s'], x['u'], Vm)) # M/s }} # ------------------------------ Steady states ------------------------------ def Cinf(self, Cai, Vm): ''' Steady-state O-gate closed-probability ''' return self.betao(Vm) / self.alphao(Vm) * self.Oinf(Cai, Vm) def Oinf(self, Cai, Vm): ''' Steady-state O-gate open-probability ''' return self.k4 / (self.k3 * (1 - self.P0inf(Cai)) + self.k4 * ( 1 + self.betao(Vm) / self.alphao(Vm))) def P0inf(self, Cai): ''' Steady-state unbound probability of Ih regulating factor ''' return self.k2 / (self.k2 + self.k1 * Cai**self.nCa) def Caiinf(self, Vm): ''' Steady-state intracellular Calcium concentration ''' return self.Cai_min - self.taur_Cai * self.iCa_to_Cai_rate * self.iCaT( self.sinf(Vm), self.uinf(Vm), Vm) # M def steadyStates(self): return {**super().steadyStates(), **{ 'O': lambda Vm: self.Oinf(self.Caiinf(Vm), Vm), 'C': lambda Vm: self.Cinf(self.Caiinf(Vm), Vm), 'P0': lambda Vm: self.P0inf(self.Caiinf(Vm)), 'Cai': lambda Vm: self.Caiinf(Vm) }} # def quasiSteadyStates(self, lkp): # qsstates = super().quasiSteadyStates(lkp) # qsstates.update({ # 'O': self.Oinf(self.Caiinf(lkp['V']), lkp['V']), # 'C': self.Cinf(lkp['V']), # 'P0': self.P0inf(self.Caiinf(lkp['V'])), # 'Cai': self.Caiinf(lkp['V']) # }) # return qsstates # ------------------------------ Membrane currents ------------------------------ def iKLeak(self, Vm): ''' Potassium leakage current ''' return self.gKLeak * (Vm - self.EK) # mA/m2 def iH(self, O, C, Vm): ''' outward mixed cationic current ''' return self.gHbar * (O + 2 * self.OL(O, C)) * (Vm - self.EH) # mA/m2 def currents(self): return {**super().currents(), **{ 'iKLeak': lambda Vm, states: self.iKLeak(Vm), 'iH': lambda Vm, states: self.iH(states['O'], states['C'], Vm) }}