Model of peripheral nerve with ephaptic coupling (Capllonch-Juan & Sepulveda 2020)

 Download zip file 
Help downloading and running models
Accession:263988
We built a computational model of a peripheral nerve trunk in which the interstitial space between the fibers and the tissues is modelled using a resistor network, thus enabling distance-dependent ephaptic coupling between myelinated axons and between fascicles as well. We used the model to simulate a) the stimulation of a nerve trunk model with a cuff electrode, and b) the propagation of action potentials along the axons. Results were used to investigate the effect of ephaptic interactions on recruitment and selectivity stemming from artificial (i.e., neural implant) stimulation and on the relative timing between action potentials during propagation.
Reference:
1 . Capllonch-Juan M, Sepulveda F (2020) Modelling the effects of ephaptic coupling on selectivity and response patterns during artificial stimulation of peripheral nerves. PLoS Comput Biol 16:e1007826 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Extracellular; Axon;
Brain Region(s)/Organism:
Cell Type(s): Myelinated neuron;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; Python;
Model Concept(s): Ephaptic coupling; Stimulus selectivity;
Implementer(s):
/
publication_data
dataset_01__fields
code
data
models
settings
src
x86_64
gaines_sensory_flut.mod *
gaines_sensory_mysa.mod *
gaines_sensory_node.mod *
gaines_sensory_stin.mod *
MRG_AXNODE.mod *
algebra.py *
analysis.py *
anatomy.py *
biophysics.py *
circlepacker.py *
contourhandler.py *
convert_results_to_xyz.py
fill_nerve.py *
geometry.py *
get_extstim.py *
PaperFig_0.8_textwidth.mplstyle *
read_results.py *
show_results.py
sim_launcher.py *
simcontrol.py
tessellations.py *
tools.py *
visualisation.py *
workspace.py *
                            
import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
from datetime import datetime
from matplotlib.colors import LogNorm
from matplotlib.ticker import LogFormatter
from matplotlib import ticker, cm

import workspace as ws
import read_results



