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;
}
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!