diff --git a/PySONIC/neurons/stn.py b/PySONIC/neurons/stn.py index 0c85b07..2268399 100644 --- a/PySONIC/neurons/stn.py +++ b/PySONIC/neurons/stn.py @@ -1,612 +1,613 @@ 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) 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 # 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 get_dCCa(self, C_Ca, Vm): ''' Return the time derivative of intracellular Calcium concentration given its current value at a specific membrane potential. :param C_Ca: Calcium concentration in submembranal space (M) :param Vm: membrane potential (mV) :return: time derivative of Calcium concentration in submembranal space (M/s) ''' self.VCa = nernst(Z_Ca, C_Ca, self.CCa_out, self.T) # mV iCaT = self.iCaT(self.pinf(Vm), self.qinf(Vm), Vm) iCaL = self.iCaL(self.cinf(Vm), self.d1inf(Vm), self.d2inf(C_Ca), Vm) return self.derC_Ca(C_Ca, iCaT, iCaL) def findCaeq(self, Vm): ''' Find the equilibrium intracellular Calcium concentration for a specific membrane potential. :param Vm: membrane potential (mV) :return: equilibrium Calcium concentration in submembranal space (M) ''' Ca_eq = brentq( lambda x: self.get_dCCa(x, Vm), self.CCa_in0 * 1e-4, self.CCa_in0 * 1e3, xtol=1e-16 ) return Ca_eq 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 currents(self, Vm, states): ''' Compute all membrane currents 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 { 'iNa': self.iNa(m, h, Vm), 'iKd': self.iKd(n, Vm), 'iA': self.iA(a, b, Vm), 'iCaT': self.iCaT(p, q, Vm), 'iCaL': self.iCaL(c, d1, d2, Vm), '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) + # C_Ca_eq = self.CCa_in0 C_Ca_eq = self.findCaeq(Vm) d2eq = self.d2inf(C_Ca_eq) req = self.rinf(C_Ca_eq) return np.array([aeq, beq, ceq, d1eq, d2eq, meq, heq, neq, peq, qeq, req, C_Ca_eq]) 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) + 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, 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/scripts/STN_regime_transition.py b/scripts/STN_regime_transition.py index 36d56b2..1547b87 100644 --- a/scripts/STN_regime_transition.py +++ b/scripts/STN_regime_transition.py @@ -1,295 +1,304 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2018-09-28 16:13:34 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-03-12 01:24:26 +# @Last Modified time: 2019-03-12 09:33:08 ''' Script to study STN transitions between different behavioral regimesl. ''' import os import pickle import numpy as np import matplotlib.pyplot as plt import matplotlib from argparse import ArgumentParser import logging from PySONIC.core import NeuronalBilayerSonophore from PySONIC.utils import * from PySONIC.postpro import getStableFixedPoints from PySONIC.neurons import getNeuronsDict # Plot parameters matplotlib.rcParams['pdf.fonttype'] = 42 matplotlib.rcParams['ps.fonttype'] = 42 matplotlib.rcParams['font.family'] = 'arial' # Set logging level logger.setLevel(logging.INFO) def getStableQmQSS(neuron, a, Fdrive, amps): # Compute net current profile for each amplitude, from QSS states and Vmeff profiles nbls = NeuronalBilayerSonophore(a, neuron, Fdrive) _, Qref, Vmeff, QS_states = nbls.getQSSvars(Fdrive, amps=amps) iNet = neuron.iNet(Vmeff, QS_states) # Find stable fixed points in iNet(Qref) profile Qstab = [] for i, Adrive in enumerate(amps): Qstab.append(getStableFixedPoints(Qref, -iNet[i, :])) return Qstab def getChargeStabilizationFromSims(inputdir, neuron, a, Fdrive, amps, tstim, PRF=100, DC=1.0): # Get filenames fnames = ['{}.pkl'.format(ASTIM_filecode(neuron.name, a, Fdrive, A, tstim, PRF, DC, 'sonic')) for A in amps] # Initialize output arrays - tstab = np.empty(amps.size) - Qstab = np.empty(amps.size) + t_stab = np.empty(amps.size) + Q_stab = np.empty(amps.size) + Ca_stab = np.empty(amps.size) # For each file for i, fn in enumerate(fnames): # Extract charge temporal profile during stimulus - fp = os.path.join(inputdir, 'STN', fn) + fp = os.path.join(inputdir, fn) # logger.info('loading data from file "{}"'.format(fn)) with open(fp, 'rb') as fh: frame = pickle.load(fh) df = frame['data'] t = df['t'].values Qm = df['Qm'].values + Ca = df['C_Ca'].values Qm = Qm[t < tstim] + Ca = Ca[t < tstim] t = t[t < tstim] dt = np.diff(t) # If charge signal is stable during last 100 ms of stimulus if np.ptp(Qm[-int(100e-3 // dt[0]):]) < 5e-5: # Compute instant of stabilization by iNet thresholding iNet_abs = np.abs(np.diff(Qm)) / dt - Qstab[i] = Qm[-1] - tstab[i] = t[np.where(iNet_abs > 1e-3)[0][-1] + 2] - logger.info('A = %.2f kPa: Qm stabilization around %.2f nC/cm2 from t = %.0f ms onward', - amps[i] * 1e-3, Qstab[i] * 1e5, tstab[i] * 1e3) + t_stab[i] = t[np.where(iNet_abs > 1e-3)[0][-1] + 2] + + # Get steady-state charge and Calcium concentration values + Q_stab[i] = Qm[-1] + Ca_stab[i] = Ca[-1] + + logger.debug('A = %.2f kPa: Qm stabilization around %.2f nC/cm2 from t = %.0f ms onward', + amps[i] * 1e-3, Q_stab[i] * 1e5, t_stab[i] * 1e3) # Otherwise, populate arrays with NaN else: - Qstab[i] = np.nan - tstab[i] = np.nan - logger.info('A = %.2f kPa: no Qm stabilization', amps[i] * 1e-3) + t_stab[i] = np.nan + Q_stab[i] = np.nan + Ca_stab[i] = np.nan + logger.debug('A = %.2f kPa: no Qm stabilization', amps[i] * 1e-3) - return Qstab, tstab + return t_stab, Q_stab, Ca_stab def plotQSSvars_vs_Qm(neuron, a, Fdrive, Adrive, fs=12): # Get quasi-steady states and effective membrane potential profiles at this amplitude nbls = NeuronalBilayerSonophore(a, neuron, Fdrive) _, Qref, Vmeff, QS_states = nbls.getQSSvars(Fdrive, amps=Adrive) # Compute QSS currents currents = neuron.currents(Vmeff, QS_states) iNet = sum(currents.values()) # Create figure fig, axes = plt.subplots(3, 1, figsize=(7, 9)) 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('$V_m^*$ (mV)', fontsize=fs) ax.plot(Qref * 1e5, Vmeff, color='C0') ax.axhline(neuron.Vm0, linewidth=0.5, color='k') # Subplot 2: quasi-steady states + colors = plt.get_cmap('tab10').colors + plt.get_cmap('Dark2').colors ax = axes[1] ax.set_ylabel('$X_\infty$', fontsize=fs) ax.set_yticks([0, 0.5, 1]) ax.set_ylim([-0.05, 1.05]) - for label, qsstate in zip(neuron.states_names, QS_states): - ax.plot(Qref * 1e5, qsstate, label=label) + for i, label in enumerate(neuron.states_names): + ax.plot(Qref * 1e5, QS_states[i], label=label, c=colors[i]) # Subplot 3: currents ax = axes[2] ax.set_ylabel('QSS currents (A/m2)', fontsize=fs) for k, I in currents.items(): ax.plot(Qref * 1e5, I * 1e-3, label=k) ax.plot(Qref * 1e5, iNet * 1e-3, color='k', label='iNet') ax.axhline(0, color='k', linewidth=0.5) fig.tight_layout() fig.subplots_adjust(right=0.8) for ax in axes[1:]: ax.legend(loc='center right', fontsize=fs, frameon=False, bbox_to_anchor=(1.3, 0.5)) fig.canvas.set_window_title( '{}_QSS_states_vs_Qm_{:.2f}kPa'.format(neuron.name, Adrive * 1e-3)) return fig def plotInetQSS_vs_Qm(neuron, a, Fdrive, amps, fs=12, cmap='viridis', zscale='lin'): # Compute net current profile for each amplitude, from QSS states and Vmeff profiles nbls = NeuronalBilayerSonophore(a, neuron, Fdrive) _, Qref, Vmeff, QS_states = nbls.getQSSvars(Fdrive, amps=amps) iNet = neuron.iNet(Vmeff, QS_states) # Define color code mymap = plt.get_cmap(cmap) zref = amps * 1e-3 if zscale == 'lin': norm = matplotlib.colors.Normalize(zref.min(), zref.max()) elif zscale == 'log': norm = matplotlib.colors.LogNorm(zref.min(), zref.max()) sm = matplotlib.cm.ScalarMappable(norm=norm, cmap=mymap) sm._A = [] # Create figure fig, ax = plt.subplots(figsize=(6, 4)) ax.set_xlabel('$\\rm Q_m\ (nC/cm^2)$', fontsize=fs) ax.set_ylabel('$\\rm I_{net, QSS}\ (A/m^2)$', fontsize=fs) for skey in ['top', 'right']: ax.spines[skey].set_visible(False) for item in ax.get_xticklabels() + ax.get_yticklabels(): item.set_fontsize(fs) figname = '{} neuron - QSS current imbalance vs. amplitude'.format(neuron.name) ax.set_title(figname, fontsize=fs) ax.axhline(0, color='k', linewidth=0.5) # Plot iNet profiles for each US amplitude (with specific color code) for i, Adrive in enumerate(amps): lbl = '{:.2f} kPa'.format(Adrive * 1e-3) c = sm.to_rgba(Adrive * 1e-3) ax.plot(Qref * 1e5, iNet[i] * 1e-3, label=lbl, c=c) for i, Adrive in enumerate(amps): Qstab = getStableFixedPoints(Qref, -iNet[i, :]) if Qstab is not None: ax.plot(Qstab * 1e5, np.zeros(Qstab.size), '.', c='k') fig.tight_layout() # Plot US amplitude colorbar fig.subplots_adjust(bottom=0.15, top=0.9, right=0.80, hspace=0.5) cbarax = fig.add_axes([0.85, 0.15, 0.03, 0.75]) fig.colorbar(sm, cax=cbarax) cbarax.set_ylabel('Amplitude (kPa)', fontsize=fs) for item in cbarax.get_yticklabels(): item.set_fontsize(fs) fig.canvas.set_window_title( '{}_iNet_QSS_vs_amp'.format(neuron.name)) return fig def compareEqChargesQSSvsSim(inputdir, neuron, a, Fdrive, amps, tstim, fs=12): # Get charge value that cancels out net current in QSS approx. and sim Qeq_QSS = getStableQmQSS(neuron, a, Fdrive, amps) - Qeq_sim, _ = getChargeStabilizationFromSims(inputdir, neuron, a, Fdrive, amps, tstim) + _, Qeq_sim, _ = getChargeStabilizationFromSims(inputdir, neuron, a, Fdrive, amps, tstim) # Plot Qm balancing net current as function of amplitude fig, ax = plt.subplots(figsize=(6, 4)) figname = '{} neuron - equilibrium charge vs. amplitude'.format(neuron.name) ax.set_title(figname) ax.set_xlabel('Amplitude (kPa)', fontsize=fs) ax.set_ylabel('$\\rm Q_{eq}\ (nC/cm^2)$', fontsize=fs) for skey in ['top', 'right']: ax.spines[skey].set_visible(False) for item in ax.get_xticklabels() + ax.get_yticklabels(): item.set_fontsize(fs) lgd = True for Adrive, Qstab in zip(amps, Qeq_QSS): if Qstab is not None: if lgd: lbl = 'QSS approximation' lgd = False else: lbl = None ax.plot(np.ones(Qstab.size) * Adrive * 1e-3, Qstab * 1e5, '.', c='C0', label=lbl) - # ax.plot(amps * 1e-3, Qeq_QSS * 1e5, label='QSS approximation') - ax.plot(amps * 1e-3, Qeq_sim * 1e5, c='C1', + ax.plot(amps * 1e-3, Qeq_sim * 1e5, '.', c='C1', label='end of {:.2f} s stimulus (simulation)'.format(tstim)) ax.legend(frameon=False, fontsize=fs) fig.tight_layout() fig.canvas.set_window_title( '{}_Qeq_vs_amp'.format(neuron.name)) return fig def main(): ap = ArgumentParser() # Stimulation parameters + ap.add_argument('-n', '--neuron', type=str, default='STN', help='Neuron type') ap.add_argument('-i', '--inputdir', type=str, default=None, help='Input directory') ap.add_argument('-o', '--outputdir', type=str, default=None, help='Output directory') ap.add_argument('-f', '--figset', type=str, nargs='+', help='Figure set', default='all') ap.add_argument('-c', '--cmap', type=str, default='viridis', help='Colormap name') ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity') ap.add_argument('-s', '--save', default=False, action='store_true', help='Save output figures as pdf') # Parse arguments args = ap.parse_args() logger.setLevel(logging.DEBUG if args.verbose else logging.INFO) figset = args.figset if figset == 'all': figset = ['a', 'b', 'c'] - neuron = getNeuronsDict()['STN']() + neuron = getNeuronsDict()[args.neuron]() a = 32e-9 # m Fdrive = 500e3 # Hz intensities = getLowIntensitiesSTN() # W/m2 amps = Intensity2Pressure(intensities) # Pa tstim = 1.0 # s figs = [] if 'a' in figset: - for Adrive in [amps[0], amps[-1]]: + for Adrive in [amps[0], amps[amps.size // 2], amps[-1]]: figs.append(plotQSSvars_vs_Qm(neuron, a, Fdrive, Adrive)) if 'b' in figset: figs.append(plotInetQSS_vs_Qm(neuron, a, Fdrive, amps)) if 'c' in figset: inputdir = args.inputdir if args.inputdir is not None else selectDirDialog() if inputdir == '': logger.error('no input directory') else: figs.append(compareEqChargesQSSvsSim(inputdir, neuron, a, Fdrive, amps, tstim)) if args.save: outputdir = args.outputdir if args.outputdir is not None else selectDirDialog() if outputdir == '': logger.error('no output directory') else: for fig in figs: s = fig.canvas.get_window_title() s = s.replace('(', '- ').replace('/', '_').replace(')', '') figname = '{}.pdf'.format(s) fig.savefig(os.path.join(outputdir, figname), transparent=True) else: plt.show() if __name__ == '__main__': main()