def show_slice(
		ax, 
		filename, 
		axis, 
		where, 
		method='nearest', 
		rsx=100, 
		rsy=100, 
		rsz=100, 
		window=(0., 1.), 
		nlevels=100, 
		vmin=None, 
		vmax=None, 
		extend='neither', 
		logscale=False, 
		save=False, 
		zorder=0, 
		**kwargs
	):
	""" Opens a file 'filename' and shows the data at a slice perpendicular
	to the 'axis' axis at a relative position 'where', where 'where' ranges
	from 0. to 1. """

	# Keyword arguments
	if len(kwargs) >= 2:
		xyzunits = kwargs['xyzunits']
		vunits = kwargs['vunits']
	else:
		xyzunits = ''
		vunits = ''
	try:
		cmap = kwargs['cmap']
	except KeyError:
		cmap = plt.cm.jet
	try:
		show_data_points = kwargs['show_data_points']
	except KeyError:
		show_data_points = False
	try:
		show_intp_points = kwargs['show_intp_points']
	except KeyError:
		show_intp_points = False

	# Set colors for bad values
	if False:
		cmap.set_under(color='black')
		cmap.set_over(color='black')
	if False:
		cmap.set_bad(color='black')

	# Load data
	data = np.loadtxt(filename, delimiter=',')
	x = data[:,0]
	y = data[:,1]
	z = data[:,2]
	v = data[:,3]

	# Sides
	xspan = x.max() - x.min()
	yspan = y.max() - y.min()
	zspan = z.max() - z.min()

	# Window
	wl, wr = window
	# Crop the window limits just in case
	wl = max(wl, 0.)
	wl = min(wl, 1.)
	wr = max(wr, 0.)
	wr = min(wr, 1.)

	# Regular mesh
	xi_ = np.linspace(x.min(), x.max(), rsx)
	yi_ = np.linspace(y.min(), y.max(), rsy)
	zi_ = np.linspace(z.min(), z.max(), rsz)
	
	if axis == 'x':
		# Indices for the windows
		wl_i = int(rsx * wl)
		wr_i = int(rsx * wr)
		# Create the mask for the data
		xleft = x.min() + xspan * wl
		xrght = x.min() + xspan * wr
		mask = np.ma.masked_where((x >= xleft) & (x <= xrght), x).mask
		# Crop the corresponding array for the regular mesh
		xi_ = xi_[wl_i:wr_i]
	elif axis == 'y':
		# Indices for the windows
		wl_i = int(rsy * wl)
		wr_i = int(rsy * wr)
		# Create the mask for the data
		yleft = y.min() + yspan * wl
		yrght = y.min() + yspan * wr
		mask = np.ma.masked_where((y >= yleft) & (y <= yrght), y).mask
		# Crop the corresponding array for the regular mesh
		yi_ = yi_[wl_i:wr_i]
	elif axis == 'z':
		# Indices for the windows
		wl_i = int(rsz * wl)
		wr_i = int(rsz * wr)
		# Create the mask for the data
		zleft = z.min() + zspan * wl
		zrght = z.min() + zspan * wr
		mask = np.ma.masked_where((z >= zleft) & (z <= zrght), z).mask
		# Crop the corresponding array for the regular mesh
		zi_ = zi_[wl_i:wr_i]

	# 'mesh' the regular mesh
	xi, yi, zi = np.meshgrid(xi_, yi_, zi_)

	# Mask data
	x = x[mask]
	y = y[mask]
	z = z[mask]
	v = v[mask]

	# Interpolation
	t0 = datetime.now()
	vi = griddata((x,y,z), v, (xi,yi,zi), method=method)
	t1 = datetime.now()
	print('Time needed: %f s'%((t1 - t0).total_seconds()))

	# Figure

	# Levels for colouring
	vi_mkd_cpd = np.ma.masked_where(np.isnan(vi), vi).compressed()
	vmin_, vmax_ = vi_mkd_cpd.min(), vi_mkd_cpd.max()
	if vmax is None:
		vmax = vmax_
	if vmin is None:
		vmin = vmin_
	# levels_i = np.linspace(vmin, vmax, nlevels)
	# Avoid an error due to a flat level list
	if vmin == vmax:
		vmin -= 0.1
		vmax += 0.1
	# Re-build the levels
	print('Limits:', vmin, vmax)
	levels_i = np.linspace(vmin, vmax, nlevels)
	loglevels = np.logspace(
				np.log10(min(np.abs(vmin), np.abs(vmax))), 
				np.log10(max(np.abs(vmin), np.abs(vmax))), 
				nlevels
			)
	# Minimum and maximum exponents
	min_exp = int(np.log10(min(np.abs(vmin), np.abs(vmax))))
	max_exp = int(np.log10(max(np.abs(vmin), np.abs(vmax)))) + 1
	# Ticks for the colorbar in case it's logscale
	cbticks = 10. ** np.arange(min_exp, max_exp + 1, 1)
	cbticks_ = cbticks.tolist()
	for i in range(2, 10, 1):
		cbticks_ = cbticks_ + (float(i) * cbticks).tolist()
	cbticks_.sort()
	# print(cbticks_, min_exp, max_exp)
	cbticklabels = []
	for c in cbticks_:
		if c in (20, 50, 100, 200, 500, 1000):
			cbticklabels.append('%i'%c)
		else:
			cbticklabels.append('')

	if axis == 'x':
		s = int(where * rsx) - wl_i
		if not logscale:
			cf = ax.contourf(
					yi[:, s, :], 
					zi[:, s, :], 
					vi[:, s, :], 
					levels_i, 
					cmap=cmap, 
					extend=extend, 
					zorder=zorder
				)
		else:
			cf = ax.contourf(
					yi[:, s, :], 
					zi[:, s, :], 
					np.abs(vi[:, s, :]), 
					loglevels, 
					cmap=cmap, 
					norm=LogNorm(
						vmin=min(np.abs(vmin), np.abs(vmax)), 
						vmax=max(np.abs(vmin), np.abs(vmax))
						), 
					# locator=ticker.LogLocator(), 
					zorder=zorder
				)
		if show_data_points:
			ax.plot(y, z, 'k.', markersize=0.6)
		if show_intp_points:
			ax.scatter(yi, zi, c='k', s=1)
		ax.set_title('x = %f'%xi[s, s, s] + ' ' + xyzunits)
		ax.set_xlabel('y (' + xyzunits + ')')
		ax.set_ylabel('z (' + xyzunits + ')')
	elif axis == 'y':
		s = int(where * rsy) - wl_i
		if not logscale:
			cf = ax.contourf(
					xi[s, :, :], 
					zi[s, :, :], 
					vi[s, :, :], 
					levels_i, 
					cmap=cmap, 
					extend=extend, 
					zorder=zorder
				)
		else:
			cf = ax.contourf(
					xi[s, :, :], 
					zi[s, :, :], 
					np.abs(vi[s, :, :]), 
					loglevels, 
					cmap=cmap, 
					norm=LogNorm(
						vmin=min(np.abs(vmin), np.abs(vmax)), 
						vmax=max(np.abs(vmin), np.abs(vmax))
						), 
					# locator=ticker.LogLocator(), 
					zorder=zorder
				)
		if show_data_points:
			ax.plot(x, z, 'k.', markersize=0.6)
		if show_intp_points:
			ax.scatter(xi, zi, c='k', s=1)
		ax.set_title('y = %f'%yi[s, s, s] + ' ' + xyzunits)
		ax.set_xlabel('x (' + xyzunits + ')')
		ax.set_ylabel('z (' + xyzunits + ')')
	elif axis == 'z':
		s = int(where * rsz) - wl_i
		if not logscale:
			cf = ax.contourf(
					xi[:, :, s], 
					yi[:, :, s], 
					vi[:, :, s], 
					levels_i, 
					cmap=cmap, 
					extend=extend, 
					zorder=zorder
				)
		else:
			cf = ax.contourf(
					xi[:, :, s], 
					yi[:, :, s], 
					np.abs(vi[:, :, s]), 
					loglevels, 
					cmap=cmap, 
					norm=LogNorm(
						vmin=min(np.abs(vmin), np.abs(vmax)), 
						vmax=max(np.abs(vmin), np.abs(vmax))
						), 
					# locator=ticker.LogLocator(), 
					zorder=zorder
				)
		if show_data_points:
			ax.plot(x, y, 'k.', markersize=0.6)
		if show_intp_points:
			ax.scatter(xi, yi, c='k', s=1)
		ax.set_title('z = %f'%zi[s, s, s] + ' ' + xyzunits)
		ax.set_xlabel('x (' + xyzunits + ')')
		ax.set_ylabel('y (' + xyzunits + ')')
		ax.set_aspect('equal')

	if logscale:
		l_f = LogFormatter(10, labelOnlyBase=False)
		# cb = plt.colorbar(cf, ticks=cbticks_, format=l_f)
		cb = plt.colorbar(cf, format=l_f)
		# cb = plt.colorbar(cf, ticks=cbticks)
		# cb = plt.colorbar(cf)
		cb.set_ticks(cbticks_)
		# cb.set_ticklabels(cbticklabels)
		# cb.set_norm(LogNorm())
		# print(cbticklabels)
	else:
		cb = plt.colorbar(cf, ax=ax)

	cb.set_label(vunits)

