diff --git a/PointNICE/channels/base.py b/PointNICE/channels/base.py index 1839889..0715e16 100644 --- a/PointNICE/channels/base.py +++ b/PointNICE/channels/base.py @@ -1,126 +1,126 @@ #!/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: 2017-08-03 14:41:42 +# @Last Modified time: 2017-08-24 17:16:19 ''' Module standard API for all neuron mechanisms. Each mechanism class can use different methods to define the membrane dynamics of a specific neuron type. However, they must contain some mandatory attributes and methods in order to be properly imported in other PointNICE modules and used in NICE simulations. ''' import abc class BaseMech(metaclass=abc.ABCMeta): ''' Abstract class defining the common API (i.e. mandatory attributes and methods) of all subclasses implementing the channels mechanisms of specific 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: - **currNet**: 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 2-dimensional linear interpolators of "effective" coefficients. This method must return a list of derivatives ordered identically as in the states0 attribute. ''' @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' - @property - @abc.abstractmethod - def states_names(self): - return 'Should never reach here' + # @property + # @abc.abstractmethod + # def states_names(self): + # return 'Should never reach here' - @property - @abc.abstractmethod - def states0(self): - return 'Should never reach here' + # @property + # @abc.abstractmethod + # def states0(self): + # return 'Should never reach here' - @property - @abc.abstractmethod - def coeff_names(self): - return 'Should never reach here' + # @property + # @abc.abstractmethod + # def coeff_names(self): + # return 'Should never reach here' @abc.abstractmethod def currNet(self, Vm, states): ''' Compute the net ionic current per unit area. :param Vm: membrane potential (mV) :states: state probabilities of the ion channels :return: current per unit area (mA/m2) ''' @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, Adrive, Qm, states, interpolators): ''' Compute the effective derivatives of channel states, based on 2-dimensional linear interpolation of "effective" coefficients that summarize the system's behaviour over an acoustic cycle. :param Vm_eff: effective membrane potential (mV) :states: state probabilities of the ion channels :param interpolators: dictionary of 2-dimensional linear interpolators of "effective" rates over the 2D amplitude x charge input domain. ''' diff --git a/PointNICE/channels/cortical.py b/PointNICE/channels/cortical.py index 3125f19..3d55a6d 100644 --- a/PointNICE/channels/cortical.py +++ b/PointNICE/channels/cortical.py @@ -1,572 +1,576 @@ #!/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: 2017-08-22 14:49:00 +# @Last Modified time: 2017-08-24 17:22:05 ''' Channels mechanisms for thalamic neurons. ''' +import logging import numpy as np from .base import BaseMech +# Get package logger +logger = logging.getLogger('PointNICE') + class Cortical(BaseMech): ''' 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) - # Names and initial states of the channels state probabilities - states_names = ['m', 'h', 'n', 'p'] - states0 = np.array([]) - - # Names of the different coefficients to be averaged in a lookup table. - coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphap', 'betap'] def __init__(self): ''' Constructor of the class ''' - pass + + # 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'] 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 * (Vdiff - 13) / (np.exp(- (Vdiff - 13) / 4) - 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) ''' Vdiff = Vm - self.VT beta = (0.28 * (Vdiff - 40) / (np.exp((Vdiff - 40) / 5) - 1)) # 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 * (Vdiff - 15) / (np.exp(-(Vdiff - 15) / 5) - 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) ''' 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 currNa(self, m, h, Vm): ''' Compute the inward Sodium current per unit area. :param m: open-probability of Sodium channels :param h: inactivation-probability of Sodium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GNa = self.GNaMax * m**3 * h return GNa * (Vm - self.VNa) def currK(self, n, Vm): ''' Compute the outward, delayed-rectifier Potassium current per unit area. :param n: open-probability of delayed-rectifier Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GKMax * n**4 return GK * (Vm - self.VK) def currM(self, p, Vm): ''' Compute the outward, slow non-inactivating Potassium current per unit area. :param p: open-probability of the slow non-inactivating Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GM = self.GMMax * p return GM * (Vm - self.VK) def currL(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.GL * (Vm - self.VL) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, p = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currM(p, Vm) + self.currL(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, Adrive, Qm, states, interpolators): ''' Concrete implementation of the abstract API method. ''' rates = np.array([interpolators[rn](Adrive, Qm) 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) GL = 0.205 # Conductance of non-specific leakage current (S/m^2) VL = -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) 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) GL = 0.38 # Conductance of non-specific leakage current (S/m^2) VL = -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) 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) GL = 0.19 # Conductance of non-specific leakage current (S/m^2) VCa = 120.0 # # Calcium Nernst potential (mV) VL = -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) - 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 currCa(self, s, u, Vm): ''' Compute the inward, low-threshold Calcium current per unit area. :param s: open-probability of the S-type activation gate of Calcium channels :param u: open-probability of the U-type inactivation gate of Calcium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GT = self.GTMax * s**2 * u return GT * (Vm - self.VCa) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, p, s, u = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currM(p, Vm) + self.currCa(s, u, Vm) + self.currL(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 Calcium channel gates steady-states seq = self.sinf(Vm) ueq = self.uinf(Vm) Ca_eqstates = np.array([seq, ueq]) # 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) sinf = self.sinf(Vm) as_avg = np.mean(sinf / Ts) bs_avg = np.mean(1 / Ts) - as_avg Tu = np.array([self.tauu(v) for v in Vm]) uinf = self.uinf(Vm) au_avg = np.mean(uinf / 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, Adrive, Qm, states, interpolators): ''' 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(Adrive, Qm, NaK_states, interpolators) # Compute Calcium channels states derivatives Ca_rates = np.array([interpolators[rn](Adrive, Qm) 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] diff --git a/PointNICE/channels/leech.py b/PointNICE/channels/leech.py index dcdf78e..9c5de48 100644 --- a/PointNICE/channels/leech.py +++ b/PointNICE/channels/leech.py @@ -1,377 +1,381 @@ #!/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: 2017-08-22 10:45:35 +# @Last Modified time: 2017-08-24 17:20:49 ''' Channels mechanisms for leech ganglion neurons. ''' +import logging from functools import partialmethod import numpy as np from .base import BaseMech +# Get package logger +logger = logging.getLogger('PointNICE') + class LeechTouch(BaseMech): ''' 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.0 # 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) VKCa = -62.0 # Calcium-dependent, Potassium current Nernst potential (mV) VL = -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) GL = 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) surface = 6434.0 # surface of cell assumed as a single soma (um2) curr_factor = 1e6 # 1/nA to 1/mA tau_Na_removal = 16.0 # Time constant for the removal of Sodium ions from the pool (s) tau_Ca_removal = 1.25 # Time constant for the removal of Calcium ions from the pool (s) tau_PumpNa_act = 0.1 # Time constant for the PumpNa current activation from Sodium ions (s) tau_KCa_act = 0.01 # Time constant for the KCa current activation from Calcium ions (s) - # Names of the channels state probabilities - states_names = ['m', 'h', 'n', 's', 'C_Na', 'A_Na', 'C_Ca', 'A_Ca'] - - # Names of the channels effective coefficients - coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas'] - - # initial channel probabilities - states0 = np.array([]) 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 = 0.016 * self.surface / self.curr_factor self.K_Ca = 0.1 * self.surface / self.curr_factor # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) # ----------------- 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 currNa(self, m, h, Vm): ''' Sodium inward current. ''' return self.GNaMax * m**3 * h * (Vm - self.VNa) def currK(self, n, Vm): ''' Potassium outward current. ''' return self.GKMax * n**2 * (Vm - self.VK) def currCa(self, s, Vm): ''' Calcium inward current. ''' return self.GCaMax * s * (Vm - self.VCa) def currKCa(self, A_Ca, Vm): ''' Calcium-activated Potassium outward current. ''' return self.GKCaMax * A_Ca * (Vm - self.VKCa) def currPumpNa(self, A_Na, Vm): ''' Outward current mimicking the activity of the NaK-ATPase pump. ''' return self.GPumpNa * A_Na * (Vm - self.VPumpNa) def currL(self, Vm): ''' Leakage current. ''' return self.GL * (Vm - self.VL) def currNet(self, Vm, states): ''' Compute the net ionic current per unit area. :param Vm: membrane potential (mV) :states: state probabilities of the ion channels :return: current per unit area (mA/m2) ''' m, h, n, s, _, A_Na, _, A_Ca = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currCa(s, Vm) + self.currL(Vm) + self.currPumpNa(A_Na, Vm) + self.currKCa(A_Ca, Vm)) # mA/m2 def steadyStates(self, Vm): ''' Compute the Sodium, Potassium and Calcium channel gates steady-states for a specific membrane potential value. :param Vm: membrane potential (mV) :return: array of steady-states ''' # 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.currNa(meq, heq, Vm) CNa_eq = self.K_Na * (-INa_eq) ANa_eq = CNa_eq # print('initial Na current: {:.2f} mA/m2'.format(INa_eq)) # print('initial Na concentration in pool: {:.2f} arb. unit'.format(CNa_eq)) # KCa current pool concentration and activation steady-state ICa_eq = self.currCa(seq, Vm) CCa_eq = self.K_Ca * (-ICa_eq) ACa_eq = CCa_eq # print('initial Ca current: {:.2f} mA/m2'.format(ICa_eq)) # print('initial Ca concentration in pool: {:.2f} arb. unit'.format(CCa_eq)) return np.array([meq, heq, neq, seq, CNa_eq, ANa_eq, CCa_eq, ACa_eq]) 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) ''' # 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.currNa(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.currCa(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): ''' 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) ''' # 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, Adrive, Qm, states, interpolators): ''' Compute the effective derivatives of channel states, based on 2-dimensional linear interpolation of "effective" coefficients that summarize the system's behaviour over an acoustic cycle. :param Vm_eff: effective membrane potential (mV) :states: state probabilities of the ion channels :param interpolators: dictionary of 2-dimensional linear interpolators of "effective" rates over the 2D amplitude x charge input domain. ''' rates = np.array([interpolators[rn](Adrive, Qm) for rn in self.coeff_names]) Vmeff = interpolators['V'](Adrive, Qm) # 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 - m) - rates[7] * s # PumpNa current pool concentration and activation state I_Na = self.currNa(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.currCa(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] diff --git a/PointNICE/channels/thalamic.py b/PointNICE/channels/thalamic.py index 14fd16e..22b4a89 100644 --- a/PointNICE/channels/thalamic.py +++ b/PointNICE/channels/thalamic.py @@ -1,771 +1,773 @@ #!/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: 2017-08-22 14:48:08 +# @Last Modified time: 2017-08-24 17:20:30 ''' Channels mechanisms for thalamic neurons. ''' +import logging import numpy as np from .base import BaseMech +# Get package logger +logger = logging.getLogger('PointNICE') + class Thalamic(BaseMech): ''' 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) - # Names and initial states of the channels state probabilities - states_names = ['m', 'h', 'n', 's', 'u'] - states0 = np.array([]) - - # Names of the different coefficients to be averaged in a lookup table. - coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas', - 'alphau', 'betau'] - - def __init__(self): ''' Constructor of the class ''' - pass + + # 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'] 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 * (Vdiff - 13) / (np.exp(- (Vdiff - 13) / 4) - 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) ''' Vdiff = Vm - self.VT beta = (0.28 * (Vdiff - 40) / (np.exp((Vdiff - 40) / 5) - 1)) # 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 * (Vdiff - 15) / (np.exp(-(Vdiff - 15) / 5) - 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) ''' 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 currNa(self, m, h, Vm): ''' Compute the inward Sodium current per unit area. :param m: open-probability of Sodium channels :param h: inactivation-probability of Sodium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GNa = self.GNaMax * m**3 * h return GNa * (Vm - self.VNa) def currK(self, n, Vm): ''' Compute the outward delayed-rectifier Potassium current per unit area. :param n: open-probability of delayed-rectifier Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GKMax * n**4 return GK * (Vm - self.VK) def currCa(self, s, u, Vm): ''' Compute the inward Calcium current per unit area. :param s: open-probability of the S-type activation gate of Calcium channels :param u: open-probability of the U-type inactivation gate of Calcium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GT = self.GTMax * s**2 * u return GT * (Vm - self.VCa) def currL(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.GL * (Vm - self.VL) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, u = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currCa(s, u, Vm) + self.currL(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) sinf = self.sinf(Vm) as_avg = np.mean(sinf / Ts) bs_avg = np.mean(1 / Ts) - as_avg Tu = np.array([self.tauu(v) for v in Vm]) uinf = self.uinf(Vm) au_avg = np.mean(uinf / 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, Adrive, Qm, states, interpolators): ''' Concrete implementation of the abstract API method. ''' rates = np.array([interpolators[rn](Adrive, Qm) 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) GL = 0.5 # Conductance of non-specific leakage current (S/m^2) VL = -90.0 # Non-specific leakage Nernst potential (mV) VT = -67.0 # Spike threshold adjustment parameter (mV) 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) GL = 0.1 # Conductance of non-specific leakage current (S/m^2) Vh = -40.0 # Mixed cationic current reversal potential (mV) VL = -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) deff = 100e-9 # effective depth beneath membrane for intracellular [Ca2+] calculation F_Ca = 1.92988e5 # Faraday constant for bivalent ion (Coulomb / mole) 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) 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 * self.F_Ca) # Define names of the channels state probabilities self.states_names += ['O', 'C', 'P0', 'C_Ca'] # 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): ''' 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) :return: derivative of proportion w.r.t. time (s-1) ''' return self.k2 * (1 - P0) - self.k1 * P0 * C_Ca**self.nCa def derC_Ca(self, C_Ca, ICa): ''' Compute the evolution of the Calcium concentration in submembranal space. Model of Ca2+ buffering and contribution from iCa 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 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 (self.CCa_min - C_Ca) / self.tau_Ca_removal - self.iT_2_CCa * ICa def currKL(self, Vm): ''' Compute the voltage-dependent leak Potassium current per unit area. :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GKL * (Vm - self.VK) def currH(self, O, C, Vm): ''' Compute the outward mixed cationic current per unit area. :param O: proportion of the channels in open form :param OL: proportion of the channels in locked-open form :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' OL = 1 - O - C return self.GhMax * (O + 2 * OL) * (Vm - self.Vh) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, u, O, C, _, _ = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currCa(s, u, Vm) + self.currKL(Vm) + self.currH(O, C, Vm) + self.currL(Vm)) # mA/m2 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.currCa(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 BA = self.betao(Vm) / self.alphao(Vm) P0_eq = self.k2 / (self.k2 + self.k1 * CCa_eq**self.nCa) 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]) # 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 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.currCa(s, u, Vm) dCCa_dt = self.derC_Ca(C_Ca, ICa) return NaKCa_derstates + [dO_dt, dC_dt, dP0_dt, dCCa_dt] 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, Adrive, Qm, states, interpolators): ''' Concrete implementation of the abstract API method. ''' rates = np.array([interpolators[rn](Adrive, Qm) for rn in self.coeff_names]) Vmeff = interpolators['V'](Adrive, Qm) # Unpack states m, h, n, s, u, O, C, P0, C_Ca = states # INa, IK, ICa 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.currCa(s, u, Vmeff) dCCa_dt = self.derC_Ca(C_Ca, ICa_eff) # Merge derivatives and return return [dmdt, dhdt, dndt, dsdt, dudt, dO_dt, dC_dt, dP0_dt, dCCa_dt] diff --git a/PointNICE/constants.py b/PointNICE/constants.py index 6e7b00f..b4fa456 100644 --- a/PointNICE/constants.py +++ b/PointNICE/constants.py @@ -1,45 +1,45 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-11-04 13:23:31 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-08-22 17:29:02 +# @Last Modified time: 2017-08-24 13:41:31 ''' Algorithmic constants used in the core modules. ''' # Fitting and pre-processing LJFIT_PM_MAX = 1e8 # intermolecular pressure at the deflection lower bound for LJ fitting (Pa) PNET_EQ_MAX = 1e-1 # error threshold for net pressure at computed equilibrium position (Pa) PMAVG_STD_ERR_MAX = 1500 # error threshold in nonlinear fit of molecular pressure (Pa) # Generic integration constants NPC_FULL = 1000 # nb of samples per acoustic period in full system SOLVER_NSTEPS = 1000 # maximum number of steps allowed during one call to the LSODA/DOP853 solvers CLASSIC_DS_FACTOR = 3 # time downsampling factor applied to output arrays of classic simulations # Effective integration DT_EFF = 5e-5 # time step for effective integration (s) # DT_EFF = 1e-6 # time step for effective integration (s) # Mechanical simulations Z_ERR_MAX = 1e-11 # periodic convergence threshold for deflection gas content (m) NG_ERR_MAX = 1e-24 # periodic convergence threshold for gas content (mol) # Hybrid integration NPC_HH = 40 # nb of samples per acoustic period in HH system DQ_UPDATE = 1e-5 # charge evolution threshold between two hybrid integrations (C/m2) DT_UPDATE = 5e-4 # time interval between two hybrid integrations (s) # Titrations TITRATION_AMAX = 2e5 # initial acoustic pressure upper bound for titration procedure (Pa) TITRATION_TMAX = 2e-1 # initial stimulus duration upper bound for titration procedure (Pa) TITRATION_DA_THR = 1e3 # acoustic pressure search range threshold for titration procedure (Pa) TITRATION_DT_THR = 1e-3 # stimulus duration search range threshold for titration procedure (s) # Spike detection SPIKE_MIN_QAMP = 10e-5 # threshold amplitude for spike detection on charge signal (C/m2) +SPIKE_MIN_VAMP = 10.0 # threshold amplitude for spike detection on potential signal (mV) SPIKE_MIN_DT = 1e-3 # minimal time interval for spike detection on charge signal (s) - diff --git a/PointNICE/plt/pltutils.py b/PointNICE/plt/pltutils.py index 3d40833..0e881fd 100644 --- a/PointNICE/plt/pltutils.py +++ b/PointNICE/plt/pltutils.py @@ -1,449 +1,471 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-23 14:55:37 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-08-23 15:41:49 +# @Last Modified time: 2017-08-24 17:45:07 ''' Plotting utilities ''' import pickle import ntpath import re import inspect import tkinter as tk from tkinter import filedialog import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Rectangle from .. import channels from ..bls import BilayerSonophore from .pltvars import pltvars +# Define global variables +neuron = None +bls = None + +# Regular expression for input files +rgxp = re.compile('([E,A])STIM_([A-Za-z]*)_(.*).pkl') + + class InteractiveLegend(object): """ Class defining an interactive matplotlib legend, where lines visibility can be toggled by simply clicking on the corresponding legend label. Other graphic objects can also be associated to the toggle of a specific line Adapted from: http://stackoverflow.com/questions/31410043/hiding-lines-after-showing-a-pyplot-figure """ def __init__(self, legend, aliases): self.legend = legend self.fig = legend.axes.figure self.lookup_artist, self.lookup_handle = self._build_lookups(legend) self._setup_connections() self.handles_aliases = aliases self.update() def _setup_connections(self): for artist in self.legend.texts + self.legend.legendHandles: artist.set_picker(10) # 10 points tolerance self.fig.canvas.mpl_connect('pick_event', self.on_pick) def _build_lookups(self, legend): ''' Method of the InteractiveLegend class building the legend lookups. ''' labels = [t.get_text() for t in legend.texts] handles = legend.legendHandles label2handle = dict(zip(labels, handles)) handle2text = dict(zip(handles, legend.texts)) lookup_artist = {} lookup_handle = {} for artist in legend.axes.get_children(): if artist.get_label() in labels: handle = label2handle[artist.get_label()] lookup_handle[artist] = handle lookup_artist[handle] = artist lookup_artist[handle2text[handle]] = artist lookup_handle.update(zip(handles, handles)) lookup_handle.update(zip(legend.texts, handles)) return lookup_artist, lookup_handle def on_pick(self, event): handle = event.artist if handle in self.lookup_artist: artist = self.lookup_artist[handle] artist.set_visible(not artist.get_visible()) self.update() def update(self): for artist in self.lookup_artist.values(): handle = self.lookup_handle[artist] if artist.get_visible(): handle.set_visible(True) if artist in self.handles_aliases: for al in self.handles_aliases[artist]: al.set_visible(True) else: handle.set_visible(False) if artist in self.handles_aliases: for al in self.handles_aliases[artist]: al.set_visible(False) self.fig.canvas.draw() def show(self): ''' showing the interactive legend ''' plt.show() def getPatchesLoc(t, states): ''' Determine the location of stimulus patches. :param t: simulation time vector (s). :param states: a vector of stimulation state (ON/OFF) at each instant in time. :return: 3-tuple with number of patches, timing of STIM-ON an STIM-OFF instants. ''' # Compute states derivatives and identify bounds indexes of pulses dstates = np.diff(states) ipatch_on = np.insert(np.where(dstates > 0.0)[0] + 1, 0, 0) ipatch_off = np.where(dstates < 0.0)[0] # Get time instants for pulses ON and OFF npatches = ipatch_on.size tpatch_on = t[ipatch_on] tpatch_off = t[ipatch_off] # return 3-tuple with #patches, pulse ON and pulse OFF instants return (npatches, tpatch_on, tpatch_off) def SaveFigDialog(dirname, filename): """ Open a FileSaveDialogBox to set the directory and name of the figure to be saved. The default directory and filename are given, and the default extension is ".pdf" :param dirname: default directory :param filename: default filename :return: full path to the chosen filename """ root = tk.Tk() root.withdraw() filename_out = filedialog.asksaveasfilename(defaultextension=".pdf", initialdir=dirname, initialfile=filename) return filename_out def plotComp(yvars, filepaths, fs=15, show_patches=True): ''' Compare profiles of several specific output variables of NICE simulations. :param yvars: list of variables names to extract and compare :param filepaths: list of full paths to output data files to be compared :param fs: labels fontsize :param show_patches: boolean indicating whether to indicate periods of stimulation with colored rectangular patches ''' nvars = len(yvars) # x and y variables plotting information t_plt = pltvars['t'] y_pltvars = [pltvars[key] for key in yvars] # Dictionary of neurons neurons = {} for _, obj in inspect.getmembers(channels): if inspect.isclass(obj) and isinstance(obj.name, str): neurons[obj.name] = obj - # Regular expression for input files - rgxp = re.compile('sim_([A-Za-z]*)_(.*).pkl') - - # Initialize figure and axes if nvars == 1: _, ax = plt.subplots(figsize=(11, 4)) axes = [ax] else: _, axes = plt.subplots(nvars, 1, figsize=(11, min(3 * nvars, 9))) labels = [ntpath.basename(fp)[4:-4].replace('_', ' ') for fp in filepaths] for i in range(nvars): ax = axes[i] pltvar = y_pltvars[i] if 'min' in pltvar and 'max' in pltvar: ax.set_ylim(pltvar['min'], pltvar['max']) if pltvar['unit']: ax.set_ylabel('${}\ ({})$'.format(pltvar['label'], pltvar['unit']), fontsize=fs) else: ax.set_ylabel('${}$'.format(pltvar['label']), fontsize=fs) if i < nvars - 1: ax.get_xaxis().set_ticklabels([]) else: ax.set_xlabel('${}\ ({})$'.format(t_plt['label'], t_plt['unit']), fontsize=fs) for tick in ax.xaxis.get_major_ticks(): tick.label.set_fontsize(fs) ax.locator_params(axis='y', nbins=2) for tick in ax.yaxis.get_major_ticks(): tick.label.set_fontsize(fs) # Loop through data files tstim_ref = 0.0 nstim = 0 j = 0 aliases = {} for filepath in filepaths: pkl_filename = ntpath.basename(filepath) # Retrieve neuron name mo = rgxp.fullmatch(pkl_filename) if not mo: print('Error: PKL file does not match regular expression pattern') quit() - neuron_name = mo.group(1) + sim_type = mo.group(1) + neuron_name = mo.group(2) # Load data print('Loading data from "' + pkl_filename + '"') with open(filepath, 'rb') as pkl_file: data = pickle.load(pkl_file) # Extract useful variables t = data['t'] states = data['states'] tstim = data['tstim'] - Fdrive = data['Fdrive'] - params = data['params'] - a = data['a'] - d = data['d'] - geom = {"a": a, "d": d} # Initialize BLS and channels mechanism - global neuron, bls + global neuron neuron = neurons[neuron_name]() - Qm0 = neuron.Cm0 * neuron.Vm0 * 1e-3 - bls = BilayerSonophore(geom, params, Fdrive, neuron.Cm0, Qm0) + neuron_states = [data[sn] for sn in neuron.states_names] + + # Initialize BLS + if sim_type == 'A': + global bls + params = data['params'] + Fdrive = data['Fdrive'] + a = data['a'] + d = data['d'] + geom = {"a": a, "d": d} + Qm0 = neuron.Cm0 * neuron.Vm0 * 1e-3 + bls = BilayerSonophore(geom, params, Fdrive, neuron.Cm0, Qm0) # Get data of variables to plot vrs = [] for i in range(nvars): pltvar = y_pltvars[i] if 'alias' in pltvar: var = eval(pltvar['alias']) elif 'key' in pltvar: var = data[pltvar['key']] else: var = data[yvars[i]] vrs.append(var) # Determine patches location npatches, tpatch_on, tpatch_off = getPatchesLoc(t, states) # Adding onset to all signals if t_plt['onset'] > 0.0: t = np.insert(t + t_plt['onset'], 0, 0.0) for i in range(nvars): vrs[i] = np.insert(vrs[i], 0, vrs[i][0]) tpatch_on += t_plt['onset'] tpatch_off += t_plt['onset'] # Plotting handles = [axes[i].plot(t * t_plt['factor'], vrs[i] * y_pltvars[i]['factor'], linewidth=2, label=labels[j]) for i in range(nvars)] plt.tight_layout() if show_patches: k = 0 # stimulation patches for ax in axes: handle = handles[k] (ybottom, ytop) = ax.get_ylim() la = [] for i in range(npatches): la.append(ax.add_patch(Rectangle((tpatch_on[i] * t_plt['factor'], ybottom), (tpatch_off[i] - tpatch_on[i]) * t_plt['factor'], ytop - ybottom, color=handle[0].get_color(), alpha=0.2))) aliases[handle[0]] = la k += 1 if tstim != tstim_ref: if nstim == 0: nstim += 1 tstim_ref = tstim else: print('Warning: comparing different stimulation durations') j += 1 iLegends = [] for k in range(nvars): axes[k].legend(loc='upper left', fontsize=fs) iLegends.append(InteractiveLegend(axes[k].legend_, aliases)) plt.show() -def plotBatch(yvars, positions, directory, filepaths, plt_show=True, plt_save=False, - ask_before_save=1, fig_ext='png', tag='fig', fs=15, show_patches=True): + +def plotBatch(vars_dict, directory, filepaths, plt_show=True, plt_save=False, + ask_before_save=1, fig_ext='png', tag='fig', fs=15, lw=4, show_patches=True): ''' Plot a figure with profiles of several specific NICE output variables, for several NICE simulations. - :param yvars: list of variables names to extract and compare + :param vars_dict: dict of lists of variables names to extract and plot together :param positions: subplot indexes of each variable :param filepaths: list of full paths to output data files to be compared :param fs: labels fontsize :param show_patches: boolean indicating whether to indicate periods of stimulation with colored rectangular patches ''' - nvars = len(yvars) - naxes = np.unique(positions).size + labels = list(vars_dict.keys()) + naxes = len(vars_dict) # x and y variables plotting information t_plt = pltvars['t'] - y_pltvars = [pltvars[key] for key in yvars] # Dictionary of neurons neurons = {} for _, obj in inspect.getmembers(channels): if inspect.isclass(obj) and isinstance(obj.name, str): neurons[obj.name] = obj - # Regular expression for input files - rgxp = re.compile('sim_([A-Za-z]*)_(.*).pkl') # Loop through data files for filepath in filepaths: # Get code from file name pkl_filename = ntpath.basename(filepath) filecode = pkl_filename[0:-4] # Retrieve neuron name mo = rgxp.fullmatch(pkl_filename) if not mo: print('Error: PKL file does not match regular expression pattern') quit() - neuron_name = mo.group(1) + sim_type = mo.group(1) + neuron_name = mo.group(2) # Load data print('Loading data from "' + pkl_filename + '"') with open(filepath, 'rb') as pkl_file: data = pickle.load(pkl_file) # Extract variables print('Extracting variables') t = data['t'] states = data['states'] - Fdrive = data['Fdrive'] - params = data['params'] - a = data['a'] - d = data['d'] - geom = {"a": a, "d": d} - # Initialize BLS and channels mechanism - global bls, neuron + # Initialize channel mechanism + global neuron neuron = neurons[neuron_name]() - Qm0 = neuron.Cm0 * neuron.Vm0 * 1e-3 - bls = BilayerSonophore(geom, params, Fdrive, neuron.Cm0, Qm0) - - # Get data of variables to plot - vrs = [] - for i in range(nvars): - pltvar = y_pltvars[i] - if 'alias' in pltvar: - var = eval(pltvar['alias']) - elif 'key' in pltvar: - var = data[pltvar['key']] - elif 'constant' in pltvar: - var = eval(pltvar['constant']) * np.ones(t.size) - else: - # var = data[varlist[i]] - var = data[yvars[i]] - vrs.append(var) + neuron_states = [data[sn] for sn in neuron.states_names] + + # Initialize BLS + if sim_type == 'A': + global bls + params = data['params'] + Fdrive = data['Fdrive'] + a = data['a'] + d = data['d'] + geom = {"a": a, "d": d} + Qm0 = neuron.Cm0 * neuron.Vm0 * 1e-3 + bls = BilayerSonophore(geom, params, Fdrive, neuron.Cm0, Qm0) # Determine patches location npatches, tpatch_on, tpatch_off = getPatchesLoc(t, states) - # Adding onset to all signals + # Adding onset to time vector and patches if t_plt['onset'] > 0.0: t = np.insert(t + t_plt['onset'], 0, 0.0) - for i in range(nvars): - vrs[i] = np.insert(vrs[i], 0, vrs[i][0]) tpatch_on += t_plt['onset'] tpatch_off += t_plt['onset'] # Plotting if naxes == 1: _, ax = plt.subplots(figsize=(11, 4)) axes = [ax] else: _, axes = plt.subplots(naxes, 1, figsize=(11, min(3 * naxes, 9))) - # Axes for i in range(naxes): + ax = axes[i] - if positions[i] < naxes - 1: + ax_pltvars = [pltvars[j] for j in vars_dict[labels[i]]] + nvars = len(ax_pltvars) + + # X-axis + if i < naxes - 1: ax.get_xaxis().set_ticklabels([]) else: ax.set_xlabel('${}\ ({})$'.format(t_plt['label'], t_plt['unit']), fontsize=fs) for tick in ax.xaxis.get_major_ticks(): tick.label.set_fontsize(fs) + + # Y-axis + if ax_pltvars[0]['unit']: + ax.set_ylabel('${}\ ({})$'.format(labels[i], ax_pltvars[0]['unit']), + fontsize=fs) + else: + ax.set_ylabel('${}$'.format(labels[i]), fontsize=fs) + if 'min' in ax_pltvars[0] and 'max' in ax_pltvars[0]: + ax_min = min([ap['min'] for ap in ax_pltvars]) + ax_max = max([ap['max'] for ap in ax_pltvars]) + ax.set_ylim(ax_min, ax_max) ax.locator_params(axis='y', nbins=2) for tick in ax.yaxis.get_major_ticks(): tick.label.set_fontsize(fs) - # Time series - icolor = 0 - for i in range(nvars): - pltvar = y_pltvars[i] - ax = axes[positions[i]] - if 'constant' in pltvar: - ax.plot(t * t_plt['factor'], vrs[i] * pltvar['factor'], '--', c='black', lw=4) - else: - ax.plot(t * t_plt['factor'], vrs[i] * pltvar['factor'], - c='C{}'.format(icolor), lw=4) - if 'min' in pltvar and 'max' in pltvar: - ax.set_ylim(pltvar['min'], pltvar['max']) - if pltvar['unit']: - ax.set_ylabel('${}\ ({})$'.format(pltvar['label'], pltvar['unit']), - fontsize=fs) + # Time series + icolor = 0 + for j in range(nvars): + + # Extract variable + pltvar = ax_pltvars[j] + if 'alias' in pltvar: + var = eval(pltvar['alias']) + elif 'key' in pltvar: + var = data[pltvar['key']] + elif 'constant' in pltvar: + var = eval(pltvar['constant']) * np.ones(t.size) + else: + var = data[vars_dict[labels[i]][j]] + if t_plt['onset'] > 0.0: + var = np.insert(var, 0, var[0]) + + # Plot variable + if 'constant' in pltvar: + ax.plot(t * t_plt['factor'], var * pltvar['factor'], '--', c='black', lw=lw, + label='${}$'.format(pltvar['label'])) else: - ax.set_ylabel('${}$'.format(pltvar['label']), fontsize=fs) - icolor += 1 + ax.plot(t * t_plt['factor'], var * pltvar['factor'], + c='C{}'.format(icolor), lw=lw, label='${}$'.format(pltvar['label'])) + icolor += 1 - # Patches - if show_patches == 1: - for ax in axes: + # Patches + if show_patches == 1: (ybottom, ytop) = ax.get_ylim() for j in range(npatches): ax.add_patch(Rectangle((tpatch_on[j] * t_plt['factor'], ybottom), (tpatch_off[j] - tpatch_on[j]) * t_plt['factor'], ytop - ybottom, color='#8A8A8A', alpha=0.1)) + # Legend + if nvars > 1: + ax.legend(fontsize=fs, loc=7, ncol=nvars // 5 + 1) + plt.tight_layout() # Save figure if needed (automatic or checked) if plt_save == 1: if ask_before_save == 1: plt_filename = SaveFigDialog(directory, '{}_{}.{}'.format(filecode, tag, fig_ext)) else: plt_filename = '{}/{}_{}.{}'.format(directory, filecode, tag, fig_ext) if plt_filename: plt.savefig(plt_filename) print('Saving figure as "{}"'.format(plt_filename)) plt.close() # Show all plots if needed if plt_show == 1: plt.show() diff --git a/PointNICE/plt/pltvars.py b/PointNICE/plt/pltvars.py index 552cf32..77c4aff 100644 --- a/PointNICE/plt/pltvars.py +++ b/PointNICE/plt/pltvars.py @@ -1,321 +1,392 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-21 14:33:36 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-08-23 14:52:26 +# @Last Modified time: 2017-08-24 17:34:28 ''' Dictionary of plotting settings for output variables of the model. ''' pltvars = { 't': { 'desc': 'time', 'label': 'time', 'unit': 'ms', 'factor': 1e3, 'onset': 3e-3 }, 'Z': { 'desc': 'leaflets deflection', 'label': 'Z', 'unit': 'nm', 'factor': 1e9, 'min': -1.0, 'max': 10.0 }, 'ng': { 'desc': 'gas content', 'label': 'gas', 'unit': '10^{-22}\ mol', 'factor': 1e22, 'min': 1.0, 'max': 15.0 }, 'Pac': { 'desc': 'acoustic pressure', 'label': 'P_{AC}', 'unit': 'kPa', 'factor': 1e-3, 'alias': 'bls.Pacoustic(t, data["Adrive"] * states, Fdrive)' }, 'Pmavg': { 'desc': 'average intermolecular pressure', 'label': 'P_M', 'unit': 'kPa', 'factor': 1e-3, 'alias': 'bls.PMavgpred(data["Z"])' }, 'Telastic': { 'desc': 'leaflet elastic tension', 'label': 'T_E', 'unit': 'mN/m', 'factor': 1e3, 'alias': 'bls.TEleaflet(data["Z"])' }, 'Qm': { 'desc': 'charge density', 'label': 'Q_m', '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': 'np.array([bls.Capct(ZZ) for ZZ in data["Z"]])' }, 'Vm': { 'desc': 'membrane potential', 'label': 'V_m', 'unit': 'mV', - 'factor': 1e3, - 'alias': 'data["Qm"] / np.array([bls.Capct(ZZ) for ZZ in data["Z"]])' - }, - - - 'iL': { - 'desc': 'leakage current', - 'label': 'I_L', - 'unit': 'mA/cm^2', 'factor': 1, - 'alias': 'neuron.currL(data["Qm"] * 1e3 / np.array([bls.Capct(ZZ) for ZZ in data["Z"]]))' }, - '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 }, 'n': { 'desc': 'iK activation gate opening', 'label': 'n-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, '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 }, '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 - data["O"] - data["C"]' }, + 'O + 2OL': { + 'desc': 'iH activation gate relative conductance', + 'label': 'O\ +\ 2O_L', + 'unit': None, + 'factor': 1, + 'min': -0.1, + 'max': 2.1, + 'alias': 'data["O"] + 2 * (1 - data["O"] - data["C"])' + }, + '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_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_{Na^{2+}}', 'unit': 'arb', 'factor': 1 }, 'VL': { 'constant': 'neuron.VL', 'desc': 'non-specific leakage current resting potential', 'label': 'A_{Na^{2+}}', 'unit': 'mV', 'factor': 1e0 }, + 'iL': { + 'desc': 'leakage current', + 'label': 'I_L', + 'unit': 'A/m^2', + 'factor': 1e-3, + 'alias': 'neuron.currL(data["Vm"])' + }, + + 'iNa': { + 'desc': 'Sodium current', + 'label': 'I_{Na}', + 'unit': 'A/m^2', + 'factor': 1e-3, + 'alias': 'neuron.currNa(data["m"], data["h"], data["Vm"])' + }, + + 'iK': { + 'desc': 'delayed-rectifier Potassium current', + 'label': 'I_K', + 'unit': 'A/m^2', + 'factor': 1e-3, + 'alias': 'neuron.currK(data["n"], data["Vm"])' + }, + + 'iM': { + 'desc': 'slow non-inactivating Potassium current', + 'label': 'I_M', + 'unit': 'A/m^2', + 'factor': 1e-3, + 'alias': 'neuron.currM(data["p"], data["Vm"])' + }, + + 'iT': { + 'desc': 'low-threshold Calcium current', + 'label': 'I_T', + 'unit': 'A/m^2', + 'factor': 1e-3, + 'alias': 'neuron.currCa(data["s"], data["u"], data["Vm"])' + }, + + 'iTs': { + 'desc': 'low-threshold Calcium current', + 'label': 'I_{TS}', + 'unit': 'A/m^2', + 'factor': 1e-3, + 'alias': 'neuron.currCa(data["s"], data["u"], data["Vm"])' + }, + + 'iH': { + 'desc': 'hyperpolarization-activated cationic current', + 'label': 'I_h', + 'unit': 'A/m^2', + 'factor': 1e-3, + 'alias': 'neuron.currH(data["O"], data["C"], data["Vm"])' + }, + + 'iKL': { + 'desc': 'leakage Potassium current', + 'label': 'I_{KL}', + 'unit': 'A/m^2', + 'factor': 1e-3, + 'alias': 'neuron.currKL(data["Vm"])' + }, + + 'iNet': { + 'desc': 'net current', + 'label': 'I_{net}', + 'unit': 'A/m^2', + 'factor': 1e-3, + 'alias': 'neuron.currNet(data["Vm"], neuron_states)' + }, + 'Veff': { 'key': 'V', 'desc': 'effective membrane potential', 'label': 'V_{m, eff}', 'unit': 'mV', 'factor': 1e0 }, 'alpham': { 'desc': 'iNa m-gate activation rate', 'label': '\\alpha_{m,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'betam': { 'desc': 'iNa m-gate inactivation rate', 'label': '\\beta_{m,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'alphah': { 'desc': 'iNa h-gate activation rate', 'label': '\\alpha_{h,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'betah': { 'desc': 'iNa h-gate inactivation rate', 'label': '\\beta_{h,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'alphan': { 'desc': 'iK n-gate activation rate', 'label': '\\alpha_{n,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'betan': { 'desc': 'iK n-gate inactivation rate', 'label': '\\beta_{n,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'alphap': { 'desc': 'iM p-gate activation rate', 'label': '\\alpha_{p,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'betap': { 'desc': 'iM p-gate inactivation rate', 'label': '\\beta_{p,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'alphas': { 'desc': 'iT s-gate activation rate', 'label': '\\alpha_{s,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'betas': { 'desc': 'iT s-gate inactivation rate', 'label': '\\beta_{s,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'alphau': { 'desc': 'iT u-gate activation rate', 'label': '\\alpha_{u,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'betau': { 'desc': 'iT u-gate inactivation rate', 'label': '\\beta_{u,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 } } diff --git a/PointNICE/solvers/SolverElec.py b/PointNICE/solvers/SolverElec.py index de7e248..7ad4fa8 100644 --- a/PointNICE/solvers/SolverElec.py +++ b/PointNICE/solvers/SolverElec.py @@ -1,147 +1,147 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-09-29 16:16:19 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-08-03 15:09:53 +# @Last Modified time: 2017-08-24 13:14:22 import logging import numpy as np import scipy.integrate as integrate # Get package logger logger = logging.getLogger('PointNICE') class SolverElec: def __init__(self): # Do nothing logger.info('Elec solver initialization') def eqHH(self, _, y, channel_mech, Iinj): ''' Compute the derivatives of a HH system variables for a specific value of injected current. :param t: time value (s, unused) :param y: vector of HH system variables at time t :param channel_mech: channels mechanism object :param Iinj: injected current (mA/m2) :return: vector of HH system derivatives at time t ''' Vm, *states = y Iionic = channel_mech.currNet(Vm, states) # mA/m2 dVmdt = (- Iionic + Iinj) / channel_mech.Cm0 # mV/s dstates = channel_mech.derStates(Vm, states) return [dVmdt, *dstates] def runSim(self, channel_mech, Astim, tstim, toffset, tonset=10e-3): ''' Compute solutions of a neuron's HH system for a specific set of electrical stimulation parameters, using a classic integration scheme. :param channel_mech: channels mechanism object :param Astim: pulse amplitude (mA/m2) :param tstim: pulse duration (s) :param toffset: offset duration (s) :param tonset: onset duration (s) :return: 2-tuple with the time profile and solution matrix ''' # Set time vector ttot = tonset + tstim + toffset dt = 1e-4 # s nsamples = int(np.round(ttot / dt)) t = np.linspace(0.0, ttot, nsamples) - tonset # Set pulse vector n_onset = int(np.round(tonset / dt)) n_stim = int(np.round(tstim / dt)) n_offset = int(np.round(toffset / dt)) pulse = np.concatenate((np.zeros(n_onset), Astim * np.ones(n_stim), np.zeros(n_offset))) # Create solver solver = integrate.ode(self.eqHH) solver.set_integrator('lsoda', nsteps=1000) # Set initial conditions y0 = [channel_mech.Vm0, *channel_mech.states0] nvar = len(y0) # Run simulation - y = np.empty((nsamples - 1, nvar)) + y = np.empty((nvar, nsamples - 1)) solver.set_initial_value(y0, t[0]) k = 1 while solver.successful() and k <= nsamples - 1: solver.set_f_params(channel_mech, pulse[k]) solver.integrate(t[k]) - y[k - 1, :] = solver.y + y[:, k - 1] = solver.y k += 1 - y = np.concatenate((np.atleast_2d(y0), y), axis=0) + y = np.concatenate((np.atleast_2d(y0).T, y), axis=1) return (t, y) def eqHH_VClamp(self, _, y, channel_mech, Vc): ''' Compute the derivatives of a HH system variables for a specific value of clamped voltage. :param t: time value (s, unused) :param y: vector of HH system variables at time t :param channel_mech: channels mechanism object :param Vc: clamped voltage (mV) :return: vector of HH system derivatives at time t ''' return channel_mech.derStates(Vc, y) def runVClamp(self, channel_mech, Vclamp, tclamp, toffset, tonset=10e-3): ''' Compute solutions of a neuron's HH system for a specific set of voltage clamp parameters, using a classic integration scheme. :param channel_mech: channels mechanism object :param Vclamp: clamped voltage (mV) :param toffset: offset duration (s) :param tclamp: clamp duration (s) :param tonset: onset duration (s) :return: 2-tuple with the time profile and solution matrix ''' # Set time vector ttot = tonset + tclamp + toffset dt = 1e-4 # s nsamples = int(np.round(ttot / dt)) t = np.linspace(0.0, ttot, nsamples) - tonset # Set clamp vector n_onset = int(np.round(tonset / dt)) n_clamp = int(np.round(tclamp / dt)) n_offset = int(np.round(toffset / dt)) clamp = np.concatenate((np.zeros(n_onset), Vclamp * np.ones(n_clamp), np.zeros(n_offset))) # Create solver solver = integrate.ode(self.eqHH_VClamp) solver.set_integrator('lsoda', nsteps=1000) # Set initial conditions y0 = channel_mech.states0 nvar = len(y0) # Run simulation y = np.empty((nsamples - 1, nvar)) solver.set_initial_value(y0, t[0]) k = 1 while solver.successful() and k <= nsamples - 1: solver.set_f_params(channel_mech, clamp[k]) solver.integrate(t[k]) y[k - 1, :] = solver.y k += 1 y = np.concatenate((np.atleast_2d(y0), y), axis=0) return (t, y) diff --git a/PointNICE/solvers/simutils.py b/PointNICE/solvers/simutils.py index ca030a8..0f14f42 100644 --- a/PointNICE/solvers/simutils.py +++ b/PointNICE/solvers/simutils.py @@ -1,409 +1,564 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-22 14:33:04 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-08-23 16:23:51 +# @Last Modified time: 2017-08-24 16:49:21 """ Utility functions used in simulations """ +import os import time import logging import pickle +import shutil +import tkinter as tk +from tkinter import filedialog import numpy as np from openpyxl import load_workbook # from . import SolverUS from .SolverUS import SolverUS +from .SolverElec import SolverElec from ..constants import * from ..utils import detectSpikes # Get package logger logger = logging.getLogger('PointNICE') +# Naming and logging settings +ESTIM_CW_code = 'ESTIM_{}_{:.1f}mA_per_m2_{:.0f}ms' +ASTIM_CW_code = 'ASTIM_{}_CW_{:.0f}nm_{:.0f}kHz_{:.0f}kPa_{:.0f}ms_{}' +ASTIM_PW_code = 'ASTIM_{}_PW_{:.0f}nm_{:.0f}kHz_{:.0f}kPa_{:.0f}ms_PRF{:.2f}kHz_DF{:.2f}_{}' +ESTIM_CW_log = '%s neuron - E-STIM simulation %u/%u (A = %.1f mA/m2, t = %.1f ms)' +ASTIM_CW_log = ('%s neuron - CW A-STIM %s simulation %u/%u (a = %.1f nm, f = %.2f kHz, ' + 'A = %.2f kPa, t = %.2f ms') +ASTIM_PW_log = ('%s neuron - PW A-STIM %s simulation %u/%u (a = %.1f nm, f = %.2f kHz, ' + 'A = %.2f kPa, t = %.2f ms, PRF = %.2f kHz, DF = %.2f)') +ASTIM_CW_titration_log = ('%s neuron - CW A-STIM titration %u/%u (a = %.1f nm, f = %.2f kHz, ' + '%s = %.2f %s') +ASTIM_PW_titration_log = ('%s neuron - PW A-STIM titration %u/%u (a = %.1f nm, f = %.2f kHz, ' + '%s = %.2f %s, PRF = %.2f kHz, DF = %.2f)') + + + +def checkBatchLog(batch_type): + ''' Determine batch directory, and add a log file to the directory if it is absent. + + :param batch_type: name of the log file to search for + :return: 2-tuple with full paths to batch directory and log file + ''' + + # Get batch directory from user + root = tk.Tk() + root.withdraw() + batch_dir = filedialog.askdirectory() + assert batch_dir, 'No batch directory chosen' + + # Determine log template from batch type + if batch_type == 'mech': + logfile = 'log_mech.xlsx' + elif batch_type == 'A-STIM': + logfile = 'log_ASTIM.xlsx' + elif batch_type == 'E-STIM': + logfile = 'log_ESTIM.xlsx' + else: + raise ValueError('Unknown batch type', batch_type) + this_dir, _ = os.path.split(__file__) + parent_dir = os.path.abspath(os.path.join(this_dir, os.pardir)) + logsrc = parent_dir + '/templates/' + logfile + assert os.path.isfile(logsrc), 'template log file "{}" not found'.format(logsrc) + + # Copy template in batch directory if no appropriate log file + logdst = batch_dir + '/' + logfile + if not os.path.isfile(logdst): + shutil.copy2(logsrc, logdst) + + return (batch_dir, logdst) + + def createSimQueue(amps, durations, offsets, PRFs, DFs): ''' Create a serialized 2D array of all parameter combinations for a series of individual parameter sweeps, while avoiding repetition of CW protocols for a given PRF sweep. :param amps: list (or 1D-array) of acoustic amplitudes :param durations: list (or 1D-array) of stimulus durations :param offsets: list (or 1D-array) of stimulus offsets (paired with durations array) :param PRFs: list (or 1D-array) of pulse-repetition frequencies :param DFs: list (or 1D-array) of duty cycle values :return: 2D-array with (amplitude, duration, offset, PRF, DF) for each stimulation protocol ''' # Convert input to 1D-arrays amps = np.array(amps) durations = np.array(durations) offsets = np.array(offsets) PRFs = np.array(PRFs) DFs = np.array(DFs) # Create index arrays iamps = range(len(amps)) idurs = range(len(durations)) # Create empty output matrix queue = np.empty((1, 5)) # Continuous protocols if 1.0 in DFs: nCW = len(amps) * len(durations) arr1 = np.ones(nCW) iCW_queue = np.array(np.meshgrid(iamps, idurs)).T.reshape(nCW, 2) CW_queue = np.hstack((amps[iCW_queue[:, 0]], durations[iCW_queue[:, 1]], offsets[iCW_queue[:, 1]], PRFs.min() * arr1, arr1)) queue = np.vstack((queue, CW_queue)) # Pulsed protocols if np.any(DFs != 1.0): pulsed_DFs = DFs[DFs != 1.0] iPRFs = range(len(PRFs)) ipulsed_DFs = range(len(pulsed_DFs)) nPW = len(amps) * len(durations) * len(PRFs) * len(pulsed_DFs) iPW_queue = np.array(np.meshgrid(iamps, idurs, iPRFs, ipulsed_DFs)).T.reshape(nPW, 4) PW_queue = np.hstack((amps[iPW_queue[:, 0]], durations[iPW_queue[:, 1]], offsets[iPW_queue[:, 1]], PRFs[iPW_queue[:, 2]], pulsed_DFs[iPW_queue[:, 3]])) queue = np.vstack((queue, PW_queue)) # Return return queue[1:, :] def xlslog(filename, sheetname, data): """ Append log data on a new row to specific sheet of excel workbook. :param filename: absolute or relative path to the Excel workbook :param sheetname: name of the Excel spreadsheet to which data is appended :param data: data structure to be added to specific columns on a new row :return: boolean indicating success (1) or failure (0) of operation """ try: wb = load_workbook(filename) ws = wb[sheetname] keys = data.keys() i = 1 row_data = {} for k in keys: row_data[k] = data[k] i += 1 ws.append(row_data) wb.save(filename) return 1 except PermissionError: # If file cannot be accessed for writing because already opened logger.error('Cannot write to "%s". Close the file and type "Y"', filename) user_str = input() if user_str in ['y', 'Y']: return xlslog(filename, sheetname, data) else: return 0 -def runSimBatch(batch_dir, log_filepath, neurons, bls_params, geom, stim_params, sim_type): +def runEStimBatch(batch_dir, log_filepath, neurons, stim_params): + ''' Run batch E-STIM simulations of the system for various neuron types and + stimulation parameters. + + :param batch_dir: full path to output directory of batch + :param log_filepath: full path log file of batch + :param neurons: array of channel mechanisms + :param stim_params: dictionary containing sweeps for all stimulation parameters + ''' + + logger.info("Starting E-STIM simulation batch") + + # Generate simulations queue + sim_queue = createSimQueue(stim_params['amps'], stim_params['durations'], + stim_params['offsets'], stim_params['PRFs'], stim_params['DFs']) + nqueue = sim_queue.shape[0] + + # Initialize solver + solver = SolverElec() + + # Run simulations + nsims = len(neurons) * nqueue + simcount = 0 + filepaths = [] + for ch_mech in neurons: + try: + for i in range(nqueue): + simcount += 1 + Astim, tstim, toffset, _, _ = sim_queue[i, :] + + # Get date and time info + date_str = time.strftime("%Y.%m.%d") + daytime_str = time.strftime("%H:%M:%S") + + # Log and define naming + logger.info(ESTIM_CW_log, ch_mech.name, simcount, nsims, Astim, tstim * 1e3) + simcode = ESTIM_CW_code.format(ch_mech.name, Astim, tstim * 1e3) + + # Run simulation + tstart = time.time() + (t, y) = solver.runSim(ch_mech, Astim, tstim, toffset) + Vm, *channels = y + tcomp = time.time() - tstart + logger.info('completed in %.2f seconds', tcomp) + + # States vector + nsamples = t.size + nstim = int(np.round(nsamples * tstim / (tstim + toffset))) + noffset = nsamples - nstim + states = np.hstack((np.ones(nstim), np.zeros(noffset))) + + # Store data in dictionary + data = { + 'Astim': Astim, + 'tstim': tstim, + 'toffset': toffset, + 't': t, + 'states': states, + 'Vm': Vm + } + for j in range(len(ch_mech.states_names)): + data[ch_mech.states_names[j]] = channels[j] + + # Export data to PKL file + output_filepath = batch_dir + '/' + simcode + ".pkl" + with open(output_filepath, 'wb') as fh: + pickle.dump(data, fh) + logger.info('simulation data exported to "%s"', output_filepath) + filepaths.append(output_filepath) + + # Detect spikes on Vm signal + n_spikes, lat, sr = detectSpikes(t, Vm, SPIKE_MIN_VAMP, SPIKE_MIN_DT) + logger.info('%u spike%s detected', n_spikes, "s" if n_spikes > 1 else "") + + # Export key metrics to log file + log = { + 'A': date_str, + 'B': daytime_str, + 'C': ch_mech.name, + 'D': Astim, + 'E': tstim * 1e3, + 'F': 'N/A', + 'G': 'N/A', + 'H': t.size, + 'I': round(tcomp, 2), + 'J': n_spikes, + 'K': lat * 1e3 if isinstance(lat, float) else 'N/A', + 'L': sr * 1e-3 if isinstance(sr, float) else 'N/A' + } + + if xlslog(log_filepath, 'Data', log) == 1: + logger.info('log exported to "%s"', log_filepath) + else: + logger.error('log export to "%s" aborted', log_filepath) + + except AssertionError as err: + logger.error(err) + + return filepaths + + + + +def runAStimBatch(batch_dir, log_filepath, neurons, bls_params, geom, stim_params, + sim_type='effective'): ''' Run batch simulations of the system for various neuron types, sonophore and stimulation parameters. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param neurons: array of channel mechanisms :param bls_params: BLS biomechanical and biophysical parameters dictionary :param geom: BLS geometric constants dictionary :param stim_params: dictionary containing sweeps for all stimulation parameters :param sim_type: selected integration method ''' - # Define naming and logging settings - sim_str_CW = 'sim_{}_CW_{:.0f}nm_{:.0f}kHz_{:.0f}kPa_{:.0f}ms_{}' - sim_str_PW = 'sim_{}_PW_{:.0f}nm_{:.0f}kHz_{:.0f}kPa_{:.0f}ms_PRF{:.2f}kHz_DF{:.2f}_{}' - CW_log = ('%s neuron - CW %s simulation %u/%u (a = %.1f nm, f = %.2f kHz, A = %.2f kPa, ' - 't = %.1f ms)') - PW_log = ('%s neuron - PW %s simulation %u/%u (a = %.1f nm, f = %.2f kHz, A = %.2f kPa, ' - ' t = %.1f ms, PRF = %.2f kHz, DF = %.2f)') - - logger.info("Starting NICE simulation batch") + logger.info("Starting A-STIM simulation batch") + # Unpack geometrical parameters a = geom['a'] d = geom['d'] # Generate simulations queue sim_queue = createSimQueue(stim_params['amps'], stim_params['durations'], stim_params['offsets'], stim_params['PRFs'], stim_params['DFs']) nqueue = sim_queue.shape[0] # Run simulations - simcount = 0 nsims = len(neurons) * len(stim_params['freqs']) * nqueue + simcount = 0 + filepaths = [] for ch_mech in neurons: for Fdrive in stim_params['freqs']: try: # Create SolverUS instance (modulus of embedding tissue depends on frequency) solver = SolverUS(geom, bls_params, ch_mech, Fdrive) for i in range(nqueue): simcount += 1 Adrive, tstim, toffset, PRF, DF = sim_queue[i, :] # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") # Log and define naming if DF == 1.0: - logger.info(CW_log, ch_mech.name, sim_type, simcount, nsims, + logger.info(ASTIM_CW_log, ch_mech.name, sim_type, simcount, nsims, a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3) - simcode = sim_str_CW.format(ch_mech.name, a * 1e9, Fdrive * 1e-3, - Adrive * 1e-3, tstim * 1e3, sim_type) + simcode = ASTIM_CW_code.format(ch_mech.name, a * 1e9, Fdrive * 1e-3, + Adrive * 1e-3, tstim * 1e3, sim_type) else: - logger.info(PW_log, ch_mech.name, sim_type, simcount, nsims, a * 1e9, + logger.info(ASTIM_PW_log, ch_mech.name, sim_type, simcount, nsims, a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3, PRF * 1e-3, DF) - simcode = sim_str_PW.format(ch_mech.name, a * 1e9, Fdrive * 1e-3, - Adrive * 1e-3, tstim * 1e3, PRF * 1e-3, - DF, sim_type) + simcode = ASTIM_PW_code.format(ch_mech.name, a * 1e9, Fdrive * 1e-3, + Adrive * 1e-3, tstim * 1e3, PRF * 1e-3, + DF, sim_type) # Run simulation tstart = time.time() (t, y, states) = solver.runSim(ch_mech, Fdrive, Adrive, tstim, toffset, PRF, DF, sim_type) Z, ng, Qm, *channels = y U = np.insert(np.diff(Z) / np.diff(t), 0, 0.0) tcomp = time.time() - tstart logger.info('completed in %.2f seconds', tcomp) # Store data in dictionary bls_params['biophys']['Qm0'] = solver.Qm0 data = { 'a': a, 'd': d, 'params': bls_params, 'Fdrive': Fdrive, 'Adrive': Adrive, 'phi': np.pi, 'tstim': tstim, 'toffset': toffset, 'PRF': PRF, 'DF': DF, 't': t, 'states': states, 'U': U, 'Z': Z, 'ng': ng, - 'Qm': Qm + 'Qm': Qm, + 'Vm': Qm * 1e3 / np.array([solver.Capct(ZZ) for ZZ in Z]) } for j in range(len(ch_mech.states_names)): data[ch_mech.states_names[j]] = channels[j] # Export data to PKL file - datafile_name = batch_dir + '/' + simcode + ".pkl" - with open(datafile_name, 'wb') as fh: + output_filepath = batch_dir + '/' + simcode + ".pkl" + with open(output_filepath, 'wb') as fh: pickle.dump(data, fh) - logger.info('simulation data exported to "%s"', datafile_name) + logger.info('simulation data exported to "%s"', output_filepath) + filepaths.append(output_filepath) # Detect spikes on Qm signal n_spikes, lat, sr = detectSpikes(t, Qm, SPIKE_MIN_QAMP, SPIKE_MIN_DT) logger.info('%u spike%s detected', n_spikes, "s" if n_spikes > 1 else "") # Export key metrics to log file log = { 'A': date_str, 'B': daytime_str, 'C': ch_mech.name, 'D': a * 1e9, 'E': d * 1e6, 'F': Fdrive * 1e-3, 'G': Adrive * 1e-3, 'H': tstim * 1e3, 'I': PRF * 1e-3 if DF < 1 else 'N/A', 'J': DF, 'K': sim_type, 'L': t.size, 'M': round(tcomp, 2), 'N': n_spikes, 'O': lat * 1e3 if isinstance(lat, float) else 'N/A', 'P': sr * 1e-3 if isinstance(sr, float) else 'N/A' } if xlslog(log_filepath, 'Data', log) == 1: logger.info('log exported to "%s"', log_filepath) else: logger.error('log export to "%s" aborted', log_filepath) except AssertionError as err: logger.error(err) + return filepaths + def runTitrationBatch(batch_dir, log_filepath, neurons, bls_params, geom, stim_params): ''' Run batch titrations of the system for various neuron types, sonophore and stimulation parameters, to determine the threshold of a specific stimulus parameter for neural excitation. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param neurons: array of channel mechanisms :param bls_params: BLS biomechanical and biophysical parameters dictionary :param geom: BLS geometric constants dictionary :param stim_params: dictionary containing sweeps for all stimulation parameters ''' - # Define naming and logging settings - sim_str_CW = 'sim_{}_CW_{:.0f}nm_{:.0f}kHz_{:.0f}kPa_{:.0f}ms_{}' - sim_str_PW = 'sim_{}_PW_{:.0f}nm_{:.0f}kHz_{:.0f}kPa_{:.0f}ms_PRF{:.2f}kHz_DF{:.2f}_{}' - CW_log = ('%s neuron - CW titration %u/%u (a = %.1f nm, f = %.2f kHz, %s = %.2f %s') - PW_log = ('%s neuron - PW titration %u/%u (a = %.1f nm, f = %.2f kHz, %s = %.2f %s, ' - 'PRF = %.2f kHz, DF = %.2f)') - - logger.info("Starting NICE titration batch") + logger.info("Starting A-STIM titration batch") # Unpack geometrical parameters a = geom['a'] d = geom['d'] # Define default parameters sim_type = 'effective' offset = 30e-3 # Determine titration parameter (x) and titrations list A = {'name': 'A', 'factor': 1e-3, 'unit': 'kPa'} t = {'name': 't', 'factor': 1e3, 'unit': 'ms'} if 'durations' not in stim_params: varin = A varout = t titr_type = 'duration' sim_queue = createSimQueue(stim_params['amps'], [0.], [offset], stim_params['PRFs'], stim_params['DFs']) sim_queue = np.delete(sim_queue, 1, axis=1) elif 'amps' not in stim_params: varin = t varout = A titr_type = 'amplitude' sim_queue = createSimQueue([0.], stim_params['durations'], [offset] * len(stim_params['durations']), stim_params['PRFs'], stim_params['DFs']) sim_queue = np.delete(sim_queue, 0, axis=1) nqueue = sim_queue.shape[0] # Run titrations - simcount = 0 nsims = len(neurons) * len(stim_params['freqs']) * nqueue + simcount = 0 + filepaths = [] for ch_mech in neurons: for Fdrive in stim_params['freqs']: try: # Create SolverUS instance (modulus of embedding tissue depends on frequency) solver = SolverUS(geom, bls_params, ch_mech, Fdrive) for i in range(nqueue): simcount += 1 input_val, toffset, PRF, DF = sim_queue[i, :] # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") # Log and define naming if DF == 1.0: - logger.info(CW_log, ch_mech.name, simcount, nsims, a * 1e9, Fdrive * 1e-3, - varin['name'], input_val * varin['factor'], varin['unit']) + logger.info(ASTIM_CW_titration_log, ch_mech.name, simcount, nsims, a * 1e9, + Fdrive * 1e-3, varin['name'], input_val * varin['factor'], + varin['unit']) else: - logger.info(PW_log, ch_mech.name, simcount, nsims, a * 1e9, Fdrive * 1e-3, - varin['name'], input_val * varin['factor'], varin['unit'], - PRF * 1e-3, DF) + logger.info(ASTIM_PW_titration_log, ch_mech.name, simcount, nsims, a * 1e9, + Fdrive * 1e-3, varin['name'], input_val * varin['factor'], + varin['unit'], PRF * 1e-3, DF) # Run titration tstart = time.time() (output_thr, t, y, states, lat) = solver.titrate(ch_mech, Fdrive, input_val, toffset, PRF, DF, titr_type) Z, ng, Qm, *channels = y U = np.insert(np.diff(Z) / np.diff(t), 0, 0.0) tcomp = time.time() - tstart logger.info('completed in %.2f s, threshold = %.2f %s', tcomp, output_thr * varout['factor'], varout['unit']) # Sort input and output as amplitude and duration if titr_type == 'amplitude': Adrive = output_thr tstim = input_val elif titr_type == 'duration': tstim = output_thr Adrive = input_val # Define output naming if DF == 1.0: - sim_str_CW = 'sim_{}_CW_{:.0f}nm_{:.0f}kHz_{:.0f}kPa_{:.0f}ms_{}' - simcode = sim_str_CW.format(ch_mech.name, a * 1e9, Fdrive * 1e-3, - Adrive * 1e-3, tstim * 1e3, sim_type) + simcode = ASTIM_CW_code.format(ch_mech.name, a * 1e9, Fdrive * 1e-3, + Adrive * 1e-3, tstim * 1e3, sim_type) else: - simcode = sim_str_PW.format(ch_mech.name, a * 1e9, Fdrive * 1e-3, - Adrive * 1e-3, tstim * 1e3, PRF * 1e-3, - DF, sim_type) + simcode = ASTIM_PW_code.format(ch_mech.name, a * 1e9, Fdrive * 1e-3, + Adrive * 1e-3, tstim * 1e3, PRF * 1e-3, + DF, sim_type) # Store data in dictionary bls_params['biophys']['Qm0'] = solver.Qm0 data = { 'a': a, 'd': d, 'params': bls_params, 'Fdrive': Fdrive, 'Adrive': Adrive, 'phi': np.pi, 'tstim': tstim, 'toffset': toffset, 'PRF': PRF, 'DF': DF, 't': t, 'states': states, 'U': U, 'Z': Z, 'ng': ng, 'Qm': Qm } for j in range(len(ch_mech.states_names)): data[ch_mech.states_names[j]] = channels[j] # Export data to PKL file - datafile_name = batch_dir + '/' + simcode + ".pkl" - with open(datafile_name, 'wb') as fh: + output_filepath = batch_dir + '/' + simcode + ".pkl" + with open(output_filepath, 'wb') as fh: pickle.dump(data, fh) - logger.info('simulation data exported to "%s"', datafile_name) + logger.info('simulation data exported to "%s"', output_filepath) + filepaths.append(output_filepath) # Detect spikes on Qm signal n_spikes, lat, sr = detectSpikes(t, Qm, SPIKE_MIN_QAMP, SPIKE_MIN_DT) logger.info('%u spike%s detected', n_spikes, "s" if n_spikes > 1 else "") # Export key metrics to log file log = { 'A': date_str, 'B': daytime_str, 'C': ch_mech.name, 'D': a * 1e9, 'E': d * 1e6, 'F': Fdrive * 1e-3, 'G': Adrive * 1e-3, 'H': tstim * 1e3, 'I': PRF * 1e-3 if DF < 1 else 'N/A', 'J': DF, 'K': sim_type, 'L': t.size, 'M': round(tcomp, 2), 'N': n_spikes, 'O': lat * 1e3 if isinstance(lat, float) else 'N/A', 'P': sr * 1e-3 if isinstance(sr, float) else 'N/A' } if xlslog(log_filepath, 'Data', log) == 1: logger.info('log exported to "%s"', log_filepath) else: logger.error('log export to "%s" aborted', log_filepath) except AssertionError as err: logger.error(err) + + return filepaths diff --git a/PointNICE/utils.py b/PointNICE/utils.py index 530b541..3d70d84 100644 --- a/PointNICE/utils.py +++ b/PointNICE/utils.py @@ -1,410 +1,372 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-09-19 22:30:46 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-08-23 16:22:25 +# @Last Modified time: 2017-08-24 14:09:42 """ Definition of generic utility functions used in other modules """ from enum import Enum from functools import partial import os -import shutil import logging import tkinter as tk from tkinter import filedialog from openpyxl import load_workbook import numpy as np import yaml # Get package logger logger = logging.getLogger('PointNICE') class PmCompMethod(Enum): """ Enum: types of computation method for the intermolecular pressure """ direct = 1 predict = 2 def LoadParamsFile(filename): """ Load a dictionary of parameters for the BLS model from an external yaml file. :param filename: name of the input file :return: parameters dictionary """ logger.info('Loading parameters from "%s"', filename) with open(filename, 'r') as f: stream = f.read() params = yaml.load(stream) return ParseNestedDict(params) LoadParams = partial(LoadParamsFile, filename=os.path.split(__file__)[0] + '/params.yaml') def getLookupDir(): """ Return the location of the directory holding lookups files. :return: absolute path to the directory """ this_dir, _ = os.path.split(__file__) return this_dir + '/lookups' def ParseNestedDict(dict_in): """ Loop through a nested dictionary object and convert all string fields to floats. """ for key, value in dict_in.items(): if isinstance(value, dict): # If value itself is dictionary dict_in[key] = ParseNestedDict(value) elif isinstance(dict_in[key], str): dict_in[key] = float(dict_in[key]) return dict_in def OpenFilesDialog(filetype, dirname=''): """ Open a FileOpenDialogBox to select one or multiple file. The default directory and file type are given. :param dirname: default directory :param filetype: default file type :return: tuple of full paths to the chosen filenames """ root = tk.Tk() root.withdraw() filenames = filedialog.askopenfilenames(filetypes=[(filetype + " files", '.' + filetype)], initialdir=dirname) if filenames: par_dir = os.path.abspath(os.path.join(filenames[0], os.pardir)) else: par_dir = None return (filenames, par_dir) def ImportExcelCol(filename, sheetname, colstr, startrow): """ Load a specific column of an excel workbook as a numpy array. :param filename: absolute or relative path to the Excel workbook :param sheetname: name of the Excel spreadsheet to which data is appended :param colstr: string of the column to import :param startrow: index of the first row to consider :return: 1D numpy array with the column data """ wb = load_workbook(filename, read_only=True) ws = wb.get_sheet_by_name(sheetname) range_start_str = colstr + str(startrow) range_stop_str = colstr + str(ws.max_row) tmp = np.array([[i.value for i in j] for j in ws[range_start_str:range_stop_str]]) return tmp[:, 0] def ConstructMatrix(serialized_inputA, serialized_inputB, serialized_output): """ Construct a 2D output matrix from serialized input. :param serialized_inputA: serialized input variable A :param serialized_inputB: serialized input variable B :param serialized_output: serialized output variable :return: 4-tuple with vectors of unique values of A (m) and B (n), output variable 2D matrix (m,n) and number of holes in the matrix """ As = np.unique(serialized_inputA) Bs = np.unique(serialized_inputB) nA = As.size nB = Bs.size output = np.zeros((nA, nB)) output[:] = np.NAN nholes = 0 for i in range(nA): iA = np.where(serialized_inputA == As[i]) for j in range(nB): iB = np.where(serialized_inputB == Bs[j]) iMatch = np.intersect1d(iA, iB) if iMatch.size == 0: nholes += 1 elif iMatch.size > 1: logger.warning('Identical serialized inputs with values (%f, %f)', As[i], Bs[j]) else: iMatch = iMatch[0] output[i, j] = serialized_output[iMatch] return (As, Bs, output, nholes) def rmse(x1, x2): """ Compute the root mean square error between two 1D arrays """ return np.sqrt(((x1 - x2) ** 2).mean()) def rsquared(x1, x2): ''' compute the R-squared coefficient between two 1D arrays ''' residuals = x1 - x2 ss_res = np.sum(residuals**2) ss_tot = np.sum((x1 - np.mean(x1))**2) return 1 - (ss_res / ss_tot) def DownSample(t_dense, y, nsparse): """ Decimate periodic signals to a specified number of samples.""" if(y.ndim) > 1: nsignals = y.shape[0] else: nsignals = 1 y = np.array([y]) # determine time step and period of input signal T = t_dense[-1] - t_dense[0] dt_dense = t_dense[1] - t_dense[0] # resample time vector linearly t_ds = np.linspace(t_dense[0], t_dense[-1], nsparse) # create MAV window nmav = int(0.03 * T / dt_dense) if nmav % 2 == 0: nmav += 1 mav = np.ones(nmav) / nmav # determine signals padding npad = int((nmav - 1) / 2) # determine indexes of sampling on convolved signals ids = np.round(np.linspace(0, t_dense.size - 1, nsparse)).astype(int) y_ds = np.empty((nsignals, nsparse)) # loop through signals for i in range(nsignals): # pad, convolve and resample pad_left = y[i, -(npad + 2):-2] pad_right = y[i, 1:npad + 1] y_ext = np.concatenate((pad_left, y[i, :], pad_right), axis=0) y_mav = np.convolve(y_ext, mav, mode='valid') y_ds[i, :] = y_mav[ids] if nsignals == 1: y_ds = y_ds[0, :] return (t_ds, y_ds) def Pressure2Intensity(p, rho, c): """ Return the spatial peak, pulse average acoustic intensity (ISPPA) associated with the specified pressure amplitude. :param p: pressure amplitude (Pa) :param rho: medium density (kg/m3) :param c: speed of sound in medium (m/s) :return: spatial peak, pulse average acoustic intensity (W/m2) """ return p**2 / (2 * rho * c) def Intensity2Pressure(I, rho, c): """ Return the pressure amplitude associated with the specified spatial peak, pulse average acoustic intensity (ISPPA). :param I: spatial peak, pulse average acoustic intensity (W/m2) :param rho: medium density (kg/m3) :param c: speed of sound in medium (m/s) :return: pressure amplitude (Pa) """ return np.sqrt(2 * rho * c * I) def find_nearest(array, value): ''' Find nearest element in 1D array. ''' idx = (np.abs(array - value)).argmin() return (idx, array[idx]) def rescale(x, lb, ub, lb_new=0, ub_new=1): ''' Rescale a value to a specific interval by linear transformation. ''' xnorm = (x - lb) / (ub - lb) return xnorm * (ub_new - lb_new) + lb_new def printPct(pct, precision): print(('{:.' + str(precision) + 'f}%').format(pct), end='', flush=True) print('\r' * (precision + 3), end='') def LennardJones(x, beta, alpha, C, m, n): """ Generic expression of a Lennard-Jones function, adapted for the context of symmetric deflection (distance = 2x). :param x: deflection (i.e. half-distance) :param beta: x-shifting factor :param alpha: x-scaling factor :param C: y-scaling factor :param m: exponent of the repulsion term :param n: exponent of the attraction term :return: Lennard-Jones potential at given distance (2x) """ return C * (np.power((alpha / (2 * x + beta)), m) - np.power((alpha / (2 * x + beta)), n)) -def CheckBatchLog(batch_type): - ''' Determine batch directory, and add a log file to the directory if it is absent. - - :param batch_type: name of the log file to search for - :return: 2-tuple with full paths to batch directory and log file - ''' - - # Get batch directory from user - root = tk.Tk() - root.withdraw() - batch_dir = filedialog.askdirectory() - assert batch_dir, 'No batch directory chosen' - - # Check presence of log file in batch directory - logdst = batch_dir + '/log.xlsx' - log_in_dir = os.path.isfile(logdst) - - # If no log file, copy template in directory - if not log_in_dir: - - # Determine log template from batch type - if batch_type == 'mech': - logfile = 'log_mech.xlsx' - elif batch_type == 'elec': - logfile = 'log_elec.xlsx' - else: - raise ValueError('Unknown batch type', batch_type) - this_dir, _ = os.path.split(__file__) - # par_dir = os.path.abspath(os.path.join(this_dir, os.pardir)) - logsrc = this_dir + '/templates/' + logfile - - # Copy template - shutil.copy2(logsrc, logdst) - - return (batch_dir, logdst) - - def detectPeaks(x, mph=None, mpd=1, threshold=0, edge='rising', kpsh=False, valley=False, ax=None): """ Detect peaks in data based on their amplitude and inter-peak distance. """ x = np.atleast_1d(x).astype('float64') if x.size < 3: return np.array([], dtype=int) if valley: x = -x # find indices of all peaks dx = x[1:] - x[:-1] # handle NaN's indnan = np.where(np.isnan(x))[0] if indnan.size: x[indnan] = np.inf dx[np.where(np.isnan(dx))[0]] = np.inf ine, ire, ife = np.array([[], [], []], dtype=int) if not edge: ine = np.where((np.hstack((dx, 0)) < 0) & (np.hstack((0, dx)) > 0))[0] else: if edge.lower() in ['rising', 'both']: ire = np.where((np.hstack((dx, 0)) <= 0) & (np.hstack((0, dx)) > 0))[0] if edge.lower() in ['falling', 'both']: ife = np.where((np.hstack((dx, 0)) < 0) & (np.hstack((0, dx)) >= 0))[0] ind = np.unique(np.hstack((ine, ire, ife))) # handle NaN's if ind.size and indnan.size: # NaN's and values close to NaN's cannot be peaks ind = ind[np.in1d(ind, np.unique(np.hstack((indnan, indnan - 1, indnan + 1))), invert=True)] # first and last values of x cannot be peaks if ind.size and ind[0] == 0: ind = ind[1:] if ind.size and ind[-1] == x.size - 1: ind = ind[:-1] # remove peaks < minimum peak height if ind.size and mph is not None: ind = ind[x[ind] >= mph] # remove peaks - neighbors < threshold if ind.size and threshold > 0: dx = np.min(np.vstack([x[ind] - x[ind - 1], x[ind] - x[ind + 1]]), axis=0) ind = np.delete(ind, np.where(dx < threshold)[0]) # detect small peaks closer than minimum peak distance if ind.size and mpd > 1: ind = ind[np.argsort(x[ind])][::-1] # sort ind by peak height idel = np.zeros(ind.size, dtype=bool) for i in range(ind.size): if not idel[i]: # keep peaks with the same height if kpsh is True idel = idel | (ind >= ind[i] - mpd) & (ind <= ind[i] + mpd) \ & (x[ind[i]] > x[ind] if kpsh else True) idel[i] = 0 # Keep current peak # remove the small peaks and sort back the indices by their occurrence ind = np.sort(ind[~idel]) return ind def detectPeaksTime(t, y, mph, mtd): """ Extension of the detectPeaks function to detect peaks in data based on their amplitude and time difference, with a non-uniform time vector. :param t: time vector (not necessarily uniform) :param y: signal :param mph: minimal peak height :param mtd: minimal time difference :return: array of peak indexes """ # Detect peaks on signal with no restriction on inter-peak distance raw_indexes = detectPeaks(y, mph, mpd=1) if raw_indexes.size > 0: # Filter relevant peaks with temporal distance n_raw = raw_indexes.size filtered_indexes = np.array([raw_indexes[0]]) for i in range(1, n_raw): i1 = filtered_indexes[-1] i2 = raw_indexes[i] if t[i2] - t[i1] < mtd: if y[i2] > y[i1]: filtered_indexes[-1] = i2 else: filtered_indexes = np.append(filtered_indexes, i2) # Return peak indexes return filtered_indexes else: return None def detectSpikes(t, Qm, min_amp, min_dt): ''' Detect spikes on a charge density signal, and return their number, latency and rate. :param t: time vector (s) :param Qm: charge density vector (C/m2) :param min_amp: minimal charge amplitude to detect spikes (C/m2) :param min_dt: minimal time interval between 2 spikes (s) :return: 3-tuple with number of spikes, latency (s) and spike rate (sp/s) ''' i_spikes = detectPeaksTime(t, Qm, min_amp, min_dt) if i_spikes is not None: latency = t[i_spikes[0]] # s n_spikes = i_spikes.size if n_spikes > 1: first_to_last_spike = t[i_spikes[-1]] - t[i_spikes[0]] # s spike_rate = n_spikes / first_to_last_spike # spikes/s else: spike_rate = 'N/A' else: latency = 'N/A' spike_rate = 'N/A' n_spikes = 0 return (n_spikes, latency, spike_rate) diff --git a/sim/ASTIM_batch.py b/sim/ASTIM_batch.py new file mode 100644 index 0000000..50b5a9e --- /dev/null +++ b/sim/ASTIM_batch.py @@ -0,0 +1,97 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# @Author: Theo Lemaire +# @Date: 2017-02-13 18:16:09 +# @Email: theo.lemaire@epfl.ch +# @Last Modified by: Theo Lemaire +# @Last Modified time: 2017-08-24 17:37:01 + +""" Run batch acoustic simulations of the NICE model. """ + +import os +import logging +import numpy as np +from PointNICE.solvers import checkBatchLog, runAStimBatch +from PointNICE.channels import * +from PointNICE.utils import LoadParams +from PointNICE.plt import plotBatch + +# Set logging options +logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%d/%m/%Y %H:%M:%S:') +logger = logging.getLogger('PointNICE') +logger.setLevel(logging.DEBUG) + +# BLS parameters +bls_params = LoadParams() + +# Geometry of BLS structure +a = 32e-9 # in-plane radius (m) +d = 0.0e-6 # embedding tissue thickness (m) +geom = {"a": a, "d": d} + +# Channels mechanisms +neurons = [ThalamoCortical()] + +# Stimulation parameters +stim_params = { + 'freqs': [3.5e5], # Hz + 'amps': [100e3], # Pa + 'durations': [50e-3], # s + 'PRFs': [1e2], # Hz + 'DFs': [1.0] +} +stim_params['offsets'] = [30e-3] * len(stim_params['durations']) # s + + +# Select output directory +try: + (batch_dir, log_filepath) = checkBatchLog('A-STIM') +except AssertionError as err: + logger.error(err) + quit() + +# Run A-STIM batch +pkl_filepaths = runAStimBatch(batch_dir, log_filepath, neurons, bls_params, geom, stim_params) +pkl_dir, _ = os.path.split(pkl_filepaths[0]) + + +vars_RS_FS = { + 'Z': ['Z'], + 'Q_m': ['Qm'], + 'i_{Na}\ kin.': ['m', 'h'], + 'i_K\ kin.': ['n'], + 'i_M\ kin.': ['p'], + 'curr.': ['iNa', 'iK', 'iM', 'iL', 'iNet'] +} + + +vars_LTS = { + 'Z': ['Z'], + 'Q_m': ['Qm'], + 'i_{Na}\ kin.': ['m', 'h'], + 'i_K\ kin.': ['n'], + 'i_M\ kin.': ['p'], + 'i_T\ kin.': ['s', 'u'], + 'curr.': ['iNa', 'iK', 'iM', 'iT', 'iL', 'iNet'] +} + +vars_RE = { + 'Z': ['Z'], + 'Q_m': ['Qm'], + 'i_{Na}\ kin.': ['m', 'h'], + 'i_K\ kin.': ['n'], + 'i_{TS}\ kin.': ['s', 'u'], + 'curr.': ['iNa', 'iK', 'iTs', 'iL', 'iNet'] +} + +vars_TC = { + 'Z': ['Z'], + 'Q_m': ['Qm'], + 'i_{Na}\ kin.': ['m', 'h'], + 'i_K\ kin.': ['n'], + 'i_{T}\ kin.': ['s', 'u'], + 'i_{H}\ kin.': ['O', 'OL', 'O + 2OL'], + 'curr.': ['iNa', 'iK', 'iT', 'iH', 'iKL', 'iL', 'iNet'] +} + +plotBatch(vars_TC, pkl_dir, pkl_filepaths) diff --git a/sim/ASTIM_mech_batch.py b/sim/ASTIM_mech_batch.py index 77989ea..b658ae9 100644 --- a/sim/ASTIM_mech_batch.py +++ b/sim/ASTIM_mech_batch.py @@ -1,141 +1,141 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-11-21 10:46:56 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-08-23 15:51:56 +# @Last Modified time: 2017-08-24 14:13:24 """ Run batch simulations of the NICE mechanical model with imposed charge densities """ import time import logging import pickle import numpy as np import PointNICE -from PointNICE.utils import LoadParams, CheckBatchLog -from PointNICE.solvers import xlslog +from PointNICE.utils import LoadParams +from PointNICE.solvers import xlslog, checkBatchLog # Set logging options logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%d/%m/%Y %H:%M:%S:') logger = logging.getLogger('PointNICE') logger.setLevel(logging.DEBUG) # Select output directory try: - (batch_dir, log_filepath) = CheckBatchLog('mech') + (batch_dir, log_filepath) = checkBatchLog('mech') except AssertionError as err: logger.error(err) quit() # Define naming and logging settings sim_str = 'sim_{:.0f}nm_{:.0f}kHz_{:.0f}kPa_{:.1f}nCcm2_mech' sim_log = 'simulation %u/%u (a = %.1f nm, d = %.1f um, f = %.2f kHz, A = %.2f kPa, Q = %.1f nC/cm2)' logger.info("Starting BLS simulation batch") # Load NICE parameters params = LoadParams() biomech = params['biomech'] ac_imp = biomech['rhoL'] * biomech['c'] # Rayl # Set geometry of NBLS structure a = 32e-9 # in-plane radius (m) d = 0.0e-6 # embedding tissue thickness (m) geom = {"a": a, "d": d} Cm0 = 1e-2 # membrane resting capacitance (F/m2) Qm0 = -54.0e-5 # membrane resting charge density (C/m2) # Set stimulation parameters freqs = [6.9e5] # Hz amps = [1.43e3] # Pa charges = np.linspace(0.0, 80.0, 81) * 1e-5 # C/m2 # Run simulations nsims = len(freqs) * len(amps) * len(charges) simcount = 0 for Fdrive in freqs: try: bls = PointNICE.BilayerSonophore(geom, params, Fdrive, Cm0, Qm0) # Create SolverUS instance (compression modulus of embedding tissue depends on frequency) for Adrive in amps: for Qm in charges: simcount += 1 # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") # Log to console logger.info(sim_log, simcount, nsims, a * 1e9, d * 1e6, Fdrive * 1e-3, Adrive * 1e-3, Qm * 1e5) # Run simulation tstart = time.time() (t, y, states) = bls.runMech(Fdrive, Adrive, Qm) (Z, ng) = y U = np.insert(np.diff(Z) / np.diff(t), 0, 0.0) tcomp = time.time() - tstart logger.info('completed in %.2f seconds', tcomp) # Export data to PKL file simcode = sim_str.format(a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, Qm * 1e5) datafile_name = batch_dir + '/' + simcode + ".pkl" data = {'a': a, 'd': d, 'params': params, 'Fdrive': Fdrive, 'Adrive': Adrive, 'phi': np.pi, 'Qm': Qm, 't': t, 'states': states, 'U': U, 'Z': Z, 'ng': ng} with open(datafile_name, 'wb') as fh: pickle.dump(data, fh) # Compute key output metrics Zmax = np.amax(Z) Zmin = np.amin(Z) Zabs_max = np.amax(np.abs([Zmin, Zmax])) eAmax = bls.arealstrain(Zabs_max) Tmax = bls.TEtot(Zabs_max) Pmmax = bls.PMavgpred(Zmin) ngmax = np.amax(ng) dUdtmax = np.amax(np.abs(np.diff(U) / np.diff(t)**2)) # Export key metrics to log xls file log = { 'A': date_str, 'B': daytime_str, 'C': a * 1e9, 'D': d * 1e6, 'E': Fdrive * 1e-3, 'F': Adrive * 1e-3, 'G': Qm * 1e5, 'H': t.size, 'I': tcomp, 'J': bls.kA + bls.kA_tissue, 'K': Zmax * 1e9, 'L': eAmax, 'M': Tmax * 1e3, 'N': (ngmax - bls.ng0) / bls.ng0, 'O': Pmmax * 1e-3, 'P': dUdtmax } success = xlslog(log_filepath, 'Data', log) if success == 1: logger.info('log exported to "%s"', log_filepath) else: logger.error('log export to "%s" aborted', log_filepath) except AssertionError as err: logger.error(err) diff --git a/sim/ASTIM_sim_batch.py b/sim/ASTIM_sim_batch.py deleted file mode 100644 index 2388c65..0000000 --- a/sim/ASTIM_sim_batch.py +++ /dev/null @@ -1,63 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -# @Author: Theo Lemaire -# @Date: 2017-02-13 18:16:09 -# @Email: theo.lemaire@epfl.ch -# @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-08-23 16:25:04 - -""" Run batch acoustic simulations of the NICE model. """ - -import logging -import numpy as np -from PointNICE.solvers import runSimBatch -from PointNICE.channels import * -from PointNICE.utils import LoadParams, CheckBatchLog - -# Set logging options -logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%d/%m/%Y %H:%M:%S:') -logger = logging.getLogger('PointNICE') -logger.setLevel(logging.DEBUG) - -# BLS parameters -bls_params = LoadParams() - -# Geometry of NBLS structure -a = 32e-9 # in-plane radius (m) -d = 0.0e-6 # embedding tissue thickness (m) -geom = {"a": a, "d": d} - -# Channels mechanisms -neurons = [CorticalRS()] - -# Stimulation parameters -stim_params = { - 'freqs': [3.5e5], # Hz - 'amps': [100e3], # Pa - 'durations': [50e-3], # s - 'PRFs': [1e2], # Hz - 'DFs': [1.0] -} -stim_params['offsets'] = [30e-3] * len(stim_params['durations']) # s - -# stim_params = { -# 'freqs': np.array([200, 400, 600, 800, 1000]) * 1e3, # Hz -# 'amps': np.array([10, 20, 40, 80, 150, 300, 600]) * 1e3, # Pa -# 'durs': np.array([20, 40, 60, 80, 100, 150, 200, 250, 300]) * 1e-3, # s -# 'PRFs': np.array([0.1, 0.2, 0.5, 1, 2, 5, 10]) * 1e3, # Hz -# 'DFs': np.array([1, 2, 5, 10, 25, 50, 75, 100]) / 100 -# } -# stim_params['offsets'] = 350e-3 - stim_params['durations'] # s - -# Simulation type -sim_type = 'effective' - -# Select output directory -try: - (batch_dir, log_filepath) = CheckBatchLog('elec') -except AssertionError as err: - logger.error(err) - quit() - -# Run simulation batch -runSimBatch(batch_dir, log_filepath, neurons, bls_params, geom, stim_params, sim_type) diff --git a/sim/ASTIM_titration_batch.py b/sim/ASTIM_titration_batch.py index 5d260c3..97fc120 100644 --- a/sim/ASTIM_titration_batch.py +++ b/sim/ASTIM_titration_batch.py @@ -1,50 +1,55 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-02-13 18:16:09 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-08-22 17:15:29 +# @Last Modified time: 2017-08-24 14:15:12 """ Run batch parameter titrations of the NICE model. """ +import os import logging import numpy as np -from PointNICE.solvers import runTitrationBatch +from PointNICE.solvers import checkBatchLog, runTitrationBatch from PointNICE.channels import * -from PointNICE.utils import LoadParams, CheckBatchLog +from PointNICE.utils import LoadParams +from PointNICE.plt import plotBatch # Set logging options logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%d/%m/%Y %H:%M:%S:') logger = logging.getLogger('PointNICE') logger.setLevel(logging.DEBUG) # BLS parameters bls_params = LoadParams() # Geometry of NBLS structure a = 32e-9 # in-plane radius (m) d = 0.0e-6 # embedding tissue thickness (m) geom = {"a": a, "d": d} # Channels mechanisms neurons = [CorticalRS()] # Stimulation parameters stim_params = { 'freqs': [3.5e5], # Hz # 'amps': [100e3], # Pa 'durations': [50e-3], # s 'PRFs': [1e2], # Hz 'DFs': [1.0] } # Select output directory try: - (batch_dir, log_filepath) = CheckBatchLog('elec') + (batch_dir, log_filepath) = checkBatchLog('A-STIM') except AssertionError as err: logger.error(err) quit() # Run titration batch -runTitrationBatch(batch_dir, log_filepath, neurons, bls_params, geom, stim_params) +pkl_filepaths = runTitrationBatch(batch_dir, log_filepath, neurons, bls_params, geom, stim_params) +pkl_dir, _ = os.path.split(pkl_filepaths[0]) + +plotBatch(['Qm'], [0], pkl_dir, pkl_filepaths) diff --git a/sim/ESTIM_batch.py b/sim/ESTIM_batch.py new file mode 100644 index 0000000..3d563a2 --- /dev/null +++ b/sim/ESTIM_batch.py @@ -0,0 +1,78 @@ +# -*- coding: utf-8 -*- +# @Author: Theo Lemaire +# @Date: 2017-08-24 11:55:07 +# @Last Modified by: Theo Lemaire +# @Last Modified time: 2017-08-24 17:45:53 + +""" Run batch electrical simulations of point-neuron models. """ + +import os +import logging +from PointNICE.solvers import checkBatchLog, runEStimBatch +from PointNICE.channels import * +from PointNICE.plt import plotBatch + +# Set logging options +logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%d/%m/%Y %H:%M:%S:') +logger = logging.getLogger('PointNICE') +logger.setLevel(logging.DEBUG) + +# Channels mechanisms +neurons = [ThalamoCortical()] + +# Stimulation parameters +stim_params = { + 'amps': [20.0], # mA/m2 + 'durations': [0.5], # s + 'PRFs': [1e2], # Hz + 'DFs': [1.0] +} +stim_params['offsets'] = [1.0] * len(stim_params['durations']) # s + +# Select output directory +try: + (batch_dir, log_filepath) = checkBatchLog('E-STIM') +except AssertionError as err: + logger.error(err) + quit() + +# Run E-STIM batch +pkl_filepaths = runEStimBatch(batch_dir, log_filepath, neurons, stim_params) +pkl_dir, _ = os.path.split(pkl_filepaths[0]) + +vars_RS_FS = { + 'V_m': ['Vm'], + 'i_{Na}\ kin.': ['m', 'h'], + 'i_L\ kin.': ['n'], + 'i_M\ kin.': ['p'], + 'curr.': ['iNa', 'iK', 'iM', 'iL', 'iNet'] +} + +vars_LTS = { + 'V_m': ['Vm'], + 'i_{Na}\ kin.': ['m', 'h'], + 'i_K\ kin.': ['n'], + 'i_M\ kin.': ['p'], + 'i_T\ kin.': ['s', 'u'], + 'curr.': ['iNa', 'iK', 'iM', 'iT', 'iL', 'iNet'] +} + +vars_RE = { + 'V_m': ['Vm'], + 'i_{Na}\ kin.': ['m', 'h'], + 'i_K\ kin.': ['n'], + 'i_{TS}\ kin.': ['s', 'u'], + 'curr.': ['iNa', 'iK', 'iTs', 'iL', 'iNet'] +} + +vars_TC = { + 'V_m': ['Vm'], + 'i_{Na}\ kin.': ['m', 'h'], + 'i_K\ kin.': ['n'], + 'i_{T}\ kin.': ['s', 'u'], + 'i_{H}\ kin.': ['O', 'OL', 'O + 2OL'], + 'curr.': ['iNa', 'iK', 'iT', 'iH', 'iKL', 'iL', 'iNet'] +} + + +plotBatch(vars_TC, pkl_dir, pkl_filepaths, lw=2)