rxd.Reaction(...) causes finitialize to not work

Extending NEURON to handle reaction-diffusion problems.

Moderators: hines, wwlytton, ramcdougal

Post Reply
neuron121!
Posts: 8
Joined: Thu Mar 31, 2022 11:00 pm

rxd.Reaction(...) causes finitialize to not work

Post by neuron121! »

Hello all,

I have been debugging some code sent to me from a MacOS for a while (I use Windows). Basically, I've traced the error to every time I call an instance of self.'anything' = rxd.Reaction('anything'), and I get the following error when I run f.initialize(-65):

Code: Select all

C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s: Assembler messages:
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:7: Error: bad register name `%rbp'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:8: Warning: .seh_pushreg ignored for this target
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:9: Error: bad register name `%rsp'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:10: Warning: .seh_setframe ignored for this target
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:11: Error: bad register name `%rsp'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:12: Warning: .seh_stackalloc ignored for this target
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:14: Error: bad register name `%rcx'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:15: Error: bad register name `%rdx'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:16: Error: bad register name `%r8'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:17: Error: bad register name `%r9'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:18: Error: bad register name `%rbp)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:19: Error: bad register name `%rax)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:20: Error: bad register name `%rax)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:21: Error: bad register name `%rip)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:23: Error: bad register name `%rbp)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:24: Error: bad register name `%rax'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:25: Error: bad register name `%rax)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:26: Error: bad register name `%rax)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:28: Error: bad register name `%rbp)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:29: Error: bad register name `%rax'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:30: Error: bad register name `%rax)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:31: Error: bad register name `%rax)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:36: Error: bad register name `%rbp)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:37: Error: bad register name `%rbp)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:38: Error: bad register name `%rax)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:39: Error: bad register name `%rbp)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:40: Error: bad register name `%rip)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:42: Error: bad register name `%rax)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:43: Error: bad register name `%rbp)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:44: Error: bad register name `%rax'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:45: Error: bad register name `%rax)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:46: Error: bad register name `%rbp)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:47: Error: bad register name `%rip)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:49: Error: bad register name `%rax)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:50: Error: bad register name `%rbp)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:51: Error: bad register name `%rax'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:52: Error: bad register name `%rax)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:53: Error: bad register name `%rbp)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:54: Error: bad register name `%rax)'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:56: Error: bad register name `%rsp'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:57: Error: bad register name `%rbp'
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:59: Warning: zero assumed for missing expression
C:\Users\jdami\AppData\Local\Temp\ccKC1hT6.s:59: 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__
    self._handle = _dlopen(self._name, mode)
FileNotFoundError: Could not find module 'C:\Users\jdami\Documents\Year4\Lazzi\Colab\ITEMS_AxonGrowth\rxddll1269e012-0c73-11ed-9c11-7470fdcdedcf.so' (or one of its dependencies). Try using the full path with constructor syntax.

NEURON: Python Callback failed
 near line 0
 ^
        finitialize(-60)
Traceback (most recent call last):
  File "c:\Users\jdami\Documents\Year4\Lazzi\Colab\ITEMS_AxonGrowth\test3_recreation.py", line 893, in <module>
    vgcc.stimulate(vgcc.somaList[0], vgcc.somaList[0], 0.5, 0, 0.4, 50, 100)
  File "c:\Users\jdami\Documents\Year4\Lazzi\Colab\ITEMS_AxonGrowth\test3_recreation.py", line 838, in stimulate
    h.finitialize(-60)
RuntimeError: hocobj_call error


It also creates the following file called rxddll1269e012-0c73-11ed-9c11-7470fdcdedcf.c:

Code: Select all

#include <math.h>
/*Some functions supported by numpy that aren't included in math.h
 * names and arguments match the wrappers used in rxdmath.py
 */