def z_slice(ax, filename, folder, z, t, tarray):
	""" This slicing function uses a much more efficient 
	method to interpolate along the z-axiz """

	# This is from a xyz file

	# Load data
	data = np.loadtxt(filename, delimiter=',')
	x = data[:,0]
	y = data[:,1]
	z_ = data[:,2]
	v = data[:,3]

	# Set the cables
	xy_all = np.array(list(zip(x, y)))
	xy_unique = np.unique(xy_all, axis=0)
	nc = len(xy_unique)

	data_intp = np.zeros(nc)

	# Construct zprofile
	for i, xy in enumerate(xy_unique):

		# Select the corresponding rows
		rows = np.where((xy == xy_all).all(axis=1))[0]

		# Sort z_ just in case
		# argsort = np.argsort(z_[rows])
		# z_.sort()
		# Sort v accordingly
		# Not yet
		# v.sort(argsort)

		# Select z_
		z__ = z_[rows]
		z_lower = z__[np.where(z__ < z)]
		z__minus_z = np.abs(z__ - z)
		z_lower_minus_z = np.abs(z_lower - z).min()
		here = np.where(z__minus_z == z_lower_minus_z)[0][0]


		# Interpolate
		v__ = v[rows]
		data_intp[i] = (z__minus_z[here + 1] * v__[here] + z__minus_z[here] * v__[here + 1]) / (z__[here + 1] - z__[here])
		# if v[here] < -100.:

	# Show result
	ax.scatter(xy_unique[:, 0], xy_unique[:, 1], c=data_intp, s=10, cmap=plt.cm.jet_r)
	ax.set_aspect('equal')

	# KEEP THE CODE BELOW, DON'T DELETE IT.
	# It works only when there are results
	if False:
		# Read results
		data, zprofile, xyr = read_results.read_results(folder)
		xx_, yy_, rr_ = xyr

		nc = len(data)
		data_intp = np.zeros(nc)

		for i in range(nc):

			# Choose data
			try:
				data_ = data[i]['vext[0]']
			except KeyError:
				data_ = data[i]['vext[1]']

			# Find the zprofile[i] value lying just below z
			zpr = np.array(zprofile[i].values())
			zpr_minus_z = np.abs(zpr - z)
			here = np.where(zpr_minus_z == zpr_minus_z.min())

			# Interpolate v
			it = np.where(tarray == t)[0][0]
			v = data_[:, it]
			# data_intp[i] = (zpr_minus_z[here] * v[here] + zpr_minus_z[here + 1] * v[here + 1]) / (zpr[here + 1] - zpr[here])
			data_intp[i] = (zpr_minus_z[here + 1] * v[here] + zpr_minus_z[here] * v[here + 1]) / (zpr[here + 1] - zpr[here])

		# Show result
		ax.scatter(xx_, yy_, c=data_intp)

	return xy_unique, data_intp

