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

 Download zip file 
Help downloading and running models
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.
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;
Gap Junctions:
Simulation Environment: NEURON; Python;
Model Concept(s): Ephaptic coupling; Stimulus selectivity;
gaines_sensory_flut.mod *
gaines_sensory_mysa.mod *
gaines_sensory_node.mod *
gaines_sensory_stin.mod *
aaa_info_dataset * * * * * * * * * * * * * * * * * *
Miguel Capllonch Juan
This module performs analysis of the results

import numpy as np
import os
import csv
from neuron import h

import workspace as ws

def remove_previous_results():
	""" Remove all previous results, either text or figures """
	resultsdir = ws.folders["data/results"]
	# Remove all .csv in the results
	for item in os.listdir(resultsdir):
		if item.endswith(".csv"):
			os.remove(os.path.join(resultsdir, item))
	# Remove all images
	imgdir = os.path.join(resultsdir, 'images')
	for item in os.listdir(imgdir):
		os.remove(os.path.join(imgdir, item))

def set_time_recording():
	""" Set the recording for the time vector of the simulation """

	# Record time
	ws.tvec = h.Vector()

def set_recording(seg, var):
	""" Set up a recording vector in hoc for any given segment 
	in the system """
	recording = h.Vector()
	# recording.buffer_size(ws.nt)
	ws.counters['recordings'] += 1

def record_cables():
	""" Set up the recording vectors for the desired cables and the
	desired variables """
	dr = ws.settings["data records"]
	rc = dr["record cables"]
	ra = dr["record axons"]
	if rc == "all":
		# Record all types of cables
		for cable in ws.cables:
			if cable.cabletype == "NAELC":
			elif cable.cabletype == "Axon":
				if ra == "all":
					# Record all axons
	elif rc == "NAELC":
		# Record only NAELC
		for cable in ws.cables:
			if cable.cabletype == "NAELC":
	elif rc == "axons":
		# Record only axons
		if ra == "all":
			# Record all axons
			for cable in ws.cables:
				if cable.cabletype == "Axon":
	elif rc == "None":
		# Nothing to record

def remove_artefacts(time, recs):
	""" Remove the stimulation artefacts in the electrode recordings """

	# Start the analysis

	# Turn time and recordings into arrays to be able to work with them
	time = np.array(time)
	recs = np.array(recs)
	# Make everything positive
	rpos = np.abs(recs)
	# Average
	ravg = rpos.mean()
	# Identify where rpos is higher than its mean: there should be the artefacts
	rhigh = np.where(rpos > ravg)

	# Identify the separation between pulses
	dt = np.gradient(time).mean()
	# The 1.0000001 factor is just to make sure we get > and not >=
		cuts = np.where(np.gradient(time[rhigh]) > 1.000001 * dt)[0][0::2]
	except ValueError:
		# This may happen, for instance, if np.gradient can't be 
		# performed. In case of error, just say there's no cuts
		cuts = []
	# Split the pulses
	segments = []
	for cut in cuts:
		ttrr = []
		ttrr.append(time[rhigh][:cut + 1])
		ttrr.append(recs[rhigh][:cut + 1])

	# Last pulse
	ttrr = []
		ttrr.append(time[rhigh][cut + 1:])
		ttrr.append(recs[rhigh][cut + 1:])
	except UnboundLocalError:
		# There are no cuts

	# Identify the jumps at the cuts

	# Choose which data are outisde the pulses. Use rhigh as a mask
	toutsp = > ravg, time)
	routsp = > ravg, recs)
	# Get the jumps from that
	jumps = np.where(routsp.mask[1:] != routsp.mask[:-1])[0]
	# Pulse counter
	k = 0
	# (Average) jump heights of each pulse
	pulseheights = []
	# New values for the recordings during the pulses
	rnews = []
	# Iterate over jumps to get the avg. jump height for each pulse
	for i, jump in enumerate(jumps):
		jheight = abs(recs[jump + 1] - recs[jump])
		if i % 2 == 0:
			jumpup = jheight
		elif len(segments) > 0:
			rawrec = segments[k][1]
			avgjh = 0.5 * (jumpup + jheight) * np.sign(rawrec)
			# Subtract this value from the corresponding segment in recs
			rnew = rawrec - avgjh
			# Add one to the pulse counter
			k += 1
			rnews = recs.copy()

	# Finally, get the entire thing
	finalrecs = recs.copy()
	if len(segments) > 0:
		finalrecs[rhigh] = [item for sublist in rnews for item in sublist]
	return finalrecs