double factorial(const double);
double degrees(const double);
void radians(const double, double*);
double log1p(const double);
double vtrap(const double, const double);
void reaction(double** species, double** params, double** rhs, double* mult, double* species_3d, double* params_3d, double* rhs_3d, double** flux, double v)
{
	double rate;
	rate = ((25000.0)*(species[0][0]))*(species[1][0])-((2.0)*(species[2][0]));
	rhs[0][0] = (-1) * rate;
	rhs[1][0] = (-1) * rate;
	rhs[2][0] = (1) * rate;
	rate = (500.0)*(species[2][0])-(0);
	rhs[1][0] += (1) * rate;
	rhs[2][0] += (-1) * rate;
	rate = ((93.827)*(pow(species[0][0], 2)))*(species[3][0])-((4.0)*(species[4][0]));
	rhs[0][0] += (-2) * rate;
	rhs[3][0] = (-1) * rate;
	rhs[4][0] = (1) * rate;
	rate = (1000.0)*(species[4][0])-(0);
	rhs[3][0] += (1) * rate;
	rhs[4][0] += (-1) * rate;
	rate = ((11.0)*(species[0][0]))*(species[5][0])-((0.002607)*(species[6][0]));
	rhs[0][0] += (-1) * rate;
	rhs[5][0] = (-1) * rate;
	rhs[6][0] = (1) * rate;
	rate = ((8.7)*(species[0][0]))*(species[5][0])-((0.0035757)*(species[8][0]));
	rhs[0][0] += (-1) * rate;
	rhs[5][0] += (-1) * rate;
	rhs[8][0] = (1) * rate;
	rate = ((11.0)*(species[0][0]))*(species[8][0])-((0.002607)*(species[9][0]));
	rhs[0][0] += (-1) * rate;
	rhs[8][0] += (-1) * rate;
	rhs[9][0] = (1) * rate;
	rate = ((8.7)*(species[0][0]))*(species[8][0])-((0.0035757)*(species[11][0]));
	rhs[0][0] += (-1) * rate;
	rhs[8][0] += (-1) * rate;
	rhs[11][0] = (1) * rate;
	rate = ((11.0)*(species[0][0]))*(species[11][0])-((0.002607)*(species[12][0]));
	rhs[0][0] += (-1) * rate;
	rhs[11][0] += (-1) * rate;
	rhs[12][0] = (1) * rate;
	rate = ((8.7)*(species[0][0]))*(species[6][0])-((0.0035757)*(species[9][0]));
	rhs[0][0] += (-1) * rate;
	rhs[6][0] += (-1) * rate;
	rhs[9][0] += (1) * rate;
	rate = ((11.0)*(species[0][0]))*(species[6][0])-((0.002607)*(species[7][0]));
	rhs[0][0] += (-1) * rate;
	rhs[6][0] += (-1) * rate;
	rhs[7][0] = (1) * rate;
	rate = ((8.7)*(species[0][0]))*(species[9][0])-((0.0035757)*(species[12][0]));
	rhs[0][0] += (-1) * rate;
	rhs[9][0] += (-1) * rate;
	rhs[12][0] += (1) * rate;
	rate = ((11.0)*(species[0][0]))*(species[9][0])-((0.002607)*(species[10][0]));
	rhs[0][0] += (-1) * rate;
	rhs[9][0] += (-1) * rate;
	rhs[10][0] = (1) * rate;
	rate = ((11.0)*(species[0][0]))*(species[12][0])-((0.002607)*(species[12][0]));
	rhs[0][0] += (-1) * rate;
	rhs[12][0] += (0) * rate;
	rate = ((8.7)*(species[0][0]))*(species[7][0])-((0.0035757)*(species[10][0]));
	rhs[0][0] += (-1) * rate;
	rhs[7][0] += (-1) * rate;
	rhs[10][0] += (1) * rate;
	rate = ((8.7)*(species[0][0]))*(species[10][0])-((0.0035757)*(species[13][0]));
	rhs[0][0] += (-1) * rate;
	rhs[10][0] += (-1) * rate;
	rhs[13][0] = (1) * rate;
	rate = ((0.01)*(species[0][1]))*(species[14][1])-((0.02)*(species[15][1]));
	rhs[0][1] = (-1) * rate;
	rhs[14][1] = (-1) * rate;
	rhs[15][1] = (1) * rate;
	rate = ((18.5)*(species[0][0]))*(species[16][0])-((9.832749999999999)*(species[17][0]));
	rhs[0][0] += (-1) * rate;
	rhs[16][0] = (-1) * rate;
	rhs[17][0] = (1) * rate;
	rhs[39][2] = (((((tanh(20000*(species[0][0])-20000*(0.0002))+1)/(2))*(species[0][0]))/(0.0002+species[0][0]))*(species[0][1]-(species[0][0]))-(species[39][2]))/(1.2);
	rate = ((params[1][2])*(species[39][2]))*(species[0][1]);
	rhs[0][1] += mult[0] * rate;
	rhs[0][0] += mult[1] * rate;
	rhs[37][2] = ((pow(species[0][0], 2))/(pow(species[0][0], 2)+1e-06)-(species[37][2]))/(1.0);
	rate = ((params[0][2])*(species[37][2]))*(species[0][0]);
	rhs[0][0] += mult[2] * rate;
	if(flux) flux[0][0] = 2.0 * rate;
	rhs[0][1] += mult[3] * rate;
	if(flux) flux[0][1] = 0.0 * rate;
	rate = ((species[18][0])*(2.78))*(species[19][0])-((0.00215)*(species[22][0]));
	rhs[19][0] = (-1) * rate;
	rhs[22][0] = (1) * rate;
	rate = ((species[18][0])*(2.78))*(species[23][0])-((2.15e-05)*(species[25][0]));
	rhs[23][0] = (-1) * rate;
	rhs[25][0] = (1) * rate;
	rate = ((species[18][0])*(2.78))*(species[24][0])-((2.15e-05)*(species[26][0]));
	rhs[24][0] = (-1) * rate;
	rhs[26][0] = (1) * rate;
	rate = ((2.7e-07)*(species[20][0]))*(species[19][0])-((0.0068)*(species[23][0]));
	rhs[19][0] += (-1) * rate;
	rhs[20][0] = (-1) * rate;
	rhs[23][0] += (1) * rate;
	rate = ((2.7e-06)*(species[20][0]))*(species[22][0])-((0.00068)*(species[25][0]));
	rhs[20][0] += (-1) * rate;
	rhs[22][0] += (-1) * rate;
	rhs[25][0] += (1) * rate;
	rate = ((2.7e-07)*(species[21][0]))*(species[19][0])-((0.0068)*(species[24][0]));
	rhs[19][0] += (-1) * rate;
	rhs[21][0] = (-1) * rate;
	rhs[24][0] += (1) * rate;
	rate = ((2.7e-06)*(species[21][0]))*(species[22][0])-((0.00068)*(species[26][0]));
	rhs[21][0] += (-1) * rate;
	rhs[22][0] += (-1) * rate;
	rhs[26][0] += (1) * rate;
	rate = (1.5e-08)*(species[20][0])-(0);
	rhs[20][0] += (-1) * rate;
	rhs[21][0] += (1) * rate;
	rhs[28][0] = (1) * rate;
	rate = (1.5e-08)*(species[23][0])-(0);
	rhs[23][0] += (-1) * rate;
	rhs[24][0] += (1) * rate;
	rhs[28][0] += (1) * rate;
	rate = (0.00065)*(species[25][0])-(0);
	rhs[25][0] += (-1) * rate;
	rhs[26][0] += (1) * rate;
	rhs[28][0] += (1) * rate;
	rate = (0.0047)*(species[30][0])-(0);
	rhs[30][0] = (-1) * rate;
	rhs[31][0] = (1) * rate;
	rate = (2.6e-05)*(species[28][0])-(0);
	rhs[27][0] = (1) * rate;
	rhs[28][0] += (-1) * rate;
	rate = (0.015)*(species[31][0])-(0);
	rhs[30][0] += (1) * rate;
	rhs[31][0] += (-1) * rate;
	rate = ((0.01)*(species[29][0]))*(species[28][0])-(0);
	rhs[28][0] += (-1) * rate;
	rhs[29][0] = (-1) * rate;
	rhs[31][0] += (1) * rate;
	rate = (0.0071)*(species[30][0])-(0);
	rhs[27][0] += (1) * rate;
	rhs[29][0] += (1) * rate;
	rhs[30][0] += (-1) * rate;
	rate = ((0.001)*(species[21][0]))*(species[27][0])-(0);
	rhs[20][0] += (1) * rate;
	rhs[21][0] += (-1) * rate;
	rhs[27][0] += (-1) * rate;
	rhs[33][0] = (species[32][0])*(1.992e-05);
	rhs[33][0] += (species[33][0])*(-0.0009960000000000001);
	rate = (0.00016600000000000002)*(species[33][0])-((0.0002324)*(species[34][0]));
	rhs[33][0] += (-1) * rate;
	rhs[34][0] = (1) * rate;
	rate = ((species[31][0])*(0.09))*(species[34][0])-(0);
	rhs[34][0] += (-1) * rate;
	rhs[35][0] = (1) * rate;
}
The name of the file seems random, and I'm not sure why it makes it or how to fix it.