def axons_activity(time, data, tlims):
	""" Show axons' activity """		
	labels = {'Sti': 'Stimulation point', 'Rec': 'Recording point'}
	colors = {'Sti': 'k', 'Rec': 'b'}
	fig, ax = plt.subplots()
	# Potential
	for vecs in data:
		for k, vv in vecs.items():
			ax.plot(time, vv, colors[k])
	ax.set_xlim(tlims[0], tlims[1])
	ax.set_xlabel('t (ms)')
	ax.set_ylabel('Vm (mV)')
	ax.set_title('Axons Vm (mV)')
	plt.savefig(os.path.join(ws.folders["data/results"], 'images/Axons.png'))

def cross_sec_results(data):
	""" Cross-section of the nerve and results """
	# Map of axons
	act_axs = []
	k = 0
	# Count it as having an AP if Vm is above 20 mV
	for i, v in enumerate(data):
		# Don't use the stimulation point to check if there is an AP, because
		# the stimulation artefact can give a false positive.
		# v = v['Sti']
		# Use instead the recording position. This is in better accordance to 
		# Raspopovic & al. (2011), who say that the AP needs to travel the whole
		# length for the axon to be tagged as activated
		v = v['Rec']
		if v.max() > -30.:
			act_axs.append(True)
		else:
			act_axs.append(False)

	# Figure
	fig, ax = plt.subplots()

	zorder = 0

	# Triangulation
	ws.nvt.draw_contours(ax, c='k', zorder=zorder)
	ws.nvt.draw_axons_and_points(ax, c='k', act_axons=act_axs, zorder=zorder)
	ws.nvt.draw_triangulation(ax, c='darkgrey', lw=0.4, zorder=zorder)
	zorder += 1

	# Scatter all points in black
	x, y = ws.nvt.x, ws.nvt.y
	ax.scatter(x, y, c='k', s=1, zorder=zorder)

df2

	# Arrange figure
	xleft = x.min()
	xrght = x.max()
	ybot = y.min()
	ytop = y.max()
	margin = 100

	ax.set_xlim(xleft - margin, xrght + margin)
	ax.set_ylim(ybot - margin, ytop + margin)
	ax.set_aspect('equal')
	ax.set_xlabel('x')
	ax.set_ylabel('y')
	title_ups = 'Tessellated nerve and current injections'
	title_len = 'Length (z-axis): %0.00f um'%ws.length
	title_inj = 'Injected current(s): %0.00f uA'%(ws.icamp * 1.e-3)
	title = title_ups + '\n' + title_len + '\n' + title_inj
	ax.set_title(title)

	# Field
	ws.nvt.pd.draw(ax, colour='dimgrey', alpha=1., linestyle='-', linewidth=.5,
			markers=None, title=title, zorder=zorder)
	fig.tight_layout()
	fig.savefig(os.path.join(ws.folders["data/results"], 'images/Nerve_and_tessellation.png'))

Loading data, please wait...