def save_vmem(data):
	""" Save the membrane potentials """
	f_sti = open(os.path.join(ws.folders["data/results"], 'vmemb_stm_site.csv'), 'w')
	f_rec = open(os.path.join(ws.folders["data/results"], 'vmemb_rec_site.csv'), 'w')
	writers = {
		'Sti': csv.writer(f_sti),
		'Rec': csv.writer(f_rec)
	for i, vecs in enumerate(data):
		for k, vv in vecs.items():
			writers[k].writerow([i] + vv.as_numpy().tolist())

def save_electrode_recordings():
	""" Save the recordings from all electrodes into files """

	# Time
	tvarr = ws.tvec.as_numpy()

	# Iterate over electrodes:
	for name, electrode in ws.electrodes.items():
		if electrode.role == 'recording':
			# Iterate over rings on the electrode
			for ring in electrode.rings:
				ring_number = ring.ring_number

				# Iterate over pads on the ring
				for pad, rec in ring.recordings.items():
					# Recording ID
					rec_id = rec['recording id']
					# Fetch the vector with the data from ws
					v = ws.recordings[rec_id]

					# Save with artefacts
					with open(os.path.join(ws.folders["data/results"], 
						%(name.replace(' ', '_'), ring_number, pad)), 'w') as f:
						fw = csv.writer(f)
						fw.writerow(['Time (ms)', 'V (mV)'])
						for tt, vv in zip(tvarr, v):
							fw.writerow([tt, vv])

					# Remove the artefacts
						v = remove_artefacts(tvarr, v.as_numpy())
					except Exception as e:
						message = "ERROR in %s.save_electrode_recordings:\n"%__name__ + type(e).__name__ + ':\n' + str(e)
						# Don't terminate, just go on
						# ws.terminate()

					# Save
					with open(os.path.join(ws.folders["data/results"], 
						%(name.replace(' ', '_'), ring_number, pad)), 'w') as f:
						fw = csv.writer(f)
						fw.writerow(['Time (ms)', 'V (mV)'])
						for tt, vv in zip(tvarr, v):
							fw.writerow([tt, vv])

def save_cable_recordings():
	""" Save the recordings from cables into files """

	# Time
	tvarr = ws.tvec.as_numpy()

	# Iterate over cables
	for cable in ws.cables:

		if cable.has_recordings:

			# Open file
			with open(os.path.join(ws.folders["data/results"], '%s%04i.csv'\
				%(cable.cabletype, cable.specific_number)), 'w') as f:
				fw = csv.writer(f)

				# Write the most basic properties of this cable:
				# Position and radius
				fw.writerow([cable.x, cable.y, cable.r])

				# Write the anatomy of the axon along the z-axis
				lp = ["segment positions:"]
				for i, seg in enumerate(ws.segments[cable.cable_number]):

				# Iterate over recorded variables
				for var, recs in cable.recordings.items():

					# Write the name of the variable
					# Iterate over recording vectors (or segments)
					for rec in recs:

						# Get numerical data
						data = ws.recordings[rec["number"]].as_numpy()

						# Store results
						# fw.writerow([ws.sectypes[cable.cable_number][i]] + [x for x in data])
						fw.writerow([rec["segment"]] + [x for x in data])

def save_memcurrents(i_mem_, i_mem_ly2):
	""" Save the membrane currents into files """
	x, y = ws.nvt.x, ws.nvt.y
	with open(os.path.join(ws.folders["data/results"], 'membranecurrents.csv'), 'w') as f:
		fw = csv.writer(f)
		fw.writerow(['Fiber', 'secname', 'x', 'y', 'z', 'i_membrane time series', 'i_membrane on layer 2'])
		for i, fiber in i_mem_.items():
			for j, segment in enumerate(fiber):
				fw.writerow([i, ws.segments[i][j], x[i], 
					y[i], ws.zprofs[i][j]] + segment.tolist() + \

def store_i_mem(i_t):
	""" Store the membrane currents of all the cells.
	So far, this only works for axons with 2 extracellular layers """
	for i in np.arange(ws.nNAELC,, 1):
		for j, seg in enumerate(ws.segments[i]):
			if 'node' in
				ws.i_mem_[i][j, i_t] = seg.i_membrane * seg.area()
				# Outwards is positive
				ws.i_mem_ly2[i][j, i_t] = (seg.vext[0] - seg.vext[1]) * seg.xg[0] * seg.area()
				ws.i_mem_[i][j, i_t] = 0.
				ws.i_mem_[i][j, i_t] = seg.i_membrane * seg.area()
				# Outwards is positive
				ws.i_mem_ly2[i][j, i_t] = (seg.vext[0] - seg.vext[1]) * seg.xg[0] * seg.area()

def postprocessing():
	""" Do all the post-processing tasks after a simulation """

	# Save the recordings from cables into files

	# Save the electrode recordings into files
	# if that's what we want
	if ws.settings["electrode recordings"]:
		# and if the nerve model we're using is compatible with that
		if ws.settings["nerve model"] != "simple":

Loading data, please wait...