I've included the source code for reference. The rxd.Reaction codes are around lines 580-740:

Code: Select all

from xml.etree.ElementInclude import include
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 os
#include <extra_numpy.h>
# from ctypes import *
# import extra_numpy.dll
#load the shared object file
# CDLL('extra_numpy.dll')
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)

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.apic]
        # self.basalList = [sec for sec in neuronMorph.basal]
        # self.somaList = [sec for sec in neuronMorph.soma]


        self.all = self.axonList + self.somaList #+ self.basalList #+ self.apicalList


        self.findExtrema()
        self.padding = 100

        # 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.create_neuron()
        # self.all = self.soma



    def createSoma(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

    def loadMorph(self):
        cell = h.Import3d_SWC_read()
        cell.input("8.swc")
        i3d = h.Import3d_GUI(cell, False)
        i3d.instantiate(self)

        return cell



    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 insertChannels(self):
        # HH Model by Fohlmeister
        for sec in self.axonList:
            # sec.insert('cad')
            # sec.depth_cad = 3 # (micron)
            # sec.taur_cad = 10 # (ms)

            sec.insert('pas')
            sec.e_pas = -60
            sec.g_pas = 0.05e-3

            sec.insert('spike')
            sec.gnabar_spike = 0.2
            sec.gkbar_spike = 0.211
            sec.gabar_spike = 3 * sec.gkbar_spike
            sec.gcabar_spike = 0.013
            sec.gkcbar_spike = 0.004 * sec.gkbar_spike

        for sec in self.somaList:
            # sec.insert('cad')
            # sec.depth_cad = 3  # (micron)
            # sec.taur_cad = 10  # (ms)

            sec.insert('pas')
            sec.e_pas = -60
            sec.g_pas = 0.05e-3
            sec.Ra = 110

            sec.insert('spike')
            sec.gnabar_spike = 0.2
            sec.gkbar_spike = 0.211
            sec.gabar_spike = 3 * sec.gkbar_spike
            sec.gcabar_spike = 0.001
            sec.gkcbar_spike = 0.004 * sec.gkbar_spike

        # for sec in self.apicalList:
        #
        #     sec.insert('cad')
        #     sec.depth_cad = 3  # (micron)
        #     sec.taur_cad = 10  # (ms)
        #
        #     sec.insert('pas')
        #     sec.e_pas = -60
        #     sec.g_pas = 0.05e-3
        #
        #     sec.insert('spike')
        #     sec.gnabar_spike = 0.08
        #     sec.gkbar_spike = 0.08
        #     sec.gabar_spike = 3 * sec.gkbar_spike
        #     sec.gcabar_spike = 0.01
        #     sec.gkcbar_spike = 0.004 * sec.gkbar_spike
        #
        # for sec in self.basalList:
        #     sec.insert('cad')
        #     sec.depth_cad = 3  # (micron)
        #     sec.taur_cad = 10  # (ms)
        #
        #     sec.insert('pas')
        #     sec.e_pas = -60
        #     sec.g_pas = 0.05e-3
        #
        #     sec.insert('spike')
        #     sec.gnabar_spike = 0.08
        #     sec.gkbar_spike = 0.08
        #     sec.gabar_spike = 3 * sec.gkbar_spike
        #     sec.gcabar_spike = 0.01
        #     sec.gkcbar_spike = 0.004 * sec.gkbar_spike

        # self.all.insert('spike')
        # self.all.insert('cad')
        # self.all.insert('pas')
        # self.all.e_pas = -60
        # self.all.g_pas = 0.05e-3

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

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

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

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

        IP3ConcInit = self.rxdDict["species"]['IP3']
        IP3Diffusion = 0.3

        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 = 8.7
        kCbdLR = KdLow * kCbdLF

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

        # Calcium-induced Ca release
        self.VCICR = 100000.0  # cm^-2 ms^-1
        self.KCICR = 0.0002  # mM
        tauCICR = 1.2  # ms
        KtCICR = 0.0002  # mM

        ## Code and kinetic models regarding AChR and Ip3-gated Ca release come from Mergenthal 2020.
        # Muscarinic AChR
        kL1F = 2.78  # (mM^-1 msec^-1) 0.00278
        kL1R = 2.15e-3  # (msec^-1) Altered from original to match ACh binding
        kL2F = 2.78  # (mM^-1 msec^-1)
        kL2R = 2.15e-5  # (msec^-1) Altered from original to match ACh binding
        kG1F = 0.00000027  # (um^2 ms^-1)
        kG1R = 0.0068  # (msec^-1)
        kG2F = 0.0000027  # (um^2 msec^-1)
        kG2R = 0.00068  # (msec^-1)
        k_NX_RLG = 0.00065  # *10.0  # (msec^-1)
        k_NX_G = 0.000000015  # (msec^-1)
        k_NX_P = 0.0047  # (msec^-1)
        kGTPase1 = 0.000026  # (msec^-1)
        kGTPase2 = 0.015  # (msec^-1)
        kPLCAssoc = 0.001 * 10.0  # 10.0  # (um^2 ms^-1) Alterred to increase rate of PLC
        kPLCDiss = 0.00071 * 10.0  # 10.0  # 0.00071  # (ms^-1)
        kReconst = 0.001  # (um^2 msec^-1)
        kPLC = 0.0003 * 100.0  # (um^2 msec^-1)

        # IP3R Gating, numbering based on Doi (2005) Supplemental Figure 1D
        kIP3R1F = 1000.0 # mM^-1 msec^-1
        kIP3R1R = 25.8 # msec^-1
        kIP3R2F = 8000.0
        kIP3R2R = 2.0
        kIP3R3F = 8.889
        kIP3R3R = 0.005
        kIP3R4F = 20.0
        kIP3R4R = 0.01
        kIP3R5F = 40.0
        kIP3R5R = 0.015
        kIP3R6F = 60.0
        kIP3R6R = 0.02
        gIP3r = self.rxdDict['channelDensities']['gIP3r']

        # Phosphatidylinositol Synthesis
        k4k = 0.0000008 * 24.9  # 0.0000228  # 0.0000008*28.5  # (msec^-1)
        k4p = 0.00012 * 8.3  # 0.00114  # 0.00012*9.5  # (msec^-1)
        k5k = 0.00002 * 8.3  # 0.00019  # 0.00002*9.5 # (msec^-1)
        k5p = 0.000028 * 8.3  # 0.000266  # 0.000028*9.5  # (msec^-1)

        kPIP2F = 0.001  # *10.0  # (msec^-1) [4]
        foldPIP2 = 3.0  # [4,5]
        kDAGase = 0.0002  # 0.0002 (ms^-1)

        KdKCNQPIP2 = 2000
        kKCNQF = 0.00005
        kKCNQR = KdKCNQPIP2 * kKCNQF


        # IP3 Decay
        kIP5PF = 59.0  # (mM^-1*msec^-1)
        kIP5PR = 0.072  # (msec^-1)
        KdIP5P = kIP5PR / kIP5PF
        kIP2 = 0.018  # (msec^-1)

        kIP3KCaF = 1111.1  # (mM^-1*msec^-1)
        kIP3KCaR = 0.1  # (msec^-1)
        KdIP3KCa = np.sqrt(kIP3KCaR / kIP3KCaF)

        kIP3KIP3F = 500  # (mM^-1*msec^-1)
        kIP3KIP3R = 0.08  # (msec^-1)
        KdIP3KIP3 = kIP3KIP3R / kIP3KIP3F
        kIP4 = 0.02  # (msec^-1)

        # SERCA
        self.VSERCA = 4.0  #
        self.KSERCA = 0.001  # mM
        self.tauSERCA = 1.0  # ms

        # 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"]["calbindin"]
        h0l0Init = totalCbd * 0.319
        h0l1Init = totalCbd * 0.157
        h0l2Init = totalCbd * 0.019
        h1l0Init = totalCbd * 0.236
        h1l1Init = totalCbd * 0.164
        h1l2Init = totalCbd * 0.0197
        h2l0Init = totalCbd * 0.0459
        h2l1Init = totalCbd * 0.0344
        h2l2Init = totalCbd * 0.00414

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

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

        # IP5P
        totalIP5P = self.rxdDict["species"]['IP5P']
        b_ip5p = totalIP5P * IP3ConcInit / (KdIP5P + IP3ConcInit)
        ub_ip5p = totalIP5P - b_ip5p

        pip2_soma_init = self.rxdDict["species"]['PIP2']
        pip2_axon_init = self.rxdDict["species"]['PIP2']
        b_kcnq_soma = 1724.1 * 0.005449  # (channels*um^-2) (based on 5.8 pS/channel [1] with 0.005449 S/cm^2)
        b_kcnq_axon = 1724.1 * 0.02647  # (channels*um^-2) (based on 5.8 pS S/channel with 0.02647 S/cm^2
        tot_kcnq_soma = b_kcnq_soma * (KdKCNQPIP2 + pip2_soma_init) / pip2_soma_init
        tot_kcnq_axon = b_kcnq_axon * (KdKCNQPIP2 + pip2_axon_init) / pip2_axon_init
        ub_kcnq_soma = tot_kcnq_soma - b_kcnq_soma
        ub_kcnq_axon = tot_kcnq_axon - b_kcnq_axon

        perc_kcnq_base_soma = (b_kcnq_soma / tot_kcnq_soma) ** 2
        perc_kcnq_base_axon = (b_kcnq_axon / tot_kcnq_axon) ** 2

        # IP3 Kinase
        total_ip3k = self.rxdDict['species']['IP3K']
        b_ip3k_ca = total_ip3k * caConcInitCyt / (KdIP3KCa + caConcInitCyt)
        ub_ip3k_ca = total_ip3k - b_ip3k_ca

        b_ip3k_ip3 = b_ip3k_ca * IP3ConcInit / (KdIP3KIP3 + IP3ConcInit)


        ########################################
        # 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, surface_fraction=1.0))
        # 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.ecs = rxd.Extracellular(-100, -100, -100, 100, 100, 100, dx=33)
        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))



        ####### 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=h0l0Init)
        self.cbdh1l0 = rxd.Species([self.cyt], name='cbdh1l0', initial=h1l0Init)
        self.cbdh2l0 = rxd.Species([self.cyt], name='cbdh2l0', initial=h2l0Init)
        self.cbdh0l1 = rxd.Species([self.cyt], name='cbdh0l1', initial=h0l1Init)
        self.cbdh1l1 = rxd.Species([self.cyt], name='cbdh1l1', initial=h1l1Init)
        self.cbdh2l1 = rxd.Species([self.cyt], name='cbdh2l1', initial=h2l1Init)
        self.cbdh0l2 = rxd.Species([self.cyt], name='cbdh0l2', initial=h0l2Init)
        self.cbdh1l2 = rxd.Species([self.cyt], name='cbdh1l2', initial=h1l2Init)
        self.cbdh2l2 = rxd.Species([self.cyt], name='cbdh2l2', initial=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)

        # Species declaration comes from Mergenthal et al., 2020. Mergenthal's model is publicly available on GitHub
        # Inositol Triphosphate
        self.ip3 = rxd.Species(self.cyt, d=IP3Diffusion, charge=0,
                               initial=IP3ConcInit)

        # Acetylcholine
        self.ach = rxd.Species(self.cyt, charge=0, initial=0.0)

        # Muscarinic mAChR
        self.rM1 = rxd.Species(self.cyt, initial=self.rxdDict['species']['M1'])
        self.gM1 = rxd.Species(self.cyt, initial=self.rxdDict['species']['GProtein'])
        self.gbgM1 = rxd.Species(self.cyt, initial=0.0)
        self.rlM1 = rxd.Species(self.cyt, initial=0.0)
        self.rgM1 = rxd.Species(self.cyt, initial=0.0)
        self.rgbgM1 = rxd.Species(self.cyt, initial=0.0)
        self.rlgM1 = rxd.Species(self.cyt, initial=0.0)
        self.rlgbgM1 = rxd.Species(self.cyt, initial=0.0)

        self.gaGdpM1 = rxd.Species(self.cyt, initial=0.0)
        self.gaGtpM1 = rxd.Species(self.cyt, initial=3.84e-6)  # 9.64e-6)
        self.plcM1 = rxd.Species(self.cyt, initial=self.rxdDict['species']['PLC'])
        self.gaGdpPlcM1 = rxd.Species(self.cyt, initial=0.0)
        self.gaGtpPlcM1 = rxd.Species(self.cyt, initial=3.097e-5)  # 1.32e-5)

        # Phosphatidylinositol Synthesis
        self.pi = rxd.Species(self.cyt, initial=self.rxdDict['species']['PI'])
        self.pi4p = rxd.Species(self.cyt, initial=self.rxdDict['species']['PI4P'])
        self.pip2 = rxd.Species(self.cyt, charge=0,
                                initial=self.rxdDict["species"]['PIP2'])

        self.pip2Bound = rxd.Species(self.cyt,
                                      initial=self.rxdDict["species"]['PIP2'] * (foldPIP2 - 1.0))
        self.dag = rxd.Species(self.cyt, initial=13.0)

        # PIP2 and KCNQ Binding
        # for sec in self.somaList:
        #     for seg in sec:
        #         sec(seg.x).base_perc_kmb_inh = perc_kcnq_base_soma
        #
        # for sec in self.axonList:
        #     for seg in sec:
        #         self.axonList[0](seg.x).base_perc_kmb_inh = perc_kcnq_base_axon
        #
        # self.kcnq = rxd.Species(self.cyt, charge=0, name='kcnq',
        #                         initial=lambda nd: ub_kcnq_axon if (nd.satisfies(self.axonList[0]) or nd.satisfies(
        #                             self.somaList[0])) else ub_kcnq_soma)
        # self.pip2Kcnq = rxd.Species(self.cyt, charge=0, name='pip2Kcnq',
        #                              initial=lambda nd: b_kcnq_axon if (nd.satisfies(self.axonList[0]) or nd.satisfies(
        #                                  self.somaList[0])) else b_kcnq_soma)

        # Inositol Triphosphate
        self.ip3 = rxd.Species(self.cyt, name='IP3', initial=IP3ConcInit)

        

        # IP3 Receptor
        self.ip3r = rxd.State(self.erCytMem, initial=0.814)
        self.ip3rc = rxd.State(self.erCytMem, initial=8.913e-5)
        self.ip3ro = rxd.State(self.erCytMem, initial=3.594e-5)
        self.ip3r1Ca = rxd.State(self.erCytMem, initial=0.146)
        self.ip3r2Ca = rxd.State(self.erCytMem, initial=0.0294)
        self.ip3r3Ca = rxd.State(self.erCytMem, initial=0.0079)
        self.ip3r4Ca = rxd.State(self.erCytMem, initial=0.00239)

        # IP3 Breakdown
        self.ip5p = rxd.Species(self.cyt, initial=lambda nd: ub_ip5p * nd.surface_area / nd.volume)  # 9.97e-4)
        self.ip5p_ip3 = rxd.Species(self.cyt, initial=lambda nd: b_ip5p * nd.surface_area / nd.volume)  # 2.32e-6)

        self.ip3k = rxd.Species(self.cyt, initial=lambda nd: ub_ip3k_ca * nd.surface_area / nd.volume)  # 0.0128)
        self.ip3k_2ca = rxd.Species(self.cyt, initial=lambda nd: b_ip3k_ca * nd.surface_area / nd.volume)  # 1.458e-7)
        self.ip3k_2ca_ip3 = rxd.Species(self.cyt,
                                        initial=lambda nd: b_ip3k_ip3 * nd.surface_area / nd.volume)  # 2.582e-9)

        # 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 or nd.segment in self.axonList)
        #                          else 0)

        # SERCA + ip3r_dense and param
        # self.ip3r_dense = rxd.Parameter(self.erCytMem,
        #                         value=lambda nd: 1.2 * nd.segment.diam / 4.0 if nd.satisfies(
        #                             self.somaList[0]) else nd.segment.diam / 4.0)

        # # AT BOTTOM:
        # # k_ip3r = gIP3r * (4*(1.0-self.x10_ip3r)*self.x10_ip3r**3 + self.x10_ip3r**4) * self.ip3r_dense  # For Bicknell and Goodhill 2016
        # # k_ip3r = gIP3r * self.ip3r_dense * self.ip3ro  # For Doi et al., 2005
        # # self.ip3r_flux = rxd.MultiCompartmentReaction(self.ca[self.er], self.ca[self.cyt], k_ip3r,
        # #                                               membrane=self.erCytMem)

        # scale = 602214.129
        # self.ip3_param = rxd.Parameter(self.cyt, value=lambda nd: nd.surface_area / scale / nd.volume)

        self.scale_serca = rxd.Parameter(self.erCytMem,
                                         value=lambda nd: 63525.0 * nd.segment.diam / 4.0)
        self.SERCAFluxInf = self.cai ** 2 / (self.KSERCA ** 2 + self.cai ** 2)
        self.SERCAState = rxd.State([self.cyt, self.erCytMem], name='SERCAstate', initial=3.96039604e-8)


        # Ca-induced Ca release
        CICRThresh = (1 + tanh(2 * 10000 *(self.cai - KtCICR))) / 2
        self.scale_cicr = rxd.Parameter([self.erCytMem],
                                        value=lambda nd: self.VCICR * nd.segment.diam / 4.0)
        CICRFluxInf = CICRThresh * self.cai / (self.KCICR + self.cai) * (self.caer - self.cai)
        self.CICRState = rxd.State([self.cyt, self.erCytMem], name='CICRstate', initial=1.749e-9)



        """
        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.
        """

        # SOMEWHERE BELOW IT BREAKS
        # BREAKS ON EVERY INSTANCE OF REACTION





        ####### WHAT? #######
        # Ca Extrusion by PMCA

        # print(self.ca, self.pmca, self.pmcaCa, kPMCABindF, kPMCABindR, self.cyt)
        self. = rxd.Reaction(self.ca + self.pmca, self.pmcaCa, kPMCABindF, kPMCABindR,
                                         region=[self.cyt])
        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, kCbdHF, kCbdHR, region = self.cyt)
        self.h0l0LBind = rxd.Reaction(self.ca + self.cbdh0l0, self.cbdh0l1, kCbdLF, kCbdLR, region = self.cyt)
        self.h0l1HBind = rxd.Reaction(self.ca + self.cbdh0l1, self.cbdh1l1, kCbdHF, kCbdHR,region = self.cyt)
        self.h0l1LBind = rxd.Reaction(self.ca + self.cbdh0l1, self.cbdh0l2, kCbdLF, kCbdLR,region = self.cyt)
        self.h0l2HBind = rxd.Reaction(self.ca + self.cbdh0l2, self.cbdh1l2, kCbdHF, kCbdHR,region = self.cyt)

        # Both above and below break it?
        

        self.h1l0LBind = rxd.Reaction(self.ca + self.cbdh1l0, self.cbdh1l1, kCbdLF, kCbdLR,region = self.cyt)
        self.h1l0HBind = rxd.Reaction(self.ca + self.cbdh1l0, self.cbdh2l0, kCbdHF, kCbdHR,region = self.cyt)
        self.h1l1LBind = rxd.Reaction(self.ca + self.cbdh1l1, self.cbdh1l2, kCbdLF, kCbdLR,region = self.cyt)
        self.h1l1HBind = rxd.Reaction(self.ca + self.cbdh1l1, self.cbdh2l1, kCbdHF, kCbdHR,region = self.cyt)
        self.h1l2HBind = rxd.Reaction(self.ca + self.cbdh1l2, self.cbdh1l2, kCbdHF, kCbdHR,region = self.cyt)

        self.h2l0LBind = rxd.Reaction(self.ca + self.cbdh2l0, self.cbdh2l1, kCbdLF, kCbdLR,region = self.cyt)
        self.h2l1LBind = rxd.Reaction(self.ca + self.cbdh2l1, self.cbdh2l2, kCbdLF, kCbdLR,region = self.cyt)
        # 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)

        

     

        # CICR
        self.dCICR = rxd.Rate(self.CICRState, (CICRFluxInf - self.CICRState) / tauCICR, regions=[self.erCytMem])
        self.CICRFlux = rxd.MultiCompartmentReaction(self.caer > self.cai, self.scale_cicr * self.CICRState,
                                                     membrane=self.erCytMem)

        # SERCA
        self.dSERCA = rxd.Rate(self.SERCAState, (self.SERCAFluxInf-self.SERCAState)/self.tauSERCA, regions=[self.erCytMem])
        self.SERCAFlux = rxd.MultiCompartmentReaction(self.cai > self.caer, self.scale_serca * self.SERCAState,
                                                      membrane=self.erCytMem, membrane_flux=True)

        # Muscarinic mAChR
        self.L1 = rxd.Reaction(self.rM1[self.cyt], self.rlM1[self.cyt],
                               kL1F * self.ach[self.cyt], kL1R, membrane=self.cyt)
        self.L2 = rxd.Reaction(self.rgM1[self.cyt], self.rlgM1[self.cyt],
                               kL2F * self.ach[self.cyt], kL2R, membrane=self.cyt)
        self.L2b = rxd.Reaction(self.rgbgM1[self.cyt], self.rlgbgM1[self.cyt],
                                kL2F * self.ach[self.cyt], kL2R, membrane=self.cyt)
        self.G1 = rxd.Reaction(self.gM1 + self.rM1, self.rgM1, kG1F, kG1R, region=[self.cyt])
        self.G2 = rxd.Reaction(self.gM1 + self.rlM1, self.rlgM1, kG2F, kG2R, region=[self.cyt])
        self.G1b = rxd.Reaction(self.gbgM1 + self.rM1, self.rgbgM1, kG1F, kG1R, region=[self.cyt])
        self.G2b = rxd.Reaction(self.gbgM1 + self.rlM1, self.rlgbgM1, kG2F, kG2R, region=[self.cyt])

        self.nx_g = rxd.Reaction(self.gM1 > self.gbgM1 + self.gaGtpM1, k_NX_G, region=[self.cyt])
        self.nx_rg = rxd.Reaction(self.rgM1 > self.rgbgM1 + self.gaGtpM1, k_NX_G, region=[self.cyt])
        self.nx_rlg = rxd.Reaction(self.rlgM1 > self.rlgbgM1 + self.gaGtpM1, k_NX_RLG, region=[self.cyt])
        self.nx_p = rxd.Reaction(self.gaGdpPlcM1 > self.gaGtpPlcM1, k_NX_P, region=[self.cyt])

        self.gtpase1 = rxd.Reaction(self.gaGtpM1 > self.gaGdpM1, kGTPase1, region=[self.cyt])
        self.gtpase2 = rxd.Reaction(self.gaGtpPlcM1 > self.gaGdpPlcM1, kGTPase2, region=[self.cyt])

        self.plc_assoc = rxd.Reaction(self.plcM1 + self.gaGtpM1 > self.gaGtpPlcM1, kPLCAssoc, region=[self.cyt])
        self.plc_diss = rxd.Reaction(self.gaGdpPlcM1 > self.gaGdpM1 + self.plcM1, kPLCDiss, region=[self.cyt])
        self.g_reconst = rxd.Reaction(self.gbgM1 + self.gaGdpM1 > self.gM1, kReconst, region=[self.cyt])
        # Phosphatidylinositol Synthesis
        self.k4 = rxd.Rate(self.pi4p, k4k * self.pi, regions=[self.cyt])
        self.p4 = rxd.Rate(self.pi4p, -k4p * self.pi4p, regions=[self.cyt])
        self.k5_p5 = rxd.Reaction(self.pi4p, self.pip2, k5k, k5p, region=[self.cyt])

        # IP3 Generation and Decay
        self.pip2_hydrolysis_dag = rxd.Reaction(self.pip2 > self.dag,
                                                kPLC * foldPIP2 * self.gaGtpPlcM1[self.cyt],
                                                region=[self.cyt])


        scale = 602214.129
        # COMMENT 1 moved up to SERCA
        # print(nd.volume)
        # print(value = lambda nd: nd.volume)
        self.ip3_param = rxd.Parameter(self.cyt, value=lambda nd: nd + 1)
        self.ip3_param = rxd.Parameter(self.cyt, value=lambda nd: nd.surface_area / scale / nd.volume)

        self.ip3_flux = rxd.Rate(self.ip3, self.pip2 * kPLC * foldPIP2 * self.gaGtpPlcM1 * self.ip3_param,
                                 regions=[self.cyt])

        self.ip5p_binding = rxd.Reaction(self.ip3 + self.ip5p, self.ip5p_ip3, kIP5PF, kIP5PR, region=self.cyt)
        self.ip2_form = rxd.Reaction(self.ip5p_ip3 > self.ip5p, kIP2, region=self.cyt)

        self.ip3k_cal_bind = rxd.Reaction(self.ip3k + 2 * self.ca, self.ip3k_2ca, kIP3KCaF, kIP3KCaR,
                                          region=self.cyt)
        self.ip3k_ip3_bind = rxd.Reaction(self.ip3k_2ca + self.ip3, self.ip3k_2ca_ip3, kIP3KIP3F, kIP3KIP3R,
                                          region=self.cyt)
        self.ip4_form = rxd.Reaction(self.ip3k_2ca_ip3 > self.ip3k_2ca, kIP4, region=self.cyt)

        self.dagase = rxd.Rate(self.dag, -kDAGase * self.dag, regions=[self.cyt])
        # PIP2 Buffering and Binding to KCNQ
        # self.pip2_buffer = rxd.Reaction(self.pip2, self.pip2_bound, (foldPIP2 - 1) * kPIP2F, kPIP2F,
        #                                 regions=[self.cyt])
        # self.kcnq_bind_pip2 = rxd.Reaction(self.pip2 + self.kcnq, self.pip2_kcnq, kKCNQF,
        #                                    kKCNQR, regions=[self.cyt])

        # IP3 Gating
        # d1 = kIP3R1F * self.ip3r * self.IP3[self.cyt] - kIP3R1R * self.ip3rc
        # d2 = kIP3R2F * self.ip3rc * self.cai - kIP3R2R * self.ip3ro
        # d3 = kIP3R3F * self.ip3r * self.cai - kIP3R3R * self.ip3r1Ca
        # d4 = kIP3R4F * self.ip3r1Ca * self.cai - kIP3R4R * self.ip3r2Ca
        # d5 = kIP3R5F * self.ip3r2Ca * self.cai - kIP3R5R * self.ip3r3Ca
        # d6 = kIP3R6F * self.ip3r3Ca * self.cai - kIP3R6R * self.ip3r4Ca
        #
        # self.dIP3R = rxd.Rate(self.ip3r, -d1-d3, regions=[self.erCytMem])
        # self.dIP3Rc = rxd.Rate(self.ip3rc, d1-d2, regions=[self.erCytMem])
        # self.dIP3Ro = rxd.Rate(self.ip3ro, -d2, regions=[self.erCytMem])
        # self.dIP3R1Ca = rxd.Rate(self.ip3r1Ca, -d4+d3, regions=[self.erCytMem])
        # self.dIP3R2Ca = rxd.Rate(self.ip3r2Ca, d4-d5, regions=[self.erCytMem])
        # self.dIP3R3Ca = rxd.Rate(self.ip3r3Ca, d5-d6, regions=[self.erCytMem])
        # self.dIP3R4Ca = rxd.Rate(self.ip3r4Ca, d6, regions=[self.erCytMem])

        # From Mergenthal 2020, based on Doi et al., 2005
        react1 = self.ip3r * self.ca[self.cyt] * kIP3R2F - kIP3R2R * self.ip3ro
        react2 = self.ip3rc * self.ip3[self.cyt] * kIP3R1F - kIP3R1R * self.ip3rc
        react3 = self.ip3rc * self.ca[self.cyt] * kIP3R3F - self.ip3r1Ca * kIP3R3R
        react4 = self.ip3r1Ca * self.ca[self.cyt] * kIP3R4F - self.ip3r2Ca * kIP3R4R
        react5 = self.ip3r2Ca * self.ca[self.cyt] * kIP3R5F - self.ip3r3Ca * kIP3R5R
        react6 = self.ip3r3Ca * self.ca[self.cyt] * kIP3R6F - self.ip3r4Ca * kIP3R6R

        self.dip3r = rxd.Rate(self.ip3r, - react2 - react3, regions=[self.erCytMem])
        self.dip3rc = rxd.Rate(self.ip3rc, -react1 + react2, regions=[self.erCytMem])
        self.dip3ro = rxd.Rate(self.ip3ro, react1, regions=[self.erCytMem])
        self.dip3r1Ca = rxd.Rate(self.ip3r1Ca, react3 - react4, regions=[self.erCytMem])
        self.dip3r2Ca = rxd.Rate(self.ip3r2Ca, react4 - react5, regions=[self.erCytMem])
        self.dip3r3Ca = rxd.Rate(self.ip3r3Ca, react5 - react6, regions=[self.erCytMem])
        self.dip3r4Ca = rxd.Rate(self.ip3r4Ca, react6, regions=[self.erCytMem])
        # COMMENT2
        # self.ip3r_dense = rxd.Parameter(self.erCytMem,
        #                                 value=lambda nd: 1.2 * nd.segment.diam / 4.0 if nd.satisfies(
        #                                     self.somaList[0]) else nd.segment.diam / 4.0)
        k_ip3r = gIP3r * (4*(1.0-self.ip3r)*self.ip3r**3 + self.ip3r**4) * self.ip3r_dense
        k_ip3r = gIP3r * self.ip3r_dense * self.ip3ro  # For Doi et al., 2005
        self.ip3r_flux = rxd.MultiCompartmentReaction(self.ca[self.er], self.ca[self.cyt], k_ip3r,
                                                      membrane=self.erCytMem)


        # SOMEWHERE ABOVE IT BREAKS 


    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.
        '''

        '''
        START
        self.tvec = h.Vector().record(h._ref_t)
        self.vvec = h.Vector().record(measureNeuron(location)._ref_v)

        self.mvec = h.Vector().record(measureNeuron(location).spike._ref_m)
        self.nvec = h.Vector().record(measureNeuron(location).spike._ref_n)
        self.hvec = h.Vector().record(measureNeuron(location).spike._ref_h)
        self.pvec = h.Vector().record(measureNeuron(location).spike._ref_p)
        self.qvec = h.Vector().record(measureNeuron(location).spike._ref_q)
        self.cvec = h.Vector().record(measureNeuron(location).spike._ref_c)


        # 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.ip3vec = h.Vector().record(self.ip3.nodes(measureNeuron(location))._ref_concentration)

        self.achvec = h.Vector().record(self.ach.nodes(measureNeuron(location))._ref_concentration)


        self.icakvec = h.Vector().record(measureNeuron(location)._ref_ica)
        self.ikvec = h.Vector().record(measureNeuron(location)._ref_ik)



        # 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)

        # self.cicrthreshvec = h.Vector().record(self.CICRThresh[self.cyt].nodes(measureNeuron(location))._ref_value)
        #
        #
        # self.cicrvec = h.Vector().record(self.CICRState[self.cyt].nodes(measureNeuron(location))._ref_value)

        self.sercavec = h.Vector().record(self.SERCAState[self.cyt].nodes(measureNeuron(location))._ref_value)
        # h.cvode.active(1)
        #
        # h.v_init = -70.0
        #
        # h.tstop = sim_dur
        # h.celsius = 34.0


        END
        '''


        # fih = h.FInitializeHandler("finitialize()", [h])
        # fih.allprint()
        print("cp0")
        h.finitialize(-60)
        h.continuerun(sim_dur)

        # self.CICRInfvec = self.VCICR * self.caivec / (self.KCICR ** 2 + self.caivec) * (self.caervec - self.caivec)





    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 = RGCMorphology()
    vgcc = VGCC_2(rgc1, 'rxd_file.json')
    vgcc.insertChannels()
    vgcc.implement_RXD(rgc1)
    # vgcc.stimulate(vgcc.somaList[0], vgcc.axonList[0], 0.5, 0, 0.0001, 50, 100)
    # print(vgcc.somaList[0].diam)
    vgcc.stimulate(vgcc.somaList[0], vgcc.somaList[0], 0.5, 0, 0.4, 50, 100)

    print("ENDED STIMULATION I FUCKING DID IT")



    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.pvec.as_numpy(), '-b', label='p-gate')
    pyplot.plot(vgcc.tvec.as_numpy(), vgcc.qvec.as_numpy(), '-r', label='q-gate')
    pyplot.plot(vgcc.tvec.as_numpy(), vgcc.cvec.as_numpy(), '-g', label='c-gate')
    pyplot.plot(vgcc.tvec.as_numpy(), vgcc.mvec.as_numpy(),  label='m-gate')
    pyplot.plot(vgcc.tvec.as_numpy(), vgcc.hvec.as_numpy(),  label='h-gate')
    pyplot.plot(vgcc.tvec.as_numpy(), vgcc.nvec.as_numpy(),  label='n-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.ikvec.as_numpy(), '-g', label='k')
    # pyplot.xlabel('t (ms)')
    # pyplot.ylabel('Current (mA/cm^2)')
    # pyplot.savefig(os.path.join('CaKCurrentTest.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.as_numpy(), vgcc.ip3vec.as_numpy(), '-g', label='k')
    pyplot.xlabel('time (ms)')
    pyplot.ylabel('IP3 Concentration (mM)')
    pyplot.savefig(os.path.join('IP3Test.svg'), format="svg")

    # fig = pyplot.figure()
    # pyplot.plot(vgcc.tvec.as_numpy(), vgcc.achvec.as_numpy(), '-g', label='k')
    # pyplot.xlabel('time (ms)')
    # pyplot.ylabel('ACh Concentration (mM)')
    # pyplot.savefig(os.path.join('AChTest.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")


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

    # 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")

    # fig = pyplot.figure()
    # pyplot.plot(vgcc.tvec.as_numpy(), vgcc.cicrthreshvec.as_numpy(), '-g', label='k')
    # pyplot.xlabel('time (ms)')
    # pyplot.ylabel('CICR Flux mA/cm^2)')
    # pyplot.savefig(os.path.join('CICRTest.svg'), format="svg")



    # print(vgcc.caivec[-1])
    # print(vgcc.soma(0.5)._ref_cai)


    rgc1 = vgcc = None

Thank you for any help!
neuron121!
Posts: 8
Joined: Thu Mar 31, 2022 11:00 pm

Re: rxd.Reaction(...) causes finitialize to not work

Post by neuron121! »

Just for a follow-up, after downloading Neuron 8.2.0, there is the SAME error every time that there is ANY rxd.Reaction() call followed by h.finitialize(). This occurs straight after downloading, and I'm not sure what to do about this error. If there are any suggestions, I would be grateful. I have a Windows 64-bit using numpy 1.22.2 and Neuron 8.2.0 (which works on my collaborators Mac).
Post Reply