diff --git a/LICENSE b/LICENSE index c2b6e9a..cbe4c87 100644 --- a/LICENSE +++ b/LICENSE @@ -1,21 +1,21 @@ MIT License -Copyright (c) 2017 TNE lab, EPFL +Copyright (c) 2018 TNE lab, EPFL Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. \ No newline at end of file diff --git a/PointNICE/channels/leech.py b/PointNICE/channels/leech.py index 2170f6d..a2e8064 100644 --- a/PointNICE/channels/leech.py +++ b/PointNICE/channels/leech.py @@ -1,1074 +1,1074 @@ #!/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 -# @Last Modified time: 2018-03-02 12:46:33 +# @Last Modified by: Theo Lemaire +# @Last Modified time: 2018-03-07 14:46:54 ''' Channels mechanisms for leech ganglion neurons. ''' import logging from functools import partialmethod import numpy as np from .base import BaseMech # Get package logger logger = logging.getLogger('PointNICE') class LeechTouch(BaseMech): ''' Class defining the membrane channel dynamics of a leech touch sensory neuron. with 4 different current types: - Inward Sodium current - Outward Potassium current - Inward Calcium current - Non-specific leakage current - Calcium-dependent, outward Potassium current - Outward, Sodium pumping current Reference: *Cataldo, E., Brunelli, M., Byrne, J.H., Av-Ron, E., Cai, Y., and Baxter, D.A. (2005). Computational model of touch sensory cells (T Cells) of the leech: role of the afterhyperpolarization (AHP) in activity-dependent conduction failure. J Comput Neurosci 18, 5–24.* ''' # Name of channel mechanism name = 'LeechT' # Cell-specific biophysical parameters Cm0 = 1e-2 # Cell membrane resting capacitance (F/m2) Vm0 = -53.58 # Cell membrane resting potential (mV) VNa = 45.0 # Sodium Nernst potential (mV) VK = -62.0 # Potassium Nernst potential (mV) VCa = 60.0 # Calcium Nernst potential (mV) VL = -48.0 # Non-specific leakage Nernst potential (mV) VPumpNa = -300.0 # Sodium pump current reversal potential (mV) GNaMax = 3500.0 # Max. conductance of Sodium current (S/m^2) GKMax = 900.0 # Max. conductance of Potassium current (S/m^2) GCaMax = 20.0 # Max. conductance of Calcium current (S/m^2) GKCaMax = 236.0 # Max. conductance of Calcium-dependent Potassium current (S/m^2) GL = 1.0 # Conductance of non-specific leakage current (S/m^2) GPumpNa = 20.0 # Max. conductance of Sodium pump current (S/m^2) taum = 0.1e-3 # Sodium activation time constant (s) taus = 0.6e-3 # Calcium activation time constant (s) # Original conversion constants from inward ion current (nA) to build-up of # intracellular ion concentration (arb.) K_Na_original = 0.016 # iNa to intracellular [Na+] K_Ca_original = 0.1 # iCa to intracellular [Ca2+] # Constants needed to convert K from original model (soma compartment) # to current model (point-neuron) surface = 6434.0e-12 # surface of cell assumed as a single soma (m2) curr_factor = 1e6 # mA to nA # Time constants for the removal of ions from intracellular pools (s) tau_Na_removal = 16.0 # Na+ removal tau_Ca_removal = 1.25 # Ca2+ removal # Time constants for the iPumpNa and iKCa currents activation # from specific intracellular ions (s) tau_PumpNa_act = 0.1 # iPumpNa activation from intracellular Na+ tau_KCa_act = 0.01 # iKCa activation from intracellular Ca2+ # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h', 'm3h'], 'i_K\ kin.': ['n'], 'i_{Ca}\ kin.': ['s'], 'pools': ['C_Na_arb', 'C_Na_arb_activation', 'C_Ca_arb', 'C_Ca_arb_activation'], 'I': ['iNa', 'iK', 'iCa', 'iKCa', 'iPumpNa', 'iL', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'C_Na', 'A_Na', 'C_Ca', 'A_Ca'] self.states0 = np.array([]) # Names of the channels effective coefficients self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas'] self.K_Na = self.K_Na_original * self.surface * self.curr_factor self.K_Ca = self.K_Ca_original * self.surface * self.curr_factor # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) # ----------------- Generic ----------------- def _xinf(self, Vm, halfmax, slope, power): ''' Generic function computing the steady-state activation/inactivation of a particular ion channel at a given voltage. :param Vm: membrane potential (mV) :param halfmax: half-(in)activation voltage (mV) :param slope: slope parameter of (in)activation function (mV) :param power: power exponent multiplying the exponential expression (integer) :return: steady-state (in)activation (-) ''' return 1 / (1 + np.exp((Vm - halfmax) / slope))**power def _taux(self, Vm, halfmax, slope, tauMax, tauMin): ''' Generic function computing the voltage-dependent, activation/inactivation time constant of a particular ion channel at a given voltage. :param Vm: membrane potential (mV) :param halfmax: voltage at which (in)activation time constant is half-maximal (mV) :param slope: slope parameter of (in)activation time constant function (mV) :return: steady-state (in)activation (-) ''' return (tauMax - tauMin) / (1 + np.exp((Vm - halfmax) / slope)) + tauMin def _derC_ion(self, Cion, Iion, Kion, tau): ''' Generic function computing the time derivative of the concentration of a specific ion in its intracellular pool. :param Cion: ion concentration in the pool (arbitrary unit) :param Iion: ionic current (mA/m2) :param Kion: scaling factor for current contribution to pool (arb. unit / nA???) :param tau: time constant for removal of ions from the pool (s) :return: variation of ionic concentration in the pool (arbitrary unit /s) ''' return (Kion * (-Iion) - Cion) / tau def _derA_ion(self, Aion, Cion, tau): ''' Generic function computing the time derivative of the concentration and time dependent activation function, for a specific pool-dependent ionic current. :param Aion: concentration and time dependent activation function (arbitrary unit) :param Cion: ion concentration in the pool (arbitrary unit) :param tau: time constant for activation function variation (s) :return: variation of activation function (arbitrary unit / s) ''' return (Cion - Aion) / tau # ------------------ Na ------------------- minf = partialmethod(_xinf, halfmax=-35.0, slope=-5.0, power=1) hinf = partialmethod(_xinf, halfmax=-50.0, slope=9.0, power=2) tauh = partialmethod(_taux, halfmax=-36.0, slope=3.5, tauMax=14.0e-3, tauMin=0.2e-3) def derM(self, Vm, m): ''' Instantaneous derivative of Sodium activation. ''' return (self.minf(Vm) - m) / self.taum # s-1 def derH(self, Vm, h): ''' Instantaneous derivative of Sodium inactivation. ''' return (self.hinf(Vm) - h) / self.tauh(Vm) # s-1 # ------------------ K ------------------- ninf = partialmethod(_xinf, halfmax=-22.0, slope=-9.0, power=1) taun = partialmethod(_taux, halfmax=-10.0, slope=10.0, tauMax=6.0e-3, tauMin=1.0e-3) def derN(self, Vm, n): ''' Instantaneous derivative of Potassium activation. ''' return (self.ninf(Vm) - n) / self.taun(Vm) # s-1 # ------------------ Ca ------------------- sinf = partialmethod(_xinf, halfmax=-10.0, slope=-2.8, power=1) def derS(self, Vm, s): ''' Instantaneous derivative of Calcium activation. ''' return (self.sinf(Vm) - s) / self.taus # s-1 # ------------------ Pools ------------------- def derC_Na(self, C_Na, I_Na): ''' Derivative of Sodium concentration in intracellular pool. ''' return self._derC_ion(C_Na, I_Na, self.K_Na, self.tau_Na_removal) def derA_Na(self, A_Na, C_Na): ''' Derivative of Sodium pool-dependent activation function for iPumpNa. ''' return self._derA_ion(A_Na, C_Na, self.tau_PumpNa_act) def derC_Ca(self, C_Ca, I_Ca): ''' Derivative of Calcium concentration in intracellular pool. ''' return self._derC_ion(C_Ca, I_Ca, self.K_Ca, self.tau_Ca_removal) def derA_Ca(self, A_Ca, C_Ca): ''' Derivative of Calcium pool-dependent activation function for iKCa. ''' return self._derA_ion(A_Ca, C_Ca, self.tau_KCa_act) # ------------------ Currents ------------------- def currNa(self, m, h, Vm): ''' Sodium inward current. ''' return self.GNaMax * m**3 * h * (Vm - self.VNa) def currK(self, n, Vm): ''' Potassium outward current. ''' return self.GKMax * n**2 * (Vm - self.VK) def currCa(self, s, Vm): ''' Calcium inward current. ''' return self.GCaMax * s * (Vm - self.VCa) def currKCa(self, A_Ca, Vm): ''' Calcium-activated Potassium outward current. ''' return self.GKCaMax * A_Ca * (Vm - self.VK) def currPumpNa(self, A_Na, Vm): ''' Outward current mimicking the activity of the NaK-ATPase pump. ''' return self.GPumpNa * A_Na * (Vm - self.VPumpNa) def currL(self, Vm): ''' Leakage current. ''' return self.GL * (Vm - self.VL) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, _, A_Na, _, A_Ca = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currCa(s, Vm) + self.currL(Vm) + self.currPumpNa(A_Na, Vm) + self.currKCa(A_Ca, Vm)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Standard gating dynamics: Solve the equation dx/dt = 0 at Vm for each x-state meq = self.minf(Vm) heq = self.hinf(Vm) neq = self.ninf(Vm) seq = self.sinf(Vm) # PumpNa pool concentration and activation steady-state INa_eq = self.currNa(meq, heq, Vm) CNa_eq = self.K_Na * (-INa_eq) ANa_eq = CNa_eq # KCa current pool concentration and activation steady-state ICa_eq = self.currCa(seq, Vm) CCa_eq = self.K_Ca * (-ICa_eq) ACa_eq = CCa_eq return np.array([meq, heq, neq, seq, CNa_eq, ANa_eq, CCa_eq, ACa_eq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack states m, h, n, s, C_Na, A_Na, C_Ca, A_Ca = states # Standard gating states derivatives dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) # PumpNa current pool concentration and activation state I_Na = self.currNa(m, h, Vm) dCNa_dt = self.derC_Na(C_Na, I_Na) dANa_dt = self.derA_Na(A_Na, C_Na) # KCa current pool concentration and activation state I_Ca = self.currCa(s, Vm) dCCa_dt = self.derC_Ca(C_Ca, I_Ca) dACa_dt = self.derA_Ca(A_Ca, C_Ca) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dCNa_dt, dANa_dt, dCCa_dt, dACa_dt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants Tm = self.taum minf = self.minf(Vm) am_avg = np.mean(minf / Tm) bm_avg = np.mean(1 / Tm) - am_avg Th = self.tauh(Vm) hinf = self.hinf(Vm) ah_avg = np.mean(hinf / Th) bh_avg = np.mean(1 / Th) - ah_avg Tn = self.taun(Vm) ninf = self.ninf(Vm) an_avg = np.mean(ninf / Tn) bn_avg = np.mean(1 / Tn) - an_avg Ts = self.taus sinf = self.sinf(Vm) as_avg = np.mean(sinf / Ts) bs_avg = np.mean(1 / Ts) - as_avg # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) Vmeff = np.interp(Qm, interp_data['Q'], interp_data['V']) # Unpack states m, h, n, s, C_Na, A_Na, C_Ca, A_Ca = states # Standard gating states derivatives dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s # PumpNa current pool concentration and activation state I_Na = self.currNa(m, h, Vmeff) dCNa_dt = self.derC_Na(C_Na, I_Na) dANa_dt = self.derA_Na(A_Na, C_Na) # KCa current pool concentration and activation state I_Ca_eff = self.currCa(s, Vmeff) dCCa_dt = self.derC_Ca(C_Ca, I_Ca_eff) dACa_dt = self.derA_Ca(A_Ca, C_Ca) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dCNa_dt, dANa_dt, dCCa_dt, dACa_dt] class LeechMech(BaseMech): ''' Class defining the basic dynamics of Sodium, Potassium and Calcium channels for several neurons of the leech. Reference: *Baccus, S.A. (1998). Synaptic facilitation by reflected action potentials: enhancement of transmission when nerve impulses reverse direction at axon branch points. Proc. Natl. Acad. Sci. U.S.A. 95, 8345–8350.* ''' alphaC_sf = 1e-5 # Calcium activation rate constant scaling factor (M) betaC = 0.1e3 # beta rate for the open-probability of Ca2+-dependent Potassium channels (s-1) T = 293.15 # Room temperature (K) Rg = 8.314 # Universal gas constant (J.mol^-1.K^-1) Faraday = 9.6485e4 # Faraday constant (C/mol) def alpham(self, Vm): ''' Compute the alpha rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = -0.03 * (Vm + 28) / (np.exp(- (Vm + 28) / 15) - 1) # ms-1 return alpha * 1e3 # s-1 def betam(self, Vm): ''' Compute the beta rate for the open-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 2.7 * np.exp(-(Vm + 53) / 18) # ms-1 return beta * 1e3 # s-1 def alphah(self, Vm): ''' Compute the alpha rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = 0.045 * np.exp(-(Vm + 58) / 18) # ms-1 return alpha * 1e3 # s-1 def betah(self, Vm): ''' Compute the beta rate for the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) .. warning:: the original paper contains an error (multiplication) in the expression of this rate constant, corrected in the mod file on ModelDB (division). ''' beta = 0.72 / (np.exp(-(Vm + 23) / 14) + 1) # ms-1 return beta * 1e3 # s-1 def alphan(self, Vm): ''' Compute the alpha rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = -0.024 * (Vm - 17) / (np.exp(-(Vm - 17) / 8) - 1) # ms-1 return alpha * 1e3 # s-1 def betan(self, Vm): ''' Compute the beta rate for the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 0.2 * np.exp(-(Vm + 48) / 35) # ms-1 return beta * 1e3 # s-1 def alphas(self, Vm): ''' Compute the alpha rate for the open-probability of Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' alpha = -1.5 * (Vm - 20) / (np.exp(-(Vm - 20) / 5) - 1) # ms-1 return alpha * 1e3 # s-1 def betas(self, Vm): ''' Compute the beta rate for the open-probability of Calcium channels. :param Vm: membrane potential (mV) :return: rate constant (s-1) ''' beta = 1.5 * np.exp(-(Vm + 25) / 10) # ms-1 return beta * 1e3 # s-1 def alphaC(self, C_Ca_in): ''' Compute the alpha rate for the open-probability of Calcium-dependent Potassium channels. :param C_Ca_in: intracellular Calcium concentration (M) :return: rate constant (s-1) ''' alpha = 0.1 * C_Ca_in / self.alphaC_sf # ms-1 return alpha * 1e3 # s-1 def derM(self, Vm, m): ''' Compute the evolution of the open-probability of Sodium channels. :param Vm: membrane potential (mV) :param m: open-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alpham(Vm) * (1 - m) - self.betam(Vm) * m def derH(self, Vm, h): ''' Compute the evolution of the inactivation-probability of Sodium channels. :param Vm: membrane potential (mV) :param h: inactivation-probability of Sodium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphah(Vm) * (1 - h) - self.betah(Vm) * h def derN(self, Vm, n): ''' Compute the evolution of the open-probability of delayed-rectifier Potassium channels. :param Vm: membrane potential (mV) :param n: open-probability of delayed-rectifier Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphan(Vm) * (1 - n) - self.betan(Vm) * n def derS(self, Vm, s): ''' Compute the evolution of the open-probability of Calcium channels. :param Vm: membrane potential (mV) :param s: open-probability of Calcium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphas(Vm) * (1 - s) - self.betas(Vm) * s def derC(self, c, C_Ca_in): ''' Compute the evolution of the open-probability of Calcium-dependent Potassium channels. :param c: open-probability of Calcium-dependent Potassium channels (prob) :param C_Ca_in: intracellular Calcium concentration (M) :return: derivative of open-probability w.r.t. time (prob/s) ''' return self.alphaC(C_Ca_in) * (1 - c) - self.betaC * c def currNa(self, m, h, Vm, C_Na_in): ''' Compute the inward Sodium current per unit area. :param m: open-probability of Sodium channels :param h: inactivation-probability of Sodium channels :param Vm: membrane potential (mV) :param C_Na_in: intracellular Sodium concentration (M) :return: current per unit area (mA/m2) ''' GNa = self.GNaMax * m**4 * h VNa = self.nernst(self.Z_Na, C_Na_in, self.C_Na_out) # Sodium Nernst potential return GNa * (Vm - VNa) def currK(self, n, Vm): ''' Compute the outward, delayed-rectifier Potassium current per unit area. :param n: open-probability of delayed-rectifier Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GKMax * n**2 return GK * (Vm - self.VK) def currCa(self, s, Vm, C_Ca_in): ''' Compute the inward Calcium current per unit area. :param s: open-probability of Calcium channels :param Vm: membrane potential (mV) :param C_Ca_in: intracellular Calcium concentration (M) :return: current per unit area (mA/m2) ''' GCa = self.GCaMax * s VCa = self.nernst(self.Z_Ca, C_Ca_in, self.C_Ca_out) # Calcium Nernst potential return GCa * (Vm - VCa) def currKCa(self, c, Vm): ''' Compute the outward Calcium-dependent Potassium current per unit area. :param c: open-probability of Calcium-dependent Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GKCa = self.GKCaMax * c return GKCa * (Vm - self.VK) def currL(self, Vm): ''' Compute the non-specific leakage current per unit area. :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' return self.GL * (Vm - self.VL) class LeechPressure(LeechMech): ''' Class defining the membrane channel dynamics of a leech pressure sensory neuron. with 7 different current types: - Inward Sodium current - Outward Potassium current - Inward high-voltage-activated Calcium current - Non-specific leakage current - Calcium-dependent, outward Potassium current - Sodium pump current - Calcium pump current Reference: *Baccus, S.A. (1998). Synaptic facilitation by reflected action potentials: enhancement of transmission when nerve impulses reverse direction at axon branch points. Proc. Natl. Acad. Sci. U.S.A. 95, 8345–8350.* ''' # Name of channel mechanism name = 'LeechP' # Cell-specific biophysical parameters Cm0 = 1e-2 # Cell membrane resting capacitance (F/m2) Vm0 = -48.865 # Cell membrane resting potential (mV) C_Na_out = 0.11 # Sodium extracellular concentration (M) C_Ca_out = 1.8e-3 # Calcium extracellular concentration (M) C_Na_in0 = 0.01 # Initial Sodium intracellular concentration (M) C_Ca_in0 = 1e-7 # Initial Calcium intracellular concentration (M) # VNa = 60 # Sodium Nernst potential, from MOD file on ModelDB (mV) # VCa = 125 # Calcium Nernst potential, from MOD file on ModelDB (mV) VK = -68.0 # Potassium Nernst potential (mV) VL = -49.0 # Non-specific leakage Nernst potential (mV) INaPmax = 70.0 # Maximum pump rate of the NaK-ATPase (mA/m2) khalf_Na = 0.012 # Sodium concentration at which NaK-ATPase is at half its maximum rate (M) ksteep_Na = 1e-3 # Sensitivity of NaK-ATPase to varying Sodium concentrations (M) iCaS = 0.1 # Calcium pump current parameter (mA/m2) GNaMax = 3500.0 # Max. conductance of Sodium current (S/m^2) GKMax = 60.0 # Max. conductance of Potassium current (S/m^2) GCaMax = 0.02 # Max. conductance of Calcium current (S/m^2) GKCaMax = 8.0 # Max. conductance of Calcium-dependent Potassium current (S/m^2) GL = 5.0 # Conductance of non-specific leakage current (S/m^2) diam = 50e-6 # Cell soma diameter (m) Z_Na = 1 # Sodium valence Z_Ca = 2 # Calcium valence # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h', 'm4h'], 'i_K\ kin.': ['n'], 'i_{Ca}\ kin.': ['s'], 'i_{KCa}\ kin.': ['c'], 'pools': ['C_Na', 'C_Ca'], 'I': ['iNa2', 'iK', 'iCa2', 'iKCa2', 'iPumpNa2', 'iPumpCa2', 'iL', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' SV_ratio = 6 / self.diam # surface to volume ratio of the (spherical) cell soma # Conversion constant from membrane ionic currents into # change rate of intracellular ionic concentrations self.K_Na = SV_ratio / (self.Z_Na * self.Faraday) * 1e-6 # Sodium (M/s) self.K_Ca = SV_ratio / (self.Z_Ca * self.Faraday) * 1e-6 # Calcium (M/s) # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'c', 'C_Na', 'C_Ca'] self.states0 = np.array([]) # Names of the channels effective coefficients self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas'] # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) def nernst(self, z_ion, C_ion_in, C_ion_out): ''' Return the Nernst potential of a specific ion given its intra and extracellular concentrations. :param z_ion: ion valence :param C_ion_in: intracellular ion concentration (M) :param C_ion_out: extracellular ion concentration (M) :return: ion Nernst potential (mV) ''' return (self.Rg * self.T) / (z_ion * self.Faraday) * np.log(C_ion_out / C_ion_in) * 1e3 def currPumpNa(self, C_Na_in): ''' Outward current mimicking the activity of the NaK-ATPase pump. :param C_Na_in: intracellular Sodium concentration (M) :return: current per unit area (mA/m2) ''' return self.INaPmax / (1 + np.exp((self.khalf_Na - C_Na_in) / self.ksteep_Na)) def currPumpCa(self, C_Ca_in): ''' Outward current representing the activity of a Calcium pump. :param C_Ca_in: intracellular Calcium concentration (M) :return: current per unit area (mA/m2) ''' return self.iCaS * (C_Ca_in - self.C_Ca_in0) / 1.5 def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, c, C_Na_in, C_Ca_in = states return (self.currNa(m, h, Vm, C_Na_in) + self.currK(n, Vm) + self.currCa(s, Vm, C_Ca_in) + self.currKCa(c, Vm) + self.currL(Vm) + (self.currPumpNa(C_Na_in) / 3.) + self.currPumpCa(C_Ca_in)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Intracellular concentrations C_Na_eq = self.C_Na_in0 C_Ca_eq = self.C_Ca_in0 # Standard gating dynamics: Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) seq = self.alphas(Vm) / (self.alphas(Vm) + self.betas(Vm)) ceq = self.alphaC(C_Ca_eq) / (self.alphaC(C_Ca_eq) + self.betaC) return np.array([meq, heq, neq, seq, ceq, C_Na_eq, C_Ca_eq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack states m, h, n, s, c, C_Na_in, C_Ca_in = states # Standard gating states derivatives dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) dcdt = self.derC(c, C_Ca_in) # Intracellular concentrations dCNa_dt = - (self.currNa(m, h, Vm, C_Na_in) + self.currPumpNa(C_Na_in)) * self.K_Na # M/s dCCa_dt = -(self.currCa(s, Vm, C_Ca_in) + self.currPumpCa(C_Ca_in)) * self.K_Ca # M/s # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dCNa_dt, dCCa_dt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) as_avg = np.mean(self.alphas(Vm)) bs_avg = np.mean(self.betas(Vm)) # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) Vmeff = np.interp(Qm, interp_data['Q'], interp_data['V']) # Unpack states m, h, n, s, c, C_Na_in, C_Ca_in = states # Standard gating states derivatives dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s # KCa current gating state derivative dcdt = self.derC(c, C_Ca_in) # Intracellular concentrations dCNa_dt = - (self.currNa(m, h, Vmeff, C_Na_in) + self.currPumpNa(C_Na_in)) * self.K_Na # M/s dCCa_dt = -(self.currCa(s, Vmeff, C_Ca_in) + self.currPumpCa(C_Ca_in)) * self.K_Ca # M/s # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dCNa_dt, dCCa_dt] class LeechRetzius(LeechMech): ''' Class defining the membrane channel dynamics of a leech Retzius neuron. with 5 different current types: - Inward Sodium current - Outward Potassium current - Inward high-voltage-activated Calcium current - Non-specific leakage current - Calcium-dependent, outward Potassium current References: *Vazquez, Y., Mendez, B., Trueta, C., and De-Miguel, F.F. (2009). Summation of excitatory postsynaptic potentials in electrically-coupled neurones. Neuroscience 163, 202–212.* *ModelDB link: https://senselab.med.yale.edu/modeldb/ShowModel.cshtml?model=120910* ''' # Name of channel mechanism - name = 'LeechR' + # name = 'LeechR' # Cell-specific biophysical parameters Cm0 = 5e-2 # Cell membrane resting capacitance (F/m2) Vm0 = -44.45 # Cell membrane resting potential (mV) VNa = 50.0 # Sodium Nernst potential, from retztemp.ses file on ModelDB (mV) VCa = 125.0 # Calcium Nernst potential, from cachdend.mod file on ModelDB (mV) VK = -79.0 # Potassium Nernst potential, from retztemp.ses file on ModelDB (mV) VL = -30.0 # Non-specific leakage Nernst potential, from leakdend.mod file on ModelDB (mV) GNaMax = 1250.0 # Max. conductance of Sodium current (S/m^2) GKMax = 10.0 # Max. conductance of Potassium current (S/m^2) GAMax = 100.0 # Max. conductance of transient Potassium current (S/m^2) GCaMax = 4.0 # Max. conductance of Calcium current (S/m^2) GKCaMax = 130.0 # Max. conductance of Calcium-dependent Potassium current (S/m^2) GL = 1.25 # Conductance of non-specific leakage current (S/m^2) Vhalf = -73.1 # mV C_Ca_in = 5e-8 # Calcium intracellular concentration, from retztemp.ses file (M) # Default plotting scheme pltvars_scheme = { 'i_{Na}\ kin.': ['m', 'h', 'm4h'], 'i_K\ kin.': ['n'], 'i_A\ kin.': ['a', 'b', 'ab'], 'i_{Ca}\ kin.': ['s'], 'i_{KCa}\ kin.': ['c'], 'I': ['iNa', 'iK', 'iCa', 'iKCa2', 'iA', 'iL', 'iNet'] } def __init__(self): ''' Constructor of the class. ''' # Names and initial states of the channels state probabilities self.states_names = ['m', 'h', 'n', 's', 'c', 'a', 'b'] self.states0 = np.array([]) # Names of the channels effective coefficients self.coeff_names = ['alpham', 'betam', 'alphah', 'betah', 'alphan', 'betan', 'alphas', 'betas', 'alphac', 'betac', 'alphaa', 'betaa' 'alphab', 'betab'] # Define initial channel probabilities (solving dx/dt = 0 at resting potential) self.states0 = self.steadyStates(self.Vm0) def ainf(self, Vm): ''' Steady-state activation probability of transient Potassium channels. Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' Vth = -55.0 # mV return 0 if Vm <= Vth else min(1, 2 * (Vm - Vth)**3 / ((11 - Vth)**3 + (Vm - Vth)**3)) def taua(self, Vm): ''' Activation time constant of transient Potassium channels. (assuming T = 20°C). Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' x = -1.5 * (Vm - self.Vhalf) * 1e-3 * self.Faraday / (self.Rg * self.T) # [-] alpha = np.exp(x) # ms-1 beta = np.exp(0.7 * x) # ms-1 return max(0.5, beta / (0.3 * (1 + alpha))) * 1e-3 # s def binf(self, Vm): ''' Steady-state inactivation probability of transient Potassium channels. Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' return 1. / (1 + np.exp((self.Vhalf - Vm) / -6.3)) def taub(self, Vm): ''' Inactivation time constant of transient Potassium channels. (assuming T = 20°C). Source: *Beck, H., Ficker, E., and Heinemann, U. (1992). Properties of two voltage-activated potassium currents in acutely isolated juvenile rat dentate gyrus granule cells. J. Neurophysiol. 68, 2086–2099.* :param Vm: membrane potential (mV) :return: time constant (s) ''' x = 2 * (Vm - self.Vhalf) * 1e-3 * self.Faraday / (self.Rg * self.T) alpha = np.exp(x) beta = np.exp(0.65 * x) return max(7.5, beta / (0.02 * (1 + alpha))) * 1e-3 # s def derA(self, Vm, a): ''' Compute the evolution of the activation-probability of transient Potassium channels. :param Vm: membrane potential (mV) :param a: activation-probability of transient Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.ainf(Vm) - a) / self.taua(Vm) def derB(self, Vm, b): ''' Compute the evolution of the inactivation-probability of transient Potassium channels. :param Vm: membrane potential (mV) :param b: inactivation-probability of transient Potassium channels (prob) :return: derivative of open-probability w.r.t. time (prob/s) ''' return (self.binf(Vm) - b) / self.taub(Vm) def currA(self, a, b, Vm): ''' Compute the outward, transient Potassium current per unit area. :param a: open-probability of transient Potassium channels :param b: inactivation-probability of transient Potassium channels :param Vm: membrane potential (mV) :return: current per unit area (mA/m2) ''' GK = self.GAMax * a * b return GK * (Vm - self.VK) def currNet(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' m, h, n, s, c, a, b = states return (self.currNa(m, h, Vm) + self.currK(n, Vm) + self.currCa(s, Vm) + self.currL(Vm) + self.currKCa(c, Vm) + self.currA(a, b, Vm)) # mA/m2 def steadyStates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Standard gating dynamics: Solve the equation dx/dt = 0 at Vm for each x-state meq = self.alpham(Vm) / (self.alpham(Vm) + self.betam(Vm)) heq = self.alphah(Vm) / (self.alphah(Vm) + self.betah(Vm)) neq = self.alphan(Vm) / (self.alphan(Vm) + self.betan(Vm)) seq = self.alphas(Vm) / (self.alphas(Vm) + self.betas(Vm)) ceq = self.alphaC(self.C_Ca_in) / (self.alphaC(self.C_Ca_in) + self.betaC) aeq = self.ainf(Vm) beq = self.binf(Vm) return np.array([meq, heq, neq, seq, ceq, aeq, beq]) def derStates(self, Vm, states): ''' Concrete implementation of the abstract API method. ''' # Unpack states m, h, n, s, c, a, b = states # Standard gating states derivatives dmdt = self.derM(Vm, m) dhdt = self.derH(Vm, h) dndt = self.derN(Vm, n) dsdt = self.derS(Vm, s) dcdt = self.derC(c, self.C_Ca_in) dadt = self.derA(Vm, a) dbdt = self.derB(Vm, b) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dadt, dbdt] def getEffRates(self, Vm): ''' Concrete implementation of the abstract API method. ''' # Compute average cycle value for rate constants am_avg = np.mean(self.alpham(Vm)) bm_avg = np.mean(self.betam(Vm)) ah_avg = np.mean(self.alphah(Vm)) bh_avg = np.mean(self.betah(Vm)) an_avg = np.mean(self.alphan(Vm)) bn_avg = np.mean(self.betan(Vm)) as_avg = np.mean(self.alphas(Vm)) bs_avg = np.mean(self.betas(Vm)) Ta = self.taua(Vm) ainf = self.ainf(Vm) aa_avg = np.mean(ainf / Ta) ba_avg = np.mean(1 / Ta) - aa_avg Tb = self.taub(Vm) binf = self.binf(Vm) ab_avg = np.mean(binf / Tb) bb_avg = np.mean(1 / Tb) - ab_avg # Return array of coefficients return np.array([am_avg, bm_avg, ah_avg, bh_avg, an_avg, bn_avg, as_avg, bs_avg, aa_avg, ba_avg, ab_avg, bb_avg]) def derStatesEff(self, Qm, states, interp_data): ''' Concrete implementation of the abstract API method. ''' rates = np.array([np.interp(Qm, interp_data['Q'], interp_data[rn]) for rn in self.coeff_names]) # Unpack states m, h, n, s, c, a, b = states # Standard gating states derivatives dmdt = rates[0] * (1 - m) - rates[1] * m dhdt = rates[2] * (1 - h) - rates[3] * h dndt = rates[4] * (1 - n) - rates[5] * n dsdt = rates[6] * (1 - s) - rates[7] * s dadt = rates[8] * (1 - a) - rates[9] * a dbdt = rates[10] * (1 - b) - rates[11] * b # KCa current gating state derivative dcdt = self.derC(c, self.C_Ca_in) # Pack derivatives and return return [dmdt, dhdt, dndt, dsdt, dcdt, dadt, dbdt] diff --git a/PointNICE/constants.py b/PointNICE/constants.py index 2be9687..bccc8f7 100644 --- a/PointNICE/constants.py +++ b/PointNICE/constants.py @@ -1,52 +1,52 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-11-04 13:23:31 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-01-26 12:07:41 +# @Last Modified time: 2018-03-07 14:55:04 ''' Algorithmic constants used in the package. ''' # Fitting and pre-processing LJFIT_PM_MAX = 1e8 # intermolecular pressure at the deflection lower bound for LJ fitting (Pa) PNET_EQ_MAX = 1e-1 # error threshold for net pressure at computed equilibrium position (Pa) PMAVG_STD_ERR_MAX = 1500 # error threshold in nonlinear fit of molecular pressure (Pa) # Mechanical simulations Z_ERR_MAX = 1e-11 # periodic convergence threshold for deflection gas content (m) NG_ERR_MAX = 1e-24 # periodic convergence threshold for gas content (mol) INPUT_CHARGE_RANGE = (-120e-5, 60e-5) # physiological charge range constraining the membrane (C/m2) # E-STIM simulations DT_ESTIM = 1e-4 # A-STIM simulations SOLVER_NSTEPS = 1000 # maximum number of steps allowed during one call to the LSODA/DOP853 solvers CLASSIC_DS_FACTOR = 3 # time downsampling factor applied to output arrays of classic simulations NPC_FULL = 1000 # nb of samples per acoustic period in full system NPC_HH = 40 # nb of samples per acoustic period in HH system DQ_UPDATE = 1e-5 # charge evolution threshold between two hybrid integrations (C/m2) DT_UPDATE = 5e-4 # time interval between two hybrid integrations (s) DT_EFF = 5e-5 # time step for effective integration (s) # Spike detection SPIKE_MIN_QAMP = 10e-5 # threshold amplitude for spike detection on charge signal (C/m2) SPIKE_MIN_VAMP = 10.0 # threshold amplitude for spike detection on potential signal (mV) SPIKE_MIN_DT = 1e-3 # minimal time interval for spike detection on charge signal (s) # Titrations TITRATION_T_OFFSET = 50e-3 # offset period for titration procedures (s) TITRATION_ASTIM_A_MAX = 3e5 # initial acoustic pressure upper bound for titration (Pa) TITRATION_ASTIM_DA_MAX = 1e3 # acoustic pressure search range threshold for titration (Pa) TITRATION_ESTIM_A_MAX = 50.0 # initial current density upper bound for titration (mA/m2) TITRATION_ESTIM_DA_MAX = 0.1 # current density search range threshold for titration (mA/m2) TITRATION_T_MAX = 2e-1 # initial stimulus duration upper bound for titration (s) TITRATION_DT_THR = 1e-3 # stimulus duration search range threshold for titration (s) TITRATION_DDF_THR = 0.01 # stimulus duty cycle search range threshold for titration (-) TITRATION_DF_MAX = 1.0 # initial stimulus duty cycle upper bound for titration (-) diff --git a/PointNICE/lookups/BLS_lookups_a100.0nm.json b/PointNICE/lookups/BLS_lookups_a100.0nm.json new file mode 100644 index 0000000..a6f83e2 --- /dev/null +++ b/PointNICE/lookups/BLS_lookups_a100.0nm.json @@ -0,0 +1 @@ +{"-80.00": {"LJ_approx": {"x0": 1.7817059603461257e-09, "C": 14691.28061115401, "nrep": 3.9113369016354067, "nattr": 0.9360192402929897}, "Delta_eq": 1.2344323203763867e-09}} \ No newline at end of file diff --git a/PointNICE/lookups/BLS_lookups_a20.0nm.json b/PointNICE/lookups/BLS_lookups_a20.0nm.json new file mode 100644 index 0000000..ed30ab3 --- /dev/null +++ b/PointNICE/lookups/BLS_lookups_a20.0nm.json @@ -0,0 +1 @@ +{"-80.00": {"LJ_approx": {"x0": 1.7921352427645028e-09, "C": 14362.363951301933, "nrep": 3.911239140144219, "nattr": 0.9577057656995782}, "Delta_eq": 1.2344323203763867e-09}} \ No newline at end of file diff --git a/PointNICE/lookups/BLS_lookups_a50.0nm.json b/PointNICE/lookups/BLS_lookups_a50.0nm.json new file mode 100644 index 0000000..7337d0d --- /dev/null +++ b/PointNICE/lookups/BLS_lookups_a50.0nm.json @@ -0,0 +1 @@ +{"-80.00": {"LJ_approx": {"x0": 1.7845793287878114e-09, "C": 14600.797866604744, "nrep": 3.9112835250651723, "nattr": 0.9433189447961238}, "Delta_eq": 1.2344323203763867e-09}} \ No newline at end of file diff --git a/PointNICE/lookups/BLS_lookups_a500.0nm.json b/PointNICE/lookups/BLS_lookups_a500.0nm.json new file mode 100644 index 0000000..08d6be3 --- /dev/null +++ b/PointNICE/lookups/BLS_lookups_a500.0nm.json @@ -0,0 +1 @@ +{"-80.00": {"LJ_approx": {"x0": 1.7789681895702171e-09, "C": 14776.517324235609, "nrep": 3.9114315801143102, "nattr": 0.9262993775092675}, "Delta_eq": 1.2344323203763867e-09}} \ No newline at end of file diff --git a/PointNICE/solvers/simutils.py b/PointNICE/solvers/simutils.py index b37744b..9a5aa14 100644 --- a/PointNICE/solvers/simutils.py +++ b/PointNICE/solvers/simutils.py @@ -1,1210 +1,1219 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-22 14:33:04 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-01-26 12:07:01 +# @Last Modified time: 2018-03-07 15:22:18 """ Utility functions used in simulations """ import os import time import logging import pickle import shutil import tkinter as tk from tkinter import filedialog import numpy as np from openpyxl import load_workbook import lockfile from ..bls import BilayerSonophore from .SolverUS import SolverUS from .SolverElec import SolverElec from ..constants import * from ..utils import getNeuronsDict # Get package logger logger = logging.getLogger('PointNICE') # Naming nomenclature for output files MECH_code = 'MECH_{:.0f}nm_{:.0f}kHz_{:.0f}kPa_{:.1f}nCcm2' ESTIM_CW_code = 'ESTIM_{}_CW_{:.1f}mA_per_m2_{:.0f}ms' ESTIM_PW_code = 'ESTIM_{}_PW_{:.1f}mA_per_m2_{:.0f}ms_PRF{:.2f}kHz_DF{:.2f}' ASTIM_CW_code = 'ASTIM_{}_CW_{:.0f}nm_{:.0f}kHz_{:.0f}kPa_{:.0f}ms_{}' ASTIM_PW_code = 'ASTIM_{}_PW_{:.0f}nm_{:.0f}kHz_{:.0f}kPa_{:.0f}ms_PRF{:.2f}kHz_DF{:.2f}_{}' # Parameters units ASTIM_params = { 'f': {'index': 0, 'factor': 1e-3, 'unit': 'kHz'}, 'A': {'index': 1, 'factor': 1e-3, 'unit': 'kPa'}, 't': {'index': 2, 'factor': 1e3, 'unit': 'ms'}, 'PRF': {'index': 4, 'factor': 1e-3, 'unit': 'kHz'}, 'DF': {'index': 5, 'factor': 1e2, 'unit': '%'} } ESTIM_params = { 'A': {'index': 0, 'factor': 1e0, 'unit': 'mA/m2'}, 't': {'index': 1, 'factor': 1e3, 'unit': 'ms'}, 'PRF': {'index': 3, 'factor': 1e-3, 'unit': 'kHz'}, 'DF': {'index': 4, 'factor': 1e2, 'unit': '%'} } # Default geometry default_diam = 32e-9 default_embedding = 0.0e-6 def setBatchDir(): ''' Select batch directory for output files.α :return: full path to batch directory ''' root = tk.Tk() root.withdraw() batch_dir = filedialog.askdirectory() assert batch_dir, 'No batch directory chosen' return batch_dir def checkBatchLog(batch_dir, batch_type): ''' Check for appropriate log file in batch directory, and create one if it is absent. :param batch_dir: full path to batch directory :param batch_type: type of simulation batch :return: 2 tuple with full path to log file and boolean stating if log file was created ''' # Determine log template from batch type if batch_type == 'MECH': logfile = 'log_MECH.xlsx' elif batch_type == 'A-STIM': logfile = 'log_ASTIM.xlsx' elif batch_type == 'E-STIM': logfile = 'log_ESTIM.xlsx' else: raise ValueError('Unknown batch type', batch_type) # Get template in package subdirectory this_dir, _ = os.path.split(__file__) parent_dir = os.path.abspath(os.path.join(this_dir, os.pardir)) logsrc = parent_dir + '/templates/' + logfile assert os.path.isfile(logsrc), 'template log file "{}" not found'.format(logsrc) # Copy template in batch directory if no appropriate log file logdst = batch_dir + '/' + logfile is_log = os.path.isfile(logdst) if not is_log: shutil.copy2(logsrc, logdst) return (logdst, not is_log) def createSimQueue(amps, durations, offsets, PRFs, DFs): ''' Create a serialized 2D array of all parameter combinations for a series of individual parameter sweeps, while avoiding repetition of CW protocols for a given PRF sweep. :param amps: list (or 1D-array) of acoustic amplitudes :param durations: list (or 1D-array) of stimulus durations :param offsets: list (or 1D-array) of stimulus offsets (paired with durations array) :param PRFs: list (or 1D-array) of pulse-repetition frequencies :param DFs: list (or 1D-array) of duty cycle values :return: 2D-array with (amplitude, duration, offset, PRF, DF) for each stimulation protocol ''' # Convert input to 1D-arrays amps = np.array(amps) durations = np.array(durations) offsets = np.array(offsets) PRFs = np.array(PRFs) DFs = np.array(DFs) # Create index arrays iamps = range(len(amps)) idurs = range(len(durations)) # Create empty output matrix queue = np.empty((1, 5)) # Continuous protocols if 1.0 in DFs: nCW = len(amps) * len(durations) arr1 = np.ones(nCW) iCW_queue = np.array(np.meshgrid(iamps, idurs)).T.reshape(nCW, 2) CW_queue = np.vstack((amps[iCW_queue[:, 0]], durations[iCW_queue[:, 1]], offsets[iCW_queue[:, 1]], PRFs.min() * arr1, arr1)).T queue = np.vstack((queue, CW_queue)) # Pulsed protocols if np.any(DFs != 1.0): pulsed_DFs = DFs[DFs != 1.0] iPRFs = range(len(PRFs)) ipulsed_DFs = range(len(pulsed_DFs)) nPW = len(amps) * len(durations) * len(PRFs) * len(pulsed_DFs) iPW_queue = np.array(np.meshgrid(iamps, idurs, iPRFs, ipulsed_DFs)).T.reshape(nPW, 4) PW_queue = np.vstack((amps[iPW_queue[:, 0]], durations[iPW_queue[:, 1]], offsets[iPW_queue[:, 1]], PRFs[iPW_queue[:, 2]], pulsed_DFs[iPW_queue[:, 3]])).T queue = np.vstack((queue, PW_queue)) # Return return queue[1:, :] def xlslog(filename, sheetname, data): """ Append log data on a new row to specific sheet of excel workbook, using a lockfile to avoid read/write errors between concurrent processes. :param filename: absolute or relative path to the Excel workbook :param sheetname: name of the Excel spreadsheet to which data is appended :param data: data structure to be added to specific columns on a new row :return: boolean indicating success (1) or failure (0) of operation """ try: lock = lockfile.FileLock(filename) lock.acquire() wb = load_workbook(filename) ws = wb[sheetname] keys = data.keys() i = 1 row_data = {} for k in keys: row_data[k] = data[k] i += 1 ws.append(row_data) wb.save(filename) lock.release() return 1 except PermissionError: # If file cannot be accessed for writing because already opened logger.warning('Cannot write to "%s". Close the file and type "Y"', filename) user_str = input() if user_str in ['y', 'Y']: return xlslog(filename, sheetname, data) else: return 0 def detectPeaks(x, mph=None, mpd=1, threshold=0, edge='rising', kpsh=False, valley=False, ax=None): """ Detect peaks in data based on their amplitude and inter-peak distance. """ x = np.atleast_1d(x).astype('float64') if x.size < 3: return np.array([], dtype=int) if valley: x = -x # find indices of all peaks dx = x[1:] - x[:-1] # handle NaN's indnan = np.where(np.isnan(x))[0] if indnan.size: x[indnan] = np.inf dx[np.where(np.isnan(dx))[0]] = np.inf ine, ire, ife = np.array([[], [], []], dtype=int) if not edge: ine = np.where((np.hstack((dx, 0)) < 0) & (np.hstack((0, dx)) > 0))[0] else: if edge.lower() in ['rising', 'both']: ire = np.where((np.hstack((dx, 0)) <= 0) & (np.hstack((0, dx)) > 0))[0] if edge.lower() in ['falling', 'both']: ife = np.where((np.hstack((dx, 0)) < 0) & (np.hstack((0, dx)) >= 0))[0] ind = np.unique(np.hstack((ine, ire, ife))) # handle NaN's if ind.size and indnan.size: # NaN's and values close to NaN's cannot be peaks ind = ind[np.in1d(ind, np.unique(np.hstack((indnan, indnan - 1, indnan + 1))), invert=True)] # first and last values of x cannot be peaks if ind.size and ind[0] == 0: ind = ind[1:] if ind.size and ind[-1] == x.size - 1: ind = ind[:-1] # remove peaks < minimum peak height if ind.size and mph is not None: ind = ind[x[ind] >= mph] # remove peaks - neighbors < threshold if ind.size and threshold > 0: dx = np.min(np.vstack([x[ind] - x[ind - 1], x[ind] - x[ind + 1]]), axis=0) ind = np.delete(ind, np.where(dx < threshold)[0]) # detect small peaks closer than minimum peak distance if ind.size and mpd > 1: ind = ind[np.argsort(x[ind])][::-1] # sort ind by peak height idel = np.zeros(ind.size, dtype=bool) for i in range(ind.size): if not idel[i]: # keep peaks with the same height if kpsh is True idel = idel | (ind >= ind[i] - mpd) & (ind <= ind[i] + mpd) \ & (x[ind[i]] > x[ind] if kpsh else True) idel[i] = 0 # Keep current peak # remove the small peaks and sort back the indices by their occurrence ind = np.sort(ind[~idel]) return ind def detectPeaksTime(t, y, mph, mtd): """ Extension of the detectPeaks function to detect peaks in data based on their amplitude and time difference, with a non-uniform time vector. :param t: time vector (not necessarily uniform) :param y: signal :param mph: minimal peak height :param mtd: minimal time difference :return: array of peak indexes """ # Detect peaks on signal with no restriction on inter-peak distance raw_indexes = detectPeaks(y, mph, mpd=1) if raw_indexes.size > 0: # Filter relevant peaks with temporal distance n_raw = raw_indexes.size filtered_indexes = np.array([raw_indexes[0]]) for i in range(1, n_raw): i1 = filtered_indexes[-1] i2 = raw_indexes[i] if t[i2] - t[i1] < mtd: if y[i2] > y[i1]: filtered_indexes[-1] = i2 else: filtered_indexes = np.append(filtered_indexes, i2) # Return peak indexes return filtered_indexes else: return None def detectSpikes(t, Qm, min_amp, min_dt): ''' Detect spikes on a charge density signal, and return their number, latency and rate. :param t: time vector (s) :param Qm: charge density vector (C/m2) :param min_amp: minimal charge amplitude to detect spikes (C/m2) :param min_dt: minimal time interval between 2 spikes (s) :return: 3-tuple with number of spikes, latency (s) and spike rate (sp/s) ''' i_spikes = detectPeaksTime(t, Qm, min_amp, min_dt) if i_spikes is not None: latency = t[i_spikes[0]] # s n_spikes = i_spikes.size if n_spikes > 1: first_to_last_spike = t[i_spikes[-1]] - t[i_spikes[0]] # s spike_rate = (n_spikes - 1) / first_to_last_spike # spikes/s else: spike_rate = 'N/A' else: latency = 'N/A' spike_rate = 'N/A' n_spikes = 0 return (n_spikes, latency, spike_rate) def runEStim(batch_dir, log_filepath, solver, neuron, Astim, tstim, toffset, PRF, DF): ''' Run a single E-STIM simulation a given neuron for specific stimulation parameters, and save the results in a PKL file. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param solver: SolverElec instance :param Astim: pulse amplitude (mA/m2) :param tstim: pulse duration (s) :param toffset: offset duration (s) :param PRF: pulse repetition frequency (Hz) :param DF: pulse duty factor (-) :return: full path to the output file ''' if DF == 1.0: simcode = ESTIM_CW_code.format(neuron.name, Astim, tstim * 1e3) else: simcode = ESTIM_PW_code.format(neuron.name, Astim, tstim * 1e3, PRF * 1e-3, DF) # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") # Run simulation tstart = time.time() (t, y, states) = solver.run(neuron, Astim, tstim, toffset, PRF, DF) Vm, *channels = y tcomp = time.time() - tstart logger.debug('completed in %.2f seconds', tcomp) # Store data in dictionary data = { 'Astim': Astim, 'tstim': tstim, 'toffset': toffset, 'PRF': PRF, 'DF': DF, 't': t, 'states': states, 'Vm': Vm } for j in range(len(neuron.states_names)): data[neuron.states_names[j]] = channels[j] # Export data to PKL file output_filepath = batch_dir + '/' + simcode + ".pkl" with open(output_filepath, 'wb') as fh: pickle.dump(data, fh) logger.debug('simulation data exported to "%s"', output_filepath) # Detect spikes on Vm signal n_spikes, lat, sr = detectSpikes(t, Vm, SPIKE_MIN_VAMP, SPIKE_MIN_DT) logger.debug('%u spike%s detected', n_spikes, "s" if n_spikes > 1 else "") # Export key metrics to log file log = { 'A': date_str, 'B': daytime_str, 'C': neuron.name, 'D': Astim, 'E': tstim * 1e3, 'F': PRF * 1e-3 if DF < 1 else 'N/A', 'G': DF, 'H': t.size, 'I': round(tcomp, 2), 'J': n_spikes, 'K': lat * 1e3 if isinstance(lat, float) else 'N/A', 'L': sr * 1e-3 if isinstance(sr, float) else 'N/A' } if xlslog(log_filepath, 'Data', log) == 1: logger.debug('log exported to "%s"', log_filepath) else: logger.error('log export to "%s" aborted', log_filepath) return output_filepath def runEStimBatch(batch_dir, log_filepath, neurons, stim_params): ''' Run batch E-STIM simulations of the system for various neuron types and stimulation parameters. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param neurons: list of neurons names :param stim_params: dictionary containing sweeps for all stimulation parameters :return: list of full paths to the output files ''' mandatory_params = ['amps', 'durations', 'offsets', 'PRFs', 'DFs'] for mp in mandatory_params: assert mp in stim_params, 'stim_params dictionary must contain "{}" field'.format(mp) # Define logging format ESTIM_CW_log = 'E-STIM simulation %u/%u: %s neuron, A = %.1f mA/m2, t = %.1f ms' ESTIM_PW_log = ('E-STIM simulation %u/%u: %s neuron, A = %.1f mA/m2, t = %.1f ms, ' 'PRF = %.2f kHz, DF = %.2f') logger.info("Starting E-STIM simulation batch") # Generate simulations queue sim_queue = createSimQueue(stim_params['amps'], stim_params['durations'], stim_params['offsets'], stim_params['PRFs'], stim_params['DFs']) nqueue = sim_queue.shape[0] # Initialize solver solver = SolverElec() # Run simulations nsims = len(neurons) * nqueue simcount = 0 filepaths = [] for nname in neurons: neuron = getNeuronsDict()[nname]() for i in range(nqueue): simcount += 1 Astim, tstim, toffset, PRF, DF = sim_queue[i, :] if DF == 1.0: logger.info(ESTIM_CW_log, simcount, nsims, neuron.name, Astim, tstim * 1e3) else: logger.info(ESTIM_PW_log, simcount, nsims, neuron.name, Astim, tstim * 1e3, PRF * 1e-3, DF) output_filepath = runEStim(batch_dir, log_filepath, solver, neuron, Astim, tstim, toffset, PRF, DF) filepaths.append(output_filepath) return filepaths def runAStim(batch_dir, log_filepath, solver, neuron, Fdrive, Adrive, tstim, toffset, PRF, DF, int_method='effective'): ''' Run a single A-STIM simulation a given neuron for specific stimulation parameters, and save the results in a PKL file. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param solver: SolverUS instance :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 int_method: selected integration method :return: full path to the output file ''' if DF == 1.0: simcode = ASTIM_CW_code.format(neuron.name, solver.a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3, int_method) else: simcode = ASTIM_PW_code.format(neuron.name, solver.a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3, PRF * 1e-3, DF, int_method) # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") # Run simulation tstart = time.time() (t, y, states) = solver.run(neuron, Fdrive, Adrive, tstim, toffset, PRF, DF, int_method) Z, ng, Qm, *channels = y U = np.insert(np.diff(Z) / np.diff(t), 0, 0.0) tcomp = time.time() - tstart logger.debug('completed in %.2f seconds', tcomp) # Store data in dictionary data = { 'a': solver.a, 'd': solver.d, 'Fdrive': Fdrive, 'Adrive': Adrive, 'phi': np.pi, 'tstim': tstim, 'toffset': toffset, 'PRF': PRF, 'DF': DF, 't': t, 'states': states, 'U': U, 'Z': Z, 'ng': ng, 'Qm': Qm, 'Vm': Qm * 1e3 / np.array([solver.Capct(ZZ) for ZZ in Z]) } for j in range(len(neuron.states_names)): data[neuron.states_names[j]] = channels[j] # Export data to PKL file output_filepath = batch_dir + '/' + simcode + ".pkl" with open(output_filepath, 'wb') as fh: pickle.dump(data, fh) logger.debug('simulation data exported to "%s"', output_filepath) # Detect spikes on Qm signal n_spikes, lat, sr = detectSpikes(t, Qm, SPIKE_MIN_QAMP, SPIKE_MIN_DT) logger.debug('%u spike%s detected', n_spikes, "s" if n_spikes > 1 else "") # Export key metrics to log file log = { 'A': date_str, 'B': daytime_str, 'C': neuron.name, 'D': solver.a * 1e9, 'E': solver.d * 1e6, 'F': Fdrive * 1e-3, 'G': Adrive * 1e-3, 'H': tstim * 1e3, 'I': PRF * 1e-3 if DF < 1 else 'N/A', 'J': DF, 'K': int_method, 'L': t.size, 'M': round(tcomp, 2), 'N': n_spikes, 'O': lat * 1e3 if isinstance(lat, float) else 'N/A', 'P': sr * 1e-3 if isinstance(sr, float) else 'N/A' } if xlslog(log_filepath, 'Data', log) == 1: logger.debug('log exported to "%s"', log_filepath) else: logger.error('log export to "%s" aborted', log_filepath) return output_filepath def runAStimBatch(batch_dir, log_filepath, neurons, stim_params, a=default_diam, int_method='effective'): ''' Run batch simulations of the system for various neuron types, sonophore and stimulation parameters. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param neurons: list of neurons names :param stim_params: dictionary containing sweeps for all stimulation parameters :param a: BLS structure diameter (m) :param int_method: selected integration method :return: list of full paths to the output files ''' mandatory_params = ['freqs', 'amps', 'durations', 'offsets', 'PRFs', 'DFs'] for mp in mandatory_params: assert mp in stim_params, 'stim_params dictionary must contain "{}" field'.format(mp) # Define logging format ASTIM_CW_log = ('A-STIM %s simulation %u/%u: %s neuron, a = %.1f nm, f = %.2f kHz, ' 'A = %.2f kPa, t = %.2f ms') ASTIM_PW_log = ('A-STIM %s simulation %u/%u: %s neuron, a = %.1f nm, f = %.2f kHz, ' 'A = %.2f kPa, t = %.2f ms, PRF = %.2f kHz, DF = %.3f') logger.info("Starting A-STIM simulation batch") # Generate simulations queue sim_queue = createSimQueue(stim_params['amps'], stim_params['durations'], stim_params['offsets'], stim_params['PRFs'], stim_params['DFs']) nqueue = sim_queue.shape[0] # Run simulations nsims = len(neurons) * len(stim_params['freqs']) * nqueue simcount = 0 filepaths = [] for nname in neurons: neuron = getNeuronsDict()[nname]() for Fdrive in stim_params['freqs']: # Initialize SolverUS solver = SolverUS(a, neuron, Fdrive) for i in range(nqueue): simcount += 1 Adrive, tstim, toffset, PRF, DF = sim_queue[i, :] # Log and define naming if DF == 1.0: logger.info(ASTIM_CW_log, int_method, simcount, nsims, neuron.name, a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3) else: logger.info(ASTIM_PW_log, int_method, simcount, nsims, neuron.name, a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3, PRF * 1e-3, DF) # Run simulation output_filepath = runAStim(batch_dir, log_filepath, solver, neuron, Fdrive, Adrive, tstim, toffset, PRF, DF, int_method) filepaths.append(output_filepath) return filepaths def titrateEStim(solver, ch_mech, Astim, tstim, toffset, PRF=1.5e3, DF=1.0): """ Use a dichotomic recursive search to determine the threshold value of a specific electric stimulation parameter needed to obtain neural excitation, keeping all other parameters fixed. The titration parameter can be stimulation amplitude, duration or any variable for which the number of spikes is a monotonically increasing function. This function is called recursively until an accurate threshold is found. :param solver: solver instance :param ch_mech: channels mechanism object :param Astim: injected current density amplitude (mA/m2) :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 (-) :return: 5-tuple with the determined amplitude threshold, time profile, solution matrix, state vector and response latency """ # Determine titration type if isinstance(Astim, tuple): t_type = 'A' interval = Astim thr = TITRATION_ESTIM_DA_MAX + maxval = TITRATION_ESTIM_A_MAX elif isinstance(tstim, tuple): t_type = 't' interval = tstim thr = TITRATION_DT_THR + maxval = TITRATION_T_MAX elif isinstance(DF, tuple): t_type = 'DF' interval = DF thr = TITRATION_DDF_THR + maxval = TITRATION_DF_MAX else: logger.error('Invalid titration type') return 0. t_var = ESTIM_params[t_type] # Check amplitude interval and define current value assert interval[0] < interval[1], '{} interval must be defined as (lb, ub)'.format(t_type) value = (interval[0] + interval[1]) / 2 # Define stimulation parameters if t_type == 'A': stim_params = [value, tstim, toffset, PRF, DF] elif t_type == 't': stim_params = [Astim, value, toffset, PRF, DF] elif t_type == 'DF': stim_params = [Astim, tstim, toffset, PRF, value] # Run simulation and detect spikes (t, y, states) = solver.run(ch_mech, *stim_params) n_spikes, latency, _ = detectSpikes(t, y[0, :], SPIKE_MIN_VAMP, SPIKE_MIN_DT) logger.debug('%.2f %s ---> %u spike%s detected', value * t_var['factor'], t_var['unit'], n_spikes, "s" if n_spikes > 1 else "") # If accurate threshold is found, return simulation results if (interval[1] - interval[0]) <= thr and n_spikes == 1: return (value, t, y, states, latency) # Otherwise, refine titration interval and iterate recursively else: if n_spikes == 0: + if (maxval - interval[1]) <= thr: # if upper bound too close to max then stop + logger.warning('no spikes detected within titration interval') + return (np.nan, t, y, states, latency) new_interval = (value, interval[1]) else: new_interval = (interval[0], value) stim_params[t_var['index']] = new_interval return titrateEStim(solver, ch_mech, *stim_params) def titrateAStim(solver, ch_mech, Fdrive, Adrive, tstim, toffset, PRF=1.5e3, DF=1.0, int_method='effective'): """ Use a dichotomic recursive search to determine the threshold value of a specific acoustic stimulation parameter needed to obtain neural excitation, keeping all other parameters fixed. The titration parameter can be stimulation amplitude, duration or any variable for which the number of spikes is a monotonically increasing function. This function is called recursively until an accurate threshold is found. :param solver: solver instance :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 int_method: selected integration method :return: 5-tuple with the determined amplitude threshold, time profile, solution matrix, state vector and response latency """ # Determine titration type if isinstance(Adrive, tuple): t_type = 'A' interval = Adrive thr = TITRATION_ASTIM_DA_MAX + maxval = TITRATION_ASTIM_A_MAX elif isinstance(tstim, tuple): t_type = 't' interval = tstim thr = TITRATION_DT_THR + maxval = TITRATION_T_MAX elif isinstance(DF, tuple): t_type = 'DF' interval = DF thr = TITRATION_DDF_THR + maxval = TITRATION_DF_MAX else: logger.error('Invalid titration type') return 0. t_var = ASTIM_params[t_type] # Check amplitude interval and define current value assert interval[0] < interval[1], '{} interval must be defined as (lb, ub)'.format(t_type) value = (interval[0] + interval[1]) / 2 # Define stimulation parameters if t_type == 'A': stim_params = [Fdrive, value, tstim, toffset, PRF, DF] elif t_type == 't': stim_params = [Fdrive, Adrive, value, toffset, PRF, DF] elif t_type == 'DF': stim_params = [Fdrive, Adrive, tstim, toffset, PRF, value] # Run simulation and detect spikes (t, y, states) = solver.run(ch_mech, *stim_params, int_method) n_spikes, latency, _ = detectSpikes(t, y[2, :], SPIKE_MIN_QAMP, SPIKE_MIN_DT) logger.debug('%.2f %s ---> %u spike%s detected', value * t_var['factor'], t_var['unit'], n_spikes, "s" if n_spikes > 1 else "") # If accurate threshold is found, return simulation results - if (interval[1] - interval[0]) <= thr: - if n_spikes == 0: - value = np.nan - logger.warning('no spikes detected within titration interval') + if (interval[1] - interval[0]) <= thr and n_spikes == 1: return (value, t, y, states, latency) # Otherwise, refine titration interval and iterate recursively else: if n_spikes == 0: + if (maxval - interval[1]) <= thr: # if upper bound too close to max then stop + logger.warning('no spikes detected within titration interval') + return (np.nan, t, y, states, latency) new_interval = (value, interval[1]) else: new_interval = (interval[0], value) stim_params[t_var['index']] = new_interval return titrateAStim(solver, ch_mech, *stim_params, int_method) def titrateAStimBatch(batch_dir, log_filepath, neurons, stim_params, a=default_diam, int_method='effective'): ''' Run batch acoustic titrations of the system for various neuron types, sonophore and stimulation parameters, to determine the threshold of a specific stimulus parameter for neural excitation. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param neurons: list of neurons names :param stim_params: dictionary containing sweeps for all stimulation parameters :param a: BLS structure diameter (m) :param int_method: selected integration method :return: list of full paths to the output files ''' # Define logging format ASTIM_titration_log = '%s neuron - A-STIM titration %u/%u (a = %.1f nm, %s)' logger.info("Starting A-STIM titration batch") # Define default parameters int_method = 'effective' # Determine titration parameter and titrations list if 'durations' not in stim_params: t_type = 't' sim_queue = createSimQueue(stim_params['amps'], [None], [TITRATION_T_OFFSET], stim_params['PRFs'], stim_params['DFs']) # sim_queue = np.delete(sim_queue, 1, axis=1) elif 'amps' not in stim_params: t_type = 'A' sim_queue = createSimQueue([None], stim_params['durations'], [TITRATION_T_OFFSET] * len(stim_params['durations']), stim_params['PRFs'], stim_params['DFs']) elif 'DF' not in stim_params: t_type = 'DF' sim_queue = createSimQueue(stim_params['amps'], stim_params['durations'], [TITRATION_T_OFFSET] * len(stim_params['durations']), stim_params['PRFs'], [None]) nqueue = sim_queue.shape[0] t_var = ASTIM_params[t_type] # Run titrations nsims = len(neurons) * len(stim_params['freqs']) * nqueue simcount = 0 filepaths = [] for nname in neurons: neuron = getNeuronsDict()[nname]() for Fdrive in stim_params['freqs']: try: # Create SolverUS instance (modulus of embedding tissue depends on frequency) solver = SolverUS(a, neuron, Fdrive) for i in range(nqueue): simcount += 1 # Extract parameters Adrive, tstim, toffset, PRF, DF = sim_queue[i, :] if Adrive is None: Adrive = (0., 2 * TITRATION_ASTIM_A_MAX) elif tstim is None: tstim = (0., 2 * TITRATION_T_MAX) elif DF is None: DF = (0., 2 * TITRATION_DF_MAX) curr_params = [Fdrive, Adrive, tstim, PRF, DF] # Generate log str log_str = '' pnames = list(ASTIM_params.keys()) j = 0 for cp in curr_params: pn = pnames[j] pi = ASTIM_params[pn] if not isinstance(cp, tuple): if log_str: log_str += ', ' log_str += '{} = {:.2f} {}'.format(pn, pi['factor'] * cp, pi['unit']) j += 1 # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") # Log logger.info(ASTIM_titration_log, neuron.name, simcount, nsims, a * 1e9, log_str) # Run titration tstart = time.time() (output_thr, t, y, states, lat) = titrateAStim(solver, neuron, Fdrive, Adrive, tstim, toffset, PRF, DF) Z, ng, Qm, *channels = y U = np.insert(np.diff(Z) / np.diff(t), 0, 0.0) tcomp = time.time() - tstart logger.info('completed in %.2f s, threshold = %.2f %s', tcomp, output_thr * t_var['factor'], t_var['unit']) # Determine output variable if t_type == 'A': Adrive = output_thr elif t_type == 't': tstim = output_thr elif t_type == 'DF': DF = output_thr # Define output naming if DF == 1.0: simcode = ASTIM_CW_code.format(neuron.name, a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3, int_method) else: simcode = ASTIM_PW_code.format(neuron.name, a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3, PRF * 1e-3, DF, int_method) # Store data in dictionary data = { 'a': solver.a, 'd': solver.d, 'Fdrive': Fdrive, 'Adrive': Adrive, 'phi': np.pi, 'tstim': tstim, 'toffset': toffset, 'PRF': PRF, 'DF': DF, 't': t, 'states': states, 'U': U, 'Z': Z, 'ng': ng, 'Qm': Qm, 'Vm': Qm * 1e3 / np.array([solver.Capct(ZZ) for ZZ in Z]) } for j in range(len(neuron.states_names)): data[neuron.states_names[j]] = channels[j] # Export data to PKL file output_filepath = batch_dir + '/' + simcode + ".pkl" with open(output_filepath, 'wb') as fh: pickle.dump(data, fh) logger.info('simulation data exported to "%s"', output_filepath) filepaths.append(output_filepath) # Detect spikes on Qm signal n_spikes, lat, sr = detectSpikes(t, Qm, SPIKE_MIN_QAMP, SPIKE_MIN_DT) logger.info('%u spike%s detected', n_spikes, "s" if n_spikes > 1 else "") # Export key metrics to log file log = { 'A': date_str, 'B': daytime_str, 'C': neuron.name, 'D': solver.a * 1e9, 'E': solver.d * 1e6, 'F': Fdrive * 1e-3, 'G': Adrive * 1e-3, 'H': tstim * 1e3, 'I': PRF * 1e-3 if DF < 1 else 'N/A', 'J': DF, 'K': int_method, 'L': t.size, 'M': round(tcomp, 2), 'N': n_spikes, 'O': lat * 1e3 if isinstance(lat, float) else 'N/A', 'P': sr * 1e-3 if isinstance(sr, float) else 'N/A' } if xlslog(log_filepath, 'Data', log) == 1: logger.info('log exported to "%s"', log_filepath) else: logger.error('log export to "%s" aborted', log_filepath) except AssertionError as err: logger.error(err) return filepaths def titrateEStimBatch(batch_dir, log_filepath, neurons, stim_params): ''' Run batch electrical titrations of the system for various neuron types and stimulation parameters, to determine the threshold of a specific stimulus parameter for neural excitation. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param neurons: list of neurons names :param stim_params: dictionary containing sweeps for all stimulation parameters :return: list of full paths to the output files ''' # Define logging format ESTIM_titration_log = '%s neuron - E-STIM titration %u/%u (%s)' logger.info("Starting E-STIM titration batch") # Determine titration parameter and titrations list if 'durations' not in stim_params: t_type = 't' sim_queue = createSimQueue(stim_params['amps'], [None], [TITRATION_T_OFFSET], stim_params['PRFs'], stim_params['DFs']) elif 'amps' not in stim_params: t_type = 'A' sim_queue = createSimQueue([None], stim_params['durations'], [TITRATION_T_OFFSET] * len(stim_params['durations']), stim_params['PRFs'], stim_params['DFs']) elif 'DF' not in stim_params: t_type = 'DF' sim_queue = createSimQueue(stim_params['amps'], stim_params['durations'], [TITRATION_T_OFFSET] * len(stim_params['durations']), stim_params['PRFs'], [None]) nqueue = sim_queue.shape[0] t_var = ESTIM_params[t_type] # Create SolverElec instance solver = SolverElec() # Run titrations nsims = len(neurons) * nqueue simcount = 0 filepaths = [] for nname in neurons: neuron = getNeuronsDict()[nname]() try: for i in range(nqueue): simcount += 1 # Extract parameters Astim, tstim, toffset, PRF, DF = sim_queue[i, :] if Astim is None: Astim = (0., 2 * TITRATION_ESTIM_A_MAX) elif tstim is None: tstim = (0., 2 * TITRATION_T_MAX) elif DF is None: DF = (0., 2 * TITRATION_DF_MAX) curr_params = [Astim, tstim, PRF, DF] # Generate log str log_str = '' pnames = list(ESTIM_params.keys()) j = 0 for cp in curr_params: pn = pnames[j] pi = ESTIM_params[pn] if not isinstance(cp, tuple): if log_str: log_str += ', ' log_str += '{} = {:.2f} {}'.format(pn, pi['factor'] * cp, pi['unit']) j += 1 # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") # Log logger.info(ESTIM_titration_log, neuron.name, simcount, nsims, log_str) # Run titration tstart = time.time() (output_thr, t, y, states, lat) = titrateEStim(solver, neuron, Astim, tstim, toffset, PRF, DF) Vm, *channels = y tcomp = time.time() - tstart logger.info('completed in %.2f s, threshold = %.2f %s', tcomp, output_thr * t_var['factor'], t_var['unit']) # Determine output variable if t_type == 'A': Astim = output_thr elif t_type == 't': tstim = output_thr elif t_type == 'DF': DF = output_thr # Define output naming if DF == 1.0: simcode = ESTIM_CW_code.format(neuron.name, Astim, tstim * 1e3) else: simcode = ESTIM_PW_code.format(neuron.name, Astim, tstim * 1e3, PRF * 1e-3, DF) # Store data in dictionary data = { 'Astim': Astim, 'tstim': tstim, 'toffset': toffset, 'PRF': PRF, 'DF': DF, 't': t, 'states': states, 'Vm': Vm } for j in range(len(neuron.states_names)): data[neuron.states_names[j]] = channels[j] # Export data to PKL file output_filepath = batch_dir + '/' + simcode + ".pkl" with open(output_filepath, 'wb') as fh: pickle.dump(data, fh) logger.info('simulation data exported to "%s"', output_filepath) filepaths.append(output_filepath) # Detect spikes on Qm signal n_spikes, lat, sr = detectSpikes(t, Vm, SPIKE_MIN_VAMP, SPIKE_MIN_DT) logger.info('%u spike%s detected', n_spikes, "s" if n_spikes > 1 else "") # Export key metrics to log file log = { 'A': date_str, 'B': daytime_str, 'C': neuron.name, 'D': Astim, 'E': tstim * 1e3, 'F': PRF * 1e-3 if DF < 1 else 'N/A', 'G': DF, 'H': t.size, 'I': round(tcomp, 2), 'J': n_spikes, 'K': lat * 1e3 if isinstance(lat, float) else 'N/A', 'L': sr * 1e-3 if isinstance(sr, float) else 'N/A' } if xlslog(log_filepath, 'Data', log) == 1: logger.info('log exported to "%s"', log_filepath) else: logger.error('log export to "%s" aborted', log_filepath) except AssertionError as err: logger.error(err) return filepaths def runMechBatch(batch_dir, log_filepath, Cm0, Qm0, stim_params, a=default_diam, d=default_embedding): ''' Run batch simulations of the mechanical system with imposed values of charge density, for various sonophore spans and stimulation parameters. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param Cm0: membrane resting capacitance (F/m2) :param Qm0: membrane resting charge density (C/m2) :param stim_params: dictionary containing sweeps for all stimulation parameters :param a: BLS in-plane diameter (m) :param d: depth of embedding tissue around plasma membrane (m) ''' # Define logging format MECH_log = ('Mechanical simulation %u/%u (a = %.1f nm, d = %.1f um, f = %.2f kHz, ' 'A = %.2f kPa, Q = %.1f nC/cm2)') logger.info("Starting mechanical simulation batch") # Unpack stimulation parameters amps = np.array(stim_params['amps']) charges = np.array(stim_params['charges']) # Generate simulations queue nA = len(amps) nQ = len(charges) sim_queue = np.array(np.meshgrid(amps, charges)).T.reshape(nA * nQ, 2) nqueue = sim_queue.shape[0] # Run simulations nsims = len(stim_params['freqs']) * nqueue simcount = 0 filepaths = [] for Fdrive in stim_params['freqs']: try: # Create BilayerSonophore instance (modulus of embedding tissue depends on frequency) bls = BilayerSonophore(a, Fdrive, Cm0, Qm0, d) for i in range(nqueue): simcount += 1 Adrive, Qm = sim_queue[i, :] # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") # Log and define naming logger.info(MECH_log, simcount, nsims, a * 1e9, d * 1e6, Fdrive * 1e-3, Adrive * 1e-3, Qm * 1e5) simcode = MECH_code.format(a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, Qm * 1e5) # Run simulation tstart = time.time() (t, y, states) = bls.runMech(Fdrive, Adrive, Qm) (Z, ng) = y U = np.insert(np.diff(Z) / np.diff(t), 0, 0.0) tcomp = time.time() - tstart logger.info('completed in %.2f seconds', tcomp) # Store data in dictionary data = { 'a': a, 'd': d, 'Cm0': Cm0, 'Qm0': Qm0, 'Fdrive': Fdrive, 'Adrive': Adrive, 'phi': np.pi, 'Qm': Qm, 't': t, 'states': states, 'U': U, 'Z': Z, 'ng': ng } # Export data to PKL file output_filepath = batch_dir + '/' + simcode + ".pkl" with open(output_filepath, 'wb') as fh: pickle.dump(data, fh) logger.info('simulation data exported to "%s"', output_filepath) filepaths.append(output_filepath) # Compute key output metrics Zmax = np.amax(Z) Zmin = np.amin(Z) Zabs_max = np.amax(np.abs([Zmin, Zmax])) eAmax = bls.arealstrain(Zabs_max) Tmax = bls.TEtot(Zabs_max) Pmmax = bls.PMavgpred(Zmin) ngmax = np.amax(ng) dUdtmax = np.amax(np.abs(np.diff(U) / np.diff(t)**2)) # Export key metrics to log file log = { 'A': date_str, 'B': daytime_str, 'C': a * 1e9, 'D': d * 1e6, 'E': Fdrive * 1e-3, 'F': Adrive * 1e-3, 'G': Qm * 1e5, 'H': t.size, 'I': tcomp, 'J': bls.kA + bls.kA_tissue, 'K': Zmax * 1e9, 'L': eAmax, 'M': Tmax * 1e3, 'N': (ngmax - bls.ng0) / bls.ng0, 'O': Pmmax * 1e-3, 'P': dUdtmax } if xlslog(log_filepath, 'Data', log) == 1: logger.info('log exported to "%s"', log_filepath) else: logger.error('log export to "%s" aborted', log_filepath) except AssertionError as err: logger.error(err) return filepaths diff --git a/README.md b/README.md index 32fd67f..d192685 100644 --- a/README.md +++ b/README.md @@ -1,65 +1,68 @@ Description ============ PointNICE is a Python implementation of the **Neuronal Intramembrane Cavitation Excitation** (NICE) model introduced by Plaksin et. al in 2014 and initially developed in MATLAB by its authors. It contains optimized methods to predict the electrical response of point-neuron models to both acoustic and electrical stimuli. This package contains several core modules: - **bls** defines the underlying biomechanical model of intramembrane cavitation (**BilayerSonophore** class), and provides an integration method to predict compute the mechanical oscillations of the plasma membrane subject to a periodic acoustic perturbation. - **solvers** contains a simple solver for electrical stimuli (**SolverElec** class) as well as a tailored solver for acoustic stimuli (**SolverUS** class). The latter directly inherits from the BilayerSonophore class upon instantiation, and is hooked to a specific "channel mechanism" in order to link the mechanical model to an electrical model of membrane dynamics. It also provides several integration methods (detailed below) to compute the behaviour of the full electro-mechanical model subject to a continuous or pulsed ultrasonic stimulus. - **channels** contains the definitions of the different channels mechanisms inherent to specific neurons, including several types of **cortical** and **thalamic** neurons. - **plt** defines plotting utilities to load results of several simulations and display/compare temporal profiles of multiple variables of interest across simulations. - **utils** defines generic utilities used across the different modules The **SolverUS** class incorporates optimized numerical integration methods to perform dynamic simulations of the model subject to acoustic perturbation, and compute the evolution of its mechanical and electrical variables: - a **classic** method that solves all variables for the entire duration of the simulation. This method uses very small time steps and is computationally expensive (simulation time: several hours) - a **hybrid** method (initially developed by Plaskin et al.) in which integration is performed in consecutive “slices” of time, during which the full system is solved until mechanical stabilization, at which point the electrical system is solely solved with predicted mechanical variables until the end of the slice. This method is more efficient (simulation time: several minutes) and provides accurate results. - a newly developed **effective** method that neglects the high amplitude oscillations of mechanical and electrical variables during each acoustic cycle, to instead grasp the net effect of the acoustic stimulus on the electrical system. To do so, the sole electrical system is solved using pre-computed coefficients that depend on membrane charge and acoustic amplitude. This method allows to run simulations of the electrical system in only a few seconds, with very accurate results of the net membrane charge density evolution. This package is meant to be easy to use as a predictive and comparative tool for researchers investigating ultrasonic and/or electrical neuro-stimulation experimentally. -Installation & +Installation ================== Install Python 3 if not already done. Open a terminal. -Activate a Python3 environment if needed, e.g. on the tnesrv5 machine: -`>> source /opt/apps/anaconda3/bin activate` +Activate a Python3 environment if needed, e.g. on the tnesrv5 machine: + + source /opt/apps/anaconda3/bin activate Check that the appropriate version of pip is activated: -`>>> pip --version -pip x.x.x from /site-packages (python 3.x)` + + pip --version Go to the PointNICE directory (where the setup.py file is located) and install it as a package: -`>> cd ->> pip install -e .` + + cd + pip install -e . + PointNICE and all its dependencies will be installed. Usage ======= Command line scripts --------------------- To run single simulations of a given point-neuron model under specific stimulation parameters, you can use the `ASTIM_run.py` and `ESTIM_run.py` command-line scripts provided by the package. For instance, to simulate a regular-spiking neuron under continuous wave ultrasonic stimulation at 500kHz and 100kPa, for 150 ms: - `python ASTIM_run.py -n=RS -f=500 -A=100 -t=150` + python ASTIM_run.py -n=RS -f=500 -A=100 -t=150 Similarly, to simulate the electrical stimulation of a thalamo-cortical neuron at 10 mA/m2 for 150 ms: - `python ESTIM_run.py -n=TC -A=10 -t=150` + python ESTIM_run.py -n=TC -A=10 -t=150 The simulation results will be save in an output PKL file in the current working directory. To view these results, you can use the dedicated Batch scripts --------------- To run a batch of simulations on different neuron types and spanning ranges of several stimulation parameters, you can run the `ASTIM_batch.py` and `ESTIM_batch.py` scripts. To do so, simply modify the **stim_params** and **neurons** variables with your own neuron types and parameter sweeps, and then run the scripts (without command-line arguments). diff --git a/plot/plot_batch.py b/plot/plot_batch.py index f502c71..b50e741 100644 --- a/plot/plot_batch.py +++ b/plot/plot_batch.py @@ -1,33 +1,33 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-03-20 12:19:55 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-01-26 09:17:25 +# @Last Modified time: 2018-03-05 19:06:17 """ Batch plot profiles of several specific output variables of NICE simulations. """ import sys import logging from PointNICE.utils import logger, OpenFilesDialog from PointNICE.plt import plotBatch # Set logging level logger.setLevel(logging.INFO) # Select data files pkl_filepaths, pkl_dir = OpenFilesDialog('pkl') if not pkl_filepaths: logger.error('No input file') sys.exit(1) # Plot profiles try: # yvars = {'Q_m': ['Qm'], 'i_{Ca}\ kin.': ['s', 'u', 's2u'], 'I': ['iNa', 'iK', 'iT', 'iL']} yvars = {'Q_m': ['Qm']} - plotBatch(pkl_dir, pkl_filepaths, title=False, vars_dict=yvars) + plotBatch(pkl_dir, pkl_filepaths, title=False) except AssertionError as err: logger.error(err) sys.exit(1) diff --git a/plot/plot_comp.py b/plot/plot_comp.py index 015b05f..fca8d90 100644 --- a/plot/plot_comp.py +++ b/plot/plot_comp.py @@ -1,35 +1,35 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-02-13 12:41:26 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-02-27 15:14:54 +# @Last Modified time: 2018-03-06 22:40:02 """ Compare profiles of several specific output variables of NICE simulations. """ import sys import logging from PointNICE.utils import logger, OpenFilesDialog from PointNICE.plt import plotComp # Set logging level logger.setLevel(logging.INFO) # Select data files pkl_filepaths, _ = OpenFilesDialog('pkl') if not pkl_filepaths: logger.error('No input file') sys.exit(1) # Comparative plot try: yvars = ['Qm'] - # labels = ['classic', 'effective'] + labels = ['classic', 'effective'] # labels = ['FS neuron', 'LTS neuron', 'RE neuron', 'RS neuron', 'TC neuron'] - plotComp(yvars, pkl_filepaths) + plotComp(yvars, pkl_filepaths, labels=labels) except AssertionError as err: logger.error(err) sys.exit(1) diff --git a/sim/ASTIM_batch.py b/sim/ASTIM_batch.py index 6ef9656..efedf86 100644 --- a/sim/ASTIM_batch.py +++ b/sim/ASTIM_batch.py @@ -1,52 +1,52 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-02-13 18:16:09 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-02-27 15:15:46 +# @Last Modified time: 2018-03-05 17:54:35 """ Run batch acoustic simulations of specific "point-neuron" models. """ import sys import os import logging import numpy as np from PointNICE.utils import logger from PointNICE.solvers import setBatchDir, checkBatchLog, runAStimBatch from PointNICE.plt import plotBatch # Set logging level logger.setLevel(logging.DEBUG) # Neurons neurons = ['LeechP'] # Stimulation parameters stim_params = { 'freqs': [350e3], # Hz 'amps': [100e3], # Pa 'durations': [150e-3], # s 'PRFs': [100.0], # Hz 'DFs': [1] } stim_params['offsets'] = [100e-3] * len(stim_params['durations']) # s try: # Select output directory batch_dir = setBatchDir() log_filepath, _ = checkBatchLog(batch_dir, 'A-STIM') # Run A-STIM batch pkl_filepaths = runAStimBatch(batch_dir, log_filepath, neurons, stim_params, - int_method='effective') + int_method='classic') pkl_dir, _ = os.path.split(pkl_filepaths[0]) # Plot resulting profiles yvars = {'Q_m': ['Qm']} plotBatch(pkl_dir, pkl_filepaths, yvars) except AssertionError as err: logger.error(err) sys.exit(1) diff --git a/sim/MECH_batch.py b/sim/MECH_batch.py index 76d1db7..5be1001 100644 --- a/sim/MECH_batch.py +++ b/sim/MECH_batch.py @@ -1,48 +1,49 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-11-21 10:46:56 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-09-14 16:06:20 +# @Last Modified time: 2018-03-07 11:03:12 """ Run batch simulations of the NICE mechanical model with imposed charge densities """ import sys import os import logging import numpy as np from PointNICE.utils import logger from PointNICE.solvers import setBatchDir, checkBatchLog, runMechBatch from PointNICE.plt import plotBatch # Set logging level logger.setLevel(logging.DEBUG) # Electrical properties of the membrane Cm0 = 1e-2 # membrane resting capacitance (F/m2) Qm0 = -80e-5 # membrane resting charge density (C/m2) +a = 500e-9 # in-plane diameter (m) # Stimulation parameters stim_params = { 'freqs': [3.5e5], # Hz - 'amps': [100e3, 200e3, 500e3], # Pa - 'charges': [50e-5] # C/m2 + 'amps': [500e3], # Pa + 'charges': [-80e-5] # C/m2 } # Select output directory try: batch_dir = setBatchDir() log_filepath, _ = checkBatchLog(batch_dir, 'MECH') # Run MECH batch - pkl_filepaths = runMechBatch(batch_dir, log_filepath, Cm0, Qm0, stim_params) + pkl_filepaths = runMechBatch(batch_dir, log_filepath, Cm0, Qm0, stim_params, a) pkl_dir, _ = os.path.split(pkl_filepaths[0]) # Plot resulting profiles plotBatch(pkl_dir, pkl_filepaths) except AssertionError as err: logger.error(err) sys.exit(1) diff --git a/tests/test_values.py b/tests/test_values.py index d24961f..bc3ca47 100644 --- a/tests/test_values.py +++ b/tests/test_values.py @@ -1,244 +1,248 @@ #!/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-11-27 13:11:08 +# @Last Modified time: 2018-03-07 15:21:22 ''' Run functionalities of the package and test validity of outputs. ''' import sys import logging from argparse import ArgumentParser import numpy as np from PointNICE.utils import logger, getNeuronsDict from PointNICE import BilayerSonophore, SolverElec, SolverUS from PointNICE.solvers import detectSpikes, titrateEStim, titrateAStim from PointNICE.constants import * # Set logging level logger.setLevel(logging.INFO) # List of implemented neurons neurons = getNeuronsDict() neurons = list(neurons.values()) 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') # Create BLS instance a = 32e-9 # m Fdrive = 350e3 # Hz Cm0 = 1e-2 # membrane resting capacitance (F/m2) Qm0 = -80e-5 # membrane resting charge density (C/m2) bls = BilayerSonophore(a, Fdrive, Cm0, Qm0) # Run mechanical simulation Adrive = 100e3 # Pa Qm = 50e-5 # C/m2 _, 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') + logger.info('Starting test: neurons resting potential') # Initialize solver solver = SolverElec() for neuron_class in neurons: # Simulate each neuron in free conditions neuron = neuron_class() + + logger.info('%s neuron simulation in free conditions', neuron.name) + _, 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') + logger.info('Passed test: neurons resting potential') 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': 4.66, 'RE': 3.61, 'TC': 4.05} - latency_refs = {'FS': 101.00e-3, 'LTS': 128.56e-3, 'RS': 103.81e-3, 'LeechT': 21.32e-3, - 'RE': 148.50e-3, 'TC': 63.46e-3} + Athr_refs = {'FS': 6.91, 'LTS': 1.54, 'RS': 5.03, 'RE': 3.61, 'TC': 4.05, + 'LeechT': 4.66, 'LeechP': 13.72} + latency_refs = {'FS': 101.00e-3, 'LTS': 128.56e-3, 'RS': 103.81e-3, 'RE': 148.50e-3, + 'TC': 63.46e-3, 'LeechT': 21.32e-3, 'LeechP': 36.84e-3} for neuron_class in neurons: # Perform titration for each neuron neuron = neuron_class() logger.info('%s neuron titration', neuron.name) (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 diameter a = 32e-9 # m # 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': 21.48e3} - latency_refs = {'FS': 63.72e-3, 'LTS': 61.92e-3, 'RS': 62.52e-3, 'RE': 79.20e-3, - 'TC': 67.38e-3, 'LeechT': 23.9e-3} + Athr_refs = {'FS': 38.66e3, 'LTS': 24.90e3, 'RS': 50.90e3, 'RE': 46.36e3, 'TC': 23.29e3, + 'LeechT': 21.39e3, 'LeechP': 22.56e3} + latency_refs = {'FS': 69.58e-3, 'LTS': 57.56e-3, 'RS': 71.59e-3, 'RE': 79.20e-3, + 'TC': 63.67e-3, 'LeechT': 25.45e-3, 'LeechP': 54.76e-3} # Titration for each neuron for neuron_class in neurons: # Initialize neuron neuron = neuron_class() logger.info('%s neuron titration', neuron.name) # Initialize solver solver = SolverUS(a, 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_all(): logger.info('Starting tests') test_MECH() test_resting_potential() test_ESTIM() test_ASTIM() logger.info('All tests successfully passed') def main(): # Define valid test sets valid_testsets = [ 'MECH', 'resting_potential', 'ESTIM', 'ASTIM', 'all' ] # Define argument parser ap = ArgumentParser() ap.add_argument('-t', '--testset', type=str, default='all', choices=valid_testsets, help='Specific test set') ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity') # Parse arguments args = ap.parse_args() if args.verbose: logger.setLevel(logging.DEBUG) else: logger.setLevel(logging.INFO) # Run test if args.testset == 'all': test_all() else: possibles = globals().copy() possibles.update(locals()) method = possibles.get('test_{}'.format(args.testset)) method() sys.exit(0) if __name__ == '__main__': main()