diff --git a/PointNICE/channels/base.py b/PointNICE/channels/base.py index 0715e16..d359b20 100644 --- a/PointNICE/channels/base.py +++ b/PointNICE/channels/base.py @@ -1,126 +1,126 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-03 11:53:04 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-08-24 17:16:19 +# @Last Modified time: 2017-08-28 17:50:40 ''' Module standard API for all neuron mechanisms. Each mechanism class can use different methods to define the membrane dynamics of a specific neuron type. However, they must contain some mandatory attributes and methods in order to be properly imported in other PointNICE modules and used in NICE simulations. ''' import abc class BaseMech(metaclass=abc.ABCMeta): ''' Abstract class defining the common API (i.e. mandatory attributes and methods) of all subclasses implementing the channels mechanisms of specific neurons. The mandatory attributes are: - **name**: a string defining the name of the mechanism. - **Cm0**: a float defining the membrane resting capacitance (in F/m2) - **Vm0**: a float defining the membrane resting potential (in mV) - **states_names**: a list of strings defining the names of the different state probabilities governing the channels behaviour (i.e. the differential HH variables). - **states0**: a 1D array of floats (NOT integers !!!) defining the initial values of the different state probabilities. - **coeff_names**: a list of strings defining the names of the different coefficients to be used in effective simulations. The mandatory methods are: - **currNet**: compute the net ionic current density (in mA/m2) across the membrane, given a specific membrane potential (in mV) and channel states. - **steadyStates**: compute the channels steady-state values for a specific membrane potential value (in mV). - **derStates**: compute the derivatives of channel states, given a specific membrane potential (in mV) and channel states. This method must return a list of derivatives ordered identically as in the states0 attribute. - **getEffRates**: get the effective rate constants of ion channels to be used in effective simulations. This method must return an array of effective rates ordered identically as in the coeff_names attribute. - **derStatesEff**: compute the effective derivatives of channel states, based on - 2-dimensional linear interpolators of "effective" coefficients. This method must + 1-dimensional linear interpolators of "effective" coefficients. This method must return a list of derivatives ordered identically as in the states0 attribute. ''' @property @abc.abstractmethod def name(self): return 'Should never reach here' @property @abc.abstractmethod def Cm0(self): return 'Should never reach here' @property @abc.abstractmethod def Vm0(self): return 'Should never reach here' # @property # @abc.abstractmethod # def states_names(self): # return 'Should never reach here' # @property # @abc.abstractmethod # def states0(self): # return 'Should never reach here' # @property # @abc.abstractmethod # def coeff_names(self): # return 'Should never reach here' @abc.abstractmethod def currNet(self, Vm, states): ''' Compute the net ionic current per unit area. :param Vm: membrane potential (mV) :states: state probabilities of the ion channels :return: current per unit area (mA/m2) ''' @abc.abstractmethod def steadyStates(self, Vm): ''' Compute the channels steady-state values for a specific membrane potential value. :param Vm: membrane potential (mV) :return: array of steady-states ''' @abc.abstractmethod def derStates(self, Vm, states): ''' Compute the derivatives of channel states. :param Vm: membrane potential (mV) :states: state probabilities of the ion channels :return: current per unit area (mA/m2) ''' @abc.abstractmethod def getEffRates(self, Vm): ''' Get the effective rate constants of ion channels, averaged along an acoustic cycle, for future use in effective simulations. :param Vm: array of membrane potential values for an acoustic cycle (mV) :return: an array of rate average constants (s-1) ''' @abc.abstractmethod - def derStatesEff(self, Adrive, Qm, states, interpolators): + def derStatesEff(self, Qm, states, interp_data): ''' Compute the effective derivatives of channel states, based on - 2-dimensional linear interpolation of "effective" coefficients + 1-dimensional linear interpolation of "effective" coefficients that summarize the system's behaviour over an acoustic cycle. - :param Vm_eff: effective membrane potential (mV) + :param Qm: membrane charge density (C/m2) :states: state probabilities of the ion channels - :param interpolators: dictionary of 2-dimensional linear interpolators - of "effective" rates over the 2D amplitude x charge input domain. + :param interp_data: dictionary of 1D vectors of "effective" coefficients + over the charge domain, for specific frequency and amplitude values. ''' diff --git a/PointNICE/channels/cortical.py b/PointNICE/channels/cortical.py index b8ee8ad..912b015 100644 --- a/PointNICE/channels/cortical.py +++ b/PointNICE/channels/cortical.py @@ -1,603 +1,605 @@ #!/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-25 18:03:08 +# @Last Modified time: 2017-08-28 17:54:21 ''' 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'] def alpham(self, Vm): ''' Compute the alpha rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = (-0.32 * (Vdiff - 13) / (np.exp(- (Vdiff - 13) / 4) - 1)) # ms-1 return alpha * 1e3 # s-1 def betam(self, Vm): ''' Compute the beta rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = (0.28 * (Vdiff - 40) / (np.exp((Vdiff - 40) / 5) - 1)) # ms-1 return beta * 1e3 # s-1 def alphah(self, Vm): ''' Compute the alpha rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = (0.128 * np.exp(-(Vdiff - 17) / 18)) # ms-1 return alpha * 1e3 # s-1 def betah(self, Vm): ''' Compute the beta rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = (4 / (1 + np.exp(-(Vdiff - 40) / 5))) # ms-1 return beta * 1e3 # s-1 def alphan(self, Vm): ''' Compute the alpha rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = (-0.032 * (Vdiff - 15) / (np.exp(-(Vdiff - 15) / 5) - 1)) # ms-1 return alpha * 1e3 # s-1 def betan(self, Vm): ''' Compute the beta rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = (0.5 * np.exp(-(Vdiff - 10) / 40)) # ms-1 return beta * 1e3 # s-1 def pinf(self, Vm): ''' Compute the asymptotic value of the open-probability of slow non-inactivating Potassium channels. :param Vm: membrane potential (mV) :return: asymptotic probability (-) ''' return 1.0 / (1 + np.exp(-(Vm + 35) / 10)) # prob def taup(self, Vm): ''' Compute the decay time constant for adaptation of slow non-inactivating Potassium channels. :param Vm: membrane potential (mV) :return: decayed time constant (s) ''' return self.TauMax / (3.3 * np.exp((Vm + 35) / 20) + np.exp(-(Vm + 35) / 20)) # s def derM(self, Vm, m): ''' Compute the evolution of the open-probability of Sodium channels. :param Vm: membrane potential (mV) :param m: open-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alpham(Vm) * (1 - m) - self.betam(Vm) * m def derH(self, Vm, h): ''' Compute the evolution of the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :param h: inactivation-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphah(Vm) * (1 - h) - self.betah(Vm) * h def derN(self, Vm, n): ''' Compute the evolution of the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :param n: open-probability of delayed-rectifier Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphan(Vm) * (1 - n) - self.betan(Vm) * n def derP(self, Vm, p): ''' Compute the evolution of the open-probability of slow non-inactivating Potassium channels. :param Vm: membrane potential (mV) :param p: open-probability of slow non-inactivating Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.pinf(Vm) - p) / self.taup(Vm) def currNa(self, m, h, Vm): ''' Compute the inward Sodium current per unit area. :param m: open-probability of Sodium channels :param h: inactivation-probability of Sodium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GNa = self.GNaMax * m**3 * h return GNa * (Vm - self.VNa) def currK(self, n, Vm): ''' Compute the outward, delayed-rectifier Potassium current per unit area. :param n: open-probability of delayed-rectifier Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GKMax * n**4 return GK * (Vm - self.VK) def currM(self, p, Vm): ''' Compute the outward, slow non-inactivating Potassium current per unit area. :param p: open-probability of the slow non-inactivating Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GM = self.GMMax * p return GM * (Vm - self.VK) def currL(self, Vm): ''' Compute the non-specific leakage current per unit area. :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GL * (Vm - self.VL) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, p = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currM(p, Vm) + self.currL(Vm)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) peq = self.pinf(Vm) return np.array([meq, heq, neq, peq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, p = states dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dpdt = self.derP(Vm, p) return [dmdt, dhdt, dndt, dpdt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) Tp = self.taup(Vm) pinf = self.pinf(Vm) ap_avg = np.mean(pinf / Tp) bp_avg = np.mean(1 / Tp) - ap_avg # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, ap_avg, bp_avg]) - def derStatesEff(self, Adrive, Qm, states, interpolators): + def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' - rates = np.array([interpolators[rn](Adrive, Qm) for rn in self.coeff_names]) + 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, Adrive, Qm, states, interpolators): + 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(Adrive, Qm, NaK_states, interpolators) + NaK_dstates = super().derStatesEff(Qm, NaK_states, interp_data) # Compute Calcium channels states derivatives - Ca_rates = np.array([interpolators[rn](Adrive, Qm) for rn in self.coeff_names[8:]]) + 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/channels/leech.py b/PointNICE/channels/leech.py index 43751d7..9df80d6 100644 --- a/PointNICE/channels/leech.py +++ b/PointNICE/channels/leech.py @@ -1,390 +1,363 @@ #!/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-25 18:21:20 +# @Last Modified time: 2017-08-28 17:57:49 ''' Channels mechanisms for leech ganglion neurons. ''' import logging from functools import partialmethod import numpy as np from .base import BaseMech # Get package logger logger = logging.getLogger('PointNICE') class LeechTouch(BaseMech): ''' Class defining the membrane channel dynamics of a leech touch sensory neuron. with 4 different current types: - Inward Sodium current - Outward Potassium current - Inward Calcium current - Non-specific leakage current - Calcium-dependent, outward Potassium current - Outward, Sodium pumping current Reference: *Cataldo, E., Brunelli, M., Byrne, J.H., Av-Ron, E., Cai, Y., and Baxter, D.A. (2005). Computational model of touch sensory cells (T Cells) of the leech: role of the afterhyperpolarization (AHP) in activity-dependent conduction failure. J Comput Neurosci 18, 5–24.* ''' # Name of channel mechanism name = 'LeechT' # Cell-specific biophysical parameters Cm0 = 1e-2 # Cell membrane resting capacitance (F/m2) Vm0 = -53.0 # Cell membrane resting potential (mV) VNa = 45.0 # Sodium Nernst potential (mV) VK = -62.0 # Potassium Nernst potential (mV) VCa = 60.0 # Calcium Nernst potential (mV) VKCa = -62.0 # Calcium-dependent, Potassium current Nernst potential (mV) VL = -48.0 # Non-specific leakage Nernst potential (mV) VPumpNa = -300.0 # Sodium pump current reversal potential (mV) GNaMax = 3500.0 # Max. conductance of Sodium current (S/m^2) GKMax = 900.0 # Max. conductance of Potassium current (S/m^2) GCaMax = 20.0 # Max. conductance of Calcium current (S/m^2) GKCaMax = 236.0 # Max. conductance of Calcium-dependent Potassium current (S/m^2) GL = 1.0 # Conductance of non-specific leakage current (S/m^2) GPumpNa = 20.0 # Max. conductance of Sodium pump current (S/m^2) taum = 0.1e-3 # Sodium activation time constant (s) taus = 0.6e-3 # Calcium activation time constant (s) surface = 6434.0 # surface of cell assumed as a single soma (um2) curr_factor = 1e6 # 1/nA to 1/mA tau_Na_removal = 16.0 # Time constant for the removal of Sodium ions from the pool (s) tau_Ca_removal = 1.25 # Time constant for the removal of Calcium ions from the pool (s) tau_PumpNa_act = 0.1 # Time constant for the PumpNa current activation from Sodium ions (s) tau_KCa_act = 0.01 # Time constant for the KCa current activation from Calcium ions (s) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h'], 'i_K\ kin.': ['n'], 'i_{Ca}\ kin.': ['s'], 'pools': ['C_Na_arb', 'C_Na_arb_activation', 'C_Ca_arb', 'C_Ca_arb_activation'], '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 = 0.016 * self.surface / self.curr_factor self.K_Ca = 0.1 * self.surface / self.curr_factor # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) # ----------------- Generic ----------------- def _xinf(self, Vm, halfmax, slope, power): ''' Generic function computing the steady-state activation/inactivation of a particular ion channel at a given voltage. :param Vm: membrane potential (mV) :param halfmax: half-(in)activation voltage (mV) :param slope: slope parameter of (in)activation function (mV) :param power: power exponent multiplying the exponential expression (integer) :return: steady-state (in)activation (-) ''' return 1 / (1 + np.exp((Vm - halfmax) / slope))**power def _taux(self, Vm, halfmax, slope, tauMax, tauMin): ''' Generic function computing the voltage-dependent, activation/inactivation time constant of a particular ion channel at a given voltage. :param Vm: membrane potential (mV) :param halfmax: voltage at which (in)activation time constant is half-maximal (mV) :param slope: slope parameter of (in)activation time constant function (mV) :return: steady-state (in)activation (-) ''' return (tauMax - tauMin) / (1 + np.exp((Vm - halfmax) / slope)) + tauMin def _derC_ion(self, Cion, Iion, Kion, tau): ''' Generic function computing the time derivative of the concentration of a specific ion in its intracellular pool. :param Cion: ion concentration in the pool (arbitrary unit) :param Iion: ionic current (mA/m2) :param Kion: scaling factor for current contribution to pool (arb. unit / nA???) :param tau: time constant for removal of ions from the pool (s) :return: variation of ionic concentration in the pool (arbitrary unit /s) ''' return (Kion * (-Iion) - Cion) / tau def _derA_ion(self, Aion, Cion, tau): ''' Generic function computing the time derivative of the concentration and time dependent activation function, for a specific pool-dependent ionic current. :param Aion: concentration and time dependent activation function (arbitrary unit) :param Cion: ion concentration in the pool (arbitrary unit) :param tau: time constant for activation function variation (s) :return: variation of activation function (arbitrary unit / s) ''' return (Cion - Aion) / tau # ------------------ Na ------------------- minf = partialmethod(_xinf, halfmax=-35.0, slope=-5.0, power=1) hinf = partialmethod(_xinf, halfmax=-50.0, slope=9.0, power=2) tauh = partialmethod(_taux, halfmax=-36.0, slope=3.5, tauMax=14.0e-3, tauMin=0.2e-3) def derM(self, Vm, m): ''' Instantaneous derivative of Sodium activation. ''' return (self.minf(Vm) - m) / self.taum # s-1 def derH(self, Vm, h): ''' Instantaneous derivative of Sodium inactivation. ''' return (self.hinf(Vm) - h) / self.tauh(Vm) # s-1 # ------------------ K ------------------- ninf = partialmethod(_xinf, halfmax=-22.0, slope=-9.0, power=1) taun = partialmethod(_taux, halfmax=-10.0, slope=10.0, tauMax=6.0e-3, tauMin=1.0e-3) def derN(self, Vm, n): ''' Instantaneous derivative of Potassium activation. ''' return (self.ninf(Vm) - n) / self.taun(Vm) # s-1 # ------------------ Ca ------------------- sinf = partialmethod(_xinf, halfmax=-10.0, slope=-2.8, power=1) def derS(self, Vm, s): ''' Instantaneous derivative of Calcium activation. ''' return (self.sinf(Vm) - s) / self.taus # s-1 # ------------------ Pools ------------------- def derC_Na(self, C_Na, I_Na): ''' Derivative of Sodium concentration in intracellular pool. ''' return self._derC_ion(C_Na, I_Na, self.K_Na, self.tau_Na_removal) def derA_Na(self, A_Na, C_Na): ''' Derivative of Sodium pool-dependent activation function for iPumpNa. ''' return self._derA_ion(A_Na, C_Na, self.tau_PumpNa_act) def derC_Ca(self, C_Ca, I_Ca): ''' Derivative of Calcium concentration in intracellular pool. ''' return self._derC_ion(C_Ca, I_Ca, self.K_Ca, self.tau_Ca_removal) def derA_Ca(self, A_Ca, C_Ca): ''' Derivative of Calcium pool-dependent activation function for iKCa. ''' return self._derA_ion(A_Ca, C_Ca, self.tau_KCa_act) # ------------------ Currents ------------------- def currNa(self, m, h, Vm): ''' Sodium inward current. ''' return self.GNaMax * m**3 * h * (Vm - self.VNa) def currK(self, n, Vm): ''' Potassium outward current. ''' return self.GKMax * n**2 * (Vm - self.VK) def currCa(self, s, Vm): ''' Calcium inward current. ''' return self.GCaMax * s * (Vm - self.VCa) def currKCa(self, A_Ca, Vm): ''' Calcium-activated Potassium outward current. ''' return self.GKCaMax * A_Ca * (Vm - self.VKCa) def currPumpNa(self, A_Na, Vm): ''' Outward current mimicking the activity of the NaK-ATPase pump. ''' return self.GPumpNa * A_Na * (Vm - self.VPumpNa) def currL(self, Vm): ''' Leakage current. ''' return self.GL * (Vm - self.VL) def currNet(self, Vm, states): - ''' Compute the net ionic current per unit area. - - :param Vm: membrane potential (mV) - :states: state probabilities of the ion channels - :return: current per unit area (mA/m2) - ''' + ''' 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): - ''' Compute the Sodium, Potassium and Calcium channel gates steady-states - for a specific membrane potential value. - - :param Vm: membrane potential (mV) - :return: array of steady-states - ''' + ''' 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 # print('initial Na current: {:.2f} mA/m2'.format(INa_eq)) # print('initial Na concentration in pool: {:.2f} arb. unit'.format(CNa_eq)) # KCa current pool concentration and activation steady-state ICa_eq = self.currCa(seq, Vm) CCa_eq = self.K_Ca * (-ICa_eq) ACa_eq = CCa_eq # print('initial Ca current: {:.2f} mA/m2'.format(ICa_eq)) # print('initial Ca concentration in pool: {:.2f} arb. unit'.format(CCa_eq)) return np.array([meq, heq, neq, seq, CNa_eq, ANa_eq, CCa_eq, ACa_eq]) def derStates(self, Vm, states): - ''' Compute the derivatives of channel states. - - :param Vm: membrane potential (mV) - :states: state probabilities of the ion channels - :return: current per unit area (mA/m2) - ''' + ''' 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): - ''' Get the effective rate constants of ion channels, averaged along an acoustic cycle, - for future use in effective simulations. - - :param Vm: array of membrane potential values for an acoustic cycle (mV) - :return: an array of rate average constants (s-1) - ''' + ''' 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, Adrive, Qm, states, interpolators): - ''' Compute the effective derivatives of channel states, based on - 2-dimensional linear interpolation of "effective" coefficients - that summarize the system's behaviour over an acoustic cycle. - - :param Vm_eff: effective membrane potential (mV) - :states: state probabilities of the ion channels - :param interpolators: dictionary of 2-dimensional linear interpolators - of "effective" rates over the 2D amplitude x charge input domain. - ''' + def derStatesEff(self, Qm, states, interp_data): + ''' Concrete implementation of the abstract API method. ''' - rates = np.array([interpolators[rn](Adrive, Qm) for rn in self.coeff_names]) - Vmeff = interpolators['V'](Adrive, Qm) + 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 - m) - rates[7] * s # PumpNa current pool concentration and activation state I_Na = self.currNa(m, h, Vmeff) dCNa_dt = self.derC_Na(C_Na, I_Na) dANa_dt = self.derA_Na(A_Na, C_Na) # KCa current pool concentration and activation state I_Ca_eff = self.currCa(s, Vmeff) dCCa_dt = self.derC_Ca(C_Ca, I_Ca_eff) dACa_dt = self.derA_Ca(A_Ca, C_Ca) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dCNa_dt, dANa_dt, dCCa_dt, dACa_dt] diff --git a/PointNICE/channels/thalamic.py b/PointNICE/channels/thalamic.py index a83ea99..bc31692 100644 --- a/PointNICE/channels/thalamic.py +++ b/PointNICE/channels/thalamic.py @@ -1,790 +1,793 @@ #!/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-25 18:04:26 +# @Last Modified time: 2017-08-28 17:55:32 ''' 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'] def alpham(self, Vm): ''' Compute the alpha rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = (-0.32 * (Vdiff - 13) / (np.exp(- (Vdiff - 13) / 4) - 1)) # ms-1 return alpha * 1e3 # s-1 def betam(self, Vm): ''' Compute the beta rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = (0.28 * (Vdiff - 40) / (np.exp((Vdiff - 40) / 5) - 1)) # ms-1 return beta * 1e3 # s-1 def alphah(self, Vm): ''' Compute the alpha rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = (0.128 * np.exp(-(Vdiff - 17) / 18)) # ms-1 return alpha * 1e3 # s-1 def betah(self, Vm): ''' Compute the beta rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = (4 / (1 + np.exp(-(Vdiff - 40) / 5))) # ms-1 return beta * 1e3 # s-1 def alphan(self, Vm): ''' Compute the alpha rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT alpha = (-0.032 * (Vdiff - 15) / (np.exp(-(Vdiff - 15) / 5) - 1)) # ms-1 return alpha * 1e3 # s-1 def betan(self, Vm): ''' Compute the beta rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' Vdiff = Vm - self.VT beta = (0.5 * np.exp(-(Vdiff - 10) / 40)) # ms-1 return beta * 1e3 # s-1 def derM(self, Vm, m): ''' Compute the evolution of the open-probability of Sodium channels. :param Vm: membrane potential (mV) :param m: open-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alpham(Vm) * (1 - m) - self.betam(Vm) * m def derH(self, Vm, h): ''' Compute the evolution of the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :param h: inactivation-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphah(Vm) * (1 - h) - self.betah(Vm) * h def derN(self, Vm, n): ''' Compute the evolution of the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :param n: open-probability of delayed-rectifier Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphan(Vm) * (1 - n) - self.betan(Vm) * n def derS(self, Vm, s): ''' Compute the evolution of the open-probability of the S-type, activation gate of Calcium channels. :param Vm: membrane potential (mV) :param s: open-probability of S-type Calcium activation gates (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.sinf(Vm) - s) / self.taus(Vm) def derU(self, Vm, u): ''' Compute the evolution of the open-probability of the U-type, inactivation gate of Calcium channels. :param Vm: membrane potential (mV) :param u: open-probability of U-type Calcium inactivation gates (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.uinf(Vm) - u) / self.tauu(Vm) def currNa(self, m, h, Vm): ''' Compute the inward Sodium current per unit area. :param m: open-probability of Sodium channels :param h: inactivation-probability of Sodium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GNa = self.GNaMax * m**3 * h return GNa * (Vm - self.VNa) def currK(self, n, Vm): ''' Compute the outward delayed-rectifier Potassium current per unit area. :param n: open-probability of delayed-rectifier Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GKMax * n**4 return GK * (Vm - self.VK) def currCa(self, s, u, Vm): ''' Compute the inward Calcium current per unit area. :param s: open-probability of the S-type activation gate of Calcium channels :param u: open-probability of the U-type inactivation gate of Calcium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GT = self.GTMax * s**2 * u return GT * (Vm - self.VCa) def currL(self, Vm): ''' Compute the non-specific leakage current per unit area. :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GL * (Vm - self.VL) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, u = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currCa(s, u, Vm) + self.currL(Vm)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) seq = self.sinf(Vm) ueq = self.uinf(Vm) return np.array([meq, heq, neq, seq, ueq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, u = states dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) dudt = self.derU(Vm, u) return [dmdt, dhdt, dndt, dsdt, dudt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) Ts = self.taus(Vm) sinf = self.sinf(Vm) as_avg = np.mean(sinf / Ts) bs_avg = np.mean(1 / Ts) - as_avg Tu = np.array([self.tauu(v) for v in Vm]) uinf = self.uinf(Vm) au_avg = np.mean(uinf / Tu) bu_avg = np.mean(1 / Tu) - au_avg # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg, au_avg, bu_avg]) - def derStatesEff(self, Adrive, Qm, states, interpolators): + + def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' - rates = np.array([interpolators[rn](Adrive, Qm) for rn in self.coeff_names]) + 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 vars_RE = { '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, Adrive, Qm, states, interpolators): + def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' - rates = np.array([interpolators[rn](Adrive, Qm) for rn in self.coeff_names]) - Vmeff = interpolators['V'](Adrive, Qm) + 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 71b865d..cd04dcb 100644 --- a/PointNICE/solvers/SolverUS.py +++ b/PointNICE/solvers/SolverUS.py @@ -1,706 +1,710 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-09-29 16:16:19 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-08-28 15:20:07 +# @Last Modified time: 2017-08-28 17:47:33 import os import warnings import pickle import logging 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 ..channels 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, geom, params, ch_mech, Fdrive): """ Constructor of the class. :param geom: BLS geometric constants dictionary :param params: BLS biomechanical and biophysical parameters dictionary :param ch_mech: channels mechanism object :param Fdrive: frequency of acoustic perturbation (Hz) """ # Check validity of input parameters assert isinstance(ch_mech, BaseMech), ('channel 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 = ch_mech.Cm0 Vm0 = ch_mech.Vm0 BilayerSonophore.__init__(self, geom, params, Fdrive, Cm0, Cm0 * Vm0 * 1e-3) logger.info('US solver initialization with %s channel mechanism', ch_mech.name) def eqHH(self, t, y, ch_mech, 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 ch_mech: channels mechanism 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 = - ch_mech.currNet(Vm, states) * 1e-3 # A/m2 dstates = ch_mech.derStates(Vm, states) # Return derivatives vector return [dQm, *dstates] - def eqHHeff(self, t, y, ch_mech, A, interpolators): - """ Compute the derivatives of the n-ODE effective HH system variables, - based on 2-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 ch_mech: channels mechanism object - :param A: acoustic drive amplitude (Pa) - :param channels: Channel object to compute a specific electrical membrane dynamics - :param interpolators: dictionary of 2-dimensional linear interpolators - of "effective" coefficients over the 2D amplitude x charge input domain. - :return: vector of effective system derivatives at time t - """ - - # Split input vector explicitly - Qm, *states = y - - # Compute charge and channel states variation - Vm = interpolators['V'](A, Qm) # mV - dQmdt = - ch_mech.currNet(Vm, states) * 1e-3 - dstates = ch_mech.derStatesEff(A, Qm, states, interpolators) - - # Return derivatives vector - return [dQmdt, *dstates] - - def eqFull(self, t, y, ch_mech, 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 ch_mech: channels mechanism 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:], ch_mech, self.Capct(y[1])) # return concatenated output return dydt_mech + dydt_elec + def eqHHeff(self, t, y, ch_mech, 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 ch_mech: channels mechanism object + :param channels: Channel object to compute a specific electrical membrane dynamics + :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 = - ch_mech.currNet(Vm, states) * 1e-3 + dstates = ch_mech.derStatesEff(Qm, states, interp_data) + + # Return derivatives vector + return [dQmdt, *dstates] + + def getEffCoeffs(self, ch_mech, 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 ch_mech: channels mechanism 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, rates, and gas content coefficients """ # 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 = ch_mech.getEffRates(Vm) # Take final cycle value for gas content ng_eff = ng[-1] # mole return (Vm_eff, rates_eff, ng_eff) def createLookup(self, ch_mech, Fdrive, amps, charges, 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 ch_mech: channels mechanism object :param Fdrive: acoustic drive frequency (Hz) :param amps: array of acoustic drive amplitudes (Pa) :param charges: array of charge densities (C/m2) :param phi: acoustic drive phase (rad) """ # Check validity of stimulation parameters assert Fdrive > 0, 'Driving frequency must be strictly positive' assert np.amin(amps) >= 0, 'Acoustic pressure amplitudes must be positive' logger.info('Creating lookup table for f = %.2f kHz', Fdrive * 1e-3) # Initialize 3D array to store effective coefficients nA = amps.size nQ = charges.size Vm = np.empty((nA, nQ)) ng = np.empty((nA, nQ)) nrates = len(ch_mech.coeff_names) rates = np.empty((nA, nQ, nrates)) # Loop through all (A, Q) combinations isim = 0 for i in range(nA): for j in range(nQ): isim += 1 # Run short simulation and store effective coefficients logger.info('sim %u/%u (A = %.2f kPa, Q = %.2f nC/cm2)', isim, nA * nQ, amps[i] * 1e-3, charges[j] * 1e5) (Vm[i, j], rates[i, j, :], ng[i, j]) = self.getEffCoeffs(ch_mech, Fdrive, amps[i], charges[j], phi) # Convert coefficients array into dictionary with specific names lookup_dict = {ch_mech.coeff_names[k]: rates[:, :, k] for k in range(nrates)} lookup_dict['V'] = Vm # mV lookup_dict['ng'] = ng # mole # Add input amplitude and charge arrays to dictionary lookup_dict['A'] = amps # Pa lookup_dict['Q'] = charges # C/m2 # Save dictionary in lookup file lookup_file = '{}_lookups_a{:.1f}nm_f{:.1f}kHz.pkl'.format(ch_mech.name, self.a * 1e9, Fdrive * 1e-3) logger.info('Saving effective coefficients arrays in lookup file: "%s"', lookup_file) lookup_filepath = '{0}/{1}/{2}'.format(getLookupDir(), ch_mech.name, lookup_file) with open(lookup_filepath, 'wb') as fh: pickle.dump(lookup_dict, fh) def __runClassic(self, ch_mech, Fdrive, Adrive, tstim, toffset, PRF, DF, 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 ch_mech: channels mechanism 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 DF: pulse duty factor (-) :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 # Determine proportion of tstim in total integration stim_prop = tstim / (tstim + toffset) # if CW stimulus: divide integration during stimulus into 100 intervals if DF == 1.0: PRF = 100 / tstim # Compute vector sizes npulses = int(np.round(PRF * tstim)) Tpulse_on = DF / PRF Tpulse_off = (1 - DF) / 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(ch_mech.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))) # Loop through all pulse (ON and OFF) intervals for i in range(npulses): # logger.debug('pulse %u/%u', i + 1, npulses) printPct(100 * stim_prop * (i + 1) / npulses, 1) # 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(ch_mech, 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(ch_mech, 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) # 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(ch_mech, 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) # 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, ch_mech, Fdrive, Adrive, tstim, toffset, PRF, DF, 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 ch_mech: channels mechanism 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 DF: pulse duty factor (-) :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_f{:.1f}kHz.pkl'.format(ch_mech.name, self.a * 1e9, Fdrive * 1e-3) lookup_path = '{}/{}/{}'.format(getLookupDir(), ch_mech.name, lookup_file) assert os.path.isfile(lookup_path), 'No lookup file for this stimulation frequency' # Load coefficients with open(lookup_path, 'rb') as fh: - coeffs = pickle.load(fh) + coeffs2d = pickle.load(fh) # Check that pressure amplitude is within lookup range - Amax = np.amax(coeffs['A']) + 1e-9 # adding margin to compensate for eventual round error + Amax = np.amax(coeffs2d['A']) + 1e-9 # adding margin to compensate for eventual round error assert Adrive <= Amax, 'Amplitude must be within [0, {:.1f}] kPa'.format(Amax * 1e-3) - # Initialize interpolators - interpolators = {cn: interp2d(coeffs['A'], coeffs['Q'], np.transpose(coeffs[cn])) - for cn in ch_mech.coeff_names} - interpolators['V'] = interp2d(coeffs['A'], coeffs['Q'], np.transpose(coeffs['V'])) - interpolators['ng'] = interp2d(coeffs['A'], coeffs['Q'], np.transpose(coeffs['ng'])) + # Derive 1D charge-based interpolation dataset for specific amplitude + coeffs_list = ['V', 'ng', *ch_mech.coeff_names] + interpolators2d = {cn: interp2d(coeffs2d['A'], coeffs2d['Q'], np.transpose(coeffs2d[cn])) + for cn in coeffs_list} + coeffs1d = {cn: interpolators2d[cn](Adrive, coeffs2d['Q']) for cn in coeffs_list} + coeffs1d['Q'] = coeffs2d['Q'] + coeffs1d['ng0'] = interpolators2d['ng'](0.0, coeffs2d['Q']) + for key in coeffs1d.keys(): + coeffs1d[key] = np.squeeze(coeffs1d[key]) + # Initialize system solvers solver_on = integrate.ode(self.eqHHeff) solver_on.set_integrator('lsoda', nsteps=SOLVER_NSTEPS) - solver_on.set_f_params(ch_mech, Adrive, interpolators) + solver_on.set_f_params(ch_mech, 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 DF == 1.0: PRF = 1 / tstim # Compute vector sizes npulses = int(np.round(PRF * tstim)) Tpulse_on = DF / PRF Tpulse_off = (1 - DF) / PRF n_pulse_on = int(np.round(Tpulse_on / dt)) + 1 n_pulse_off = int(np.round(Tpulse_off / dt)) 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(ch_mech.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] = interpolators['ng'](Adrive, y_pulse[0, k]) # mole + 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(ch_mech, 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] = interpolators['ng'](0.0, y_pulse[0, k]) # mole + 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(ch_mech, 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(ch_mech, 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] = interpolators['ng'](0.0, y_off[0, k]) # mole + 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(ch_mech, 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, ch_mech, 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 ch_mech: channels mechanism 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(ch_mech, 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(ch_mech.states0, (2, 1)).T y = np.vstack((y_membrane, y_channels)) nvar = y.shape[0] # 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(ch_mech, 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 assert (y_full[1, k] > -0.5 * self.Delta), 'Deflection out of range' k += 1 except (Warning, AssertionError) as inst: sim_error = True logger.error('Full system integration error at step %u', k) print(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:] # print('convergence after {} cycles'.format(j)) # 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(ch_mech, 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) print(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 # Print progress printPct(100 * (t[-1] / (tstim + toffset)), 1) irep += 1 # Return output return (t, y[1:, :], states) def run(self, ch_mech, Fdrive, Adrive, tstim, toffset, PRF=None, DF=1.0, sim_type='effective'): """ Run simulation of the system for a specific set of US stimulation parameters. :param ch_mech: channels mechanism 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 DF: pulse duty factor (-) :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 + # # 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 - # Check validity of input parameters assert isinstance(ch_mech, BaseMech), ('channel mechanism must be inherited ' 'from the BaseMech class') for param in [Fdrive, Adrive, tstim, toffset, DF]: 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 strictly positive' assert tstim > 0, 'Stimulus duration must be strictly positive' assert toffset >= 0, 'Stimulus offset must be positive or null' assert DF > 0 and DF <= 1, 'Duty cycle must be within [0; 1)' if DF < 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(ch_mech, Fdrive, Adrive, tstim, toffset, PRF, DF) elif sim_type == 'effective': return self.__runEffective(ch_mech, Fdrive, Adrive, tstim, toffset, PRF, DF) elif sim_type == 'hybrid': assert DF == 1.0, 'Hybrid method can only handle continuous wave stimuli' return self.__runHybrid(ch_mech, Fdrive, Adrive, tstim, toffset) + diff --git a/tests/test_basic.py b/tests/test_basic.py index a65623d..30e1d07 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -1,271 +1,271 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-06-14 18:37:45 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-08-28 14:18:56 +# @Last Modified time: 2017-08-28 18:00:05 ''' Test the basic functionalities of the package. ''' import logging import inspect import numpy as np from PointNICE.utils import load_BLS_params from PointNICE import BilayerSonophore, channels, SolverElec, SolverUS from PointNICE.solvers import detectSpikes, titrateEStim, titrateAStim from PointNICE.constants import * # Set logging options logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%d/%m/%Y %H:%M:%S:') logger = logging.getLogger('PointNICE') logger.setLevel(logging.DEBUG) # List of implemented neurons neurons = [] for _, obj in inspect.getmembers(channels): if inspect.isclass(obj) and isinstance(obj.name, str): if obj.name != 'LeechT': neurons.append(obj) def test_MECH(): ''' Maximal negative and positive deflections of the BLS structure for a specific sonophore size, resting membrane properties and stimulation parameters. ''' logger.info('Starting test: Mechanical simulation') # BLS geometry and parameters geom = {"a": 32e-9, "d": 0.0e-6} params = load_BLS_params() # Create BLS instance Fdrive = 350e3 # Hz Cm0 = 1e-2 # membrane resting capacitance (F/m2) Qm0 = -80e-5 # membrane resting charge density (C/m2) bls = BilayerSonophore(geom, params, Fdrive, Cm0, Qm0) # Run mechanical simulation Adrive = 100e3 # Pa Qm = 50e-5 # C/m2 bls.runMech(Fdrive, Adrive, Qm) _, y, _ = bls.runMech(Fdrive, Adrive, Qm) # Check validity of deflection extrema Z, _ = y Zlast = Z[-NPC_FULL:] Zmin, Zmax = (Zlast.min(), Zlast.max()) logger.info('Zmin = %.2f nm, Zmax = %.2f nm', Zmin * 1e9, Zmax * 1e9) Zmin_ref, Zmax_ref = (-0.116e-9, 5.741e-9) assert np.abs(Zmin - Zmin_ref) < 1e-12, 'Unexpected BLS compression amplitude' assert np.abs(Zmax - Zmax_ref) < 1e-12, 'Unexpected BLS expansion amplitude' logger.info('Passed test: Mechanical simulation') def test_resting_potential(): ''' Neurons membrane potential in free conditions should stabilize to their specified resting potential value. ''' conv_err_msg = ('{} neuron membrane potential in free conditions does not converge to ' 'stable value (gap after 20s: {:.2e} mV)') value_err_msg = ('{} neuron steady-state membrane potential in free conditions differs ' 'significantly from specified resting potential (gap = {:.2f} mV)') logger.info('Starting test: neurons in free conditions') # Initialize solver solver = SolverElec() for neuron_class in neurons: # Simulate each neuron in free conditions neuron = neuron_class() _, y, _ = solver.run(neuron, Astim=0.0, tstim=20.0, toffset=0.0) Vm_free, *_ = y # Check membrane potential convergence Vm_free_last, Vm_free_beforelast = (Vm_free[-1], Vm_free[-2]) Vm_free_conv = Vm_free_last - Vm_free_beforelast assert np.abs(Vm_free_conv) < 1e-5, conv_err_msg.format(neuron.name, Vm_free_conv) # Check membrane potential convergence to resting potential Vm_free_diff = Vm_free_last - neuron.Vm0 assert np.abs(Vm_free_diff) < 0.1, value_err_msg.format(neuron.name, Vm_free_diff) logger.info('Passed test: neurons in free conditions') def test_ESTIM(): ''' Threshold E-STIM amplitude and needed to obtain an action potential and response latency should match reference values. ''' Athr_err_msg = ('{} neuron threshold amplitude for excitation does not match reference value' '(gap = {:.2f} mA/m2)') latency_err_msg = ('{} neuron latency for excitation at threshold amplitude does not match ' 'reference value (gap = {:.2f} ms)') logger.info('Starting test: E-STIM titration') # Initialize solver solver = SolverElec() Arange = (0.0, 2 * TITRATION_ESTIM_A_MAX) # mA/m2 # Stimulation parameters tstim = 100e-3 # s toffset = 50e-3 # s # Reference values Athr_refs = {'FS': 6.91, 'LTS': 1.54, 'RS': 5.03, 'LeechT': 5.54, 'RE': 3.61, 'TC': 4.05} latency_refs = {'FS': 101.00e-3, 'LTS': 128.56e-3, 'RS': 103.81e-3, 'LeechT': 20.22e-3, 'RE': 148.50e-3, 'TC': 63.46e-3} for neuron_class in neurons: # Perform titration for each neuron neuron = neuron_class() print(neuron.name, 'neuron titration') (Athr, t, y, _, latency) = titrateEStim(solver, neuron, Arange, tstim, toffset, PRF=None, DF=1.0) Vm = y[0, :] # Check that final number of spikes is 1 n_spikes, _, _ = detectSpikes(t, Vm, SPIKE_MIN_VAMP, SPIKE_MIN_DT) assert n_spikes == 1, 'Number of spikes after titration should be exactly 1' # Check threshold amplitude Athr_diff = Athr - Athr_refs[neuron.name] assert np.abs(Athr_diff) < 0.1, Athr_err_msg.format(neuron.name, Athr_diff) # Check response latency lat_diff = (latency - latency_refs[neuron.name]) * 1e3 assert np.abs(lat_diff) < 1.0, latency_err_msg.format(neuron.name, lat_diff) logger.info('Passed test: E-STIM titration') def test_ASTIM(): ''' Threshold A-STIM amplitude and needed to obtain an action potential and response latency should match reference values. ''' Athr_err_msg = ('{} neuron threshold amplitude for excitation does not match reference value' '(gap = {:.2f} kPa)') latency_err_msg = ('{} neuron latency for excitation at threshold amplitude does not match ' 'reference value (gap = {:.2f} ms)') logger.info('Starting test: A-STIM titration') # BLS geometry and parameters geom = {"a": 32e-9, "d": 0.0e-6} params = load_BLS_params() # Stimulation parameters Fdrive = 350e3 # Hz tstim = 50e-3 # s toffset = 30e-3 # s Arange = (0.0, 2 * TITRATION_ASTIM_A_MAX) # Pa # Reference values Athr_refs = {'FS': 38.67e3, 'LTS': 24.80e3, 'RS': 51.17e3, 'RE': 46.36e3, 'TC': 23.24e3, 'LeechT': None} latency_refs = {'FS': 63.72e-3, 'LTS': 61.92e-3, 'RS': 62.52e-3, 'RE': 79.20e-3, 'TC': 68.53e-3, 'LeechT': None} # Titration for each neuron for neuron_class in neurons: # Initialize neuron neuron = neuron_class() print(neuron.name, 'neuron titration') # Initialize solver solver = SolverUS(geom, params, neuron, Fdrive) # Perform titration (Athr, t, y, _, latency) = titrateAStim(solver, neuron, Fdrive, Arange, tstim, toffset, PRF=None, DF=1.0) Qm = y[2] # Check that final number of spikes is 1 n_spikes, _, _ = detectSpikes(t, Qm, SPIKE_MIN_QAMP, SPIKE_MIN_DT) assert n_spikes == 1, 'Number of spikes after titration should be exactly 1' # Check threshold amplitude Athr_diff = (Athr - Athr_refs[neuron.name]) * 1e-3 assert np.abs(Athr_diff) < 0.1, Athr_err_msg.format(neuron.name, Athr_diff) # Check response latency lat_diff = (latency - latency_refs[neuron.name]) * 1e3 assert np.abs(lat_diff) < 1.0, latency_err_msg.format(neuron.name, lat_diff) logger.info('Passed test: A-STIM titration') def test_ASTIM_classic(): ''' Short classic A-STIM simulation should not raise any errors. ''' logger.info('Starting test: A-STIM classic simulation') # BLS geometry and parameters geom = {"a": 32e-9, "d": 0.0e-6} params = load_BLS_params() # Stimulation parameters Fdrive = 350e3 # Hz Adrive = 100e3 # Pa tstim = 1e-3 # s toffset = 1e-3 # s # Initialize RS neuron rs_neuron = channels.CorticalRS() # Initialize solver solver = SolverUS(geom, params, rs_neuron, Fdrive) # Run short classic simulation of the system solver.run(rs_neuron, Fdrive, Adrive, tstim, toffset, sim_type='classic') # If no error is raised, test is passed successfully logger.info('Passed test: A-STIM classic simulation') def test_ASTIM_hybrid(): ''' Short hybrid A-STIM simulation should not raise any errors. ''' logger.info('Starting test: A-STIM hybrid simulation') # BLS geometry and parameters geom = {"a": 32e-9, "d": 0.0e-6} params = load_BLS_params() # Stimulation parameters Fdrive = 350e3 # Hz Adrive = 100e3 # Pa tstim = 10e-3 # s toffset = 1e-3 # s # Initialize RS neuron rs_neuron = channels.CorticalRS() # Initialize solver solver = SolverUS(geom, params, rs_neuron, Fdrive) # Run short classic simulation of the system solver.run(rs_neuron, Fdrive, Adrive, tstim, toffset, sim_type='hybrid') # If no error is raised, test is passed successfully logger.info('Passed test: A-STIM hybrid simulation') if __name__ == '__main__': logger.info('Starting tests') test_MECH() test_resting_potential() test_ESTIM() test_ASTIM() test_ASTIM_classic() test_ASTIM_hybrid() logger.info('All tests successfully passed')