diff --git a/PySONIC/neurons/stn.py b/PySONIC/neurons/stn.py index f4f3595..838c180 100644 --- a/PySONIC/neurons/stn.py +++ b/PySONIC/neurons/stn.py @@ -1,562 +1,562 @@ import numpy as np 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) CCa_in0 = 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 + 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': ['C_Ca'], 'i_{KCa}\ kin.': ['r'], 'I': ['iLeak', 'iNa', 'iKd', 'iA', 'iCaT2', 'iCaL2', 'iKCa', 'iNet'] } def __init__(self): ''' Constructor of the class ''' # Names and initial states of the channels state probabilities self.states_names = ['a', 'b', 'c', 'd1', 'd2', 'm', 'h', 'n', 'p', 'q', 'r', 'C_Ca'] # Names of the different coefficients to be averaged in a lookup table. self.coeff_names = [ 'alphaa', 'betaa', 'alphab', 'betab', 'alphac', 'betac', 'alphad1', 'betad1', 'alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphap', 'betap', 'alphaq', 'betaq', ] # Compute Calcium reversal potential for Cai = 5 nM self.VCa = nernst(Z_Ca, self.CCa_in0, self.CCa_out, self.T) # mV # Compute deff for that reversal potential iCaT = self.iCaT( self.pinf(self.Vm0), self.qinf(self.Vm0), self.Vm0) # mA/m2 iCaL = self.iCaL( self.cinf(self.Vm0), self.d1inf(self.Vm0), self.d2inf(self.CCa_in0), self.Vm0) # mA/m2 self.deff = -(iCaT + iCaL) / (Z_Ca * FARADAY * self.KCa * self.CCa_in0) * 1e-6 # m # Compute conversion factor from electrical current (mA/m2) to Calcium concentration (M) self.i2CCa = 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 _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 derC_Ca(self, C_Ca, iCaT, iCaL): ''' Compute the evolution of the Calcium concentration in submembranal space. :param Vm: membrane potential (mV) :param C_Ca: Calcium concentration in submembranal space (M) :param iCaT: inward, low-threshold Calcium current (mA/m2) :param iCaL: inward, high-threshold Calcium current (mA/m2) :return: derivative of Calcium concentration in submembranal space w.r.t. time (M/s) ''' return - self.i2CCa * (iCaT + iCaL) - C_Ca * self.KCa def iNa(self, m, h, Vm): ''' Compute the inward Sodium current per unit area. :param m: open-probability of m-gate :param h: inactivation-probability of h-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GNaMax * m**3 * h * (Vm - self.VNa) def iKd(self, n, Vm): ''' Compute the outward delayed-rectifier Potassium current per unit area. :param n: open-probability of n-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GKMax * n**4 * (Vm - self.VK) def iA(self, a, b, Vm): ''' Compute the outward A-type Potassium current per unit area. :param a: open-probability of a-gate :param b: open-probability of b-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GAMax * a**2 * b * (Vm - self.VK) def iCaT(self, p, q, Vm): ''' Compute the inward low-threshold Calcium current per unit area. :param p: open-probability of p-gate :param q: open-probability of q-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GTMax * p**2 * q * (Vm - self.VCa) def iCaL(self, c, d1, d2, Vm): ''' Compute the inward high-threshold Calcium current per unit area. :param c: open-probability of c-gate :param d1: open-probability of d1-gate :param d2: open-probability of d2-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GLMax * c**2 * d1 * d2 * (Vm - self.VCa) def iKCa(self, r, Vm): ''' Compute the outward, Calcium activated Potassium current per unit area. :param r: open-probability of r-gate :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GKCaMax * r**2 * (Vm - self.VK) def iLeak(self, Vm): ''' Compute the non-specific leakage current per unit area. :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GLeak * (Vm - self.VLeak) def iNet(self, Vm, states): ''' Compute net membrane current per unit area. ''' a, b, c, d1, d2, m, h, n, p, q, r, CCa_in = states # update VCa based on intracellular Calcium concentration self.VCa = nernst(Z_Ca, CCa_in, self.CCa_out, self.T) # mV return ( self.iNa(m, h, Vm) + self.iKd(n, Vm) + self.iA(a, b, Vm) + self.iCaT(p, q, Vm) + self.iCaL(c, d1, d2, Vm) + self.iKCa(r, Vm) + 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) d2eq = self.d2inf(self.CCa_in0) req = self.rinf(self.CCa_in0) return np.array([aeq, beq, ceq, d1eq, d2eq, meq, heq, neq, peq, qeq, req, self.CCa_in0]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' a, b, c, d1, d2, m, h, n, p, q, r, CCa_in = states dadt = self.derA(Vm, a) dbdt = self.derB(Vm, b) dcdt = self.derC(Vm, c) dd1dt = self.derD1(Vm, d1) dd2dt = self.derD2(CCa_in, d2) dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dpdt = self.derP(Vm, p) dqdt = self.derQ(Vm, q) drdt = self.derR(CCa_in, r) iCaT = self.iCaT(p, q, Vm) iCaL = self.iCaL(c, d1, d2, Vm) dCCaindt = self.derC_Ca(CCa_in, iCaT, iCaL) return [dadt, dbdt, dcdt, dd1dt, dd2dt, dmdt, dhdt, dndt, dpdt, dqdt, drdt, dCCaindt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants Ta = self.taua(Vm) alphaa_avg = np.mean(self.ainf(Vm) / Ta) betaa_avg = np.mean(1 / Ta) - alphaa_avg Tb = self.taub(Vm) alphab_avg = np.mean(self.binf(Vm) / Tb) betab_avg = np.mean(1 / Tb) - alphab_avg Tc = self.tauc(Vm) alphac_avg = np.mean(self.cinf(Vm) / Tc) betac_avg = np.mean(1 / Tc) - alphac_avg Td1 = self.taud1(Vm) alphad1_avg = np.mean(self.ainf(Vm) / Td1) betad1_avg = np.mean(1 / Td1) - alphad1_avg 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, CCa_in = states dadt = rates[0] * (1 - a) - rates[1] * a dbdt = rates[2] * (1 - b) - rates[3] * b dcdt = rates[4] * (1 - c) - rates[5] * c dd1dt = rates[6] * (1 - d1) - rates[7] * d1 dd2dt = self.derD2(CCa_in, 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(CCa_in, r) iCaT = self.iCaT(p, q, Vmeff) iCaL = self.iCaL(c, d1, d2, Vmeff) dCCaindt = self.derC_Ca(CCa_in, iCaT, iCaL) return [dadt, dbdt, dcdt, dd1dt, dd2dt, dmdt, dhdt, dndt, dpdt, dqdt, drdt, dCCaindt] diff --git a/PySONIC/plt/pltvars.py b/PySONIC/plt/pltvars.py index 08d0bb3..a41a7c1 100644 --- a/PySONIC/plt/pltvars.py +++ b/PySONIC/plt/pltvars.py @@ -1,803 +1,803 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-21 14:33:36 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-03-06 17:51:11 +# @Last Modified time: 2019-03-07 14:50:31 ''' Dictionary of plotting settings for output variables of the model. ''' pltvars = { 't_ms': { 'desc': 'time', 'label': 'time', 'unit': 'ms', 'factor': 1e3, 'onset': 1e-3 }, 't_us': { 'desc': 'time', 'label': 'time', 'unit': 'us', 'factor': 1e6, 'onset': 1e-6 }, 'Z': { 'desc': 'leaflets deflection', 'label': 'Z', 'unit': 'nm', 'factor': 1e9, 'min': -1.0, 'max': 10.0 }, 'ng': { 'desc': 'gas content', 'label': 'n_g', 'unit': '10^{-22}\ mol', 'factor': 1e22, 'min': 1.0, 'max': 15.0 # 'unit': 'ymol', # 'factor': 1e24, # 'min': 100.0, # 'max': 1500.0 }, 'Pac': { 'desc': 'acoustic pressure', 'label': 'P_{AC}', 'unit': 'kPa', 'factor': 1e-3, 'alias': 'bls.Pacoustic(t, meta["Adrive"] * states, Fdrive)' }, 'Pmavg': { 'desc': 'average intermolecular pressure', 'label': 'P_M', 'unit': 'kPa', 'factor': 1e-3, 'alias': 'bls.PMavgpred(df["Z"].values)' }, 'Telastic': { 'desc': 'leaflet elastic tension', 'label': 'T_E', 'unit': 'mN/m', 'factor': 1e3, 'alias': 'bls.TEleaflet(df["Z"].values)' }, 'Qm': { 'desc': 'charge density', # 'label': 'Q_m', 'label': 'charge', 'unit': 'nC/cm^2', 'factor': 1e5, 'min': -100, 'max': 50 }, 'Cm': { 'desc': 'membrane capacitance', 'label': 'C_m', 'unit': 'uF/cm^2', 'factor': 1e2, 'min': 0.0, 'max': 1.5, 'alias': 'bls.v_Capct(df["Z"].values)' }, 'Vm': { 'desc': 'membrane potential', 'label': 'V_m', 'unit': 'mV', 'factor': 1, }, 'Veff': { 'key': 'V', 'desc': 'effective membrane potential', 'label': 'V_{m, eff}', 'unit': 'mV', 'factor': 1e0 }, 'm': { 'desc': 'iNa activation gate opening', 'label': 'm-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'h': { 'desc': 'iNa inactivation gate opening', 'label': 'h-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'm2h': { 'desc': 'iNa relative conductance', 'label': 'm^2h', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': 'df["m"].values**2 * df["h"].values' }, 'm3h': { 'desc': 'iNa relative conductance', 'label': 'm^3h', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': 'df["m"].values**3 * df["h"].values' }, 'm4h': { 'desc': 'iNa relative conductance', 'label': 'm^4h', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': 'df["m"].values**4 * df["h"].values' }, 'n': { 'desc': 'iK activation gate opening', 'label': 'n-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'n4': { 'desc': 'iK relative conductance', 'label': 'n^4', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': 'df["n"].values**4' }, 'p': { 'desc': 'iM activation gate opening', 'label': 'p-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 's': { 'desc': 'iCa activation gates opening', 'label': 's-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'u': { 'desc': 'iCa inactivation gates opening', 'label': 'u-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 's2u': { 'desc': 'iT relative conductance', 'label': 's^2u', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': 'df["s"].values**2 * df["u"].values' }, 'p': { 'desc': 'iT activation gates opening', 'label': 'p-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'q': { 'desc': 'iT inactivation gates opening', 'label': 'q-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'p2q': { 'desc': 'iT relative conductance', 'label': 'p^2q', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': 'df["p"].values**2 * df["q"].values' }, 'r': { 'desc': 'iCaK Ca2+-gated activation gates opening', 'label': 'r-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'r2': { 'desc': 'iCaK relative conductance', 'label': 'r^2', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': 'df["r"].values**2' }, # 'q': { # 'desc': 'iCaL activation gates opening', # 'label': 'q-gate', # 'unit': None, # 'factor': 1, # 'min': -0.1, # 'max': 1.1 # }, # 'r': { # 'desc': 'iCaL inactivation gates opening', # 'label': 'r-gate', # 'unit': None, # 'factor': 1, # 'min': -0.1, # 'max': 1.1 # }, # 'q2r': { # 'desc': 'iCaL relative conductance', # 'label': 'q^2r', # 'unit': None, # 'factor': 1, # 'min': -0.1, # 'max': 1.1, # 'alias': 'df["q"].values**2 * df["r"].values' # }, # 'c': { # 'desc': 'iKCa activation gates opening', # 'label': 'c-gate', # 'unit': None, # 'factor': 1, # 'min': -0.1, # 'max': 1.1 # }, 'a': { 'desc': 'iA activation gates opening', 'label': 'a-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'b': { 'desc': 'iA inactivation gates opening', 'label': 'b-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'a2b': { 'desc': 'iA relative conductance', 'label': 'ab', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': 'df["a"].values**2 * df["b"].values' }, 'c': { 'desc': 'iL activation gates opening', 'label': 'c-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'd1': { 'desc': 'iL inactivation gates opening', 'label': 'd1-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'd2': { 'desc': 'iL Ca2+-gated inactivation gates opening', 'label': 'd2-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'c2d1d2': { 'desc': 'iL relative conductance', 'label': 'c^2d_1d_2', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': 'df["c"].values**2 * df["d1"].values * df["d2"].values' }, 'O': { 'desc': 'iH activation gate opening', 'label': 'O', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'OL': { 'desc': 'iH activation gate locked-opening', 'label': 'O_L', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, 'alias': '1 - df["O"].values - df["C"].values' }, 'O + 2OL': { 'desc': 'iH activation gate relative conductance', 'label': 'O\ +\ 2O_L', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 2.1, 'alias': 'df["O"].values + 2 * (1 - df["O"].values - df["C"].values)' }, 'P0': { 'desc': 'iH regulating factor activation', 'label': 'P_0', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'C_Ca': { 'desc': 'sumbmembrane Ca2+ concentration', 'label': '[Ca^{2+}]_i', 'unit': 'uM', 'factor': 1e6 # 'min': 0, # 'max': 150.0 }, 'C_Na': { 'desc': 'sumbmembrane Na+ concentration', 'label': '[Na^+]_i', 'unit': 'uM', 'factor': 1e6 # 'min': 0, # 'max': 150.0 }, 'C_Na_arb': { 'key': 'C_Na', 'desc': 'submembrane Na+ concentration', 'label': '[Na^+]', 'unit': 'arb.', 'factor': 1 }, 'C_Na_arb_activation': { 'key': 'A_Na', 'desc': 'Na+ dependent PumpNa current activation', 'label': 'A_{Na^+}', 'unit': 'arb', 'factor': 1 }, 'C_Ca_arb': { 'key': 'C_Ca', 'desc': 'submembrane Ca2+ concentration', 'label': '[Ca^{2+}]', 'unit': 'arb.', 'factor': 1 }, 'C_Ca_arb_activation': { 'key': 'A_Ca', 'desc': 'Ca2+ dependent Potassium current activation', 'label': 'A_{Ca^{2+}}', 'unit': 'arb', 'factor': 1 }, 'VLeak': { 'constant': 'neuron.VLeak', 'desc': 'non-specific leakage current resting potential', 'label': 'V_{leak}', 'unit': 'mV', 'factor': 1e0 }, 'iLeak': { 'desc': 'leakage current', 'label': 'I_{Leak}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iLeak(df["Vm"].values)' }, 'iNa': { 'desc': 'Sodium current', 'label': 'I_{Na}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iNa(df["m"].values, df["h"].values, df["Vm"].values)' }, 'iNa2': { 'desc': 'Sodium current', 'label': 'I_{Na}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iNa(df["m"].values, df["h"].values, df["Vm"].values, df["C_Na"].values)' }, 'iKd': { 'desc': 'delayed-rectifier Potassium current', 'label': 'I_K', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iKd(df["n"].values, df["Vm"].values)' }, 'iM': { 'desc': 'slow non-inactivating Potassium current', 'label': 'I_M', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iM(df["p"].values, df["Vm"].values)' }, 'iP': { 'desc': 'non-specific delayed current', 'label': 'I_P', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iP(df["p"].values, df["Vm"].values)' }, 'iA': { 'desc': 'transient Potassium current', 'label': 'I_A', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iA(df["a"].values, df["b"].values, df["Vm"].values)' }, 'iCaT': { 'desc': 'low-threshold (T-type) Calcium current', 'label': 'I_{CaT}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iCaT(df["s"].values, df["u"].values, df["Vm"].values)' }, 'iCaT2': { 'desc': 'low-threshold (T-type) Calcium current', 'label': 'I_{CaT}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iCaT(df["p"].values, df["q"].values, df["Vm"].values)' }, 'iCaL': { 'desc': 'high-threshold (L-type) Calcium current', 'label': 'I_{CaL}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iCaL(df["q"].values, df["r"].values, df["Vm"].values)' }, 'iCaL2': { 'desc': 'high-threshold (L-type) Calcium current', 'label': 'I_{CaL}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iCaL(df["c"].values, df["d1"].values, df["d2"].values, df["Vm"].values)' }, 'iCaK': { 'desc': 'Calcium-activated Potassium current', 'label': 'I_{CaK}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iCaK(df["r"].values, df["Vm"].values)' }, 'iCa': { 'desc': 'leech Calcium current', 'label': 'I_{Ca}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iCa(df["s"].values, df["Vm"].values)' }, 'iCa2': { 'desc': 'leech Calcium current', 'label': 'I_{Ca}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iCa(df["s"].values, df["Vm"].values, df["C_Ca"].values)' }, 'iH': { 'desc': 'hyperpolarization-activated cationic current', 'label': 'I_h', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iH(df["O"].values, df["C"].values, df["Vm"].values)' }, 'iKLeak': { 'desc': 'leakage Potassium current', 'label': 'I_{KLeak}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iKLeak(df["Vm"].values)' }, 'iKCa': { 'desc': 'Calcium-activated Potassium current', 'label': 'I_{KCa}', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.iKCa(df["A_Ca"].values, df["Vm"].values)' + 'alias': 'neuron.iKCa(df["r"].values, df["Vm"].values)' }, 'iKCa2': { 'desc': 'Calcium-activated Potassium current', 'label': 'I_{KCa}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iKCa(df["c"].values, df["Vm"].values)' }, 'iPumpNa': { 'desc': 'Outward current mimicking the activity of the NaK-ATPase pump', 'label': 'I_{PumpNa}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iPumpNa(df["A_Na"].values, df["Vm"].values)' }, 'iPumpNa2': { 'desc': 'Outward current mimicking the activity of the NaK-ATPase pump', 'label': 'I_{PumpNa}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iPumpNa(df["C_Na"].values)' }, 'iPumpCa2': { 'desc': 'Outward current describing the removal of Ca2+ from the intracellular space', 'label': 'I_{PumpCa}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iPumpCa(df["C_Ca"].values)' }, 'iNet': { 'desc': 'net current', 'label': 'I_{net}', 'unit': 'A/m^2', 'factor': 1e-3, 'alias': 'neuron.iNet(df["Vm"].values, neuron_states)' }, 'alphaa': { 'desc': 'iA a-gate activation rate', 'label': '\\alpha_a', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betaa': { 'desc': 'iA a-gate inactivation rate', 'label': '\\beta_a', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphab': { 'desc': 'iA b-gate activation rate', 'label': '\\alpha_b', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betab': { 'desc': 'iA b-gate inactivation rate', 'label': '\\beta_b', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphac': { 'desc': 'iL c-gate activation rate', 'label': '\\alpha_c', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betac': { 'desc': 'iL c-gate inactivation rate', 'label': '\\beta_c', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphad1': { 'desc': 'iL d1-gate activation rate', 'label': '\\alpha_d1', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betad1': { 'desc': 'iL d1-gate inactivation rate', 'label': '\\beta_d1', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alpham': { 'desc': 'iNa m-gate activation rate', 'label': '\\alpha_m', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betam': { 'desc': 'iNa m-gate inactivation rate', 'label': '\\beta_m', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphah': { 'desc': 'iNa h-gate activation rate', 'label': '\\alpha_h', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betah': { 'desc': 'iNa h-gate inactivation rate', 'label': '\\beta_h', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphan': { 'desc': 'iK n-gate activation rate', 'label': '\\alpha_n', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betan': { 'desc': 'iK n-gate inactivation rate', 'label': '\\beta_n', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphap': { 'desc': 'iM p-gate activation rate', 'label': '\\alpha_p', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betap': { 'desc': 'iM p-gate inactivation rate', 'label': '\\beta_p', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphas': { 'desc': 'iT s-gate activation rate', 'label': '\\alpha_s', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betas': { 'desc': 'iT s-gate inactivation rate', 'label': '\\beta_s', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphau': { 'desc': 'iT u-gate activation rate', 'label': '\\alpha_u', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betau': { 'desc': 'iT u-gate inactivation rate', 'label': '\\beta_u', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphaq': { 'desc': 'iT q-gate activation rate', 'label': '\\alpha_q', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betaq': { 'desc': 'iT q-gate inactivation rate', 'label': '\\beta_q', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'alphao': { 'desc': 'iH channels activation rate (between closed and open forms)', 'label': '\\alpha_O', 'unit': 'ms^{-1}', 'factor': 1e-3 }, 'betao': { 'desc': 'iH channels inactivation rate (between closed and open forms)', 'label': '\\beta_O', 'unit': 'ms^{-1}', 'factor': 1e-3 } } diff --git a/paper figures/deprecated/figQSS.py b/paper figures/deprecated/figQSS.py index 34ba853..0b404f8 100644 --- a/paper figures/deprecated/figQSS.py +++ b/paper figures/deprecated/figQSS.py @@ -1,355 +1,355 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2018-09-28 16:13:34 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-03-06 17:37:28 +# @Last Modified time: 2019-03-07 11:57:55 ''' Subpanels of the QSS approximation figure. ''' import os import logging import numpy as np from scipy.interpolate import interp1d import matplotlib.pyplot as plt import matplotlib import matplotlib.cm as cm from argparse import ArgumentParser from PySONIC.core import NeuronalBilayerSonophore from PySONIC.utils import logger, getLookups2D, selectDirDialog from PySONIC.neurons import getNeuronsDict # Plot parameters matplotlib.rcParams['pdf.fonttype'] = 42 matplotlib.rcParams['ps.fonttype'] = 42 matplotlib.rcParams['font.family'] = 'arial' # Figure basename figbase = os.path.splitext(__file__)[0] def plotQuasiSteadySystem(neuron, a, Fdrive, PRF, DC, fs=8, markers=['-', '--', '.-'], title=None): neuron = getNeuronsDict()[neuron]() # Determine spiking threshold Vthr = neuron.VT # mV Qthr = neuron.Cm0 * Vthr * 1e-3 # C/m2 # Get lookups amps, Qref, lookups2D, _ = getLookups2D(neuron.name, a=a, Fdrive=Fdrive) amps *= 1e-3 lookups1D = {key: interp1d(Qref, y2D, axis=1)(Qthr) for key, y2D in lookups2D.items()} # Remove unnecessary items ot get ON rates and effective potential at threshold charge rates_on = lookups1D rates_on.pop('ng') Vm_on = rates_on.pop('V') Vm_off = Qthr / neuron.Cm0 * 1e3 # Compute neuron OFF rates at current charge value rates_off = neuron.getRates(Vm_off) # Compute pulse-average quasi-steady states qsstates_pulse = np.empty((len(neuron.states_names), amps.size)) for j, x in enumerate(neuron.states_names): # If channel state, compute pulse-average steady-state values if x in neuron.getGates(): x = x.lower() alpha_str, beta_str = ['{}{}'.format(s, x) for s in ['alpha', 'beta']] alphax_pulse = rates_on[alpha_str] * DC + rates_off[alpha_str] * (1 - DC) betax_pulse = rates_on[beta_str] * DC + rates_off[beta_str] * (1 - DC) qsstates_pulse[j, :] = alphax_pulse / (alphax_pulse + betax_pulse) # Otherwise assume the state has reached a steady-state value for Vthr else: qsstates_pulse[j, :] = np.ones(amps.size) * neuron.steadyStates(Vthr)[j] # Compute quasi-steady ON and OFF currents iLeak_on = neuron.iLeak(Vm_on) iLeak_off = np.ones(amps.size) * neuron.iLeak(Vm_off) m = qsstates_pulse[0, :] h = qsstates_pulse[1, :] iNa_on = neuron.iNa(m, h, Vm_on) iNa_off = neuron.iNa(m, h, Vm_off) n = qsstates_pulse[2, :] iKd_on = neuron.iKd(n, Vm_on) iKd_off = neuron.iKd(n, Vm_off) p = qsstates_pulse[3, :] iM_on = neuron.iM(p, Vm_on) iM_off = neuron.iM(p, Vm_off) if neuron.name == 'LTS': s = qsstates_pulse[4, :] u = qsstates_pulse[5, :] - iCa_on = neuron.iCa(s, u, Vm_on) - iCa_off = neuron.iCa(s, u, Vm_off) + iCaT_on = neuron.iCaT(s, u, Vm_on) + iCaT_off = neuron.iCaT(s, u, Vm_off) iNet_on = neuron.iNet(Vm_on, qsstates_pulse) iNet_off = neuron.iNet(Vm_off, qsstates_pulse) # Compute quasi-steady ON, OFF and net charge variations, and threshold amplitude dQ_on = -iNet_on * DC / PRF dQ_off = -iNet_off * (1 - DC) / PRF dQ_net = dQ_on + dQ_off Athr = np.interp(0, dQ_net, amps, left=0., right=np.nan) # Create figure fig, axes = plt.subplots(4, 1, figsize=(4, 6)) axes[-1].set_xlabel('Amplitude (kPa)', fontsize=fs) for ax in axes: for skey in ['top', 'right']: ax.spines[skey].set_visible(False) ax.set_xscale('log') ax.set_xlim(1e1, 1e2) ax.set_xticks([1e1, 1e2]) for item in ax.get_xticklabels() + ax.get_yticklabels(): item.set_fontsize(fs) for item in ax.get_xticklabels(minor=True): item.set_visible(False) figname = '{} neuron thr dynamics {:.1f}nC_cm2 {:.0f}% DC'.format( neuron.name, Qthr * 1e5, DC * 1e2) fig.suptitle(figname, fontsize=fs) # Subplot 1: Vmeff ax = axes[0] ax.set_ylabel('Effective potential (mV)', fontsize=fs) Vbounds = (-120, -40) ax.set_ylim(Vbounds) ax.set_yticks([Vbounds[0], neuron.Vm0, Vbounds[1]]) ax.set_yticklabels(['{:.0f}'.format(Vbounds[0]), '$V_{m0}$', '{:.0f}'.format(Vbounds[1])]) ax.plot(amps, Vm_on, color='C0', label='ON') ax.plot(amps, Vm_off * np.ones(amps.size), '--', color='C0', label='OFF') ax.axhline(neuron.Vm0, linewidth=0.5, color='k') # Subplot 2: quasi-steady states ax = axes[1] ax.set_ylabel('Quasi-steady states', fontsize=fs) ax.set_yticks([0, 0.5, 0.6]) ax.set_yticklabels(['0', '0.5', '1']) ax.set_ylim([-0.05, 0.65]) d = .01 f = 1.03 xcut = ax.get_xlim()[0] for ycut in [0.54, 0.56]: ax.plot([xcut / f, xcut * f], [ycut - d, ycut + d], color='k', clip_on=False) for label, qsstate in zip(neuron.states_names, qsstates_pulse): if label == 'h': qsstate -= 0.4 ax.plot(amps, qsstate, label=label) # Subplot 3: currents ax = axes[2] ax.set_ylabel('QS Currents (mA/m2)', fontsize=fs) Ibounds = (-10, 10) ax.set_ylim(Ibounds) ax.set_yticks([Ibounds[0], 0.0, Ibounds[1]]) ax.plot(amps, iLeak_on, color='C0', label='$I_{Leak}$') ax.plot(amps, iLeak_off, '--', color='C0') ax.plot(amps, iNa_on, '-', color='C1', label='$I_{Na}$') ax.plot(amps, iNa_off, '--', color='C1') ax.plot(amps, iKd_on, '-', color='C2', label='$I_{K}$') ax.plot(amps, iKd_off, '--', color='C2') ax.plot(amps, iM_on, '-', color='C3', label='$I_{M}$') ax.plot(amps, iM_off, '--', color='C3') if neuron.name == 'LTS': - ax.plot(amps, iCa_on, color='C5', label='$I_{Ca}$') - ax.plot(amps, iCa_off, '--', color='C5') + ax.plot(amps, iCaT_on, color='C5', label='$I_{CaT}$') + ax.plot(amps, iCaT_off, '--', color='C5') ax.plot(amps, iNet_on, '-', color='k', label='$I_{Net}$') ax.plot(amps, iNet_off, '--', color='k') # Subplot 4: charge variations and activation threshold ax = axes[3] ax.set_ylabel('$\\rm \Delta Q_{QS}\ (nC/cm^2)$', fontsize=fs) dQbounds = (-0.06, 0.1) ax.set_ylim(dQbounds) ax.set_yticks([dQbounds[0], 0.0, dQbounds[1]]) ax.plot(amps, dQ_on, color='C0', label='ON') ax.plot(amps, dQ_off, '--', color='C0', label='OFF') ax.plot(amps, dQ_net, '-.', color='C0', label='Net') ax.plot([Athr] * 2, [ax.get_ylim()[0], 0], linestyle='--', color='k') ax.plot([Athr], [0], 'o', c='k') ax.axhline(0, color='k', linewidth=0.5) fig.tight_layout() fig.subplots_adjust(right=0.8) for ax in axes: ax.legend(loc='center right', fontsize=fs, frameon=False, bbox_to_anchor=(1.3, 0.5)) if title is not None: fig.canvas.set_window_title(title) return fig def plotdQvsDC(neuron, a, Fdrive, PRF, DCs, fs=8, title=None): neuron = getNeuronsDict()[neuron]() # Determine spiking threshold Vthr = neuron.VT # mV Qthr = neuron.Cm0 * Vthr * 1e-3 # C/m2 # Get lookups amps, Qref, lookups2D, _ = getLookups2D(neuron.name, a=a, Fdrive=Fdrive) amps *= 1e-3 lookups1D = {key: interp1d(Qref, y2D, axis=1)(Qthr) for key, y2D in lookups2D.items()} # Remove unnecessary items ot get ON rates and effective potential at threshold charge rates_on = lookups1D rates_on.pop('ng') Vm_on = rates_on.pop('V') rates_off = neuron.getRates(Vthr) # For each duty cycle, compute net charge variation at Qthr along the amplitude range, # and identify rheobase amplitude Athr = np.empty_like(DCs) dQnet = np.empty((DCs.size, amps.size)) for i, DC in enumerate(DCs): # Compute pulse-average quasi-steady states qsstates_pulse = np.empty((len(neuron.states_names), amps.size)) for j, x in enumerate(neuron.states_names): # If channel state, compute pulse-average steady-state values if x in neuron.getGates(): x = x.lower() alpha_str, beta_str = ['{}{}'.format(s, x) for s in ['alpha', 'beta']] alphax_pulse = rates_on[alpha_str] * DC + rates_off[alpha_str] * (1 - DC) betax_pulse = rates_on[beta_str] * DC + rates_off[beta_str] * (1 - DC) qsstates_pulse[j, :] = alphax_pulse / (alphax_pulse + betax_pulse) # Otherwise assume the state has reached a steady-state value for Vthr else: qsstates_pulse[j, :] = np.ones(amps.size) * neuron.steadyStates(Vthr)[j] # Compute the pulse average net current along the amplitude space iNet_on = neuron.iNet(Vm_on, qsstates_pulse) iNet_off = neuron.iNet(Vthr, qsstates_pulse) iNet_avg = iNet_on * DC + iNet_off * (1 - DC) dQnet[i, :] = -iNet_avg / PRF # Compute threshold amplitude Athr[i] = np.interp(0, dQnet[i, :], amps, left=0., right=np.nan) # Create figure fig, ax = plt.subplots(figsize=(4, 2)) figname = '{} neuron thr vs DC'.format(neuron.name, Qthr * 1e5) fig.suptitle(figname, fontsize=fs) for key in ['top', 'right']: ax.spines[key].set_visible(False) ax.set_xscale('log') for item in ax.get_xticklabels() + ax.get_yticklabels(): item.set_fontsize(fs) for item in ax.get_xticklabels(minor=True): item.set_visible(False) ax.set_xlabel('Amplitude (kPa)', fontsize=fs) ax.set_ylabel('$\\rm \Delta Q_{QS}\ (nC/cm^2)$', fontsize=fs) ax.set_xlim(1e1, 1e2) ax.axhline(0., linewidth=0.5, color='k') ax.set_ylim(-0.06, 0.12) ax.set_yticks([-0.05, 0.0, 0.10]) ax.set_yticklabels(['-0.05', '0', '0.10']) norm = matplotlib.colors.LogNorm(DCs.min(), DCs.max()) sm = cm.ScalarMappable(norm=norm, cmap='viridis') sm._A = [] for i, DC in enumerate(DCs): ax.plot(amps, dQnet[i, :], c=sm.to_rgba(DC), label='{:.0f}% DC'.format(DC * 1e2)) ax.plot([Athr[i]] * 2, [ax.get_ylim()[0], 0], linestyle='--', c=sm.to_rgba(DC)) ax.plot([Athr[i]], [0], 'o', c=sm.to_rgba(DC)) fig.tight_layout() fig.subplots_adjust(right=0.8) ax.legend(loc='center right', fontsize=fs, frameon=False, bbox_to_anchor=(1.3, 0.5)) if title is not None: fig.canvas.set_window_title(title) return fig def plotRheobaseAmps(neurons, a, Fdrive, DCs_dense, DCs_sparse, fs=8, title=None): fig, ax = plt.subplots(figsize=(3, 3)) ax.set_title('Rheobase amplitudes', fontsize=fs) ax.set_xlabel('Duty cycle (%)', fontsize=fs) ax.set_ylabel('$\\rm A_T\ (kPa)$', fontsize=fs) for key in ['top', 'right']: ax.spines[key].set_visible(False) for item in ax.get_xticklabels() + ax.get_yticklabels(): item.set_fontsize(fs) ax.set_xticks([25, 50, 75, 100]) ax.set_yscale('log') ax.set_ylim([10, 600]) norm = matplotlib.colors.LogNorm(DCs_sparse.min(), DCs_sparse.max()) sm = cm.ScalarMappable(norm=norm, cmap='viridis') sm._A = [] for i, neuron in enumerate(neurons): neuron = getNeuronsDict()[neuron]() nbls = NeuronalBilayerSonophore(a, neuron) Athrs_dense = nbls.findRheobaseAmps(DCs_dense, Fdrive, neuron.VT)[0] * 1e-3 # kPa Athrs_sparse = nbls.findRheobaseAmps(DCs_sparse, Fdrive, neuron.VT)[0] * 1e-3 # kPa ax.plot(DCs_dense * 1e2, Athrs_dense, label='{} neuron'.format(neuron.name)) for DC, Athr in zip(DCs_sparse, Athrs_sparse): ax.plot(DC * 1e2, Athr, 'o', label='{:.0f}% DC'.format(DC * 1e2) if i == len(neurons) - 1 else None, c=sm.to_rgba(DC)) ax.legend(fontsize=fs, frameon=False) fig.tight_layout() if title is not None: fig.canvas.set_window_title(title) return fig def main(): ap = ArgumentParser() # Runtime options ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity') ap.add_argument('-o', '--outdir', type=str, help='Output directory') ap.add_argument('-f', '--figset', type=str, nargs='+', help='Figure set', default='all') ap.add_argument('-s', '--save', default=False, action='store_true', help='Save output figures as pdf') args = ap.parse_args() loglevel = logging.DEBUG if args.verbose is True else logging.INFO logger.setLevel(loglevel) figset = args.figset if figset == 'all': figset = ['a', 'b', 'c', 'e'] logger.info('Generating panels {} of {}'.format(figset, figbase)) # Parameters a = 32e-9 # m Fdrive = 500e3 # Hz PRF = 100.0 # Hz DC = 0.5 DCs_sparse = np.array([5, 15, 50, 75, 95]) / 1e2 DCs_dense = np.arange(1, 101) / 1e2 # Figures figs = [] if 'a' in figset: figs += [ plotQuasiSteadySystem('RS', a, Fdrive, PRF, DC, title=figbase + 'a RS'), plotQuasiSteadySystem('LTS', a, Fdrive, PRF, DC, title=figbase + 'a LTS') ] if 'b' in figset: figs += [ plotdQvsDC('RS', a, Fdrive, PRF, DCs_sparse, title=figbase + 'b RS'), plotdQvsDC('LTS', a, Fdrive, PRF, DCs_sparse, title=figbase + 'b LTS') ] if 'c' in figset: figs.append(plotRheobaseAmps(['RS', 'LTS'], a, Fdrive, DCs_dense, DCs_sparse, title=figbase + 'c')) if args.save: outdir = selectDirDialog() if args.outdir is None else args.outdir if outdir == '': logger.error('No input directory chosen') return for fig in figs: figname = '{}.pdf'.format(fig.canvas.get_window_title()) fig.savefig(os.path.join(outdir, figname), transparent=True) else: plt.show() if __name__ == '__main__': main() diff --git a/scripts/testQSS.py b/scripts/testQSS.py new file mode 100644 index 0000000..ae0c8ba --- /dev/null +++ b/scripts/testQSS.py @@ -0,0 +1,100 @@ +# -*- coding: utf-8 -*- +# @Author: Theo Lemaire +# @Date: 2018-09-28 16:13:34 +# @Last Modified by: Theo Lemaire +# @Last Modified time: 2019-03-07 14:47:20 + + +import numpy as np +from scipy.interpolate import interp1d +import matplotlib.pyplot as plt +import matplotlib + +from PySONIC.utils import getLookups2D +from PySONIC.neurons import getNeuronsDict + + +# Plot parameters +matplotlib.rcParams['pdf.fonttype'] = 42 +matplotlib.rcParams['ps.fonttype'] = 42 +matplotlib.rcParams['font.family'] = 'arial' + + +def plotQSS(neuron, a, Fdrive, Adrive, fs=12): + + # Get lookups for specific (a, f, A) combination + Aref, Qref, lookups2D, _ = getLookups2D(neuron.name, a=a, Fdrive=Fdrive) + lookups1D = {key: interp1d(Aref, y2D, axis=0)(Adrive) for key, y2D in lookups2D.items()} + + # Remove unnecessary items ot get ON rates and effective potential + rates = lookups1D + rates.pop('ng') + Vm = rates.pop('V') + + # Compute quasi-steady states for each charge value + qsstates = np.empty((len(neuron.states_names), Qref.size)) + for j, x in enumerate(neuron.states_names): + # If channel state, compute steady-state values from rate constants + if x in neuron.getGates(): + x = x.lower() + alpha_str, beta_str = ['{}{}'.format(s, x) for s in ['alpha', 'beta']] + alphax = rates[alpha_str] + betax = rates[beta_str] + qsstates[j, :] = alphax / (alphax + betax) + # Otherwise assume the state has reached a steady-state value at the specific charge value + else: + qsstates[j, :] = np.array([neuron.steadyStates(Q / neuron.Cm0 * 1e3)[j] for Q in Qref]) + + # Compute quasi-steady currents + iNet = neuron.iNet(Vm, qsstates) + + # Create figure + fig, axes = plt.subplots(3, 1, figsize=(6, 7)) + axes[-1].set_xlabel('Charge Density (nC/cm2)', fontsize=fs) + for ax in axes: + for skey in ['top', 'right']: + ax.spines[skey].set_visible(False) + for item in ax.get_xticklabels() + ax.get_yticklabels(): + item.set_fontsize(fs) + for item in ax.get_xticklabels(minor=True): + item.set_visible(False) + figname = '{} neuron QSS dynamics @ {:.2f}kPa'.format(neuron.name, Adrive * 1e-3) + fig.suptitle(figname, fontsize=fs) + + # Subplot 1: Vmeff + ax = axes[0] + ax.set_ylabel('Effective potential (mV)', fontsize=fs) + ax.plot(Qref * 1e5, Vm, color='C0') + ax.axhline(neuron.Vm0, linewidth=0.5, color='k') + + # Subplot 2: quasi-steady states + ax = axes[1] + ax.set_ylabel('Quasi-steady states', fontsize=fs) + ax.set_yticks([0, 0.5, 1]) + ax.set_ylim([-0.05, 1.05]) + for label, qsstate in zip(neuron.states_names, qsstates): + ax.plot(Qref * 1e5, qsstate, label=label) + + # Subplot 3: currents + ax = axes[2] + ax.set_ylabel('QS Currents (mA/m2)', fontsize=fs) + ax.plot(Qref * 1e5, iNet, '-', color='C0', label='$I_{Net}$') + ax.axhline(0, color='k', linewidth=0.5) + + fig.tight_layout() + fig.subplots_adjust(right=0.8) + for ax in axes: + ax.legend(loc='center right', fontsize=fs, frameon=False, bbox_to_anchor=(1.3, 0.5)) + + return fig + + +neuron = getNeuronsDict()['STN']() +a = 32e-9 # m +Fdrive = 500e3 # Hz +amps = np.array([10, 20, 60]) * 1e3 # Pa + +for Adrive in amps: + plotQSS(neuron, a, Fdrive, Adrive) + +plt.show()