diff --git a/PointNICE/bls.py b/PointNICE/bls.py index c26a56a..e35a346 100644 --- a/PointNICE/bls.py +++ b/PointNICE/bls.py @@ -1,685 +1,685 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-09-29 16:16:19 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-07-23 14:03:27 +# @Last Modified time: 2018-08-21 14:27:37 import inspect import logging import warnings import numpy as np import scipy.integrate as integrate from scipy.optimize import brentq, curve_fit from .utils import * from .constants import * # Get package logger logger = logging.getLogger('PointNICE') class BilayerSonophore: """ This class contains the geometric and mechanical parameters of the Bilayer Sonophore Model, as well as all the core functions needed to compute the dynamics (kinetics and kinematics) of the bilayer membrane cavitation, and run dynamic BLS simulations. """ # BIOMECHANICAL PARAMETERS T = 309.15 # Temperature (K) Rg = 8.314 # Universal gas constant (Pa.m^3.mol^-1.K^-1) delta0 = 2.0e-9 # Thickness of the leaflet (m) Delta_ = 1.4e-9 # Initial gap between the two leaflets on a non-charged membrane at equil. (m) pDelta = 1.0e5 # Attraction/repulsion pressure coefficient (Pa) m = 5.0 # Exponent in the repulsion term (dimensionless) n = 3.3 # Exponent in the attraction term (dimensionless) rhoL = 1075.0 # Density of the surrounding fluid (kg/m^3) muL = 7.0e-4 # Dynamic viscosity of the surrounding fluid (Pa.s) muS = 0.035 # Dynamic viscosity of the leaflet (Pa.s) kA = 0.24 # Area compression modulus of the leaflet (N/m) alpha = 7.56 # Tissue shear loss modulus frequency coefficient (Pa.s) C0 = 0.62 # Initial gas molar concentration in the surrounding fluid (mol/m^3) kH = 1.613e5 # Henry's constant (Pa.m^3/mol) P0 = 1.0e5 # Static pressure in the surrounding fluid (Pa) Dgl = 3.68e-9 # Diffusion coefficient of gas in the fluid (m^2/s) xi = 0.5e-9 # Boundary layer thickness for gas transport across leaflet (m) c = 1515.0 # Speed of sound in medium (m/s) # BIOPHYSICAL PARAMETERS epsilon0 = 8.854e-12 # Vacuum permittivity (F/m) epsilonR = 1.0 # Relative permittivity of intramembrane cavity (dimensionless) def __init__(self, diameter, Fdrive, Cm0, Qm0, embedding_depth=0.0): """ Constructor of the class. :param diameter: in-plane diameter of the sonophore structure within the membrane (m) :param Fdrive: frequency of acoustic perturbation (Hz) :param Cm0: membrane resting capacitance (F/m2) :param Qm0: membrane resting charge density (C/m2) :param embedding_depth: depth of the embedding tissue around the membrane (m) """ logger.debug('%.1f nm BLS initialization at %.2f kHz, %.2f nC/cm2', diameter * 1e9, Fdrive * 1e-3, Qm0 * 1e5) # Extract resting constants and geometry self.Cm0 = Cm0 self.Qm0 = Qm0 self.a = diameter self.d = embedding_depth self.S0 = np.pi * self.a**2 # Derive frequency-dependent tissue elastic modulus G_tissue = self.alpha * Fdrive # G'' (Pa) self.kA_tissue = 2 * G_tissue * self.d # kA of the tissue layer (N/m) # Check existence of lookups for derived parameters lookups = get_BLS_lookups(self.a) Qkey = '{:.2f}'.format(Qm0 * 1e5) # If no lookup, compute parameters and store them in lookup if not lookups or Qkey not in lookups: # Find Delta that cancels out Pm + Pec at Z = 0 (m) if self.Qm0 == 0.0: D_eq = self.Delta_ else: (D_eq, Pnet_eq) = self.findDeltaEq(self.Qm0) assert Pnet_eq < PNET_EQ_MAX, 'High Pnet at Z = 0 with ∆ = %.2f nm' % (D_eq * 1e9) self.Delta = D_eq # Find optimal Lennard-Jones parameters to approximate PMavg (LJ_approx, std_err, _) = self.LJfitPMavg() assert std_err < PMAVG_STD_ERR_MAX, 'High error in PmAvg nonlinear fit:'\ ' std_err = %.2f Pa' % std_err self.LJ_approx = LJ_approx lookups[Qkey] = {'LJ_approx': LJ_approx, 'Delta_eq': D_eq} logger.debug('Saving BLS derived parameters to lookup file') save_BLS_lookups(self.a, lookups) # If lookup exists, load parameters from it else: logger.debug('Loading BLS derived parameters from lookup file') self.LJ_approx = lookups[Qkey]['LJ_approx'] self.Delta = lookups[Qkey]['Delta_eq'] # Compute initial volume and gas content self.V0 = np.pi * self.Delta * self.a**2 self.ng0 = self.gasPa2mol(self.P0, self.V0) def __str__(self): s = '-------- Bilayer Sonophore --------\n' s += 'class attributes:\n' class_attrs = inspect.getmembers(self.__class__, lambda a: not(inspect.isroutine(a))) class_attrs = [a for a in class_attrs if not(a[0].startswith('__') and a[0].endswith('__'))] for ca in class_attrs: s += '{} = {}\n'.format(ca[0], ca[1]) s += 'instance attributes:\n' inst_attrs = inspect.getmembers(self, lambda a: not(inspect.isroutine(a))) inst_attrs = [a for a in inst_attrs if not(a[0].startswith('__') and a[0].endswith('__')) and a not in class_attrs] for ia in inst_attrs: s += '{} = {}\n'.format(ia[0], ia[1]) return s def reinit(self): logger.debug('Re-initializing BLS object') # Find Delta that cancels out Pm + Pec at Z = 0 (m) if self.Qm0 == 0.0: D_eq = self.Delta_ else: (D_eq, Pnet_eq) = self.findDeltaEq(self.Qm0) assert Pnet_eq < PNET_EQ_MAX, 'High Pnet at Z = 0 with ∆ = %.2f nm' % (D_eq * 1e9) self.Delta = D_eq # Compute initial volume and gas content self.V0 = np.pi * self.Delta * self.a**2 self.ng0 = self.gasPa2mol(self.P0, self.V0) def curvrad(self, Z): """ Return the (signed) instantaneous curvature radius of the leaflet. :param Z: leaflet apex outward deflection value (m) :return: leaflet curvature radius (m) """ if Z == 0.0: return np.inf else: return (self.a**2 + Z**2) / (2 * Z) def v_curvrad(self, Z): ''' Vectorized curvrad function ''' - return np.array([self.curvrad(z) for z in Z]) + return np.array(list(map(self.curvrad, Z))) def surface(self, Z): """ Return the surface area of the stretched leaflet (spherical cap). :param Z: leaflet apex outward deflection value (m) :return: surface of the stretched leaflet (m^2) """ return np.pi * (self.a**2 + Z**2) def volume(self, Z): """ Return the total volume of the inter-leaflet space (cylinder +/- spherical cap). :param Z: leaflet apex outward deflection value (m) :return: inner volume of the bilayer sonophore structure (m^3) """ return np.pi * self.a**2 * self.Delta\ * (1 + (Z / (3 * self.Delta) * (3 + Z**2 / self.a**2))) def arealstrain(self, Z): """ Compute the areal strain of the stretched leaflet. epsilon = (S - S0)/S0 = (Z/a)^2 :param Z: leaflet apex outward deflection value (m) :return: areal strain (dimensionless) """ return (Z / self.a)**2 def Capct(self, Z): """ Compute the membrane capacitance per unit area, under the assumption of parallel-plate capacitor with average inter-layer distance. :param Z: leaflet apex outward deflection value (m) :return: capacitance per unit area (F/m2) """ if Z == 0.0: return self.Cm0 else: return ((self.Cm0 * self.Delta / self.a**2) * (Z + (self.a**2 - Z**2 - Z * self.Delta) / (2 * Z) * np.log((2 * Z + self.Delta) / self.Delta))) def v_Capct(self, Z): ''' Vectorized Capct function ''' - return np.array([self.Capct(z) for z in Z]) + return np.array(list(map(self.Capct, Z))) def derCapct(self, Z, U): """ Compute the derivative of the membrane capacitance per unit area with respect to time, under the assumption of parallel-plate capacitor. :param Z: leaflet apex outward deflection value (m) :param U: leaflet apex outward deflection velocity (m/s) :return: derivative of capacitance per unit area (F/m2.s) """ dCmdZ = ((self.Cm0 * self.Delta / self.a**2) * ((Z**2 + self.a**2) / (Z * (2 * Z + self.Delta)) - ((Z**2 + self.a**2) * np.log((2 * Z + self.Delta) / self.Delta)) / (2 * Z**2))) return dCmdZ * U def localdef(self, r, Z, R): """ Compute the (signed) local transverse leaflet deviation at a distance r from the center of the dome. :param r: in-plane distance from center of the sonophore (m) :param Z: leaflet apex outward deflection value (m) :param R: leaflet curvature radius (m) :return: local transverse leaflet deviation (m) """ if np.abs(Z) == 0.0: return 0.0 else: return np.sign(Z) * (np.sqrt(R**2 - r**2) - np.abs(R) + np.abs(Z)) def Pacoustic(self, t, Adrive, Fdrive, phi=np.pi): """ Compute the acoustic pressure at a specific time, given the amplitude, frequency and phase of the acoustic stimulus. :param t: time of interest :param Adrive: acoustic drive amplitude (Pa) :param Fdrive: acoustic drive frequency (Hz) :param phi: acoustic drive phase (rad) """ return Adrive * np.sin(2 * np.pi * Fdrive * t - phi) def PMlocal(self, r, Z, R): """ Compute the local intermolecular pressure. :param r: in-plane distance from center of the sonophore (m) :param Z: leaflet apex outward deflection value (m) :param R: leaflet curvature radius (m) :return: local intermolecular pressure (Pa) """ z = self.localdef(r, Z, R) relgap = (2 * z + self.Delta) / self.Delta_ return self.pDelta * ((1 / relgap)**self.m - (1 / relgap)**self.n) def PMavg(self, Z, R, S): """ Compute the average intermolecular pressure felt across the leaflet by quadratic integration. :param Z: leaflet apex outward deflection value (m) :param R: leaflet curvature radius (m) :param S: surface of the stretched leaflet (m^2) :return: averaged intermolecular resultant pressure across the leaflet (Pa) .. warning:: quadratic integration is computationally expensive. """ # Intermolecular force over an infinitely thin ring of radius r fMring = lambda r, Z, R: 2 * np.pi * r * self.PMlocal(r, Z, R) # Integrate from 0 to a fTotal, _ = integrate.quad(fMring, 0, self.a, args=(Z, R)) return fTotal / S def v_PMavg(self, Z, R, S): ''' Vectorized PMavg function ''' - return np.array([self.PMavg(Z[i], R[i], S[i]) for i in range(Z.size)]) + return np.array(list(map(self.PMavg, Z, R, S))) def LJfitPMavg(self): """ Determine optimal parameters of a Lennard-Jones expression approximating the average intermolecular pressure. These parameters are obtained by a nonlinear fit of the Lennard-Jones function for a range of deflection values between predetermined Zmin and Zmax. :return: 3-tuple with optimized LJ parameters for PmAvg prediction (Map) and the standard and max errors of the prediction in the fitting range (in Pascals) """ # Determine lower bound of deflection range: when Pm = Pmmax PMmax = LJFIT_PM_MAX # Pa Zminlb = -0.49 * self.Delta Zminub = 0.0 f = lambda Z, Pmmax: self.PMavg(Z, self.curvrad(Z), self.surface(Z)) - PMmax Zmin = brentq(f, Zminlb, Zminub, args=(PMmax), xtol=1e-16) # Create vectors for geometric variables Zmax = 2 * self.a Z = np.arange(Zmin, Zmax, 1e-11) Pmavg = self.v_PMavg(Z, self.v_curvrad(Z), self.surface(Z)) # Compute optimal nonlinear fit of custom LJ function with initial guess x0_guess = self.delta0 C_guess = 0.1 * self.pDelta nrep_guess = self.m nattr_guess = self.n pguess = (x0_guess, C_guess, nrep_guess, nattr_guess) popt, _ = curve_fit(lambda x, x0, C, nrep, nattr: LennardJones(x, self.Delta, x0, C, nrep, nattr), Z, Pmavg, p0=pguess, maxfev=10000) (x0_opt, C_opt, nrep_opt, nattr_opt) = popt Pmavg_fit = LennardJones(Z, self.Delta, x0_opt, C_opt, nrep_opt, nattr_opt) # Compute prediction error residuals = Pmavg - Pmavg_fit ss_res = np.sum(residuals**2) N = residuals.size std_err = np.sqrt(ss_res / N) max_err = max(np.abs(residuals)) logger.debug('LJ approx: x0 = %.2f nm, C = %.2f kPa, m = %.2f, n = %.2f', x0_opt * 1e9, C_opt * 1e-3, nrep_opt, nattr_opt) LJ_approx = {"x0": x0_opt, "C": C_opt, "nrep": nrep_opt, "nattr": nattr_opt} return (LJ_approx, std_err, max_err) def PMavgpred(self, Z): """ Return the predicted intermolecular pressure based on a specific Lennard-Jones function fitted on the deflection physiological range. :param Z: leaflet apex outward deflection value (m) :return: predicted average intermolecular pressure (Pa) """ return LennardJones(Z, self.Delta, self.LJ_approx['x0'], self.LJ_approx['C'], self.LJ_approx['nrep'], self.LJ_approx['nattr']) def Pelec(self, Z, Qm): """ Compute the electric equivalent pressure term. :param Z: leaflet apex outward deflection value (m) :param Qm: membrane charge density (C/m2) :return: electric equivalent pressure (Pa) """ relS = self.S0 / self.surface(Z) abs_perm = self.epsilon0 * self.epsilonR # F/m return -relS * Qm**2 / (2 * abs_perm) # Pa def findDeltaEq(self, Qm): """ Compute the Delta that cancels out the (Pm + Pec) equation at Z = 0 for a given membrane charge density, using the Brent method to refine the pressure root iteratively. :param Qm: membrane charge density (C/m2) :return: equilibrium value (m) and associated pressure (Pa) """ f = lambda Delta: (self.pDelta * ((self.Delta_ / Delta)**self.m - (self.Delta_ / Delta)**self.n) + self.Pelec(0.0, Qm)) Delta_lb = 0.1 * self.Delta_ Delta_ub = 2.0 * self.Delta_ Delta_eq = brentq(f, Delta_lb, Delta_ub, xtol=1e-16) logger.debug('∆eq = %.2f nm', Delta_eq * 1e9) return (Delta_eq, f(Delta_eq)) def gasflux(self, Z, P): """ Compute the gas molar flux through the BLS boundary layer for an unsteady system. :param Z: leaflet apex outward deflection value (m) :param P: internal gas pressure in the inter-leaflet space (Pa) :return: gas molar flux (mol/s) """ dC = self.C0 - P / self.kH return 2 * self.surface(Z) * self.Dgl * dC / self.xi def gasmol2Pa(self, ng, V): """ Compute the gas pressure in the inter-leaflet space for an unsteady system, from the value of gas molar content. :param ng: internal molar content (mol) :param V: inner volume of the bilayer sonophore structure (m^3) :return: internal gas pressure (Pa) """ return ng * self.Rg * self.T / V def gasPa2mol(self, P, V): """ Compute the gas molar content in the inter-leaflet space for an unsteady system, from the value of internal gas pressure. :param P: internal gas pressure in the inter-leaflet space (Pa) :param V: inner volume of the bilayer sonophore structure (m^3) :return: internal gas molar content (mol) """ return P * V / (self.Rg * self.T) def PtotQS(self, Z, ng, Qm, Pac, Pm_comp_method): """ Compute the balance pressure of the quasi-steady system, upon application of an external perturbation on a charged membrane: Ptot = Pm + Pg + Pec - P0 - Pac. :param Z: leaflet apex outward deflection value (m) :param ng: internal molar content (mol) :param Qm: membrane charge density (C/m2) :param Pac: external acoustic perturbation (Pa) :param Pm_comp_method: type of method used to compute average intermolecular pressure :return: total balance pressure (Pa) """ if Pm_comp_method is PmCompMethod.direct: Pm = self.PMavg(Z, self.curvrad(Z), self.surface(Z)) elif Pm_comp_method is PmCompMethod.predict: Pm = self.PMavgpred(Z) return Pm + self.gasmol2Pa(ng, self.volume(Z)) - self.P0 - Pac + self.Pelec(Z, Qm) def balancedefQS(self, ng, Qm, Pac=0.0, Pm_comp_method=PmCompMethod.predict): """ Compute the leaflet deflection upon application of an external perturbation to a quasi-steady system with a charged membrane. This function uses the Brent method (progressive approximation of function root) to solve the following transcendental equation for Z: Pm + Pg + Pec - P0 - Pac = 0. :param ng: internal molar content (mol) :param Qm: membrane charge density (C/m2) :param Pac: external acoustic perturbation (Pa) :param Pm_comp_method: type of method used to compute average intermolecular pressure :return: leaflet deflection (Z) canceling out the balance equation """ lb = -0.49 * self.Delta ub = self.a Plb = self.PtotQS(lb, ng, Qm, Pac, Pm_comp_method) Pub = self.PtotQS(ub, ng, Qm, Pac, Pm_comp_method) assert (Plb > 0 > Pub), '[%d, %d] is not a sign changing interval for PtotQS' % (lb, ub) return brentq(self.PtotQS, lb, ub, args=(ng, Qm, Pac, Pm_comp_method), xtol=1e-16) def TEleaflet(self, Z): """ Compute the circumferential elastic tension felt across the entire leaflet upon stretching. :param Z: leaflet apex outward deflection value (m) :return: circumferential elastic tension (N/m) """ return self.kA * self.arealstrain(Z) def TEtissue(self, Z): """ Compute the circumferential elastic tension felt across the embedding viscoelastic tissue layer upon stretching. :param Z: leaflet apex outward deflection value (m) :return: circumferential elastic tension (N/m) """ return self.kA_tissue * self.arealstrain(Z) def TEtot(self, Z): """ Compute the total circumferential elastic tension (leaflet and embedding tissue) felt upon stretching. :param Z: leaflet apex outward deflection value (m) :return: circumferential elastic tension (N/m) """ return self.TEleaflet(Z) + self.TEtissue(Z) def PEtot(self, Z, R): """ Compute the total elastic tension pressure (leaflet + embedding tissue) felt upon stretching. :param Z: leaflet apex outward deflection value (m) :param R: leaflet curvature radius (m) :return: elastic tension pressure (Pa) """ return - self.TEtot(Z) / R def PVleaflet(self, U, R): """ Compute the viscous stress felt across the entire leaflet upon stretching. :param U: leaflet apex outward deflection velocity (m/s) :param R: leaflet curvature radius (m) :return: leaflet viscous stress (Pa) """ return - 12 * U * self.delta0 * self.muS / R**2 def PVfluid(self, U, R): """ Compute the viscous stress felt across the entire fluid upon stretching. :param U: leaflet apex outward deflection velocity (m/s) :param R: leaflet curvature radius (m) :return: fluid viscous stress (Pa) """ return - 4 * U * self.muL / np.abs(R) def accP(self, Pres, R): """ Compute the pressure-driven acceleration of the leaflet in the unsteady system, upon application of an external perturbation. :param Pres: net resultant pressure (Pa) :param R: leaflet curvature radius (m) :return: pressure-driven acceleration (m/s^2) """ return Pres / (self.rhoL * np.abs(R)) def accNL(self, U, R): """ Compute the non-linear term of the leaflet acceleration in the unsteady system, upon application of an external perturbation. :param U: leaflet apex outward deflection velocity (m/s) :param R: leaflet curvature radius (m) :return: nonlinear acceleration (m/s^2) .. note:: A simplified version of nonlinear acceleration (neglecting dR/dH) is used here. """ # return - (3/2 - 2*R/H) * U**2 / R return -(3 * U**2) / (2 * R) def eqMech(self, y, t, Adrive, Fdrive, Qm, phi, Pm_comp_method=PmCompMethod.predict): """ Compute the derivatives of the 3-ODE mechanical system variables, with an imposed constant charge density. :param y: vector of HH system variables at time t :param t: specific instant in time (s) :param Adrive: acoustic drive amplitude (Pa) :param Fdrive: acoustic drive frequency (Hz) :param Qm: membrane charge density (F/m2) :param phi: acoustic drive phase (rad) :param Pm_comp_method: type of method used to compute average intermolecular pressure :return: vector of mechanical system derivatives at time t """ # Split input vector explicitly (U, Z, ng) = y # Correct deflection value is below critical compression if Z < -0.5 * self.Delta: logger.warning('Deflection out of range: Z = %.2f nm', Z * 1e9) Z = -0.49 * self.Delta # Compute curvature radius R = self.curvrad(Z) # Compute total pressure Pg = self.gasmol2Pa(ng, self.volume(Z)) if Pm_comp_method is PmCompMethod.direct: Pm = self.PMavg(Z, self.curvrad(Z), self.surface(Z)) elif Pm_comp_method is PmCompMethod.predict: Pm = self.PMavgpred(Z) Ptot = (Pm + Pg - self.P0 - self.Pacoustic(t, Adrive, Fdrive, phi) + self.PEtot(Z, R) + self.PVleaflet(U, R) + self.PVfluid(U, R) + self.Pelec(Z, Qm)) # Compute derivatives dUdt = self.accP(Ptot, R) + self.accNL(U, R) dZdt = U dngdt = self.gasflux(Z, Pg) # Return derivatives vector return [dUdt, dZdt, dngdt] def run(self, Fdrive, Adrive, Qm, phi=np.pi, Pm_comp_method=PmCompMethod.predict): """ Compute short solutions of the mechanical system for specific US stimulation parameters and with an imposed membrane charge density. :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param phi: acoustic drive phase (rad) :param Qm: imposed membrane charge density (C/m2) :param Pm_comp_method: type of method used to compute average intermolecular pressure :return: 3-tuple with the time profile, the solution matrix and a state vector """ # Check validity of stimulation parameters if not all(isinstance(param, float) for param in [Fdrive, Adrive, Qm, phi]): raise InputError('Invalid stimulation parameters (must be float typed)') if Fdrive <= 0: raise InputError('Invalid US driving frequency: {} kHz (must be strictly positive)' .format(Fdrive * 1e-3)) if Adrive < 0: raise InputError('Invalid US pressure amplitude: {} kPa (must be positive or null)' .format(Adrive * 1e-3)) if Qm < CHARGE_RANGE[0] or Qm > CHARGE_RANGE[1]: raise InputError('Invalid applied charge: {} nC/cm2 (must be within [{}, {}] interval' .format(Qm * 1e5, CHARGE_RANGE[0] * 1e5, CHARGE_RANGE[1] * 1e5)) if phi < 0 or phi >= 2 * np.pi: raise InputError('Invalid US pressure phase: {:.2f} rad (must be within [0, 2 PI[ rad' .format(phi)) # Raise warnings as error warnings.filterwarnings('error') # Determine mechanical system time step Tdrive = 1 / Fdrive dt_mech = Tdrive / NPC_FULL t_mech_cycle = np.linspace(0, Tdrive - dt_mech, NPC_FULL) # Initialize system variables t0 = 0.0 Z0 = 0.0 U0 = 0.0 ng0 = self.ng0 # Solve quasi-steady equation to compute first deflection value Pac1 = self.Pacoustic(t0 + dt_mech, Adrive, Fdrive, phi) Z1 = self.balancedefQS(ng0, Qm, Pac1, Pm_comp_method) U1 = (Z1 - Z0) / dt_mech # Construct arrays to hold system variables states = np.array([1, 1]) t = np.array([t0, t0 + dt_mech]) y = np.array([[U0, U1], [Z0, Z1], [ng0, ng0]]) # Integrate mechanical system for a few acoustic cycles until stabilization j = 0 ng_last = None Z_last = None periodic_conv = False while not periodic_conv and j < NCYCLES_MAX: t_mech = t_mech_cycle + t[-1] + dt_mech y_mech = integrate.odeint(self.eqMech, y[:, -1], t_mech, args=(Adrive, Fdrive, Qm, phi, Pm_comp_method)).T # Compare Z and ng signals over the last 2 acoustic periods if j > 0: Z_rmse = rmse(Z_last, y_mech[1, :]) ng_rmse = rmse(ng_last, y_mech[2, :]) logger.debug('step %u: Z_rmse = %.2e m, ng_rmse = %.2e mol', j, Z_rmse, ng_rmse) if Z_rmse < Z_ERR_MAX and ng_rmse < NG_ERR_MAX: periodic_conv = True # Update last vectors for next comparison Z_last = y_mech[1, :] ng_last = y_mech[2, :] # Concatenate time and solutions to global vectors states = np.concatenate([states, np.ones(NPC_FULL)], axis=0) t = np.concatenate([t, t_mech], axis=0) y = np.concatenate([y, y_mech], axis=1) # Increment loop index j += 1 if j == NCYCLES_MAX: logger.warning('No convergence: stopping after %u cycles', j) else: logger.debug('Periodic convergence after %u cycles', j) states[-1] = 0 # return output variables return (t, y[1:, :], states) diff --git a/PointNICE/lookups/BLS_lookups_a50.0nm.json b/PointNICE/lookups/BLS_lookups_a50.0nm.json index 7337d0d..38ee4b2 100644 --- a/PointNICE/lookups/BLS_lookups_a50.0nm.json +++ b/PointNICE/lookups/BLS_lookups_a50.0nm.json @@ -1 +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 +{"-80.00": {"LJ_approx": {"x0": 1.7845793287878114e-09, "C": 14600.797866604744, "nrep": 3.9112835250651723, "nattr": 0.9433189447961238}, "Delta_eq": 1.2344323203763867e-09}, "-71.90": {"LJ_approx": {"x0": 1.7114794411958874e-09, "C": 16732.841575829825, "nrep": 3.914827883392826, "nattr": 0.9781401389550711}, "Delta_eq": 1.2553492695740507e-09}} \ No newline at end of file diff --git a/PointNICE/neurons/nmodl/FS.mod b/PointNICE/neurons/nmodl/FS.mod deleted file mode 100644 index cbd9964..0000000 --- a/PointNICE/neurons/nmodl/FS.mod +++ /dev/null @@ -1,189 +0,0 @@ -TITLE FS neuron -: -: Equations governing the behavior of a fast spiking neuron -: during ultrasonic stimulation, according to the NICE model. -: -: Written by Theo Lemaire and Simon Narduzzi, EPFL, 2018 -: Contact: theo.lemaire@epfl.ch - - -INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)} - -UNITS { - (mA) = (milliamp) - (mV) = (millivolt) - (uF) = (microfarad) -} - -NEURON { - : Definition of NEURON mechanism - - : mechanism name and constituting currents - SUFFIX FS - NONSPECIFIC_CURRENT iNa - NONSPECIFIC_CURRENT iKd - NONSPECIFIC_CURRENT iM - NONSPECIFIC_CURRENT iLeak - - : simulation parameters (stimulus and dt) - RANGE duration, PRF, DC, stimon - - : physiological variables - RANGE Q, Vmeff, gnabar, gkdbar, gmbar, gleak, eleak, ek, ena -} - - -PARAMETER { - : Parameters set by the user upon initialization - - : simulation parameters (stimulus and dt) - duration (ms) - PRF (Hz) - DC - dt (ms) - - : membrane properties - cm = 1 (uF/cm2) - ena = 50 (mV) - ek = -90 (mV) - eleak = -70.4 (mV) - gnabar = 0.058 (mho/cm2) - gkdbar = 0.0039 (mho/cm2) - gmbar = 7.87e-5 (mho/cm2) - gleak = 3.8e-5 (mho/cm2) -} - -STATE { - : Differential variables other than v, i.e. the ion channels gating states - - m : iNa activation gate - h : iNa inactivation gate - n : iKd activation gate - p : iM activation gate -} - -ASSIGNED { - : Variables computed during the simulation and whose value can be retrieved - Q (nC/cm2) - Vmeff (mV) - v (mV) - iNa (mA/cm2) - iKd (mA/cm2) - iM (mA/cm2) - iLeak (mA/cm2) - alpha_h (/ms) - beta_h (/ms) - alpha_m (/ms) - beta_m (/ms) - alpha_n (/ms) - beta_n (/ms) - alpha_p (/ms) - beta_p (/ms) - stimon - tint (ms) -} - -INCLUDE "stimonoff.hoc" - -BREAKPOINT { - : Main computation block - - : compute Q - Q = v * cm - - : update stimulation boolean - stimonoff() - - : integrate states - SOLVE states METHOD cnexp - - : compute Vmeff - if(stimon) {Vmeff = veff_on(Q)} else {Vmeff = veff_off(Q)} - - : compute ionic currents - iNa = gnabar * m * m * m * h * (Vmeff - ena) - iKd = gkdbar * n * n * n * n * (Vmeff - ek) - iM = gmbar * p * (Vmeff - ek) - iLeak = gleak * (Vmeff - eleak) -} - -DERIVATIVE states { - : Compute states derivatives - - : compute ON or OFF rate constants accordingly to stimon value - if(stimon) {interpolate_on(Q)} else {interpolate_off(Q)} - - : compute gating states derivatives - m' = alpha_m * (1 - m) - beta_m * m - h' = alpha_h * (1 - h) - beta_h * h - n' = alpha_n * (1 - n) - beta_n * n - p' = alpha_p * (1 - p) - beta_p * p -} - -UNITSOFF : turning off units checking for functions not working with standard SI units - -INITIAL { - : Initialize variables and parameters - - : compute Q - Q = v * cm - - : set initial states values - m = alpham_off(Q) / (alpham_off(Q) + betam_off(Q)) - h = alphah_off(Q) / (alphah_off(Q) + betah_off(Q)) - n = alphan_off(Q) / (alphan_off(Q) + betan_off(Q)) - p = alphap_off(Q) / (alphap_off(Q) + betap_off(Q)) - - : initialize tint, set stimulation state and correct PRF - tint = 0 - stimon = 1 - PRF = PRF / 2 : for some reason PRF must be halved in order to ensure proper on/off switch -} - -: Define function tables for variables to be interpolated during ON and OFF periods -FUNCTION_TABLE veff_on(x) (mV) -FUNCTION_TABLE veff_off(x) (mV) -FUNCTION_TABLE alpham_on(x) (/ms) -FUNCTION_TABLE alpham_off(x) (/ms) -FUNCTION_TABLE betam_on(x) (/ms) -FUNCTION_TABLE betam_off(x) (/ms) -FUNCTION_TABLE alphah_on(x) (/ms) -FUNCTION_TABLE alphah_off(x) (/ms) -FUNCTION_TABLE betah_on(x) (/ms) -FUNCTION_TABLE betah_off(x) (/ms) -FUNCTION_TABLE alphan_on(x) (/ms) -FUNCTION_TABLE alphan_off(x) (/ms) -FUNCTION_TABLE betan_on(x) (/ms) -FUNCTION_TABLE betan_off(x) (/ms) -FUNCTION_TABLE alphap_on(x) (/ms) -FUNCTION_TABLE alphap_off(x) (/ms) -FUNCTION_TABLE betap_on(x) (/ms) -FUNCTION_TABLE betap_off(x) (/ms) - -PROCEDURE interpolate_on(Q){ - : Interpolate rate constants durong US-ON periods from appropriate tables - - alpha_m = alpham_on(Q) - beta_m = betam_on(Q) - alpha_h = alphah_on(Q) - beta_h = betah_on(Q) - alpha_n = alphan_on(Q) - beta_n = betan_on(Q) - alpha_p = alphap_on(Q) - beta_p = betap_on(Q) -} - -PROCEDURE interpolate_off(Q){ - : Interpolate rate constants durong US-OFF periods from appropriate tables - - alpha_m = alpham_off(Q) - beta_m = betam_off(Q) - alpha_h = alphah_off(Q) - beta_h = betah_off(Q) - alpha_n = alphan_off(Q) - beta_n = betan_off(Q) - alpha_p = alphap_off(Q) - beta_p = betap_off(Q) -} - -UNITSON : turn back units checking \ No newline at end of file diff --git a/PointNICE/neurons/nmodl/LTS.mod b/PointNICE/neurons/nmodl/LTS.mod deleted file mode 100644 index 6c4cc9c..0000000 --- a/PointNICE/neurons/nmodl/LTS.mod +++ /dev/null @@ -1,220 +0,0 @@ -TITLE LTS neuron -: -: Equations governing the behavior of a low-threshold spiking neuron -: during ultrasonic stimulation, according to the NICE model. -: -: Written by Theo Lemaire and Simon Narduzzi, EPFL, 2018 -: Contact: theo.lemaire@epfl.ch - - -INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)} - -UNITS { - (mA) = (milliamp) - (mV) = (millivolt) - (uF) = (microfarad) -} - -NEURON { - : Definition of NEURON mechanism - - : mechanism name and constituting currents - SUFFIX LTS - NONSPECIFIC_CURRENT iNa - NONSPECIFIC_CURRENT iKd - NONSPECIFIC_CURRENT iM - NONSPECIFIC_CURRENT iCa - NONSPECIFIC_CURRENT iLeak - - : simulation parameters (stimulus and dt) - RANGE duration, PRF, DC, stimon - - : physiological variables - RANGE Q, Vmeff, gnabar, gkdbar, gmbar, gcabar, gleak, eleak, ek, ena, eca -} - - -PARAMETER { - : Parameters set by the user upon initialization - - : simulation parameters (stimulus and dt) - duration (ms) - PRF (Hz) - DC - dt (ms) - - : membrane properties - cm = 1 (uF/cm2) - ena = 50 (mV) - eca = 120 (mV) - ek = -90 (mV) - eleak = -50.0 (mV) - gnabar = 0.050 (mho/cm2) - gkdbar = 0.004 (mho/cm2) - gmbar = 2.8e-5 (mho/cm2) - gcabar = 0.0004 (mho/cm2) - gleak = 1.9e-5 (mho/cm2) -} - -STATE { - : Differential variables other than v, i.e. the ion channels gating states - - m : iNa activation gate - h : iNa inactivation gate - n : iKd activation gate - p : iM activation gate - s : iCa activation gate - u : iCa inactivation gate -} - -ASSIGNED { - : Variables computed during the simulation and whose value can be retrieved - Q (nC/cm2) - Vmeff (mV) - v (mV) - iNa (mA/cm2) - iKd (mA/cm2) - iM (mA/cm2) - iCa (mA/cm2) - iLeak (mA/cm2) - alpha_h (/ms) - beta_h (/ms) - alpha_m (/ms) - beta_m (/ms) - alpha_n (/ms) - beta_n (/ms) - alpha_p (/ms) - beta_p (/ms) - alpha_s (/ms) - beta_s (/ms) - alpha_u (/ms) - beta_u (/ms) - stimon - tint (ms) -} - -INCLUDE "stimonoff.hoc" - -BREAKPOINT { - : Main computation block - - : compute Q - Q = v * cm - - : update stimulation boolean - stimonoff() - - : integrate states - SOLVE states METHOD cnexp - - : compute Vmeff - if(stimon) {Vmeff = veff_on(Q)} else {Vmeff = veff_off(Q)} - - : compute ionic currents - iNa = gnabar * m * m * m * h * (Vmeff - ena) - iKd = gkdbar * n * n * n * n * (Vmeff - ek) - iM = gmbar * p * (Vmeff - ek) - iCa = gcabar * s * s * u * (Vmeff - eca) - iLeak = gleak * (Vmeff - eleak) -} - -DERIVATIVE states { - : Compute states derivatives - - : compute ON or OFF rate constants accordingly to stimon value - if(stimon) {interpolate_on(Q)} else {interpolate_off(Q)} - - : compute gating states derivatives - m' = alpha_m * (1 - m) - beta_m * m - h' = alpha_h * (1 - h) - beta_h * h - n' = alpha_n * (1 - n) - beta_n * n - p' = alpha_p * (1 - p) - beta_p * p - s' = alpha_s * (1 - s) - beta_s * s - u' = alpha_u * (1 - u) - beta_u * u -} - -UNITSOFF : turning off units checking for functions not working with standard SI units - -INITIAL { - : Initialize variables and parameters - - : compute Q - Q = v * cm - - : set initial states values - m = alpham_off(Q) / (alpham_off(Q) + betam_off(Q)) - h = alphah_off(Q) / (alphah_off(Q) + betah_off(Q)) - n = alphan_off(Q) / (alphan_off(Q) + betan_off(Q)) - p = alphap_off(Q) / (alphap_off(Q) + betap_off(Q)) - s = alphas_off(Q) / (alphas_off(Q) + betas_off(Q)) - u = alphau_off(Q) / (alphau_off(Q) + betau_off(Q)) - - : initialize tint, set stimulation state and correct PRF - tint = 0 - stimon = 1 - PRF = PRF / 2 : for some reason PRF must be halved in order to ensure proper on/off switch -} - -: Define function tables for variables to be interpolated during ON and OFF periods -FUNCTION_TABLE veff_on(x) (mV) -FUNCTION_TABLE veff_off(x) (mV) -FUNCTION_TABLE alpham_on(x) (/ms) -FUNCTION_TABLE alpham_off(x) (/ms) -FUNCTION_TABLE betam_on(x) (/ms) -FUNCTION_TABLE betam_off(x) (/ms) -FUNCTION_TABLE alphah_on(x) (/ms) -FUNCTION_TABLE alphah_off(x) (/ms) -FUNCTION_TABLE betah_on(x) (/ms) -FUNCTION_TABLE betah_off(x) (/ms) -FUNCTION_TABLE alphan_on(x) (/ms) -FUNCTION_TABLE alphan_off(x) (/ms) -FUNCTION_TABLE betan_on(x) (/ms) -FUNCTION_TABLE betan_off(x) (/ms) -FUNCTION_TABLE alphap_on(x) (/ms) -FUNCTION_TABLE alphap_off(x) (/ms) -FUNCTION_TABLE betap_on(x) (/ms) -FUNCTION_TABLE betap_off(x) (/ms) -FUNCTION_TABLE alphas_on(x) (/ms) -FUNCTION_TABLE alphas_off(x) (/ms) -FUNCTION_TABLE betas_on(x) (/ms) -FUNCTION_TABLE betas_off(x) (/ms) -FUNCTION_TABLE alphau_on(x) (/ms) -FUNCTION_TABLE alphau_off(x) (/ms) -FUNCTION_TABLE betau_on(x) (/ms) -FUNCTION_TABLE betau_off(x) (/ms) - -PROCEDURE interpolate_on(Q){ - : Interpolate rate constants durong US-ON periods from appropriate tables - - alpha_m = alpham_on(Q) - beta_m = betam_on(Q) - alpha_h = alphah_on(Q) - beta_h = betah_on(Q) - alpha_n = alphan_on(Q) - beta_n = betan_on(Q) - alpha_p = alphap_on(Q) - beta_p = betap_on(Q) - alpha_s = alphas_on(Q) - beta_s = betas_on(Q) - alpha_u = alphau_on(Q) - beta_u = betau_on(Q) -} - -PROCEDURE interpolate_off(Q){ - : Interpolate rate constants durong US-OFF periods from appropriate tables - - alpha_m = alpham_off(Q) - beta_m = betam_off(Q) - alpha_h = alphah_off(Q) - beta_h = betah_off(Q) - alpha_n = alphan_off(Q) - beta_n = betan_off(Q) - alpha_p = alphap_off(Q) - beta_p = betap_off(Q) - alpha_s = alphas_off(Q) - beta_s = betas_off(Q) - alpha_u = alphau_off(Q) - beta_u = betau_off(Q) -} - -UNITSON : turn back units checking \ No newline at end of file diff --git a/PointNICE/neurons/nmodl/RE.mod b/PointNICE/neurons/nmodl/RE.mod deleted file mode 100644 index 5def3ae..0000000 --- a/PointNICE/neurons/nmodl/RE.mod +++ /dev/null @@ -1,203 +0,0 @@ -TITLE RE neuron -: -: Equations governing the behavior of a thalamic reticular neuron -: during ultrasonic stimulation, according to the NICE model. -: -: Written by Theo Lemaire and Simon Narduzzi, EPFL, 2018 -: Contact: theo.lemaire@epfl.ch - - -INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)} - -UNITS { - (mA) = (milliamp) - (mV) = (millivolt) - (uF) = (microfarad) -} - -NEURON { - : Definition of NEURON mechanism - - : mechanism name and constituting currents - SUFFIX RE - NONSPECIFIC_CURRENT iNa - NONSPECIFIC_CURRENT iKd - NONSPECIFIC_CURRENT iCa - NONSPECIFIC_CURRENT iLeak - - : simulation parameters (stimulus and dt) - RANGE duration, PRF, DC, stimon - - : physiological variables - RANGE Q, Vmeff, gnabar, gkdbar, gcabar, gleak, eleak, ek, ena, eca -} - - -PARAMETER { - : Parameters set by the user upon initialization - - : simulation parameters (stimulus and dt) - duration (ms) - PRF (Hz) - DC - dt (ms) - - : membrane properties - cm = 1 (uF/cm2) - ena = 50 (mV) - eca = 120 (mV) - ek = -90 (mV) - eleak = -90.0 (mV) - gnabar = 0.2 (mho/cm2) - gkdbar = 0.02 (mho/cm2) - gcabar = 0.003 (mho/cm2) - gleak = 5e-5 (mho/cm2) -} - -STATE { - : Differential variables other than v, i.e. the ion channels gating states - - m : iNa activation gate - h : iNa inactivation gate - n : iKd activation gate - s : iCa activation gate - u : iCa inactivation gate -} - -ASSIGNED { - : Variables computed during the simulation and whose value can be retrieved - Q (nC/cm2) - Vmeff (mV) - v (mV) - iNa (mA/cm2) - iKd (mA/cm2) - iCa (mA/cm2) - iLeak (mA/cm2) - alpha_h (/ms) - beta_h (/ms) - alpha_m (/ms) - beta_m (/ms) - alpha_n (/ms) - beta_n (/ms) - alpha_s (/ms) - beta_s (/ms) - alpha_u (/ms) - beta_u (/ms) - stimon - tint (ms) -} - -INCLUDE "stimonoff.hoc" - -BREAKPOINT { - : Main computation block - - : compute Q - Q = v * cm - - : update stimulation boolean - stimonoff() - - : integrate states - SOLVE states METHOD cnexp - - : compute Vmeff - if(stimon) {Vmeff = veff_on(Q)} else {Vmeff = veff_off(Q)} - - : compute ionic currents - iNa = gnabar * m * m * m * h * (Vmeff - ena) - iKd = gkdbar * n * n * n * n * (Vmeff - ek) - iCa = gcabar * s * s * u * (Vmeff - eca) - iLeak = gleak * (Vmeff - eleak) -} - -DERIVATIVE states { - : Compute states derivatives - - : compute ON or OFF rate constants accordingly to stimon value - if(stimon) {interpolate_on(Q)} else {interpolate_off(Q)} - - : compute gating states derivatives - m' = alpha_m * (1 - m) - beta_m * m - h' = alpha_h * (1 - h) - beta_h * h - n' = alpha_n * (1 - n) - beta_n * n - s' = alpha_s * (1 - s) - beta_s * s - u' = alpha_u * (1 - u) - beta_u * u -} - -UNITSOFF : turning off units checking for functions not working with standard SI units - -INITIAL { - : Initialize variables and parameters - - : compute Q - Q = v * cm - - : set initial states values - m = alpham_off(Q) / (alpham_off(Q) + betam_off(Q)) - h = alphah_off(Q) / (alphah_off(Q) + betah_off(Q)) - n = alphan_off(Q) / (alphan_off(Q) + betan_off(Q)) - s = alphas_off(Q) / (alphas_off(Q) + betas_off(Q)) - u = alphau_off(Q) / (alphau_off(Q) + betau_off(Q)) - - : initialize tint, set stimulation state and correct PRF - tint = 0 - stimon = 1 - PRF = PRF / 2 : for some reason PRF must be halved in order to ensure proper on/off switch -} - -: Define function tables for variables to be interpolated during ON and OFF periods -FUNCTION_TABLE veff_on(x) (mV) -FUNCTION_TABLE veff_off(x) (mV) -FUNCTION_TABLE alpham_on(x) (/ms) -FUNCTION_TABLE alpham_off(x) (/ms) -FUNCTION_TABLE betam_on(x) (/ms) -FUNCTION_TABLE betam_off(x) (/ms) -FUNCTION_TABLE alphah_on(x) (/ms) -FUNCTION_TABLE alphah_off(x) (/ms) -FUNCTION_TABLE betah_on(x) (/ms) -FUNCTION_TABLE betah_off(x) (/ms) -FUNCTION_TABLE alphan_on(x) (/ms) -FUNCTION_TABLE alphan_off(x) (/ms) -FUNCTION_TABLE betan_on(x) (/ms) -FUNCTION_TABLE betan_off(x) (/ms) -FUNCTION_TABLE alphas_on(x) (/ms) -FUNCTION_TABLE alphas_off(x) (/ms) -FUNCTION_TABLE betas_on(x) (/ms) -FUNCTION_TABLE betas_off(x) (/ms) -FUNCTION_TABLE alphau_on(x) (/ms) -FUNCTION_TABLE alphau_off(x) (/ms) -FUNCTION_TABLE betau_on(x) (/ms) -FUNCTION_TABLE betau_off(x) (/ms) - -PROCEDURE interpolate_on(Q){ - : Interpolate rate constants durong US-ON periods from appropriate tables - - alpha_m = alpham_on(Q) - beta_m = betam_on(Q) - alpha_h = alphah_on(Q) - beta_h = betah_on(Q) - alpha_n = alphan_on(Q) - beta_n = betan_on(Q) - alpha_s = alphas_on(Q) - beta_s = betas_on(Q) - alpha_u = alphau_on(Q) - beta_u = betau_on(Q) -} - -PROCEDURE interpolate_off(Q){ - : Interpolate rate constants durong US-OFF periods from appropriate tables - - alpha_m = alpham_off(Q) - beta_m = betam_off(Q) - alpha_h = alphah_off(Q) - beta_h = betah_off(Q) - alpha_n = alphan_off(Q) - beta_n = betan_off(Q) - alpha_s = alphas_off(Q) - beta_s = betas_off(Q) - alpha_u = alphau_off(Q) - beta_u = betau_off(Q) -} - -UNITSON : turn back units checking \ No newline at end of file diff --git a/PointNICE/neurons/nmodl/RS.mod b/PointNICE/neurons/nmodl/RS.mod deleted file mode 100644 index a21b2b1..0000000 --- a/PointNICE/neurons/nmodl/RS.mod +++ /dev/null @@ -1,189 +0,0 @@ -TITLE RS neuron -: -: Equations governing the behavior of a regular spiking neuron -: during ultrasonic stimulation, according to the NICE model. -: -: Written by Theo Lemaire and Simon Narduzzi, EPFL, 2018 -: Contact: theo.lemaire@epfl.ch - - -INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)} - -UNITS { - (mA) = (milliamp) - (mV) = (millivolt) - (uF) = (microfarad) -} - -NEURON { - : Definition of NEURON mechanism - - : mechanism name and constituting currents - SUFFIX RS - NONSPECIFIC_CURRENT iNa - NONSPECIFIC_CURRENT iKd - NONSPECIFIC_CURRENT iM - NONSPECIFIC_CURRENT iLeak - - : simulation parameters (stimulus and dt) - RANGE duration, PRF, DC, stimon - - : physiological variables - RANGE Q, Vmeff, gnabar, gkdbar, gmbar, gleak, eleak, ek, ena -} - - -PARAMETER { - : Parameters set by the user upon initialization - - : simulation parameters (stimulus and dt) - duration (ms) - PRF (Hz) - DC - dt (ms) - - : membrane properties - cm = 1 (uF/cm2) - ena = 50 (mV) - ek = -90 (mV) - eleak = -70.3 (mV) - gnabar = 0.056 (mho/cm2) - gkdbar = 0.006 (mho/cm2) - gmbar = 7.5e-5 (mho/cm2) - gleak = 2.05e-5 (mho/cm2) -} - -STATE { - : Differential variables other than v, i.e. the ion channels gating states - - m : iNa activation gate - h : iNa inactivation gate - n : iKd activation gate - p : iM activation gate -} - -ASSIGNED { - : Variables computed during the simulation and whose value can be retrieved - Q (nC/cm2) - Vmeff (mV) - v (mV) - iNa (mA/cm2) - iKd (mA/cm2) - iM (mA/cm2) - iLeak (mA/cm2) - alpha_h (/ms) - beta_h (/ms) - alpha_m (/ms) - beta_m (/ms) - alpha_n (/ms) - beta_n (/ms) - alpha_p (/ms) - beta_p (/ms) - stimon - tint (ms) -} - -INCLUDE "stimonoff.hoc" - -BREAKPOINT { - : Main computation block - - : compute Q - Q = v * cm - - : update stimulation boolean - stimonoff() - - : integrate states - SOLVE states METHOD cnexp - - : compute Vmeff - if(stimon) {Vmeff = veff_on(Q)} else {Vmeff = veff_off(Q)} - - : compute ionic currents - iNa = gnabar * m * m * m * h * (Vmeff - ena) - iKd = gkdbar * n * n * n * n * (Vmeff - ek) - iM = gmbar * p * (Vmeff - ek) - iLeak = gleak * (Vmeff - eleak) -} - -DERIVATIVE states { - : Compute states derivatives - - : compute ON or OFF rate constants accordingly to stimon value - if(stimon) {interpolate_on(Q)} else {interpolate_off(Q)} - - : compute gating states derivatives - m' = alpha_m * (1 - m) - beta_m * m - h' = alpha_h * (1 - h) - beta_h * h - n' = alpha_n * (1 - n) - beta_n * n - p' = alpha_p * (1 - p) - beta_p * p -} - -UNITSOFF : turning off units checking for functions not working with standard SI units - -INITIAL { - : Initialize variables and parameters - - : compute Q - Q = v * cm - - : set initial states values - m = alpham_off(Q) / (alpham_off(Q) + betam_off(Q)) - h = alphah_off(Q) / (alphah_off(Q) + betah_off(Q)) - n = alphan_off(Q) / (alphan_off(Q) + betan_off(Q)) - p = alphap_off(Q) / (alphap_off(Q) + betap_off(Q)) - - : initialize tint, set stimulation state and correct PRF - tint = 0 - stimon = 1 - PRF = PRF / 2 : for some reason PRF must be halved in order to ensure proper on/off switch -} - -: Define function tables for variables to be interpolated during ON and OFF periods -FUNCTION_TABLE veff_on(x) (mV) -FUNCTION_TABLE veff_off(x) (mV) -FUNCTION_TABLE alpham_on(x) (/ms) -FUNCTION_TABLE alpham_off(x) (/ms) -FUNCTION_TABLE betam_on(x) (/ms) -FUNCTION_TABLE betam_off(x) (/ms) -FUNCTION_TABLE alphah_on(x) (/ms) -FUNCTION_TABLE alphah_off(x) (/ms) -FUNCTION_TABLE betah_on(x) (/ms) -FUNCTION_TABLE betah_off(x) (/ms) -FUNCTION_TABLE alphan_on(x) (/ms) -FUNCTION_TABLE alphan_off(x) (/ms) -FUNCTION_TABLE betan_on(x) (/ms) -FUNCTION_TABLE betan_off(x) (/ms) -FUNCTION_TABLE alphap_on(x) (/ms) -FUNCTION_TABLE alphap_off(x) (/ms) -FUNCTION_TABLE betap_on(x) (/ms) -FUNCTION_TABLE betap_off(x) (/ms) - -PROCEDURE interpolate_on(Q){ - : Interpolate rate constants durong US-ON periods from appropriate tables - - alpha_m = alpham_on(Q) - beta_m = betam_on(Q) - alpha_h = alphah_on(Q) - beta_h = betah_on(Q) - alpha_n = alphan_on(Q) - beta_n = betan_on(Q) - alpha_p = alphap_on(Q) - beta_p = betap_on(Q) -} - -PROCEDURE interpolate_off(Q){ - : Interpolate rate constants durong US-OFF periods from appropriate tables - - alpha_m = alpham_off(Q) - beta_m = betam_off(Q) - alpha_h = alphah_off(Q) - beta_h = betah_off(Q) - alpha_n = alphan_off(Q) - beta_n = betan_off(Q) - alpha_p = alphap_off(Q) - beta_p = betap_off(Q) -} - -UNITSON : turn back units checking \ No newline at end of file diff --git a/PointNICE/neurons/nmodl/TC.mod b/PointNICE/neurons/nmodl/TC.mod deleted file mode 100644 index 8383b6d..0000000 --- a/PointNICE/neurons/nmodl/TC.mod +++ /dev/null @@ -1,283 +0,0 @@ -TITLE TC neuron -: -: Equations governing the behavior of a thalamo-cortical neuron -: during ultrasonic stimulation, according to the NICE model. -: -: Written by Theo Lemaire and Simon Narduzzi, EPFL, 2018 -: Contact: theo.lemaire@epfl.ch - - -INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)} - -UNITS { - :(mA) = (milliamp) - :(mV) = (millivolt) - :(uF) = (microfarad) - :(nC) = (nanocoulomb) - - (molar) = (1/liter) : moles do not appear in units - (M) = (molar) - (mM) = (millimolar) - (um) = (micron) - (mA) = (milliamp) - (msM) = (ms mM) -} - -NEURON { - : Definition of NEURON mechanism - - : mechanism name and constituting currents - SUFFIX TC - NONSPECIFIC_CURRENT iNa - NONSPECIFIC_CURRENT iKd - NONSPECIFIC_CURRENT iKl - NONSPECIFIC_CURRENT iH - NONSPECIFIC_CURRENT iCa - NONSPECIFIC_CURRENT iLeak - - : simulation parameters (stimulus and dt) - RANGE duration, PRF, DC, stimon - - : physiological variables - RANGE Q, Vmeff, gnabar, gkdbar, gkl, gcabar, gleak, ghbar, eleak, ek, ena, eca, eh - RANGE k1, k2, k3, k4, nca - RANGE depth, taur, camin -} - -CONSTANT { - FARADAY = 96494 (coul) : moles do not appear in units -} - - -PARAMETER { - : Parameters set by the user upon initialization - - : simulation parameters (stimulus and dt) - duration (ms) - PRF (Hz) - DC - dt (ms) - - : membrane properties - cm = 1 (uF/cm2) - ena = 50 (mV) - eca = 120 (mV) - ek = -90 (mV) - eh = -40 (mV) - eleak = -70 (mV) - gnabar = 0.09 (mho/cm2) - gkdbar = 0.01 (mho/cm2) - gkl = 1.38e-5 (mho/cm2) - gcabar = 0.002 (mho/cm2) - ghbar = 1.75e-5 (mho/cm2) - gleak = 1e-5 (mho/cm2) - - : iH Calcium dependence properties - k1 = 2.5e19 (1/M*M*M*M*ms) : CB protein Calcium-driven activation rate - k2 = 0.0004 (1/ms) : CB protein inactivation rate - k3 = 0.1 (1/ms) : CB protein iH channel binding rate - k4 = 0.001 (1/ms) : CB protein iH channel unbinding rate - nca = 4 : number of Calcium binding sites on CB protein - - : submembrane Calcium evolution properties - depth = 1e-7 (m) : depth of shell - taur = 5 (ms) : rate of calcium removal - camin = 5e-8 (M) : minimal intracellular Calcium concentration - -} - -STATE { - : Differential variables other than v - - m : iNa activation gate - h : iNa inactivation gate - n : iKd activation gate - s : iCa activation gate - u : iCa inactivation gate - C1 : iH channel closed state - O1 : iH channel open state - P0 : proportion of unbound CB protein - C_Ca (M) : submembrane Calcium concentration -} - -ASSIGNED { - : Variables computed during the simulation and whose value can be retrieved - Q (nC/cm2) - Vmeff (mV) - v (mV) - iNa (mA/cm2) - iKd (mA/cm2) - iKl (mA/cm2) - iCa (mA/cm2) - iH (mA/cm2) - iLeak (mA/cm2) - alpha_h (/ms) - beta_h (/ms) - alpha_m (/ms) - beta_m (/ms) - alpha_n (/ms) - beta_n (/ms) - alpha_s (/ms) - beta_s (/ms) - alpha_u (/ms) - beta_u (/ms) - alpha_o (/ms) - beta_o (/ms) - iCadrive (M/ms) - stimon - tint (ms) -} - -INCLUDE "stimonoff.hoc" - -BREAKPOINT { - : Main computation block - - : compute Q - Q = v * cm - - : update stimulation boolean - stimonoff() - - : integrate states - SOLVE states METHOD cnexp - - : check iH states and restrict them if needed - if(O1 < 0.) {O1 = 0.} - if(O1 > 1.) {O1 = 1.} - if(C1 < 0.) {C1 = 0.} - if(C1 > 1.) {C1 = 1.} - - : compute Vmeff - if(stimon) {Vmeff = veff_on(Q)} else {Vmeff = veff_off(Q)} - - : compute ionic currents - iNa = gnabar * m * m * m * h * (Vmeff - ena) - iKd = gkdbar * n * n * n * n * (Vmeff - ek) - iKl = gkl * (Vmeff - ek) - iCa = gcabar * s * s * u * (Vmeff - eca) - iLeak = gleak * (Vmeff - eleak) - iH = ghbar * (O1 + 2 * (1 - O1 - C1)) * (Vmeff - eh) -} - -UNITSOFF : turning off units checking for functions not working with standard SI units - -DERIVATIVE states { - : Compute states derivatives - - : compute ON or OFF rate constants accordingly to stimon value - if(stimon) {interpolate_on(Q)} else {interpolate_off(Q)} - - : compute gating states derivatives - m' = alpha_m * (1 - m) - beta_m * m - h' = alpha_h * (1 - h) - beta_h * h - n' = alpha_n * (1 - n) - beta_n * n - s' = alpha_s * (1 - s) - beta_s * s - u' = alpha_u * (1 - u) - beta_u * u - - : compute derivatives of variables for the kinetics system of Ih - iCadrive = -1e-5 * iCa / (2 * FARADAY * depth) - C_Ca' = (camin - C_Ca) / taur + iCadrive - P0' = k2 * (1 - P0) - k1 * P0 * capow(C_Ca) - C1' = beta_o * O1 - alpha_o * C1 - O1' = alpha_o * C1 - beta_o * O1 - k3 * O1 * (1 - P0) + k4 * (1 - O1 - C1) -} - -INITIAL { - : Initialize variables and parameters - - : compute Q - Q = v * cm - - : set initial states values - m = alpham_off(Q) / (alpham_off(Q) + betam_off(Q)) - h = alphah_off(Q) / (alphah_off(Q) + betah_off(Q)) - n = alphan_off(Q) / (alphan_off(Q) + betan_off(Q)) - s = alphas_off(Q) / (alphas_off(Q) + betas_off(Q)) - u = alphau_off(Q) / (alphau_off(Q) + betau_off(Q)) - - : compute steady-state Calcium concentration - iCa = gcabar * s * s * u * (veff_off(Q) - eca) - iCadrive = -1e-5 * iCa / (2 * FARADAY * depth) - C_Ca = camin + taur * iCadrive - - : compute steady values for the kinetics system of Ih - P0 = k2 / (k2 + k1 * C_Ca^nca) - O1 = k4 / (k3 * (1 - P0) + k4 * (1 + betao_off(Q) / alphao_off(Q))) - C1 = betao_off(Q) / alphao_off(Q) * O1 - - : initialize tint, set stimulation state and correct PRF - tint = 0 - stimon = 1 - PRF = PRF / 2 : for some reason PRF must be halved in order to ensure proper on/off switch -} - -: Define function tables for variables to be interpolated during ON and OFF periods -FUNCTION_TABLE veff_on(x) (mV) -FUNCTION_TABLE veff_off(x) (mV) -FUNCTION_TABLE alpham_on(x) (/ms) -FUNCTION_TABLE alpham_off(x) (/ms) -FUNCTION_TABLE betam_on(x) (/ms) -FUNCTION_TABLE betam_off(x) (/ms) -FUNCTION_TABLE alphah_on(x) (/ms) -FUNCTION_TABLE alphah_off(x) (/ms) -FUNCTION_TABLE betah_on(x) (/ms) -FUNCTION_TABLE betah_off(x) (/ms) -FUNCTION_TABLE alphan_on(x) (/ms) -FUNCTION_TABLE alphan_off(x) (/ms) -FUNCTION_TABLE betan_on(x) (/ms) -FUNCTION_TABLE betan_off(x) (/ms) -FUNCTION_TABLE alphas_on(x) (/ms) -FUNCTION_TABLE alphas_off(x) (/ms) -FUNCTION_TABLE betas_on(x) (/ms) -FUNCTION_TABLE betas_off(x) (/ms) -FUNCTION_TABLE alphau_on(x) (/ms) -FUNCTION_TABLE alphau_off(x) (/ms) -FUNCTION_TABLE betau_on(x) (/ms) -FUNCTION_TABLE betau_off(x) (/ms) -FUNCTION_TABLE alphao_on(x) (/ms) -FUNCTION_TABLE alphao_off(x) (/ms) -FUNCTION_TABLE betao_on(x) (/ms) -FUNCTION_TABLE betao_off(x) (/ms) - -PROCEDURE interpolate_on(Q){ - : Interpolate rate constants durong US-ON periods from appropriate tables - - alpha_m = alpham_on(Q) - beta_m = betam_on(Q) - alpha_h = alphah_on(Q) - beta_h = betah_on(Q) - alpha_n = alphan_on(Q) - beta_n = betan_on(Q) - alpha_s = alphas_on(Q) - beta_s = betas_on(Q) - alpha_u = alphau_on(Q) - beta_u = betau_on(Q) - alpha_o = alphao_on(Q) - beta_o = betao_on(Q) -} - -PROCEDURE interpolate_off(Q){ - : Interpolate rate constants durong US-OFF periods from appropriate tables - - alpha_m = alpham_off(Q) - beta_m = betam_off(Q) - alpha_h = alphah_off(Q) - beta_h = betah_off(Q) - alpha_n = alphan_off(Q) - beta_n = betan_off(Q) - alpha_s = alphas_off(Q) - beta_s = betas_off(Q) - alpha_u = alphau_off(Q) - beta_u = betau_off(Q) - alpha_o = alphao_off(Q) - beta_o = betao_off(Q) -} - - -FUNCTION capow(C_Ca (M)) (M*M*M*M) { - capow = C_Ca^nca -} - - -UNITSON : turn back units checking \ No newline at end of file diff --git a/PointNICE/neurons/nmodl/stimonoff.hoc b/PointNICE/neurons/nmodl/stimonoff.hoc deleted file mode 100644 index 477cfe9..0000000 --- a/PointNICE/neurons/nmodl/stimonoff.hoc +++ /dev/null @@ -1,22 +0,0 @@ -: Procedure handling the automatic switch between the ON and OFF systems during the -: numerical integration of ASTIM simulations. -: -: Written by Theo Lemaire, EPFL, 2018 -: Contact: theo.lemaire@epfl.ch - - -PROCEDURE stimonoff(){ - : Switch integration to ON or OFF system according to stimulus - - if(t < duration - dt){ - if(tint >= 1000 / PRF){ : OFF-> ON at pulse onset - stimon = 1 - tint = 0 - }else if(tint > DC * 1000 / PRF && stimon == 1){ : ON -> OFF at pulse offset - stimon = 0 - } - tint = tint + dt : increment by dt - }else if(stimon == 1){ : ON -> OFF at stimulus offset - stimon = 0 - } -} \ No newline at end of file diff --git a/PointNICE/solvers/SolverUS.py b/PointNICE/solvers/SolverUS.py index 3deab25..5ba6031 100644 --- a/PointNICE/solvers/SolverUS.py +++ b/PointNICE/solvers/SolverUS.py @@ -1,926 +1,747 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-09-29 16:16:19 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-07-20 18:21:35 +# @Last Modified time: 2018-08-21 14:19:43 import os import warnings import pickle import logging import progressbar as pb import numpy as np import scipy.integrate as integrate from scipy.interpolate import interp2d -from neuron import h from ..bls import BilayerSonophore from ..utils import * from ..constants import * from ..neurons import BaseMech # Get package logger logger = logging.getLogger('PointNICE') -# global list of paths already loaded by load_mechanisms -nrn_dll_loaded = [] - class SolverUS(BilayerSonophore): """ This class extends the BilayerSonophore class by adding a biophysical Hodgkin-Huxley model on top of the mechanical BLS model. """ def __init__(self, diameter, neuron, Fdrive, embedding_depth=0.0): """ Constructor of the class. :param diameter: in-plane diameter of the sonophore structure within the membrane (m) :param neuron: neuron object :param Fdrive: frequency of acoustic perturbation (Hz) :param embedding_depth: depth of the embedding tissue around the membrane (m) """ # Check validity of input parameters if not isinstance(neuron, BaseMech): raise InputError('Invalid neuron type: "{}" (must inherit from BaseMech class)' .format(neuron.name)) if not isinstance(Fdrive, float): raise InputError('Invalid US driving frequency (must be float typed)') if Fdrive < 0: raise InputError('Invalid US driving frequency: {} kHz (must be positive or null)' .format(Fdrive * 1e-3)) # TODO: check parameters dictionary (float type, mandatory members) # Initialize BLS object Cm0 = neuron.Cm0 Vm0 = neuron.Vm0 BilayerSonophore.__init__(self, diameter, Fdrive, Cm0, Cm0 * Vm0 * 1e-3, embedding_depth) logger.debug('US solver initialization with %s neuron', neuron.name) def eqHH(self, y, t, neuron, Cm): """ Compute the derivatives of the n-ODE HH system variables, based on a value of membrane capacitance. :param y: vector of HH system variables at time t :param t: specific instant in time (s) :param neuron: neuron object :param Cm: membrane capacitance (F/m2) :return: vector of HH system derivatives at time t """ # Split input vector explicitly Qm, *states = y # Compute membrane potential Vm = Qm / Cm * 1e3 # mV # Compute derivatives dQm = - neuron.currNet(Vm, states) * 1e-3 # A/m2 dstates = neuron.derStates(Vm, states) # Return derivatives vector return [dQm, *dstates] def eqHH2(self, t, y, neuron, Cm): return self.eqHH(y, t, neuron, Cm) def eqFull(self, y, t, neuron, Adrive, Fdrive, phi): """ Compute the derivatives of the (n+3) ODE full NBLS system variables. :param y: vector of state variables :param t: specific instant in time (s) :param neuron: neuron object :param Adrive: acoustic drive amplitude (Pa) :param Fdrive: acoustic drive frequency (Hz) :param phi: acoustic drive phase (rad) :return: vector of derivatives """ # Compute derivatives of mechanical and electrical systems dydt_mech = self.eqMech(y[:3], t, Adrive, Fdrive, y[3], phi) dydt_elec = self.eqHH(y[3:], t, neuron, self.Capct(y[1])) # return concatenated output return dydt_mech + dydt_elec def eqFull2(self, t, y, neuron, Adrive, Fdrive, phi): return self.eqFull(y, t, neuron, Adrive, Fdrive, phi) def eqHHeff(self, t, y, neuron, interp_data): """ Compute the derivatives of the n-ODE effective HH system variables, based on 1-dimensional linear interpolation of "effective" coefficients that summarize the system's behaviour over an acoustic cycle. :param t: specific instant in time (s) :param y: vector of HH system variables at time t :param neuron: neuron object :param interp_data: dictionary of 1D data points of "effective" coefficients over the charge domain, for specific frequency and amplitude values. :return: vector of effective system derivatives at time t """ # Split input vector explicitly Qm, *states = y # Compute charge and channel states variation Vm = np.interp(Qm, interp_data['Q'], interp_data['V']) # mV dQmdt = - neuron.currNet(Vm, states) * 1e-3 dstates = neuron.derStatesEff(Qm, states, interp_data) # Return derivatives vector return [dQmdt, *dstates] def __runClassic(self, neuron, Fdrive, Adrive, tstim, toffset, PRF, DC, phi=np.pi): """ Compute solutions of the system for a specific set of US stimulation parameters, using a classic integration scheme. The first iteration uses the quasi-steady simplification to compute the initiation of motion from a flat leaflet configuration. Afterwards, the ODE system is solved iteratively until completion. :param neuron: neuron object :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param tstim: duration of US stimulation (s) :param toffset: duration of the offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) :param phi: acoustic drive phase (rad) :return: 3-tuple with the time profile, the effective solution matrix and a state vector """ # Raise warnings as error warnings.filterwarnings('error') # Determine system time step Tdrive = 1 / Fdrive dt = Tdrive / NPC_FULL # if CW stimulus: divide integration during stimulus into 100 intervals if DC == 1.0: PRF = 100 / tstim # Compute vector sizes npulses = int(np.round(PRF * tstim)) Tpulse_on = DC / PRF Tpulse_off = (1 - DC) / PRF n_pulse_on = int(np.round(Tpulse_on / dt)) n_pulse_off = int(np.round(Tpulse_off / dt)) n_off = int(np.round(toffset / dt)) # Solve quasi-steady equation to compute first deflection value Z0 = 0.0 ng0 = self.ng0 Qm0 = self.Qm0 Pac1 = self.Pacoustic(dt, Adrive, Fdrive, phi) Z1 = self.balancedefQS(ng0, Qm0, Pac1) # Initialize global arrays states = np.array([1, 1]) t = np.array([0., dt]) y_membrane = np.array([[0., (Z1 - Z0) / dt], [Z0, Z1], [ng0, ng0], [Qm0, Qm0]]) y_channels = np.tile(neuron.states0, (2, 1)).T y = np.vstack((y_membrane, y_channels)) nvar = y.shape[0] # Initialize pulse time and states vectors t_pulse0 = np.linspace(0, Tpulse_on + Tpulse_off, n_pulse_on + n_pulse_off) states_pulse = np.concatenate((np.ones(n_pulse_on), np.zeros(n_pulse_off))) # Initialize progress bar if logger.getEffectiveLevel() <= logging.INFO: widgets = ['Running: ', pb.Percentage(), ' ', pb.Bar(), ' ', pb.ETA()] pbar = pb.ProgressBar(widgets=widgets, max_value=int(npulses * (toffset + tstim) / tstim)) pbar.start() # Loop through all pulse (ON and OFF) intervals for i in range(npulses): # Construct and initialize arrays t_pulse = t_pulse0 + t[-1] y_pulse = np.empty((nvar, n_pulse_on + n_pulse_off)) # Integrate ON system y_pulse[:, :n_pulse_on] = integrate.odeint(self.eqFull, y[:, -1], t_pulse[:n_pulse_on], args=(neuron, Adrive, Fdrive, phi)).T # Integrate OFF system if n_pulse_off > 0: y_pulse[:, n_pulse_on:] = integrate.odeint(self.eqFull, y_pulse[:, n_pulse_on - 1], t_pulse[n_pulse_on:], args=(neuron, 0.0, 0.0, 0.0)).T # Append pulse arrays to global arrays states = np.concatenate([states, states_pulse[1:]]) t = np.concatenate([t, t_pulse[1:]]) y = np.concatenate([y, y_pulse[:, 1:]], axis=1) # Update progress bar if logger.getEffectiveLevel() <= logging.INFO: pbar.update(i) # Integrate offset interval if n_off > 0: t_off = np.linspace(0, toffset, n_off) + t[-1] states_off = np.zeros(n_off) y_off = integrate.odeint(self.eqFull, y[:, -1], t_off, args=(neuron, 0.0, 0.0, 0.0)).T # Concatenate offset arrays to global arrays states = np.concatenate([states, states_off[1:]]) t = np.concatenate([t, t_off[1:]]) y = np.concatenate([y, y_off[:, 1:]], axis=1) # Terminate progress bar if logger.getEffectiveLevel() <= logging.INFO: pbar.finish() # Downsample arrays in time-domain accordgin to target temporal resolution ds_factor = int(np.round(CLASSIC_TARGET_DT / dt)) if ds_factor > 1: Fs = 1 / (dt * ds_factor) logger.info('Downsampling output arrays by factor %u (Fs = %.2f MHz)', ds_factor, Fs * 1e-6) t = t[::ds_factor] y = y[:, ::ds_factor] states = states[::ds_factor] # Compute membrane potential vector (in mV) Vm = y[3, :] / self.v_Capct(y[1, :]) * 1e3 # mV # Return output variables with Vm # return (t, y[1:, :], states) return (t, np.vstack([y[1:4, :], Vm, y[4:, :]]), states) def __runEffective(self, neuron, Fdrive, Adrive, tstim, toffset, PRF, DC, dt=DT_EFF): """ Compute solutions of the system for a specific set of US stimulation parameters, using charge-predicted "effective" coefficients to solve the HH equations at each step. :param neuron: neuron object :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param tstim: duration of US stimulation (s) :param toffset: duration of the offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) :param dt: integration time step (s) :return: 3-tuple with the time profile, the effective solution matrix and a state vector """ # Raise warnings as error warnings.filterwarnings('error') - # Check lookup file existence - lookup_file = '{}_lookups_a{:.1f}nm.pkl'.format(neuron.name, self.a * 1e9) - lookup_path = '{}/{}'.format(getLookupDir(), lookup_file) - if not os.path.isfile(lookup_path): - raise InputError('Missing lookup file: "{}"'.format(lookup_file)) - - # Load lookups dictionary - with open(lookup_path, 'rb') as fh: - lookups3D = pickle.load(fh) - - # Retrieve 1D inputs from lookups dictionary - freqs = lookups3D.pop('f') - amps = lookups3D.pop('A') - charges = lookups3D.pop('Q') + # Load appropriate 2D lookups + Aref, Qref, lookups2D = getLookups2D(neuron.name, self.a, Fdrive) - # Check that stimulation parameters are within lookup range + # Check that acoustic amplitude is within lookup range margin = 1e-9 # adding margin to compensate for eventual round error - frange = (freqs.min() - margin, freqs.max() + margin) - Arange = (amps.min() - margin, amps.max() + margin) - - if Fdrive < frange[0] or Fdrive > frange[1]: - raise InputError(('Invalid frequency: {:.2f} kHz (must be within ' + - '{:.1f} kHz - {:.1f} MHz lookup interval)') - .format(Fdrive * 1e-3, frange[0] * 1e-3, frange[1] * 1e-6)) + Arange = (Aref.min() - margin, Aref.max() + margin) if Adrive < Arange[0] or Adrive > Arange[1]: - raise InputError(('Invalid amplitude: {:.2f} kPa (must be within ' + - '{:.1f} - {:.1f} kPa lookup interval)') - .format(Adrive * 1e-3, Arange[0] * 1e-3, Arange[1] * 1e-3)) - - # Interpolate 3D lookups at US frequency - lookups2D = itrpLookupsFreq(lookups3D, freqs, Fdrive) + raise InputError('Invalid amplitude: {}Pa (must be within {}Pa - {} Pa lookup interval)' + .format(*si_format([Adrive, *Arange], precision=2, space=' '))) # Interpolate 2D lookups at US amplitude (along with "ng" at zero amplitude) - lookups1D = {key: np.squeeze(interp2d(amps, charges, lookups2D[key].T)(Adrive, charges)) + lookups1D = {key: np.squeeze(interp2d(Aref, Qref, lookups2D[key].T)(Adrive, Qref)) for key in lookups2D.keys()} - lookups1D['ng0'] = np.squeeze(interp2d(amps, charges, lookups2D['ng'].T)(0.0, charges)) + lookups1D['ng0'] = np.squeeze(interp2d(Aref, Qref, lookups2D['ng'].T)(0.0, Qref)) # Add reference charge vector to 1D lookup dictionary - lookups1D['Q'] = charges + lookups1D['Q'] = Qref # Initialize system solvers solver_on = integrate.ode(self.eqHHeff) solver_on.set_integrator('lsoda', nsteps=SOLVER_NSTEPS) solver_on.set_f_params(neuron, lookups1D) solver_off = integrate.ode(self.eqHH2) solver_off.set_integrator('lsoda', nsteps=SOLVER_NSTEPS) # if CW stimulus: change PRF to have exactly one integration interval during stimulus if DC == 1.0: PRF = 1 / tstim # Compute vector sizes npulses = int(np.round(PRF * tstim)) Tpulse_on = DC / PRF Tpulse_off = (1 - DC) / PRF # For high-PRF pulsed protocols: adapt time step to ensure minimal # number of samples during TON or TOFF dt_warning_msg = 'high-PRF protocol: lowering time step to %.2e s to properly integrate %s' for key, Tpulse in {'TON': Tpulse_on, 'TOFF': Tpulse_off}.items(): if Tpulse > 0 and Tpulse / dt < MIN_SAMPLES_PER_PULSE_INT: dt = Tpulse / MIN_SAMPLES_PER_PULSE_INT logger.warning(dt_warning_msg, dt, key) n_pulse_on = int(np.round(Tpulse_on / dt)) + 1 n_pulse_off = int(np.round(Tpulse_off / dt)) # Compute ofset size n_off = int(np.round(toffset / dt)) # Initialize global arrays states = np.array([1]) t = np.array([0.0]) y = np.atleast_2d(np.insert(neuron.states0, 0, self.Qm0)).T nvar = y.shape[0] Zeff = np.array([0.0]) ngeff = np.array([self.ng0]) # Initializing accurate pulse time vector t_pulse_on = np.linspace(0, Tpulse_on, n_pulse_on) t_pulse_off = np.linspace(dt, Tpulse_off, n_pulse_off) + Tpulse_on t_pulse0 = np.concatenate([t_pulse_on, t_pulse_off]) states_pulse = np.concatenate((np.ones(n_pulse_on), np.zeros(n_pulse_off))) # Loop through all pulse (ON and OFF) intervals for i in range(npulses): # Construct and initialize arrays t_pulse = t_pulse0 + t[-1] y_pulse = np.empty((nvar, n_pulse_on + n_pulse_off)) ngeff_pulse = np.empty(n_pulse_on + n_pulse_off) Zeff_pulse = np.empty(n_pulse_on + n_pulse_off) y_pulse[:, 0] = y[:, -1] ngeff_pulse[0] = ngeff[-1] Zeff_pulse[0] = Zeff[-1] # Initialize iterator k = 0 # Integrate ON system solver_on.set_initial_value(y_pulse[:, k], t_pulse[k]) while solver_on.successful() and k < n_pulse_on - 1: k += 1 solver_on.integrate(t_pulse[k]) y_pulse[:, k] = solver_on.y ngeff_pulse[k] = np.interp(y_pulse[0, k], lookups1D['Q'], lookups1D['ng']) # mole Zeff_pulse[k] = self.balancedefQS(ngeff_pulse[k], y_pulse[0, k]) # m # Integrate OFF system if n_pulse_off > 0: solver_off.set_initial_value(y_pulse[:, k], t_pulse[k]) solver_off.set_f_params(neuron, self.Capct(Zeff_pulse[k])) while solver_off.successful() and k < n_pulse_on + n_pulse_off - 1: k += 1 solver_off.integrate(t_pulse[k]) y_pulse[:, k] = solver_off.y ngeff_pulse[k] = np.interp(y_pulse[0, k], lookups1D['Q'], lookups1D['ng0']) # mole Zeff_pulse[k] = self.balancedefQS(ngeff_pulse[k], y_pulse[0, k]) # m solver_off.set_f_params(neuron, self.Capct(Zeff_pulse[k])) # Append pulse arrays to global arrays states = np.concatenate([states[:-1], states_pulse]) t = np.concatenate([t, t_pulse[1:]]) y = np.concatenate([y, y_pulse[:, 1:]], axis=1) Zeff = np.concatenate([Zeff, Zeff_pulse[1:]]) ngeff = np.concatenate([ngeff, ngeff_pulse[1:]]) # Integrate offset interval if n_off > 0: t_off = np.linspace(0, toffset, n_off) + t[-1] states_off = np.zeros(n_off) y_off = np.empty((nvar, n_off)) ngeff_off = np.empty(n_off) Zeff_off = np.empty(n_off) y_off[:, 0] = y[:, -1] ngeff_off[0] = ngeff[-1] Zeff_off[0] = Zeff[-1] solver_off.set_initial_value(y_off[:, 0], t_off[0]) solver_off.set_f_params(neuron, self.Capct(Zeff_pulse[-1])) k = 0 while solver_off.successful() and k < n_off - 1: k += 1 solver_off.integrate(t_off[k]) y_off[:, k] = solver_off.y ngeff_off[k] = np.interp(y_off[0, k], lookups1D['Q'], lookups1D['ng0']) # mole Zeff_off[k] = self.balancedefQS(ngeff_off[k], y_off[0, k]) # m solver_off.set_f_params(neuron, self.Capct(Zeff_off[k])) # Concatenate offset arrays to global arrays states = np.concatenate([states, states_off[1:]]) t = np.concatenate([t, t_off[1:]]) y = np.concatenate([y, y_off[:, 1:]], axis=1) Zeff = np.concatenate([Zeff, Zeff_off[1:]]) ngeff = np.concatenate([ngeff, ngeff_off[1:]]) # Compute membrane potential vector (in mV) Vm = np.zeros(states.size) Vm[states == 0] = y[0, states == 0] / self.v_Capct(Zeff[states == 0]) * 1e3 # mV Vm[states == 1] = np.interp(y[0, states == 1], lookups1D['Q'], lookups1D['V']) # mV # Add Zeff, ngeff and Vm to solution matrix y = np.vstack([Zeff, ngeff, y[0, :], Vm, y[1:, :]]) # return output variables return (t, y, states) - def __runNEURON(self, neuron, Fdrive, Adrive, tstim, toffset, PRF, DC, dt=DT_EFF): - - # Raise warnings as error - warnings.filterwarnings('error') - - # Define aliases for NMODL-protected variable names - neuron_aliases = {'O': 'O1', 'C': 'C1'} - - # Load mechanisms DLL file - nmodl_dir = getNmodlDir() # + '/{}'.format(neuron.name) - mod_file = nmodl_dir + '/{}.mod'.format(neuron.name) - dll_file = nmodl_dir + '/nrnmech.dll' - - # Raise warning if source MOD file is more recent than compiled DLL file - if not os.path.isfile(dll_file) or os.path.getmtime(mod_file) > os.path.getmtime(dll_file): - logger.warning('"%s.mod" file more recent than dll -> needs recompiling', neuron.name) - - # switch to NMODL directory, invoke mknrndll command and go back to original directory - # print('saving current working directory') - # cwd = os.getcwd() - # print('switching to NMODL directory') - # os.system('cd {}'.format(nmodl_dir)) - # print('invoking mknrndll command') - # exit_code = os.system('c:/nrn/mingw/bin/sh c:/nrn/lib/mknrndll.sh /c/nrn') - # print(exit_code) - # print('simulating return ') - # print('switching back to original directory') - # os.system('cd {}'.format(cwd)) - - # Load mechanisms DLL if not already loaded - if dll_file not in nrn_dll_loaded: - h.nrn_load_dll(dll_file) - nrn_dll_loaded.append(dll_file) - - # Check lookup file existence - lookup_file = '{}_lookups_a{:.1f}nm.pkl'.format(neuron.name, self.a * 1e9) - lookup_path = '{}/{}'.format(getLookupDir(), lookup_file) - if not os.path.isfile(lookup_path): - raise InputError('Missing lookup file: "{}"'.format(lookup_file)) - - # Load lookups dictionary - with open(lookup_path, 'rb') as fh: - lookups3D = pickle.load(fh) - - # Retrieve 1D inputs from lookups dictionary - freqs = lookups3D.pop('f') - amps = lookups3D.pop('A') - charges = lookups3D.pop('Q') - - # Check that stimulation parameters are within lookup range - margin = 1e-9 # adding margin to compensate for eventual round error - frange = (freqs.min() - margin, freqs.max() + margin) - Arange = (amps.min() - margin, amps.max() + margin) - - if Fdrive < frange[0] or Fdrive > frange[1]: - raise InputError(('Invalid frequency: {:.2f} kHz (must be within ' + - '{:.1f} kHz - {:.1f} MHz lookup interval)') - .format(Fdrive * 1e-3, frange[0] * 1e-3, frange[1] * 1e-6)) - if Adrive < Arange[0] or Adrive > Arange[1]: - raise InputError(('Invalid amplitude: {:.2f} kPa (must be within ' + - '{:.1f} - {:.1f} kPa lookup interval)') - .format(Adrive * 1e-3, Arange[0] * 1e-3, Arange[1] * 1e-3)) - - # Interpolate 3D lookups at US frequency - lookups2D = itrpLookupsFreq(lookups3D, freqs, Fdrive) - - # Interpolate 2D lookups at US amplitude and zero - lookupsA = {key: np.squeeze(interp2d(amps, charges, lookups2D[key].T)(Adrive, charges)) - for key in lookups2D.keys()} - lookups0 = {key: np.squeeze(interp2d(amps, charges, lookups2D[key].T)(0.0, charges)) - for key in lookups2D.keys()} - - # Rescale rate constants to ms-1 - for k in lookupsA.keys(): - if 'alpha' in k or 'beta' in k: - lookupsA[k] *= 1e-3 - lookups0[k] *= 1e-3 - - # Rename V lookup - lookupsA['veff'] = lookupsA.pop('V') - lookups0['veff'] = lookups0.pop('V') - - # Translate 1D lookups to h-vectors - qvec = h.Vector(charges * 1e5) - lookupsA_hoc = {k: h.Vector(lookupsA[k]) for k in lookupsA.keys()} - lookups0_hoc = {k: h.Vector(lookups0[k]) for k in lookups0.keys()} - - # Assign h-vectors lookups to mechanism tables - add_table_on = "h.table_{0}_on_{1}(lookupsA_hoc['{0}']._ref_x[0], qvec.size(), qvec._ref_x[0])" - add_table_off = "h.table_{0}_off_{1}(lookups0_hoc['{0}']._ref_x[0], qvec.size(), qvec._ref_x[0])" - eval(add_table_on.format('veff', neuron.name)) - eval(add_table_off.format('veff', neuron.name)) - for gate in neuron.getGates(): - gate = gate.lower() - eval(add_table_on.format('alpha{}'.format(gate), neuron.name)) - eval(add_table_off.format('alpha{}'.format(gate), neuron.name)) - eval(add_table_on.format('beta{}'.format(gate), neuron.name)) - eval(add_table_off.format('beta{}'.format(gate), neuron.name)) - - # Create point-neuron compartment - pointneuron = h.Section(name='pointneuron') - pointneuron.nseg = 1 - pointneuron.Ra = 1 # dummy axoplasmic resistance - pointneuron.diam = 1 # unit diameter - pointneuron.L = 1 # unit length - - # Add mechanisms - pointneuron.insert(neuron.name) - - # Set parameters of acoustic stimulus - setattr(pointneuron, 'duration_{}'.format(neuron.name), tstim * 1e3) - if DC < 1: - setattr(pointneuron, 'PRF_{}'.format(neuron.name), PRF) - else: - setattr(pointneuron, 'PRF_{}'.format(neuron.name), tstim * 1e3) - setattr(pointneuron, 'DC_{}'.format(neuron.name), DC) - - # Record evolution of time and output variables - t = h.Vector() - states = h.Vector() - Qm = h.Vector() - Vm = h.Vector() - t.record(h._ref_t) - Qm.record(getattr(pointneuron(0.5), '_ref_Q_{}'.format(neuron.name))) - Vm.record(getattr(pointneuron(0.5), '_ref_Vmeff_{}'.format(neuron.name))) - states.record(getattr(pointneuron(0.5), '_ref_stimon_{}'.format(neuron.name))) - probes = {} - for state in neuron.states_names: - probes[state] = h.Vector() - if state in neuron_aliases.keys(): - alias = neuron_aliases[state] - else: - alias = state - probes[state].record(getattr(pointneuron(0.5), '_ref_{}_{}'.format(alias, neuron.name))) - - # Run simulation - h.t = 0 - h.finitialize(neuron.Vm0) - tstop = (tstim + toffset) * 1e3 - h.dt = dt * 1e3 - while h.t < tstop: - h.fadvance() - - t = np.array(t.to_python()) * 1e-3 # s - Qm = np.array(Qm.to_python()) * 1e-5 # C/m2 - Vm = np.array(Vm.to_python()) # mV - channels = [np.array(probes[state].to_python()) for state in neuron.states_names] - y = np.array([np.zeros(t.size), np.zeros(t.size), Qm, Vm, *channels]) - states = np.array(states.to_python()) - - return (t, y, states) - - def __runHybrid(self, neuron, Fdrive, Adrive, tstim, toffset, phi=np.pi): """ Compute solutions of the system for a specific set of US stimulation parameters, using a hybrid integration scheme. The first iteration uses the quasi-steady simplification to compute the initiation of motion from a flat leaflet configuration. Afterwards, the NBLS ODE system is solved iteratively for "slices" of N-microseconds, in a 2-steps scheme: - First, the full (n+3) ODE system is integrated for a few acoustic cycles until Z and ng reach a stable periodic solution (limit cycle) - Second, the signals of the 3 mechanical variables over the last acoustic period are selected and resampled to a far lower sampling rate - Third, the HH n-ODE system is integrated for the remaining time of the slice, using periodic expansion of the mechanical signals to precompute the values of capacitance. :param neuron: neuron object :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param tstim: duration of US stimulation (s) :param toffset: duration of the offset (s) :param phi: acoustic drive phase (rad) :return: 3-tuple with the time profile, the solution matrix and a state vector .. warning:: This method cannot handle pulsed stimuli """ # Raise warnings as error warnings.filterwarnings('error') # Initialize full and HH systems solvers solver_full = integrate.ode(self.eqFull2) solver_full.set_f_params(neuron, Adrive, Fdrive, phi) solver_full.set_integrator('lsoda', nsteps=SOLVER_NSTEPS) solver_hh = integrate.ode(self.eqHH2) solver_hh.set_integrator('dop853', nsteps=SOLVER_NSTEPS, atol=1e-12) # Determine full and HH systems time steps Tdrive = 1 / Fdrive dt_full = Tdrive / NPC_FULL dt_hh = Tdrive / NPC_HH n_full_per_hh = int(NPC_FULL / NPC_HH) t_full_cycle = np.linspace(0, Tdrive - dt_full, NPC_FULL) t_hh_cycle = np.linspace(0, Tdrive - dt_hh, NPC_HH) # Determine number of samples in prediction vectors npc_pred = NPC_FULL - n_full_per_hh + 1 # Solve quasi-steady equation to compute first deflection value Z0 = 0.0 ng0 = self.ng0 Qm0 = self.Qm0 Pac1 = self.Pacoustic(dt_full, Adrive, Fdrive, phi) Z1 = self.balancedefQS(ng0, Qm0, Pac1) # Initialize global arrays states = np.array([1, 1]) t = np.array([0., dt_full]) y_membrane = np.array([[0., (Z1 - Z0) / dt_full], [Z0, Z1], [ng0, ng0], [Qm0, Qm0]]) y_channels = np.tile(neuron.states0, (2, 1)).T y = np.vstack((y_membrane, y_channels)) nvar = y.shape[0] # Initialize progress bar if logger.getEffectiveLevel() == logging.DEBUG: widgets = ['Running: ', pb.Percentage(), ' ', pb.Bar(), ' ', pb.ETA()] pbar = pb.ProgressBar(widgets=widgets, max_value=1000) pbar.start() # For each hybrid integration interval irep = 0 sim_error = False while not sim_error and t[-1] < tstim + toffset: # Integrate full system for a few acoustic cycles until stabilization periodic_conv = False j = 0 ng_last = None Z_last = None while not sim_error and not periodic_conv: if t[-1] > tstim: solver_full.set_f_params(neuron, 0.0, 0.0, 0.0) t_full = t_full_cycle + t[-1] + dt_full y_full = np.empty((nvar, NPC_FULL)) y0_full = y[:, -1] solver_full.set_initial_value(y0_full, t[-1]) k = 0 while solver_full.successful() and k <= NPC_FULL - 1: solver_full.integrate(t_full[k]) y_full[:, k] = solver_full.y k += 1 # Compare Z and ng signals over the last 2 acoustic periods if j > 0 and rmse(Z_last, y_full[1, :]) < Z_ERR_MAX \ and rmse(ng_last, y_full[2, :]) < NG_ERR_MAX: periodic_conv = True # Update last vectors for next comparison Z_last = y_full[1, :] ng_last = y_full[2, :] # Concatenate time and solutions to global vectors states = np.concatenate([states, np.ones(NPC_FULL)], axis=0) t = np.concatenate([t, t_full], axis=0) y = np.concatenate([y, y_full], axis=1) # Increment loop index j += 1 # Retrieve last period of the 3 mechanical variables to propagate in HH system t_last = t[-npc_pred:] mech_last = y[0:3, -npc_pred:] # Downsample signals to specified HH system time step (_, mech_pred) = DownSample(t_last, mech_last, NPC_HH) # Integrate HH system until certain dQ or dT is reached Q0 = y[3, -1] dQ = 0.0 t0_interval = t[-1] dt_interval = 0.0 j = 0 if t[-1] < tstim: tlim = tstim else: tlim = tstim + toffset while (not sim_error and t[-1] < tlim and (np.abs(dQ) < DQ_UPDATE or dt_interval < DT_UPDATE)): t_hh = t_hh_cycle + t[-1] + dt_hh y_hh = np.empty((nvar - 3, NPC_HH)) y0_hh = y[3:, -1] solver_hh.set_initial_value(y0_hh, t[-1]) k = 0 while solver_hh.successful() and k <= NPC_HH - 1: solver_hh.set_f_params(neuron, self.Capct(mech_pred[1, k])) solver_hh.integrate(t_hh[k]) y_hh[:, k] = solver_hh.y k += 1 # Concatenate time and solutions to global vectors states = np.concatenate([states, np.zeros(NPC_HH)], axis=0) t = np.concatenate([t, t_hh], axis=0) y = np.concatenate([y, np.concatenate([mech_pred, y_hh], axis=0)], axis=1) # Compute charge variation from interval beginning dQ = y[3, -1] - Q0 dt_interval = t[-1] - t0_interval # Increment loop index j += 1 # Update progress bar if logger.getEffectiveLevel() == logging.DEBUG: pbar.update(int(1000 * (t[-1] / (tstim + toffset)))) irep += 1 # Terminate progress bar if logger.getEffectiveLevel() == logging.DEBUG: pbar.finish() # Compute membrane potential vector (in mV) Vm = y[3, :] / self.v_Capct(y[1, :]) * 1e3 # mV # Return output variables with Vm # return (t, y[1:, :], states) return (t, np.vstack([y[1:4, :], Vm, y[4:, :]]), states) def run(self, neuron, Fdrive, Adrive, tstim, toffset, PRF=None, DC=1.0, sim_type='effective'): """ Run simulation of the system for a specific set of US stimulation parameters. :param neuron: neuron object :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param tstim: duration of US stimulation (s) :param toffset: duration of the offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) :param sim_type: selected integration method :return: 3-tuple with the time profile, the solution matrix and a state vector """ # Check validity of simulation type if sim_type not in ('classic', 'effective', 'NEURON', 'hybrid'): raise InputError('Invalid integration method: "{}"'.format(sim_type)) # Check validity of stimulation parameters if not isinstance(neuron, BaseMech): raise InputError('Invalid neuron type: "{}" (must inherit from BaseMech class)' .format(neuron.name)) if not all(isinstance(param, float) for param in [Fdrive, Adrive, tstim, toffset, DC]): raise InputError('Invalid stimulation parameters (must be float typed)') if Fdrive <= 0: raise InputError('Invalid US driving frequency: {} kHz (must be strictly positive)' .format(Fdrive * 1e-3)) if Adrive < 0: raise InputError('Invalid US pressure amplitude: {} kPa (must be positive or null)' .format(Adrive * 1e-3)) if tstim <= 0: raise InputError('Invalid stimulus duration: {} ms (must be strictly positive)' .format(tstim * 1e3)) if toffset < 0: raise InputError('Invalid stimulus offset: {} ms (must be positive or null)' .format(toffset * 1e3)) if DC <= 0.0 or DC > 1.0: raise InputError('Invalid duty cycle: {} (must be within ]0; 1])'.format(DC)) if DC < 1.0: if not isinstance(PRF, float): raise InputError('Invalid PRF value (must be float typed)') if PRF is None: raise InputError('Missing PRF value (must be provided when DC < 1)') if PRF < 1 / tstim: raise InputError('Invalid PRF: {} Hz (PR interval exceeds stimulus duration' .format(PRF)) if PRF >= Fdrive: raise InputError('Invalid PRF: {} Hz (must be smaller than driving frequency)' .format(PRF)) # Call appropriate simulation function if sim_type == 'classic': return self.__runClassic(neuron, Fdrive, Adrive, tstim, toffset, PRF, DC) elif sim_type == 'effective': return self.__runEffective(neuron, Fdrive, Adrive, tstim, toffset, PRF, DC) elif sim_type == 'NEURON': return self.__runNEURON(neuron, Fdrive, Adrive, tstim, toffset, PRF, DC) elif sim_type == 'hybrid': if DC < 1.0: raise InputError('Pulsed protocol incompatible with hybrid integration method') return self.__runHybrid(neuron, Fdrive, Adrive, tstim, toffset) def findRheobaseAmps(self, neuron, Fdrive, DCs, Vthr, curr='net'): ''' Find the rheobase amplitudes (i.e. threshold acoustic amplitudes of infinite duration that would result in excitation) of a specific neuron for various stimulation duty cycles. :param neuron: neuron object :param Fdrive: acoustic drive frequency (Hz) :param DCs: duty cycles vector (-) :param Vthr: threshold membrane potential above which the neuron necessarily fires (mV) :return: rheobase amplitudes vector (Pa) ''' # Check lookup file existence lookup_file = '{}_lookups_a{:.1f}nm.pkl'.format(neuron.name, self.a * 1e9) lookup_path = '{}/{}'.format(getLookupDir(), lookup_file) if not os.path.isfile(lookup_path): raise InputError('Missing lookup file: "{}"'.format(lookup_file)) # Load lookups dictionary with open(lookup_path, 'rb') as fh: lookups3D = pickle.load(fh) # Retrieve 1D inputs from lookups dictionary freqs = lookups3D.pop('f') amps = lookups3D.pop('A') charges = lookups3D.pop('Q') # Check that stimulation parameters are within lookup range margin = 1e-9 # adding margin to compensate for eventual round error frange = (freqs.min() - margin, freqs.max() + margin) if Fdrive < frange[0] or Fdrive > frange[1]: raise InputError(('Invalid frequency: {:.2f} kHz (must be within ' + '{:.1f} kHz - {:.1f} MHz lookup interval)') .format(Fdrive * 1e-3, frange[0] * 1e-3, frange[1] * 1e-6)) # Interpolate 3D lookpus at given frequency and threshold charge lookups2D = itrpLookupsFreq(lookups3D, freqs, Fdrive) Qthr = neuron.Cm0 * Vthr * 1e-3 # C/m2 lookups1D = {key: np.squeeze(interp2d(amps, charges, lookups2D[key].T)(amps, Qthr)) for key in lookups2D.keys()} # Remove unnecessary items ot get ON rates and effective potential at threshold charge rates_on = lookups1D rates_on.pop('ng') Vm_on = rates_on.pop('V') # Compute neuron OFF rates at threshold potential rates_off = neuron.getRates(Vthr) # Compute rheobase amplitudes rheboase_amps = np.empty(DCs.size) for i, DC in enumerate(DCs): sstates_pulse = np.empty((len(neuron.states_names), amps.size)) for j, x in enumerate(neuron.states_names): # If channel state, compute pulse-average steady-state values if x in neuron.getGates(): x = x.lower() alpha_str, beta_str = ['{}{}'.format(s, x) for s in ['alpha', 'beta']] alphax_pulse = rates_on[alpha_str] * DC + rates_off[alpha_str] * (1 - DC) betax_pulse = rates_on[beta_str] * DC + rates_off[beta_str] * (1 - DC) sstates_pulse[j, :] = alphax_pulse / (alphax_pulse + betax_pulse) # Otherwise assume the state has reached a steady-state value for Vthr else: sstates_pulse[j, :] = np.ones(amps.size) * neuron.steadyStates(Vthr)[j] # Compute the pulse average net (or leakage) current along the amplitude space if curr == 'net': iNet_on = neuron.currNet(Vm_on, sstates_pulse) iNet_off = neuron.currNet(Vthr, sstates_pulse) elif curr == 'leak': iNet_on = neuron.currL(Vm_on) iNet_off = neuron.currL(Vthr) iNet_avg = iNet_on * DC + iNet_off * (1 - DC) # Find the threshold amplitude that cancels the pulse average net current rheboase_amps[i] = np.interp(0, -iNet_avg, amps, left=0., right=np.nan) inan = np.where(np.isnan(rheboase_amps))[0] if inan.size > 0: if inan.size == rheboase_amps.size: logger.error('No rheobase amplitudes within [%s - %sPa] for the provided duty cycles', *si_format((amps.min(), amps.max()))) else: minDC = DCs[inan.max() + 1] logger.warning('No rheobase amplitudes within [%s - %sPa] below %.1f%% duty cycle', *si_format((amps.min(), amps.max())), minDC * 1e2) return rheboase_amps diff --git a/PointNICE/utils.py b/PointNICE/utils.py index 875257f..36bc89b 100644 --- a/PointNICE/utils.py +++ b/PointNICE/utils.py @@ -1,472 +1,511 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-09-19 22:30:46 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-07-20 17:22:24 +# @Last Modified time: 2018-08-21 14:08:06 """ Definition of generic utility functions used in other modules """ from enum import Enum import operator import os import pickle import tkinter as tk from tkinter import filedialog import inspect import json import yaml from openpyxl import load_workbook import numpy as np import colorlog from . import neurons def setLogger(): log_formatter = colorlog.ColoredFormatter( '%(log_color)s %(asctime)s %(message)s', datefmt='%d/%m/%Y %H:%M:%S:', reset=True, log_colors={ 'DEBUG': 'green', 'INFO': 'white', 'WARNING': 'yellow', 'ERROR': 'red', 'CRITICAL': 'red,bg_white', }, style='%' ) log_handler = colorlog.StreamHandler() log_handler.setFormatter(log_formatter) color_logger = colorlog.getLogger('PointNICE') color_logger.addHandler(log_handler) return color_logger # Get package logger logger = setLogger() class InputError(Exception): ''' Exception raised for errors in the input. ''' pass class PmCompMethod(Enum): """ Enum: types of computation method for the intermolecular pressure """ direct = 1 predict = 2 def loadYAML(filepath): """ Load a dictionary of parameters from an external YAML file. :param filepath: full path to the YAML file :return: parameters dictionary """ _, filename = os.path.split(filepath) logger.debug('Loading parameters from "%s"', filename) with open(filepath, 'r') as f: stream = f.read() params = yaml.load(stream) return ParseNestedDict(params) def getLookupDir(): """ Return the location of the directory holding lookups files. :return: absolute path to the directory """ this_dir, _ = os.path.split(__file__) return os.path.join(this_dir, 'lookups') -def getNmodlDir(): - """ Return the location of the directory holding the NEURON mechanisms dll. - - :return: absolute path to the directory - """ - this_dir, _ = os.path.split(__file__) - return os.path.join(this_dir, 'neurons\\nmodl') - - def ParseNestedDict(dict_in): """ Loop through a nested dictionary object and convert all string fields to floats. """ for key, value in dict_in.items(): if isinstance(value, dict): # If value itself is dictionary dict_in[key] = ParseNestedDict(value) elif isinstance(dict_in[key], str): dict_in[key] = float(dict_in[key]) return dict_in def OpenFilesDialog(filetype, dirname=''): """ Open a FileOpenDialogBox to select one or multiple file. The default directory and file type are given. :param dirname: default directory :param filetype: default file type :return: tuple of full paths to the chosen filenames """ root = tk.Tk() root.withdraw() filenames = filedialog.askopenfilenames(filetypes=[(filetype + " files", '.' + filetype)], initialdir=dirname) if filenames: par_dir = os.path.abspath(os.path.join(filenames[0], os.pardir)) else: par_dir = None return (filenames, par_dir) def ImportExcelCol(filename, sheetname, colstr, startrow): """ Load a specific column of an excel workbook as a numpy array. :param filename: absolute or relative path to the Excel workbook :param sheetname: name of the Excel spreadsheet to which data is appended :param colstr: string of the column to import :param startrow: index of the first row to consider :return: 1D numpy array with the column data """ wb = load_workbook(filename, read_only=True) ws = wb[sheetname] range_start_str = colstr + str(startrow) range_stop_str = colstr + str(ws.max_row) tmp = np.array([[i.value for i in j] for j in ws[range_start_str:range_stop_str]]) return tmp[:, 0] def ConstructMatrix(serialized_inputA, serialized_inputB, serialized_output): """ Construct a 2D output matrix from serialized input. :param serialized_inputA: serialized input variable A :param serialized_inputB: serialized input variable B :param serialized_output: serialized output variable :return: 4-tuple with vectors of unique values of A (m) and B (n), output variable 2D matrix (m,n) and number of holes in the matrix """ As = np.unique(serialized_inputA) Bs = np.unique(serialized_inputB) nA = As.size nB = Bs.size output = np.zeros((nA, nB)) output[:] = np.NAN nholes = 0 for i in range(nA): iA = np.where(serialized_inputA == As[i]) for j in range(nB): iB = np.where(serialized_inputB == Bs[j]) iMatch = np.intersect1d(iA, iB) if iMatch.size == 0: nholes += 1 elif iMatch.size > 1: logger.warning('Identical serialized inputs with values (%f, %f)', As[i], Bs[j]) else: iMatch = iMatch[0] output[i, j] = serialized_output[iMatch] return (As, Bs, output, nholes) def rmse(x1, x2): """ Compute the root mean square error between two 1D arrays """ return np.sqrt(((x1 - x2) ** 2).mean()) def rsquared(x1, x2): ''' compute the R-squared coefficient between two 1D arrays ''' residuals = x1 - x2 ss_res = np.sum(residuals**2) ss_tot = np.sum((x1 - np.mean(x1))**2) return 1 - (ss_res / ss_tot) def DownSample(t_dense, y, nsparse): """ Decimate periodic signals to a specified number of samples.""" if(y.ndim) > 1: nsignals = y.shape[0] else: nsignals = 1 y = np.array([y]) # determine time step and period of input signal T = t_dense[-1] - t_dense[0] dt_dense = t_dense[1] - t_dense[0] # resample time vector linearly t_ds = np.linspace(t_dense[0], t_dense[-1], nsparse) # create MAV window nmav = int(0.03 * T / dt_dense) if nmav % 2 == 0: nmav += 1 mav = np.ones(nmav) / nmav # determine signals padding npad = int((nmav - 1) / 2) # determine indexes of sampling on convolved signals ids = np.round(np.linspace(0, t_dense.size - 1, nsparse)).astype(int) y_ds = np.empty((nsignals, nsparse)) # loop through signals for i in range(nsignals): # pad, convolve and resample pad_left = y[i, -(npad + 2):-2] pad_right = y[i, 1:npad + 1] y_ext = np.concatenate((pad_left, y[i, :], pad_right), axis=0) y_mav = np.convolve(y_ext, mav, mode='valid') y_ds[i, :] = y_mav[ids] if nsignals == 1: y_ds = y_ds[0, :] return (t_ds, y_ds) def Pressure2Intensity(p, rho=1075.0, c=1515.0): """ Return the spatial peak, pulse average acoustic intensity (ISPPA) associated with the specified pressure amplitude. :param p: pressure amplitude (Pa) :param rho: medium density (kg/m3) :param c: speed of sound in medium (m/s) :return: spatial peak, pulse average acoustic intensity (W/m2) """ return p**2 / (2 * rho * c) def Intensity2Pressure(I, rho=1075.0, c=1515.0): """ Return the pressure amplitude associated with the specified spatial peak, pulse average acoustic intensity (ISPPA). :param I: spatial peak, pulse average acoustic intensity (W/m2) :param rho: medium density (kg/m3) :param c: speed of sound in medium (m/s) :return: pressure amplitude (Pa) """ return np.sqrt(2 * rho * c * I) def find_nearest(array, value): ''' Find nearest element in 1D array. ''' idx = (np.abs(array - value)).argmin() return (idx, array[idx]) def rescale(x, lb=None, ub=None, lb_new=0, ub_new=1): ''' Rescale a value to a specific interval by linear transformation. ''' if lb is None: lb = x.min() if ub is None: ub = x.max() xnorm = (x - lb) / (ub - lb) return xnorm * (ub_new - lb_new) + lb_new def LennardJones(x, beta, alpha, C, m, n): """ Generic expression of a Lennard-Jones function, adapted for the context of symmetric deflection (distance = 2x). :param x: deflection (i.e. half-distance) :param beta: x-shifting factor :param alpha: x-scaling factor :param C: y-scaling factor :param m: exponent of the repulsion term :param n: exponent of the attraction term :return: Lennard-Jones potential at given distance (2x) """ return C * (np.power((alpha / (2 * x + beta)), m) - np.power((alpha / (2 * x + beta)), n)) def getNeuronsDict(): ''' Return dictionary of neurons classes that can be instantiated. ''' neurons_dict = {} for _, obj in inspect.getmembers(neurons): if inspect.isclass(obj) and isinstance(obj.name, str): neurons_dict[obj.name] = obj return neurons_dict def get_BLS_lookups(a): lookup_path = getLookupDir() + '/BLS_lookups_a{:.1f}nm.json'.format(a * 1e9) try: with open(lookup_path) as fh: sample = json.load(fh) return sample except FileNotFoundError: return {} def save_BLS_lookups(a, lookups): """ Save BLS parameter into specific lookup file :return: absolute path to the directory """ lookup_path = getLookupDir() + '/BLS_lookups_a{:.1f}nm.json'.format(a * 1e9) with open(lookup_path, 'w') as fh: json.dump(lookups, fh) def extractCompTimes(filenames): ''' Extract computation times from a list of simulation files. ''' tcomps = np.empty(len(filenames)) for i, fn in enumerate(filenames): logger.info('Loading data from "%s"', fn) with open(fn, 'rb') as fh: frame = pickle.load(fh) meta = frame['meta'] tcomps[i] = meta['tcomp'] return tcomps def computeMeshEdges(x, scale='lin'): ''' Compute the appropriate edges of a mesh that quads a linear or logarihtmic distribution. :param x: the input vector :param scale: the type of distribution ('lin' for linear, 'log' for logarihtmic) :return: the edges vector ''' if scale is 'log': x = np.log10(x) dx = x[1] - x[0] if scale is 'lin': y = np.linspace(x[0] - dx / 2, x[-1] + dx / 2, x.size + 1) elif scale is 'log': y = np.logspace(x[0] - dx / 2, x[-1] + dx / 2, x.size + 1) return y si_prefixes = { 'y': 1e-24, # yocto 'z': 1e-21, # zepto 'a': 1e-18, # atto 'f': 1e-15, # femto 'p': 1e-12, # pico 'n': 1e-9, # nano 'u': 1e-6, # micro 'm': 1e-3, # mili '': 1e0, # None 'k': 1e3, # kilo 'M': 1e6, # mega 'G': 1e9, # giga 'T': 1e12, # tera 'P': 1e15, # peta 'E': 1e18, # exa 'Z': 1e21, # zetta 'Y': 1e24, # yotta } def si_format(x, precision=0, space=''): ''' Format a float according to the SI unit system, with the appropriate prefix letter. ''' if isinstance(x, float) or isinstance(x, int) or isinstance(x, np.float) or\ isinstance(x, np.int32) or isinstance(x, np.int64): if x == 0: factor = 1e0 prefix = '' else: sorted_si_prefixes = sorted(si_prefixes.items(), key=operator.itemgetter(1)) vals = [tmp[1] for tmp in sorted_si_prefixes] # vals = list(si_prefixes.values()) ix = np.searchsorted(vals, np.abs(x)) - 1 if np.abs(x) == vals[ix + 1]: ix += 1 factor = vals[ix] prefix = sorted_si_prefixes[ix][0] # prefix = list(si_prefixes.keys())[ix] return '{{:.{}f}}{}{}'.format(precision, space, prefix).format(x / factor) elif isinstance(x, list) or isinstance(x, tuple): return [si_format(item, precision, space) for item in x] elif isinstance(x, np.ndarray) and x.ndim == 1: return [si_format(float(item), precision, space) for item in x] else: print(type(x)) def getCycleAverage(t, y, T): ''' Compute the cycle-averaged profile of a signal given a specific periodicity. :param t: time vector (s) :param y: signal vector :param T: period (s) :return: cycle-averaged signal vector ''' nsamples = y.size ncycles = int((t[-1] - t[0]) // T) npercycle = int(nsamples // ncycles) return np.mean(np.reshape(y[:ncycles * npercycle], (ncycles, npercycle)), axis=1) def itrpLookupsFreq(lookups3D, freqs, Fdrive): """ Interpolate a dictionary of 3D lookups at a given frequency. :param lookups3D: dictionary of 3D lookups :param freqs: array of lookup frequencies (Hz) :param Fdrive: acoustic drive frequency (Hz) :return: a dictionary of 2D lookups interpolated a the given frequency """ # If Fdrive in lookup frequencies, simply take (A, Q) slice at Fdrive index if Fdrive in freqs: iFdrive = np.searchsorted(freqs, Fdrive) - logger.debug('Using lookups directly at %.2f kHz', freqs[iFdrive] * 1e-3) + # logger.debug('Using lookups directly at %.2f kHz', freqs[iFdrive] * 1e-3) lookups2D = {key: np.squeeze(lookups3D[key][iFdrive, :, :]) for key in lookups3D.keys()} # Otherwise, interpolate linearly between 2 (A, Q) slices at Fdrive bounding values indexes else: ilb = np.searchsorted(freqs, Fdrive) - 1 iub = ilb + 1 - logger.debug('Interpolating lookups between %.2f kHz and %.2f kHz', - freqs[ilb] * 1e-3, freqs[iub] * 1e-3) + # logger.debug('Interpolating lookups between %.2f kHz and %.2f kHz', + # freqs[ilb] * 1e-3, freqs[iub] * 1e-3) lookups2D_lb = {key: np.squeeze(lookups3D[key][ilb, :, :]) for key in lookups3D.keys()} lookups2D_ub = {key: np.squeeze(lookups3D[key][iub, :, :]) for key in lookups3D.keys()} Fratio = (Fdrive - freqs[ilb]) / (freqs[iub] - freqs[ilb]) lookups2D = {key: lookups2D_lb[key] + (lookups2D_ub[key] - lookups2D_lb[key]) * Fratio for key in lookups3D.keys()} return lookups2D +def getLookups2D(mechname, a, Fdrive): + ''' Retrieve appropriate 2D lookup tables and reference vectors + for a given membrane mechanism, sonophore diameter and US frequency. + + :param mechname: name of membrane density mechanism + :param a: sonophore diameter (m) + :param Fdrive: US frequency (Hz) + :return: 3-tuple with 1D numpy arrays of reference acoustic amplitudes and charge densities, + and a dictionary of 2D lookup numpy arrays + ''' + + # Check lookup file existence + lookup_file = '{}_lookups_a{:.1f}nm.pkl'.format(mechname, a * 1e9) + lookup_path = '{}/{}'.format(getLookupDir(), lookup_file) + if not os.path.isfile(lookup_path): + raise InputError('Missing lookup file: "{}"'.format(lookup_file)) + + # Load lookups dictionary + with open(lookup_path, 'rb') as fh: + lookups3D = pickle.load(fh) + + # Retrieve 1D inputs from lookups dictionary + Fref = lookups3D.pop('f') + Aref = lookups3D.pop('A') + Qref = lookups3D.pop('Q') + + # Check that US frequency is within lookup range + margin = 1e-9 # adding margin to compensate for eventual round error + Frange = (Fref.min() - margin, Fref.max() + margin) + if Fdrive < Frange[0] or Fdrive > Frange[1]: + raise InputError('Invalid frequency: {}Hz (must be within {}Hz - {}Hz lookup interval)' + .format(*si_format([Fdrive, *Frange], precision=2, space=' '))) + + # Interpolate 3D lookups at US frequency + lookups2D = itrpLookupsFreq(lookups3D, Fref, Fdrive) + + return Aref, Qref, lookups2D + + def nDindexes(dims, index): ''' Find index positions in a n-dimensional array. :param dims: dimensions of the n-dimensional array (tuple or list) :param index: index position in the flattened n-dimensional array :return: list of indexes along each array dimension ''' dims = list(dims) # Find size of each array dimension dims.reverse() dimsizes = [1] r = 1 for x in dims[:-1]: r *= x dimsizes.append(r) dims.reverse() dimsizes.reverse() # Find indexes indexes = [] remainder = index for dimsize in dimsizes[:-1]: idim, remainder = divmod(remainder, dimsize) indexes.append(idim) indexes.append(remainder) - return indexes \ No newline at end of file + return indexes + + +def pow10_format(number, precision=2): + ''' Format a number in power of 10 notation. ''' + ret_string = '{0:.{1:d}e}'.format(number, precision) + a, b = ret_string.split("e") + a = float(a) + b = int(b) + return '{}10^{{{}}}'.format('{} * '.format(a) if a != 1. else '', b) diff --git a/plot/plot_eff_coeffs.py b/plot/plot_eff_coeffs.py index b3f15e5..6ef2b00 100644 --- a/plot/plot_eff_coeffs.py +++ b/plot/plot_eff_coeffs.py @@ -1,34 +1,34 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-02-15 15:59:37 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-05-25 22:04:17 +# @Last Modified time: 2018-08-06 11:01:27 import numpy as np import matplotlib.pyplot as plt from PointNICE.plt import plotEffVars from PointNICE.neurons import * ''' Plot the profiles of effective variables as a function of charge density with amplitude color code. ''' # Set parameters neuron = CorticalRS() Fdrive = 500e3 amps = np.logspace(np.log10(1), np.log10(600), 10) * 1e3 charges = np.linspace(neuron.Vm0, 50, 100) * 1e-5 # Define variables to plot gates = ['m', 'h', 'n', 'p'] keys = ['V', 'ng'] for x in gates: keys += ['alpha{}'.format(x), 'beta{}'.format(x)] # Plot effective variables fig = plotEffVars(neuron, Fdrive, amps=amps, charges=charges, keys=keys, ncolmax=2, fs=8) plt.show() diff --git a/tests/comp_traces_Python_NEURON.py b/tests/comp_traces_Python_NEURON.py deleted file mode 100644 index 9568598..0000000 --- a/tests/comp_traces_Python_NEURON.py +++ /dev/null @@ -1,87 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -# @Author: Theo Lemaire -# @Date: 2017-02-13 18:16:09 -# @Email: theo.lemaire@epfl.ch -# @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-07-24 14:43:59 - -""" Compare output of pure-Python and NEURON acoustic simulations. """ - -import logging -import pickle -import matplotlib.pyplot as plt - -from PointNICE.utils import logger, getNeuronsDict, si_format -from PointNICE.solvers import SolverUS, checkBatchLog, AStimWorker, setBatchDir -from PointNICE.plt import plotComp - -# Set logging level -logger.setLevel(logging.INFO) - -# Set neurons -neurons = ['RS', 'FS', 'LTS', 'RE', 'TC'] - -# Set stimulation parameters -a = 32e-9 # m -Adrive = 100e3 # Pa -Fdrive = 500e3 # Hz -PRF = 100. # Hz -DC = .5 -tstim = 150e-3 # s -toffset = 100e-3 # s - -ASTIM_CW_code = 'ASTIM_{}_CW_{:.0f}nm_{:.0f}kHz_{:.1f}kPa_{:.0f}ms_{}' -ASTIM_PW_code = 'ASTIM_{}_PW_{:.0f}nm_{:.0f}kHz_{:.1f}kPa_{:.0f}ms_PRF{:.2f}Hz_DC{:.2f}%_{}' -suffix = 'Python_vs_NEURON' - -# Select output directory -batch_dir = setBatchDir() -mu_acc_factor = 0 - -# For each neuron -for nname in neurons: - - if DC == 1: - simcode = ASTIM_CW_code.format(nname, a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3, - suffix) - else: - simcode = ASTIM_PW_code.format(nname, a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3, - PRF, DC * 1e2, suffix) - - neuron = getNeuronsDict()[nname]() - - # Initialize solver - solver = SolverUS(a, neuron, Fdrive) - - log_filepath, _ = checkBatchLog(batch_dir, 'A-STIM') - - # Run NEURON and Python A-STIM point-neuron simulations - logger.info('Running simulations for %s neuron', nname) - NEURONsim_file = AStimWorker(1, batch_dir, log_filepath, solver, neuron, Fdrive, Adrive, - tstim, toffset, PRF, DC, 'NEURON', 2).__call__() - effsim_file = AStimWorker(1, batch_dir, log_filepath, solver, neuron, Fdrive, Adrive, - tstim, toffset, PRF, DC, 'effective', 2).__call__() - - # Compare computation times - with open(effsim_file, 'rb') as fh: - tcomp_eff = pickle.load(fh)['meta']['tcomp'] - logger.info('Python simulation completed in %.1f s', tcomp_eff) - with open(NEURONsim_file, 'rb') as fh: - tcomp_NEURON = pickle.load(fh)['meta']['tcomp'] - logger.info('NEURON simulation completed in %.1f ms (%.1f times faster)', - tcomp_NEURON * 1e3, tcomp_eff / tcomp_NEURON) - mu_acc_factor += tcomp_eff / tcomp_NEURON - - # Plot resulting profiles - filepaths = [effsim_file, NEURONsim_file] - fig = plotComp('Qm', filepaths, labels=['Python', 'NEURON'], showfig=False) - ax = fig.gca() - ax.set_title('{} neuron - {}Hz, {}Pa'.format(nname, *si_format([Fdrive, Adrive], space=' '))) - fig.tight_layout() - # fig.savefig('{}/{}.png'.format(batch_dir, simcode)) - -mu_acc_factor /= len(neurons) -logger.info('average acceleration factor = %.1f', mu_acc_factor) - -plt.show() diff --git a/tests/test_basic.py b/tests/test_basic.py index 769b15d..fcaa889 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -1,239 +1,260 @@ #!/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: 2018-03-19 19:59:39 +# @Last Modified time: 2018-08-21 14:24:51 ''' Test the basic functionalities of the package. ''' import os import sys import logging import time import cProfile import pstats from argparse import ArgumentParser -from PointNICE.utils import logger +from PointNICE.utils import logger, InputError from PointNICE import BilayerSonophore, SolverElec, SolverUS from PointNICE.neurons import * def test_MECH(is_profiled=False): ''' Mechanical simulation. ''' logger.info('Test: running MECH 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) # Stimulation parameters Adrive = 100e3 # Pa Qm = 50e-5 # C/m2 # Run simulation if is_profiled: pfile = 'tmp.stats' cProfile.runctx('bls.run(Fdrive, Adrive, Qm)', globals(), locals(), pfile) stats = pstats.Stats(pfile) os.remove(pfile) stats.strip_dirs() stats.sort_stats('cumulative') stats.print_stats() else: bls.run(Fdrive, Adrive, Qm) def test_ESTIM(is_profiled=False): ''' Electrical simulation ''' logger.info('Test: running ESTIM simulation') # Initialize neuron neuron = CorticalRS() # Initialize solver solver = SolverElec() # Stimulation parameters Astim = 10.0 # mA/m2 tstim = 100e-3 # s toffset = 50e-3 # s # Run simulation if is_profiled: pfile = 'tmp.stats' cProfile.runctx('solver.run(neuron, Astim, tstim, toffset)', globals(), locals(), pfile) stats = pstats.Stats(pfile) os.remove(pfile) stats.strip_dirs() stats.sort_stats('cumulative') stats.print_stats() else: solver.run(neuron, Astim, tstim, toffset, PRF=None, DC=1.0) solver.run(neuron, Astim, tstim, toffset, PRF=100.0, DC=0.05) def test_ASTIM_effective(is_profiled=False): ''' Effective acoustic simulation ''' - logger.info('Test: running ASTIM effective simulation') + logger.info('Test: ASTIM effective simulation') - # Initialize neuron + # Default parameters neuron = CorticalRS() - - # Stimulation parameters + a = 32e-9 # m Fdrive = 500e3 # Hz Adrive = 100e3 # Pa tstim = 50e-3 # s toffset = 10e-3 # s - # Initialize solver - a = 32e-9 # m - solver = SolverUS(a, neuron, Fdrive) # Run simulation if is_profiled: + solver = SolverUS(a, neuron, Fdrive) pfile = 'tmp.stats' cProfile.runctx("solver.run(neuron, Fdrive, Adrive, tstim, toffset, sim_type='effective')", globals(), locals(), pfile) stats = pstats.Stats(pfile) os.remove(pfile) stats.strip_dirs() stats.sort_stats('cumulative') stats.print_stats() else: + # test error 1: no lookups exist for sonophore diameter + try: + solver = SolverUS(50e-9, neuron, Fdrive) + solver.run(neuron, Fdrive, Adrive, tstim, toffset, sim_type='effective') + except InputError as err: + logger.debug('No lookups for sonophore diameter: OK') + + # test error 2: frequency outside of lookups range + Foutside = 10e3 # Hz + try: + solver = SolverUS(a, neuron, Foutside) + solver.run(neuron, Foutside, Adrive, tstim, toffset, sim_type='effective') + except InputError as err: + logger.debug('Out of range frequency: OK') + + # test error 3: amplitude outside of lookups range + Aoutside = 1e6 # Pa + try: + solver = SolverUS(a, neuron, Fdrive) + solver.run(neuron, Fdrive, Aoutside, tstim, toffset, sim_type='effective') + except InputError as err: + logger.debug('Out of range amplitude: OK') + + # test: normal stimulation completion (CW and PW) + solver = SolverUS(a, neuron, Fdrive) solver.run(neuron, Fdrive, Adrive, tstim, toffset, sim_type='effective') solver.run(neuron, Fdrive, Adrive, tstim, toffset, PRF=100.0, DC=0.05, sim_type='effective') - def test_ASTIM_classic(is_profiled=False): ''' Classic acoustic simulation ''' logger.info('Test: running ASTIM classic simulation') # Initialize neuron neuron = CorticalRS() # Stimulation parameters Fdrive = 500e3 # Hz Adrive = 100e3 # Pa tstim = 1e-6 # s toffset = 1e-6 # s # Initialize solver a = 32e-9 # m solver = SolverUS(a, neuron, Fdrive) # Run simulation if is_profiled: pfile = 'tmp.stats' cProfile.runctx("solver.run(neuron, Fdrive, Adrive, tstim, toffset, sim_type='classic')", globals(), locals(), pfile) stats = pstats.Stats(pfile) os.remove(pfile) stats.strip_dirs() stats.sort_stats('cumulative') stats.print_stats() else: solver.run(neuron, Fdrive, Adrive, tstim, toffset, sim_type='classic') def test_ASTIM_hybrid(is_profiled=False): ''' Hybrid acoustic simulation ''' logger.info('Test: running ASTIM hybrid simulation') # Initialize neuron neuron = CorticalRS() # Stimulation parameters Fdrive = 350e3 # Hz Adrive = 100e3 # Pa tstim = 1e-3 # s toffset = 1e-3 # s # Initialize solver a = 32e-9 # m solver = SolverUS(a, neuron, Fdrive) # Run simulation if is_profiled: pfile = 'tmp.stats' cProfile.runctx("solver.run(neuron, Fdrive, Adrive, tstim, toffset, sim_type='hybrid')", globals(), locals(), pfile) stats = pstats.Stats(pfile) os.remove(pfile) stats.strip_dirs() stats.sort_stats('cumulative') stats.print_stats() else: solver.run(neuron, Fdrive, Adrive, tstim, toffset, sim_type='hybrid') def test_all(): t0 = time.time() test_MECH() test_ESTIM() test_ASTIM_effective() test_ASTIM_classic() test_ASTIM_hybrid() tcomp = time.time() - t0 logger.info('All tests completed in %.0f s', tcomp) def main(): # Define valid test sets valid_testsets = [ 'MECH', 'ESTIM', 'ASTIM_effective', 'ASTIM_classic', 'ASTIM_hybrid', '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') ap.add_argument('-p', '--profile', default=False, action='store_true', help='Profile test set') # Parse arguments args = ap.parse_args() if args.verbose: logger.setLevel(logging.DEBUG) else: logger.setLevel(logging.INFO) if args.profile and args.testset == 'all': logger.error('profiling can only be run on individual tests') sys.exit(2) # Run test if args.testset == 'all': test_all() else: possibles = globals().copy() possibles.update(locals()) method = possibles.get('test_{}'.format(args.testset)) method(args.profile) sys.exit(0) if __name__ == '__main__': main()