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 was created by Miguel Capllonch Juan on 13 November 2018.
It handles most aspects of the system's biophysics

import os
import csv
import numpy as np
from neuron import h
from datetime import datetime
import planar as pl

import algebra as alg
import workspace as ws
import anatomy
import electrodes
import analysis as als
import tools as tls

class AxonGeneralModel():
	This is a class to put in all the possible properties of an 
	axon model type, including anatomical and biophysical general properties
	This class is only instantiated once per model type, so each 
	instance of a physical axon does not have to do everything 
	from scratch
	def __init__(self, model):
		self.model = model
		self.properties_ = ws.axonmodel_settings[model]

	def build_regressions(self):
		Build the necessary regression functions of the axon model from the 
		json files in the work space
		if 'regressions' in self.properties_.keys():
			regfile = os.path.join(ws.folders['models'], 
			self.regressions = alg.parameter_regressions(regfile)

class Cable:
	This is the basic class for any cable, be this an axon or a NAELC
	def __init__(self):

		# Counter. Add one to the class instances counter and save it 
		# as an id number as well
		ws.counters['cables'] += 1
		self.cable_number = ws.counters['cables'] - 1

		# Flag indicating if this cable has recordings
		self.has_recordings = False

		# Basic geometrical properties
		if ws.settings["nerve model"] == "resistor network":
			self.x = ws.nvt.x[self.cable_number]
			self.y = ws.nvt.y[self.cable_number]
			self.r = ws.nvt.r[self.cable_number]
		elif ws.settings["nerve model"] == "simple":
			self.x = ws.nvt.x[self.cable_number + ws.nvt.nNAELC]
			self.y = ws.nvt.y[self.cable_number + ws.nvt.nNAELC]
			self.r = ws.nvt.r[self.cable_number + ws.nvt.nNAELC]

		# Create properties dictionary = {}

		# Find the fascicle to which this cable belongs["fascicle"] = ws.nvt.cables_fascicles[self.cable_number]
	def finish_build(self):
		Perform the last steps in building this cable model

		# Append some properties to the corresponding dictionaries
		i = self.cable_number

		# Append the cell to the list of cables in the work space

		# Sections and lengths
		secs = list(self.cable_cell.all)
		ws.sections[i] = secs[:] 
		ws.sectypes[i] =['sctyps']
		ws.seclens[i] = [sec.L for sec in secs]
		# print('created seclens:', i)
		# Segments
		for j, sec in enumerate(secs):
			ws.secsegs[(i, j)] = list(sec.allseg())[1:-1]
		segs = tls.flattenlist([list(sec.allseg())[1:-1] for sec in ws.sections[i]])
		ws.segments[i] = segs[:]
		# Longitudinal profile (z-profiles)
		ws.zprofs[i] = anatomy.zprofile(self.cable_cell)
		# Maximum number of segments
		if len(ws.zprofs[i]) > ws.maxnseg:
			ws.maxnseg = len(ws.zprofs[i])

		# Store attributes['zprofile'] = ws.zprofs[i]

		# Intracellular area
			self.intracellular_area = np.pi * (0.5 *["axonD"]) ** 2
		except KeyError:
			# This is a NAELC or is an axon not configured to have axonD (which will be tackled if it happens)
			self.intracellular_area = 0.

	def set_recordings(self):
		Prepare the recording vectors for the desired variables
		# ws.log("Setting recordings for cable %i"%self.cable_number)

		# Update the flag indicating if this cable has recordings
		self.has_recordings = True

		# Variables to record from the settings
		rvars = ws.settings["data records"]["variables"]

		# Dictionary with all the recordings:
		# One entry per variable. Each entry contains a list of 
		# recording identifiers (one for each variable and segment of  
		# the cable)
		self.recordings = {}

		# Which sections to record
		where_to_record = ws.settings["data records"]["record sections"]

		# If this cable is a NAELC, record only vext[0]
		if self.cabletype == 'NAELC':
			rvars = ['vext[0]']
			where_to_record = ["all"]

		# Record the variables
		timings = []
		for var in rvars:

			# Create the list
			self.recordings[var] = []

			# Create the recording vectors
			for i, seg in enumerate(ws.segments[self.cable_number]):

				record = False
				t0 =

				if where_to_record[0] == "all":
					record = True
					# for sectypename in where_to_record:
					# 	if["sctyps"][i] == sectypename:
							# record = True
					if["sctyps"][i] in where_to_record:
						record = True

				if record:
						"number": ws.counters['recordings'], 
						"segment": str(seg)
					als.set_recording(seg, var)

				t1 =
				timings.append((t1 - t0).total_seconds())

		ws.log("\tAverage time needed by als.set_recording for cable %i: %f ms. Max: %f ms"%(self.cable_number, 1e3 * np.array(timings).mean(), 1e3 * np.array(timings).max()))

class Axon(Cable):
	This is a class to put in all the possible properties of an 
	axon model, including anatomy and biophysics
	This class inherits its properties anatomical and biophysical 
	general properties from THE INSTANCE of AxonGeneralModel which 
	describes the model type
	def __init__(self, model):

		# Inheritance
		# Type of cable
		self.cabletype = 'Axon'
		# print('***Creating axon with the model %s'%model)
		# Counter
		ws.counters['axons'] += 1
		# Specific identifier; specific for this type of cable
		self.specific_number = ws.counters['axons'] - 1
		# Properties
		# Put AxonGeneralModel.properties_ into['fiber_radius'] = self.r['fiberD'] = 2. *['fiber_radius']
		# Axon model
		#['model'] = ws.nvt.models[self.cable_number]

	def build(self):
		General class that chooses between the build methods for 
		myelinated or unmyelinated axons
		build_methods = {
			'myelinated': self.build_myelinated, 
			'unmyelinated': self.build_unmyelinated

		# Finish building this cable and saving its properties

	def build_unmyelinated(self):
		Build the axon as unmyelinated
		nseg = anatomy.nseg_dlambda_rule(
			)['nseg'] = nseg
		self.cable_cell = h.UnMyelAxon(['sctyps'] = np.zeros(nseg,

	def build_myelinated(self):
		Build the necessary properties of the axon model from the 
		json files in the work space

		# Variables that depend on a regression from some stored data
		iv = 'independent variable'
		rgs = self.regressions
		for k, func in rgs.items():
			if k != iv:
				x =[rgs[iv]][k] = func(x)

		# Other properties:

		# Variable name shortenings
		deltax =['deltax']

		# 1. STIN_length
		# This is still model-specific hard-coding
		if self.model == 'MRG' or self.model == 'gaines_sensory':
			nodl =['nodelength']
			pl1 =['paralength1']
			pl2 =['paralength2']['STIN_length'] = deltax - nodl - 2. * (pl1 + pl2)
		elif self.model == 'MRG_STINonly':
			nodl =['nodelength']['STIN_length'] = deltax - nodl

		# Start point along the z-axis
		if ws.settings["nerve model"] == "resistor network":["start"] = deltax * ws.nvt.start_positions[self.cable_number]
		elif ws.settings["nerve model"] == "simple":["start"] = deltax * ws.nvt.start_positions[self.cable_number + ws.nvt.nNAELC]

		# 3. Total internodal length['in_l'] = deltax -['nodelength']

		# 4. Update the minimum length of one chunk of sections
		# and of the node to fiber ratio, g
		ws.deltax_min = min((ws.deltax_min, deltax))
		ws.g_min = min((ws.g_min,['g']))

		# Build the axon in hoc
		sctyps, section_counter, n_max, lengths = anatomy.all_vars(self)

		# Lengths to hoc
		lengths_hoc = {}
		for key, value in lengths.items():
			lengths_hoc[key] = h.Vector(value)

		# Number of sections of each type and section types
		nsecs = sum(section_counter.values())
		sectypes = h.List()
		for st in sctyps:

		# Store variables as attributes inside the 'properties' dictionary['nsecs'] = nsecs['sectypes'] = sectypes['section_counter'] = section_counter['mx_n'] = n_max['lengths'] = lengths_hoc

		# Axon in hoc
		if self.model == 'MRG':
			# print(' ----------- Setting up %s'%self.model)
			self.cable_cell = h.MRGMyelAxon(
		elif self.model == 'gaines_sensory':
			# print(' ----------- Setting up %s'%self.model)
			self.cable_cell = h.Gaines2018Sensory(
		# z-profile
		zp = anatomy.zprofile(self.cable_cell)

		# If this is the smallest axon, save its longitudinal profile
		if['fiber_radius'] == ws.raxmin:
			ws.zlayers_min = zp

		# Number of segments on each section of this axon
		this_nseg = [sec.nseg for sec in self.cable_cell.all]

		# I need sctyps to reflect all the segments, not just the sections
		sctyps_str = [this_nseg[ii] * [sctyps[ii].sectype] for ii in range(len(sctyps))]
		sctyps_str = np.array(tls.flattenlist(sctyps_str))

		# Get the z-profiles of the nodes of Ranvier only
		zprofs_RN = zp[np.where(sctyps_str == 'node')[0]]
		ws.zprofs_RN[self.cable_number] = zprofs_RN

		# Store these properties['zprofs_RN'] = zprofs_RN['sctyps'] = sctyps_str

class NAELC(Cable):
	This is a Non-Axonal Cable 
	def __init__(self, *dummyargs):
		# Inheritance
		# Type of cable
		self.cabletype = 'NAELC'
		# Counter
		ws.counters['NAELC'] += 1
		# Specific identifier; specific for this type of cable
		self.specific_number = ws.counters['NAELC'] - 1
		# Properties ={
					'diam': ws.cdiam, 
					'length': ws.length, 
					'extlayers': ws.NAELC_extlayers

	def build(self):
		Build the cable into the hoc space
		"""['nseg'] = ws.NAELC_nseg['sctyps'] = np.zeros(ws.NAELC_nseg,
		self.cable_cell = h.Wire(

		# Finish building this cable and saving its properties

class ResistorNetwork():
	This is a class for the Resistor Network that describes the 
	interactions between all axons
	def __init__(self):
	def get_resistances(self):
		""" Calculate the values of the transverse resistances """

		self.resistors = {}
		self.res_lengs = {}
		for pair in ws.pairs:
			i, j = pair
			# Calculation of the resistor's value: 
			# shape factor of the resistor (adimensional)
			# and value in MOhm*cm
			l = ws.nvt.len_con[pair]
			if ws.EC["resistor network"]["interaxonal connections"] == "membranes":
				l -= (ws.nvt.r[i] + ws.nvt.r[j])
			elif ws.EC["resistor network"]["interaxonal connections"] == "centers":
			# Width
			w = ws.nvt.len_seg[pair]
			# Shape factor (adimensional (um / um))
			shape_factor = l / w
			# Large shape factors mean large resistances; not interested
			if shape_factor < 100:

				# Resistance value
				res = ws.rho['endo_trans'] * shape_factor

				# Perineurium
				resist_perineurium = 0.

				# If the pair is at opposite sides of the perineurium, 
				# add the perineurium's resistance p.u.l or 
				# resistivity (Ohm*cm)
				if ('NAELC' in ws.nvt.cables[i] and 'xon' in ws.nvt.cables[j]) \
				or ('NAELC' in ws.nvt.cables[j] and 'xon' in ws.nvt.cables[i]):
					resist_perineurium += ws.resist_perineurium

				# Also, that happens when two axons are in connected but 
				# in different fascicles
				# if 'xon' in ws.nvt.cables[i] and 'xon' in ws.nvt.cables[j]:
				if ws.cables[i].cabletype == "Axon" and ws.cables[j].cabletype == "Axon":
					# ws.log("Checking if there\'s double layer of perineurium for axons %i and %i:"%(ws.cables[i].specific_number, ws.cables[j].specific_number))
					# ws.log("%s, %s"%(ws.cables[i].properties["fascicle"], ws.cables[j].properties["fascicle"]))
					if ws.cables[i].properties["fascicle"] != ws.cables[j].properties["fascicle"]:
						resist_perineurium += 2. * ws.resist_perineurium
						# ws.log("\tAdding double perineurium between axons %i and %i"%(ws.cables[i].specific_number, ws.cables[j].specific_number))
						# ws.log("\tThese axons are placed at (%f, %f) um and (%f, %f) um, respectively"%(ws.cables[i].x, ws.cables[i].y, ws.cables[j].x, ws.cables[j].y))

				# Finally, add it up
				res += resist_perineurium / w
				# Store in the dictionary
				self.resistors[pair] = res

	def connect_resistors(self):
		""" Connect the transverse resistors with ephapses in NEURON """

		def locate_and_create(pair, positions):
			""" For a pair of cables and a set of resistor positions, 
			locate, create and connect those resistors to the cables """

			# Pair indices
			i, j = pair

			# Get resistors' lengths
			rl = get_rl(pair, positions)

			# Iterate over resistors positions
			for zr, l in zip(positions, rl):

				# For each cable, see where this falls

				# On cable i:
				i_sec, z_on_i = anatomy.locate(ws.seclens[i], zr)
				seci = ws.sections[i][i_sec]
				sgi = seci(z_on_i)

				# On cable j:
				j_sec, z_on_j = anatomy.locate(ws.seclens[j], zr)
				secj = ws.sections[j][j_sec]
				sgj = secj(z_on_j)

				# Value of the resistor (MOhm)
				ephres = res / l
				# Create the resistor
				create_eph(ephres, sgi, sgj, seci, secj, i, j, 
				# ws.log("Created a transverse resistor at z = %f um with value %f MOhm"%(zr, ephres))

		ws.eph_counter = 0
		ws.nrn_res = {}
		for pair, res in self.resistors.items():
			i, j = pair

			# Cable
			cable_cell_i = ws.cables_hoc[i]
			cable_cell_j = ws.cables_hoc[j]

			# I have to connect the cables in a staggered manner

			# If I am only taking the nodes of Ranvier into account
			# and the connection is to be set between two axons:
			if ws.EC["resistor network"]["transverse resistor locations"] == "regular":

				# Separation between transverse resistors
				sep = float(ws.EC["resistor network"]["resistors separation"])

				# Positions of the resistors
				# res_positions = np.arange(0, ws.length, sep)
				res_positions = np.arange(0.5, ws.length, sep)

				# Create and connect the resistors
				locate_and_create(pair, res_positions)

			elif ws.EC["resistor network"]["transverse resistor locations"] == "only nodes of Ranvier" and \
			('xon' in ws.cables[i].cabletype and 'xon' in ws.cables[j].cabletype):

				# Get resistors' positions and lengths (um)

				# Merge the ws.zprofs_RN of the two axons
				# This array will contain the postions of the resistors 
				# over the z-axis
				zp_pair = np.append(ws.zprofs_RN[i], ws.zprofs_RN[j])
				# Remove repeated values (can --will, actually, happen if 
				# there's aligned nodes)
				zp_pair = np.unique(zp_pair)

				# Create and connect the resistors
				locate_and_create(pair, zp_pair)

				# If I am taking all nodes into account, the compartments of only one
				# of the two cables (arbitrarily chosen; cable 'i' by default) 
				# are plugged to transverse resistors. These are then connected 
				# to the second cable wherever they must, regardless of the 
				# topology of this second cable.

				# Get resistors' lengths (um)
				# From what is explained above, in this case the z-profile of one
				# of the cables only is used as zp_pair

				rl = get_rl(pair, ws.zprofs[i])
				# Iterate over segments
				for ii, sgi in enumerate(ws.segments[i]):

					# Get the section and segment of cable j which sgi lies against
					zi = ws.zprofs[i][ii]
					j_sec, z_ioverj = anatomy.locate(ws.seclens[j], zi)
					seci = sgi.sec
					secj = ws.sections[j][j_sec]
					sgj = secj(z_ioverj)

					# Value of the resistor (MOhm)
					ephres = res / rl[ii]
					# Create the resistor
					create_eph(ephres, sgi, sgj, seci, secj, i, j, 

# Functions

def finish_fascicles():
	Finish giving basic information to the fascicles once the cables have been set up and the resistivities have been adjusted (if they have)
	for k, fas in ws.nvt.fascicles.items():
		# Intracellular area
		# Print properties
	# Now print the properties of the nerve

def get_rl(pair, zp_pair):
	Get the lengths (all over the z-axis) of all the transverse
	resistors for a pair of cables
	They will be used to obtain the values of each transverse resistor in MOhm.

	# Distances between nodes
	dz = zp_pair[1:] - zp_pair[:-1]
	# Get the lengths of the resistors (um)
	rl = 0.5 * (dz[1:] + dz[:-1])
	# Add one element on the left: first resistor
	# and one on the right: last resistor
	rl = np.append(np.array([zp_pair[0] + 0.5 * dz[0]]), rl)
	rl = np.append(rl, np.array([0.5 * dz[-1] + ws.length - zp_pair[-1]]))
	# Store its vale
	ws.rnet.res_lengs[pair] = rl
	return rl

# def create_eph(ephres, sgi, sgj, seci, secj, i, j, nrn_res, zr):
def create_eph(ephres, sgi, sgj, seci, secj, i, j, zr):
	Create a resistor that acts as an ephaptic connection

	# Choose the layer depending on the type of fiber or cable
	# Layer where to point to
	layeri = ws.cables[i].properties['extlayers']
	layerj = ws.cables[j].properties['extlayers']

	# If zr is not in a cuff, we will make the ephaptic connection
	# super resistive and hence disappear
	if not ws.EC["resistor network"]["outside cuffs"]:
		ephres = disconnect_under_cuff(zr, ephres)

	# If the connection is too large, do not implement it
	if ephres < 1.e2:

		eph = h.EphapDBA(2)
		eph.ex0_loc(0, sgi.x, layeri, sec=seci)
		eph.ex1_loc(1, sgj.x, layerj, sec=secj)
		eph.add_submatrix(0, 1, ephres)
		# Disconnect segments from ground
		# This is redundant, but just in case
		sgi.xg[layeri - 1] = 0.
		sgj.xg[layerj - 1] = 0.

		# Add to dictionary
		ws.nrn_res = tls.append_items(ws.nrn_res, (i, j), [eph])
		ws.eph_counter += 1

def disconnect_under_cuff(zi, ephres):
	Check if a transverse resistor is not under the region of a cuff 
	electrode and if so, give it a highest value so it's not connected
	if electrodes.check_cuff_cover(zi) is None:
		ephres *= ws.rho['SuperResistive'] / ws.rho['endo_trans']
	return ephres

def build_cables():
	Build all the cables, NAELC and fibers

	ws.log("Building the cables in NEURON")
	t0 =

	# Dictionaries where to store information in the workspace
	ws.zprofs_RN = {}
	ws.sections = {}
	ws.sectypes = {}
	ws.secsegs = {}
	ws.seclens = {}
	ws.zprofs = {}
	ws.segments = {}
	cables = []
	ws.cables_hoc = h.List()

	# New building method
	# First, axon model properties
	# Axon general model(s)
	ws.ax_gral_models = {}
	# ws.ax_gral_models[ws.settings['axon model']] = AxonGeneralModel(ws.settings['axon model'])
	for m in ws.nvt.models_set:
		if 'NAELC' not in m:
			ws.ax_gral_models[m] = AxonGeneralModel(m)

	# Classes of cables
	ws.cable_classes = {
		"Axon": Axon

	# First, instantiate all the cables.
	# Then, build axons. From them, obtain the desired nseg for NAELC.
	# Finally, build the NAELC using this value

	# Instantiate everything
	nvt = ws.nvt
	for i in range(

		# Create the cable(s)

		if ws.settings["nerve model"] == "resistor network":
			# Use NAELC
			# cable = ws.cable_classes[nvt.cables[i]](ws.settings['axon model'])
			cable = ws.cable_classes[nvt.cables[i]](nvt.models[i])
			# Build it only if it is an axon
			# The NAELC will be built later on
			if nvt.cables[i] == 'Axon':
				# print('building', i, nvt.cables[i])
			# Append it

		elif ws.settings["nerve model"] == "simple":
			# Don't use NAELC
			if nvt.cables[i] == 'Axon':
				# cable = ws.cable_classes[nvt.cables[i]](ws.settings['axon model'])
				cable = ws.cable_classes[nvt.cables[i]](nvt.models[i])
				# Build the axon
				# Append it
	# Obtain nseg for the NAELC if there are axons
	if ws.counters['axons'] > 0:
		# There's axons
		daxmin = 2. * ws.raxmin
		if ws.EC["resistor network"]["transverse resistor locations"] == "only nodes of Ranvier":
			nseg = int(ws.length / ws.deltax_min) + 1
		elif ws.EC["resistor network"]["transverse resistor locations"] == "regular":
			nseg = int(ws.length / float(ws.EC["resistor network"]["resistors separation"])) + 1
			# nseg for an internode
			nseg = anatomy.nseg_dlambda_rule(ws.deltax_min, ws.g_min * daxmin, 
				cables[ws.counters['NAELC']].properties['rhoa'], cm=cables[ws.counters['NAELC']].properties['cm'])
			# nseg for the whole nerve length
			nseg = int(ws.length / (ws.deltax_min / nseg)) + 1

		# Make sure it's an odd number
		# This adds 2 if it's odd, but it doesn't matter; it's even better
		nseg += nseg % 2 + 1

		# If we're here, there are no axons
		nseg = ws.nseg

	ws.NAELC_nseg = nseg
	ws.log("NAELC segment length = %f"%(ws.length / nseg))
	ws.log("nseg: %i"%nseg)
	if ws.counters['axons'] > 0:
		ws.log("Min. deltax: %f"%ws.deltax_min)

	# Now build the NAELC according to the value of nseg 
	# obtained from the axons
	for i in range(ws.counters['NAELC']):
		# print('building', i, cables[i])

	# print(ws.nvt.cables)

	# Store the cables in the workspace
	ws.cables = cables

	# Maximum number of segments. This is a number to be found
	ws.maxnseg = 0

	t1 =
	ws.log("Created %i cables in %i seconds"%(, (t1 - t0).total_seconds()))

def electrics():
	In this function, xraxial is set for every axons and all 
	connections to ground are made

	# The following loop over cables and segments and calculates:
	# - The cross-sectional area of each cable and therefore:
	#  - its value for xraxial
	#  - its connections to ground at the ends of the nerve

	for i, seglist in ws.segments.items():
		# Level or layer where to point to
		level = ws.cables[i].properties['extlayers'] - 1

		# Cross-sectional area and xraxial
		# Cross-sectional area assigned to the cable (cm2)
		carea = ws.nvt.free_areas[i] * 1.e-8

		# Attribute xraxial
		for j, segment in enumerate(seglist):
			# First and foremost, disconnect ALL segments from ground
			segment.xg[level] = 0.

			# If there's no EC outside the cuff-covered regions,
			# check if this point is under a cuff and connect it to ground
			if not ws.EC["resistor network"]["outside cuffs"]:
				cuff_cover = electrodes.check_cuff_cover(ws.zprofs[i][j])
				if cuff_cover is None:
					segment.xg[level] = 1.e9

			# Attribute xraxial (MOhm/cm; so I have to convert rho (in MOhm*um) 
			# to MOhm*cm)

			# Longitudinal resistivity for this cable as a weighted average between epineurium and endoneurium (the endoneurium is, in turn, a weighted average between all the fascicles intersecting the cross-sectional area of this cable)
			ctendo = ws.nvt.cables_tissues[i]['endoneurium']
			ctepi = ws.nvt.cables_tissues[i]['epineurium']
			fas_ = ws.nvt.fascicles
			resistivity = ws.rho['endo_longt'] * sum([ctendo[k] for k in fas_]) + ctepi * ws.rho['epi_longt']

			# Apply the resistivity in xraxial
			xrx =  1.e-4 * resistivity / carea
			segment.xraxial[level] = xrx

		# Connect to ground all cables at both ends of the nerve
		if ws.EC["resistor network"]["outside cuffs"]:
			segleft = seglist[0]
			segrght = seglist[-1]
			# Note: The areas of each segment need to be in cm2
			segleft.xg[level] = ws.xg_distal[i] / (segleft.area() * 1.e-8)
			segrght.xg[level] = ws.xg_proxim[i] / (segrght.area() * 1.e-8)
			# In this case, just ground the ends
			# This is regardless of the ends being cuffed or not; 
			# these are at the base of the cylinder
			seglist[0].sec(0.).xg[level] = 1.e9
			seglist[-1].sec(1.).xg[level] = 1.e9

	# Connect grounds to the nerve's walls (not the ends; these are already sorted)
	# The xg[level] in NEURON for each segment is in mho/cm2
	# All points in the contour are connected to a far ground
	for i in np.arange(ws.npc):
		# Level or layer where to point to
		level = ws.cables[i].properties['extlayers'] - 1
		for j, segment in enumerate(ws.segments[i]):

			# Convert area to cm2 (it's in um2)
			# Cable's cylindrical wall area
			area = segment.area() * 1.e-8

			# Set xg
			# This depends on whether a segment is under the cuff or not
			cuff_cover = electrodes.check_cuff_cover(ws.zprofs[i][j])
			if cuff_cover is not None:
				# xg must be in mho/cm2, so xg_rwall should be in mho; it is

				# If this segment is at one end of the nerve, it already 
				# has one connection to ground that needs to be added
				if segment in (ws.segments[i][0], ws.segments[i][-1]):
					segment.xg[level] += \
						ws.xg_rwall['with_cuff'][cuff_cover] / area
					segment.xg[level] = \
						ws.xg_rwall['with_cuff'][cuff_cover] / area
				if ws.EC["resistor network"]["outside cuffs"]:
					# If this segment is at one end of the nerve, it already 
					# has one connection to ground that needs to be added
					if segment in (ws.segments[i][0], ws.segments[i][-1]):
						segment.xg[level] += \
							ws.xg_rwall['wout_cuff'] / area
						segment.xg[level] = \
							ws.xg_rwall['wout_cuff'] / area
					# If this is not meant to be used, just ground it
					# (in a future version, I will simply crop these
					# bits of NAELC outside the cuffs)
					segment.xg[level] = 1.e9

def set_direct_extracellular_stimulation():
	""" Set up the stimulation over e_extracellular on the axons """

	# File to read the data from
	espath = os.path.join(ws.folders['data/load'], 'extstim.csv')

	# Create the time window
	# All ones for now
	# It will be cropped when the information for each pulse is fetched
	time_window = np.ones(ws.nt)

	# Time series of the extracellular field for the axons
	stim_tseries = []

	# Dictionaries with the arrays and more info
	pulses = {}

	# Open and read
	with open(espath, 'r') as f:
		reader = csv.reader(f)

		# Read line by line
		for row in reader:

			# First element of each row
			key = row[0]

			# New pulse
			if key == 'Pulse':
				k = int(row[1])
				pulses[k] = {}
				# Reset the time window
				time_window[:] = 1.


				# Parse the information in the line
					# See if this is an axon
					i = int(key)
				except ValueError:
					# If not, see if row[1] is a numerical value
						pulses[k][key] = float(row[1])
					except ValueError:
						# If not, it's a string; just save it as such 
						pulses[k][key] = row[1]
					# If it is an axon, save the array
					field = np.array([float(x) for x in row[1:]])

					# If we're here, this means we already have all the information
					# Give values to the time window (only once)
					if i == 0:
						if pulses[k]['Type'] == 'square pulse':
							time_window[np.where(ws.tarray < pulses[k]['Delay'])] = 0.
							time_window[np.where(ws.tarray >= pulses[k]['Delay'] + pulses[k]['Duration'])] = 0.

					# If this is the first pulse, create the time series for 
					# this axon
					if k == 0:
						stim_tseries.append(np.zeros((len(field), len(time_window))))

					# Now create the time series of the extracellular field 
					# for each segment of this axon and play it in NEURON
					for j, seg in enumerate(ws.segments[i]):
						stim_tseries[i][j] += field[j] * time_window
	# Now that each axon has its full stimulation time 
	# series, play it (fts: 'field time series')

	ws.stim_field_vectors = []

	for i in list(ws.segments.keys())[ws.counters['NAELC']:]:

		for j, seg in enumerate(ws.segments[i]):

			# Create and play the time series
			fts = h.Vector(stim_tseries[i][j]), h.dt)

			# Add it to the list
			# This is necessary for this to work
			# Otherwise, fts is lost

Loading data, please wait...