Error at implementing the ephap.zip example in Python

Moderator: wwlytton

Post Reply
mcapllonchj

Error at implementing the ephap.zip example in Python

Post by mcapllonchj »

Hello,

I have translated part of the code from ephap.zip into Python.
If I create a class ephap like:

Code: Select all

h.load_file("ephap.hoc")
ephap = h.Ephap()
, everything works fine, of course, but I wanted to have as much code as I could in Python, ideally all of it, so I have tried to create a Python class equivalent to h.Ephap. I don't know how to implement it entirely in Python, so I created the simplest class in hoc I could create and loaded it into neuron.h in Python. The hoc file is (ephaptic_forpython.hoc):

Code: Select all

begintemplate EphapPython
proc init() {}
endtemplate EphapPython
and the python script is:

Code: Select all

"""
Ephaptic coupling example by Ted Carnevale implemented almost exclusively 
in Python by Miguel Capllonch Juan on August 2017
"""

import matplotlib.pyplot as plt
from neuron import h, gui, hclass

# The written class EphapPython needs to inherit from a class previously defined
# in hoc, so I have created the simplest possible class I managed in this
 # hoc file below
h.load_file("ephaptic_forpython.hoc")

class EphapPython(hclass(h.EphapPython)):
  """
  docstring for EphapPython
  """
  def __init__(self):
    self._n = 7
    self._c = h.Matrix(self._n, self._n, 2)
    self._g = h.Matrix(self._n, self._n, 2)
    self._b = h.Vector(self._n)
    self._y = h.Vector(self._n)
    self._y0 = h.Vector(self._n)
    self._ncellport = 4
    self._xloc = h.Vector(self._ncellport)
    self._layer = h.Vector(self._ncellport)
    self._afac = h.Vector(self._ncellport)
    self._sf = h.StringFunctions()

    # Element parameters
    # Transfer resistance (MOhm)
    self.Rx = 10.
    # Named voltage nodes and initial conditions
    # I think I don't need this, for the moment... and it gives me an error
    # self._sf.alias(self, "ex0", self._y.x[0])

    # I think I have to add this:
    self._srs = self._ncellport * [None]

  def exx_loc_basics(self, ix, x, sec, layer):
    self._srs[ix] = h.SectionRef(sec)
    self._xloc.x[ix] = x
    self._layer.x[ix] = layer

  def ex0_loc(self, x, sec):
    self.exx_loc_basics(3, x, sec, 0)
    self.exx_loc_basics(0, x, sec, 1)

  def ex1_loc(self, x, sec):
    self.exx_loc_basics(2, x, sec, 0)
    self.exx_loc_basics(1, x, sec, 1)

  def refill(self):
    for i in xrange(self._ncellport):
      a = self._srs[i].sec(self._xloc[i]).area()
      if a == 0: a = 100
      self._afac.x[i] = -100. / a

    self._g.setval(0, 5, (-1.0) * self._afac.x[0])
    self._g.setval(1, 6, (1.0) * self._afac.x[1])
    self._g.setval(4, 4, -1.0 / self.Rx)
    self._g.setval(4, 5, 1.0)
    self._g.setval(4, 6, -1.0)
    self._g.setval(5, 0, -1.0)
    self._g.setval(5, 4, 1.0)
    self._g.setval(6, 1, 1.0)
    self._g.setval(6, 4, -1.0)

  def install(self):
    self._sl = h.SectionList()
    for i in xrange(self._ncellport):
      self._sl.append(self._srs[i].sec)
    self.refill()
    self._y.copy(self._y0)
    self._lm = h.LinearMechanism(self._c, self._g, self._y, self._b, self._sl,
      self._xloc, self._layer)

###############################################################################

def initialize():
	"""Initialize the simulation for the first time step"""
	h.finitialize(v_init)
	h.fcurrent()

def go():
	"""Main driver for the simulation"""
	initialize()
	h.run()

# Define a BScell class
h.load_file("bscell.hoc")

# Create the two cells and add the ephaptic coupling mechanism
cells = h.List()
for i in xrange(2):
	cells.append(h.BScell())
	cells[i].position(0, 100*i, 0)

	# IMPORTANT!!!
	# Disconnect axon.xg(0.5) from ground
	cells[i].axon(0.5).xg[0] = 0

# Set up ephaptic coupling

# Option 1: Using ehpap.hoc from Carnevale
# h.load_file("ephap.hoc")
# ephap = h.Ephap()

# Option 2: Using my implemented class in Python
ephap = EphapPython()

ephap.ex0_loc(0.5, sec=cells[0].axon)
ephap.ex1_loc(0.5, sec=cells[1].axon)
ephap.install()
ephap.Rx = 18.81e3 # Threshold value for trigerring a spike in cell 1
ephap.install()

# Current injection
stim = h.IClamp(cells[0].soma(0))
stim.amp = 1
stim.delay = 10
stim.dur = 0.1

# Initial conditions and time parameters
v_init = -64.974
h.tstop = 20.
h.dt = 0.025

# Vectors for recording
tvec = h.Vector()
v0s = h.Vector()
v1s = h.Vector()
v0a = h.Vector()
v1a = h.Vector()

# Record
tvec.record(h._ref_t)
v0s.record(cells[0].soma(0.5)._ref_v)
v1s.record(cells[1].soma(0.5)._ref_v)
v0a.record(cells[0].axon(0.5)._ref_v)
v1a.record(cells[1].axon(0.5)._ref_v)

# Initialize and run
go()

# Plot results
fig, ax = plt.subplots()
ax.plot(tvec, v0s, 'red', label='Cell 0 soma')
ax.plot(tvec, v0a, 'red', ls='--', label='Cell 0 axon')
ax.plot(tvec, v1s, 'blue', label='Cell 1 soma')
ax.plot(tvec, v1a, 'blue', ls='--', label='Cell 1 axon')
ax.legend()
plt.show()
With this, I get the following error:

Code: Select all

NEURON: spFactor error: Singular
 near line 0
 ^
        fadvance()
      advance()
    step()
  continuerun(20)
and others
The error comes at h.run(), and it is caused by the presence of the linear mechanism. If I comment the line with LinearMechanism, the simulation (without ephaptic coupling, obviously) works successfully.
Another odd thing that came up to happen is that defining a dummy linear mechanism outside the class or as global yields no error message and NEURON says "finished", but the program actually does not go further than h.finitialize().
Would you mind point me out to what I did wrong? It is the first time I handle the hclass factory, so I might have messed something up, maybe?
Thank you!
Post Reply