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()
Code: Select all
begintemplate EphapPython
proc init() {}
endtemplate EphapPython
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()
Code: Select all
NEURON: spFactor error: Singular
near line 0
^
fadvance()
advance()
step()
continuerun(20)
and others
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!