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;
AXNODE.mod *
aaa_info_dataset * * * * * * * * * * * * * * * * *
This module handles virtually all the necessary processes for setting 
up stimulation and recording from electrodes.
One more basic function (analysis.set_recording) that is used in this 
module is defined in the module analysis. This is so because this 
function is not necessarily used by an electrode

import os
import numpy as np
from collections import OrderedDict
from neuron import h

import anatomy
import geometry as geo
import analysis as als
import workspace as ws

class CuffRing():
	Ring of a cuff electrode

	Input parameters:

	z: Position of the ring along the z-axis
	d: Diameter
	npads: Number of pads
	thts: Angular positions of the pads

	Note about the stimulation set-up: this class has no method 
	set_stimulation because some rings won't have any stimulation at 
	all if it is not defined in the configuration files
	def __init__(self, ring_number, z, center, d, npads, thts):

		self.ring_number = ring_number
		self.z = z = center
		self.d = d
		self.npads = npads
		self.thts = thts
		# Dictionary with the information about all the injections
		self.injections = {}
		# Dictionary with the information about all the recordings
		self.recordings = {}

		# Locate the contact points with the nerve

	def locate_pads_on_nerve(self):
		""" Find the locations in the nerve corresponding 
		to the pads """

		# Iterate over the ring's pads
		self.contact_points = []

		for i_pad in range(self.npads):
			# The pad's angular position determines which cable cell(s) 
			# is (are) stimulated

			# Select which of the targets is in contact 
			# with a pad 'i_pad' according to its angular position
			angdiff = np.abs(ws.contour_angles - self.thts[i_pad])
			which_cables = np.where(angdiff == angdiff.min())[0]
			# Iterate over possible cables
			ncables_here = len(which_cables)
			for i_cable in which_cables:
				# For cable_cell i, see where this is
					i_sec, zpos = anatomy.locate(ws.seclens[i_cable], 
				except TypeError:
					msg = 'ERROR: A ring of an electrode could not be placed on the nerve'
				# Add the point
				point = {
						'cable': i_cable, 
						'section': i_sec, 
						'z': zpos
				except IndexError:

	def assign_injections(self, pad, protocol):
		""" Assign current injections to the pads """
		# Dictionary with the information about all the injections
		self.injections[pad] = []

		# Iterate over the contact points of this pad
		# NO! Just one pad is allowed. Otherwise, the current with 
		# intensity 'amp', which is meant for only one injection point, 
		# would be added to all the self.contact_points, so we would 
		# have (len(self.contact_points) - 1) more current injections 
		# than we intended
		point = self.contact_points[pad][0]

		# IClamp object(s)
		# Only square pulses are supported for now
		if protocol['type'] == 'square pulses':
			square_pulses(self.injections[pad], point, protocol)
		if protocol['type'] == 'ground':

	def set_recordings(self):
		""" Assign the recording vectors to each pad """

		# Iterate over all pads on this ring
		for pad in range(self.npads):

			# Dictionary with the information about all the recordings
			self.recordings[pad] = {}

			# Only one is allowed for each pad in this case
			point = self.contact_points[pad][0]
			# Set the recording for this pad
			set_recording(self.recordings[pad], point, 'vext[0]')

class CuffElectrode():
	Cuff electrode with as many rings as one wants

	Input parameters:

	z: Position (center) along the z-axis
	nrings: Number of rings
	rng_sp: Separation between rings
	nppr: Number of pads per ring
	l: Length of the cuff
	d_in: Inner diameter
	thickness: Insulator thickness
	rho: Insulator resistivity (Ohm*cm)
	def __init__(self, name): = name
		settings = ws.electrodes_settings[name]
		self.type = settings['type']
		self.role = settings['role']
		self.z = settings['z-position']
		self.nrings = settings['number of rings']
		self.rng_sp = settings['inter-ring distance']
		self.nppr = settings['pads per ring']
		self.length = settings['length']
		self.thickness = settings['thickness']
		self.rho  = settings['insulator resistivity'] = ws.nvt.circumcenter
		self.d_in = ws.nvt.circumdiameter
		self.d_out = self.d_in + 2. * self.thickness
		if self.role == 'stimulation':
			self.stimulation_protocol = \
				settings['stimulation protocol']

		# Necessary by-products

		# Ends of the cuff
		z = self.z
		l = self.length
		nrings = self.nrings
		left_end = z - 0.5 * l
		rght_end = z + 0.5 * l

		# Positions of the rings
		rings_halfspan = 0.5 * (nrings - 1) * self.rng_sp
		z_ring_left = z - rings_halfspan
		z_ring_rght = z + rings_halfspan
		self.rings_poss = np.linspace(z_ring_left, z_ring_rght, nrings)

		# Angular positions of the pads on the rings
		self.thts = settings['angular positions']
		# Consistency check
		if len(self.thts) != self.nppr:
			msg = 'ERROR: Different number of pads and angles for them'

		# Warning if size mismatch
		if z_ring_left < left_end or z_ring_rght > rght_end:
			msg = 'ERROR: Cuff electrode rings outside the cuff.'

		# Store variables into attributes
		self.left_end = left_end
		self.rght_end = rght_end
		self.z_ring_left = z_ring_left
		self.z_ring_rght = z_ring_rght

		# Create rings

	def create_rings(self):
		""" Create the rings """
		self.rings = []
		for i in range(self.nrings):
			z = self.rings_poss[i]
			self.rings.append(CuffRing(i, z,, self.d_in, 
				self.nppr, self.thts))

	def set_stimulation(self):
		""" Set up the current injections according to the assigned 
		protocol """
		sp = self.stimulation_protocol
		for pad, protocol in sp.items():
			ring_number, pad_no = [int(x) for x in pad.split(',')]
			self.rings[ring_number].assign_injections(pad_no, protocol)

	def set_recordings(self):
		""" Set up the recording vectors on the pads """
		for ring in self.rings:

class PointElectrodeSet():
	Point electrode set for intracellular current injections, extracellular current sources, recordings or connections to ground
	def __init__(self, name, index): = name
		self.index = index
		self.settings = ws.electrodes_settings[name]
		self.type = self.settings['type']
		self.role = self.settings['role']
	def set_stim_info(self):
		""" Set up just a basic thing about stimulation info """

		# Stimulation protocol
		if self.role == 'stimulation':
				self.stimulation_protocol = \
					self.settings['stimulation protocol'][self.index]
			except KeyError:
				self.stimulation_protocol = \
					self.settings['stimulation protocol']
		if self.role == 'ground':
			self.stimulation_protocol = {
				'type': 'ground'

	def set_stimulation(self):
		""" Set up the current injections according to the assigned protocol """
		# Set up the injection if we are using the RN for all types of electrodes, or if we are using a "simple" model, only for intracellular injections
		if (ws.settings["nerve model"] == "resistor network") | (ws.settings["nerve model"] == "simple" and self.type == "intracellular"):
			# Dictionary containing the information about all 
			# the current injections
			self.injections = {}

			# Injections list
			# There's only one 'pad' for an intracellular electrode
			self.injections[0] = []

			# Fetch protocol
			protocol = self.stimulation_protocol

			# Iterate over all the points in this set
			for point in self.points:
				# Create IClamp object(s)
				# Only square pulses are supported for now
				if protocol['type'] == 'square pulses':
					ws.log('Setting up injection on %s'%str(point))
					square_pulses(self.injections[0], point, protocol)
				# Connections to ground
				if protocol['type'] == 'ground':
					ws.log('Setting up ground on %s'%str(point))

	def set_recordings(self):
		""" Set up the recording vectors in the points of this electrode set """
		self.recordings = {}
		self.recordings[0] = {}
		for point in self.points:
			set_recording(self.recordings[0], point, 'v')

class IntracellularElectrodeSet(PointElectrodeSet):
	Intracellular electrode for intracellular current injections or recordings
	def __init__(self, name, index=0):

		PointElectrodeSet.__init__(self, name, 0)
		settings = self.settings

			self.z = settings['z-position']
		except KeyError:
			self.axons = settings['axons']
				self.node = settings['node of Ranvier']
			except KeyError:
				self.node = settings['internode']
				self.position = settings['position']
				# For a node of Ranvier, inject current in the middle
				self.position = 0.5
			# Warning if position mismatch
			if self.z < 0. or self.z > ws.length:
				ws.log('ERROR: Electrode is outside the limits of the nerve.')
		# Note: the warnings or errors for targetting a non-existing 
		# node will be automatically raised by NEURON

		# Set points
		self.points = []
		# Setting up these points is necessary to set up stimulation and recording
		# Besides, it gives the process uniformity across different types of electrodes
		if self.axons[0] == "all":
			self.axons = range(ws.naxons_total)

		for axon in self.axons:
			axon = int(axon)
			if ws.settings["nerve model"] == "resistor network":
				cable_i = ws.nNAELC + axon
			elif ws.settings["nerve model"] == "simple":
				cable_i = axon
			k = 0
			if "node of Ranvier" in settings.keys():
				for j, sec in enumerate(ws.sections[cable_i]):
					secname =
					if "NODE" in secname or "node" in secname:
						if k == self.node:
						k += 1

						'cable': cable_i, 
						'section': j, 
						'z': self.position


class ExtracellularPointElectrode(PointElectrodeSet):
	Extracellular point electrode for current sources, recordings or connections to ground
	def __init__(self, name, index=0):

		PointElectrodeSet.__init__(self, name, index)
		settings = self.settings

		try: = settings['point']
		except KeyError: = settings['points'][self.index]
		self.x, self.y, self.z =
		# Find the cable this corresponds to
		dists = geo.dist((self.x, self.y), (ws.nvt.x, ws.nvt.y))
		these = np.where(dists == dists.min())[0]
		# print('Points:', these)
		# Choose only one cable
		cable_i = these[0]

		# print('Cable:', ws.nvt.cables[cable_i])

		# Locate the section and segment (z value) of the cable
		j_sec, z_on_j = anatomy.locate(ws.seclens[cable_i], self.z)

		# print('Section and segment:', j_sec, z_on_j)

		# Cable and point of the cable to work on
		self.points = [
						'cable': cable_i, 
						'section': j_sec, 
						'z': z_on_j


class ExtracellularPointElectrodeSet():
	Set of ExtracellularPointElectrode instances
	def __init__(self, name): = name
		self.role = ws.electrodes_settings[name]["role"]
		self.type = "extracellular points set"
		self.electrodes = []

	def set_stimulation(self):
		""" Create the list of extracellular point electrodes """
		for i, (x, y, z) in ws.electrodes_settings[]["points"].items():
			# Create the ExtracellularPointElectrode instance and store 
			# it
			pe = ExtracellularPointElectrode(, i)

class GroundPointSet():
	Set of points which are connected to ground
	def __init__(self, name): = name
		self.role = "ground"
		self.type = "ground points set"

	def set_stimulation(self):
		""" Connect grounds """

		if ws.settings["nerve model"] == "resistor network":

			# Go straigth to connect the wanted points to ground
			for (x, y, z) in ws.electrodes_settings[]["points"]:

				# Find the cable this corresponds to
				dists = geo.dist((x, y), (ws.nvt.x, ws.nvt.y))
				these = np.where(dists == dists.min())[0]
				# print('Points:', these)
				# Choose only one cable
				cable_i = these[0]

				# print('Cable:', ws.nvt.cables[cable_i])

				# Locate the section and segment (z value) of the cable
				j_sec, z_on_j = anatomy.locate(ws.seclens[cable_i], z)

				# Cable and point of the cable to work on
				point = {
						'cable': cable_i, 
						'section': j_sec, 
						'z': z_on_j

				# Ground it

# Once the classes are defined, define the types of electrodes
electrode_types = {
	"cuff": CuffElectrode, 
	"intracellular": IntracellularElectrodeSet, 
	"extracellular point": ExtracellularPointElectrode, 
	"extracellular points set": ExtracellularPointElectrodeSet, 
	"ground points set": GroundPointSet

# Functions

def create_electrodes():
	""" Create the electrodes which are specified in the 
	electrodes.json file """
	ws.electrodes = {}	
	for name, el in ws.electrodes_settings.items():
		print('Electrode: %s'%name)
		ws.electrodes[name] = electrode_types[el['type']](name)

def set_stimulation():
	""" Set up the stimulation protocols described in the json files. 
	This information is already in the electrodes attributes """
	for electrode in ws.electrodes.values():
		if electrode.role == "stimulation":

def set_recordings():
	""" Set up the stimulation protocols described in the json files. 
	This information is already in the electrodes attributes """

	# # Time
	# set_time_recording()

	# Electrode recordings
	for electrode in ws.electrodes.values():
		if electrode.role == 'recording':
			# Tell ws that there are electrode recordings
			ws.settings["electrode recordings"] = True
			# Set the recordings

def check_cuff_cover(zi):
	Check if a value over the z-axis is under the region of a cuff 
	for name, electrode in ws.electrodes.items():
		if electrode.type == 'cuff':
			if electrode.left_end <= zi <= electrode.rght_end:
				return name
	return None

def prepare_ground_paths():
	Calculate the variables that need to be used to connect the outer 
	axons to the distant ground, through the cuff insulators and the 

	# Shortennings
	nvt = ws.nvt
	container_rho = ws.container_settings['resistivity']
	container_diam = ws.container_settings['diameter']
	# Convert the resistivity to Ohm*um
	container_rho *= 1.e4

	# Bases of the cylinder: conductivity to ground
	# Conductances p.u.area in both cases (mho/cm2)
	# Just direct grounding
	ws.xg_distal = 1.e9
	ws.xg_proxim = 1.e9
	# Actually, I want to think that there's medium at the bases
	# This base will have the cross-sectional area of the container, 
	# not the nerve, so I need to use that to compute rg
	# Let's assume its thickness is the same as in the radial direction
	ws.xg_distal = []
	ws.xg_proxim = []
	# Iterate over cables
	for i in range(ws.counters['cables']):
		area = nvt.free_areas[i]
		# Area that corresponds to this cable for the base of the 
		# container
		# NOTE: This is just an approximation that would be ideal 
		# should the nerve be circular in its cross-section
		# If a nerve is not, it's still OK, because I use its 
		# circumcircle
		area_augmented = area * (container_diam / nvt.circumdiameter) ** 2
		# Conductance in mho
		# Note that the 'depth' (distance from segment to ground along 
		# the z-axis, which means through the saline bath) is the 
		# container's radius, according to our assumption above (ASSUMPTION *1)
		ws.xg_distal.append(area_augmented / (container_rho * 0.5 * container_diam))

	# Thickness of the containing medium. There's one for each electrode
	ws.container_settings['thickness'] = {}

	# Round wall of the cylinder

	# Area that the current from one segment to ground has to cross (um2)
	# Length of a segment times its portion of the perimeter
	area = (ws.length / ws.NAELC_nseg) \
	* np.pi * 0.5 * (nvt.circumdiameter + container_diam) / ws.npc
	# Resistance to ground p.u.area (Ohm*um2)
	# The distances must be in cm (they have to be converted from um)
	# No, it now must be in um again, becuase I've put container_rho in Ohm*um
	rg = {}
	rg['wout_cuff'] = container_rho * 0.5 * (container_diam - nvt.circumdiameter)
	rg['with_cuff'] = {}

	# Conductance to ground in mho (p.u.area (mho/um2) * area (um2))
	ws.xg_rwall = {}

	# Without cuff

	ws.xg_rwall['wout_cuff'] = area / rg['wout_cuff']

	# With cuff
	ws.xg_rwall['with_cuff'] = {}

	# Iterate over cuff electrodes
	for name, el in ws.electrodes.items():
		if el.type == 'cuff':

			# Add the thickness from the electrodes thickness
			ws.container_settings['thickness'][name] = \
				0.5 * container_diam - (nvt.circumradius \
				 + ws.electrodes_settings[name]['thickness'])

			# Round wall of the cylinder

			# Resistance to ground p.u.area (Ohm*um2)
			# The resistivity has to be converted from Ohm*cm to Ohm*um
			rg['with_cuff'][name] = \
				1.e4 * ws.electrodes_settings[name]['insulator resistivity'] \
				* ws.electrodes_settings[name]['thickness'] \
				+ container_rho * ws.container_settings['thickness'][name]

			# Conductance to ground in mho (p.u.area (mho/um2) * area (um2))

			# Conductance
			ws.xg_rwall['with_cuff'][name] = area / rg['with_cuff'][name]

def square_pulses(injections_, point, protocol):
	""" Assign current injections to the pads """
	# Injection segment in hoc
	seg = ws.sections[point['cable']][point['section']](point['z'])

	# IClamp object(s)
	for i, amp in enumerate(protocol['currents']):
		# Gather the data
		dur = protocol['pulse durations'][i]
		delay = protocol['pulse onset times'][i]
		# Create IClamp object and add its properties
		injection = h.IClamp(seg)
		# Intensity from uA to nA
		injection.amp = amp * 1.e3
		injection.dur = dur
		injection.delay = delay
		# Add this information as an attribute
			'injection id': ws.counters['current injections'], 
			'point': point, 
			'amp': amp, 
			'dur': dur, 
			'delay': delay
		# Update the counter
		ws.counters['current injections'] += 1
		# Add the duration and delay to ws.totaldur
		ws.totaldur = max(ws.totaldur, delay + dur)
		# Append it to ws
		ws.injections_info = ws.injections_info + injections_
		ws.log('Injection added to: %s, %s, %s, %s'%(str(point), str(amp), str(dur), str(delay)))

def ground_segment(point):
	""" This connects a segment directly to ground. It is obviously not 
	a current injection protocol, but still is a necessary kind of 
	protocol in case it's needed, I think. And this is a good place to 
	define it.
	Because the paths or connections to ground have already been 
	defined in this module previously and declared in h in 
	biophysics.electrics, this replaces whatever path to ground this 
	segment had, if any. If it didn't have it, it creates it. """

	# Corresponding segment in hoc
	seg = ws.sections[point['cable']][point['section']](point['z'])

	# Level or layer where to point to
	level = ws.cables[point['cable']].properties['extlayers'] - 1

	# Ground it directly
	seg.xg[level] = 1.e9

def set_recording(recordings_, point, var):
	""" Assign the recording vectors to each pad """

	# This is the contact segment in hoc
	seg = ws.sections[point['cable']][point['section']](point['z'])

	# Set up the recording information
			'recording id': ws.counters['recordings'], 
			'point': point, 
			'variable': var

	# Set up the recording object
	als.set_recording(seg, var)

Loading data, please wait...