diff --git a/PySONIC/plt/pltutils.py b/PySONIC/plt/pltutils.py index c808a1f..ddb314c 100644 --- a/PySONIC/plt/pltutils.py +++ b/PySONIC/plt/pltutils.py @@ -1,1288 +1,1309 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-23 14:55:37 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-08-28 15:36:12 +# @Last Modified time: 2018-08-29 18:02:07 ''' Plotting utilities ''' import sys import os import pickle import ntpath import re import logging import tkinter as tk from tkinter import filedialog import numpy as np from scipy.interpolate import interp2d # from scipy.optimize import brentq import matplotlib import matplotlib.pyplot as plt from matplotlib.patches import Rectangle import matplotlib.cm as cm from matplotlib.ticker import FormatStrFormatter import pandas as pd from .. import neurons from ..constants import DT_EFF from ..utils import getNeuronsDict, getLookupDir, rescale, InputError, computeMeshEdges, si_format, itrpLookupsFreq from ..bls import BilayerSonophore from .pltvars import pltvars matplotlib.rcParams['pdf.fonttype'] = 42 matplotlib.rcParams['ps.fonttype'] = 42 matplotlib.rcParams['font.family'] = 'arial' # Get package logger logger = logging.getLogger('PySONIC') # Define global variables neuron = None bls = None timeunits = {'ASTIM': 't_ms', 'ESTIM': 't_ms', 'MECH': 't_us'} # Regular expression for input files rgxp = re.compile('(ESTIM|ASTIM)_([A-Za-z]*)_(.*).pkl') rgxp_mech = re.compile('(MECH)_(.*).pkl') # nb = '[0-9]*[.]?[0-9]+' # rgxp_ASTIM = re.compile('(ASTIM)_(\w+)_(PW|CW)_({0})nm_({0})kHz_({0})kPa_({0})ms(.*)_(\w+).pkl'.format(nb)) # rgxp_ESTIM = re.compile('(ESTIM)_(\w+)_(PW|CW)_({0})mA_per_m2_({0})ms(.*).pkl'.format(nb)) # rgxp_PW = re.compile('_PRF({0})kHz_DC({0})_(PW|CW)_(\d+)kHz_(\d+)kPa_(\d+)ms_(.*).pkl'.format(nb)) # Figure naming conventions ESTIM_CW_title = '{} neuron: CW E-STIM {:.2f}mA/m2, {:.0f}ms' ESTIM_PW_title = '{} neuron: PW E-STIM {:.2f}mA/m2, {:.0f}ms, {:.2f}Hz PRF, {:.0f}% DC' ASTIM_CW_title = '{} neuron: CW A-STIM {:.0f}kHz, {:.0f}kPa, {:.0f}ms' ASTIM_PW_title = '{} neuron: PW A-STIM {:.0f}kHz, {:.0f}kPa, {:.0f}ms, {:.2f}Hz PRF, {:.2f}% DC' MECH_title = '{:.0f}nm BLS structure: MECH-STIM {:.0f}kHz, {:.0f}kPa' def cm2inch(*tupl): inch = 2.54 if isinstance(tupl[0], tuple): return tuple(i / inch for i in tupl[0]) else: return tuple(i / inch for i in tupl) class InteractiveLegend(object): """ Class defining an interactive matplotlib legend, where lines visibility can be toggled by simply clicking on the corresponding legend label. Other graphic objects can also be associated to the toggle of a specific line Adapted from: http://stackoverflow.com/questions/31410043/hiding-lines-after-showing-a-pyplot-figure """ def __init__(self, legend, aliases): self.legend = legend self.fig = legend.axes.figure self.lookup_artist, self.lookup_handle = self._build_lookups(legend) self._setup_connections() self.handles_aliases = aliases self.update() def _setup_connections(self): for artist in self.legend.texts + self.legend.legendHandles: artist.set_picker(10) # 10 points tolerance self.fig.canvas.mpl_connect('pick_event', self.on_pick) def _build_lookups(self, legend): ''' Method of the InteractiveLegend class building the legend lookups. ''' labels = [t.get_text() for t in legend.texts] handles = legend.legendHandles label2handle = dict(zip(labels, handles)) handle2text = dict(zip(handles, legend.texts)) lookup_artist = {} lookup_handle = {} for artist in legend.axes.get_children(): if artist.get_label() in labels: handle = label2handle[artist.get_label()] lookup_handle[artist] = handle lookup_artist[handle] = artist lookup_artist[handle2text[handle]] = artist lookup_handle.update(zip(handles, handles)) lookup_handle.update(zip(legend.texts, handles)) return lookup_artist, lookup_handle def on_pick(self, event): handle = event.artist if handle in self.lookup_artist: artist = self.lookup_artist[handle] artist.set_visible(not artist.get_visible()) self.update() def update(self): for artist in self.lookup_artist.values(): handle = self.lookup_handle[artist] if artist.get_visible(): handle.set_visible(True) if artist in self.handles_aliases: for al in self.handles_aliases[artist]: al.set_visible(True) else: handle.set_visible(False) if artist in self.handles_aliases: for al in self.handles_aliases[artist]: al.set_visible(False) self.fig.canvas.draw() def show(self): ''' showing the interactive legend ''' plt.show() def getPatchesLoc(t, states): ''' Determine the location of stimulus patches. :param t: simulation time vector (s). :param states: a vector of stimulation state (ON/OFF) at each instant in time. :return: 3-tuple with number of patches, timing of STIM-ON an STIM-OFF instants. ''' # Compute states derivatives and identify bounds indexes of pulses dstates = np.diff(states) ipatch_on = np.insert(np.where(dstates > 0.0)[0] + 1, 0, 0) ipatch_off = np.where(dstates < 0.0)[0] + 1 if ipatch_off.size < ipatch_on.size: ioff = t.size - 1 if ipatch_off.size == 0: ipatch_off = np.array([ioff]) else: ipatch_off = np.insert(ipatch_off, ipatch_off.size - 1, ioff) # Get time instants for pulses ON and OFF npatches = ipatch_on.size tpatch_on = t[ipatch_on] tpatch_off = t[ipatch_off] # return 3-tuple with #patches, pulse ON and pulse OFF instants return (npatches, tpatch_on, tpatch_off) def SaveFigDialog(dirname, filename): """ Open a FileSaveDialogBox to set the directory and name of the figure to be saved. The default directory and filename are given, and the default extension is ".pdf" :param dirname: default directory :param filename: default filename :return: full path to the chosen filename """ root = tk.Tk() root.withdraw() filename_out = filedialog.asksaveasfilename(defaultextension=".pdf", initialdir=dirname, initialfile=filename) return filename_out def plotComp(varname, filepaths, labels=None, fs=15, lw=2, colors=None, lines=None, patches='one', xticks=None, yticks=None, blacklegend=False, straightlegend=False, showfig=True, inset=None, figsize=(11, 4)): ''' Compare profiles of several specific output variables of NICE simulations. :param varname: name of variable to extract and compare :param filepaths: list of full paths to output data files to be compared :param labels: list of labels to use in the legend :param fs: labels fontsize :param patches: string indicating whether to indicate periods of stimulation with colored rectangular patches ''' # Input check 1: variable name if varname not in pltvars: raise InputError('Unknown plot variable: "{}"'.format(varname)) pltvar = pltvars[varname] # Input check 2: labels if labels is not None: if len(labels) != len(filepaths): raise InputError('Invalid labels ({}): not matching number of compared files ({})' .format(len(labels), len(filepaths))) if not all(isinstance(x, str) for x in labels): raise InputError('Invalid labels: must be string typed') # Input check 3: line styles and colors if colors is None: colors = ['C{}'.format(j) for j in range(len(filepaths))] if lines is None: lines = ['-'] * len(filepaths) # Input check 4: STIM-ON patches greypatch = False if patches == 'none': patches = [False] * len(filepaths) elif patches == 'all': patches = [True] * len(filepaths) elif patches == 'one': patches = [True] + [False] * (len(filepaths) - 1) greypatch = True elif isinstance(patches, list): if len(patches) != len(filepaths): raise InputError('Invalid patches ({}): not matching number of compared files ({})' .format(len(patches), len(filepaths))) if not all(isinstance(p, bool) for p in patches): raise InputError('Invalid patch sequence: all list items must be boolean typed') else: raise InputError('Invalid patches: must be either "none", all", "one", or a boolean list') # Initialize figure and axis fig, ax = plt.subplots(figsize=figsize) ax.set_zorder(0) for item in ['top', 'right']: ax.spines[item].set_visible(False) if 'min' in pltvar and 'max' in pltvar: # optional min and max on y-axis ax.set_ylim(pltvar['min'], pltvar['max']) if pltvar['unit']: # y-label with optional unit ax.set_ylabel('$\\rm {}\ ({})$'.format(pltvar['label'], pltvar['unit']), fontsize=fs) else: ax.set_ylabel('$\\rm {}$'.format(pltvar['label']), fontsize=fs) if xticks is not None: # optional x-ticks ax.set_xticks(xticks) if yticks is not None: # optional y-ticks ax.set_yticks(yticks) else: ax.locator_params(axis='y', nbins=2) if any(ax.get_yticks() < 0): ax.yaxis.set_major_formatter(FormatStrFormatter('%+.0f')) for tick in ax.xaxis.get_major_ticks() + ax.yaxis.get_major_ticks(): tick.label.set_fontsize(fs) # Optional inset axis if inset is not None: inset_ax = fig.add_axes(ax.get_position()) inset_ax.set_xlim(inset['xlims'][0], inset['xlims'][1]) inset_ax.set_ylim(inset['ylims'][0], inset['ylims'][1]) inset_ax.set_xticks([]) inset_ax.set_yticks([]) # inset_ax.patch.set_alpha(1.0) inset_ax.set_zorder(1) inset_ax.add_patch(Rectangle((inset['xlims'][0], inset['ylims'][0]), inset['xlims'][1] - inset['xlims'][0], inset['ylims'][1] - inset['ylims'][0], color='w')) # Retrieve neurons dictionary neurons_dict = getNeuronsDict() # Loop through data files aliases = {} for j, filepath in enumerate(filepaths): # Retrieve sim type pkl_filename = ntpath.basename(filepath) mo1 = rgxp.fullmatch(pkl_filename) mo2 = rgxp_mech.fullmatch(pkl_filename) if mo1: mo = mo1 elif mo2: mo = mo2 else: logger.error('Error: "%s" file does not match regexp pattern', pkl_filename) sys.exit(1) sim_type = mo.group(1) if sim_type not in ('MECH', 'ASTIM', 'ESTIM'): raise InputError('Invalid simulation type: {}'.format(sim_type)) if j == 0: sim_type_ref = sim_type t_plt = pltvars[timeunits[sim_type]] elif sim_type != sim_type_ref: raise InputError('Invalid comparison: different simulation types') # Load data logger.info('Loading data from "%s"', pkl_filename) with open(filepath, 'rb') as fh: frame = pickle.load(fh) df = frame['data'] meta = frame['meta'] # Extract variables t = df['t'].values states = df['states'].values nsamples = t.size # Initialize neuron object if ESTIM or ASTIM sim type if sim_type in ['ASTIM', 'ESTIM']: neuron_name = mo.group(2) global neuron neuron = neurons_dict[neuron_name]() Cm0 = neuron.Cm0 Qm0 = Cm0 * neuron.Vm0 * 1e-3 # Extract neuron states if needed if 'alias' in pltvar and 'neuron_states' in pltvar['alias']: neuron_states = [df[sn].values for sn in neuron.states_names] else: Cm0 = meta['Cm0'] Qm0 = meta['Qm0'] # Initialize BLS if needed if sim_type in ['MECH', 'ASTIM'] and 'alias' in pltvar and 'bls' in pltvar['alias']: global bls bls = BilayerSonophore(meta['a'], meta['Fdrive'], Cm0, Qm0) # Determine patches location npatches, tpatch_on, tpatch_off = getPatchesLoc(t, states) # Add onset to time vectors if t_plt['onset'] > 0.0: tonset = np.array([-t_plt['onset'], -t[0] - t[1]]) t = np.hstack((tonset, t)) states = np.hstack((states, np.zeros(2))) # Set x-axis label ax.set_xlabel('$\\rm {}\ ({})$'.format(t_plt['label'], t_plt['unit']), fontsize=fs) # Extract variable to plot if 'alias' in pltvar: var = eval(pltvar['alias']) elif 'key' in pltvar: var = df[pltvar['key']].values elif 'constant' in pltvar: var = eval(pltvar['constant']) * np.ones(nsamples) else: var = df[varname].values if var.size == t.size - 2: if varname is 'Vm': var = np.hstack((np.array([neuron.Vm0] * 2), var)) else: var = np.hstack((np.array([var[0]] * 2), var)) # var = np.insert(var, 0, var[0]) # Determine legend label if labels is not None: label = labels[j] else: if sim_type == 'ESTIM': if meta['DC'] == 1.0: label = ESTIM_CW_title.format(neuron_name, meta['Astim'], meta['tstim'] * 1e3) else: label = ESTIM_PW_title.format(neuron_name, meta['Astim'], meta['tstim'] * 1e3, meta['PRF'], meta['DC'] * 1e2) elif sim_type == 'ASTIM': if meta['DC'] == 1.0: label = ASTIM_CW_title.format(neuron_name, meta['Fdrive'] * 1e-3, meta['Adrive'] * 1e-3, meta['tstim'] * 1e3) else: label = ASTIM_PW_title.format(neuron_name, meta['Fdrive'] * 1e-3, meta['Adrive'] * 1e-3, meta['tstim'] * 1e3, meta['PRF'], meta['DC'] * 1e2) elif sim_type == 'MECH': label = MECH_title.format(meta['a'] * 1e9, meta['Fdrive'] * 1e-3, meta['Adrive'] * 1e-3) # Plot trace handle = ax.plot(t * t_plt['factor'], var * pltvar['factor'], linewidth=lw, linestyle=lines[j], color=colors[j], label=label) if inset is not None: inset_window = np.logical_and(t > (inset['xlims'][0] / t_plt['factor']), t < (inset['xlims'][1] / t_plt['factor'])) inset_ax.plot(t[inset_window] * t_plt['factor'], var[inset_window] * pltvar['factor'], linewidth=lw, linestyle=lines[j], color=colors[j]) # Add optional STIM-ON patches if patches[j]: (ybottom, ytop) = ax.get_ylim() la = [] color = '#8A8A8A' if greypatch else handle[0].get_color() for i in range(npatches): la.append(ax.axvspan(tpatch_on[i] * t_plt['factor'], tpatch_off[i] * t_plt['factor'], edgecolor='none', facecolor=color, alpha=0.2)) aliases[handle[0]] = la if inset is not None: cond_on = np.logical_and(tpatch_on > (inset['xlims'][0] / t_plt['factor']), tpatch_on < (inset['xlims'][1] / t_plt['factor'])) cond_off = np.logical_and(tpatch_off > (inset['xlims'][0] / t_plt['factor']), tpatch_off < (inset['xlims'][1] / t_plt['factor'])) cond_glob = np.logical_and(tpatch_on < (inset['xlims'][0] / t_plt['factor']), tpatch_off > (inset['xlims'][1] / t_plt['factor'])) cond_onoff = np.logical_or(cond_on, cond_off) cond = np.logical_or(cond_onoff, cond_glob) npatches_inset = np.sum(cond) for i in range(npatches_inset): inset_ax.add_patch(Rectangle((tpatch_on[cond][i] * t_plt['factor'], ybottom), (tpatch_off[cond][i] - tpatch_on[cond][i]) * t_plt['factor'], ytop - ybottom, color=color, alpha=0.1)) fig.tight_layout() # Optional operations on inset: if inset is not None: # Re-position inset axis axpos = ax.get_position() left, right, = rescale(inset['xcoords'], ax.get_xlim()[0], ax.get_xlim()[1], axpos.x0, axpos.x0 + axpos.width) bottom, top, = rescale(inset['ycoords'], ax.get_ylim()[0], ax.get_ylim()[1], axpos.y0, axpos.y0 + axpos.height) inset_ax.set_position([left, bottom, right - left, top - bottom]) for i in inset_ax.spines.values(): i.set_linewidth(2) # Materialize inset target region with contour frame ax.plot(inset['xlims'], [inset['ylims'][0]] * 2, linestyle='-', color='k') ax.plot(inset['xlims'], [inset['ylims'][1]] * 2, linestyle='-', color='k') ax.plot([inset['xlims'][0]] * 2, inset['ylims'], linestyle='-', color='k') ax.plot([inset['xlims'][1]] * 2, inset['ylims'], linestyle='-', color='k') # Link target and inset with dashed lines if possible if inset['xcoords'][1] < inset['xlims'][0]: ax.plot([inset['xcoords'][1], inset['xlims'][0]], [inset['ycoords'][0], inset['ylims'][0]], linestyle='--', color='k') ax.plot([inset['xcoords'][1], inset['xlims'][0]], [inset['ycoords'][1], inset['ylims'][1]], linestyle='--', color='k') elif inset['xcoords'][0] > inset['xlims'][1]: ax.plot([inset['xcoords'][0], inset['xlims'][1]], [inset['ycoords'][0], inset['ylims'][0]], linestyle='--', color='k') ax.plot([inset['xcoords'][0], inset['xlims'][1]], [inset['ycoords'][1], inset['ylims'][1]], linestyle='--', color='k') else: logger.warning('Inset x-coordinates intersect with those of target region') # Create interactive legend leg = ax.legend(loc=1, fontsize=fs, frameon=False) if blacklegend: for l in leg.get_lines(): l.set_color('k') if straightlegend: for l in leg.get_lines(): l.set_linestyle('-') interactive_legend = InteractiveLegend(ax.legend_, aliases) if showfig: plt.show() return fig def plotBatch(directory, filepaths, vars_dict=None, plt_show=True, plt_save=False, ask_before_save=True, fig_ext='png', tag='fig', fs=15, lw=2, title=True, show_patches=True): ''' Plot a figure with profiles of several specific NICE output variables, for several NICE simulations. :param positions: subplot indexes of each variable :param filepaths: list of full paths to output data files to be compared :param vars_dict: dict of lists of variables names to extract and plot together :param plt_show: boolean stating whether to show the created figures :param plt_save: boolean stating whether to save the created figures :param ask_before_save: boolean stating whether to show the created figures :param fig_ext: file extension for the saved figures :param tag: suffix added to the end of the figures name :param fs: labels font size :param lw: curves line width :param title: boolean stating whether to display a general title on the figures :param show_patches: boolean indicating whether to indicate periods of stimulation with colored rectangular patches ''' # Check validity of plot variables if vars_dict: yvars = list(sum(list(vars_dict.values()), [])) for key in yvars: if key not in pltvars: raise InputError('Unknown plot variable: "{}"'.format(key)) # Dictionary of neurons neurons_dict = getNeuronsDict() # Loop through data files for filepath in filepaths: # Get code from file name pkl_filename = ntpath.basename(filepath) filecode = pkl_filename[0:-4] # Retrieve sim type mo1 = rgxp.fullmatch(pkl_filename) mo2 = rgxp_mech.fullmatch(pkl_filename) if mo1: mo = mo1 elif mo2: mo = mo2 else: logger.error('Error: "%s" file does not match regexp pattern', pkl_filename) sys.exit(1) sim_type = mo.group(1) if sim_type not in ('MECH', 'ASTIM', 'ESTIM'): raise InputError('Invalid simulation type: {}'.format(sim_type)) # Load data logger.info('Loading data from "%s"', pkl_filename) with open(filepath, 'rb') as fh: frame = pickle.load(fh) df = frame['data'] meta = frame['meta'] # Extract variables logger.info('Extracting variables') t = df['t'].values states = df['states'].values nsamples = t.size # Initialize channel mechanism if sim_type in ['ASTIM', 'ESTIM']: neuron_name = mo.group(2) global neuron neuron = neurons_dict[neuron_name]() neuron_states = [df[sn].values for sn in neuron.states_names] Cm0 = neuron.Cm0 Qm0 = Cm0 * neuron.Vm0 * 1e-3 t_plt = pltvars['t_ms'] else: Cm0 = meta['Cm0'] Qm0 = meta['Qm0'] t_plt = pltvars['t_us'] # Initialize BLS if sim_type in ['MECH', 'ASTIM']: global bls Fdrive = meta['Fdrive'] a = meta['a'] bls = BilayerSonophore(a, Fdrive, Cm0, Qm0) # Determine patches location npatches, tpatch_on, tpatch_off = getPatchesLoc(t, states) # Adding onset to time vector if t_plt['onset'] > 0.0: tonset = np.array([-t_plt['onset'], -t[0] - t[1]]) t = np.hstack((tonset, t)) states = np.hstack((states, np.zeros(2))) # Determine variables to plot if not provided if not vars_dict: if sim_type == 'ASTIM': vars_dict = {'Z': ['Z'], 'Q_m': ['Qm']} elif sim_type == 'ESTIM': vars_dict = {'V_m': ['Vm']} elif sim_type == 'MECH': vars_dict = {'P_{AC}': ['Pac'], 'Z': ['Z'], 'n_g': ['ng']} if sim_type in ['ASTIM', 'ESTIM'] and hasattr(neuron, 'pltvars_scheme'): vars_dict.update(neuron.pltvars_scheme) labels = list(vars_dict.keys()) naxes = len(vars_dict) # Plotting if naxes == 1: _, ax = plt.subplots(figsize=(11, 4)) axes = [ax] else: _, axes = plt.subplots(naxes, 1, figsize=(11, min(3 * naxes, 9))) for i in range(naxes): ax = axes[i] for item in ['top', 'right']: ax.spines[item].set_visible(False) ax_pltvars = [pltvars[j] for j in vars_dict[labels[i]]] nvars = len(ax_pltvars) # X-axis if i < naxes - 1: ax.get_xaxis().set_ticklabels([]) else: ax.set_xlabel('${}\ ({})$'.format(t_plt['label'], t_plt['unit']), fontsize=fs) for tick in ax.xaxis.get_major_ticks(): tick.label.set_fontsize(fs) # Y-axis if ax_pltvars[0]['unit']: ax.set_ylabel('${}\ ({})$'.format(labels[i], ax_pltvars[0]['unit']), fontsize=fs) else: ax.set_ylabel('${}$'.format(labels[i]), fontsize=fs) if 'min' in ax_pltvars[0] and 'max' in ax_pltvars[0]: ax_min = min([ap['min'] for ap in ax_pltvars]) ax_max = max([ap['max'] for ap in ax_pltvars]) ax.set_ylim(ax_min, ax_max) ax.locator_params(axis='y', nbins=2) for tick in ax.yaxis.get_major_ticks(): tick.label.set_fontsize(fs) # Time series icolor = 0 for j in range(nvars): # Extract variable pltvar = ax_pltvars[j] if 'alias' in pltvar: var = eval(pltvar['alias']) elif 'key' in pltvar: var = df[pltvar['key']].values elif 'constant' in pltvar: var = eval(pltvar['constant']) * np.ones(nsamples) else: var = df[vars_dict[labels[i]][j]].values if var.size == t.size - 2: if pltvar['desc'] == 'membrane potential': var = np.hstack((np.array([neuron.Vm0] * 2), var)) else: var = np.hstack((np.array([var[0]] * 2), var)) # var = np.insert(var, 0, var[0]) # Plot variable if 'constant' in pltvar or pltvar['desc'] in ['net current']: ax.plot(t * t_plt['factor'], var * pltvar['factor'], '--', c='black', lw=lw, label='${}$'.format(pltvar['label'])) else: ax.plot(t * t_plt['factor'], var * pltvar['factor'], c='C{}'.format(icolor), lw=lw, label='${}$'.format(pltvar['label'])) icolor += 1 # Patches if show_patches == 1: (ybottom, ytop) = ax.get_ylim() for j in range(npatches): ax.axvspan(tpatch_on[j] * t_plt['factor'], tpatch_off[j] * t_plt['factor'], edgecolor='none', facecolor='#8A8A8A', alpha=0.2) # Legend if nvars > 1: ax.legend(fontsize=fs, loc=7, ncol=nvars // 4 + 1) # Title if title: if sim_type == 'ESTIM': if meta['DC'] == 1.0: fig_title = ESTIM_CW_title.format(neuron.name, meta['Astim'], meta['tstim'] * 1e3) else: fig_title = ESTIM_PW_title.format(neuron.name, meta['Astim'], meta['tstim'] * 1e3, meta['PRF'], meta['DC'] * 1e2) elif sim_type == 'ASTIM': if meta['DC'] == 1.0: fig_title = ASTIM_CW_title.format(neuron.name, Fdrive * 1e-3, meta['Adrive'] * 1e-3, meta['tstim'] * 1e3) else: fig_title = ASTIM_PW_title.format(neuron.name, Fdrive * 1e-3, meta['Adrive'] * 1e-3, meta['tstim'] * 1e3, meta['PRF'], meta['DC'] * 1e2) elif sim_type == 'MECH': fig_title = MECH_title.format(a * 1e9, Fdrive * 1e-3, meta['Adrive'] * 1e-3) axes[0].set_title(fig_title, fontsize=fs) plt.tight_layout() # Save figure if needed (automatic or checked) if plt_save: if ask_before_save: plt_filename = SaveFigDialog(directory, '{}_{}.{}'.format(filecode, tag, fig_ext)) else: plt_filename = '{}/{}_{}.{}'.format(directory, filecode, tag, fig_ext) if plt_filename: plt.savefig(plt_filename) logger.info('Saving figure as "{}"'.format(plt_filename)) plt.close() # Show all plots if needed if plt_show: plt.show() def plotGatingKinetics(neuron, fs=15): ''' Plot the voltage-dependent steady-states and time constants of activation and inactivation gates of the different ionic currents involved in a specific neuron's membrane. :param neuron: specific channel mechanism object :param fs: labels and title font size ''' # Input membrane potential vector Vm = np.linspace(-100, 50, 300) xinf_dict = {} taux_dict = {} logger.info('Computing %s neuron gating kinetics', neuron.name) names = neuron.states_names print(names) for xname in names: Vm_state = True # Names of functions of interest xinf_func_str = xname.lower() + 'inf' taux_func_str = 'tau' + xname.lower() alphax_func_str = 'alpha' + xname.lower() betax_func_str = 'beta' + xname.lower() # derx_func_str = 'der' + xname.upper() # 1st choice: use xinf and taux function if hasattr(neuron, xinf_func_str) and hasattr(neuron, taux_func_str): xinf_func = getattr(neuron, xinf_func_str) taux_func = getattr(neuron, taux_func_str) xinf = np.array([xinf_func(v) for v in Vm]) if isinstance(taux_func, float): taux = taux_func * np.ones(len(Vm)) else: taux = np.array([taux_func(v) for v in Vm]) # 2nd choice: use alphax and betax functions elif hasattr(neuron, alphax_func_str) and hasattr(neuron, betax_func_str): alphax_func = getattr(neuron, alphax_func_str) betax_func = getattr(neuron, betax_func_str) alphax = np.array([alphax_func(v) for v in Vm]) if isinstance(betax_func, float): betax = betax_func * np.ones(len(Vm)) else: betax = np.array([betax_func(v) for v in Vm]) taux = 1.0 / (alphax + betax) xinf = taux * alphax # # 3rd choice: use derX choice # elif hasattr(neuron, derx_func_str): # derx_func = getattr(neuron, derx_func_str) # xinf = brentq(lambda x: derx_func(neuron.Vm, x), 0, 1) else: Vm_state = False if not Vm_state: logger.error('no function to compute %s-state gating kinetics', xname) else: xinf_dict[xname] = xinf taux_dict[xname] = taux fig, axes = plt.subplots(2) fig.suptitle('{} neuron: gating dynamics'.format(neuron.name)) ax = axes[0] ax.get_xaxis().set_ticklabels([]) ax.set_ylabel('$X_{\infty}$', fontsize=fs) for xname in names: if xname in xinf_dict: ax.plot(Vm, xinf_dict[xname], lw=2, label='$' + xname + '_{\infty}$') ax.legend(fontsize=fs, loc=7) ax = axes[1] ax.set_xlabel('$V_m\ (mV)$', fontsize=fs) ax.set_ylabel('$\\tau_X\ (ms)$', fontsize=fs) for xname in names: if xname in taux_dict: ax.plot(Vm, taux_dict[xname] * 1e3, lw=2, label='$\\tau_{' + xname + '}$') ax.legend(fontsize=fs, loc=7) plt.show() def plotRateConstants(neuron, fs=15): ''' Plot the voltage-dependent activation and inactivation rate constants for each gate of all ionic currents involved in a specific neuron's membrane. :param neuron: specific channel mechanism object :param fs: labels and title font size ''' # Input membrane potential vector Vm = np.linspace(neuron.Vm0 - 10, 50, 100) alphax_dict = {} betax_dict = {} logger.info('Computing %s neuron gating kinetics', neuron.name) names = neuron.states_names for xname in names: Vm_state = True # Names of functions of interest xinf_func_str = xname.lower() + 'inf' taux_func_str = 'tau' + xname.lower() alphax_func_str = 'alpha' + xname.lower() betax_func_str = 'beta' + xname.lower() # 1st choice: use alphax and betax functions if hasattr(neuron, alphax_func_str) and hasattr(neuron, betax_func_str): alphax_func = getattr(neuron, alphax_func_str) betax_func = getattr(neuron, betax_func_str) alphax = np.array([alphax_func(v) for v in Vm]) betax = np.array([betax_func(v) for v in Vm]) # 2nd choice: use xinf and taux function elif hasattr(neuron, xinf_func_str) and hasattr(neuron, taux_func_str): xinf_func = getattr(neuron, xinf_func_str) taux_func = getattr(neuron, taux_func_str) xinf = np.array([xinf_func(v) for v in Vm]) taux = np.array([taux_func(v) for v in Vm]) alphax = xinf / taux betax = 1.0 / taux - alphax else: Vm_state = False if not Vm_state: logger.error('no function to compute %s-state gating kinetics', xname) else: alphax_dict[xname] = alphax betax_dict[xname] = betax naxes = len(alphax_dict) _, axes = plt.subplots(naxes, figsize=(11, min(3 * naxes, 9))) for i, xname in enumerate(alphax_dict.keys()): ax1 = axes[i] if i == 0: ax1.set_title('{} neuron: rate constants'.format(neuron.name)) if i == naxes - 1: ax1.set_xlabel('$V_m\ (mV)$', fontsize=fs) else: ax1.get_xaxis().set_ticklabels([]) ax1.set_ylabel('$\\alpha_{' + xname + '}\ (ms^{-1})$', fontsize=fs, color='C0') for label in ax1.get_yticklabels(): label.set_color('C0') ax1.plot(Vm, alphax_dict[xname] * 1e-3, lw=2) ax2 = ax1.twinx() ax2.set_ylabel('$\\beta_{' + xname + '}\ (ms^{-1})$', fontsize=fs, color='C1') for label in ax2.get_yticklabels(): label.set_color('C1') ax2.plot(Vm, betax_dict[xname] * 1e-3, lw=2, color='C1') plt.tight_layout() plt.show() def setGrid(n, ncolmax=3): ''' Determine number of rows and columns in figure grid, based on number of variables to plot. ''' if n <= ncolmax: return (1, n) else: return ((n - 1) // ncolmax + 1, ncolmax) def plotEffVars(neuron, Fdrive, a=32e-9, amps=None, charges=None, keys=None, fs=12, ncolmax=2): ''' Plot the profiles of effective variables of a specific neuron for a given frequency. For each variable, one line chart per amplitude is plotted, using charge as the input variable on the abscissa and a linear color code for the amplitude value. :param neuron: channel mechanism object :param Fdrive: acoustic drive frequency (Hz) :param a: sonophore diameter (m) :param amps: vector of amplitudes at which variables must be plotted (Pa) :param charges: vector of charges at which variables must be plotted (C/m2) :param keys: list of variables to plot :param fs: figure fontsize :param ncolmax: max number of columns on the figure :return: handle to the created figure ''' # Check lookup file existence lookup_file = '{}_lookups_a{:.1f}nm.pkl'.format(neuron.name, a * 1e9) lookup_path = '{}/{}'.format(getLookupDir(), lookup_file) if not os.path.isfile(lookup_path): raise InputError('Missing lookup file: "{}"'.format(lookup_file)) # Load coefficients with open(lookup_path, 'rb') as fh: lookups3D = pickle.load(fh) # Retrieve 1D inputs from lookup dictionary freqs = lookups3D.pop('f') amps_ref = lookups3D.pop('A') charges_ref = lookups3D.pop('Q') # Filter lookups keys if provided if keys is not None: lookups3D = {key: lookups3D[key] for key in keys} # Interpolate 3D lookups at US frequency lookups2D = itrpLookupsFreq(lookups3D, freqs, Fdrive) if 'V' in lookups2D: lookups2D['Vm'] = lookups2D.pop('V') keys[keys.index('V')] = 'Vm' # Define log-amplitude color code if amps is None: amps = amps_ref mymap = cm.get_cmap('Oranges') norm = matplotlib.colors.LogNorm(amps.min(), amps.max()) sm = cm.ScalarMappable(norm=norm, cmap=mymap) sm._A = [] # Plot logger.info('plotting') nrows, ncols = setGrid(len(lookups2D), ncolmax=ncolmax) xvar = pltvars['Qm'] if charges is None: charges = charges_ref Qbounds = np.array([charges.min(), charges.max()]) * xvar['factor'] fig, _ = plt.subplots(figsize=(3 * ncols, 1 * nrows), squeeze=False) for j, key in enumerate(keys): ax = plt.subplot2grid((nrows, ncols), (j // ncols, j % ncols)) for s in ['right', 'top']: ax.spines[s].set_visible(False) yvar = pltvars[key] if j // ncols == nrows - 1: ax.set_xlabel('$\\rm {}\ ({})$'.format(xvar['label'], xvar['unit']), fontsize=fs) ax.set_xticks(Qbounds) else: ax.set_xticks([]) ax.spines['bottom'].set_visible(False) ax.xaxis.set_label_coords(0.5, -0.1) ax.yaxis.set_label_coords(-0.02, 0.5) for item in ax.get_xticklabels() + ax.get_yticklabels(): item.set_fontsize(fs) ymin = np.inf ymax = -np.inf y0 = np.squeeze(interp2d(amps_ref, charges_ref, lookups2D[key].T)(0, charges)) # Plot effective variable for each selected amplitude for Adrive in amps: y = np.squeeze(interp2d(amps_ref, charges_ref, lookups2D[key].T)(Adrive, charges)) if 'alpha' in key or 'beta' in key: y[y > y0.max() * 2] = np.nan ax.plot(charges * xvar['factor'], y * yvar['factor'], c=sm.to_rgba(Adrive)) ymin = min(ymin, y.min()) ymax = max(ymax, y.max()) # Plot reference variable ax.plot(charges * xvar['factor'], y0 * yvar['factor'], '--', c='k') ymax = max(ymax, y0.max()) ymin = min(ymin, y0.min()) # Set axis y-limits if 'alpha' in key or 'beta' in key: ymax = y0.max() * 2 ylim = [ymin * yvar['factor'], ymax * yvar['factor']] if key == 'ng': ylim = [np.floor(ylim[0] * 1e2) / 1e2, np.ceil(ylim[1] * 1e2) / 1e2] else: factor = 1 / np.power(10, np.floor(np.log10(ylim[1]))) print(key, ylim[1], factor) ylim = [np.floor(ylim[0] * factor) / factor, np.ceil(ylim[1] * factor) / factor] dy = ylim[1] - ylim[0] ax.set_yticks(ylim) ax.set_ylim(ylim) # ax.set_ylim([ylim[0] - 0.05 * dy, ylim[1] + 0.05 * dy]) # Annotate variable and unit xlim = ax.get_xlim() if np.argmax(y0) < np.argmin(y0): xtext = xlim[0] + 0.6 * (xlim[1] - xlim[0]) else: xtext = xlim[0] + 0.01 * (xlim[1] - xlim[0]) if key in ['Vm', 'ng']: ytext = ylim[0] + 0.85 * dy else: ytext = ylim[0] + 0.15 * dy ax.text(xtext, ytext, '$\\rm {}\ ({})$'.format(yvar['label'], yvar['unit']), fontsize=fs) fig.suptitle('{} neuron: original vs. effective variables @ {:.0f} kHz'.format( neuron.name, Fdrive * 1e-3)) # Plot colorbar fig.subplots_adjust(left=0.10, bottom=0.05, top=0.9, right=0.85) cbarax = fig.add_axes([0.87, 0.05, 0.04, 0.85]) fig.colorbar(sm, cax=cbarax) cbarax.set_ylabel('amplitude (Pa)', fontsize=fs) for item in cbarax.get_yticklabels(): item.set_fontsize(fs) return fig def plotActivationMap(DCs, amps, actmap, FRlims, title=None, Ascale='log', FRscale='log', fs=8): ''' Plot a neuron's activation map over the amplitude x duty cycle 2D space. :param DCs: duty cycle vector :param amps: amplitude vector :param actmap: 2D activation matrix :param FRlims: lower and upper bounds of firing rate color-scale :param title: figure title :param Ascale: scale to use for the amplitude dimension ('lin' or 'log') :param FRscale: scale to use for the firing rate coloring ('lin' or 'log') :param fs: fontsize to use for the title and labels :return: 3-tuple with the handle to the generated figure and the mesh x and y coordinates ''' # Check firing rate bounding minFR, maxFR = (actmap[actmap > 0].min(), actmap.max()) logger.info('FR range: %.0f - %.0f Hz', minFR, maxFR) if minFR < FRlims[0]: logger.warning('Minimal firing rate (%.0f Hz) is below defined lower bound (%.0f Hz)', minFR, FRlims[0]) if maxFR > FRlims[1]: logger.warning('Maximal firing rate (%.0f Hz) is above defined upper bound (%.0f Hz)', maxFR, FRlims[1]) # Plot activation map if FRscale == 'lin': norm = matplotlib.colors.Normalize(*FRlims) elif FRscale == 'log': norm = matplotlib.colors.LogNorm(*FRlims) fig, ax = plt.subplots(figsize=cm2inch(8, 5.8)) fig.subplots_adjust(left=0.15, bottom=0.15, right=0.8, top=0.92) if title is not None: ax.set_title(title, fontsize=fs) if Ascale == 'log': ax.set_yscale('log') ax.set_xlabel('Duty cycle (%)', fontsize=fs, labelpad=-0.5) ax.set_ylabel('Amplitude (kPa)', fontsize=fs) ax.set_xlim(np.array([DCs.min(), DCs.max()]) * 1e2) for item in ax.get_xticklabels() + ax.get_yticklabels(): item.set_fontsize(fs) xedges = computeMeshEdges(DCs) yedges = computeMeshEdges(amps, scale='log') actmap[actmap == -1] = np.nan actmap[actmap == 0] = 1e-3 cmap = plt.get_cmap('viridis') cmap.set_bad('silver') cmap.set_under('k') ax.pcolormesh(xedges * 1e2, yedges * 1e-3, actmap, cmap=cmap, norm=norm) # Plot firing rate colorbar sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) sm._A = [] pos1 = ax.get_position() # get the map axis position cbarax = fig.add_axes([pos1.x1 + 0.02, pos1.y0, 0.03, pos1.height]) fig.colorbar(sm, cax=cbarax) cbarax.set_ylabel('Firing rate (Hz)', fontsize=fs) for item in cbarax.get_yticklabels(): item.set_fontsize(fs) return (fig, xedges, yedges) def plotDualMaxMap(DCs, amps, maxmap, factor, actmap, title, lbl='xy', Ascale='log', fs=8): ''' Plot a variable maximum map over the amplitude x duty cycle 2D space, with different color codes for the sub and supra-threshold regions. :param DCs: duty cycle vector :param amps: amplitude vector :param maxmap: 2D variable maximum matrix :param factor: unit factor to use for the colorbars labels :param actmap: 2D activation matrix :param title: figure title :param lbl: indicates whether to label the x and y axes :param Ascale: scale to use for the amplitude dimension ('lin' or 'log') :param fs: fontsize to use for the title and labels :return: a handle to the generated figure ''' # Split variable max map into sub-threshold and supra-threshold max maps maxmap_sub, maxmap_supra = maxmap.copy(), maxmap.copy() maxmap_sub[actmap >= 0] = np.nan maxmap_supra[actmap < 0] = np.nan # Plot dual max map fig, ax = plt.subplots(figsize=cm2inch(8, 5.8)) fig.subplots_adjust(left=0.15, bottom=0.15, right=0.8, top=0.92) ax.set_title(title, fontsize=fs) if Ascale == 'log': ax.set_yscale('log') if 'x' in lbl: ax.set_xlabel('Duty cycle (%)', fontsize=fs) else: ax.set_xticklabels([]) if 'y' in lbl: ax.set_ylabel('Amplitude (kPa)', fontsize=fs) else: ax.set_yticklabels([]) for item in ax.get_xticklabels() + ax.get_yticklabels(): item.set_fontsize(fs) xedges = computeMeshEdges(DCs) yedges = computeMeshEdges(amps, scale=Ascale) # Plot the 2 corresponding colorbars sm_sub = ax.pcolormesh(xedges * 1e2, yedges * 1e-3, maxmap_sub * factor, cmap='Blues') sm_supra = ax.pcolormesh(xedges * 1e2, yedges * 1e-3, maxmap_supra * factor, cmap='Reds') pos1 = ax.get_position() # get the map axis position height = (pos1.height - 0.05) / 2.0 cbarax_sub = fig.add_axes([pos1.x1 + 0.02, pos1.y0, 0.03, height]) cbar_sub = fig.colorbar(sm_sub, cax=cbarax_sub, format='%.1f') cbar_sub.set_ticks(np.array([np.nanmin(maxmap_sub), np.nanmax(maxmap_sub)]) * factor) cbarax_supra = fig.add_axes([pos1.x1 + 0.02, pos1.y1 - height, 0.03, height]) cbar_supra = fig.colorbar(sm_supra, cax=cbarax_supra, format='%.1f') cbar_supra.set_ticks(np.array([np.nanmin(maxmap_supra), np.nanmax(maxmap_supra)]) * factor) for item in cbarax_sub.get_yticklabels() + cbarax_supra.get_yticklabels(): item.set_fontsize(fs) return fig def plotRawTrace(fpath, key, ybounds): ''' Plot the raw signal of a given variable within specified bounds. :param foath: full path to the data file :param key: key to the target variable :param ybounds: y-axis bounds :return: handle to the generated figure ''' # Check file existence fname = ntpath.basename(fpath) if not os.path.isfile(fpath): raise InputError('Error: "{}" file does not exist'.format(fname)) # Load data logger.debug('Loading data from "%s"', fname) with open(fpath, 'rb') as fh: frame = pickle.load(fh) df = frame['data'] t = df['t'].values y = df[key].values * pltvars[key]['factor'] Δy = y.max() - y.min() logger.info('d%s = %.1f %s', key, Δy, pltvars[key]['unit']) # Plot trace fig, ax = plt.subplots(figsize=cm2inch(12.5, 5.8)) fig.canvas.set_window_title(fname) for s in ['top', 'bottom', 'left', 'right']: ax.spines[s].set_visible(False) ax.set_xticks([]) ax.set_yticks([]) ax.set_ylim(ybounds) ax.plot(t, y, color='k', linewidth=1) fig.tight_layout() return fig def plotTraces(fpath, keys, tbounds): ''' Plot the raw signal of sevral variables within specified bounds. :param foath: full path to the data file :param key: key to the target variable :param tbounds: x-axis bounds :return: handle to the generated figure ''' # Check file existence fname = ntpath.basename(fpath) if not os.path.isfile(fpath): raise InputError('Error: "{}" file does not exist'.format(fname)) # Load data logger.debug('Loading data from "%s"', fname) with open(fpath, 'rb') as fh: frame = pickle.load(fh) df = frame['data'] t = df['t'].values * 1e3 # Plot trace fs = 8 fig, ax = plt.subplots(figsize=cm2inch(7, 3)) fig.canvas.set_window_title(fname) plt.subplots_adjust(left=0.2, bottom=0.2, right=0.95, top=0.95) for s in ['top', 'right']: ax.spines[s].set_visible(False) for s in ['bottom', 'left']: ax.spines[s].set_position(('axes', -0.03)) ax.spines[s].set_linewidth(2) ax.yaxis.set_tick_params(width=2) # ax.spines['bottom'].set_linewidth(2) ax.set_xlim(tbounds) ax.set_xticks([]) ymin = np.nan ymax = np.nan dt = tbounds[1] - tbounds[0] ax.set_xlabel('{}s'.format(si_format(dt * 1e-3, space=' ')), fontsize=fs) ax.set_ylabel('mV - $\\rm nC/cm^2$', fontsize=fs, labelpad=-15) colors = {'Vm': 'darkgrey', 'Qm': 'k'} for key in keys: y = df[key].values * pltvars[key]['factor'] ymin = np.nanmin([ymin, y.min()]) ymax = np.nanmax([ymax, y.max()]) # if key == 'Qm': # y0 = y[0] # ax.plot(t, y0 * np.ones(t.size), '--', color='k', linewidth=1) Δy = y.max() - y.min() logger.info('d%s = %.1f %s', key, Δy, pltvars[key]['unit']) ax.plot(t, y, color=colors[key], linewidth=1) ax.yaxis.set_major_formatter(FormatStrFormatter('%.0f')) # ax.set_yticks([ymin, ymax]) ax.set_ylim([-200, 100]) ax.set_yticks([-200, 100]) for item in ax.get_yticklabels(): item.set_fontsize(fs) # fig.tight_layout() return fig -def plotSignals(t, signals, states=None, ax=None, onset=None, lbls=None, fs=10): +def plotSignals(t, signals, states=None, ax=None, onset=None, lbls=None, fs=10, cmode='qual'): ''' Plot several signals on one graph. :param t: time vector :param signals: list of signal vectors :param states (optional): stimulation state vector :param ax (optional): handle to figure axis :param onset (optional): onset to add to signals on graph :param lbls (optional): list of legend labels :param fs (optional): font size to use on graph + :param cmode: color mode ('seq' for sequentiual or 'qual' for qualitative) :return: figure handle (if no axis provided as input) ''' + # If no axis provided, create figure if ax is None: fig, ax = plt.subplots(figsize=(8, 3)) argout = fig else: argout = None # Set axis aspect ax.spines['top'].set_visible(False) ax.spines['right'].set_visible(False) for item in ax.get_xticklabels() + ax.get_xticklabels(): item.set_fontsize(fs) # Compute number of signals nsignals = len(signals) + # Adapt labels for sequential color mode + if cmode == 'seq' and lbls is not None: + lbls[1:-1] = ['.'] * (nsignals - 2) + # Add stimulation patches if states provided if states is not None: npatches, tpatch_on, tpatch_off = getPatchesLoc(t, states) for i in range(npatches): ax.axvspan(tpatch_on[i], tpatch_off[i], edgecolor='none', facecolor='#8A8A8A', alpha=0.2) # Add onset of provided if onset is not None: t0, y0 = onset t = np.hstack((np.array([t0, 0.]), t)) signals = np.hstack((np.ones((nsignals, 2)) * y0, signals)) + # Determine colorset + nlevels = nsignals + if cmode == 'seq': + norm = matplotlib.colors.Normalize(0, nlevels - 1) + sm = cm.ScalarMappable(norm=norm, cmap=cm.get_cmap('viridis')) + sm._A = [] + colors = [sm.to_rgba(i) for i in range(nlevels)] + elif cmode == 'qual': + nlevels_max = 10 + if nlevels > nlevels_max: + raise Warning('Number of signals higher than number of color levels') + colors = ['C{}'.format(i) for i in range(nlevels)] + else: + raise InputError('Unknown color mode') + # Plot signals for i, var in enumerate(signals): - ax.plot(t, var, label=lbls[i] if lbls is not None else None, c='C{}'.format(i)) + ax.plot(t, var, label=lbls[i] if lbls is not None else None, c=colors[i]) # Add legend if lbls is not None: ax.legend(fontsize=fs, frameon=False) return argout diff --git a/README.md b/README.md index 28aecaf..fded99e 100644 --- a/README.md +++ b/README.md @@ -1,70 +1,65 @@ -Description -============ +## Description Python implementation of the **multi-Scale Optimized Neuronal Intramembrane Cavitation** (SONIC) model to compute individual neural responses to acoustic stimuli, as predicted by the *intramembrane cavitation* hypothesis. This is an optmimized variant of the original **Neuronal Intramembrane Cavitation Excitation** (NICE) model introduced by Plaksin et. al in 2014. This package contains several core modules: - - **bls** defines the underlying biomechanical model of intramembrane cavitation (**BilayerSonophore** class), and provides an integration method to predict compute the mechanical oscillations of the plasma membrane subject to a periodic acoustic perturbation. - - **solvers** contains a simple solver for electrical stimuli (**SolverElec** class) as well as a tailored solver for acoustic stimuli (**SolverUS** class). The latter directly inherits from the BilayerSonophore class upon instantiation, and is hooked to a specific "channel mechanism" in order to link the mechanical model to an electrical model of membrane dynamics. It also provides several integration methods (detailed below) to compute the behaviour of the full electro-mechanical model subject to a continuous or pulsed ultrasonic stimulus. - - **neurons** contains the definitions of the different channels mechanisms inherent to specific neurons, including several types of **cortical** and **thalamic** neurons. - - **plt** defines plotting utilities to load results of several simulations and display/compare temporal profiles of multiple variables of interest across simulations. - - **utils** defines generic utilities used across the different modules +- **bls** defines the underlying biomechanical model of intramembrane cavitation (**BilayerSonophore** class), and provides an integration method to predict compute the mechanical oscillations of the plasma membrane subject to a periodic acoustic perturbation. +- **solvers** contains a simple solver for electrical stimuli (**SolverElec** class) as well as a tailored solver for acoustic stimuli (**SolverUS** class). The latter directly inherits from the BilayerSonophore class upon instantiation, and is hooked to a specific "channel mechanism" in order to link the mechanical model to an electrical model of membrane dynamics. It also provides several integration methods (detailed below) to compute the behaviour of the full electro-mechanical model subject to a continuous or pulsed ultrasonic stimulus. +- **neurons** contains the definitions of the different channels mechanisms inherent to specific neurons, including several types of **cortical** and **thalamic** neurons. +- **plt** defines plotting utilities to load results of several simulations and display/compare temporal profiles of multiple variables of interest across simulations. +- **utils** defines generic utilities used across the different modules The **SolverUS** class incorporates optimized numerical integration methods to perform dynamic simulations of the model subject to acoustic perturbation, and compute the evolution of its mechanical and electrical variables: - - a **classic** method that solves all variables for the entire duration of the simulation. This method uses very small time steps and is computationally expensive (simulation time: several hours) - - a **hybrid** method (initially developed by Plaskin et al.) in which integration is performed in consecutive “slices” of time, during which the full system is solved until mechanical stabilization, at which point the electrical system is solely solved with predicted mechanical variables until the end of the slice. This method is more efficient (simulation time: several minutes) and provides accurate results. - - a newly developed **effective** method that neglects the high amplitude oscillations of mechanical and electrical variables during each acoustic cycle, to instead grasp the net effect of the acoustic stimulus on the electrical system. To do so, the sole electrical system is solved using pre-computed coefficients that depend on membrane charge and acoustic amplitude. This method allows to run simulations of the electrical system in only a few seconds, with very accurate results of the net membrane charge density evolution. +- a **classic** method that solves all variables for the entire duration of the simulation. This method uses very small time steps and is computationally expensive (simulation time: several hours) +- a **hybrid** method (initially developed by Plaskin et al.) in which integration is performed in consecutive “slices” of time, during which the full system is solved until mechanical stabilization, at which point the electrical system is solely solved with predicted mechanical variables until the end of the slice. This method is more efficient (simulation time: several minutes) and provides accurate results. +- a newly developed **effective** method that neglects the high amplitude oscillations of mechanical and electrical variables during each acoustic cycle, to instead grasp the net effect of the acoustic stimulus on the electrical system. To do so, the sole electrical system is solved using pre-computed coefficients that depend on membrane charge and acoustic amplitude. This method allows to run simulations of the electrical system in only a few seconds, with very accurate results of the net membrane charge density evolution. This package is meant to be easy to use as a predictive and comparative tool for researchers investigating ultrasonic and/or electrical neuro-stimulation experimentally. -Installation -================== +## Installation Install Python 3 if not already done. Open a terminal. Activate a Python3 environment if needed, e.g. on the tnesrv5 machine: - source /opt/apps/anaconda3/bin activate +```$ source /opt/apps/anaconda3/bin activate``` Check that the appropriate version of pip is activated: - pip --version +```$ pip --version``` Go to the package directory (where the setup.py file is located) and install it: - cd - pip install -e . +``` +$ cd +$ pip install -e . +``` *PySONIC* and all its dependencies will be installed. -Usage -======= +## Usage -Command line scripts ---------------------- +### Command line scripts To run single simulations of a given point-neuron model under specific stimulation parameters, you can use the `ASTIM_run.py` and `ESTIM_run.py` command-line scripts provided by the package. For instance, to simulate a regular-spiking neuron under continuous wave ultrasonic stimulation at 500kHz and 100kPa, for 150 ms: - python ASTIM_run.py -n=RS -f=500 -A=100 -t=150 +```python run_ASTIM.py -n=RS -f=500 -A=100 -t=150``` Similarly, to simulate the electrical stimulation of a thalamo-cortical neuron at 10 mA/m2 for 150 ms: - python ESTIM_run.py -n=TC -A=10 -t=150 +```python run_ESTIM.py -n=TC -A=10 -t=150``` The simulation results will be save in an output PKL file in the current working directory. To view these results, you can use the dedicated -Batch scripts ---------------- - -To run a batch of simulations on different neuron types and spanning ranges of several stimulation parameters, you can run the `ASTIM_batch.py` and `ESTIM_batch.py` scripts. To do so, simply modify the **stim_params** and **neurons** variables with your own neuron types and parameter sweeps, and then run the scripts (without command-line arguments). - +### Batch scripts +To run a batch of simulations on different neuron types and spanning ranges of several stimulation parameters, you can run the `batch_ASTIM.py` and `batch_ESTIM.py` scripts. To do so, simply modify the **stim_params** and **neurons** variables with your own neuron types and parameter sweeps, and then run the scripts (without command-line arguments).