diff --git a/PySONIC/core/pneuron.py b/PySONIC/core/pneuron.py index fe6f75d..a83262c 100644 --- a/PySONIC/core/pneuron.py +++ b/PySONIC/core/pneuron.py @@ -1,511 +1,517 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-03 11:53:04 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-03-13 14:50:53 +# @Last Modified time: 2019-03-13 15:40:06 import os import time import pickle import abc import inspect import numpy as np from scipy.integrate import odeint import pandas as pd from ..postpro import findPeaks from ..constants import * from ..utils import si_format, logger, ESTIM_filecode from ..batches import xlslog class PointNeuron(metaclass=abc.ABCMeta): ''' Abstract class defining the common API (i.e. mandatory attributes and methods) of all subclasses implementing the channels mechanisms of specific point neurons. The mandatory attributes are: - **name**: a string defining the name of the mechanism. - **Cm0**: a float defining the membrane resting capacitance (in F/m2) - **Vm0**: a float defining the membrane resting potential (in mV) - **states_names**: a list of strings defining the names of the different state probabilities governing the channels behaviour (i.e. the differential HH variables). - **states0**: a 1D array of floats (NOT integers !!!) defining the initial values of the different state probabilities. - **coeff_names**: a list of strings defining the names of the different coefficients to be used in effective simulations. The mandatory methods are: - **iNet**: compute the net ionic current density (in mA/m2) across the membrane, given a specific membrane potential (in mV) and channel states. - **steadyStates**: compute the channels steady-state values for a specific membrane potential value (in mV). - **derStates**: compute the derivatives of channel states, given a specific membrane potential (in mV) and channel states. This method must return a list of derivatives ordered identically as in the states0 attribute. - **getEffRates**: get the effective rate constants of ion channels to be used in effective simulations. This method must return an array of effective rates ordered identically as in the coeff_names attribute. - **derStatesEff**: compute the effective derivatives of channel states, based on 1-dimensional linear interpolators of "effective" coefficients. This method must return a list of derivatives ordered identically as in the states0 attribute. ''' def __repr__(self): return self.__class__.__name__ def pprint(self): return '{} neuron'.format(self.__class__.__name__) @property @abc.abstractmethod def name(self): return 'Should never reach here' @property @abc.abstractmethod def Cm0(self): return 'Should never reach here' @property @abc.abstractmethod def Vm0(self): return 'Should never reach here' @abc.abstractmethod def currents(self, Vm, states): ''' Compute all ionic currents per unit area. :param Vm: membrane potential (mV) :states: state probabilities of the ion channels :return: dictionary of ionic currents per unit area (mA/m2) ''' def iNet(self, Vm, states): ''' Net membrane current :param Vm: membrane potential (mV) :states: state probabilities of the ion channels :return: current per unit area (mA/m2) ''' return sum(self.currents(Vm, states).values()) def getCurrentsPltVars(self): ''' Return a dictionary of implemented currents with their description. ''' # Get list of all current names cnames = list(self.currents(np.nan, [np.nan] * len(self.states_names)).keys()) # Generate a dictionary of 1 plot varibale from each ionic current cpltvars = {} for cname in cnames: cfunc = getattr(self, cname) cargs = inspect.getargspec(cfunc)[0][1:] cpltvars[cname] = dict( desc=inspect.getdoc(cfunc).splitlines()[0], label='I_{{{}}}'.format(cname[1:]), unit='A/m^2', factor=1e-3, alias='neuron.{}({})'.format(cname, ', '.join(['df["{}"]'.format(a) for a in cargs])) ) # add plot variable for net membrane current cpltvars['iNet'] = dict( desc=inspect.getdoc(getattr(self, 'iNet')).splitlines()[0], label='I_{net}', unit='A/m^2', factor=1e-3, alias='neuron.iNet(df["Vm"], neuron_states)' ) return cpltvars @abc.abstractmethod def steadyStates(self, Vm): ''' Compute the channels steady-state values for a specific membrane potential value. :param Vm: membrane potential (mV) :return: array of steady-states ''' @abc.abstractmethod def derStates(self, Vm, states): ''' Compute the derivatives of channel states. :param Vm: membrane potential (mV) :states: state probabilities of the ion channels :return: current per unit area (mA/m2) ''' @abc.abstractmethod def getEffRates(self, Vm): ''' Get the effective rate constants of ion channels, averaged along an acoustic cycle, for future use in effective simulations. :param Vm: array of membrane potential values for an acoustic cycle (mV) :return: an array of rate average constants (s-1) ''' @abc.abstractmethod def derStatesEff(self, Qm, states, interp_data): ''' Compute the effective derivatives of channel states, based on 1-dimensional linear interpolation of "effective" coefficients that summarize the system's behaviour over an acoustic cycle. :param Qm: membrane charge density (C/m2) :states: state probabilities of the ion channels :param interp_data: dictionary of 1D vectors of "effective" coefficients over the charge domain, for specific frequency and amplitude values. ''' + + def 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 getGates(self): ''' Retrieve the names of the neuron's states that match an ion channel gating. ''' gates = [] for x in self.states_names: if 'alpha{}'.format(x.lower()) in self.coeff_names: gates.append(x) return gates def getRates(self, Vm): ''' Compute the ion channels rate constants for a given membrane potential. :param Vm: membrane potential (mV) :return: a dictionary of rate constants and their values at the given potential. ''' rates = {} for x in self.getGates(): x = x.lower() alpha_str, beta_str = ['{}{}'.format(s, x.lower()) for s in ['alpha', 'beta']] inf_str, tau_str = ['{}inf'.format(x.lower()), 'tau{}'.format(x.lower())] if hasattr(self, 'alpha{}'.format(x)): alphax = getattr(self, alpha_str)(Vm) betax = getattr(self, beta_str)(Vm) elif hasattr(self, '{}inf'.format(x)): xinf = getattr(self, inf_str)(Vm) taux = getattr(self, tau_str)(Vm) alphax = xinf / taux betax = 1 / taux - alphax rates[alpha_str] = alphax rates[beta_str] = betax return rates def Vderivatives(self, y, t, Iinj): ''' Compute the derivatives of a V-cast HH system for a specific value of injected current. :param y: vector of HH system variables at time t :param t: time value (s, unused) :param Iinj: injected current (mA/m2) :return: vector of HH system derivatives at time t ''' Vm, *states = y Iionic = self.iNet(Vm, states) # mA/m2 dVmdt = (- Iionic + Iinj) / self.Cm0 # mV/s dstates = self.derStates(Vm, states) return [dVmdt, *dstates] def Qderivatives(self, y, t, Cm=None): ''' Compute the derivatives of the n-ODE HH system variables, based on a value of membrane capacitance. :param y: vector of HH system variables at time t :param t: specific instant in time (s) :param Cm: membrane capacitance (F/m2) :return: vector of HH system derivatives at time t ''' if Cm is None: Cm = self.Cm0 Qm, *states = y Vm = Qm / Cm * 1e3 # mV dQm = - self.iNet(Vm, states) * 1e-3 # A/m2 dstates = self.derStates(Vm, states) return [dQm, *dstates] def checkInputs(self, Astim, tstim, toffset, PRF, DC): ''' Check validity of electrical stimulation parameters. :param Astim: pulse amplitude (mA/m2) :param tstim: pulse duration (s) :param toffset: offset duration (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) ''' # Check validity of stimulation parameters if not all(isinstance(param, float) for param in [Astim, tstim, toffset, DC]): raise TypeError('Invalid stimulation parameters (must be float typed)') if tstim <= 0: raise ValueError('Invalid stimulus duration: {} ms (must be strictly positive)' .format(tstim * 1e3)) if toffset < 0: raise ValueError('Invalid stimulus offset: {} ms (must be positive or null)' .format(toffset * 1e3)) if DC <= 0.0 or DC > 1.0: raise ValueError('Invalid duty cycle: {} (must be within ]0; 1])'.format(DC)) if DC < 1.0: if not isinstance(PRF, float): raise TypeError('Invalid PRF value (must be float typed)') if PRF is None: raise AttributeError('Missing PRF value (must be provided when DC < 1)') if PRF < 1 / tstim: raise ValueError('Invalid PRF: {} Hz (PR interval exceeds stimulus duration)' .format(PRF)) def simulate(self, Astim, tstim, toffset, PRF=None, DC=1.0): ''' Compute solutions of a neuron's HH system for a specific set of electrical stimulation parameters, using a classic integration scheme. :param Astim: pulse amplitude (mA/m2) :param tstim: pulse duration (s) :param toffset: offset duration (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) :return: 3-tuple with the time profile and solution matrix and a state vector ''' # Check validity of stimulation parameters self.checkInputs(Astim, tstim, toffset, PRF, DC) # Determine system time step dt = DT_ESTIM # if CW stimulus: divide integration during stimulus into single interval if DC == 1.0: PRF = 1 / tstim # Compute vector sizes npulses = int(np.round(PRF * tstim)) Tpulse_on = DC / PRF Tpulse_off = (1 - DC) / PRF # For high-PRF pulsed protocols: adapt time step to ensure minimal # number of samples during TON or TOFF dt_warning_msg = 'high-PRF protocol: lowering time step to %.2e s to properly integrate %s' for key, Tpulse in {'TON': Tpulse_on, 'TOFF': Tpulse_off}.items(): if Tpulse > 0 and Tpulse / dt < MIN_SAMPLES_PER_PULSE_INT: dt = Tpulse / MIN_SAMPLES_PER_PULSE_INT logger.warning(dt_warning_msg, dt, key) n_pulse_on = int(np.round(Tpulse_on / dt)) n_pulse_off = int(np.round(Tpulse_off / dt)) # Compute offset size n_off = int(np.round(toffset / dt)) # Set initial conditions y0 = [self.Vm0, *self.states0] nvar = len(y0) # Initialize global arrays t = np.array([0.]) states = np.array([1]) y = np.array([y0]).T # Initialize pulse time and states vectors t_pulse0 = np.linspace(0, Tpulse_on + Tpulse_off, n_pulse_on + n_pulse_off) states_pulse = np.concatenate((np.ones(n_pulse_on), np.zeros(n_pulse_off))) # Loop through all pulse (ON and OFF) intervals for i in range(npulses): # Construct and initialize arrays t_pulse = t_pulse0 + t[-1] y_pulse = np.empty((nvar, n_pulse_on + n_pulse_off)) # Integrate ON system y_pulse[:, :n_pulse_on] = odeint( self.Vderivatives, y[:, -1], t_pulse[:n_pulse_on], args=(Astim,)).T # Integrate OFF system if n_pulse_off > 0: y_pulse[:, n_pulse_on:] = odeint( self.Vderivatives, y_pulse[:, n_pulse_on - 1], t_pulse[n_pulse_on:], args=(0.0,)).T # Append pulse arrays to global arrays states = np.concatenate([states, states_pulse[1:]]) t = np.concatenate([t, t_pulse[1:]]) y = np.concatenate([y, y_pulse[:, 1:]], axis=1) # Integrate offset interval if n_off > 0: t_off = np.linspace(0, toffset, n_off) + t[-1] states_off = np.zeros(n_off) y_off = odeint(self.Vderivatives, y[:, -1], t_off, args=(0.0, )).T # Concatenate offset arrays to global arrays states = np.concatenate([states, states_off[1:]]) t = np.concatenate([t, t_off[1:]]) y = np.concatenate([y, y_off[:, 1:]], axis=1) # Return output variables return (t, y, states) def titrate(self, tstim, toffset, PRF=None, DC=1.0, Arange=(0., 2 * TITRATION_ESTIM_A_MAX)): ''' Use a dichotomic recursive search to determine the threshold amplitude needed to obtain neural excitation for a given duration, PRF and duty cycle. :param tstim: duration of US stimulation (s) :param toffset: duration of the offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) :param Arange: search interval for Astim, iteratively refined :return: 5-tuple with the determined threshold, time profile, solution matrix, state vector and response latency ''' Astim = (Arange[0] + Arange[1]) / 2 # Run simulation and detect spikes t0 = time.time() (t, y, states) = self.simulate(Astim, tstim, toffset, PRF, DC) tcomp = time.time() - t0 dt = t[1] - t[0] ipeaks, *_ = findPeaks(y[0, :], SPIKE_MIN_VAMP, int(np.ceil(SPIKE_MIN_DT / dt)), SPIKE_MIN_VPROM) nspikes = ipeaks.size latency = t[ipeaks[0]] if nspikes > 0 else None logger.debug('A = %sA/m2 ---> %s spike%s detected', si_format(Astim * 1e-3, 2, space=' '), nspikes, "s" if nspikes > 1 else "") # If accurate threshold is found, return simulation results if (Arange[1] - Arange[0]) <= TITRATION_ESTIM_DA_MAX and nspikes == 1: return (Astim, t, y, states, latency, tcomp) # Otherwise, refine titration interval and iterate recursively else: if nspikes == 0: # if Astim too close to max then stop if (TITRATION_ESTIM_A_MAX - Astim) <= TITRATION_ESTIM_DA_MAX: return (np.nan, t, y, states, latency, tcomp) Arange = (Astim, Arange[1]) else: Arange = (Arange[0], Astim) return self.titrate(tstim, toffset, PRF, DC, Arange=Arange) def runAndSave(self, outdir, tstim, toffset, PRF=None, DC=1.0, Astim=None): ''' Run a simulation of the point-neuron Hodgkin-Huxley system with specific parameters, and save the results in a PKL file. :param outdir: full path to output directory :param tstim: stimulus duration (s) :param toffset: stimulus offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: stimulus duty cycle (-) :param Astim: stimulus amplitude (mA/m2) ''' # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") if Astim is not None: logger.info('%s: simulation @ A = %sA/m2, t = %ss (%ss offset)%s', self, si_format(Astim * 1e-3, 2, space=' '), *si_format([tstim, toffset], 1, space=' '), (', PRF = {}Hz, DC = {:.2f}%'.format(si_format(PRF, 2, space=' '), DC * 1e2) if DC < 1.0 else '')) # Run simulation tstart = time.time() t, y, states = self.simulate(Astim, tstim, toffset, PRF, DC) Vm, *channels = y tcomp = time.time() - tstart # Detect spikes on Vm signal dt = t[1] - t[0] ipeaks, *_ = findPeaks(Vm, SPIKE_MIN_VAMP, int(np.ceil(SPIKE_MIN_DT / dt)), SPIKE_MIN_VPROM) nspikes = ipeaks.size lat = t[ipeaks[0]] if nspikes > 0 else 'N/A' outstr = '{} spike{} detected'.format(nspikes, 's' if nspikes > 1 else '') else: logger.info('%s: titration @ t = %ss%s', self, si_format(tstim, 1, space=' '), (', PRF = {}Hz, DC = {:.2f}%'.format(si_format(PRF, 2, space=' '), DC * 1e2) if DC < 1.0 else '')) # Run titration Astim, t, y, states, lat, tcomp = self.titrate(tstim, toffset, PRF, DC) Vm, *channels = y nspikes = 1 if Astim is np.nan: outstr = 'no spikes detected within titration interval' nspikes = 0 else: nspikes = 1 outstr = 'Athr = {}A/m2'.format(si_format(Astim * 1e-3, 2, space=' ')) logger.debug('completed in %s, %s', si_format(tcomp, 1), outstr) sr = np.mean(1 / np.diff(t[ipeaks])) if nspikes > 1 else None # Store dataframe and metadata df = pd.DataFrame({ 't': t, 'states': states, 'Vm': Vm, 'Qm': Vm * self.Cm0 * 1e-3 }) for j in range(len(self.states_names)): df[self.states_names[j]] = channels[j] meta = { 'neuron': self.name, 'Astim': Astim, 'tstim': tstim, 'toffset': toffset, 'PRF': PRF, 'DC': DC, 'tcomp': tcomp } # Export into to PKL file simcode = ESTIM_filecode(self.name, Astim, tstim, PRF, DC) outpath = '{}/{}.pkl'.format(outdir, simcode) with open(outpath, 'wb') as fh: pickle.dump({'meta': meta, 'data': df}, fh) logger.debug('simulation data exported to "%s"', outpath) # Export key metrics to log file logpath = os.path.join(outdir, 'log_ESTIM.xlsx') logentry = { 'Date': date_str, 'Time': daytime_str, 'Neuron Type': self.name, 'Astim (mA/m2)': Astim, 'Tstim (ms)': tstim * 1e3, 'PRF (kHz)': PRF * 1e-3 if DC < 1 else 'N/A', 'Duty factor': DC, '# samples': t.size, 'Comp. time (s)': round(tcomp, 2), '# spikes': nspikes, 'Latency (ms)': lat * 1e3 if isinstance(lat, float) else 'N/A', 'Spike rate (sp/ms)': sr * 1e-3 if isinstance(sr, float) else 'N/A' } if xlslog(logpath, logentry) == 1: logger.debug('log exported to "%s"', logpath) else: logger.error('log export to "%s" aborted', self.logpath) return outpath def findRheobaseAmps(self, DCs, Vthr, curr='net'): ''' Find the rheobase amplitudes (i.e. threshold amplitudes of infinite duration that would result in excitation) of a specific neuron for various stimulation duty cycles. :param DCs: duty cycles vector (-) :param Vthr: threshold membrane potential above which the neuron necessarily fires (mV) :return: rheobase amplitudes vector (mA/m2) ''' # Compute the pulse average net (or leakage) current along the amplitude space if curr == 'net': iNet = self.iNet(Vthr, self.steadyStates(Vthr)) elif curr == 'leak': iNet = self.iLeak(Vthr) # Compute rheobase amplitudes return iNet / np.array(DCs) diff --git a/PySONIC/neurons/cortical.py b/PySONIC/neurons/cortical.py index 8af6936..ed76952 100644 --- a/PySONIC/neurons/cortical.py +++ b/PySONIC/neurons/cortical.py @@ -1,814 +1,810 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-07-31 15:19:51 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-03-13 15:07:59 +# @Last Modified time: 2019-03-13 15:40:25 import numpy as np from ..core import PointNeuron from ..utils import vtrap class Cortical(PointNeuron): ''' Class defining the generic membrane channel dynamics of a cortical neuron with 4 different current types: - Inward Sodium current - Outward, delayed-rectifier Potassium current - Outward, slow non.inactivating Potassium current - Non-specific leakage current This generic class cannot be used directly as it does not contain any specific parameters. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* ''' # Generic biophysical parameters of cortical cells Cm0 = 1e-2 # Cell membrane resting capacitance (F/m2) Vm0 = 0.0 # Dummy value for membrane potential (mV) VNa = 50.0 # Sodium Nernst potential (mV) VK = -90.0 # Potassium Nernst potential (mV) VCa = 120.0 # # Calcium Nernst potential (mV) def __init__(self): ''' Constructor of the class ''' # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 'p'] self.states0 = np.array([]) # Names of the different coefficients to be averaged in a lookup table. self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphap', 'betap'] - # Charge interval bounds for lookup creation - self.Qbounds = np.array([np.round(self.Vm0 - 25.0), 50.0]) * self.Cm0 * 1e-3 # C/m2 - - def alpham(self, Vm): ''' Compute the alpha rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = 0.32 * vtrap(13 - Vdiff, 4) # ms-1 return alpha * 1e3 # s-1 def betam(self, Vm): ''' Compute the beta rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = 0.28 * vtrap(Vdiff - 40, 5) # ms-1 return beta * 1e3 # s-1 def alphah(self, Vm): ''' Compute the alpha rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = (0.128 * np.exp(-(Vdiff - 17) / 18)) # ms-1 return alpha * 1e3 # s-1 def betah(self, Vm): ''' Compute the beta rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = (4 / (1 + np.exp(-(Vdiff - 40) / 5))) # ms-1 return beta * 1e3 # s-1 def alphan(self, Vm): ''' Compute the alpha rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = 0.032 * vtrap(15 - Vdiff, 5) # ms-1 return alpha * 1e3 # s-1 def betan(self, Vm): ''' Compute the beta rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = (0.5 * np.exp(-(Vdiff - 10) / 40)) # ms-1 return beta * 1e3 # s-1 def pinf(self, Vm): ''' Compute the asymptotic value of the open-probability of slow non-inactivating Potassium channels. :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1 + np.exp(-(Vm + 35) / 10)) # prob def taup(self, Vm): ''' Compute the decay time constant for adaptation of slow non-inactivating Potassium channels. :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' return self.TauMax / (3.3 * np.exp((Vm + 35) / 20) + np.exp(-(Vm + 35) / 20)) # s def derM(self, Vm, m): ''' Compute the evolution of the open-probability of Sodium channels. :param Vm: membrane potential (mV) :param m: open-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alpham(Vm) * (1 - m) - self.betam(Vm) * m def derH(self, Vm, h): ''' Compute the evolution of the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :param h: inactivation-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphah(Vm) * (1 - h) - self.betah(Vm) * h def derN(self, Vm, n): ''' Compute the evolution of the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :param n: open-probability of delayed-rectifier Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphan(Vm) * (1 - n) - self.betan(Vm) * n def derP(self, Vm, p): ''' Compute the evolution of the open-probability of slow non-inactivating Potassium channels. :param Vm: membrane potential (mV) :param p: open-probability of slow non-inactivating Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.pinf(Vm) - p) / self.taup(Vm) def iNa(self, m, h, Vm): ''' Sodium current :param m: open-probability of m-gate :param h: open-probability of h-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GNa = self.GNaMax * m**3 * h return GNa * (Vm - self.VNa) def iKd(self, n, Vm): ''' Delayed-rectifier Potassium current :param n: open-probability of n-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GKMax * n**4 return GK * (Vm - self.VK) def iM(self, p, Vm): ''' Slow non-inactivating Potassium current :param p: open-probability of p-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GM = self.GMMax * p return GM * (Vm - self.VK) def iLeak(self, Vm): ''' Non-specific leakage current :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GLeak * (Vm - self.VLeak) def currents(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, p = states return { 'iNa': self.iNa(m, h, Vm), 'iKd': self.iKd(n, Vm), 'iM': self.iM(p, Vm), 'iLeak': self.iLeak(Vm) } # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) peq = self.pinf(Vm) return np.array([meq, heq, neq, peq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, p = states dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dpdt = self.derP(Vm, p) return [dmdt, dhdt, dndt, dpdt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) Tp = self.taup(Vm) pinf = self.pinf(Vm) ap_avg = np.mean(pinf / Tp) bp_avg = np.mean(1 / Tp) - ap_avg # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, ap_avg, bp_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) m, h, n, p = states dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dpdt = rates[6] * (1 - p) - rates[7] * p return [dmdt, dhdt, dndt, dpdt] class CorticalRS(Cortical): ''' Specific membrane channel dynamics of a cortical regular spiking, excitatory pyramidal neuron. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* ''' # Name of channel mechanism name = 'RS' # Cell-specific biophysical parameters Vm0 = -71.9 # Cell membrane resting potential (mV) GNaMax = 560.0 # Max. conductance of Sodium current (S/m^2) GKMax = 60.0 # Max. conductance of delayed Potassium current (S/m^2) GMMax = 0.75 # Max. conductance of slow non-inactivating Potassium current (S/m^2) GLeak = 0.205 # Conductance of non-specific leakage current (S/m^2) VLeak = -70.3 # Non-specific leakage Nernst potential (mV) VT = -56.2 # Spike threshold adjustment parameter (mV) TauMax = 0.608 # Max. adaptation decay of slow non-inactivating Potassium current (s) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_{Kd}\ kin.': ['n'], 'i_M\ kin.': ['p'] } def __init__(self): ''' Constructor of the class. ''' # Instantiate parent class super().__init__() # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) class CorticalFS(Cortical): ''' Specific membrane channel dynamics of a cortical fast-spiking, inhibitory neuron. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* ''' # Name of channel mechanism name = 'FS' # Cell-specific biophysical parameters Vm0 = -71.4 # Cell membrane resting potential (mV) GNaMax = 580.0 # Max. conductance of Sodium current (S/m^2) GKMax = 39.0 # Max. conductance of delayed Potassium current (S/m^2) GMMax = 0.787 # Max. conductance of slow non-inactivating Potassium current (S/m^2) GLeak = 0.38 # Conductance of non-specific leakage current (S/m^2) VLeak = -70.4 # Non-specific leakage Nernst potential (mV) VT = -57.9 # Spike threshold adjustment parameter (mV) TauMax = 0.502 # Max. adaptation decay of slow non-inactivating Potassium current (s) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_{Kd}\ kin.': ['n'], 'i_M\ kin.': ['p'] } def __init__(self): ''' Constructor of the class. ''' # Instantiate parent class super().__init__() # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) class CorticalLTS(Cortical): ''' Specific membrane channel dynamics of a cortical low-threshold spiking, inhibitory neuron with an additional inward Calcium current due to the presence of a T-type channel. References: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* *Huguenard, J.R., and McCormick, D.A. (1992). Simulation of the currents involved in rhythmic oscillations in thalamic relay neurons. J. Neurophysiol. 68, 1373–1383.* ''' # Name of channel mechanism name = 'LTS' # Cell-specific biophysical parameters Vm0 = -54.0 # Cell membrane resting potential (mV) GNaMax = 500.0 # Max. conductance of Sodium current (S/m^2) GKMax = 40.0 # Max. conductance of delayed Potassium current (S/m^2) GMMax = 0.28 # Max. conductance of slow non-inactivating Potassium current (S/m^2) GTMax = 4.0 # Max. conductance of low-threshold Calcium current (S/m^2) GLeak = 0.19 # Conductance of non-specific leakage current (S/m^2) VLeak = -50.0 # Non-specific leakage Nernst potential (mV) VT = -50.0 # Spike threshold adjustment parameter (mV) TauMax = 4.0 # Max. adaptation decay of slow non-inactivating Potassium current (s) Vx = -7.0 # Voltage-dependence uniform shift factor at 36°C (mV) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_{Kd}\ kin.': ['n'], 'i_M\ kin.': ['p'], 'i_{CaT}\ kin.': ['s', 'u'] } def __init__(self): ''' Constructor of the class. ''' # Instantiate parent class super().__init__() # Add names of cell-specific Calcium channel probabilities self.states_names += ['s', 'u'] # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) # Define the names of the different coefficients to be averaged in a lookup table. self.coeff_names += ['alphas', 'betas', 'alphau', 'betau'] def sinf(self, Vm): ''' Compute the asymptotic value of the open-probability of the S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1.0 + np.exp(-(Vm + self.Vx + 57.0) / 6.2)) # prob def taus(self, Vm): ''' Compute the decay time constant for adaptation of S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' tmp = np.exp(-(Vm + self.Vx + 132.0) / 16.7) + np.exp((Vm + self.Vx + 16.8) / 18.2) return 1.0 / 3.7 * (0.612 + 1.0 / tmp) * 1e-3 # s def uinf(self, Vm): ''' Compute the asymptotic value of the open-probability of the U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1.0 + np.exp((Vm + self.Vx + 81.0) / 4.0)) # prob def tauu(self, Vm): ''' Compute the decay time constant for adaptation of U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' if Vm + self.Vx < -80.0: return 1.0 / 3.7 * np.exp((Vm + self.Vx + 467.0) / 66.6) * 1e-3 # s else: return 1.0 / 3.7 * (np.exp(-(Vm + self.Vx + 22) / 10.5) + 28.0) * 1e-3 # s def derS(self, Vm, s): ''' Compute the evolution of the open-probability of the S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :param s: open-probability of S-type Calcium activation gates (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.sinf(Vm) - s) / self.taus(Vm) def derU(self, Vm, u): ''' Compute the evolution of the open-probability of the U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :param u: open-probability of U-type Calcium inactivation gates (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.uinf(Vm) - u) / self.tauu(Vm) def iCaT(self, s, u, Vm): ''' Low-threshold (T-type) Calcium current :param s: open-probability of s-gate :param u: open-probability of u-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GT = self.GTMax * s**2 * u return GT * (Vm - self.VCa) def currents(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, p, s, u = states currents = super().currents(Vm, [m, h, n, p]) currents['iCaT'] = self.iCaT(s, u, Vm) # mA/m2 return currents def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Call parent method to compute Sodium and Potassium channels gates steady-states NaK_eqstates = super().steadyStates(Vm) # Compute Calcium channel gates steady-states Ca_eqstates = np.array([self.sinf(Vm), self.uinf(Vm)]) # Merge all steady-states and return return np.concatenate((NaK_eqstates, Ca_eqstates)) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack input states *NaK_states, s, u = states # Call parent method to compute Sodium and Potassium channels states derivatives NaK_derstates = super().derStates(Vm, NaK_states) # Compute Calcium channels states derivatives dsdt = self.derS(Vm, s) dudt = self.derU(Vm, u) # Merge all states derivatives and return return NaK_derstates + [dsdt, dudt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Call parent method to compute Sodium and Potassium effective rate constants NaK_rates = super().getEffRates(Vm) # Compute Calcium effective rate constants Ts = self.taus(Vm) as_avg = np.mean(self.sinf(Vm) / Ts) bs_avg = np.mean(1 / Ts) - as_avg Tu = np.array([self.tauu(v) for v in Vm]) au_avg = np.mean(self.uinf(Vm) / Tu) bu_avg = np.mean(1 / Tu) - au_avg Ca_rates = np.array([as_avg, bs_avg, au_avg, bu_avg]) # Merge all rates and return return np.concatenate((NaK_rates, Ca_rates)) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' # Unpack input states *NaK_states, s, u = states # Call parent method to compute Sodium and Potassium channels states derivatives NaK_dstates = super().derStatesEff(Qm, NaK_states, interp_data) # Compute Calcium channels states derivatives Ca_rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names[8:]]) dsdt = Ca_rates[0] * (1 - s) - Ca_rates[1] * s dudt = Ca_rates[2] * (1 - u) - Ca_rates[3] * u # Merge all states derivatives and return return NaK_dstates + [dsdt, dudt] class CorticalIB(Cortical): ''' Specific membrane channel dynamics of a cortical intrinsically bursting neuron with an additional inward Calcium current due to the presence of a L-type channel. References: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* *Reuveni, I., Friedman, A., Amitai, Y., and Gutnick, M.J. (1993). Stepwise repolarization from Ca2+ plateaus in neocortical pyramidal cells: evidence for nonhomogeneous distribution of HVA Ca2+ channels in dendrites. J. Neurosci. 13, 4609–4621.* ''' # Name of channel mechanism name = 'IB' # Cell-specific biophysical parameters Vm0 = -71.4 # Cell membrane resting potential (mV) GNaMax = 500 # Max. conductance of Sodium current (S/m^2) GKMax = 50 # Max. conductance of delayed Potassium current (S/m^2) GMMax = 0.3 # Max. conductance of slow non-inactivating Potassium current (S/m^2) GCaLMax = 1.0 # Max. conductance of L-type Calcium current (S/m^2) GLeak = 0.1 # Conductance of non-specific leakage current (S/m^2) VLeak = -70 # Non-specific leakage Nernst potential (mV) VT = -56.2 # Spike threshold adjustment parameter (mV) TauMax = 0.608 # Max. adaptation decay of slow non-inactivating Potassium current (s) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_{Kd}\ kin.': ['n'], 'i_M\ kin.': ['p'], 'i_{CaL}\ kin.': ['q', 'r'] } def __init__(self): ''' Constructor of the class. ''' # Instantiate parent class super().__init__() # Add names of cell-specific Calcium channel probabilities self.states_names += ['q', 'r'] # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) # Define the names of the different coefficients to be averaged in a lookup table. self.coeff_names += ['alphaq', 'betaq', 'alphar', 'betar'] def alphaq(self, Vm): ''' Compute the alpha rate for the open-probability of L-type Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = 0.055 * vtrap(-(Vm + 27), 3.8) # ms-1 return alpha * 1e3 # s-1 def betaq(self, Vm): ''' Compute the beta rate for the open-probability of L-type Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 0.94 * np.exp(-(Vm + 75) / 17) # ms-1 return beta * 1e3 # s-1 def alphar(self, Vm): ''' Compute the alpha rate for the inactivation-probability of L-type Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = 0.000457 * np.exp(-(Vm + 13) / 50) # ms-1 return alpha * 1e3 # s-1 def betar(self, Vm): ''' Compute the beta rate for the inactivation-probability of L-type Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 0.0065 / (np.exp(-(Vm + 15) / 28) + 1) # ms-1 return beta * 1e3 # s-1 def derQ(self, Vm, q): ''' Compute the evolution of the open-probability of the Q (activation) gate of L-type Calcium channels. :param Vm: membrane potential (mV) :param q: open-probability of Q gate (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphaq(Vm) * (1 - q) - self.betaq(Vm) * q def derR(self, Vm, r): ''' Compute the evolution of the open-probability of the R (inactivation) gate of L-type Calcium channels. :param Vm: membrane potential (mV) :param r: open-probability of R gate (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphar(Vm) * (1 - r) - self.betar(Vm) * r def iCaL(self, q, r, Vm): ''' High-threshold (L-type) Calcium current :param q: open-probability of q-gate :param r: open-probability of r-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GCaL = self.GCaLMax * q**2 * r return GCaL * (Vm - self.VCa) def currents(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, p, q, r = states return { 'iNa': self.iNa(m, h, Vm), 'iKd': self.iKd(n, Vm), 'iM': self.iM(p, Vm), 'iCaL': self.iCaL(q, r, Vm), 'iLeak': self.iLeak(Vm) } # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Call parent method to compute Sodium and Potassium channels gates steady-states NaK_eqstates = super().steadyStates(Vm) # Compute L-type Calcium channel gates steady-states qeq = self.alphaq(Vm) / (self.alphaq(Vm) + self.betaq(Vm)) req = self.alphar(Vm) / (self.alphar(Vm) + self.betar(Vm)) CaL_eqstates = np.array([qeq, req]) # Merge all steady-states and return return np.concatenate((NaK_eqstates, CaL_eqstates)) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack input states *NaK_states, q, r = states # Call parent method to compute Sodium and Potassium channels states derivatives NaK_derstates = super().derStates(Vm, NaK_states) # Compute L-type Calcium channels states derivatives dqdt = self.derQ(Vm, q) drdt = self.derR(Vm, r) # Merge all states derivatives and return return NaK_derstates + [dqdt, drdt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Call parent method to compute Sodium and Potassium effective rate constants NaK_rates = super().getEffRates(Vm) # Compute Calcium effective rate constants aq_avg = np.mean(self.alphaq(Vm)) bq_avg = np.mean(self.betaq(Vm)) ar_avg = np.mean(self.alphar(Vm)) br_avg = np.mean(self.betar(Vm)) CaL_rates = np.array([aq_avg, bq_avg, ar_avg, br_avg]) # Merge all rates and return return np.concatenate((NaK_rates, CaL_rates)) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' # Unpack input states *NaK_states, q, r = states # Call parent method to compute Sodium and Potassium channels states derivatives NaK_dstates = super().derStatesEff(Qm, NaK_states, interp_data) # Compute Calcium channels states derivatives CaL_rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names[8:]]) dqdt = CaL_rates[0] * (1 - q) - CaL_rates[1] * q drdt = CaL_rates[2] * (1 - r) - CaL_rates[3] * r # Merge all states derivatives and return return NaK_dstates + [dqdt, drdt] diff --git a/PySONIC/neurons/fh.py b/PySONIC/neurons/fh.py index 32a0abe..75cc252 100644 --- a/PySONIC/neurons/fh.py +++ b/PySONIC/neurons/fh.py @@ -1,298 +1,295 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2019-01-07 18:41:06 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-03-13 14:36:34 +# @Last Modified time: 2019-03-13 15:40:33 import numpy as np from ..core import PointNeuron from ..utils import vtrap, ghkDrive from ..constants import Celsius2Kelvin, Z_Na, Z_K class FrankenhaeuserHuxley(PointNeuron): ''' Class defining the membrane channel dynamics of a Xenopus myelinated neuron with 4 different current types: - Inward Sodium current - Outward, delayed-rectifier Potassium current - Non-specific delayed current - Non-specific leakage current Reference: *Frankenhaeuser, B., and Huxley, A.F. (1964). The action potential in the myelinated nerve fibre of Xenopus laevis as computed on the basis of voltage clamp data. J Physiol 171, 302–315.* ''' # Name of channel mechanism name = 'FH' # Cell biophysical parameters Cm0 = 2e-2 # Cell membrane resting capacitance (F/m2) Vm0 = -70. # Membrane resting potential (mV) celsius = 20.0 # Temperature (Celsius) GLeak = 300.3 # Leakage conductance (S/m2) VLeak = -69.974 # Leakage resting potential (mV) pNabar = 8e-5 # Sodium permeability constant (m/s) pPbar = .54e-5 # Non-specific permeability constant (m/s) pKbar = 1.2e-5 # Potassium permeability constant (m/s) nai = 13.74e-3 # Sodium intracellular concentration (M) nao = 114.5e-3 # Sodium extracellular concentration (M) ki = 120e-3 # Potassium intracellular concentration (M) ko = 2.5e-3 # Potassium extracellular concentration (M) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_{Kd}\ kin.': ['n'], 'i_P\ kin.': ['p'] } def __init__(self): ''' Constructor of the class ''' self.q10 = 3**((self.celsius - 20) / 10) self.T = self.celsius + Celsius2Kelvin # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 'p'] self.states0 = self.steadyStates(self.Vm0) # Names of the different coefficients to be averaged in a lookup table. self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphap', 'betap'] - # Charge interval bounds for lookup creation - self.Qbounds = np.array([np.round(self.Vm0 - 25.0), 50.0]) * self.Cm0 * 1e-3 # C/m2 - def alpham(self, Vm): ''' Compute the alpha rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.Vm0 alpha = 0.36 * vtrap(22. - Vdiff, 3.) # ms-1 return self.q10 * alpha * 1e3 # s-1 def betam(self, Vm): ''' Compute the beta rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.Vm0 beta = 0.4 * vtrap(Vdiff - 13., 20.) # ms-1 return self.q10 * beta * 1e3 # s-1 def alphah(self, Vm): ''' Compute the alpha rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.Vm0 alpha = 0.1 * vtrap(Vdiff + 10.0, 6.) # ms-1 return self.q10 * alpha * 1e3 # s-1 def betah(self, Vm): ''' Compute the beta rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.Vm0 beta = 4.5 / (np.exp((45. - Vdiff) / 10.) + 1) # ms-1 return self.q10 * beta * 1e3 # s-1 def alphan(self, Vm): ''' Compute the alpha rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.Vm0 alpha = 0.02 * vtrap(35. - Vdiff, 10.0) # ms-1 return self.q10 * alpha * 1e3 # s-1 def betan(self, Vm): ''' Compute the beta rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.Vm0 beta = 0.05 * vtrap(Vdiff - 10., 10.) # ms-1 return self.q10 * beta * 1e3 # s-1 def alphap(self, Vm): ''' Compute the alpha rate for the open-probability of non-specific delayed current. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.Vm0 alpha = 0.006 * vtrap(40. - Vdiff, 10.0) # ms-1 return self.q10 * alpha * 1e3 # s-1 def betap(self, Vm): ''' Compute the beta rate for the open-probability of non-specific delayed current. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.Vm0 beta = 0.09 * vtrap(Vdiff + 25., 20.) # ms-1 return self.q10 * beta * 1e3 # s-1 def derM(self, Vm, m): ''' Compute the evolution of the open-probability of Sodium channels. :param Vm: membrane potential (mV) :param m: open-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alpham(Vm) * (1 - m) - self.betam(Vm) * m def derH(self, Vm, h): ''' Compute the evolution of the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :param h: inactivation-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphah(Vm) * (1 - h) - self.betah(Vm) * h def derN(self, Vm, n): ''' Compute the evolution of the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :param n: open-probability of delayed-rectifier Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphan(Vm) * (1 - n) - self.betan(Vm) * n def derP(self, Vm, p): ''' Compute the evolution of the open-probability of non-specific delayed current. :param Vm: membrane potential (mV) :param p: open-probability (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphap(Vm) * (1 - p) - self.betap(Vm) * p def iNa(self, m, h, Vm): ''' Sodium current :param m: open-probability of m-gate :param h: open-probability of h-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' iNa_drive = ghkDrive(Vm, Z_Na, self.nai, self.nao, self.T) # mC/m3 return self.pNabar * m**2 * h * iNa_drive # mA/m2 def iKd(self, n, Vm): ''' Delayed-rectifier Potassium current :param n: open-probability of n-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' iKd_drive = ghkDrive(Vm, Z_K, self.ki, self.ko, self.T) # mC/m3 return self.pKbar * n**2 * iKd_drive # mA/m2 def iP(self, p, Vm): ''' Non-specific delayed current :param p: open-probability of p-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' iP_drive = ghkDrive(Vm, Z_Na, self.nai, self.nao, self.T) # mC/m3 return self.pPbar * p**2 * iP_drive # mA/m2 def iLeak(self, Vm): ''' Non-specific leakage current :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GLeak * (Vm - self.VLeak) def currents(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, p = states return { 'iNa': self.iNa(m, h, Vm), 'iKd': self.iKd(n, Vm), 'iP': self.iP(p, Vm), 'iLeak': self.iLeak(Vm) } # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) peq = self.alphap(Vm) / (self.alphap(Vm) + self.betap(Vm)) return np.array([meq, heq, neq, peq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, p = states dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dpdt = self.derP(Vm, p) return [dmdt, dhdt, dndt, dpdt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) ap_avg = np.mean(self.alphap(Vm)) bp_avg = np.mean(self.betap(Vm)) # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, ap_avg, bp_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) m, h, n, p = states dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dpdt = rates[6] * (1 - p) - rates[7] * p return [dmdt, dhdt, dndt, dpdt] diff --git a/PySONIC/neurons/leech.py b/PySONIC/neurons/leech.py index e6e30f2..9fcc91c 100644 --- a/PySONIC/neurons/leech.py +++ b/PySONIC/neurons/leech.py @@ -1,1110 +1,1102 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-07-31 15:20:54 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-03-13 15:10:25 +# @Last Modified time: 2019-03-13 15:41:16 from functools import partialmethod import numpy as np from ..core import PointNeuron from ..constants import FARADAY, Rg, Z_Na, Z_Ca from ..utils import nernst class LeechTouch(PointNeuron): ''' Class defining the membrane channel dynamics of a leech touch sensory neuron. with 4 different current types: - Inward Sodium current - Outward Potassium current - Inward Calcium current - Non-specific leakage current - Calcium-dependent, outward Potassium current - Outward, Sodium pumping current Reference: *Cataldo, E., Brunelli, M., Byrne, J.H., Av-Ron, E., Cai, Y., and Baxter, D.A. (2005). Computational model of touch sensory cells (T Cells) of the leech: role of the afterhyperpolarization (AHP) in activity-dependent conduction failure. J Comput Neurosci 18, 5–24.* ''' # Name of channel mechanism name = 'LeechT' # Cell-specific biophysical parameters Cm0 = 1e-2 # Cell membrane resting capacitance (F/m2) Vm0 = -53.58 # Cell membrane resting potential (mV) VNa = 45.0 # Sodium Nernst potential (mV) VK = -62.0 # Potassium Nernst potential (mV) VCa = 60.0 # Calcium Nernst potential (mV) VLeak = -48.0 # Non-specific leakage Nernst potential (mV) VPumpNa = -300.0 # Sodium pump current reversal potential (mV) GNaMax = 3500.0 # Max. conductance of Sodium current (S/m^2) GKMax = 900.0 # Max. conductance of Potassium current (S/m^2) GCaMax = 20.0 # Max. conductance of Calcium current (S/m^2) GKCaMax = 236.0 # Max. conductance of Calcium-dependent Potassium current (S/m^2) GLeak = 1.0 # Conductance of non-specific leakage current (S/m^2) GPumpNa = 20.0 # Max. conductance of Sodium pump current (S/m^2) taum = 0.1e-3 # Sodium activation time constant (s) taus = 0.6e-3 # Calcium activation time constant (s) # Original conversion constants from inward ion current (nA) to build-up of # intracellular ion concentration (arb.) K_Na_original = 0.016 # iNa to intracellular [Na+] K_Ca_original = 0.1 # iCa to intracellular [Ca2+] # Constants needed to convert K from original model (soma compartment) # to current model (point-neuron) surface = 6434.0e-12 # surface of cell assumed as a single soma (m2) curr_factor = 1e6 # mA to nA # Time constants for the removal of ions from intracellular pools (s) tau_Na_removal = 16.0 # Na+ removal tau_Ca_removal = 1.25 # Ca2+ removal # Time constants for the iPumpNa and iKCa currents activation # from specific intracellular ions (s) tau_PumpNa_act = 0.1 # iPumpNa activation from intracellular Na+ tau_KCa_act = 0.01 # iKCa activation from intracellular Ca2+ # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_K\ kin.': ['n'], 'i_{Ca}\ kin.': ['s'], 'pools': ['C_Na_arb', 'C_Na_arb_activation', 'C_Ca_arb', 'C_Ca_arb_activation'] } def __init__(self): ''' Constructor of the class. ''' # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'C_Na', 'A_Na', 'C_Ca', 'A_Ca'] self.states0 = np.array([]) # Names of the channels effective coefficients self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas'] self.K_Na = self.K_Na_original * self.surface * self.curr_factor self.K_Ca = self.K_Ca_original * self.surface * self.curr_factor # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) - # Charge interval bounds for lookup creation - self.Qbounds = np.array([np.round(self.Vm0 - 10.0), 50.0]) * self.Cm0 * 1e-3 # C/m2 - # ----------------- Generic ----------------- def _xinf(self, Vm, halfmax, slope, power): ''' Generic function computing the steady-state activation/inactivation of a particular ion channel at a given voltage. :param Vm: membrane potential (mV) :param halfmax: half-(in)activation voltage (mV) :param slope: slope parameter of (in)activation function (mV) :param power: power exponent multiplying the exponential expression (integer) :return: steady-state (in)activation (-) ''' return 1 / (1 + np.exp((Vm - halfmax) / slope))**power def _taux(self, Vm, halfmax, slope, tauMax, tauMin): ''' Generic function computing the voltage-dependent, activation/inactivation time constant of a particular ion channel at a given voltage. :param Vm: membrane potential (mV) :param halfmax: voltage at which (in)activation time constant is half-maximal (mV) :param slope: slope parameter of (in)activation time constant function (mV) :return: steady-state (in)activation (-) ''' return (tauMax - tauMin) / (1 + np.exp((Vm - halfmax) / slope)) + tauMin def _derC_ion(self, Cion, Iion, Kion, tau): ''' Generic function computing the time derivative of the concentration of a specific ion in its intracellular pool. :param Cion: ion concentration in the pool (arbitrary unit) :param Iion: ionic current (mA/m2) :param Kion: scaling factor for current contribution to pool (arb. unit / nA???) :param tau: time constant for removal of ions from the pool (s) :return: variation of ionic concentration in the pool (arbitrary unit /s) ''' return (Kion * (-Iion) - Cion) / tau def _derA_ion(self, Aion, Cion, tau): ''' Generic function computing the time derivative of the concentration and time dependent activation function, for a specific pool-dependent ionic current. :param Aion: concentration and time dependent activation function (arbitrary unit) :param Cion: ion concentration in the pool (arbitrary unit) :param tau: time constant for activation function variation (s) :return: variation of activation function (arbitrary unit / s) ''' return (Cion - Aion) / tau # ------------------ Na ------------------- minf = partialmethod(_xinf, halfmax=-35.0, slope=-5.0, power=1) hinf = partialmethod(_xinf, halfmax=-50.0, slope=9.0, power=2) tauh = partialmethod(_taux, halfmax=-36.0, slope=3.5, tauMax=14.0e-3, tauMin=0.2e-3) def derM(self, Vm, m): ''' Instantaneous derivative of Sodium activation. ''' return (self.minf(Vm) - m) / self.taum # s-1 def derH(self, Vm, h): ''' Instantaneous derivative of Sodium inactivation. ''' return (self.hinf(Vm) - h) / self.tauh(Vm) # s-1 # ------------------ K ------------------- ninf = partialmethod(_xinf, halfmax=-22.0, slope=-9.0, power=1) taun = partialmethod(_taux, halfmax=-10.0, slope=10.0, tauMax=6.0e-3, tauMin=1.0e-3) def derN(self, Vm, n): ''' Instantaneous derivative of Potassium activation. ''' return (self.ninf(Vm) - n) / self.taun(Vm) # s-1 # ------------------ Ca ------------------- sinf = partialmethod(_xinf, halfmax=-10.0, slope=-2.8, power=1) def derS(self, Vm, s): ''' Instantaneous derivative of Calcium activation. ''' return (self.sinf(Vm) - s) / self.taus # s-1 # ------------------ Pools ------------------- def derC_Na(self, C_Na, I_Na): ''' Derivative of Sodium concentration in intracellular pool. ''' return self._derC_ion(C_Na, I_Na, self.K_Na, self.tau_Na_removal) def derA_Na(self, A_Na, C_Na): ''' Derivative of Sodium pool-dependent activation function for iPumpNa. ''' return self._derA_ion(A_Na, C_Na, self.tau_PumpNa_act) def derC_Ca(self, C_Ca, I_Ca): ''' Derivative of Calcium concentration in intracellular pool. ''' return self._derC_ion(C_Ca, I_Ca, self.K_Ca, self.tau_Ca_removal) def derA_Ca(self, A_Ca, C_Ca): ''' Derivative of Calcium pool-dependent activation function for iKCa. ''' return self._derA_ion(A_Ca, C_Ca, self.tau_KCa_act) # ------------------ Currents ------------------- def iNa(self, m, h, Vm): ''' Sodium current :param m: open-probability of m-gate :param h: open-probability of h-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GNaMax * m**3 * h * (Vm - self.VNa) def iK(self, n, Vm): ''' Delayed-rectifier Potassium current :param n: open-probability of n-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GKMax * n**2 * (Vm - self.VK) def iCa(self, s, Vm): ''' Calcium current :param s: open-probability of s-gate :param u: open-probability of u-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' ''' Calcium inward current. ''' return self.GCaMax * s * (Vm - self.VCa) def iKCa(self, A_Ca, Vm): ''' Calcium-activated Potassium current :param A_Ca: Calcium pool-dependent gate opening :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GKCaMax * A_Ca * (Vm - self.VK) def iPumpNa(self, A_Na, Vm): ''' NaK-ATPase pump current :param A_Na: Sodium pool-dependent gate opening. :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GPumpNa * A_Na * (Vm - self.VPumpNa) def iLeak(self, Vm): ''' Non-specific leakage current :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GLeak * (Vm - self.VLeak) def currents(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, _, A_Na, _, A_Ca = states return { 'iNa': self.iNa(m, h, Vm), 'iK': self.iK(n, Vm), 'iCa': self.iCa(s, Vm), 'iLeak': self.iLeak(Vm), 'iPumpNa': self.iPumpNa(A_Na, Vm), 'iKCa': self.iKCa(A_Ca, Vm) } # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Standard gating dynamics: Solve the equation dx/dt = 0 at Vm for each x-state meq = self.minf(Vm) heq = self.hinf(Vm) neq = self.ninf(Vm) seq = self.sinf(Vm) # PumpNa pool concentration and activation steady-state INa_eq = self.iNa(meq, heq, Vm) CNa_eq = self.K_Na * (-INa_eq) ANa_eq = CNa_eq # KCa current pool concentration and activation steady-state ICa_eq = self.iCa(seq, Vm) CCa_eq = self.K_Ca * (-ICa_eq) ACa_eq = CCa_eq return np.array([meq, heq, neq, seq, CNa_eq, ANa_eq, CCa_eq, ACa_eq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack states m, h, n, s, C_Na, A_Na, C_Ca, A_Ca = states # Standard gating states derivatives dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) # PumpNa current pool concentration and activation state I_Na = self.iNa(m, h, Vm) dCNa_dt = self.derC_Na(C_Na, I_Na) dANa_dt = self.derA_Na(A_Na, C_Na) # KCa current pool concentration and activation state I_Ca = self.iCa(s, Vm) dCCa_dt = self.derC_Ca(C_Ca, I_Ca) dACa_dt = self.derA_Ca(A_Ca, C_Ca) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dCNa_dt, dANa_dt, dCCa_dt, dACa_dt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants Tm = self.taum minf = self.minf(Vm) am_avg = np.mean(minf / Tm) bm_avg = np.mean(1 / Tm) - am_avg Th = self.tauh(Vm) hinf = self.hinf(Vm) ah_avg = np.mean(hinf / Th) bh_avg = np.mean(1 / Th) - ah_avg Tn = self.taun(Vm) ninf = self.ninf(Vm) an_avg = np.mean(ninf / Tn) bn_avg = np.mean(1 / Tn) - an_avg Ts = self.taus sinf = self.sinf(Vm) as_avg = np.mean(sinf / Ts) bs_avg = np.mean(1 / Ts) - as_avg # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) Vmeff = np.interp(Qm, interp_data['Q'], interp_data['V']) # Unpack states m, h, n, s, C_Na, A_Na, C_Ca, A_Ca = states # Standard gating states derivatives dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s # PumpNa current pool concentration and activation state I_Na = self.iNa(m, h, Vmeff) dCNa_dt = self.derC_Na(C_Na, I_Na) dANa_dt = self.derA_Na(A_Na, C_Na) # KCa current pool concentration and activation state I_Ca_eff = self.iCa(s, Vmeff) dCCa_dt = self.derC_Ca(C_Ca, I_Ca_eff) dACa_dt = self.derA_Ca(A_Ca, C_Ca) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dCNa_dt, dANa_dt, dCCa_dt, dACa_dt] class LeechMech(PointNeuron): ''' Class defining the basic dynamics of Sodium, Potassium and Calcium channels for several neurons of the leech. Reference: *Baccus, S.A. (1998). Synaptic facilitation by reflected action potentials: enhancement of transmission when nerve impulses reverse direction at axon branch points. Proc. Natl. Acad. Sci. U.S.A. 95, 8345–8350.* ''' alphaC_sf = 1e-5 # Calcium activation rate constant scaling factor (M) betaC = 0.1e3 # beta rate for the open-probability of Ca2+-dependent Potassium channels (s-1) T = 293.15 # Room temperature (K) def alpham(self, Vm): ''' Compute the alpha rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = -0.03 * (Vm + 28) / (np.exp(- (Vm + 28) / 15) - 1) # ms-1 return alpha * 1e3 # s-1 def betam(self, Vm): ''' Compute the beta rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 2.7 * np.exp(-(Vm + 53) / 18) # ms-1 return beta * 1e3 # s-1 def alphah(self, Vm): ''' Compute the alpha rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = 0.045 * np.exp(-(Vm + 58) / 18) # ms-1 return alpha * 1e3 # s-1 def betah(self, Vm): ''' Compute the beta rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) .. warning:: the original paper contains an error (multiplication) in the expression of this rate constant, corrected in the mod file on ModelDB (division). ''' beta = 0.72 / (np.exp(-(Vm + 23) / 14) + 1) # ms-1 return beta * 1e3 # s-1 def alphan(self, Vm): ''' Compute the alpha rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = -0.024 * (Vm - 17) / (np.exp(-(Vm - 17) / 8) - 1) # ms-1 return alpha * 1e3 # s-1 def betan(self, Vm): ''' Compute the beta rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 0.2 * np.exp(-(Vm + 48) / 35) # ms-1 return beta * 1e3 # s-1 def alphas(self, Vm): ''' Compute the alpha rate for the open-probability of Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = -1.5 * (Vm - 20) / (np.exp(-(Vm - 20) / 5) - 1) # ms-1 return alpha * 1e3 # s-1 def betas(self, Vm): ''' Compute the beta rate for the open-probability of Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 1.5 * np.exp(-(Vm + 25) / 10) # ms-1 return beta * 1e3 # s-1 def alphaC(self, C_Ca_in): ''' Compute the alpha rate for the open-probability of Calcium-dependent Potassium channels. :param C_Ca_in: intracellular Calcium concentration (M) :return: rate constant (s-1) ''' alpha = 0.1 * C_Ca_in / self.alphaC_sf # ms-1 return alpha * 1e3 # s-1 def derM(self, Vm, m): ''' Compute the evolution of the open-probability of Sodium channels. :param Vm: membrane potential (mV) :param m: open-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alpham(Vm) * (1 - m) - self.betam(Vm) * m def derH(self, Vm, h): ''' Compute the evolution of the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :param h: inactivation-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphah(Vm) * (1 - h) - self.betah(Vm) * h def derN(self, Vm, n): ''' Compute the evolution of the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :param n: open-probability of delayed-rectifier Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphan(Vm) * (1 - n) - self.betan(Vm) * n def derS(self, Vm, s): ''' Compute the evolution of the open-probability of Calcium channels. :param Vm: membrane potential (mV) :param s: open-probability of Calcium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphas(Vm) * (1 - s) - self.betas(Vm) * s def derC(self, c, C_Ca_in): ''' Compute the evolution of the open-probability of Calcium-dependent Potassium channels. :param c: open-probability of Calcium-dependent Potassium channels (prob) :param C_Ca_in: intracellular Calcium concentration (M) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphaC(C_Ca_in) * (1 - c) - self.betaC * c def iNa(self, m, h, Vm, C_Na_in): ''' Sodium current :param m: open-probability of m-gate :param h: open-probability of Sodium channels :param Vm: membrane potential (mV) :param C_Na_in: intracellular Sodium concentration (M) :return: current per unit area (mA/m2) ''' GNa = self.GNaMax * m**4 * h VNa = nernst(Z_Na, C_Na_in, self.C_Na_out, self.T) # mV return GNa * (Vm - VNa) def iK(self, n, Vm): ''' Delayed-rectifier Potassium current :param n: open-probability of n-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GKMax * n**2 return GK * (Vm - self.VK) def iCa(self, s, Vm, C_Ca_in): ''' Calcium current :param s: open-probability of s-gate :param Vm: membrane potential (mV) :param C_Ca_in: intracellular Calcium concentration (M) :return: current per unit area (mA/m2) ''' GCa = self.GCaMax * s VCa = nernst(Z_Ca, C_Ca_in, self.C_Ca_out, self.T) # mV return GCa * (Vm - VCa) def iKCa(self, c, Vm): ''' Calcium-activated Potassium current :param c: open-probability of c-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GKCa = self.GKCaMax * c return GKCa * (Vm - self.VK) def iLeak(self, Vm): ''' Non-specific leakage current :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GLeak * (Vm - self.VLeak) class LeechPressure(LeechMech): ''' Class defining the membrane channel dynamics of a leech pressure sensory neuron. with 7 different current types: - Inward Sodium current - Outward Potassium current - Inward high-voltage-activated Calcium current - Non-specific leakage current - Calcium-dependent, outward Potassium current - Sodium pump current - Calcium pump current Reference: *Baccus, S.A. (1998). Synaptic facilitation by reflected action potentials: enhancement of transmission when nerve impulses reverse direction at axon branch points. Proc. Natl. Acad. Sci. U.S.A. 95, 8345–8350.* ''' # Name of channel mechanism name = 'LeechP' # Cell-specific biophysical parameters Cm0 = 1e-2 # Cell membrane resting capacitance (F/m2) Vm0 = -48.865 # Cell membrane resting potential (mV) C_Na_out = 0.11 # Sodium extracellular concentration (M) C_Ca_out = 1.8e-3 # Calcium extracellular concentration (M) C_Na_in0 = 0.01 # Initial Sodium intracellular concentration (M) C_Ca_in0 = 1e-7 # Initial Calcium intracellular concentration (M) # VNa = 60 # Sodium Nernst potential, from MOD file on ModelDB (mV) # VCa = 125 # Calcium Nernst potential, from MOD file on ModelDB (mV) VK = -68.0 # Potassium Nernst potential (mV) VLeak = -49.0 # Non-specific leakage Nernst potential (mV) INaPmax = 70.0 # Maximum pump rate of the NaK-ATPase (mA/m2) khalf_Na = 0.012 # Sodium concentration at which NaK-ATPase is at half its maximum rate (M) ksteep_Na = 1e-3 # Sensitivity of NaK-ATPase to varying Sodium concentrations (M) iCaS = 0.1 # Calcium pump current parameter (mA/m2) GNaMax = 3500.0 # Max. conductance of Sodium current (S/m^2) GKMax = 60.0 # Max. conductance of Potassium current (S/m^2) GCaMax = 0.02 # Max. conductance of Calcium current (S/m^2) GKCaMax = 8.0 # Max. conductance of Calcium-dependent Potassium current (S/m^2) GLeak = 5.0 # Conductance of non-specific leakage current (S/m^2) diam = 50e-6 # Cell soma diameter (m) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_K\ kin.': ['n'], 'i_{Ca}\ kin.': ['s'], 'i_{KCa}\ kin.': ['c'], 'pools': ['C_Na', 'C_Ca'] } def __init__(self): ''' Constructor of the class. ''' SV_ratio = 6 / self.diam # surface to volume ratio of the (spherical) cell soma # Conversion constant from membrane ionic currents into # change rate of intracellular ionic concentrations self.K_Na = SV_ratio / (Z_Na * FARADAY) * 1e-6 # Sodium (M/s) self.K_Ca = SV_ratio / (Z_Ca * FARADAY) * 1e-6 # Calcium (M/s) # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'c', 'C_Na', 'C_Ca'] self.states0 = np.array([]) # Names of the channels effective coefficients self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas'] # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) - # Charge interval bounds for lookup creation - self.Qbounds = np.array([np.round(self.Vm0 - 10.0), 60.0]) * self.Cm0 * 1e-3 # C/m2 - def iPumpNa(self, C_Na_in): ''' NaK-ATPase pump current :param C_Na_in: intracellular Sodium concentration (M) :return: current per unit area (mA/m2) ''' return self.INaPmax / (1 + np.exp((self.khalf_Na - C_Na_in) / self.ksteep_Na)) def iPumpCa(self, C_Ca_in): ''' Calcium pump current :param C_Ca_in: intracellular Calcium concentration (M) :return: current per unit area (mA/m2) ''' return self.iCaS * (C_Ca_in - self.C_Ca_in0) / 1.5 def currents(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, c, C_Na_in, C_Ca_in = states return { 'iNa': self.iNa(m, h, Vm, C_Na_in), 'iK': self.iK(n, Vm), 'iCa': self.iCa(s, Vm, C_Ca_in), 'iKCa': self.iKCa(c, Vm), 'iLeak': self.iLeak(Vm), 'iPumpNa': self.iPumpNa(C_Na_in) / 3., 'iPumpCa': self.iPumpCa(C_Ca_in) } # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Intracellular concentrations C_Na_eq = self.C_Na_in0 C_Ca_eq = self.C_Ca_in0 # Standard gating dynamics: Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) seq = self.alphas(Vm) / (self.alphas(Vm) + self.betas(Vm)) ceq = self.alphaC(C_Ca_eq) / (self.alphaC(C_Ca_eq) + self.betaC) return np.array([meq, heq, neq, seq, ceq, C_Na_eq, C_Ca_eq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack states m, h, n, s, c, C_Na_in, C_Ca_in = states # Standard gating states derivatives dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) dcdt = self.derC(c, C_Ca_in) # Intracellular concentrations dCNa_dt = - (self.iNa(m, h, Vm, C_Na_in) + self.iPumpNa(C_Na_in)) * self.K_Na # M/s dCCa_dt = -(self.iCa(s, Vm, C_Ca_in) + self.iPumpCa(C_Ca_in)) * self.K_Ca # M/s # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dCNa_dt, dCCa_dt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) as_avg = np.mean(self.alphas(Vm)) bs_avg = np.mean(self.betas(Vm)) # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) Vmeff = np.interp(Qm, interp_data['Q'], interp_data['V']) # Unpack states m, h, n, s, c, C_Na_in, C_Ca_in = states # Standard gating states derivatives dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s # KCa current gating state derivative dcdt = self.derC(c, C_Ca_in) # Intracellular concentrations dCNa_dt = - (self.iNa(m, h, Vmeff, C_Na_in) + self.iPumpNa(C_Na_in)) * self.K_Na # M/s dCCa_dt = -(self.iCa(s, Vmeff, C_Ca_in) + self.iPumpCa(C_Ca_in)) * self.K_Ca # M/s # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dCNa_dt, dCCa_dt] class LeechRetzius(LeechMech): ''' Class defining the membrane channel dynamics of a leech Retzius neuron. with 5 different current types: - Inward Sodium current - Outward Potassium current - Inward high-voltage-activated Calcium current - Non-specific leakage current - Calcium-dependent, outward Potassium current References: *Vazquez, Y., Mendez, B., Trueta, C., and De-Miguel, F.F. (2009). Summation of excitatory postsynaptic potentials in electrically-coupled neurones. Neuroscience 163, 202–212.* *ModelDB link: https://senselab.med.yale.edu/modeldb/ShowModel.cshtml?model=120910* ''' # Name of channel mechanism # name = 'LeechR' # Cell-specific biophysical parameters Cm0 = 5e-2 # Cell membrane resting capacitance (F/m2) Vm0 = -44.45 # Cell membrane resting potential (mV) VNa = 50.0 # Sodium Nernst potential, from retztemp.ses file on ModelDB (mV) VCa = 125.0 # Calcium Nernst potential, from cachdend.mod file on ModelDB (mV) VK = -79.0 # Potassium Nernst potential, from retztemp.ses file on ModelDB (mV) VLeak = -30.0 # Non-specific leakage Nernst potential, from leakdend.mod file on ModelDB (mV) GNaMax = 1250.0 # Max. conductance of Sodium current (S/m^2) GKMax = 10.0 # Max. conductance of Potassium current (S/m^2) GAMax = 100.0 # Max. conductance of transient Potassium current (S/m^2) GCaMax = 4.0 # Max. conductance of Calcium current (S/m^2) GKCaMax = 130.0 # Max. conductance of Calcium-dependent Potassium current (S/m^2) GLeak = 1.25 # Conductance of non-specific leakage current (S/m^2) Vhalf = -73.1 # mV C_Ca_in = 5e-8 # Calcium intracellular concentration, from retztemp.ses file (M) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_K\ kin.': ['n'], 'i_A\ kin.': ['a', 'b'], 'i_{Ca}\ kin.': ['s'], 'i_{KCa}\ kin.': ['c'] } def __init__(self): ''' Constructor of the class. ''' # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'c', 'a', 'b'] self.states0 = np.array([]) # Names of the channels effective coefficients self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas', 'alphac', 'betac', 'alphaa', 'betaa' 'alphab', 'betab'] # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) - self.Qbounds = np.array([np.round(self.Vm0 - 10.0), 50.0]) * self.Cm0 * 1e-3 # C/m2 - def ainf(self, Vm): ''' Steady-state activation probability of transient Potassium channels. Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' Vth = -55.0 # mV return 0 if Vm <= Vth else min(1, 2 * (Vm - Vth)**3 / ((11 - Vth)**3 + (Vm - Vth)**3)) def taua(self, Vm): ''' Activation time constant of transient Potassium channels. (assuming T = 20°C). Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' x = -1.5 * (Vm - self.Vhalf) * 1e-3 * 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): ''' Steady-state inactivation probability of transient Potassium channels. Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' return 1. / (1 + np.exp((self.Vhalf - Vm) / -6.3)) def taub(self, Vm): ''' Inactivation time constant of transient Potassium channels. (assuming T = 20°C). Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' x = 2 * (Vm - self.Vhalf) * 1e-3 * FARADAY / (Rg * self.T) alpha = np.exp(x) beta = np.exp(0.65 * x) return max(7.5, beta / (0.02 * (1 + alpha))) * 1e-3 # s def derA(self, Vm, a): ''' Compute the evolution of the activation-probability of transient Potassium channels. :param Vm: membrane potential (mV) :param a: activation-probability of transient Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.ainf(Vm) - a) / self.taua(Vm) def derB(self, Vm, b): ''' Compute the evolution of the inactivation-probability of transient Potassium channels. :param Vm: membrane potential (mV) :param b: inactivation-probability of transient Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.binf(Vm) - b) / self.taub(Vm) def iA(self, a, b, Vm): ''' Transient Potassium current :param a: open-probability of a-gate :param b: open-probability of b-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GAMax * a * b return GK * (Vm - self.VK) def currents(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, c, a, b = states return { 'iNa': self.iNa(m, h, Vm), 'iK': self.iK(n, Vm), 'iCa': self.iCa(s, Vm), 'iLeak': self.iLeak(Vm), 'iKCa': self.iKCa(c, Vm), 'iA': self.iA(a, b, Vm) } # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Standard gating dynamics: Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) seq = self.alphas(Vm) / (self.alphas(Vm) + self.betas(Vm)) ceq = self.alphaC(self.C_Ca_in) / (self.alphaC(self.C_Ca_in) + self.betaC) aeq = self.ainf(Vm) beq = self.binf(Vm) return np.array([meq, heq, neq, seq, ceq, aeq, beq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack states m, h, n, s, c, a, b = states # Standard gating states derivatives dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) dcdt = self.derC(c, self.C_Ca_in) dadt = self.derA(Vm, a) dbdt = self.derB(Vm, b) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dadt, dbdt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) as_avg = np.mean(self.alphas(Vm)) bs_avg = np.mean(self.betas(Vm)) Ta = self.taua(Vm) ainf = self.ainf(Vm) aa_avg = np.mean(ainf / Ta) ba_avg = np.mean(1 / Ta) - aa_avg Tb = self.taub(Vm) binf = self.binf(Vm) ab_avg = np.mean(binf / Tb) bb_avg = np.mean(1 / Tb) - ab_avg # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg, aa_avg, ba_avg, ab_avg, bb_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) # Unpack states m, h, n, s, c, a, b = states # Standard gating states derivatives dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s dadt = rates[8] * (1 - a) - rates[9] * a dbdt = rates[10] * (1 - b) - rates[11] * b # KCa current gating state derivative dcdt = self.derC(c, self.C_Ca_in) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dadt, dbdt] diff --git a/PySONIC/neurons/stn.py b/PySONIC/neurons/stn.py index 68dd26e..9310b56 100644 --- a/PySONIC/neurons/stn.py +++ b/PySONIC/neurons/stn.py @@ -1,590 +1,595 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2018-11-29 16:56:45 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-03-13 14:31:51 +# @Last Modified time: 2019-03-13 15:41:26 import numpy as np from scipy.optimize import brentq from ..core import PointNeuron from ..constants import FARADAY, Z_Ca from ..utils import nernst class OtsukaSTN(PointNeuron): ''' Class defining the Otsuka model of sub-thalamic nucleus neuron with 5 different current types: - Inward Sodium current (iNa) - Outward, delayed-rectifer Potassium current (iKd) - Inward, A-type Potassium current (iA) - Inward, low-threshold Calcium current (iCaT) - Inward, high-threshold Calcium current (iCaL) - Outward, Calcium-dependent Potassium current (iKCa) - Non-specific leakage current (iLeak) 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.* ''' name = 'STN' # Resting parameters Cm0 = 1e-2 # Cell membrane resting capacitance (F/m2) Vm0 = -58.0 # Resting membrane potential (mV) Cai0 = 5e-9 # M (5 nM) # Reversal potentials VNa = 60.0 # Sodium Nernst potential (mV) VK = -90.0 # Potassium Nernst potential (mV) # Physical constants T = 306.15 # K (33°C) # Calcium dynamics CCa_out = 2e-3 # M (2 mM) KCa = 2e3 # s-1 # Leakage current GLeak = 3.5 # Conductance of non-specific leakage current (S/m^2) VLeak = -60.0 # Leakage reversal potential (mV) # Fast Na current GNaMax = 490.0 # Max. conductance of Sodium current (S/m^2) thetax_m = -40 # mV thetax_h = -45.5 # mV kx_m = -8 # mV kx_h = 6.4 # mV tau0_m = 0.2 * 1e-3 # s tau1_m = 3 * 1e-3 # s tau0_h = 0 * 1e-3 # s tau1_h = 24.5 * 1e-3 # s thetaT_m = -53 # mV thetaT1_h = -50 # mV thetaT2_h = -50 # mV sigmaT_m = -0.7 # mV sigmaT1_h = -15 # mV sigmaT2_h = 16 # mV # Delayed rectifier K+ current GKMax = 570.0 # Max. conductance of delayed-rectifier Potassium current (S/m^2) thetax_n = -41 # mV kx_n = -14 # mV tau0_n = 0 * 1e-3 # s tau1_n = 11 * 1e-3 # s thetaT1_n = -40 # mV thetaT2_n = -40 # mV sigmaT1_n = -40 # mV sigmaT2_n = 50 # mV # T-type Ca2+ current GTMax = 50.0 # Max. conductance of low-threshold Calcium current (S/m^2) thetax_p = -56 # mV thetax_q = -85 # mV kx_p = -6.7 # mV kx_q = 5.8 # mV tau0_p = 5 * 1e-3 # s tau1_p = 0.33 * 1e-3 # s tau0_q = 0 * 1e-3 # s tau1_q = 400 * 1e-3 # s thetaT1_p = -27 # mV thetaT2_p = -102 # mV thetaT1_q = -50 # mV thetaT2_q = -50 # mV sigmaT1_p = -10 # mV sigmaT2_p = 15 # mV sigmaT1_q = -15 # mV sigmaT2_q = 16 # mV # L-type Ca2+ current GLMax = 150.0 # Max. conductance of high-threshold Calcium current (S/m^2) thetax_c = -30.6 # mV thetax_d1 = -60 # mV thetax_d2 = 0.1 * 1e-6 # M kx_c = -5 # mV kx_d1 = 7.5 # mV kx_d2 = 0.02 * 1e-6 # M tau0_c = 45 * 1e-3 # s tau1_c = 10 * 1e-3 # s tau0_d1 = 400 * 1e-3 # s tau1_d1 = 500 * 1e-3 # s tau_d2 = 130 * 1e-3 # s thetaT1_c = -27 # mV thetaT2_c = -50 # mV thetaT1_d1 = -40 # mV thetaT2_d1 = -20 # mV sigmaT1_c = -20 # mV sigmaT2_c = 15 # mV sigmaT1_d1 = -15 # mV sigmaT2_d1 = 20 # mV # A-type K+ current GAMax = 50.0 # Max. conductance of A-type Potassium current (S/m^2) thetax_a = -45 # mV thetax_b = -90 # mV kx_a = -14.7 # mV kx_b = 7.5 # mV tau0_a = 1 * 1e-3 # s tau1_a = 1 * 1e-3 # s tau0_b = 0 * 1e-3 # s tau1_b = 200 * 1e-3 # s thetaT_a = -40 # mV thetaT1_b = -60 # mV thetaT2_b = -40 # mV sigmaT_a = -0.5 # mV sigmaT1_b = -30 # mV sigmaT2_b = 10 # mV # Ca2+-activated K+ current GKCaMax = 10.0 # Max. conductance of Calcium-dependent Potassium current (S/m^2) thetax_r = 0.17 * 1e-6 # M kx_r = -0.08 * 1e-6 # M tau_r = 2 * 1e-3 # s # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_{Kd}\ kin.': ['n'], 'i_A\ kin.': ['a', 'b'], 'i_{CaT}\ kin.': ['p', 'q'], 'i_{CaL}\ kin.': ['c', 'd1', 'd2'], 'Ca^{2+}_i': ['Cai'], 'i_{KCa}\ kin.': ['r'], } def __init__(self): ''' Constructor of the class ''' # Names and initial states of the channels state probabilities self.states_names = ['a', 'b', 'c', 'd1', 'd2', 'm', 'h', 'n', 'p', 'q', 'r', 'Cai'] # Names of the different coefficients to be averaged in a lookup table. self.coeff_names = [ 'alphaa', 'betaa', 'alphab', 'betab', 'alphac', 'betac', 'alphad1', 'betad1', 'alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphap', 'betap', 'alphaq', 'betaq', ] - # Compute deff for that reversal potential - iCaT = self.iCaT(self.pinf(self.Vm0), self.qinf(self.Vm0), self.Vm0, self.Cai0) # mA/m2 - iCaL = self.iCaL(self.cinf(self.Vm0), self.d1inf(self.Vm0), self.d2inf(self.Cai0), - self.Vm0, self.Cai0) # mA/m2 - self.deff = -(iCaT + iCaL) / (Z_Ca * FARADAY * self.KCa * self.Cai0) * 1e-6 # m + # Compute deff for neuron's resting potential and equilibrium Calcium concentration + self.deff = self.getEffectiveDepth(self.Cai0, self.Vm0) # m # Compute conversion factor from electrical current (mA/m2) to Calcium concentration (M/s) - self.i2CCa = 1e-6 / (Z_Ca * self.deff * FARADAY) + self.i2Cai = 1e-6 / (Z_Ca * self.deff * FARADAY) # Initial states self.states0 = self.steadyStates(self.Vm0) - # Charge interval bounds for lookup creation - self.Qbounds = np.array([np.round(self.Vm0 - 25.0), 50.0]) * self.Cm0 * 1e-3 # C/m2 + + 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 * self.KCa * Cai) * 1e-6 # m def _xinf(self, var, theta, k): ''' Generic function computing the steady-state activation/inactivation of a particular ion channel at a given voltage or ion concentration. :param var: membrane potential (mV) or ion concentration (mM) :param theta: half-(in)activation voltage or concentration (mV or mM) :param k: slope parameter of (in)activation function (mV or mM) :return: steady-state (in)activation (-) ''' return 1 / (1 + np.exp((var - theta) / k)) def ainf(self, Vm): return self._xinf(Vm, self.thetax_a, self.kx_a) def binf(self, Vm): return self._xinf(Vm, self.thetax_b, self.kx_b) def cinf(self, Vm): return self._xinf(Vm, self.thetax_c, self.kx_c) def d1inf(self, Vm): return self._xinf(Vm, self.thetax_d1, self.kx_d1) def d2inf(self, Cai): return self._xinf(Cai, self.thetax_d2, self.kx_d2) def minf(self, Vm): return self._xinf(Vm, self.thetax_m, self.kx_m) def hinf(self, Vm): return self._xinf(Vm, self.thetax_h, self.kx_h) def ninf(self, Vm): return self._xinf(Vm, self.thetax_n, self.kx_n) def pinf(self, Vm): return self._xinf(Vm, self.thetax_p, self.kx_p) def qinf(self, Vm): return self._xinf(Vm, self.thetax_q, self.kx_q) def rinf(self, Cai): return self._xinf(Cai, self.thetax_r, self.kx_r) def _taux1(self, Vm, theta, sigma, tau0, tau1): ''' Generic function computing the voltage-dependent, activation/inactivation time constant of a particular ion channel at a given voltage (first variant). :param Vm: membrane potential (mV) :param theta: voltage at which (in)activation time constant is half-maximal (mV) :param sigma: slope parameter of (in)activation time constant function (mV) :param tau0: minimal time constant (s) :param tau1: modulated time constant (s) :return: (in)activation time constant (s) ''' return tau0 + tau1 / (1 + np.exp(-(Vm - theta) / sigma)) def taua(self, Vm): return self._taux1(Vm, self.thetaT_a, self.sigmaT_a, self.tau0_a, self.tau1_a) def taum(self, Vm): return self._taux1(Vm, self.thetaT_m, self.sigmaT_m, self.tau0_m, self.tau1_m) def _taux2(self, Vm, theta1, theta2, sigma1, sigma2, 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) def derA(self, Vm, a): return (self.ainf(Vm) - a) / self.taua(Vm) def derB(self, Vm, b): return (self.binf(Vm) - b) / self.taub(Vm) def derC(self, Vm, c): return (self.cinf(Vm) - c) / self.tauc(Vm) def derD1(self, Vm, d1): return (self.d1inf(Vm) - d1) / self.taud1(Vm) def derD2(self, Cai, d2): return (self.d2inf(Cai) - d2) / self.tau_d2 def derM(self, Vm, m): return (self.minf(Vm) - m) / self.taum(Vm) def derH(self, Vm, h): return (self.hinf(Vm) - h) / self.tauh(Vm) def derN(self, Vm, n): return (self.ninf(Vm) - n) / self.taun(Vm) def derP(self, Vm, p): return (self.pinf(Vm) - p) / self.taup(Vm) def derQ(self, Vm, q): return (self.qinf(Vm) - q) / self.tauq(Vm) def derR(self, Cai, r): return (self.rinf(Cai) - r) / self.tau_r def derCai(self, p, q, c, d1, d2, Cai, Vm): ''' Compute the time derivative of the Calcium concentration in submembranal space. :param Vm: membrane potential (mV) :param Cai: Calcium concentration in submembranal space (M) :param p: open-probability of p-gate :param q: open-probability of q-gate :param c: open-probability of c-gate :param d1: open-probability of d1-gate :param d2: open-probability of d2-gate :return: time derivative of Calcium concentration in submembranal space (M/s) ''' iCaT = self.iCaT(p, q, Vm, Cai) iCaL = self.iCaL(c, d1, d2, Vm, Cai) - return - self.i2CCa * (iCaT + iCaL) - Cai * self.KCa + return - self.i2Cai * (iCaT + iCaL) - Cai * self.KCa def derCaiSteadyState(self, Cai, Vm): ''' Compute the time derivative of the Calcium concentration in submembranal space, assuming quasi-steady gating states. :param Vm: membrane potential (mV) :param Cai: Calcium concentration in submembranal space (M) :return: time derivative of Calcium concentration in submembranal space (M/s) ''' return self.derCai(self.pinf(Vm), self.qinf(Vm), self.cinf(Vm), self.d1inf(Vm), self.d2inf(Cai), Cai, Vm) def findCaiSteadyState(self, Vm): ''' Find the steady-state intracellular Calcium concentration for a specific membrane potential. :param Vm: membrane potential (mV) :return: steady-state Calcium concentration in submembranal space (M) ''' return brentq( lambda x: self.derCaiSteadyState(x, Vm), self.Cai0 * 1e-4, self.Cai0 * 1e3, xtol=1e-16 ) def iNa(self, m, h, Vm): ''' Sodium current :param m: open-probability of m-gate :param h: open-probability of h-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GNaMax * m**3 * h * (Vm - self.VNa) def iKd(self, n, Vm): ''' Delayed-rectifier Potassium current :param n: open-probability of n-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GKMax * n**4 * (Vm - self.VK) def iA(self, a, b, Vm): ''' A-type Potassium current :param a: open-probability of a-gate :param b: open-probability of b-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GAMax * a**2 * b * (Vm - self.VK) def iCaT(self, p, q, Vm, Cai): ''' Low-threshold (T-type) Calcium current :param p: open-probability of p-gate :param q: open-probability of q-gate :param Vm: membrane potential (mV) :param Cai: submembrane Calcium concentration (M) :return: current per unit area (mA/m2) ''' return self.GTMax * p**2 * q * (Vm - nernst(Z_Ca, Cai, self.CCa_out, self.T)) def iCaL(self, c, d1, d2, Vm, Cai): ''' High-threshold (L-type) Calcium current :param c: open-probability of c-gate :param d1: open-probability of d1-gate :param d2: open-probability of d2-gate :param Vm: membrane potential (mV) :param Cai: submembrane Calcium concentration (M) :return: current per unit area (mA/m2) ''' return self.GLMax * c**2 * d1 * d2 * (Vm - nernst(Z_Ca, Cai, self.CCa_out, self.T)) def iKCa(self, r, Vm): ''' Calcium-activated Potassium current :param r: open-probability of r-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GKCaMax * r**2 * (Vm - self.VK) def iLeak(self, Vm): ''' Non-specific leakage current :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GLeak * (Vm - self.VLeak) def currents(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' a, b, c, d1, d2, m, h, n, p, q, r, Cai = states return { 'iNa': self.iNa(m, h, Vm), 'iKd': self.iKd(n, Vm), 'iA': self.iA(a, b, Vm), 'iCaT': self.iCaT(p, q, Vm, Cai), 'iCaL': self.iCaL(c, d1, d2, Vm, Cai), 'iKCa': self.iKCa(r, Vm), 'iLeak': self.iLeak(Vm) } # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Solve the equation dx/dt = 0 at Vm for each x-state aeq = self.ainf(Vm) beq = self.binf(Vm) ceq = self.cinf(Vm) d1eq = self.d1inf(Vm) meq = self.minf(Vm) heq = self.hinf(Vm) neq = self.ninf(Vm) peq = self.pinf(Vm) qeq = self.qinf(Vm) # Cai_eq = self.Cai0 Cai_eq = self.findCaiSteadyState(Vm) d2eq = self.d2inf(Cai_eq) req = self.rinf(Cai_eq) return np.array([aeq, beq, ceq, d1eq, d2eq, meq, heq, neq, peq, qeq, req, Cai_eq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' a, b, c, d1, d2, m, h, n, p, q, r, Cai = states dadt = self.derA(Vm, a) dbdt = self.derB(Vm, b) dcdt = self.derC(Vm, c) dd1dt = self.derD1(Vm, d1) dd2dt = self.derD2(Cai, d2) dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dpdt = self.derP(Vm, p) dqdt = self.derQ(Vm, q) drdt = self.derR(Cai, r) dCaidt = self.derCai(p, q, c, d1, d2, Cai, Vm) return [dadt, dbdt, dcdt, dd1dt, dd2dt, dmdt, dhdt, dndt, dpdt, dqdt, drdt, dCaidt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants Ta = self.taua(Vm) alphaa_avg = np.mean(self.ainf(Vm) / Ta) betaa_avg = np.mean(1 / Ta) - alphaa_avg Tb = self.taub(Vm) alphab_avg = np.mean(self.binf(Vm) / Tb) betab_avg = np.mean(1 / Tb) - alphab_avg Tc = self.tauc(Vm) alphac_avg = np.mean(self.cinf(Vm) / Tc) betac_avg = np.mean(1 / Tc) - alphac_avg Td1 = self.taud1(Vm) alphad1_avg = np.mean(self.d1inf(Vm) / Td1) betad1_avg = np.mean(1 / Td1) - alphad1_avg Tm = self.taum(Vm) alpham_avg = np.mean(self.minf(Vm) / Tm) betam_avg = np.mean(1 / Tm) - alpham_avg Th = self.tauh(Vm) alphah_avg = np.mean(self.hinf(Vm) / Th) betah_avg = np.mean(1 / Th) - alphah_avg Tn = self.taun(Vm) alphan_avg = np.mean(self.ninf(Vm) / Tn) betan_avg = np.mean(1 / Tn) - alphan_avg Tp = self.taup(Vm) alphap_avg = np.mean(self.pinf(Vm) / Tp) betap_avg = np.mean(1 / Tp) - alphap_avg Tq = self.tauq(Vm) alphaq_avg = np.mean(self.qinf(Vm) / Tq) betaq_avg = np.mean(1 / Tq) - alphaq_avg # Return array of coefficients return np.array([ alphaa_avg, betaa_avg, alphab_avg, betab_avg, alphac_avg, betac_avg, alphad1_avg, betad1_avg, alpham_avg, betam_avg, alphah_avg, betah_avg, alphan_avg, betan_avg, alphap_avg, betap_avg, alphaq_avg, betaq_avg ]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) Vmeff = np.interp(Qm, interp_data['Q'], interp_data['V']) a, b, c, d1, d2, m, h, n, p, q, r, Cai = states dadt = rates[0] * (1 - a) - rates[1] * a dbdt = rates[2] * (1 - b) - rates[3] * b dcdt = rates[4] * (1 - c) - rates[5] * c dd1dt = rates[6] * (1 - d1) - rates[7] * d1 dd2dt = self.derD2(Cai, d2) dmdt = rates[8] * (1 - m) - rates[9] * m dhdt = rates[10] * (1 - h) - rates[11] * h dndt = rates[12] * (1 - n) - rates[13] * n dpdt = rates[14] * (1 - p) - rates[15] * p dqdt = rates[16] * (1 - q) - rates[17] * q drdt = self.derR(Cai, r) dCaidt = self.derCai(p, q, c, d1, d2, Cai, Vmeff) return [dadt, dbdt, dcdt, dd1dt, dd2dt, dmdt, dhdt, dndt, dpdt, dqdt, drdt, dCaidt] diff --git a/PySONIC/neurons/thalamic.py b/PySONIC/neurons/thalamic.py index 08e9239..1263c08 100644 --- a/PySONIC/neurons/thalamic.py +++ b/PySONIC/neurons/thalamic.py @@ -1,789 +1,784 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-07-31 15:20:54 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-03-13 15:08:46 +# @Last Modified time: 2019-03-13 15:41:33 import numpy as np from ..core import PointNeuron from ..constants import FARADAY, Z_Ca from ..utils import vtrap class Thalamic(PointNeuron): ''' Class defining the generic membrane channel dynamics of a thalamic neuron with 4 different current types: - Inward Sodium current - Outward Potassium current - Inward Calcium current - Non-specific leakage current This generic class cannot be used directly as it does not contain any specific parameters. Reference: *Plaksin, M., Kimmel, E., and Shoham, S. (2016). Cell-Type-Selective Effects of Intramembrane Cavitation as a Unifying Theoretical Framework for Ultrasonic Neuromodulation. eNeuro 3.* ''' # Generic biophysical parameters of thalamic cells Cm0 = 1e-2 # Cell membrane resting capacitance (F/m2) Vm0 = 0.0 # Dummy value for membrane potential (mV) VNa = 50.0 # Sodium Nernst potential (mV) VK = -90.0 # Potassium Nernst potential (mV) VCa = 120.0 # Calcium Nernst potential (mV) def __init__(self): ''' Constructor of the class ''' # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'u'] self.states0 = np.array([]) # Names of the different coefficients to be averaged in a lookup table. self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas', 'alphau', 'betau'] - # Charge interval bounds for lookup creation - self.Qbounds = np.array([np.round(self.Vm0 - 25.0), 50.0]) * self.Cm0 * 1e-3 # C/m2 - def alpham(self, Vm): ''' Compute the alpha rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = 0.32 * vtrap(13 - Vdiff, 4) # ms-1 return alpha * 1e3 # s-1 def betam(self, Vm): ''' Compute the beta rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = 0.28 * vtrap(Vdiff - 40, 5) # ms-1 return beta * 1e3 # s-1 def alphah(self, Vm): ''' Compute the alpha rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = (0.128 * np.exp(-(Vdiff - 17) / 18)) # ms-1 return alpha * 1e3 # s-1 def betah(self, Vm): ''' Compute the beta rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = (4 / (1 + np.exp(-(Vdiff - 40) / 5))) # ms-1 return beta * 1e3 # s-1 def alphan(self, Vm): ''' Compute the alpha rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = 0.032 * vtrap(15 - Vdiff, 5) # ms-1 return alpha * 1e3 # s-1 def betan(self, Vm): ''' Compute the beta rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = (0.5 * np.exp(-(Vdiff - 10) / 40)) # ms-1 return beta * 1e3 # s-1 def derM(self, Vm, m): ''' Compute the evolution of the open-probability of Sodium channels. :param Vm: membrane potential (mV) :param m: open-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alpham(Vm) * (1 - m) - self.betam(Vm) * m def derH(self, Vm, h): ''' Compute the evolution of the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :param h: inactivation-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphah(Vm) * (1 - h) - self.betah(Vm) * h def derN(self, Vm, n): ''' Compute the evolution of the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :param n: open-probability of delayed-rectifier Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphan(Vm) * (1 - n) - self.betan(Vm) * n def derS(self, Vm, s): ''' Compute the evolution of the open-probability of the S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :param s: open-probability of S-type Calcium activation gates (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.sinf(Vm) - s) / self.taus(Vm) def derU(self, Vm, u): ''' Compute the evolution of the open-probability of the U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :param u: open-probability of U-type Calcium inactivation gates (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.uinf(Vm) - u) / self.tauu(Vm) def iNa(self, m, h, Vm): ''' Sodium current :param m: open-probability of m-gate :param h: open-probability of h-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GNa = self.GNaMax * m**3 * h return GNa * (Vm - self.VNa) def iKd(self, n, Vm): ''' Delayed-rectifier Potassium current :param n: open-probability of n-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GKMax * n**4 return GK * (Vm - self.VK) def iCaT(self, s, u, Vm): ''' Low-threshold (Ts-type) Calcium current :param s: open-probability of s-gate :param u: open-probability of u-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GT = self.GTMax * s**2 * u return GT * (Vm - self.VCa) def iLeak(self, Vm): ''' Non-specific leakage current :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GLeak * (Vm - self.VLeak) def currents(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, u = states return { 'iNa': self.iNa(m, h, Vm), 'iKd': self.iKd(n, Vm), 'iCaT': self.iCaT(s, u, Vm), 'iLeak': self.iLeak(Vm) } # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) seq = self.sinf(Vm) ueq = self.uinf(Vm) return np.array([meq, heq, neq, seq, ueq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, u = states dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) dudt = self.derU(Vm, u) return [dmdt, dhdt, dndt, dsdt, dudt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) Ts = self.taus(Vm) as_avg = np.mean(self.sinf(Vm) / Ts) bs_avg = np.mean(1 / Ts) - as_avg Tu = np.array([self.tauu(v) for v in Vm]) au_avg = np.mean(self.uinf(Vm) / Tu) bu_avg = np.mean(1 / Tu) - au_avg # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg, au_avg, bu_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) m, h, n, s, u = states dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s dudt = rates[8] * (1 - u) - rates[9] * u return [dmdt, dhdt, dndt, dsdt, dudt] class ThalamicRE(Thalamic): ''' Specific membrane channel dynamics of a thalamic reticular neuron. References: *Destexhe, A., Contreras, D., Steriade, M., Sejnowski, T.J., and Huguenard, J.R. (1996). In vivo, in vitro, and computational analysis of dendritic calcium currents in thalamic reticular neurons. J. Neurosci. 16, 169–185.* *Huguenard, J.R., and Prince, D.A. (1992). A novel T-type current underlies prolonged Ca(2+)-dependent burst firing in GABAergic neurons of rat thalamic reticular nucleus. J. Neurosci. 12, 3804–3817.* ''' # Name of channel mechanism name = 'RE' # Cell-specific biophysical parameters Vm0 = -89.5 # Cell membrane resting potential (mV) GNaMax = 2000.0 # Max. conductance of Sodium current (S/m^2) GKMax = 200.0 # Max. conductance of Potassium current (S/m^2) GTMax = 30.0 # Max. conductance of low-threshold Calcium current (S/m^2) GLeak = 0.5 # Conductance of non-specific leakage current (S/m^2) VLeak = -90.0 # Non-specific leakage Nernst potential (mV) VT = -67.0 # Spike threshold adjustment parameter (mV) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_{Kd}\ kin.': ['n'], 'i_{TS}\ kin.': ['s', 'u'] } def __init__(self): ''' Constructor of the class. ''' # Instantiate parent class super().__init__() # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) def sinf(self, Vm): ''' Compute the asymptotic value of the open-probability of the S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1.0 + np.exp(-(Vm + 52.0) / 7.4)) # prob def taus(self, Vm): ''' Compute the decay time constant for adaptation of S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' return (1 + 0.33 / (np.exp((Vm + 27.0) / 10.0) + np.exp(-(Vm + 102.0) / 15.0))) * 1e-3 # s def uinf(self, Vm): ''' Compute the asymptotic value of the open-probability of the U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1.0 + np.exp((Vm + 80.0) / 5.0)) # prob def tauu(self, Vm): ''' Compute the decay time constant for adaptation of U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' return (28.3 + 0.33 / (np.exp((Vm + 48.0) / 4.0) + np.exp(-(Vm + 407.0) / 50.0))) * 1e-3 # s class ThalamoCortical(Thalamic): ''' Specific membrane channel dynamics of a thalamo-cortical neuron, with a specific hyperpolarization-activated, mixed cationic current and a leakage Potassium current. References: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* *Destexhe, A., Bal, T., McCormick, D.A., and Sejnowski, T.J. (1996). Ionic mechanisms underlying synchronized oscillations and propagating waves in a model of ferret thalamic slices. J. Neurophysiol. 76, 2049–2070.* *McCormick, D.A., and Huguenard, J.R. (1992). A model of the electrophysiological properties of thalamocortical relay neurons. J. Neurophysiol. 68, 1384–1400.* ''' # Name of channel mechanism name = 'TC' # Cell-specific biophysical parameters # Vm0 = -63.4 # Cell membrane resting potential (mV) Vm0 = -61.93 # Cell membrane resting potential (mV) GNaMax = 900.0 # Max. conductance of Sodium current (S/m^2) GKMax = 100.0 # Max. conductance of Potassium current (S/m^2) GTMax = 20.0 # Max. conductance of low-threshold Calcium current (S/m^2) GKL = 0.138 # Conductance of leakage Potassium current (S/m^2) GhMax = 0.175 # Max. conductance of mixed cationic current (S/m^2) GLeak = 0.1 # Conductance of non-specific leakage current (S/m^2) Vh = -40.0 # Mixed cationic current reversal potential (mV) VLeak = -70.0 # Non-specific leakage Nernst potential (mV) VT = -52.0 # Spike threshold adjustment parameter (mV) Vx = 0.0 # Voltage-dependence uniform shift factor at 36°C (mV) - tau_Ca_removal = 5e-3 # decay time constant for intracellular Ca2+ dissolution (s) - CCa_min = 50e-9 # minimal intracellular Calcium concentration (M) + 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) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_{Kd}\ kin.': ['n'], 'i_{T}\ kin.': ['s', 'u'], 'i_{H}\ kin.': ['O', 'OL'] } def __init__(self): ''' Constructor of the class. ''' # Instantiate parent class super().__init__() # Compute current to concentration conversion constant - self.iT_2_CCa = 1e-6 / (self.deff * Z_Ca * FARADAY) + self.i2Cai = 1e-6 / (self.deff * Z_Ca * FARADAY) # Define names of the channels state probabilities - self.states_names += ['O', 'C', 'P0', 'C_Ca'] + self.states_names += ['O', 'C', 'P0', 'Cai'] # Define the names of the different coefficients to be averaged in a lookup table. self.coeff_names += ['alphao', 'betao'] # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) def sinf(self, Vm): ''' Compute the asymptotic value of the open-probability of the S-type, activation gate of Calcium channels. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1.0 + np.exp(-(Vm + self.Vx + 57.0) / 6.2)) # prob def taus(self, Vm): ''' Compute the decay time constant for adaptation of S-type, activation gate of Calcium channels. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' tmp = np.exp(-(Vm + self.Vx + 132.0) / 16.7) + np.exp((Vm + self.Vx + 16.8) / 18.2) return 1.0 / 3.7 * (0.612 + 1.0 / tmp) * 1e-3 # s def uinf(self, Vm): ''' Compute the asymptotic value of the open-probability of the U-type, inactivation gate of Calcium channels. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1.0 + np.exp((Vm + self.Vx + 81.0) / 4.0)) # prob def tauu(self, Vm): ''' Compute the decay time constant for adaptation of U-type, inactivation gate of Calcium channels. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' if Vm + self.Vx < -80.0: return 1.0 / 3.7 * np.exp((Vm + self.Vx + 467.0) / 66.6) * 1e-3 # s else: return 1 / 3.7 * (np.exp(-(Vm + self.Vx + 22) / 10.5) + 28.0) * 1e-3 # s def derS(self, Vm, s): ''' Compute the evolution of the open-probability of the S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :param s: open-probability of S-type Calcium activation gates (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.sinf(Vm) - s) / self.taus(Vm) def derU(self, Vm, u): ''' Compute the evolution of the open-probability of the U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :param u: open-probability of U-type Calcium inactivation gates (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.uinf(Vm) - u) / self.tauu(Vm) def oinf(self, Vm): ''' Voltage-dependent steady-state activation of hyperpolarization-activated cation current channels. Reference: *Huguenard, J.R., and McCormick, D.A. (1992). Simulation of the currents involved in rhythmic oscillations in thalamic relay neurons. J. Neurophysiol. 68, 1373–1383.* :param Vm: membrane potential (mV) :return: steady-state activation (-) ''' return 1.0 / (1.0 + np.exp((Vm + 75.0) / 5.5)) def tauo(self, Vm): ''' Time constant for activation of hyperpolarization-activated cation current channels. Reference: *Huguenard, J.R., and McCormick, D.A. (1992). Simulation of the currents involved in rhythmic oscillations in thalamic relay neurons. J. Neurophysiol. 68, 1373–1383.* :param Vm: membrane potential (mV) :return: time constant (s) ''' return 1 / (np.exp(-14.59 - 0.086 * Vm) + np.exp(-1.87 + 0.0701 * Vm)) * 1e-3 def alphao(self, Vm): ''' Transition rate between closed and open form of hyperpolarization-activated cation current channels. :param Vm: membrane potential (mV) :return: transition rate (s-1) ''' return self.oinf(Vm) / self.tauo(Vm) def betao(self, Vm): ''' Transition rate between open and closed form of hyperpolarization-activated cation current channels. :param Vm: membrane potential (mV) :return: transition rate (s-1) ''' return (1 - self.oinf(Vm)) / self.tauo(Vm) def derC(self, C, O, Vm): ''' Compute the evolution of the proportion of hyperpolarization-activated cation current channels in closed state. Kinetics scheme of Calcium dependent activation derived from: *Destexhe, A., Bal, T., McCormick, D.A., and Sejnowski, T.J. (1996). Ionic mechanisms underlying synchronized oscillations and propagating waves in a model of ferret thalamic slices. J. Neurophysiol. 76, 2049–2070.* :param Vm: membrane potential (mV) :param C: proportion of Ih channels in closed state (-) :param O: proportion of Ih channels in open state (-) :return: derivative of proportion w.r.t. time (s-1) ''' return self.betao(Vm) * O - self.alphao(Vm) * C def derO(self, C, O, P0, Vm): ''' Compute the evolution of the proportion of hyperpolarization-activated cation current channels in open state. Kinetics scheme of Calcium dependent activation derived from: *Destexhe, A., Bal, T., McCormick, D.A., and Sejnowski, T.J. (1996). Ionic mechanisms underlying synchronized oscillations and propagating waves in a model of ferret thalamic slices. J. Neurophysiol. 76, 2049–2070.* :param Vm: membrane potential (mV) :param C: proportion of Ih channels in closed state (-) :param O: proportion of Ih channels in open state (-) :param P0: proportion of Ih channels regulating factor in unbound state (-) :return: derivative of proportion w.r.t. time (s-1) ''' return - self.derC(C, O, Vm) - self.k3 * O * (1 - P0) + self.k4 * (1 - O - C) - def derP0(self, P0, C_Ca): + def derP0(self, P0, Cai): ''' Compute the evolution of the proportion of Ih channels regulating factor in unbound state. Kinetics scheme of Calcium dependent activation derived from: *Destexhe, A., Bal, T., McCormick, D.A., and Sejnowski, T.J. (1996). Ionic mechanisms underlying synchronized oscillations and propagating waves in a model of ferret thalamic slices. J. Neurophysiol. 76, 2049–2070.* :param Vm: membrane potential (mV) :param P0: proportion of Ih channels regulating factor in unbound state (-) - :param C_Ca: Calcium concentration in effective submembranal space (M) + :param Cai: Calcium concentration in effective submembranal space (M) :return: derivative of proportion w.r.t. time (s-1) ''' - return self.k2 * (1 - P0) - self.k1 * P0 * C_Ca**self.nCa + return self.k2 * (1 - P0) - self.k1 * P0 * Cai**self.nCa - def derC_Ca(self, C_Ca, ICa): + def derCai(self, Cai, s, u, Vm): ''' Compute the evolution of the Calcium concentration in submembranal space. 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.* - + :param Cai: Calcium concentration in submembranal space (M) + :param s: open-probability of s-gate + :param u: open-probability of u-gate :param Vm: membrane potential (mV) - :param C_Ca: Calcium concentration in submembranal space (M) - :param ICa: inward Calcium current filling up the submembranal space with Ca2+ (mA/m2) - :return: derivative of Calcium concentration in submembranal space w.r.t. time (s-1) + :return: time derivative of Calcium concentration in submembranal space (M/s) ''' - return (self.CCa_min - C_Ca) / self.tau_Ca_removal - self.iT_2_CCa * ICa + return (self.Cai_min - Cai) / self.taur_Cai - self.i2Cai * self.iCaT(s, u, Vm) def iKLeak(self, Vm): ''' Potassium leakage current :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GKL * (Vm - self.VK) def iH(self, O, C, Vm): ''' Outward mixed cationic current :param O: proportion of the channels in open form :param C: proportion of the channels in closed form :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' OL = 1 - O - C # proportion of channels in locked-open form return self.GhMax * (O + 2 * OL) * (Vm - self.Vh) def currents(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, u, O, C, _, _ = states currents = super().currents(Vm, [m, h, n, s, u]) currents['iKLeak'] = self.iKLeak(Vm) # mA/m2 currents['iH'] = self.iH(O, C, Vm) # mA/m2 return currents def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Call parent method to compute Sodium, Potassium and Calcium channels gates steady-states NaKCa_eqstates = super().steadyStates(Vm) # Compute steady-state Calcium current seq = NaKCa_eqstates[3] ueq = NaKCa_eqstates[4] - iTeq = self.iCaT(seq, ueq, Vm) + iCaTeq = self.iCaT(seq, ueq, Vm) # Compute steady-state variables for the kinetics system of Ih - CCa_eq = self.CCa_min - self.tau_Ca_removal * self.iT_2_CCa * iTeq - P0_eq = self.k2 / (self.k2 + self.k1 * CCa_eq**self.nCa) + Cai_eq = self.Cai_min - self.taur_Cai * self.i2Cai * iCaTeq + P0_eq = self.k2 / (self.k2 + self.k1 * Cai_eq**self.nCa) BA = self.betao(Vm) / self.alphao(Vm) O_eq = self.k4 / (self.k3 * (1 - P0_eq) + self.k4 * (1 + BA)) C_eq = BA * O_eq - kin_eqstates = np.array([O_eq, C_eq, P0_eq, CCa_eq]) + kin_eqstates = np.array([O_eq, C_eq, P0_eq, Cai_eq]) # Merge all steady-states and return return np.concatenate((NaKCa_eqstates, kin_eqstates)) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' - m, h, n, s, u, O, C, P0, C_Ca = states + m, h, n, s, u, O, C, P0, Cai = states NaKCa_states = [m, h, n, s, u] NaKCa_derstates = super().derStates(Vm, NaKCa_states) - dO_dt = self.derO(C, O, P0, Vm) - dC_dt = self.derC(C, O, Vm) - dP0_dt = self.derP0(P0, C_Ca) - ICa = self.iCaT(s, u, Vm) - dCCa_dt = self.derC_Ca(C_Ca, ICa) + dOdt = self.derO(C, O, P0, Vm) + dCdt = self.derC(C, O, Vm) + dP0dt = self.derP0(P0, Cai) + dCaidt = self.derCai(Cai, s, u, Vm) - return NaKCa_derstates + [dO_dt, dC_dt, dP0_dt, dCCa_dt] + return NaKCa_derstates + [dOdt, dCdt, dP0dt, dCaidt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute effective coefficients for Sodium, Potassium and Calcium conductances NaKCa_effrates = super().getEffRates(Vm) # Compute effective coefficients for Ih conductance ao_avg = np.mean(self.alphao(Vm)) bo_avg = np.mean(self.betao(Vm)) iH_effrates = np.array([ao_avg, bo_avg]) # Return array of coefficients return np.concatenate((NaKCa_effrates, iH_effrates)) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) Vmeff = np.interp(Qm, interp_data['Q'], interp_data['V']) # Unpack states - m, h, n, s, u, O, C, P0, C_Ca = states + m, h, n, s, u, O, C, P0, Cai = states - # INa, IK, ICa effective states derivatives + # INa, IK, iCaT effective states derivatives dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s dudt = rates[8] * (1 - u) - rates[9] * u # Ih effective states derivatives - dC_dt = rates[11] * O - rates[10] * C - dO_dt = - dC_dt - self.k3 * O * (1 - P0) + self.k4 * (1 - O - C) - dP0_dt = self.derP0(P0, C_Ca) - ICa_eff = self.iCaT(s, u, Vmeff) - dCCa_dt = self.derC_Ca(C_Ca, ICa_eff) + dCdt = rates[11] * O - rates[10] * C + dOdt = - dCdt - self.k3 * O * (1 - P0) + self.k4 * (1 - O - C) + dP0dt = self.derP0(P0, Cai) + dCaidt = self.derCai(Cai, s, u, Vmeff) # Merge derivatives and return - return [dmdt, dhdt, dndt, dsdt, dudt, dO_dt, dC_dt, dP0_dt, dCCa_dt] + return [dmdt, dhdt, dndt, dsdt, dudt, dOdt, dCdt, dP0dt, dCaidt] diff --git a/scripts/run_lookups.py b/scripts/run_lookups.py index c96b02c..0f97648 100644 --- a/scripts/run_lookups.py +++ b/scripts/run_lookups.py @@ -1,170 +1,170 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-06-02 17:50:10 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-01-23 18:38:38 +# @Last Modified time: 2019-03-13 15:42:19 ''' Create lookup table for specific neuron. ''' import os import pickle import logging import numpy as np from argparse import ArgumentParser from PySONIC.utils import logger, getNeuronLookupsFile from PySONIC.batches import createQueue, runBatch from PySONIC.neurons import getNeuronsDict from PySONIC.core import NeuronalBilayerSonophore # Default parameters defaults = dict( neuron='RS', radius=np.array([16.0, 32.0, 64.0]), # nm freq=np.array([20., 100., 500., 1e3, 2e3, 3e3, 4e3]), # kHz amp=np.insert(np.logspace(np.log10(0.1), np.log10(600), num=50), 0, 0.0), # kPa ) def computeAStimLookups(neuron, aref, fref, Aref, Qref, phi=np.pi, mpi=False, loglevel=logging.INFO): ''' Run simulations of the mechanical system for a multiple combinations of imposed sonophore raduis, US frequencies, acoustic amplitudes and charge densities, compute effective coefficients and store them in a dictionary of 3D arrays. :param neuron: neuron object :param aref: array of sonophore radii (m) :param fref: array of acoustic drive frequencies (Hz) :param Aref: array of acoustic drive amplitudes (Pa) :param Qref: array of membrane charge densities (C/m2) :param phi: acoustic drive phase (rad) :param mpi: boolean statting wether or not to use multiprocessing :return: lookups dictionary ''' # Check validity of input parameters for key, values in {'radii': aref, 'frequencies': fref, 'amplitudes': Aref}.items(): if not (isinstance(values, list) or isinstance(values, np.ndarray)): raise TypeError('Invalid {} (must be provided as list or numpy array)'.format(key)) if not all(isinstance(x, float) for x in values): raise TypeError('Invalid {} (must all be float typed)'.format(key)) if len(values) == 0: raise ValueError('Empty {} array'.format(key)) if key in ('radii', 'frequencies') and min(values) <= 0: raise ValueError('Invalid {} (must all be strictly positive)'.format(key)) if key is 'amplitudes' and min(values) < 0: raise ValueError('Invalid {} (must all be positive or null)'.format(key)) # populate inputs dictionary inputs = dict( a=aref, # nm f=fref, # Hz A=Aref, # Pa Q=Qref # C/m2 ) # create simulation queue na, nf, nA, nQ = len(aref), len(fref), len(Aref), len(Qref) queue = createQueue((fref, Aref, Qref)) # run simulations and populate outputs (list of lists) logger.info('Starting simulation batch for %s neuron', neuron.name) outputs = [] for a in aref: nbls = NeuronalBilayerSonophore(a, neuron) outputs += runBatch(nbls, 'computeEffVars', queue, mpi=mpi, loglevel=loglevel) outputs = np.array(outputs).T # Split comp times and lookups tcomps = outputs[0] outputs = outputs[1:] # Reshape comp times into 4D array tcomps = tcomps.reshape(na, nf, nA, nQ) # reshape outputs into 4D arrays and add them to lookups dictionary logger.info('Reshaping output into lookup tables') keys = ['V', 'ng'] + neuron.coeff_names assert len(keys) == len(outputs), 'Lookup keys not matching array size' lookups = {} for key, output in zip(keys, outputs): lookups[key] = output.reshape(na, nf, nA, nQ) # Store inputs, lookup data and comp times in dictionary df = { 'input': inputs, 'lookup': lookups, 'tcomp': tcomps } return df def main(): ap = ArgumentParser() # Runtime options ap.add_argument('--mpi', default=False, action='store_true', help='Use multiprocessing') ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity') ap.add_argument('-t', '--test', default=False, action='store_true', help='Test configuration') # Stimulation parameters ap.add_argument('-n', '--neuron', type=str, default=defaults['neuron'], help='Neuron name (string)') ap.add_argument('-a', '--radius', nargs='+', type=float, help='Sonophore radius (nm)') ap.add_argument('-f', '--freq', nargs='+', type=float, help='US frequency (kHz)') ap.add_argument('-A', '--amp', nargs='+', type=float, help='Acoustic pressure amplitude (kPa)') ap.add_argument('-Q', '--charge', nargs='+', type=float, help='Membrane charge density (nC/cm2)') # Parse arguments args = {key: value for key, value in vars(ap.parse_args()).items() if value is not None} loglevel = logging.DEBUG if args['verbose'] is True else logging.INFO logger.setLevel(loglevel) mpi = args['mpi'] neuron_str = args['neuron'] radii = np.array(args.get('radius', defaults['radius'])) * 1e-9 # m freqs = np.array(args.get('freq', defaults['freq'])) * 1e3 # Hz amps = np.array(args.get('amp', defaults['amp'])) * 1e3 # Pa # Check neuron name validity if neuron_str not in getNeuronsDict(): logger.error('Unknown neuron type: "%s"', neuron_str) return neuron = getNeuronsDict()[neuron_str]() if 'charge' in args: charges = np.array(args['charge']) * 1e-5 # C/m2 else: - charges = np.arange(neuron.Qbounds[0], neuron.Qbounds[1] + 1e-5, 1e-5) # C/m2 + charges = np.arange(neuron.Qbounds()[0], neuron.Qbounds()[1] + 1e-5, 1e-5) # C/m2 if args['test']: radii = np.array([radii.min(), radii.max()]) freqs = np.array([freqs.min(), freqs.max()]) amps = np.array([amps.min(), amps.max()]) charges = np.array([charges.min(), 0., charges.max()]) # Check if lookup file already exists lookup_path = getNeuronLookupsFile(neuron.name) if os.path.isfile(lookup_path): logger.warning('"%s" file already exists and will be overwritten. ' + 'Continue? (y/n)', lookup_path) user_str = input() if user_str not in ['y', 'Y']: logger.error('%s Lookup creation canceled', neuron.name) return # compute lookups df = computeAStimLookups( neuron, radii, freqs, amps, charges, mpi=mpi, loglevel=loglevel) # Save dictionary in lookup file logger.info('Saving %s neuron lookup table in file: "%s"', neuron.name, lookup_path) with open(lookup_path, 'wb') as fh: pickle.dump(df, fh) if __name__ == '__main__': main() diff --git a/scripts/run_lookups_fs.py b/scripts/run_lookups_fs.py index fb12a90..29ea0d9 100644 --- a/scripts/run_lookups_fs.py +++ b/scripts/run_lookups_fs.py @@ -1,153 +1,153 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-06-02 17:50:10 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-01-23 18:12:24 +# @Last Modified time: 2019-03-13 15:42:29 ''' Create lookup table for specific neuron. ''' import os import pickle import logging import numpy as np from argparse import ArgumentParser from PySONIC.utils import logger, getNeuronLookupsFile from PySONIC.batches import createQueue, runBatch from PySONIC.neurons import getNeuronsDict from PySONIC.core import NeuronalBilayerSonophore # Default parameters defaults = dict( neuron='RS', radius=32.0, # nm freq=500, # kHz amp=np.insert(np.logspace(np.log10(0.1), np.log10(600), num=50), 0, 0.0), # kPa ) def computeAStimLookups(neuron, a, Fdrive, Aref, Qref, fsref, mpi=False, loglevel=logging.INFO): # Check validity of input parameters for key, values in {'coverage fractions': fsref, 'amplitudes': Aref}.items(): if not (isinstance(values, list) or isinstance(values, np.ndarray)): raise TypeError('Invalid {} (must be provided as list or numpy array)'.format(key)) if not all(isinstance(x, float) for x in values): raise TypeError('Invalid {} (must all be float typed)'.format(key)) if len(values) == 0: raise ValueError('Empty {} array'.format(key)) if key is 'coverage fractions' and min(values) <= 0: raise ValueError('Invalid {} (must all be strictly positive)'.format(key)) if key is 'amplitudes' and min(values) < 0: raise ValueError('Invalid {} (must all be positive or null)'.format(key)) # populate inputs dictionary inputs = dict( fs=fsref, # (-) A=Aref, # Pa Q=Qref # C/m2 ) # create simulation queue nA, nQ, nfs = len(Aref), len(Qref), len(fsref) queue = createQueue(([Fdrive], Aref, Qref, fsref)) # run simulations and populate outputs (list of lists) logger.info('Starting simulation batch for %s neuron', neuron.name) nbls = NeuronalBilayerSonophore(a, neuron) outputs = runBatch(nbls, 'computeEffVars', queue, mpi=mpi, loglevel=loglevel) outputs = np.array(outputs).T # Split comp times and lookups tcomps = outputs[0] outputs = outputs[1:] # Reshape comp times into 4D array tcomps = tcomps.reshape(nA, nQ, nfs) # reshape outputs into 4D arrays and add them to lookups dictionary logger.info('Reshaping output into lookup tables') keys = ['V', 'ng'] + neuron.coeff_names assert len(keys) == len(outputs), 'Lookup keys not matching array size' lookups = {} for key, output in zip(keys, outputs): lookups[key] = output.reshape(nA, nQ, nfs) # Store inputs, lookup data and comp times in dictionary df = { 'input': inputs, 'lookup': lookups, 'tcomp': tcomps } return df def main(): ap = ArgumentParser() # Runtime options ap.add_argument('--mpi', default=False, action='store_true', help='Use multiprocessing') ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity') ap.add_argument('-t', '--test', default=False, action='store_true', help='Test configuration') # Stimulation parameters ap.add_argument('-n', '--neuron', type=str, default=defaults['neuron'], help='Neuron name (string)') ap.add_argument('-a', '--radius', type=float, default=defaults['radius'], help='Sonophore radius (nm)') ap.add_argument('-f', '--freq', type=float, default=defaults['freq'], help='US frequency (kHz)') ap.add_argument('-A', '--amp', nargs='+', type=float, help='Acoustic pressure amplitude (kPa)') # Parse arguments args = {key: value for key, value in vars(ap.parse_args()).items() if value is not None} loglevel = logging.DEBUG if args['verbose'] is True else logging.INFO logger.setLevel(loglevel) mpi = args['mpi'] neuron_str = args['neuron'] a = args['radius'] * 1e-9 # m Fdrive = args['freq'] * 1e3 # Hz amps = np.array(args.get('amp', defaults['amp'])) * 1e3 # Pa fs = np.linspace(0, 100, 101) * 1e-2 # (-) # Check neuron name validity if neuron_str not in getNeuronsDict(): logger.error('Unknown neuron type: "%s"', neuron_str) return neuron = getNeuronsDict()[neuron_str]() - charges = np.arange(neuron.Qbounds[0], neuron.Qbounds[1] + 1e-5, 1e-5) # C/m2 + charges = np.arange(neuron.Qbounds()[0], neuron.Qbounds()[1] + 1e-5, 1e-5) # C/m2 if args['test']: fs = np.array([fs.min(), fs.max()]) amps = np.array([amps.min(), amps.max()]) charges = np.array([charges.min(), 0., charges.max()]) # Check if lookup file already exists lookup_path = getNeuronLookupsFile(neuron.name, a=a, Fdrive=Fdrive, fs=True) if os.path.isfile(lookup_path): logger.warning('"%s" file already exists and will be overwritten. ' + 'Continue? (y/n)', lookup_path) user_str = input() if user_str not in ['y', 'Y']: logger.error('%s Lookup creation canceled', neuron.name) return # compute lookups df = computeAStimLookups( neuron, a, Fdrive, amps, charges, fs, mpi=mpi, loglevel=loglevel) # Save dictionary in lookup file logger.info('Saving %s neuron lookup table in file: "%s"', neuron.name, lookup_path) with open(lookup_path, 'wb') as fh: pickle.dump(df, fh) if __name__ == '__main__': main()