diff --git a/PointNICE/lookups/FS_lookups_a32.0nm.pkl b/PointNICE/lookups/old/FS_lookups_a32.0nm.pkl similarity index 100% rename from PointNICE/lookups/FS_lookups_a32.0nm.pkl rename to PointNICE/lookups/old/FS_lookups_a32.0nm.pkl diff --git a/PointNICE/lookups/LTS_lookups_a32.0nm.pkl b/PointNICE/lookups/old/LTS_lookups_a32.0nm.pkl similarity index 100% rename from PointNICE/lookups/LTS_lookups_a32.0nm.pkl rename to PointNICE/lookups/old/LTS_lookups_a32.0nm.pkl diff --git a/PointNICE/lookups/LeechP_lookups_a32.0nm.pkl b/PointNICE/lookups/old/LeechP_lookups_a32.0nm.pkl similarity index 100% rename from PointNICE/lookups/LeechP_lookups_a32.0nm.pkl rename to PointNICE/lookups/old/LeechP_lookups_a32.0nm.pkl diff --git a/PointNICE/lookups/LeechT_lookups_a32.0nm.pkl b/PointNICE/lookups/old/LeechT_lookups_a32.0nm.pkl similarity index 100% rename from PointNICE/lookups/LeechT_lookups_a32.0nm.pkl rename to PointNICE/lookups/old/LeechT_lookups_a32.0nm.pkl diff --git a/PointNICE/lookups/RE_lookups_a32.0nm.pkl b/PointNICE/lookups/old/RE_lookups_a32.0nm.pkl similarity index 100% rename from PointNICE/lookups/RE_lookups_a32.0nm.pkl rename to PointNICE/lookups/old/RE_lookups_a32.0nm.pkl diff --git a/PointNICE/lookups/RS_lookups_a32.0nm.pkl b/PointNICE/lookups/old/RS_lookups_a32.0nm.pkl similarity index 100% rename from PointNICE/lookups/RS_lookups_a32.0nm.pkl rename to PointNICE/lookups/old/RS_lookups_a32.0nm.pkl diff --git a/PointNICE/lookups/TC_lookups_a32.0nm.pkl b/PointNICE/lookups/old/TC_lookups_a32.0nm.pkl similarity index 100% rename from PointNICE/lookups/TC_lookups_a32.0nm.pkl rename to PointNICE/lookups/old/TC_lookups_a32.0nm.pkl diff --git a/PointNICE/neurons/cortical.py b/PointNICE/neurons/cortical.py index 912b015..19b6111 100644 --- a/PointNICE/neurons/cortical.py +++ b/PointNICE/neurons/cortical.py @@ -1,605 +1,608 @@ #!/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-28 17:54:21 +# @Last Modified time: 2018-03-13 15:02:02 ''' 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) def __init__(self): ''' Constructor of the class ''' # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 'p'] self.states0 = np.array([]) # Names of the different coefficients to be averaged in a lookup table. self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphap', 'betap'] + # Charge interval bounds for lookup creation + self.Qbounds = (np.round(self.Vm0 - 10.0) * 1e-5, 50.0e-5) + 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, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) m, h, n, p = states dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dpdt = rates[6] * (1 - p) - rates[7] * p return [dmdt, dhdt, dndt, dpdt] class CorticalRS(Cortical): ''' Specific membrane channel dynamics of a cortical regular spiking, excitatory pyramidal neuron. Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* ''' # Name of channel mechanism name = 'RS' # Cell-specific biophysical parameters Vm0 = -71.9 # Cell membrane resting potential (mV) GNaMax = 560.0 # Max. conductance of Sodium current (S/m^2) GKMax = 60.0 # Max. conductance of delayed Potassium current (S/m^2) GMMax = 0.75 # Max. conductance of slow non-inactivating Potassium current (S/m^2) 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) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_K\ kin.': ['n'], 'i_M\ kin.': ['p'], 'I': ['iNa', 'iK', 'iM', 'iL', 'iNet'] } 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) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_K\ kin.': ['n'], 'i_M\ kin.': ['p'], 'I': ['iNa', 'iK', 'iM', 'iL', 'iNet'] } 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) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_K\ kin.': ['n'], 'i_M\ kin.': ['p'], 'i_T\ kin.': ['s', 'u'], 'I': ['iNa', 'iK', 'iM', 'iT', 'iL', 'iNet'] } 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, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' # Unpack input states *NaK_states, s, u = states # Call parent method to compute Sodium and Potassium channels states derivatives NaK_dstates = super().derStatesEff(Qm, NaK_states, interp_data) # Compute Calcium channels states derivatives Ca_rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names[8:]]) dsdt = Ca_rates[0] * (1 - s) - Ca_rates[1] * s dudt = Ca_rates[2] * (1 - u) - Ca_rates[3] * u # Merge all states derivatives and return return NaK_dstates + [dsdt, dudt] diff --git a/PointNICE/neurons/leech.py b/PointNICE/neurons/leech.py index a2e8064..9592c93 100644 --- a/PointNICE/neurons/leech.py +++ b/PointNICE/neurons/leech.py @@ -1,1074 +1,1080 @@ #!/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: 2018-03-07 14:46:54 +# @Last Modified time: 2018-03-13 15:04:00 ''' 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.58 # Cell membrane resting potential (mV) VNa = 45.0 # Sodium Nernst potential (mV) VK = -62.0 # Potassium Nernst potential (mV) VCa = 60.0 # Calcium Nernst potential (mV) 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) # Original conversion constants from inward ion current (nA) to build-up of # intracellular ion concentration (arb.) K_Na_original = 0.016 # iNa to intracellular [Na+] K_Ca_original = 0.1 # iCa to intracellular [Ca2+] # Constants needed to convert K from original model (soma compartment) # to current model (point-neuron) surface = 6434.0e-12 # surface of cell assumed as a single soma (m2) curr_factor = 1e6 # mA to nA # Time constants for the removal of ions from intracellular pools (s) tau_Na_removal = 16.0 # Na+ removal tau_Ca_removal = 1.25 # Ca2+ removal # Time constants for the iPumpNa and iKCa currents activation # from specific intracellular ions (s) tau_PumpNa_act = 0.1 # iPumpNa activation from intracellular Na+ tau_KCa_act = 0.01 # iKCa activation from intracellular Ca2+ # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h', 'm3h'], 'i_K\ kin.': ['n'], 'i_{Ca}\ kin.': ['s'], 'pools': ['C_Na_arb', 'C_Na_arb_activation', 'C_Ca_arb', 'C_Ca_arb_activation'], 'I': ['iNa', 'iK', 'iCa', 'iKCa', 'iPumpNa', 'iL', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'C_Na', 'A_Na', 'C_Ca', 'A_Ca'] self.states0 = np.array([]) # Names of the channels effective coefficients self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas'] self.K_Na = self.K_Na_original * self.surface * self.curr_factor self.K_Ca = self.K_Ca_original * self.surface * self.curr_factor # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) + # Charge interval bounds for lookup creation + self.Qbounds = (np.round(self.Vm0 - 10.0) * 1e-5, 50.0e-5) + # ----------------- 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.VK) 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): ''' Concrete implementation of the abstract API method. ''' 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): ''' Concrete implementation of the abstract API method. ''' # Standard gating dynamics: Solve the equation dx/dt = 0 at Vm for each x-state meq = self.minf(Vm) heq = self.hinf(Vm) neq = self.ninf(Vm) seq = self.sinf(Vm) # PumpNa pool concentration and activation steady-state INa_eq = self.currNa(meq, heq, Vm) CNa_eq = self.K_Na * (-INa_eq) ANa_eq = CNa_eq # KCa current pool concentration and activation steady-state ICa_eq = self.currCa(seq, Vm) CCa_eq = self.K_Ca * (-ICa_eq) ACa_eq = CCa_eq return np.array([meq, heq, neq, seq, CNa_eq, ANa_eq, CCa_eq, ACa_eq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack states m, h, n, s, C_Na, A_Na, C_Ca, A_Ca = states # Standard gating states derivatives dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) # PumpNa current pool concentration and activation state I_Na = self.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): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants Tm = self.taum minf = self.minf(Vm) am_avg = np.mean(minf / Tm) bm_avg = np.mean(1 / Tm) - am_avg Th = self.tauh(Vm) hinf = self.hinf(Vm) ah_avg = np.mean(hinf / Th) bh_avg = np.mean(1 / Th) - ah_avg Tn = self.taun(Vm) ninf = self.ninf(Vm) an_avg = np.mean(ninf / Tn) bn_avg = np.mean(1 / Tn) - an_avg Ts = self.taus sinf = self.sinf(Vm) as_avg = np.mean(sinf / Ts) bs_avg = np.mean(1 / Ts) - as_avg # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) Vmeff = np.interp(Qm, interp_data['Q'], interp_data['V']) # Unpack states m, h, n, s, C_Na, A_Na, C_Ca, A_Ca = states # Standard gating states derivatives dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s # PumpNa current pool concentration and activation state I_Na = self.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] class LeechMech(BaseMech): ''' Class defining the basic dynamics of Sodium, Potassium and Calcium channels for several neurons of the leech. Reference: *Baccus, S.A. (1998). Synaptic facilitation by reflected action potentials: enhancement of transmission when nerve impulses reverse direction at axon branch points. Proc. Natl. Acad. Sci. U.S.A. 95, 8345–8350.* ''' alphaC_sf = 1e-5 # Calcium activation rate constant scaling factor (M) betaC = 0.1e3 # beta rate for the open-probability of Ca2+-dependent Potassium channels (s-1) T = 293.15 # Room temperature (K) Rg = 8.314 # Universal gas constant (J.mol^-1.K^-1) Faraday = 9.6485e4 # Faraday constant (C/mol) def alpham(self, Vm): ''' Compute the alpha rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = -0.03 * (Vm + 28) / (np.exp(- (Vm + 28) / 15) - 1) # ms-1 return alpha * 1e3 # s-1 def betam(self, Vm): ''' Compute the beta rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 2.7 * np.exp(-(Vm + 53) / 18) # ms-1 return beta * 1e3 # s-1 def alphah(self, Vm): ''' Compute the alpha rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = 0.045 * np.exp(-(Vm + 58) / 18) # ms-1 return alpha * 1e3 # s-1 def betah(self, Vm): ''' Compute the beta rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) .. warning:: the original paper contains an error (multiplication) in the expression of this rate constant, corrected in the mod file on ModelDB (division). ''' beta = 0.72 / (np.exp(-(Vm + 23) / 14) + 1) # ms-1 return beta * 1e3 # s-1 def alphan(self, Vm): ''' Compute the alpha rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = -0.024 * (Vm - 17) / (np.exp(-(Vm - 17) / 8) - 1) # ms-1 return alpha * 1e3 # s-1 def betan(self, Vm): ''' Compute the beta rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 0.2 * np.exp(-(Vm + 48) / 35) # ms-1 return beta * 1e3 # s-1 def alphas(self, Vm): ''' Compute the alpha rate for the open-probability of Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = -1.5 * (Vm - 20) / (np.exp(-(Vm - 20) / 5) - 1) # ms-1 return alpha * 1e3 # s-1 def betas(self, Vm): ''' Compute the beta rate for the open-probability of Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 1.5 * np.exp(-(Vm + 25) / 10) # ms-1 return beta * 1e3 # s-1 def alphaC(self, C_Ca_in): ''' Compute the alpha rate for the open-probability of Calcium-dependent Potassium channels. :param C_Ca_in: intracellular Calcium concentration (M) :return: rate constant (s-1) ''' alpha = 0.1 * C_Ca_in / self.alphaC_sf # ms-1 return alpha * 1e3 # s-1 def derM(self, Vm, m): ''' Compute the evolution of the open-probability of Sodium channels. :param Vm: membrane potential (mV) :param m: open-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alpham(Vm) * (1 - m) - self.betam(Vm) * m def derH(self, Vm, h): ''' Compute the evolution of the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :param h: inactivation-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphah(Vm) * (1 - h) - self.betah(Vm) * h def derN(self, Vm, n): ''' Compute the evolution of the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :param n: open-probability of delayed-rectifier Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphan(Vm) * (1 - n) - self.betan(Vm) * n def derS(self, Vm, s): ''' Compute the evolution of the open-probability of Calcium channels. :param Vm: membrane potential (mV) :param s: open-probability of Calcium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphas(Vm) * (1 - s) - self.betas(Vm) * s def derC(self, c, C_Ca_in): ''' Compute the evolution of the open-probability of Calcium-dependent Potassium channels. :param c: open-probability of Calcium-dependent Potassium channels (prob) :param C_Ca_in: intracellular Calcium concentration (M) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphaC(C_Ca_in) * (1 - c) - self.betaC * c def currNa(self, m, h, Vm, C_Na_in): ''' 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) :param C_Na_in: intracellular Sodium concentration (M) :return: current per unit area (mA/m2) ''' GNa = self.GNaMax * m**4 * h VNa = self.nernst(self.Z_Na, C_Na_in, self.C_Na_out) # Sodium Nernst potential return GNa * (Vm - 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**2 return GK * (Vm - self.VK) def currCa(self, s, Vm, C_Ca_in): ''' Compute the inward Calcium current per unit area. :param s: open-probability of Calcium channels :param Vm: membrane potential (mV) :param C_Ca_in: intracellular Calcium concentration (M) :return: current per unit area (mA/m2) ''' GCa = self.GCaMax * s VCa = self.nernst(self.Z_Ca, C_Ca_in, self.C_Ca_out) # Calcium Nernst potential return GCa * (Vm - VCa) def currKCa(self, c, Vm): ''' Compute the outward Calcium-dependent Potassium current per unit area. :param c: open-probability of Calcium-dependent Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GKCa = self.GKCaMax * c return GKCa * (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) class LeechPressure(LeechMech): ''' Class defining the membrane channel dynamics of a leech pressure sensory neuron. with 7 different current types: - Inward Sodium current - Outward Potassium current - Inward high-voltage-activated Calcium current - Non-specific leakage current - Calcium-dependent, outward Potassium current - Sodium pump current - Calcium pump current Reference: *Baccus, S.A. (1998). Synaptic facilitation by reflected action potentials: enhancement of transmission when nerve impulses reverse direction at axon branch points. Proc. Natl. Acad. Sci. U.S.A. 95, 8345–8350.* ''' # Name of channel mechanism name = 'LeechP' # Cell-specific biophysical parameters Cm0 = 1e-2 # Cell membrane resting capacitance (F/m2) Vm0 = -48.865 # Cell membrane resting potential (mV) C_Na_out = 0.11 # Sodium extracellular concentration (M) C_Ca_out = 1.8e-3 # Calcium extracellular concentration (M) C_Na_in0 = 0.01 # Initial Sodium intracellular concentration (M) C_Ca_in0 = 1e-7 # Initial Calcium intracellular concentration (M) # VNa = 60 # Sodium Nernst potential, from MOD file on ModelDB (mV) # VCa = 125 # Calcium Nernst potential, from MOD file on ModelDB (mV) VK = -68.0 # Potassium Nernst potential (mV) VL = -49.0 # Non-specific leakage Nernst potential (mV) INaPmax = 70.0 # Maximum pump rate of the NaK-ATPase (mA/m2) khalf_Na = 0.012 # Sodium concentration at which NaK-ATPase is at half its maximum rate (M) ksteep_Na = 1e-3 # Sensitivity of NaK-ATPase to varying Sodium concentrations (M) iCaS = 0.1 # Calcium pump current parameter (mA/m2) GNaMax = 3500.0 # Max. conductance of Sodium current (S/m^2) GKMax = 60.0 # Max. conductance of Potassium current (S/m^2) GCaMax = 0.02 # Max. conductance of Calcium current (S/m^2) GKCaMax = 8.0 # Max. conductance of Calcium-dependent Potassium current (S/m^2) GL = 5.0 # Conductance of non-specific leakage current (S/m^2) diam = 50e-6 # Cell soma diameter (m) Z_Na = 1 # Sodium valence Z_Ca = 2 # Calcium valence # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h', 'm4h'], 'i_K\ kin.': ['n'], 'i_{Ca}\ kin.': ['s'], 'i_{KCa}\ kin.': ['c'], 'pools': ['C_Na', 'C_Ca'], 'I': ['iNa2', 'iK', 'iCa2', 'iKCa2', 'iPumpNa2', 'iPumpCa2', 'iL', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' SV_ratio = 6 / self.diam # surface to volume ratio of the (spherical) cell soma # Conversion constant from membrane ionic currents into # change rate of intracellular ionic concentrations self.K_Na = SV_ratio / (self.Z_Na * self.Faraday) * 1e-6 # Sodium (M/s) self.K_Ca = SV_ratio / (self.Z_Ca * self.Faraday) * 1e-6 # Calcium (M/s) # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'c', 'C_Na', 'C_Ca'] self.states0 = np.array([]) # Names of the channels effective coefficients self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas'] # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) + # Charge interval bounds for lookup creation + self.Qbounds = (np.round(self.Vm0 - 10.0) * 1e-5, 60.0e-5) + def nernst(self, z_ion, C_ion_in, C_ion_out): ''' Return the Nernst potential of a specific ion given its intra and extracellular concentrations. :param z_ion: ion valence :param C_ion_in: intracellular ion concentration (M) :param C_ion_out: extracellular ion concentration (M) :return: ion Nernst potential (mV) ''' return (self.Rg * self.T) / (z_ion * self.Faraday) * np.log(C_ion_out / C_ion_in) * 1e3 def currPumpNa(self, C_Na_in): ''' Outward current mimicking the activity of the NaK-ATPase pump. :param C_Na_in: intracellular Sodium concentration (M) :return: current per unit area (mA/m2) ''' return self.INaPmax / (1 + np.exp((self.khalf_Na - C_Na_in) / self.ksteep_Na)) def currPumpCa(self, C_Ca_in): ''' Outward current representing the activity of a Calcium pump. :param C_Ca_in: intracellular Calcium concentration (M) :return: current per unit area (mA/m2) ''' return self.iCaS * (C_Ca_in - self.C_Ca_in0) / 1.5 def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, c, C_Na_in, C_Ca_in = states return (self.currNa(m, h, Vm, C_Na_in) + self.currK(n, Vm) + self.currCa(s, Vm, C_Ca_in) + self.currKCa(c, Vm) + self.currL(Vm) + (self.currPumpNa(C_Na_in) / 3.) + self.currPumpCa(C_Ca_in)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Intracellular concentrations C_Na_eq = self.C_Na_in0 C_Ca_eq = self.C_Ca_in0 # Standard gating dynamics: Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) seq = self.alphas(Vm) / (self.alphas(Vm) + self.betas(Vm)) ceq = self.alphaC(C_Ca_eq) / (self.alphaC(C_Ca_eq) + self.betaC) return np.array([meq, heq, neq, seq, ceq, C_Na_eq, C_Ca_eq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack states m, h, n, s, c, C_Na_in, C_Ca_in = states # Standard gating states derivatives dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) dcdt = self.derC(c, C_Ca_in) # Intracellular concentrations dCNa_dt = - (self.currNa(m, h, Vm, C_Na_in) + self.currPumpNa(C_Na_in)) * self.K_Na # M/s dCCa_dt = -(self.currCa(s, Vm, C_Ca_in) + self.currPumpCa(C_Ca_in)) * self.K_Ca # M/s # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dCNa_dt, dCCa_dt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) as_avg = np.mean(self.alphas(Vm)) bs_avg = np.mean(self.betas(Vm)) # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) Vmeff = np.interp(Qm, interp_data['Q'], interp_data['V']) # Unpack states m, h, n, s, c, C_Na_in, C_Ca_in = states # Standard gating states derivatives dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s # KCa current gating state derivative dcdt = self.derC(c, C_Ca_in) # Intracellular concentrations dCNa_dt = - (self.currNa(m, h, Vmeff, C_Na_in) + self.currPumpNa(C_Na_in)) * self.K_Na # M/s dCCa_dt = -(self.currCa(s, Vmeff, C_Ca_in) + self.currPumpCa(C_Ca_in)) * self.K_Ca # M/s # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dCNa_dt, dCCa_dt] class LeechRetzius(LeechMech): ''' Class defining the membrane channel dynamics of a leech Retzius neuron. with 5 different current types: - Inward Sodium current - Outward Potassium current - Inward high-voltage-activated Calcium current - Non-specific leakage current - Calcium-dependent, outward Potassium current References: *Vazquez, Y., Mendez, B., Trueta, C., and De-Miguel, F.F. (2009). Summation of excitatory postsynaptic potentials in electrically-coupled neurones. Neuroscience 163, 202–212.* *ModelDB link: https://senselab.med.yale.edu/modeldb/ShowModel.cshtml?model=120910* ''' # Name of channel mechanism # name = 'LeechR' # Cell-specific biophysical parameters Cm0 = 5e-2 # Cell membrane resting capacitance (F/m2) Vm0 = -44.45 # Cell membrane resting potential (mV) VNa = 50.0 # Sodium Nernst potential, from retztemp.ses file on ModelDB (mV) VCa = 125.0 # Calcium Nernst potential, from cachdend.mod file on ModelDB (mV) VK = -79.0 # Potassium Nernst potential, from retztemp.ses file on ModelDB (mV) VL = -30.0 # Non-specific leakage Nernst potential, from leakdend.mod file on ModelDB (mV) GNaMax = 1250.0 # Max. conductance of Sodium current (S/m^2) GKMax = 10.0 # Max. conductance of Potassium current (S/m^2) GAMax = 100.0 # Max. conductance of transient Potassium current (S/m^2) GCaMax = 4.0 # Max. conductance of Calcium current (S/m^2) GKCaMax = 130.0 # Max. conductance of Calcium-dependent Potassium current (S/m^2) GL = 1.25 # Conductance of non-specific leakage current (S/m^2) Vhalf = -73.1 # mV C_Ca_in = 5e-8 # Calcium intracellular concentration, from retztemp.ses file (M) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h', 'm4h'], 'i_K\ kin.': ['n'], 'i_A\ kin.': ['a', 'b', 'ab'], 'i_{Ca}\ kin.': ['s'], 'i_{KCa}\ kin.': ['c'], 'I': ['iNa', 'iK', 'iCa', 'iKCa2', 'iA', 'iL', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'c', 'a', 'b'] self.states0 = np.array([]) # Names of the channels effective coefficients self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas', 'alphac', 'betac', 'alphaa', 'betaa' 'alphab', 'betab'] # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) def ainf(self, Vm): ''' Steady-state activation probability of transient Potassium channels. Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' Vth = -55.0 # mV return 0 if Vm <= Vth else min(1, 2 * (Vm - Vth)**3 / ((11 - Vth)**3 + (Vm - Vth)**3)) def taua(self, Vm): ''' Activation time constant of transient Potassium channels. (assuming T = 20°C). Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' x = -1.5 * (Vm - self.Vhalf) * 1e-3 * self.Faraday / (self.Rg * self.T) # [-] alpha = np.exp(x) # ms-1 beta = np.exp(0.7 * x) # ms-1 return max(0.5, beta / (0.3 * (1 + alpha))) * 1e-3 # s def binf(self, Vm): ''' Steady-state inactivation probability of transient Potassium channels. Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' return 1. / (1 + np.exp((self.Vhalf - Vm) / -6.3)) def taub(self, Vm): ''' Inactivation time constant of transient Potassium channels. (assuming T = 20°C). Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' x = 2 * (Vm - self.Vhalf) * 1e-3 * self.Faraday / (self.Rg * self.T) alpha = np.exp(x) beta = np.exp(0.65 * x) return max(7.5, beta / (0.02 * (1 + alpha))) * 1e-3 # s def derA(self, Vm, a): ''' Compute the evolution of the activation-probability of transient Potassium channels. :param Vm: membrane potential (mV) :param a: activation-probability of transient Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.ainf(Vm) - a) / self.taua(Vm) def derB(self, Vm, b): ''' Compute the evolution of the inactivation-probability of transient Potassium channels. :param Vm: membrane potential (mV) :param b: inactivation-probability of transient Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.binf(Vm) - b) / self.taub(Vm) def currA(self, a, b, Vm): ''' Compute the outward, transient Potassium current per unit area. :param a: open-probability of transient Potassium channels :param b: inactivation-probability of transient Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GAMax * a * b return GK * (Vm - self.VK) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, c, a, b = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currCa(s, Vm) + self.currL(Vm) + self.currKCa(c, Vm) + self.currA(a, b, Vm)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Standard gating dynamics: Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) seq = self.alphas(Vm) / (self.alphas(Vm) + self.betas(Vm)) ceq = self.alphaC(self.C_Ca_in) / (self.alphaC(self.C_Ca_in) + self.betaC) aeq = self.ainf(Vm) beq = self.binf(Vm) return np.array([meq, heq, neq, seq, ceq, aeq, beq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack states m, h, n, s, c, a, b = states # Standard gating states derivatives dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) dcdt = self.derC(c, self.C_Ca_in) dadt = self.derA(Vm, a) dbdt = self.derB(Vm, b) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dadt, dbdt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) as_avg = np.mean(self.alphas(Vm)) bs_avg = np.mean(self.betas(Vm)) Ta = self.taua(Vm) ainf = self.ainf(Vm) aa_avg = np.mean(ainf / Ta) ba_avg = np.mean(1 / Ta) - aa_avg Tb = self.taub(Vm) binf = self.binf(Vm) ab_avg = np.mean(binf / Tb) bb_avg = np.mean(1 / Tb) - ab_avg # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg, aa_avg, ba_avg, ab_avg, bb_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) # Unpack states m, h, n, s, c, a, b = states # Standard gating states derivatives dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s dadt = rates[8] * (1 - a) - rates[9] * a dbdt = rates[10] * (1 - b) - rates[11] * b # KCa current gating state derivative dcdt = self.derC(c, self.C_Ca_in) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dadt, dbdt] diff --git a/PointNICE/neurons/thalamic.py b/PointNICE/neurons/thalamic.py index c7d8e76..2669563 100644 --- a/PointNICE/neurons/thalamic.py +++ b/PointNICE/neurons/thalamic.py @@ -1,794 +1,797 @@ #!/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-30 19:26:17 +# @Last Modified time: 2018-03-13 15:02:14 ''' 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) def __init__(self): ''' Constructor of the class ''' # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'u'] self.states0 = np.array([]) # Names of the different coefficients to be averaged in a lookup table. self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas', 'alphau', 'betau'] + # Charge interval bounds for lookup creation + self.Qbounds = (np.round(self.Vm0 - 10.0) * 1e-5, 50.0e-5) + 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, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) m, h, n, s, u = states dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s dudt = rates[8] * (1 - u) - rates[9] * u return [dmdt, dhdt, dndt, dsdt, dudt] class ThalamicRE(Thalamic): ''' Specific membrane channel dynamics of a thalamic reticular neuron. References: *Destexhe, A., Contreras, D., Steriade, M., Sejnowski, T.J., and Huguenard, J.R. (1996). In vivo, in vitro, and computational analysis of dendritic calcium currents in thalamic reticular neurons. J. Neurosci. 16, 169–185.* *Huguenard, J.R., and Prince, D.A. (1992). A novel T-type current underlies prolonged Ca(2+)-dependent burst firing in GABAergic neurons of rat thalamic reticular nucleus. J. Neurosci. 12, 3804–3817.* ''' # Name of channel mechanism name = 'RE' # Cell-specific biophysical parameters Vm0 = -89.5 # Cell membrane resting potential (mV) GNaMax = 2000.0 # Max. conductance of Sodium current (S/m^2) GKMax = 200.0 # Max. conductance of Potassium current (S/m^2) GTMax = 30.0 # Max. conductance of low-threshold Calcium current (S/m^2) 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) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_K\ kin.': ['n'], 'i_{TS}\ kin.': ['s', 'u'], 'I': ['iNa', 'iK', 'iTs', 'iL', 'iNet'] } 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) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_K\ kin.': ['n'], 'i_{T}\ kin.': ['s', 'u'], 'i_{H}\ kin.': ['O', 'OL', 'O + 2OL'], 'I': ['iNa', 'iK', 'iT', 'iH', 'iKL', 'iL', 'iNet'] } 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, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) Vmeff = np.interp(Qm, interp_data['Q'], interp_data['V']) # Unpack states m, h, n, s, u, O, C, P0, C_Ca = states # 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/solvers/SolverUS.py b/PointNICE/solvers/SolverUS.py index 38d7ec3..62167d5 100644 --- a/PointNICE/solvers/SolverUS.py +++ b/PointNICE/solvers/SolverUS.py @@ -1,786 +1,790 @@ #!/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: 2018-03-13 14:29:47 +# @Last Modified time: 2018-03-13 15:08:52 import os import warnings import pickle import logging import progressbar as pb import numpy as np import scipy.integrate as integrate from scipy.interpolate import interp2d from ..bls import BilayerSonophore from ..utils import * from ..constants import * from ..neurons import BaseMech # Get package logger logger = logging.getLogger('PointNICE') class SolverUS(BilayerSonophore): """ This class extends the BilayerSonophore class by adding a biophysical Hodgkin-Huxley model on top of the mechanical BLS model. """ def __init__(self, diameter, neuron, Fdrive, embedding_depth=0.0): """ Constructor of the class. :param diameter: in-plane diameter of the sonophore structure within the membrane (m) :param neuron: neuron object :param Fdrive: frequency of acoustic perturbation (Hz) :param embedding_depth: depth of the embedding tissue around the membrane (m) """ # Check validity of input parameters assert isinstance(neuron, BaseMech), ('neuron mechanism must be inherited ' 'from the BaseMech class') assert Fdrive >= 0., 'Driving frequency must be positive' # TODO: check parameters dictionary (float type, mandatory members) # Initialize BLS object Cm0 = neuron.Cm0 Vm0 = neuron.Vm0 BilayerSonophore.__init__(self, diameter, Fdrive, Cm0, Cm0 * Vm0 * 1e-3, embedding_depth) logger.debug('US solver initialization with %s neuron', neuron.name) def eqHH(self, t, y, neuron, Cm): """ Compute the derivatives of the n-ODE HH system variables, based on a value of membrane capacitance. :param t: specific instant in time (s) :param y: vector of HH system variables at time t :param neuron: neuron object :param Cm: membrane capacitance (F/m2) :return: vector of HH system derivatives at time t """ # Split input vector explicitly Qm, *states = y # Compute membrane potential Vm = Qm / Cm * 1e3 # mV # Compute derivatives dQm = - neuron.currNet(Vm, states) * 1e-3 # A/m2 dstates = neuron.derStates(Vm, states) # Return derivatives vector return [dQm, *dstates] def eqFull(self, t, y, neuron, Adrive, Fdrive, phi): """ Compute the derivatives of the (n+3) ODE full NBLS system variables. :param t: specific instant in time (s) :param y: vector of state variables :param neuron: neuron object :param Adrive: acoustic drive amplitude (Pa) :param Fdrive: acoustic drive frequency (Hz) :param phi: acoustic drive phase (rad) :return: vector of derivatives """ # Compute derivatives of mechanical and electrical systems dydt_mech = self.eqMech(t, y[:3], Adrive, Fdrive, y[3], phi) dydt_elec = self.eqHH(t, y[3:], neuron, self.Capct(y[1])) # return concatenated output return dydt_mech + dydt_elec def eqHHeff(self, t, y, neuron, interp_data): """ Compute the derivatives of the n-ODE effective HH system variables, based on 1-dimensional linear interpolation of "effective" coefficients that summarize the system's behaviour over an acoustic cycle. :param t: specific instant in time (s) :param y: vector of HH system variables at time t :param neuron: neuron object :param interp_data: dictionary of 1D data points of "effective" coefficients over the charge domain, for specific frequency and amplitude values. :return: vector of effective system derivatives at time t """ # Split input vector explicitly Qm, *states = y # Compute charge and channel states variation Vm = np.interp(Qm, interp_data['Q'], interp_data['V']) # mV dQmdt = - neuron.currNet(Vm, states) * 1e-3 dstates = neuron.derStatesEff(Qm, states, interp_data) # Return derivatives vector return [dQmdt, *dstates] def getEffCoeffs(self, neuron, Fdrive, Adrive, Qm, phi=np.pi): """ Compute "effective" coefficients of the HH system for a specific combination of stimulus frequency, stimulus amplitude and charge density. A short mechanical simulation is run while imposing the specific charge density, until periodic stabilization. The HH coefficients are then averaged over the last acoustic cycle to yield "effective" coefficients. :param neuron: neuron object :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param Qm: imposed charge density (C/m2) :param phi: acoustic drive phase (rad) :return: tuple with the effective potential, gas content and channel rates """ # Run simulation and retrieve deflection and gas content vectors from last cycle (_, y, _) = self.runMech(Fdrive, Adrive, Qm, phi) (Z, ng) = y Z_last = Z[-NPC_FULL:] # m # Compute membrane potential vector Vm = np.array([Qm / self.Capct(ZZ) * 1e3 for ZZ in Z_last]) # mV # Compute average cycle value for membrane potential and rate constants Vm_eff = np.mean(Vm) # mV rates_eff = neuron.getEffRates(Vm) # Take final cycle value for gas content ng_eff = ng[-1] # mole return (Vm_eff, ng_eff, *rates_eff) - def createLookup(self, neuron, freqs, amps, charges, phi=np.pi): + def createLookup(self, neuron, freqs, amps, phi=np.pi): """ Run simulations of the mechanical system for a multiple combinations of imposed charge densities and acoustic amplitudes, compute effective coefficients and store them as 2D arrays in a lookup file. :param neuron: neuron object :param freqs: array of acoustic drive frequencies (Hz) :param amps: array of acoustic drive amplitudes (Pa) - :param charges: array of charge densities (C/m2) :param phi: acoustic drive phase (rad) """ + # Check if lookup file already exists + lookup_file = '{}_lookups_a{:.1f}nm.pkl'.format(neuron.name, self.a * 1e9) + lookup_filepath = '{0}/{1}'.format(getLookupDir(), lookup_file) + assert not os.path.isfile(lookup_filepath), '"{}" file already exists'.format(lookup_file) + # Check validity of stimulation parameters assert freqs.min() > 0, 'Driving frequencies must be strictly positive' assert amps.min() >= 0, 'Acoustic pressure amplitudes must be positive' logger.info('Creating lookup table for %s neuron', neuron.name) + # Create neuron-specific charge vector + charges = np.arange(neuron.Qbounds[0], neuron.Qbounds[1] + 1e-5, 1e-5) # C/m2 + # Initialize lookup dictionary of 3D array to store effective coefficients nf = freqs.size nA = amps.size nQ = charges.size coeffs_names = ['V', 'ng', *neuron.coeff_names] ncoeffs = len(coeffs_names) lookup_dict = {cn: np.empty((nf, nA, nQ)) for cn in coeffs_names} # Loop through all (f, A, Q) combinations nsims = nf * nA * nQ isim = 0 log_str = 'short simulation %u/%u (f = %.2f kHz, A = %.2f kPa, Q = %.2f nC/cm2)' for i in range(nf): for j in range(nA): for k in range(nQ): isim += 1 # Run short simulation and store effective coefficients logger.info(log_str, isim, nsims, freqs[i] * 1e-3, amps[j] * 1e-3, charges[k] * 1e5) sim_coeffs = self.getEffCoeffs(neuron, freqs[i], amps[j], charges[k], phi) for icoeff in range(ncoeffs): lookup_dict[coeffs_names[icoeff]][i, j, k] = sim_coeffs[icoeff] # Add input frequency, amplitude and charge arrays to lookup dictionary lookup_dict['f'] = freqs # Hz lookup_dict['A'] = amps # Pa lookup_dict['Q'] = charges # C/m2 # Save dictionary in lookup file - lookup_file = '{}_lookups_a{:.1f}nm.pkl'.format(neuron.name, self.a * 1e9) logger.info('Saving %s neuron lookup table in file: "%s"', neuron.name, lookup_file) - - lookup_filepath = '{0}/{1}'.format(getLookupDir(), lookup_file) with open(lookup_filepath, 'wb') as fh: pickle.dump(lookup_dict, fh) def __runClassic(self, neuron, Fdrive, Adrive, tstim, toffset, PRF, DC, phi=np.pi): """ Compute solutions of the system for a specific set of US stimulation parameters, using a classic integration scheme. The first iteration uses the quasi-steady simplification to compute the initiation of motion from a flat leaflet configuration. Afterwards, the ODE system is solved iteratively until completion. :param neuron: neuron object :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param tstim: duration of US stimulation (s) :param toffset: duration of the offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) :param phi: acoustic drive phase (rad) :return: 3-tuple with the time profile, the effective solution matrix and a state vector """ # Raise warnings as error warnings.filterwarnings('error') # Initialize system solver solver_full = integrate.ode(self.eqFull) solver_full.set_integrator('lsoda', nsteps=SOLVER_NSTEPS, ixpr=True) # Determine system time step Tdrive = 1 / Fdrive dt = Tdrive / NPC_FULL # if CW stimulus: divide integration during stimulus into 100 intervals if DC == 1.0: PRF = 100 / tstim # Compute vector sizes npulses = int(np.round(PRF * tstim)) Tpulse_on = DC / PRF Tpulse_off = (1 - DC) / PRF n_pulse_on = int(np.round(Tpulse_on / dt)) n_pulse_off = int(np.round(Tpulse_off / dt)) n_off = int(np.round(toffset / dt)) # Solve quasi-steady equation to compute first deflection value Z0 = 0.0 ng0 = self.ng0 Qm0 = self.Qm0 Pac1 = self.Pacoustic(dt, Adrive, Fdrive, phi) Z1 = self.balancedefQS(ng0, Qm0, Pac1) # Initialize global arrays states = np.array([1, 1]) t = np.array([0., dt]) y_membrane = np.array([[0., (Z1 - Z0) / dt], [Z0, Z1], [ng0, ng0], [Qm0, Qm0]]) y_channels = np.tile(neuron.states0, (2, 1)).T y = np.vstack((y_membrane, y_channels)) nvar = y.shape[0] # Initialize pulse time and states vectors t_pulse0 = np.linspace(0, Tpulse_on + Tpulse_off, n_pulse_on + n_pulse_off) states_pulse = np.concatenate((np.ones(n_pulse_on), np.zeros(n_pulse_off))) # Initialize progress bar if logger.getEffectiveLevel() == logging.DEBUG: widgets = ['Running: ', pb.Percentage(), ' ', pb.Bar(), ' ', pb.ETA()] pbar = pb.ProgressBar(widgets=widgets, max_value=int(npulses * (toffset + tstim) / tstim)) pbar.start() # Loop through all pulse (ON and OFF) intervals for i in range(npulses): # Construct and initialize arrays t_pulse = t_pulse0 + t[-1] y_pulse = np.empty((nvar, n_pulse_on + n_pulse_off)) y_pulse[:, 0] = y[:, -1] # Initialize iterator k = 0 # Integrate ON system solver_full.set_f_params(neuron, Adrive, Fdrive, phi) solver_full.set_initial_value(y_pulse[:, k], t_pulse[k]) while solver_full.successful() and k < n_pulse_on - 1: k += 1 solver_full.integrate(t_pulse[k]) y_pulse[:, k] = solver_full.y # Integrate OFF system solver_full.set_f_params(neuron, 0.0, 0.0, 0.0) solver_full.set_initial_value(y_pulse[:, k], t_pulse[k]) while solver_full.successful() and k < n_pulse_on + n_pulse_off - 1: k += 1 solver_full.integrate(t_pulse[k]) y_pulse[:, k] = solver_full.y # Append pulse arrays to global arrays states = np.concatenate([states, states_pulse[1:]]) t = np.concatenate([t, t_pulse[1:]]) y = np.concatenate([y, y_pulse[:, 1:]], axis=1) # Update progress bar if logger.getEffectiveLevel() == logging.DEBUG: pbar.update(i) # Integrate offset interval if n_off > 0: t_off = np.linspace(0, toffset, n_off) + t[-1] states_off = np.zeros(n_off) y_off = np.empty((nvar, n_off)) y_off[:, 0] = y[:, -1] solver_full.set_initial_value(y_off[:, 0], t_off[0]) solver_full.set_f_params(neuron, 0.0, 0.0, 0.0) k = 0 while solver_full.successful() and k < n_off - 1: k += 1 solver_full.integrate(t_off[k]) y_off[:, k] = solver_full.y # Concatenate offset arrays to global arrays states = np.concatenate([states, states_off[1:]]) t = np.concatenate([t, t_off[1:]]) y = np.concatenate([y, y_off[:, 1:]], axis=1) # Terminate progress bar if logger.getEffectiveLevel() == logging.DEBUG: pbar.finish() # Downsample arrays in time-domain to reduce overall size t = t[::CLASSIC_DS_FACTOR] y = y[:, ::CLASSIC_DS_FACTOR] states = states[::CLASSIC_DS_FACTOR] # Return output variables return (t, y[1:, :], states) def __runEffective(self, neuron, Fdrive, Adrive, tstim, toffset, PRF, DC, dt=DT_EFF): """ Compute solutions of the system for a specific set of US stimulation parameters, using charge-predicted "effective" coefficients to solve the HH equations at each step. :param neuron: neuron object :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param tstim: duration of US stimulation (s) :param toffset: duration of the offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) :param dt: integration time step (s) :return: 3-tuple with the time profile, the effective solution matrix and a state vector """ # Raise warnings as error warnings.filterwarnings('error') # Check lookup file existence lookup_file = '{}_lookups_a{:.1f}nm.pkl'.format(neuron.name, self.a * 1e9) lookup_path = '{}/{}'.format(getLookupDir(), lookup_file) assert os.path.isfile(lookup_path), ('No lookup file available for {} ' 'neuron type').format(neuron.name) # Load coefficients with open(lookup_path, 'rb') as fh: lookup_dict = pickle.load(fh) # Retrieve 1D inputs from lookup dictionary freqs = lookup_dict['f'] amps = lookup_dict['A'] charges = lookup_dict['Q'] # Check that stimulation parameters are within lookup range margin = 1e-9 # adding margin to compensate for eventual round error frange = (freqs.min() - margin, freqs.max() + margin) Arange = (amps.min() - margin, amps.max() + margin) assert frange[0] <= Fdrive <= frange[1], \ 'Fdrive must be within [{:.1f}, {:.1f}] kHz'.format(*[f * 1e-3 for f in frange]) assert Arange[0] <= Adrive <= Arange[1], \ 'Adrive must be within [{:.1f}, {:.1f}] kPa'.format(*[A * 1e-3 for A in Arange]) # Define interpolation datasets to be projected coeffs_list = ['V', 'ng', *neuron.coeff_names] # If Fdrive in lookup frequencies, simply project (A, Q) interpolation dataset # at Fdrive index onto 1D charge-based interpolation dataset if Fdrive in freqs: iFdrive = np.searchsorted(freqs, Fdrive) logger.debug('Using lookups directly at %.2f kHz', freqs[iFdrive] * 1e-3) coeffs1d = {} for cn in coeffs_list: coeff2d = np.squeeze(lookup_dict[cn][iFdrive, :, :]) itrp = interp2d(amps, charges, coeff2d.T) coeffs1d[cn] = itrp(Adrive, charges) if cn == 'ng': coeffs1d['ng0'] = itrp(0.0, charges) # Otherwise, project 2 (A, Q) interpolation datasets at Fdrive bounding values # indexes in lookup frequencies onto two 1D charge-based interpolation datasets, and # interpolate between them afterwards else: ilb = np.searchsorted(freqs, Fdrive) - 1 logger.debug('Interpolating lookups between %.2f kHz and %.2f kHz', freqs[ilb] * 1e-3, freqs[ilb + 1] * 1e-3) coeffs1d = {} for cn in coeffs_list: coeffs1d_bounds = [] ng0_bounds = [] for iFdrive in [ilb, ilb + 1]: coeff2d = np.squeeze(lookup_dict[cn][iFdrive, :, :]) itrp = interp2d(amps, charges, coeff2d.T) coeffs1d_bounds.append(itrp(Adrive, charges)) if cn == 'ng': ng0_bounds.append(itrp(0.0, charges)) coeffs1d_bounds = np.squeeze(np.array([coeffs1d_bounds])) itrp = interp2d(freqs[ilb:ilb + 2], charges, coeffs1d_bounds.T) coeffs1d[cn] = itrp(Fdrive, charges) if cn == 'ng': ng0_bounds = np.squeeze(np.array([ng0_bounds])) itrp = interp2d(freqs[ilb:ilb + 2], charges, ng0_bounds.T) coeffs1d['ng0'] = itrp(Fdrive, charges) # Squeeze interpolated vectors extra dimensions and add input charges vector coeffs1d = {key: np.squeeze(value) for key, value in coeffs1d.items()} coeffs1d['Q'] = charges # Initialize system solvers solver_on = integrate.ode(self.eqHHeff) solver_on.set_integrator('lsoda', nsteps=SOLVER_NSTEPS) solver_on.set_f_params(neuron, coeffs1d) solver_off = integrate.ode(self.eqHH) solver_off.set_integrator('lsoda', nsteps=SOLVER_NSTEPS) # if CW stimulus: change PRF to have exactly one integration interval during stimulus if DC == 1.0: PRF = 1 / tstim # Compute vector sizes npulses = int(np.round(PRF * tstim)) Tpulse_on = DC / PRF Tpulse_off = (1 - DC) / PRF n_pulse_on = int(np.round(Tpulse_on / dt)) + 1 n_pulse_off = int(np.round(Tpulse_off / dt)) # For high-PRF pulsed protocols: adapt time step if greater than TON or TOFF dt_warning_msg = 'high-PRF protocol: lowering integration time step to %.2e ms to match %s' if Tpulse_on > 0 and n_pulse_on == 0: logger.warning(dt_warning_msg, Tpulse_on * 1e3, 'TON') dt = Tpulse_on n_pulse_on = int(np.round(Tpulse_on / dt)) n_pulse_off = int(np.round(Tpulse_off / dt)) if Tpulse_off > 0 and n_pulse_off == 0: logger.warning(dt_warning_msg, Tpulse_off * 1e3, 'TOFF') dt = Tpulse_off n_pulse_on = int(np.round(Tpulse_on / dt)) n_pulse_off = int(np.round(Tpulse_off / dt)) # Compute ofset size n_off = int(np.round(toffset / dt)) # Initialize global arrays states = np.array([1]) t = np.array([0.0]) y = np.atleast_2d(np.insert(neuron.states0, 0, self.Qm0)).T nvar = y.shape[0] Zeff = np.array([0.0]) ngeff = np.array([self.ng0]) # Initializing accurate pulse time vector t_pulse_on = np.linspace(0, Tpulse_on, n_pulse_on) t_pulse_off = np.linspace(dt, Tpulse_off, n_pulse_off) + Tpulse_on t_pulse0 = np.concatenate([t_pulse_on, t_pulse_off]) states_pulse = np.concatenate((np.ones(n_pulse_on), np.zeros(n_pulse_off))) # Loop through all pulse (ON and OFF) intervals for i in range(npulses): # Construct and initialize arrays t_pulse = t_pulse0 + t[-1] y_pulse = np.empty((nvar, n_pulse_on + n_pulse_off)) ngeff_pulse = np.empty(n_pulse_on + n_pulse_off) Zeff_pulse = np.empty(n_pulse_on + n_pulse_off) y_pulse[:, 0] = y[:, -1] ngeff_pulse[0] = ngeff[-1] Zeff_pulse[0] = Zeff[-1] # Initialize iterator k = 0 # Integrate ON system solver_on.set_initial_value(y_pulse[:, k], t_pulse[k]) while solver_on.successful() and k < n_pulse_on - 1: k += 1 solver_on.integrate(t_pulse[k]) y_pulse[:, k] = solver_on.y ngeff_pulse[k] = np.interp(y_pulse[0, k], coeffs1d['Q'], coeffs1d['ng']) # mole Zeff_pulse[k] = self.balancedefQS(ngeff_pulse[k], y_pulse[0, k]) # m # Integrate OFF system solver_off.set_initial_value(y_pulse[:, k], t_pulse[k]) solver_off.set_f_params(neuron, self.Capct(Zeff_pulse[k])) while solver_off.successful() and k < n_pulse_on + n_pulse_off - 1: k += 1 solver_off.integrate(t_pulse[k]) y_pulse[:, k] = solver_off.y ngeff_pulse[k] = np.interp(y_pulse[0, k], coeffs1d['Q'], coeffs1d['ng0']) # mole Zeff_pulse[k] = self.balancedefQS(ngeff_pulse[k], y_pulse[0, k]) # m solver_off.set_f_params(neuron, self.Capct(Zeff_pulse[k])) # Append pulse arrays to global arrays states = np.concatenate([states[:-1], states_pulse]) t = np.concatenate([t, t_pulse[1:]]) y = np.concatenate([y, y_pulse[:, 1:]], axis=1) Zeff = np.concatenate([Zeff, Zeff_pulse[1:]]) ngeff = np.concatenate([ngeff, ngeff_pulse[1:]]) # Integrate offset interval if n_off > 0: t_off = np.linspace(0, toffset, n_off) + t[-1] states_off = np.zeros(n_off) y_off = np.empty((nvar, n_off)) ngeff_off = np.empty(n_off) Zeff_off = np.empty(n_off) y_off[:, 0] = y[:, -1] ngeff_off[0] = ngeff[-1] Zeff_off[0] = Zeff[-1] solver_off.set_initial_value(y_off[:, 0], t_off[0]) solver_off.set_f_params(neuron, self.Capct(Zeff_pulse[k])) k = 0 while solver_off.successful() and k < n_off - 1: k += 1 solver_off.integrate(t_off[k]) y_off[:, k] = solver_off.y ngeff_off[k] = np.interp(y_off[0, k], coeffs1d['Q'], coeffs1d['ng0']) # mole Zeff_off[k] = self.balancedefQS(ngeff_off[k], y_off[0, k]) # m solver_off.set_f_params(neuron, self.Capct(Zeff_off[k])) # Concatenate offset arrays to global arrays states = np.concatenate([states, states_off[1:]]) t = np.concatenate([t, t_off[1:]]) y = np.concatenate([y, y_off[:, 1:]], axis=1) Zeff = np.concatenate([Zeff, Zeff_off[1:]]) ngeff = np.concatenate([ngeff, ngeff_off[1:]]) # Add Zeff and ngeff to solution matrix y = np.vstack([Zeff, ngeff, y]) # return output variables return (t, y, states) def __runHybrid(self, neuron, Fdrive, Adrive, tstim, toffset, phi=np.pi): """ Compute solutions of the system for a specific set of US stimulation parameters, using a hybrid integration scheme. The first iteration uses the quasi-steady simplification to compute the initiation of motion from a flat leaflet configuration. Afterwards, the NBLS ODE system is solved iteratively for "slices" of N-microseconds, in a 2-steps scheme: - First, the full (n+3) ODE system is integrated for a few acoustic cycles until Z and ng reach a stable periodic solution (limit cycle) - Second, the signals of the 3 mechanical variables over the last acoustic period are selected and resampled to a far lower sampling rate - Third, the HH n-ODE system is integrated for the remaining time of the slice, using periodic expansion of the mechanical signals to precompute the values of capacitance. :param neuron: neuron object :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param tstim: duration of US stimulation (s) :param toffset: duration of the offset (s) :param phi: acoustic drive phase (rad) :return: 3-tuple with the time profile, the solution matrix and a state vector .. warning:: This method cannot handle pulsed stimuli """ # Raise warnings as error warnings.filterwarnings('error') # Initialize full and HH systems solvers solver_full = integrate.ode(self.eqFull) solver_full.set_f_params(neuron, Adrive, Fdrive, phi) solver_full.set_integrator('lsoda', nsteps=SOLVER_NSTEPS) solver_hh = integrate.ode(self.eqHH) solver_hh.set_integrator('dop853', nsteps=SOLVER_NSTEPS, atol=1e-12) # Determine full and HH systems time steps Tdrive = 1 / Fdrive dt_full = Tdrive / NPC_FULL dt_hh = Tdrive / NPC_HH n_full_per_hh = int(NPC_FULL / NPC_HH) t_full_cycle = np.linspace(0, Tdrive - dt_full, NPC_FULL) t_hh_cycle = np.linspace(0, Tdrive - dt_hh, NPC_HH) # Determine number of samples in prediction vectors npc_pred = NPC_FULL - n_full_per_hh + 1 # Solve quasi-steady equation to compute first deflection value Z0 = 0.0 ng0 = self.ng0 Qm0 = self.Qm0 Pac1 = self.Pacoustic(dt_full, Adrive, Fdrive, phi) Z1 = self.balancedefQS(ng0, Qm0, Pac1) # Initialize global arrays states = np.array([1, 1]) t = np.array([0., dt_full]) y_membrane = np.array([[0., (Z1 - Z0) / dt_full], [Z0, Z1], [ng0, ng0], [Qm0, Qm0]]) y_channels = np.tile(neuron.states0, (2, 1)).T y = np.vstack((y_membrane, y_channels)) nvar = y.shape[0] # Initialize progress bar if logger.getEffectiveLevel() == logging.DEBUG: widgets = ['Running: ', pb.Percentage(), ' ', pb.Bar(), ' ', pb.ETA()] pbar = pb.ProgressBar(widgets=widgets, max_value=1000) pbar.start() # For each hybrid integration interval irep = 0 sim_error = False while not sim_error and t[-1] < tstim + toffset: # Integrate full system for a few acoustic cycles until stabilization periodic_conv = False j = 0 ng_last = None Z_last = None while not sim_error and not periodic_conv: if t[-1] > tstim: solver_full.set_f_params(neuron, 0.0, 0.0, 0.0) t_full = t_full_cycle + t[-1] + dt_full y_full = np.empty((nvar, NPC_FULL)) y0_full = y[:, -1] solver_full.set_initial_value(y0_full, t[-1]) k = 0 try: # try to integrate and catch errors/warnings while solver_full.successful() and k <= NPC_FULL - 1: solver_full.integrate(t_full[k]) y_full[:, k] = solver_full.y k += 1 except (Warning, AssertionError) as inst: sim_error = True logger.error('Full system integration error at step %u', k) logger.error(inst) # Compare Z and ng signals over the last 2 acoustic periods if j > 0 and rmse(Z_last, y_full[1, :]) < Z_ERR_MAX \ and rmse(ng_last, y_full[2, :]) < NG_ERR_MAX: periodic_conv = True # Update last vectors for next comparison Z_last = y_full[1, :] ng_last = y_full[2, :] # Concatenate time and solutions to global vectors states = np.concatenate([states, np.ones(NPC_FULL)], axis=0) t = np.concatenate([t, t_full], axis=0) y = np.concatenate([y, y_full], axis=1) # Increment loop index j += 1 # Retrieve last period of the 3 mechanical variables to propagate in HH system t_last = t[-npc_pred:] mech_last = y[0:3, -npc_pred:] # Downsample signals to specified HH system time step (_, mech_pred) = DownSample(t_last, mech_last, NPC_HH) # Integrate HH system until certain dQ or dT is reached Q0 = y[3, -1] dQ = 0.0 t0_interval = t[-1] dt_interval = 0.0 j = 0 if t[-1] < tstim: tlim = tstim else: tlim = tstim + toffset while (not sim_error and t[-1] < tlim and (np.abs(dQ) < DQ_UPDATE or dt_interval < DT_UPDATE)): t_hh = t_hh_cycle + t[-1] + dt_hh y_hh = np.empty((nvar - 3, NPC_HH)) y0_hh = y[3:, -1] solver_hh.set_initial_value(y0_hh, t[-1]) k = 0 try: # try to integrate and catch errors/warnings while solver_hh.successful() and k <= NPC_HH - 1: solver_hh.set_f_params(neuron, self.Capct(mech_pred[1, k])) solver_hh.integrate(t_hh[k]) y_hh[:, k] = solver_hh.y k += 1 except (Warning, AssertionError) as inst: sim_error = True logger.error('HH system integration error at step %u', k) logger.error(inst) # Concatenate time and solutions to global vectors states = np.concatenate([states, np.zeros(NPC_HH)], axis=0) t = np.concatenate([t, t_hh], axis=0) y = np.concatenate([y, np.concatenate([mech_pred, y_hh], axis=0)], axis=1) # Compute charge variation from interval beginning dQ = y[3, -1] - Q0 dt_interval = t[-1] - t0_interval # Increment loop index j += 1 # Update progress bar if logger.getEffectiveLevel() == logging.DEBUG: pbar.update(int(1000 * (t[-1] / (tstim + toffset)))) irep += 1 # Terminate progress bar if logger.getEffectiveLevel() == logging.DEBUG: pbar.finish() # Return output return (t, y[1:, :], states) def run(self, neuron, Fdrive, Adrive, tstim, toffset, PRF=None, DC=1.0, sim_type='effective'): """ Run simulation of the system for a specific set of US stimulation parameters. :param neuron: neuron object :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param tstim: duration of US stimulation (s) :param toffset: duration of the offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) :param sim_type: selected integration method :return: 3-tuple with the time profile, the solution matrix and a state vector """ # # Check validity of simulation type sim_types = ('classic, effective, hybrid') assert sim_type in sim_types, 'Allowed simulation types are {}'.format(sim_types) # Check validity of stimulation parameters assert isinstance(neuron, BaseMech), ('neuron mechanism must be inherited ' 'from the BaseMech class') for param in [Fdrive, Adrive, tstim, toffset, DC]: assert isinstance(param, float), 'stimulation parameters must be float typed' assert Fdrive > 0, 'Driving frequency must be strictly positive' assert Adrive >= 0, 'Acoustic pressure amplitude must be positive' assert tstim > 0, 'Stimulus duration must be strictly positive' assert toffset >= 0, 'Stimulus offset must be positive or null' assert DC > 0 and DC <= 1, 'Duty cycle must be within [0; 1)' if DC < 1.0: assert isinstance(PRF, float), 'if provided, the PRF parameter must be float typed' assert PRF is not None, 'PRF must be provided when using duty cycles smaller than 1' assert PRF >= 1 / tstim, 'PR interval must be smaller than stimulus duration' assert PRF < Fdrive, 'PRF must be smaller than driving frequency' # Call appropriate simulation function if sim_type == 'classic': return self.__runClassic(neuron, Fdrive, Adrive, tstim, toffset, PRF, DC) elif sim_type == 'effective': return self.__runEffective(neuron, Fdrive, Adrive, tstim, toffset, PRF, DC) elif sim_type == 'hybrid': assert DC == 1.0, 'Hybrid method can only handle continuous wave stimuli' return self.__runHybrid(neuron, Fdrive, Adrive, tstim, toffset) diff --git a/plot/plot_P_vs_I.py b/plot/plot_P_vs_I.py index 2ff6837..025ef83 100644 --- a/plot/plot_P_vs_I.py +++ b/plot/plot_P_vs_I.py @@ -1,35 +1,34 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-07-17 11:47:50 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-07-21 18:40:04 +# @Last Modified time: 2018-03-13 15:21:18 ''' plot profile of acoustic Intensity (in W/cm^2) vs Pressure (in kPa) ''' import numpy as np import matplotlib.pyplot as plt from PointNICE.utils import Pressure2Intensity rho = 1075 # kg/m3 c = 1515 # m/s +P = np.logspace(np.log10(1e1), np.log10(1e7), num=500) # Pa +Int = Pressure2Intensity(P, rho, c) # W/m2 fig, ax = plt.subplots() ax.set_xlabel('$Pressure\ (kPa)$') ax.set_ylabel('$I_{SPPA}\ (W/cm^2)$') ax.set_xscale('log') ax.set_yscale('log') - -P = np.logspace(np.log10(1e1), np.log10(1e7), num=500) # Pa -Int = Pressure2Intensity(P, rho, c) # W/m2 ax.plot(P * 1e-3, Int * 1e-4) Psnaps = np.logspace(1, 7, 7) # Pa for Psnap in Psnaps: Isnap = Pressure2Intensity(Psnap, rho, c) # W/m2 ax.plot(np.array([Psnap, Psnap]) * 1e-3, np.array([0.0, Isnap]) * 1e-4, '--', color='black') ax.plot(np.array([0, Psnap]) * 1e-3, np.array([Isnap, Isnap]) * 1e-4, '--', color='black') plt.show() diff --git a/sim/ASTIM_lookups.py b/sim/ASTIM_lookups.py index a7264a5..5d970f5 100644 --- a/sim/ASTIM_lookups.py +++ b/sim/ASTIM_lookups.py @@ -1,42 +1,43 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-06-02 17:50:10 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-03-13 14:32:29 +# @Last Modified time: 2018-03-13 17:13:16 """ Create lookup tables for different acoustic frequencies. """ import logging import numpy as np import PointNICE from PointNICE.utils import logger from PointNICE.neurons import * # Set logging level logger.setLevel(logging.DEBUG) # BLS diameter (m) a = 32e-9 # Channel mechanisms -neurons = [LeechPressure()] +neurons = [CorticalRS()] # Stimulation parameters -freqs = np.arange(100, 1001, 100) * 1e3 # Hz +freqs = np.array([20., 100., 500., 1000., 2000., 3000., 4000.]) * 1e3 # Hz amps = np.logspace(np.log10(0.1), np.log10(600), num=50) * 1e3 # Pa amps = np.insert(amps, 0, 0.0) # adding amplitude 0 logger.info('Starting batch lookup creation') for neuron in neurons: # Create a SolverUS instance (with dummy frequency parameter) solver = PointNICE.SolverUS(a, neuron, 0.0) - charges = np.arange(np.round(neuron.Vm0 - 10.0), 50.0 + 1.0, 1.0) * 1e-5 # C/m2 # Create lookup file - solver.createLookup(neuron, freqs, amps, charges) - -logger.info('Lookup tables successfully created') + try: + solver.createLookup(neuron, freqs, amps) + logger.info('%s Lookup table successfully created', neuron.name) + except AssertionError as err: + logger.error(err)