
			import inspect
			import pandas as pd
			import numpy as np
			import matplotlib.pyplot as plt
			from matplotlib import cm, colors

			from ..postpro import getFixedPoints
			from ..core import NeuronalBilayerSonophore
			from .pltutils import *
			from ..utils import logger
			from ..batches import runBatch


			def plotVarQSSDynamics(neuron, a, Fdrive, Adrive, charges, varname, varrange, fs=12):
			''' Plot the QSS-approximated derivative of a specific variable as function of
			the variable itself, as well as equilibrium values, for various membrane
			charge densities at a given acoustic amplitude.

			:param neuron: neuron object
			:param a: sonophore radius (m)
			:param Fdrive: US frequency (Hz)
			:param Adrive: US amplitude (Pa)
			:param charges: charge density vector (C/m2)
			:param varname: name of variable to plot
			:param varrange: range over which to compute the derivative
			:return: figure handle
			'''

			# Extract information about variable to plot
			pltvar = neuron.getPltVars()[varname]

			# Get methods to compute derivative and steady-state of variable of interest
			derX_func = getattr(neuron, 'der{}{}'.format(varname[0].upper(), varname[1:]))
			Xinf_func = getattr(neuron, '{}inf'.format(varname))
			derX_args = inspect.getargspec(derX_func)[0][1:]
			Xinf_args = inspect.getargspec(Xinf_func)[0][1:]

			# Get dictionary of charge and amplitude dependent QSS variables
			nbls = NeuronalBilayerSonophore(a, neuron, Fdrive)
			_, Qref, lookups, QSS = nbls.quasiSteadyStates(
			Fdrive, amps=Adrive, charges=charges, squeeze_output=True)
			df = QSS
			df['Vm'] = lookups['V']

			# Create figure
			fig, ax = plt.subplots(figsize=(6, 4))
			ax.set_title('{} neuron - QSS {} dynamics @ {:.2f} kPa'.format(
			neuron.name, pltvar['desc'], Adrive * 1e-3), fontsize=fs)
			ax.set_xscale('log')
			for key in ['top', 'right']:
			ax.spines[key].set_visible(False)
			ax.set_xlabel('$\\rm {}\ ({})$'.format(pltvar['label'], pltvar.get('unit', '')),
			fontsize=fs)
			ax.set_ylabel('$\\rm QSS\ d{}/dt\ ({}/s)$'.format(pltvar['label'], pltvar.get('unit', '1')),
			fontsize=fs)
			ax.set_ylim(-40, 40)
			ax.axhline(0, c='k', linewidth=0.5)

			y0_str = '{}0'.format(varname)
			if hasattr(neuron, y0_str):
			ax.axvline(getattr(neuron, y0_str) * pltvar.get('factor', 1),
			label=y0_str, c='k', linewidth=0.5)

			# For each charge value
			icolor = 0
			for j, Qm in enumerate(charges):
			lbl = 'Q = {:.0f} nC/cm2'.format(Qm * 1e5)

			# Compute variable derivative as a function of its value, as well as equilibrium value,
			# keeping other variables at quasi steady-state
			derX_inputs = [varrange if arg == varname else df[arg][j] for arg in derX_args]
			Xinf_inputs = [df[arg][j] for arg in Xinf_args]
			dX_QSS = neuron.derCai(*derX_inputs)
			Xeq_QSS = neuron.Caiinf(*Xinf_inputs)

			# Plot variable derivative and its root as a function of the variable itself
			c = 'C{}'.format(icolor)
			ax.plot(varrange * pltvar.get('factor', 1), dX_QSS * pltvar.get('factor', 1), c=c, label=lbl)
			ax.axvline(Xeq_QSS * pltvar.get('factor', 1), linestyle='--', c=c)
			icolor += 1

			ax.legend(frameon=False, fontsize=fs - 3)
			for item in ax.get_xticklabels() + ax.get_yticklabels():
			item.set_fontsize(fs)
			fig.tight_layout()

			fig.canvas.set_window_title('{}_QSS_{}_dynamics_{:.2f}kPa'.format(
			neuron.name, varname, Adrive * 1e-3))

			return fig


			def plotQSSvars(neuron, a, Fdrive, Adrive, fs=12):
			''' Plot effective membrane potential, quasi-steady states and resulting membrane currents
			as a function of membrane charge density, for a given acoustic amplitude.

			:param neuron: neuron object
			:param a: sonophore radius (m)
			:param Fdrive: US frequency (Hz)
			:param Adrive: US amplitude (Pa)
			:return: figure handle
			'''

			# Get neuron-specific pltvars
			pltvars = neuron.getPltVars()

			# Compute neuron-specific charge and amplitude dependent QS states at this amplitude
			nbls = NeuronalBilayerSonophore(a, neuron, Fdrive)
			_, Qref, lookups, QSS = nbls.quasiSteadyStates(Fdrive, amps=Adrive, squeeze_output=True)
			Vmeff = lookups['V']

			# Compute QSS currents
			currents = neuron.currents(Vmeff, np.array([QSS[k] for k in neuron.states]))
			iNet = sum(currents.values())

			# Compute fixed points in dQdt profile
			dQdt = -iNet
			Q_SFPs = getFixedPoints(Qref, dQdt, filter='stable')
			Q_UFPs = getFixedPoints(Qref, dQdt, filter='unstable')

			# Extract dimensionless states
			norm_QSS = {}
			for x in neuron.states:
			if 'unit' not in pltvars[x]:
			norm_QSS[x] = QSS[x]

			# Create figure
			fig, axes = plt.subplots(3, 1, figsize=(7, 9))
			axes[-1].set_xlabel('$\\rm Q_m\ (nC/cm^2)$', fontsize=fs)
			for ax in axes:
			for skey in ['top', 'right']:
			ax.spines[skey].set_visible(False)
			for item in ax.get_xticklabels() + ax.get_yticklabels():
			item.set_fontsize(fs)
			for item in ax.get_xticklabels(minor=True):
			item.set_visible(False)
			figname = '{} neuron - QSS dynamics @ {:.2f} kPa'.format(neuron.name, Adrive * 1e-3)
			fig.suptitle(figname, fontsize=fs)

			# Subplot: Vmeff
			ax = axes[0]
			ax.set_ylabel('$V_m^*$ (mV)', fontsize=fs)
			ax.plot(Qref * 1e5, Vmeff, color='k')
			ax.axhline(neuron.Vm0, linewidth=0.5, color='k')

			# Subplot: dimensionless quasi-steady states
			cset = plt.get_cmap('Dark2').colors + plt.get_cmap('tab10').colors
			ax = axes[1]
			ax.set_ylabel('QSS gating variables (-)', fontsize=fs)
			ax.set_yticks([0, 0.5, 1])
			ax.set_ylim([-0.05, 1.05])
			for i, (label, QS_state) in enumerate(norm_QSS.items()):
			ax.plot(Qref * 1e5, QS_state, label=label, c=cset[i])

			# Subplot: currents
			ax = axes[2]
			cset = plt.get_cmap('tab10').colors
			ax.set_ylabel('QSS currents ($\\rm A/m^2$)', fontsize=fs)
			for i, (k, I) in enumerate(currents.items()):
			ax.plot(Qref * 1e5, -I * 1e-3, '--', c=cset[i],
			label='$\\rm -{}$'.format(neuron.getPltVars()[k]['label']))
			ax.plot(Qref * 1e5, -iNet * 1e-3, color='k', label='$\\rm -I_{Net}$')
			ax.axhline(0, color='k', linewidth=0.5)
			if Q_SFPs.size > 0:
			ax.plot(Q_SFPs * 1e5, np.zeros(Q_SFPs.size), 'o', c='k', markersize=5, zorder=2)
			if Q_SFPs.size > 0:
			ax.plot(Q_UFPs * 1e5, np.zeros(Q_UFPs.size), 'o', c='k', markersize=5, mfc='none', zorder=2)

			fig.tight_layout()
			fig.subplots_adjust(right=0.8)
			for ax in axes[1:]:
			ax.legend(loc='center right', fontsize=fs, frameon=False, bbox_to_anchor=(1.3, 0.5))
			for ax in axes[:-1]:
			ax.set_xticklabels([])

			fig.canvas.set_window_title(
			'{}_QSS_states_vs_Qm_{:.2f}kPa'.format(neuron.name, Adrive * 1e-3))

			return fig


			def plotQSSVarVsAmp(neuron, a, Fdrive, varname, amps=None, DC=1.,
			fs=12, cmap='viridis', yscale='lin', zscale='lin'):
			''' Plot a specific QSS variable (state or current) as a function of
			membrane charge density, for various acoustic amplitudes.

			:param neuron: neuron object
			:param a: sonophore radius (m)
			:param Fdrive: US frequency (Hz)
			:param amps: US amplitudes (Pa)
			:param DC: duty cycle (-)
			:param varname: extraction key for variable to plot
			:return: figure handle
			'''

			# Determine stimulation modality
			if a is None and Fdrive is None:
			stim_type = 'elec'
			a = 32e-9
			Fdrive = 500e3
			else:
			stim_type = 'US'

			# Extract information about variable to plot
			pltvar = neuron.getPltVars()[varname]
			Qvar = neuron.getPltVars()['Qm']
			Afactor = {'US': 1e-3, 'elec': 1.}[stim_type]
			# Q_SFPs = []
			# Q_UFPs = []

			log = 'plotting {} neuron QSS {} vs. amp for {} stimulation @ {:.0f}% DC'.format(
			neuron.name, varname, stim_type, DC * 1e2)
			logger.info(log)

			nbls = NeuronalBilayerSonophore(a, neuron, Fdrive)

			# Get reference dictionaries for zero amplitude
			_, Qref, lookups0, QSS0 = nbls.quasiSteadyStates(Fdrive, amps=0., squeeze_output=True)
			Vmeff0 = lookups0['V']
			if stim_type == 'elec': # if E-STIM case, compute steady states with constant capacitance
			Vmeff0 = Qref / neuron.Cm0 * 1e3
			QSS0 = neuron.steadyStates(Vmeff0)
			df0 = QSS0
			df0['Vm'] = Vmeff0

			# Create figure
			fig, ax = plt.subplots(figsize=(6, 4))
			title = '{} neuron - {}steady-state {}'.format(
			neuron.name, 'quasi-' if amps is not None else '', pltvar['desc'])
			if amps is not None:
			title += '\nvs. {} amplitude @ {:.0f}% DC'.format(stim_type, DC * 1e2)
			ax.set_title(title, fontsize=fs)
			ax.set_xlabel('$\\rm {}\ ({})$'.format(Qvar['label'], Qvar['unit']), fontsize=fs)
			ax.set_ylabel('$\\rm QSS\ {}\ ({})$'.format(pltvar['label'], pltvar.get('unit', '')),
			fontsize=fs)
			if yscale == 'log':
			ax.set_yscale('log')
			for key in ['top', 'right']:
			ax.spines[key].set_visible(False)

			# Plot y-variable reference line, if any
			y0 = None
			y0_str = '{}0'.format(varname)
			if hasattr(neuron, y0_str):
			y0 = getattr(neuron, y0_str) * pltvar.get('factor', 1)
			elif varname in neuron.getCurrentsNames() + ['iNet', 'dQdt']:
			y0 = 0.
			y0_str = ''
			if y0 is not None:
			ax.axhline(y0, label=y0_str, c='k', linewidth=0.5)

			# Plot reference QSS profile of variable as a function of charge density
			var0 = extractPltVar(
			neuron, pltvar, pd.DataFrame({k: df0[k] for k in df0.keys()}), name=varname)
			ax.plot(Qref * Qvar['factor'], var0, '--', c='k', zorder=1,
			label='$\\rm A_{{{}}}=0$'.format(stim_type))
			# if varname == 'dQdt':
			# Q_SFPs += getFixedPoints(Qref, var0, filter='stable').tolist()
			# Q_UFPs += getFixedPoints(Qref, var0, filter='unstable').tolist()

			# Define color code
			mymap = plt.get_cmap(cmap)
			zref = amps * Afactor
			if zscale == 'lin':
			norm = colors.Normalize(zref.min(), zref.max())
			elif zscale == 'log':
			norm = colors.LogNorm(zref.min(), zref.max())
			sm = cm.ScalarMappable(norm=norm, cmap=mymap)
			sm._A = []

			# Get amplitude-dependent QSS dictionary
			if stim_type == 'US':
			# Get dictionary of charge and amplitude dependent QSS variables
			_, Qref, lookups, QSS = nbls.quasiSteadyStates(
			Fdrive, amps=amps, DCs=DC, squeeze_output=True)
			df = QSS
			df['Vm'] = lookups['V']
			else:
			# Repeat zero-amplitude QSS dictionary for all amplitudes
			df = {k: np.tile(df0[k], (amps.size, 1)) for k in df0}

			# Plot QSS profiles for various amplitudes
			for i, A in enumerate(amps):
			var = extractPltVar(
			neuron, pltvar, pd.DataFrame({k: df[k][i] for k in df.keys()}), name=varname)
			if varname == 'dQdt' and stim_type == 'elec':
			var += A * DC * pltvar['factor']
			ax.plot(Qref * Qvar['factor'], var, c=sm.to_rgba(A * Afactor), zorder=0)
			# if varname == 'dQdt':
			# # mark eq. point if starting point provided, otherwise mark all FPs
			# Q_SFPs += getFixedPoints(Qref, var, filter='stable').tolist()
			# Q_UFPs += getFixedPoints(Qref, var, filter='unstable').tolist()

			# # Plot fixed-points, if any
			# if len(Q_SFPs) > 0:
			# ax.plot(np.array(Q_SFPs) * Qvar['factor'], np.zeros(len(Q_SFPs)), 'o', c='k',
			# markersize=5, zorder=2)
			# if len(Q_UFPs) > 0:
			# ax.plot(np.array(Q_UFPs) * Qvar['factor'], np.zeros(len(Q_UFPs)), 'x', c='k',
			# markersize=5, zorder=2)

			# Add legend and adjust layout
			ax.legend(frameon=False, fontsize=fs)
			for item in ax.get_xticklabels() + ax.get_yticklabels():
			item.set_fontsize(fs)
			fig.tight_layout()
			fig.subplots_adjust(bottom=0.15, top=0.9, right=0.80, hspace=0.5)

			# Plot amplitude colorbar
			if amps is not None:
			cbarax = fig.add_axes([0.85, 0.15, 0.03, 0.75])
			fig.colorbar(sm, cax=cbarax)
			cbarax.set_ylabel(
			'Amplitude ({})'.format({'US': 'kPa', 'elec': 'mA/m2'}[stim_type]), fontsize=fs)
			for item in cbarax.get_yticklabels():
			item.set_fontsize(fs)

			title = '{}_{}SS_{}'.format(neuron.name, 'Q' if amps is not None else '', varname)
			if amps is not None:
			title += '_vs_{}A_{}_{:.0f}%DC'.format(zscale, stim_type, DC * 1e2)
			fig.canvas.set_window_title(title)

			return fig


			def plotEqChargeVsAmp(neurons, a, Fdrive, amps=None, tstim=250e-3, toffset=50e-3, PRF=100.0,
			DCs=[1.], fs=12, xscale='lin', titrate=False, mpi=False):
			''' Plot the equilibrium membrane charge density as a function of acoustic amplitude,
			given an initial value of membrane charge density.

			:param neurons: neuron objects
			:param a: sonophore radius (m)
			:param Fdrive: US frequency (Hz)
			:param amps: US amplitudes (Pa)
			:return: figure handle
			'''

			# Determine stimulation modality
			if a is None and Fdrive is None:
			stim_type = 'elec'
			a = 32e-9
			Fdrive = 500e3
			else:
			stim_type = 'US'

			logger.info('plotting equilibrium charges for %s stimulation', stim_type)

			# Create figure
			fig, ax = plt.subplots(figsize=(6, 4))
			figname = 'charge stability vs. amplitude'
			ax.set_title(figname)
			ax.set_xlabel('Amplitude ({})'.format({'US': 'kPa', 'elec': 'mA/m2'}[stim_type]),
			fontsize=fs)
			ax.set_ylabel('$\\rm Q_m\ (nC/cm^2)$', fontsize=fs)
			if xscale == 'log':
			ax.set_xscale('log')
			for skey in ['top', 'right']:
			ax.spines[skey].set_visible(False)
			for item in ax.get_xticklabels() + ax.get_yticklabels():
			item.set_fontsize(fs)

			Qrange = (np.inf, -np.inf)

			icolor = 0
			for i, neuron in enumerate(neurons):

			nbls = NeuronalBilayerSonophore(a, neuron, Fdrive)

			# Compute reference charge variation array for zero amplitude
			_, Qref, lookups0, QSS0 = nbls.quasiSteadyStates(Fdrive, amps=0., squeeze_output=True)
			Qrange = (min(Qrange[0], Qref.min()), max(Qrange[1], Qref.max()))
			Vmeff0 = lookups0['V']
			if stim_type == 'elec': # if E-STIM case, compute steady states with constant capacitance
			Vmeff0 = Qref / neuron.Cm0 * 1e3
			QSS0 = neuron.steadyStates(Vmeff0)
			dQdt0 = -neuron.iNet(Vmeff0, np.array([QSS0[k] for k in neuron.states])) # mA/m2

			# Compute 3D QSS charge variation array
			if stim_type == 'US':
			_, _, lookups, QSS = nbls.quasiSteadyStates(Fdrive, amps=amps, DCs=DCs)
			dQdt = -neuron.iNet(lookups['V'], np.array([QSS[k] for k in neuron.states])) # mA/m2
			Afactor = 1e-3
			else:
			Afactor = 1.
			dQdt = np.empty((amps.size, Qref.size, DCs.size))
			for iA, A in enumerate(amps):
			for iDC, DC in enumerate(DCs):
			dQdt[iA, :, iDC] = dQdt0 + A * DC

			# For each duty cycle
			for iDC, DC in enumerate(DCs):
			color = 'k' if len(neurons) * len(DCs) == 1 else 'C{}'.format(icolor)

			# Initialize containers for stable and unstable fixed points
			SFPs = []
			UFPs = []

			stab_points = []

			# Generate QSS batch queue
			QSS_queue = []
			for iA, Adrive in enumerate(amps):
			lookups1D = {k: v[iA, :, iDC] for k, v in lookups.items()}
			lookups1D['Q'] = Qref
			QSS_queue.append([Fdrive, Adrive, DC, lookups1D, dQdt[iA, :, iDC]])

			# Run batch to find stable and unstable fixed points at each amplitude
			QSS_output = runBatch(nbls, 'quasiSteadyStateFixedPoints', QSS_queue, mpi=mpi)

			# Generate simulations batch queue
			sim_queue = nbls.createQueue([Fdrive], amps, [tstim], [toffset], [PRF], [DC], method)

			# Run batch to find stabilization points at each amplitude
			sim_output = runBatch(nbls, 'runIfNone', sim_queue, extra_params=[outdir], mpi=mpi)

			# Retrieve batch output
			for i, Adrive in enumerate(amps):
			SFPs += [(Adrive, Qm) for Qm in QSS_output[i][0]]
			UFPs += [(Adrive, Qm) for Qm in QSS_output[i][1]]

			# TODO: get stabilization point from simulation, if any


			# Plot charge SFPs and UFPs for each acoustic amplitude
			lbl = '{} neuron - {{}}stable fixed points @ {:.0f} % DC'.format(
			neuron.name, DC * 1e2)
			if len(SFPs) > 0:
			A_SFPs, Q_SFPs = np.array(SFPs).T
			ax.plot(np.array(A_SFPs) * Afactor, np.array(Q_SFPs) * 1e5, 'o', c=color,
			markersize=3, label=lbl.format(''))
			if len(UFPs) > 0:
			A_UFPs, Q_UFPs = np.array(UFPs).T
			ax.plot(np.array(A_UFPs) * Afactor, np.array(Q_UFPs) * 1e5, 'x', c=color,
			markersize=3, label=lbl.format('un'))

			# If specified, compute and plot the threshold excitation amplitude
			if titrate:
			if stim_type == 'US':
			Athr = nbls.titrate(Fdrive, tstim, toffset, PRF=PRF, DC=DC)
			ax.axvline(Athr * Afactor, c=color, linestyle='--')
			else:
			for Arange, ls in zip([(0., amps.max(amps.min(), 0.)), ()], ['--', '-.']):
			Athr = neuron.titrate(tstim, toffset, PRF=PRF, DC=DC, Arange=Arange)
			ax.axvline(Athr * Afactor, c=color, linestyle=ls)
			icolor += 1

			# Post-process figure
			ax.set_ylim(np.array([Qrange[0], 0]) * 1e5)
			ax.legend(frameon=False, fontsize=fs)
			fig.tight_layout()

			fig.canvas.set_window_title('QSS_Qstab_vs_{}A_{}_{}_{}%DC{}'.format(
			xscale,
			'_'.join([n.name for n in neurons]),
			stim_type,
			'_'.join(['{:.0f}'.format(DC * 1e2) for DC in DCs]),
			'_with_thresholds' if titrate else ''
			))

			return fig

R4670/PySONIC/plt/QSS.py
6cb5dac16172master

QSS.py
View Options

R4670/PySONIC/plt/QSS.py6cb5dac16172master

QSS.pyView Options

R4670/PySONIC/plt/QSS.py
6cb5dac16172master

QSS.py
View Options