h.finitialize() fails with a hocobj call error

Anything that doesn't fit elsewhere.
Post Reply
neuron121!
Posts: 8
Joined: Thu Mar 31, 2022 11:00 pm

h.finitialize() fails with a hocobj call error

Post by neuron121! »

Hello,
I got some python code sent to me that works on their end, but it gets the following error on my end:

Code: Select all

# lots of similar lines
...
...
C:\Users\jdami\AppData\Local\Temp\ccWASCyx.s:1829: Error: bad register name `%rbp)'
C:\Users\jdami\AppData\Local\Temp\ccWASCyx.s:1830: Error: bad register name `%rsp'
C:\Users\jdami\AppData\Local\Temp\ccWASCyx.s:1831: Error: bad register name `%rbp'
C:\Users\jdami\AppData\Local\Temp\ccWASCyx.s:1833: Warning: zero assumed for missing expression
C:\Users\jdami\AppData\Local\Temp\ccWASCyx.s:1833: Warning: zero assumed for missing expression
Traceback (most recent call last):
  File "c:\nrn\lib\python\neuron\rxd\rxd.py", line 1872, in _init
    _compile_reactions()
  File "c:\nrn\lib\python\neuron\rxd\rxd.py", line 1484, in _compile_reactions
    _c_compile(fxn_string),
  File "c:\nrn\lib\python\neuron\rxd\rxd.py", line 564, in _c_compile
    dll = ctypes.cdll["%s.so" % os.path.abspath(filename)]
  File "C:\Users\jdami\AppData\Local\Programs\Python\Python310\lib\ctypes\__init__.py", line 449, in __getitem__
    return getattr(self, name)
  File "C:\Users\jdami\AppData\Local\Programs\Python\Python310\lib\ctypes\__init__.py", line 444, in __getattr__
    dll = self._dlltype(name)
  File "C:\Users\jdami\AppData\Local\Programs\Python\Python310\lib\ctypes\__init__.py", line 374, in __init__
FileNotFoundError: Could not find module 'C:\Users\jdami\Documents\Year 4\Lazzi\Colab\ITEMS_AxonGrowth\rxddll9e37c954-e82a-11ec-b65e-7470fdcdedcf.so' (or one of its dependencies). Try using the full path with constructor syntax.

NEURON: Python Callback failed
 near line 0
 ^
        finitialize()
Traceback (most recent call last):
    vgcc.stimulate(vgcc.somaList[0], vgcc.somaList[0], 0.5, 0, 0.5, 60, 100)
  File "c:\Users\jdami\Documents\Year 4\Lazzi\Colab\ITEMS_AxonGrowth\CaSignalingPMCAClass.py", line 588, in stimulate
    h.finitialize(-70 * mV)
RuntimeError: hocobj_call error
Basically, the function h.finitialize() breaks because there is an hocobj call error. I have no idea how to fix this. I'm including all of the code for reference. All needed files are there, but I'm only showing the main file and the one where the morphology is declared:

Main file:

Code: Select all

from neuron import h, rxd
import numpy as np
import json
from matplotlib import pyplot
from neuron.units import mM, mV, nM, uM
from neuron.rxd import v
from neuron.rxd.rxdmath import vtrap, exp, log, tanh, fabs
import RGCMorphology as rgcMorph
# import calcium_optimization_init_conds as cbdInit
import os
h.load_file('stdrun.hoc')
h.load_file("import3d.hoc")
h.load_file("stdlib.hoc")

class VGCC_2:

    def __init__(self, neuronMorph, JSONFile):

        self.axonList = [sec for sec in neuronMorph.axon]
        self.somaList = [sec for sec in neuronMorph.soma]
        self.apicalList = [sec for sec in neuronMorph.soma]
        self.all = self.axonList + self.somaList + self.apicalList

        # self.somaList[0].pt3dadd(-300,0,0,self.somaList[0].diam)
        # self.somaList[0].pt3dadd(-300+self.somaList[0].L, 0, 0, self.somaList[0].diam)

        with open(JSONFile, 'r') as fileObj:
            self.rxdDict = json.load(fileObj)

        self.findExtrema()
        self.padding = 100

    def create_neuron(self):
        '''
        THINGS TO THINK ABOUT: Do we really need a function that creates a new neuron? Perhaps we can pass a neuron as a parameter
        and implement VGCC's in that neuron, so we may not even need this function.

        :param self:
        :return:
        '''

        self.axon = h.Section(name="axon")
        self.soma = h.Section(name="soma")
        self.soma.pt3dclear()
        self.soma.pt3dadd(-90, 0, 0, 30)
        self.soma.pt3dadd(-60, 0, 0, 30)
        self.axon.nseg = 11
        self.soma.nseg = 11
        self.axon.connect(self.soma)

        self.axon.diam = 10
        self.axon.L = 100

        return self.soma

    def findExtrema(self):
        '''
        Finds the extrema of the morphology, and stores them into variables xmin, ymin, zmin, xmax, ymax, and zmax.
        These coordinates will be used in defining the extracellular space in 3D using rxd.Extracellular()

        :return: NONE
        '''

        self.xmin = self.ymin = self.zmin = self.xmax = self.ymax = self.zmax = None
        for sec in self.all:
            n3d = sec.n3d()
            xList = [sec.x3d(i) for i in range(n3d)]
            yList = [sec.y3d(i) for i in range(n3d)]
            zList = [sec.z3d(i) for i in range(n3d)]

            xMinSec, yMinSec, zMinSec = min(xList), min(yList), min(zList)
            xMaxSec, yMaxSec, zMaxSec = max(xList), max(yList), max(zList)

            if self.xmin is None:
                self.xmin, self.ymin, self.zmin = xMinSec, yMinSec, zMinSec
                self.xmax, self.ymax, self.zmax = xMaxSec, yMaxSec, zMaxSec
            else:
                self.xmin, self.ymin, self.zmin = min(self.xmin,xMinSec), min(self.ymin,yMinSec), min(self.zmin,zMinSec)
                self.xmax, self.ymax, self.zmax = max(self.xmax,xMaxSec), max(self.ymax,yMaxSec), max(self.zmax,zMaxSec)





    def concInit(self, cytc, exc, erc):
        return lambda nd: exc if nd.region == self.ecs else (cytc if nd.region == self.cyt else erc)

    def implement_RXD(self, neuronMorph):
        ########################################
        # Parameters

        # Initial ion concentrations
        naConcInitCyt = self.rxdDict["ions"]["na"]["cyt"]
        naConcInitEcs = self.rxdDict["ions"]["na"]["ecs"]

        kConcInitCyt = self.rxdDict["ions"]["k"]["cyt"]
        kConcInitEcs = self.rxdDict["ions"]["k"]["ecs"]

        caConcInitCyt = self.rxdDict["ions"]["ca"]["cyt"]
        caConcInitEcs = self.rxdDict["ions"]["ca"]["ecs"]
        caConcInitEr = self.rxdDict["ions"]["ca"]["er"]

        h.celsius = 37
        self.k = 1.38064852e-23  # Boltzmann constant
        self.T = h.celsius + 273.15  # Temperature in Kelvin
        # Parameters between this comment and the next comment comes from Hodgkin Huxley tutorial on NEURON website.
        self.e = 1.60217662e-19
        self.scale = 1e-14 / self.e  # convert S/cm2 to molecules/um2
        # Make sure to use CHANNEL DENSITIES. Below max conductance values
        # come from Targeted Simulation of RGC in Epiretinal Prostheses... Paknahad et al. (2020) (Soma of A2 cell type)
        self.gNaBar = 0.3 * self.scale  # molecules/um2 ms mV
        self.gKBar = 0.12 * self.scale  # molecules/um2 ms mV
        self.gKABar = 0 #3 * self.gKBar
        self.gHBar = 0.012 * self.scale
        self.gCaBarL = .0025 * self.scale  ## molecules/um2 ms mV
        self.gCaBarN = 0.0025 * self.scale  # molecules/um2 ms mV
        self.gCaBarT = 0.00025 * self.scale  # molecules/um2 ms mV
        self.gL = 0.0003 * self.scale  # molecules/um2 ms mV, assuming 30 um soma diameter and transmembrane resistance of 20000 Ohm*cm^2

        self.el = -70
        self.q10 = 3.0 ** ((h.celsius - 37.0) / 10.0)

        self.erFraction = self.rxdDict["region"]["erFraction"]
        self.cytFraction = 1 - self.erFraction

        # Rate Constants
        # PMCA
        kPMCABindF = 25000.0  # msec^-1 mM^-1
        kPMCABindR = 2.0 # msec^-1
        kPMCARel = 500.0  # msec^-1

        # NCX
        kNCXBindF = 93.827 #93.827 msec^-1 mM^-1
        kNCXBindR = 4.0 # msec^-1
        kNCXRel = 1000.0 # msec^-1 mM^-1

        # SERCA
        kSERCABindF = 17147.0
        kSERCABindR = 1.0
        kSERCARel = 250.0

        # Calreticulin
        KdCalr = 2.0
        kCalrBindF = 0.01
        kCalrBindR = KdCalr * kCalrBindF

        # Calbindin-D282k
        KdHigh = 237e-6
        kCbdHF = 11.0
        kCbdHR = KdHigh * kCbdHF
        KdLow = 411e-6
        kCbdLF = 87.0
        kCbdLR = KdLow * kCbdLF

        # Parvalbumin
        KdParv = 0.5315
        kParvBindF = 18.5
        kParvBindR = KdParv * kParvBindF

        # Initialize channel densities (steady-state) in respective kinetic states according to Hill-Langmuir equation
        # PMCA
        totalPMCA = self.rxdDict["channelDensities"]["totalPMCA"]  # Channel density of PMCA from Adam's paper. NEED TO calibrate
        pmcaCaInit = totalPMCA * caConcInitCyt / ((kPMCABindR / kPMCABindF) + caConcInitCyt)
        pmcaInit = totalPMCA - pmcaCaInit

        # NCX
        totalNCX = self.rxdDict["channelDensities"]["totalNCX"]
        ncx2CaInit = totalNCX * (caConcInitCyt**2) / ((kNCXBindR/kNCXBindF) + (caConcInitCyt**2))
        ncxInit = totalNCX - ncx2CaInit

        # SERCA
        totalSERCA = self.rxdDict["channelDensities"]["totalSERCA"]
        serca2CaInit = totalSERCA * (caConcInitCyt ** 2) / ((kSERCABindR / kSERCABindF) + (caConcInitCyt ** 2))
        sercaInit = totalSERCA - serca2CaInit

        # Calbindin-D28k
        totalCbd = self.rxdDict["endoBuffers"]["totalCbd"]
        # h0l0Init, h0l1Init, h0l2Init, h1l0Init, h1l1Init, h1l2Init, h2l0Init, h2l1Init, h2l2Init = cbdInit.get_init_conds_calbindin(totalCbd, caConcInitCyt)

        # Calreticulin
        totalCalr = self.rxdDict["endoBuffers"]["totalCalr"]
        calrCaInit = totalCalr * caConcInitCyt / ((kCalrBindR / kCalrBindF) + caConcInitCyt)
        calrInit = totalCalr - calrCaInit

        # Parvalbumin
        totalParv = self.rxdDict["endoBuffers"]["totalParv"]
        parvCaInit = totalParv * caConcInitCyt / ((kCalrBindR / kCalrBindF) + caConcInitCyt)
        parvInit = totalParv - parvCaInit


        ########################################
        # Hodgkin-Huxley equations for voltage-gated ion channels

        '''
        Note: "vtrap(x,y) is 1/(exp(x/y)-1) if |x/y|>=1e-6 or y*(1.0 - x/y/2.0) otherwise"
        https://neuron.yale.edu/neuron/docs/hodgkin-huxley-using-rxd
        vtrap is used to avoid singularities in equations that define transition probabilities.
        In some equations, denominator goes to zero, so vtrap offers a workaround.
        
        Below HH equations come from Targeted Simulation of RGC in Epiretinal Prostheses... Paknahad et al. (2020)
        '''
        # Na channel
        # m-gate
        self.alpha = 0.1 * vtrap(-(v + 40.0), 10.0)
        self.beta = 4.0 * exp(-(v + 65.0) / 18.0)
        self.mtau = 1.0 / (self.q10 * (self.alpha + self.beta))
        self.minf = self.alpha / (self.alpha + self.beta)
        # self.alpha = (-0.6 * (v + 30.0)) * vtrap(-(v + 30.0), 10.0)
        # self.beta = 20.0 * exp(-(v + 55.0) / 18.0)
        # self.mtau = 1.0 / (self.q10 * (self.alpha + self.beta))
        # self.minf = self.alpha / (self.alpha + self.beta)

        # h-gate
        self.alpha = 0.07 * exp(-(v + 65.0) / 20.0)
        self.beta = 1.0 / (exp(-(v + 35.0) / 10.0) + 1.0)
        self.htau = 1.0 / (self.q10 * (self.alpha + self.beta))
        self.hinf = self.alpha / (self.alpha + self.beta)
        # self.alpha = 0.4 * exp(-(v + 50.0) / 20.0)
        # self.beta = 6.0 / (exp(-0.1 * (v + 20.0) / 10.0) + 1.0)
        # self.htau = 1.0 / (self.q10 * (self.alpha + self.beta))
        # self.hinf = self.alpha / (self.alpha + self.beta)

        # K channel (Delayed-rectifier)
        # n-gate
        self.alpha = 0.01 * vtrap(-(v + 55.0), 10.0)
        self.beta = 0.125 * exp(-(v + 65.0) / 80.0)
        self.ntau = 1.0 / (self.q10 * (self.alpha + self.beta))
        self.ninf = self.alpha / (self.alpha + self.beta)
        # self.alpha = -0.02 * (v + 40.0) * vtrap(-(v + 40.0), 10.0)
        # self.beta = 0.4 * exp(-(v + 50.0) / 80.0)
        # self.ntau = 1.0 / (self.q10 * (self.alpha + self.beta))
        # self.ninf = self.alpha / (self.alpha + self.beta)

        # K,A channel
        # m-gate
        # self.alpha = -0.006 * (v + 90.0) * vtrap(-(v + 90.0), 10.0)
        # self.beta = 0.1 * exp(-(v + 30.0) / 10.0)
        # self.mtauA = 1.0 / (self.q10 * (self.alpha + self.beta))
        # self.minfA = self.alpha / (self.alpha + self.beta)
        #
        # # h-gate
        # self.alpha = 0.04 * exp(-(v + 70.0) / 20.0)
        # self.beta = 0.6 / (exp(-(v + 40.0) / 10.0) + 1)
        # self.htauA = 1.0 / (self.q10 * (self.alpha + self.beta))
        # self.hinfA = self.alpha / (self.alpha + self.beta)
        #
        # # H channel
        # # l-gate
        # self.alpha = exp(0.08316 * (v + 75.0))
        # self.beta = exp(0.033264 * (v + 75.0))
        # self.ltauH = 1.0 / (self.q10 * (self.alpha + self.beta))
        # self.linfH = self.alpha / (self.alpha + self.beta)
        '''
        The equations come from "VDCC documentation v1.3bis" Word Document file on AXON REGROWTH Google Drive folder.
        '''
        # Ca channel

        # T-channel
        # # m-gate
        # self.alpha = 0.2 * vtrap((-v + 19.26), 10)
        # self.beta = 9.0e-3 * exp(-v / 22.03)
        # self.minf_T = self.alpha / (self.alpha + self.beta)
        # self.tau_Tm = 1.0 / (self.q10 * (self.alpha + self.beta))
        # # h-gate
        # self.alpha = 1e-6 * exp(-v / 16.26)
        # self.beta = 1 / (exp((-v + 29.79) / 10) + 1)
        # self.hinf_T = self.alpha / (self.alpha + self.beta)
        # self.tau_Th = 1.0 / (self.q10 * (self.alpha + self.beta))

        # N-channel
        # m-gate
        self.alpha = 0.19 * vtrap((-v + 19.88), 10)
        self.beta = 4.6e-2 * exp(-v / 20.73)
        self.minf_N = self.alpha / (self.alpha + self.beta)
        self.tau_Nm = 1.0 / (self.q10 * (self.alpha + self.beta))

        # h-gate
        self.alpha = 1.6e-4 * exp(-v / 48.4)
        self.beta = 1 / (exp((-v + 39) / 10) + 1)
        self.hinf_N = self.alpha / (self.alpha + self.beta)
        self.tau_Nh = 1.0 / (self.q10 * (self.alpha + self.beta))

        # L-channel
        # m-gate
        self.alpha = 15.69 * vtrap((-v + 81.5), 10)
        self.beta = 0.29 * exp(-v / 10.86)
        self.minf_L = self.alpha / (self.alpha + self.beta)
        self.tau_Lm = 1.0 / (self.q10 * (self.alpha + self.beta))

        '''
        This code uses the complex equation of Ca current. A simple if-statement cannot compare crxd.v to a number due to how
        crxd.v is defined. A workaround uses y = (1 + tanh(x - threshold))/2, which transitions from 0 to 1 fairly quickly,
        allowing the function to be used as a binary switch that outputs 1 if our value is greater than threshold or 0 if value is
        less than threshold. As a result, temp3 stores either 1 - temp2/2 or temp2/exp(temp) - 1.
        '''
        # Formulas for Non-modulated Ca currents
        self.threshold = 1.10e-4
        self.m = 10000  # steepness of switch
        self.temp1 = 0.0853 * self.T / 2
        self.temp2 = v / self.temp1
        # switch goes to 0 for abs(temp2) < threshold and 1 for abs(temp2) > threshold
        self.switch = (1 + tanh(((fabs(self.temp2) - self.threshold) * 100))) / 2
        self.temp3 = (1 - self.switch) * (1 - self.temp2 / 2) + self.switch * (self.temp2 / exp(self.temp2) - 1)



        ########################################
        # RXD
        '''
        Code that uses RXD should define the following criteria: WHERE the reaction is happening, WHO (what species) are we interested
        in, and WHAT about those species of interest? With that being said, the following code
        answers, 'where?' The reaction is happening in: the cytoplasm, cell membrane, and extracellular fluid.
        '''
        ####### WHERE? #######
        self.cyt = rxd.Region(self.all, name='cyt', nrn_region='i', geometry = rxd.FractionalVolume(self.cytFraction))
        self.mem = rxd.Region(self.all, name='mem', geometry=rxd.membrane())
        # self.ecs = rxd.Region(self.all, name='ecs', nrn_region = 'o')
        self.ecs = rxd.Extracellular(self.xmin-self.padding, self.ymin-self.padding, self.zmin-self.padding,
                                     self.xmax+self.padding, self.ymax+self.padding, self.zmax+self.padding,
                                     dx=50)
        self.er = rxd.Region(self.all, name='er', geometry=rxd.FractionalVolume(self.erFraction))
        self.erCytMem = rxd.Region(self.all, name='erCytMem', geometry=rxd.ScalableBorder(0.1, on_cell_surface=False))

        '''
        This code aims to graph Ca current in order to study Ca signaling within a neuron. Hodgkin Huxley
        model will be able to track current of Na, K, and Ca ions flowing in/out of neuron.
        '''
        ####### WHO? #######
        # Ions
        self.k = rxd.Species([self.cyt, self.ecs], name="k", d=1, charge=1,
                             initial=self.concInit(kConcInitCyt, kConcInitEcs,0))
        self.na = rxd.Species([self.cyt, self.ecs], name="na", d=1, charge=1,
                              initial=self.concInit(naConcInitCyt, naConcInitEcs, 0))
        self.ca = rxd.Species([self.cyt, self.ecs, self.er], name='ca', d=1, charge=2,
                              initial=self.concInit(caConcInitCyt, caConcInitEcs, caConcInitEr))

        # Plasma Membrane Ca ATPase
        self.pmca = rxd.Species([self.cyt], name="pmca", initial=pmcaInit)
        self.pmcaCa = rxd.Species([self.cyt], name="pmcaCa", initial=pmcaCaInit) # PMCA bound to Ca


        # Na/Ca Exchange
        self.ncx = rxd.Species([self.cyt], name='ncx', initial=ncxInit)
        self.ncx2Ca = rxd.Species([self.cyt], name='ncx2Ca', initial=ncx2CaInit)

        # Sarcoplasmic/Endoplasmic Reticulum Ca ATPase
        self.serca = rxd.Species(self.cyt, name='serca', initial=sercaInit)
        self.serca2Ca = rxd.Species(self.cyt, name='serca2Ca', initial=serca2CaInit)

        # # Calbindin
        # self.cbdh0l0 = rxd.Species([self.cyt], name='cbdh0l0', initial=self.h0l0Init)
        # self.cbdh1l0 = rxd.Species([self.cyt], name='cbdh1l0', initial=self.h1l0Init)
        # self.cbdh2l0 = rxd.Species([self.cyt], name='cbdh2l0', initial=self.h2l0Init)
        # self.cbdh0l1 = rxd.Species([self.cyt], name='cbdh0l1', initial=self.h0l1Init)
        # self.cbdh1l1 = rxd.Species([self.cyt], name='cbdh1l1', initial=self.h1l1Init)
        # self.cbdh2l1 = rxd.Species([self.cyt], name='cbdh2l1', initial=self.h2l1Init)
        # self.cbdh0l2 = rxd.Species([self.cyt], name='cbdh0l2', initial=self.h0l2Init)
        # self.cbdh1l2 = rxd.Species([self.cyt], name='cbdh1l2', initial=self.h1l2Init)
        # self.cbdh2l2 = rxd.Species([self.cyt], name='cbdh2l2', initial=self.h2l2Init)

        # Calreticulin
        self.calr = rxd.Species(self.er, name='calr', initial=calrInit)
        self.calrCa = rxd.Species(self.er, name='calrCa', initial=calrCaInit)

        # Parvalbumin
        self.parv = rxd.Species(self.cyt, name='parv', initial=parvInit)
        self.parvCa = rxd.Species(self.cyt, name='parvCa', initial=parvCaInit)

        # Leak
        self.x = rxd.Species([self.cyt, self.ecs], name='x', charge=1)

        self.ki = self.k[self.cyt]
        self.ko = self.k[self.ecs]
        self.nai = self.na[self.cyt]
        self.nao = self.na[self.ecs]
        self.xi = self.x[self.cyt]
        self.xo = self.x[self.ecs]
        self.cai = self.ca[self.cyt]
        self.cao = self.ca[self.ecs]
        self.caer = self.ca[self.er]

        # To limit reactions to soma and axon
        inRegion = rxd.Parameter([self.cyt], name = 'inRegion',
                                 value = lambda nd: 1
                                 if (nd.segment in self.somaList[0] or nd.segment in self.axonList)
                                 else 0)

        # gating states -> open probabilities
        # K+ channel (delayed rectifier)
        self.ngate = rxd.State([self.cyt, self.mem], name='ngate', initial=0.2445865)

        # # K+ Channel (A channel)
        # self.mgateA = rxd.State([self.cyt, self.mem], name='mgateA', initial=0.0247887308)
        # self.hgateA = rxd.State([self.cyt, self.mem], name='hgateA', initial=0.543209973)
        #
        # # H channel (K+)
        # self.lgateH = rxd.State([self.cyt, self.mem], name='lgateH', initial=0.6416789483)
        #
        # Na+ Channel
        self.mgate = rxd.State([self.cyt, self.mem], name='mgate', initial=0.0289055)
        self.hgate = rxd.State([self.cyt, self.mem], name='hgate', initial=0.7540796)

        # Ca+2 Channel
        # mgate_T = rxd.State([cyt, mem], name='mgate_T', initial=0.010871215)
        # hgate_T = rxd.State([cyt, mem], name='hgate_T', initial=0.615047335)
        self.mgate_N = rxd.State([self.cyt, self.mem], name='mgate_N', initial=0.001581552)
        self.hgate_N = rxd.State([self.cyt, self.mem], name='hgate_N', initial=0.973556944)
        self.mgate_L = rxd.State([self.cyt, self.mem], name='mgate_L', initial=3.42574e-6)


        """
        NOTE: 'mgate' and 'm_gate' are different. 'mgate' is an object of class rxd.State, representing
        open probability value. We are telling NEURON 'who' we are interested in, which in this case are the
        various types of gates.
        'm_gate'is an rxd.Rate object which we are telling NEURON to track over time according
        to the equation (minf-mgate)/mtau. Using the rxd.Rate class allows the rxd.State objects
        to change over time.
        """

        ####### WHAT? #######
        # Na channel
        self.m_gate = rxd.Rate(self.mgate, (self.minf - self.mgate) / self.mtau)
        self.h_gate = rxd.Rate(self.hgate, (self.hinf - self.hgate) / self.htau)

        # K Channel
        # Delayed rectifier
        self.n_gate = rxd.Rate(self.ngate, (self.ninf - self.ngate) / self.ntau)

        # # A channel
        # self.m_gateA = rxd.Rate(self.mgateA, (self.minfA - self.mgateA) / self.mtauA)
        # self.h_gateA = rxd.Rate(self.hgateA, (self.hinfA - self.hgateA) / self.htauA)

        # H channel
        # self.l_gateH = rxd.Rate(self.lgateH, (self.linfH - self.lgateH) / self.ltauH)

        # Ca Channel
        # m_gate_T = rxd.Rate(mgate_T, (minf_T - mgate_T) / tau_Tm)
        # h_gate_T = rxd.Rate(hgate_T, (hinf_T - hgate_T) / tau_Th)
        self.m_gate_N = rxd.Rate(self.mgate_N, (self.minf_N - self.mgate_N) / self.tau_Nm)
        self.h_gate_N = rxd.Rate(self.hgate_N, (self.hinf_N - self.hgate_N) / self.tau_Nh)
        self.m_gate_L = rxd.Rate(self.mgate_L, (self.minf_L - self.mgate_L) / self.tau_Lm)

        # Nernst Potentials
        self.Ena = 1e3 * h.R * self.T * log(self.nao / self.nai) / h.FARADAY
        self.Ek = 1e3 * h.R * self.T * log(self.ko / self.ki) / h.FARADAY
        self.Eca = 1e3 * h.R * self.T * log(self.cao / self.cai) / (2 * h.FARADAY)
        self.Eh = 0 # mV


        # Conductances
        self.gNa = inRegion * self.gNaBar * self.mgate ** 3 * self.hgate
        self.gK = inRegion * self.gKBar * self.ngate ** 4
        # self.gCaT = inRegion * self.gcabar_T * self.mgate_T ** 2 * self.hgate_T
        self.gCaN = inRegion * self.gCaBarN * self.mgate_N ** 2 * self.hgate_N
        self.gCaL = inRegion * self.gCaBarL * self.mgate_L ** 2
        # self.gKA = inRegion * self.gKABar * self.mgateA ** 3 * self.hgateA

        # self.gna = self.gnabar * self.mgate ** 3 * self.hgate
        # self.gk = self.gkbar * self.ngate ** 4
        # # gca_T = gcabar_T * mgate_T ** 2 * hgate_T
        # self.gca_N = self.gcabar_N * self.mgate_N ** 2 * self.hgate_N
        # self.gca_L = self.gcabar_L * self.mgate_L ** 2

        self.dvf = 0.001 / (0.001 + self.cai) * self.temp1 * self.temp3 * (1 - self.cai / self.cao * exp(self.temp2))
        self.I_non_mod_L = self.gCaL * self.dvf
        self.I_non_mod_N = self.gCaN * self.dvf

        # Ca Extrusion by PMCA
        self.pmcaBinding = rxd.Reaction(self.ca + self.pmca, self.pmcaCa, kPMCABindF, kPMCABindR,
                                         region=[self.cyt])
        self.pmcaRelease = rxd.Reaction(self.pmcaCa, self.pmca, kPMCARel, region=[self.cyt])

        # Na/Ca Exchange
        self.ncxBinding = rxd.Reaction(2*self.ca + self.ncx, self.ncx2Ca, kNCXBindF, kNCXBindR, region=[self.cyt])
        self.ncxRelease = rxd.Reaction(self.ncx2Ca, self.ncx, kNCXRel, region=[self.cyt])

        # SERCA
        self.sercaBinding = rxd.Reaction(2*self.ca + self.serca, self.serca2Ca, kSERCABindF, kSERCABindR, region=[self.cyt])
        self.sercaRelease = rxd.MultiCompartmentReaction(self.serca2Ca[self.cyt], 2*self.caer + self.serca[self.cyt], kSERCARel, membrane=self.erCytMem, membrane_flux=True)

        # Calbindin-D28k
        # self.h0l0HBind = rxd.Reaction(self.ca + self.cbdh0l0, self.cbdh1l0)
        # self.h0l0LBind
        # self.h0l1HBind
        # self.h0l1LBind
        # self.h0l2HBind
        # self.h1l0LBind
        # self.h1l0HBind
        # self.h1l1LBind
        # self.h1l1HBind
        # self.h1l2HBind
        # self.h2l0LBind
        # self.h2l1LBind

        # Calreticulin
        self.calrBuffering = rxd.Reaction(self.caer + self.calr, self.calrCa, kCalrBindF, kCalrBindR, region=self.er)

        # Parvalbumin
        self.parvBuffering = rxd.Reaction(self.cai + self.parv, self.parvCa, kParvBindF, kParvBindR, region=self.cyt)
        '''
        rxd.MultiCompartmentReaction tracks current across multiple areas, in this case, cytoplasm and extracellular fluid
        '''
        # Current through ion channels
        self.INa = rxd.MultiCompartmentReaction(self.nai, self.nao, self.gNa * (v - self.Ena), mass_action=False,
                                                membrane=self.mem, membrane_flux=True)

        totalIK = (self.gK + self.gKABar) * (v - self.Ek)
        self.IK = rxd.MultiCompartmentReaction(self.ki, self.ko, totalIK, mass_action=False,
                                               membrane=self.mem, membrane_flux=True)

        # ICaT = rxd.MultiCompartmentReaction(cai, cao, I_non_mod_T, mass_action=False,
        # membrane=mem, membrane_flux=True)

        self.ICaN = rxd.MultiCompartmentReaction(self.cai, self.cao, self.I_non_mod_N, mass_action=False,
                                                 membrane=self.mem, membrane_flux=True)

        self.ICaL = rxd.MultiCompartmentReaction(self.cai, self.cao, self.I_non_mod_L, mass_action=False,
                                                 membrane=self.mem, membrane_flux=True)

        self.ILeak = rxd.MultiCompartmentReaction(self.xi, self.xo, self.gL * (v - self.el), mass_action=False,
                                                  membrane=self.mem, membrane_flux=True)



    def stimulate(self, stimNeuron, measureNeuron, location, delay, amp, stim_dur, sim_dur):
        '''


        :param neuron: Section of neuron to place electrode
        :param location: Location within section
        :param delay: No current from start of simulation
        :param amp: Current amplitude in nA
        :param stim_dur: Duration of current injection
        :param sim_dur: Duration of simulation
        :return: NONE
        '''

        self.stim = h.IClamp(stimNeuron(location))
        self.stim.delay = delay

        self.stim.amp = amp
        self.stim.dur = stim_dur

        '''
        Following vectors are necessary in order to create graphs.
        '''
        self.tvec = h.Vector().record(h._ref_t)
        self.vvec = h.Vector().record(measureNeuron(location)._ref_v)

        # self.mvec = h.Vector().record(self.mgate[self.cyt].nodes(measureNeuron(location))._ref_value)
        # self.nvec = h.Vector().record(self.ngate[self.cyt].nodes(measureNeuron(location))._ref_value)
        # self.hvec = h.Vector().record(self.hgate[self.cyt].nodes(measureNeuron(location))._ref_value)
        self.caL_mvec = h.Vector().record(self.mgate_L[self.cyt].nodes(measureNeuron(location))._ref_value)
        self.caN_mvec = h.Vector().record(self.mgate_N[self.cyt].nodes(measureNeuron(location))._ref_value)
        self.caN_hvec = h.Vector().record(self.hgate_N[self.cyt].nodes(measureNeuron(location))._ref_value)
        #
        # self.kivec = h.Vector().record(measureNeuron(location)._ref_ki)
        # self.kovec = h.Vector().record(measureNeuron(location)._ref_ko)
        # self.naovec = h.Vector().record(measureNeuron(location)._ref_nao)
        # self.naivec = h.Vector().record(measureNeuron(location)._ref_nai)
        self.caovec = h.Vector().record(measureNeuron(location)._ref_cao)
        self.caivec = h.Vector().record(measureNeuron(location)._ref_cai)
        self.caervec = h.Vector().record(self.caer.nodes(measureNeuron(location))._ref_concentration)

        # self.icavec = h.Vector().record(measureNeuron(location)._ref_ica)

        # self.pmcavec = h.Vector().record(self.pmca.nodes(measureNeuron(location))._ref_value)
        # self.pmcaCavec = h.Vector().record(self.pmcaCa.nodes(measureNeuron(location))._ref_value)
        #
        # self.ncxvec = h.Vector().record(self.ncx.nodes(measureNeuron(location))._ref_value)
        # self.ncx2cavec = h.Vector().record(self.ncx2Ca.nodes(measureNeuron(location))._ref_value)

        self.calrvec = h.Vector().record(self.calr.nodes(measureNeuron(location))._ref_concentration)
        self.calrcavec = h.Vector().record(self.calrCa.nodes(measureNeuron(location))._ref_concentration)

        # self.sercavec = h.Vector().record(self.serca[self.cyt].nodes(measureNeuron(location))._ref_concentration)
        # self.serca2cavec = h.Vector().record(self.serca2Ca[self.cyt].nodes(measureNeuron(location))._ref_concentration)

        # h.cvode.active(1)
        #
        # h.v_init = -70.0
        #
        # h.tstop = sim_dur
        # h.celsius = 34.0

        # This line needed:
        # h.finitialize(-70 * mV)
        fih = h.FInitializeHandler(callable)
        fih.allprint()

        print(h.hname())
        h.continuerun(sim_dur)




    def interp(self, vec):
        return vec.as_numpy() - 0.5 * np.diff(vec.as_numpy(), prepend=vec.x[0])

    def flux(self, vec):
        return np.diff(vec.as_numpy(), prepend=vec.x[0]) / h.dt

    def rxd_flux_mols(self, reaction):
        self.rate_str = repr(-reaction.f_rate) if reaction.b_rate is None else repr(reaction.b_rate - reaction.f_rate)

        self.rxd_to_py = {'v': 'self.interp(self.vvec)',
                          'log': 'np.log',
                          'exp': 'np.exp',
                          'tanh': 'np.tanh',
                          'fabs': 'np.fabs',
                          'self.k[self.cyt]': 'self.interp(self.kivec)',
                          'self.k[self.ecs]': 'self.interp(self.kovec)',
                          'self.na[self.cyt]': 'self.interp(self.naivec)',
                          'self.na[self.ecs]': 'self.interp(self.naovec)',
                          'self.ngate': 'self.interp(self.nvec)',
                          'self.mgate': 'self.interp(self.mvec)',
                          'self.hgate': 'self.interp(self.hvec)',
                          'self.mgate_L': 'self.interp(self.ca_Lvec)',
                          'self.mgate_N': 'self.interp(self.caN_mvec)',
                          'self.hgate_N': 'self.interp(self.caN_hvec)',
                          'self.ca[self.cyt]': 'self.interp(self.caivec)',
                          'self.ca[self.ecs]': 'self.interp(self.caovec)'}

        self.current_rate = functools.reduce(lambda r, kv: r.replace(*kv), self.rxd_to_py.items(), self.rate_str)

        self.flux_mols = eval(self.current_rate)

        return self.flux_mols

    def rxd_flux_mM(self, reaction, node):
        self.flux_mols = self.rxd_flux_mols(reaction)

        return self.flux_mols * node.surface_area / (node.volume * rxd.constants.molecules_per_mM_um3())

if __name__ == "__main__":
    rgc1 = rgcMorph.RGCMorphology()
    vgcc = VGCC_2(rgc1, 'rxd_file.json')
    vgcc.implement_RXD(rgc1)
    vgcc.stimulate(vgcc.somaList[0], vgcc.somaList[0], 0.5, 0, 0.5, 60, 100)
    # print(vgcc.somaList[0].diam)


    fig = pyplot.figure()
    pyplot.plot(vgcc.tvec.as_numpy(), vgcc.vvec.as_numpy(), '-b', label='k')
    pyplot.xlabel('t (ms)')
    pyplot.ylabel('Voltage (mV)')
    pyplot.savefig(os.path.join('MemPotTest.svg'), format="svg")


    # fig = pyplot.figure()
    # pyplot.plot(vgcc.tvec.as_numpy(), vgcc.caL_mvec.as_numpy(), '-b', label='L m-gate')
    # # pyplot.plot(vgcc.tvec.as_numpy(), vgcc.caN_hvec.as_numpy(), '-r', label='N h-gate')
    # pyplot.plot(vgcc.tvec.as_numpy(), vgcc.caN_mvec.as_numpy(), '-g', label='N m-gate')
    # pyplot.legend()
    # pyplot.xlabel('t (ms)')
    # pyplot.ylabel('Open Probability')
    # pyplot.savefig(os.path.join('OpenProbTest.svg'), format="svg")
    #
    # fig = pyplot.figure()
    # pyplot.plot(vgcc.tvec, vgcc.icavec.as_numpy(), '-g', label='k')
    # pyplot.xlabel('t (ms)')
    # pyplot.ylabel('Current (mA/cm^2)')
    # pyplot.savefig(os.path.join('CaCurrentClass.svg'), format="svg")

    fig = pyplot.figure()
    pyplot.plot(vgcc.tvec.as_numpy(), vgcc.caivec.as_numpy(), '-g', label='k')
    pyplot.xlabel('time (ms)')
    pyplot.ylabel('Intracellular Ca Concentration (mM)')
    pyplot.savefig(os.path.join('IntraCaTest.svg'), format="svg")

    fig = pyplot.figure()
    pyplot.plot(vgcc.tvec.as_numpy(), vgcc.caervec.as_numpy(), '-g', label='k')
    pyplot.xlabel('time (ms)')
    pyplot.ylabel('ER Lumenal Ca Concentration (mM)')
    pyplot.savefig(os.path.join('LumenCaTest.svg'), format="svg")

    # fig = pyplot.figure()
    # pyplot.plot(vgcc.tvec, vgcc.pmcavec.as_numpy(), '-b', label='PMCA')
    # pyplot.plot(vgcc.tvec, vgcc.pmcaCavec.as_numpy(), '-r', label='PMCA_CA')
    # pyplot.xlabel('t (ms)')
    # pyplot.ylabel('Concentration')
    # pyplot.legend()
    # pyplot.savefig(os.path.join('PMCAClass.svg'), format="svg")

    # fig = pyplot.figure()
    # pyplot.plot(vgcc.tvec, vgcc.ncxvec.as_numpy(), '-b', label='ncx')
    # pyplot.plot(vgcc.tvec, vgcc.ncx2cavec.as_numpy(), '-r', label='ncx_2ca')
    # pyplot.xlabel('t (ms)')
    # pyplot.ylabel('Concentration')
    # pyplot.legend()
    # pyplot.savefig(os.path.join('NCXClass.svg'), format="svg")

    # fig = pyplot.figure()
    # pyplot.plot(vgcc.tvec, vgcc.sercavec.as_numpy(), '-b', label='serca')
    # pyplot.plot(vgcc.tvec, vgcc.serca2cavec.as_numpy(), '-r', label='serca_2ca')
    # pyplot.xlabel('t (ms)')
    # pyplot.ylabel('Concentration')
    # pyplot.legend()
    # pyplot.savefig(os.path.join('SERCA.svg'), format="svg")

    # fig = pyplot.figure()
    # pyplot.plot(vgcc.tvec, vgcc.calrvec.as_numpy(), '-b', label='calr')
    # pyplot.plot(vgcc.tvec, vgcc.calrcavec.as_numpy(), '-r', label='calr_ca')
    # pyplot.xlabel('t (ms)')
    # pyplot.ylabel('Concentration')
    # pyplot.legend()
    # pyplot.savefig(os.path.join('calr.svg'), format="svg")

    # print(vgcc.caivec[-1])


    rgc1 = vgcc = None

RGCMorphology.py:

Code: Select all

from neuron import h

h.load_file('stdrun.hoc')
h.load_file("import3d.hoc")
h.load_file("stdlib.hoc")

class RGCMorphology:
    def __init__(self):
        self.load_morphology()


    def load_morphology(self):
        cell = h.Import3d_SWC_read()
        cell.input("8.swc")
        i3d = h.Import3d_GUI(cell, False)
        i3d.instantiate(self)
If you have any suggestions, I will try anything. Thank you!
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: h.finitialize() fails with a hocobj call error

Post by hines »

Running the model requires the 8.swc file. Please send that to michael.hines@yale.edu as a an attachment. Or zip file if there are are other files needed to run the model.
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: h.finitialize() fails with a hocobj call error

Post by hines »

Thank you for sending me the zip file with the code. In getting back into context, I did notice above that
Could not find module 'C:\Users\jdami\Documents\Year 4\Lazzi\...
Could you rerun after modifying the location of your files so that you eliminate the space in "Year 4". Changing to "Year4" is sufficient. Let me know if that fixes your issue.
neuron121!
Posts: 8
Joined: Thu Mar 31, 2022 11:00 pm

Re: h.finitialize() fails with a hocobj call error

Post by neuron121! »

Hello. After removing the space, the same error was achieved.
neuron121!
Posts: 8
Joined: Thu Mar 31, 2022 11:00 pm

Re: h.finitialize() fails with a hocobj call error

Post by neuron121! »

Is there anything else I can try? I'm still stuck on this issue
ramcdougal
Posts: 267
Joined: Fri Nov 28, 2008 3:38 pm
Location: Yale School of Public Health

Re: h.finitialize() fails with a hocobj call error

Post by ramcdougal »

Here's a possible clue... it's looking for a .so file. Usually on Windows, the equivalent is a .dll

Are you using the regular Windows version of NEURON or a version the windows subsystem for linux (WSL)?

Either way, if you try running with python -i myprogram.py (that is, add the -i flag so it doesn't just quit), is it generating any files in that directory that look anything like the file it complains about not being able to load? Like a .c file or a .dll?
Post Reply