I got some python code sent to me that works on their end, but it gets the following error on my end:
Code: Select all
# lots of similar lines
...
...
C:\Users\jdami\AppData\Local\Temp\ccWASCyx.s:1829: Error: bad register name `%rbp)'
C:\Users\jdami\AppData\Local\Temp\ccWASCyx.s:1830: Error: bad register name `%rsp'
C:\Users\jdami\AppData\Local\Temp\ccWASCyx.s:1831: Error: bad register name `%rbp'
C:\Users\jdami\AppData\Local\Temp\ccWASCyx.s:1833: Warning: zero assumed for missing expression
C:\Users\jdami\AppData\Local\Temp\ccWASCyx.s:1833: Warning: zero assumed for missing expression
Traceback (most recent call last):
File "c:\nrn\lib\python\neuron\rxd\rxd.py", line 1872, in _init
_compile_reactions()
File "c:\nrn\lib\python\neuron\rxd\rxd.py", line 1484, in _compile_reactions
_c_compile(fxn_string),
File "c:\nrn\lib\python\neuron\rxd\rxd.py", line 564, in _c_compile
dll = ctypes.cdll["%s.so" % os.path.abspath(filename)]
File "C:\Users\jdami\AppData\Local\Programs\Python\Python310\lib\ctypes\__init__.py", line 449, in __getitem__
return getattr(self, name)
File "C:\Users\jdami\AppData\Local\Programs\Python\Python310\lib\ctypes\__init__.py", line 444, in __getattr__
dll = self._dlltype(name)
File "C:\Users\jdami\AppData\Local\Programs\Python\Python310\lib\ctypes\__init__.py", line 374, in __init__
FileNotFoundError: Could not find module 'C:\Users\jdami\Documents\Year 4\Lazzi\Colab\ITEMS_AxonGrowth\rxddll9e37c954-e82a-11ec-b65e-7470fdcdedcf.so' (or one of its dependencies). Try using the full path with constructor syntax.
NEURON: Python Callback failed
near line 0
^
finitialize()
Traceback (most recent call last):
vgcc.stimulate(vgcc.somaList[0], vgcc.somaList[0], 0.5, 0, 0.5, 60, 100)
File "c:\Users\jdami\Documents\Year 4\Lazzi\Colab\ITEMS_AxonGrowth\CaSignalingPMCAClass.py", line 588, in stimulate
h.finitialize(-70 * mV)
RuntimeError: hocobj_call error
Main file:
Code: Select all
from neuron import h, rxd
import numpy as np
import json
from matplotlib import pyplot
from neuron.units import mM, mV, nM, uM
from neuron.rxd import v
from neuron.rxd.rxdmath import vtrap, exp, log, tanh, fabs
import RGCMorphology as rgcMorph
# import calcium_optimization_init_conds as cbdInit
import os
h.load_file('stdrun.hoc')
h.load_file("import3d.hoc")
h.load_file("stdlib.hoc")
class VGCC_2:
def __init__(self, neuronMorph, JSONFile):
self.axonList = [sec for sec in neuronMorph.axon]
self.somaList = [sec for sec in neuronMorph.soma]
self.apicalList = [sec for sec in neuronMorph.soma]
self.all = self.axonList + self.somaList + self.apicalList
# self.somaList[0].pt3dadd(-300,0,0,self.somaList[0].diam)
# self.somaList[0].pt3dadd(-300+self.somaList[0].L, 0, 0, self.somaList[0].diam)
with open(JSONFile, 'r') as fileObj:
self.rxdDict = json.load(fileObj)
self.findExtrema()
self.padding = 100
def create_neuron(self):
'''
THINGS TO THINK ABOUT: Do we really need a function that creates a new neuron? Perhaps we can pass a neuron as a parameter
and implement VGCC's in that neuron, so we may not even need this function.
:param self:
:return:
'''
self.axon = h.Section(name="axon")
self.soma = h.Section(name="soma")
self.soma.pt3dclear()
self.soma.pt3dadd(-90, 0, 0, 30)
self.soma.pt3dadd(-60, 0, 0, 30)
self.axon.nseg = 11
self.soma.nseg = 11
self.axon.connect(self.soma)
self.axon.diam = 10
self.axon.L = 100
return self.soma
def findExtrema(self):
'''
Finds the extrema of the morphology, and stores them into variables xmin, ymin, zmin, xmax, ymax, and zmax.
These coordinates will be used in defining the extracellular space in 3D using rxd.Extracellular()
:return: NONE
'''
self.xmin = self.ymin = self.zmin = self.xmax = self.ymax = self.zmax = None
for sec in self.all:
n3d = sec.n3d()
xList = [sec.x3d(i) for i in range(n3d)]
yList = [sec.y3d(i) for i in range(n3d)]
zList = [sec.z3d(i) for i in range(n3d)]
xMinSec, yMinSec, zMinSec = min(xList), min(yList), min(zList)
xMaxSec, yMaxSec, zMaxSec = max(xList), max(yList), max(zList)
if self.xmin is None:
self.xmin, self.ymin, self.zmin = xMinSec, yMinSec, zMinSec
self.xmax, self.ymax, self.zmax = xMaxSec, yMaxSec, zMaxSec
else:
self.xmin, self.ymin, self.zmin = min(self.xmin,xMinSec), min(self.ymin,yMinSec), min(self.zmin,zMinSec)
self.xmax, self.ymax, self.zmax = max(self.xmax,xMaxSec), max(self.ymax,yMaxSec), max(self.zmax,zMaxSec)
def concInit(self, cytc, exc, erc):
return lambda nd: exc if nd.region == self.ecs else (cytc if nd.region == self.cyt else erc)
def implement_RXD(self, neuronMorph):
########################################
# Parameters
# Initial ion concentrations
naConcInitCyt = self.rxdDict["ions"]["na"]["cyt"]
naConcInitEcs = self.rxdDict["ions"]["na"]["ecs"]
kConcInitCyt = self.rxdDict["ions"]["k"]["cyt"]
kConcInitEcs = self.rxdDict["ions"]["k"]["ecs"]
caConcInitCyt = self.rxdDict["ions"]["ca"]["cyt"]
caConcInitEcs = self.rxdDict["ions"]["ca"]["ecs"]
caConcInitEr = self.rxdDict["ions"]["ca"]["er"]
h.celsius = 37
self.k = 1.38064852e-23 # Boltzmann constant
self.T = h.celsius + 273.15 # Temperature in Kelvin
# Parameters between this comment and the next comment comes from Hodgkin Huxley tutorial on NEURON website.
self.e = 1.60217662e-19
self.scale = 1e-14 / self.e # convert S/cm2 to molecules/um2
# Make sure to use CHANNEL DENSITIES. Below max conductance values
# come from Targeted Simulation of RGC in Epiretinal Prostheses... Paknahad et al. (2020) (Soma of A2 cell type)
self.gNaBar = 0.3 * self.scale # molecules/um2 ms mV
self.gKBar = 0.12 * self.scale # molecules/um2 ms mV
self.gKABar = 0 #3 * self.gKBar
self.gHBar = 0.012 * self.scale
self.gCaBarL = .0025 * self.scale ## molecules/um2 ms mV
self.gCaBarN = 0.0025 * self.scale # molecules/um2 ms mV
self.gCaBarT = 0.00025 * self.scale # molecules/um2 ms mV
self.gL = 0.0003 * self.scale # molecules/um2 ms mV, assuming 30 um soma diameter and transmembrane resistance of 20000 Ohm*cm^2
self.el = -70
self.q10 = 3.0 ** ((h.celsius - 37.0) / 10.0)
self.erFraction = self.rxdDict["region"]["erFraction"]
self.cytFraction = 1 - self.erFraction
# Rate Constants
# PMCA
kPMCABindF = 25000.0 # msec^-1 mM^-1
kPMCABindR = 2.0 # msec^-1
kPMCARel = 500.0 # msec^-1
# NCX
kNCXBindF = 93.827 #93.827 msec^-1 mM^-1
kNCXBindR = 4.0 # msec^-1
kNCXRel = 1000.0 # msec^-1 mM^-1
# SERCA
kSERCABindF = 17147.0
kSERCABindR = 1.0
kSERCARel = 250.0
# Calreticulin
KdCalr = 2.0
kCalrBindF = 0.01
kCalrBindR = KdCalr * kCalrBindF
# Calbindin-D282k
KdHigh = 237e-6
kCbdHF = 11.0
kCbdHR = KdHigh * kCbdHF
KdLow = 411e-6
kCbdLF = 87.0
kCbdLR = KdLow * kCbdLF
# Parvalbumin
KdParv = 0.5315
kParvBindF = 18.5
kParvBindR = KdParv * kParvBindF
# Initialize channel densities (steady-state) in respective kinetic states according to Hill-Langmuir equation
# PMCA
totalPMCA = self.rxdDict["channelDensities"]["totalPMCA"] # Channel density of PMCA from Adam's paper. NEED TO calibrate
pmcaCaInit = totalPMCA * caConcInitCyt / ((kPMCABindR / kPMCABindF) + caConcInitCyt)
pmcaInit = totalPMCA - pmcaCaInit
# NCX
totalNCX = self.rxdDict["channelDensities"]["totalNCX"]
ncx2CaInit = totalNCX * (caConcInitCyt**2) / ((kNCXBindR/kNCXBindF) + (caConcInitCyt**2))
ncxInit = totalNCX - ncx2CaInit
# SERCA
totalSERCA = self.rxdDict["channelDensities"]["totalSERCA"]
serca2CaInit = totalSERCA * (caConcInitCyt ** 2) / ((kSERCABindR / kSERCABindF) + (caConcInitCyt ** 2))
sercaInit = totalSERCA - serca2CaInit
# Calbindin-D28k
totalCbd = self.rxdDict["endoBuffers"]["totalCbd"]
# h0l0Init, h0l1Init, h0l2Init, h1l0Init, h1l1Init, h1l2Init, h2l0Init, h2l1Init, h2l2Init = cbdInit.get_init_conds_calbindin(totalCbd, caConcInitCyt)
# Calreticulin
totalCalr = self.rxdDict["endoBuffers"]["totalCalr"]
calrCaInit = totalCalr * caConcInitCyt / ((kCalrBindR / kCalrBindF) + caConcInitCyt)
calrInit = totalCalr - calrCaInit
# Parvalbumin
totalParv = self.rxdDict["endoBuffers"]["totalParv"]
parvCaInit = totalParv * caConcInitCyt / ((kCalrBindR / kCalrBindF) + caConcInitCyt)
parvInit = totalParv - parvCaInit
########################################
# Hodgkin-Huxley equations for voltage-gated ion channels
'''
Note: "vtrap(x,y) is 1/(exp(x/y)-1) if |x/y|>=1e-6 or y*(1.0 - x/y/2.0) otherwise"
https://neuron.yale.edu/neuron/docs/hodgkin-huxley-using-rxd
vtrap is used to avoid singularities in equations that define transition probabilities.
In some equations, denominator goes to zero, so vtrap offers a workaround.
Below HH equations come from Targeted Simulation of RGC in Epiretinal Prostheses... Paknahad et al. (2020)
'''
# Na channel
# m-gate
self.alpha = 0.1 * vtrap(-(v + 40.0), 10.0)
self.beta = 4.0 * exp(-(v + 65.0) / 18.0)
self.mtau = 1.0 / (self.q10 * (self.alpha + self.beta))
self.minf = self.alpha / (self.alpha + self.beta)
# self.alpha = (-0.6 * (v + 30.0)) * vtrap(-(v + 30.0), 10.0)
# self.beta = 20.0 * exp(-(v + 55.0) / 18.0)
# self.mtau = 1.0 / (self.q10 * (self.alpha + self.beta))
# self.minf = self.alpha / (self.alpha + self.beta)
# h-gate
self.alpha = 0.07 * exp(-(v + 65.0) / 20.0)
self.beta = 1.0 / (exp(-(v + 35.0) / 10.0) + 1.0)
self.htau = 1.0 / (self.q10 * (self.alpha + self.beta))
self.hinf = self.alpha / (self.alpha + self.beta)
# self.alpha = 0.4 * exp(-(v + 50.0) / 20.0)
# self.beta = 6.0 / (exp(-0.1 * (v + 20.0) / 10.0) + 1.0)
# self.htau = 1.0 / (self.q10 * (self.alpha + self.beta))
# self.hinf = self.alpha / (self.alpha + self.beta)
# K channel (Delayed-rectifier)
# n-gate
self.alpha = 0.01 * vtrap(-(v + 55.0), 10.0)
self.beta = 0.125 * exp(-(v + 65.0) / 80.0)
self.ntau = 1.0 / (self.q10 * (self.alpha + self.beta))
self.ninf = self.alpha / (self.alpha + self.beta)
# self.alpha = -0.02 * (v + 40.0) * vtrap(-(v + 40.0), 10.0)
# self.beta = 0.4 * exp(-(v + 50.0) / 80.0)
# self.ntau = 1.0 / (self.q10 * (self.alpha + self.beta))
# self.ninf = self.alpha / (self.alpha + self.beta)
# K,A channel
# m-gate
# self.alpha = -0.006 * (v + 90.0) * vtrap(-(v + 90.0), 10.0)
# self.beta = 0.1 * exp(-(v + 30.0) / 10.0)
# self.mtauA = 1.0 / (self.q10 * (self.alpha + self.beta))
# self.minfA = self.alpha / (self.alpha + self.beta)
#
# # h-gate
# self.alpha = 0.04 * exp(-(v + 70.0) / 20.0)
# self.beta = 0.6 / (exp(-(v + 40.0) / 10.0) + 1)
# self.htauA = 1.0 / (self.q10 * (self.alpha + self.beta))
# self.hinfA = self.alpha / (self.alpha + self.beta)
#
# # H channel
# # l-gate
# self.alpha = exp(0.08316 * (v + 75.0))
# self.beta = exp(0.033264 * (v + 75.0))
# self.ltauH = 1.0 / (self.q10 * (self.alpha + self.beta))
# self.linfH = self.alpha / (self.alpha + self.beta)
'''
The equations come from "VDCC documentation v1.3bis" Word Document file on AXON REGROWTH Google Drive folder.
'''
# Ca channel
# T-channel
# # m-gate
# self.alpha = 0.2 * vtrap((-v + 19.26), 10)
# self.beta = 9.0e-3 * exp(-v / 22.03)
# self.minf_T = self.alpha / (self.alpha + self.beta)
# self.tau_Tm = 1.0 / (self.q10 * (self.alpha + self.beta))
# # h-gate
# self.alpha = 1e-6 * exp(-v / 16.26)
# self.beta = 1 / (exp((-v + 29.79) / 10) + 1)
# self.hinf_T = self.alpha / (self.alpha + self.beta)
# self.tau_Th = 1.0 / (self.q10 * (self.alpha + self.beta))
# N-channel
# m-gate
self.alpha = 0.19 * vtrap((-v + 19.88), 10)
self.beta = 4.6e-2 * exp(-v / 20.73)
self.minf_N = self.alpha / (self.alpha + self.beta)
self.tau_Nm = 1.0 / (self.q10 * (self.alpha + self.beta))
# h-gate
self.alpha = 1.6e-4 * exp(-v / 48.4)
self.beta = 1 / (exp((-v + 39) / 10) + 1)
self.hinf_N = self.alpha / (self.alpha + self.beta)
self.tau_Nh = 1.0 / (self.q10 * (self.alpha + self.beta))
# L-channel
# m-gate
self.alpha = 15.69 * vtrap((-v + 81.5), 10)
self.beta = 0.29 * exp(-v / 10.86)
self.minf_L = self.alpha / (self.alpha + self.beta)
self.tau_Lm = 1.0 / (self.q10 * (self.alpha + self.beta))
'''
This code uses the complex equation of Ca current. A simple if-statement cannot compare crxd.v to a number due to how
crxd.v is defined. A workaround uses y = (1 + tanh(x - threshold))/2, which transitions from 0 to 1 fairly quickly,
allowing the function to be used as a binary switch that outputs 1 if our value is greater than threshold or 0 if value is
less than threshold. As a result, temp3 stores either 1 - temp2/2 or temp2/exp(temp) - 1.
'''
# Formulas for Non-modulated Ca currents
self.threshold = 1.10e-4
self.m = 10000 # steepness of switch
self.temp1 = 0.0853 * self.T / 2
self.temp2 = v / self.temp1
# switch goes to 0 for abs(temp2) < threshold and 1 for abs(temp2) > threshold
self.switch = (1 + tanh(((fabs(self.temp2) - self.threshold) * 100))) / 2
self.temp3 = (1 - self.switch) * (1 - self.temp2 / 2) + self.switch * (self.temp2 / exp(self.temp2) - 1)
########################################
# RXD
'''
Code that uses RXD should define the following criteria: WHERE the reaction is happening, WHO (what species) are we interested
in, and WHAT about those species of interest? With that being said, the following code
answers, 'where?' The reaction is happening in: the cytoplasm, cell membrane, and extracellular fluid.
'''
####### WHERE? #######
self.cyt = rxd.Region(self.all, name='cyt', nrn_region='i', geometry = rxd.FractionalVolume(self.cytFraction))
self.mem = rxd.Region(self.all, name='mem', geometry=rxd.membrane())
# self.ecs = rxd.Region(self.all, name='ecs', nrn_region = 'o')
self.ecs = rxd.Extracellular(self.xmin-self.padding, self.ymin-self.padding, self.zmin-self.padding,
self.xmax+self.padding, self.ymax+self.padding, self.zmax+self.padding,
dx=50)
self.er = rxd.Region(self.all, name='er', geometry=rxd.FractionalVolume(self.erFraction))
self.erCytMem = rxd.Region(self.all, name='erCytMem', geometry=rxd.ScalableBorder(0.1, on_cell_surface=False))
'''
This code aims to graph Ca current in order to study Ca signaling within a neuron. Hodgkin Huxley
model will be able to track current of Na, K, and Ca ions flowing in/out of neuron.
'''
####### WHO? #######
# Ions
self.k = rxd.Species([self.cyt, self.ecs], name="k", d=1, charge=1,
initial=self.concInit(kConcInitCyt, kConcInitEcs,0))
self.na = rxd.Species([self.cyt, self.ecs], name="na", d=1, charge=1,
initial=self.concInit(naConcInitCyt, naConcInitEcs, 0))
self.ca = rxd.Species([self.cyt, self.ecs, self.er], name='ca', d=1, charge=2,
initial=self.concInit(caConcInitCyt, caConcInitEcs, caConcInitEr))
# Plasma Membrane Ca ATPase
self.pmca = rxd.Species([self.cyt], name="pmca", initial=pmcaInit)
self.pmcaCa = rxd.Species([self.cyt], name="pmcaCa", initial=pmcaCaInit) # PMCA bound to Ca
# Na/Ca Exchange
self.ncx = rxd.Species([self.cyt], name='ncx', initial=ncxInit)
self.ncx2Ca = rxd.Species([self.cyt], name='ncx2Ca', initial=ncx2CaInit)
# Sarcoplasmic/Endoplasmic Reticulum Ca ATPase
self.serca = rxd.Species(self.cyt, name='serca', initial=sercaInit)
self.serca2Ca = rxd.Species(self.cyt, name='serca2Ca', initial=serca2CaInit)
# # Calbindin
# self.cbdh0l0 = rxd.Species([self.cyt], name='cbdh0l0', initial=self.h0l0Init)
# self.cbdh1l0 = rxd.Species([self.cyt], name='cbdh1l0', initial=self.h1l0Init)
# self.cbdh2l0 = rxd.Species([self.cyt], name='cbdh2l0', initial=self.h2l0Init)
# self.cbdh0l1 = rxd.Species([self.cyt], name='cbdh0l1', initial=self.h0l1Init)
# self.cbdh1l1 = rxd.Species([self.cyt], name='cbdh1l1', initial=self.h1l1Init)
# self.cbdh2l1 = rxd.Species([self.cyt], name='cbdh2l1', initial=self.h2l1Init)
# self.cbdh0l2 = rxd.Species([self.cyt], name='cbdh0l2', initial=self.h0l2Init)
# self.cbdh1l2 = rxd.Species([self.cyt], name='cbdh1l2', initial=self.h1l2Init)
# self.cbdh2l2 = rxd.Species([self.cyt], name='cbdh2l2', initial=self.h2l2Init)
# Calreticulin
self.calr = rxd.Species(self.er, name='calr', initial=calrInit)
self.calrCa = rxd.Species(self.er, name='calrCa', initial=calrCaInit)
# Parvalbumin
self.parv = rxd.Species(self.cyt, name='parv', initial=parvInit)
self.parvCa = rxd.Species(self.cyt, name='parvCa', initial=parvCaInit)
# Leak
self.x = rxd.Species([self.cyt, self.ecs], name='x', charge=1)
self.ki = self.k[self.cyt]
self.ko = self.k[self.ecs]
self.nai = self.na[self.cyt]
self.nao = self.na[self.ecs]
self.xi = self.x[self.cyt]
self.xo = self.x[self.ecs]
self.cai = self.ca[self.cyt]
self.cao = self.ca[self.ecs]
self.caer = self.ca[self.er]
# To limit reactions to soma and axon
inRegion = rxd.Parameter([self.cyt], name = 'inRegion',
value = lambda nd: 1
if (nd.segment in self.somaList[0] or nd.segment in self.axonList)
else 0)
# gating states -> open probabilities
# K+ channel (delayed rectifier)
self.ngate = rxd.State([self.cyt, self.mem], name='ngate', initial=0.2445865)
# # K+ Channel (A channel)
# self.mgateA = rxd.State([self.cyt, self.mem], name='mgateA', initial=0.0247887308)
# self.hgateA = rxd.State([self.cyt, self.mem], name='hgateA', initial=0.543209973)
#
# # H channel (K+)
# self.lgateH = rxd.State([self.cyt, self.mem], name='lgateH', initial=0.6416789483)
#
# Na+ Channel
self.mgate = rxd.State([self.cyt, self.mem], name='mgate', initial=0.0289055)
self.hgate = rxd.State([self.cyt, self.mem], name='hgate', initial=0.7540796)
# Ca+2 Channel
# mgate_T = rxd.State([cyt, mem], name='mgate_T', initial=0.010871215)
# hgate_T = rxd.State([cyt, mem], name='hgate_T', initial=0.615047335)
self.mgate_N = rxd.State([self.cyt, self.mem], name='mgate_N', initial=0.001581552)
self.hgate_N = rxd.State([self.cyt, self.mem], name='hgate_N', initial=0.973556944)
self.mgate_L = rxd.State([self.cyt, self.mem], name='mgate_L', initial=3.42574e-6)
"""
NOTE: 'mgate' and 'm_gate' are different. 'mgate' is an object of class rxd.State, representing
open probability value. We are telling NEURON 'who' we are interested in, which in this case are the
various types of gates.
'm_gate'is an rxd.Rate object which we are telling NEURON to track over time according
to the equation (minf-mgate)/mtau. Using the rxd.Rate class allows the rxd.State objects
to change over time.
"""
####### WHAT? #######
# Na channel
self.m_gate = rxd.Rate(self.mgate, (self.minf - self.mgate) / self.mtau)
self.h_gate = rxd.Rate(self.hgate, (self.hinf - self.hgate) / self.htau)
# K Channel
# Delayed rectifier
self.n_gate = rxd.Rate(self.ngate, (self.ninf - self.ngate) / self.ntau)
# # A channel
# self.m_gateA = rxd.Rate(self.mgateA, (self.minfA - self.mgateA) / self.mtauA)
# self.h_gateA = rxd.Rate(self.hgateA, (self.hinfA - self.hgateA) / self.htauA)
# H channel
# self.l_gateH = rxd.Rate(self.lgateH, (self.linfH - self.lgateH) / self.ltauH)
# Ca Channel
# m_gate_T = rxd.Rate(mgate_T, (minf_T - mgate_T) / tau_Tm)
# h_gate_T = rxd.Rate(hgate_T, (hinf_T - hgate_T) / tau_Th)
self.m_gate_N = rxd.Rate(self.mgate_N, (self.minf_N - self.mgate_N) / self.tau_Nm)
self.h_gate_N = rxd.Rate(self.hgate_N, (self.hinf_N - self.hgate_N) / self.tau_Nh)
self.m_gate_L = rxd.Rate(self.mgate_L, (self.minf_L - self.mgate_L) / self.tau_Lm)
# Nernst Potentials
self.Ena = 1e3 * h.R * self.T * log(self.nao / self.nai) / h.FARADAY
self.Ek = 1e3 * h.R * self.T * log(self.ko / self.ki) / h.FARADAY
self.Eca = 1e3 * h.R * self.T * log(self.cao / self.cai) / (2 * h.FARADAY)
self.Eh = 0 # mV
# Conductances
self.gNa = inRegion * self.gNaBar * self.mgate ** 3 * self.hgate
self.gK = inRegion * self.gKBar * self.ngate ** 4
# self.gCaT = inRegion * self.gcabar_T * self.mgate_T ** 2 * self.hgate_T
self.gCaN = inRegion * self.gCaBarN * self.mgate_N ** 2 * self.hgate_N
self.gCaL = inRegion * self.gCaBarL * self.mgate_L ** 2
# self.gKA = inRegion * self.gKABar * self.mgateA ** 3 * self.hgateA
# self.gna = self.gnabar * self.mgate ** 3 * self.hgate
# self.gk = self.gkbar * self.ngate ** 4
# # gca_T = gcabar_T * mgate_T ** 2 * hgate_T
# self.gca_N = self.gcabar_N * self.mgate_N ** 2 * self.hgate_N
# self.gca_L = self.gcabar_L * self.mgate_L ** 2
self.dvf = 0.001 / (0.001 + self.cai) * self.temp1 * self.temp3 * (1 - self.cai / self.cao * exp(self.temp2))
self.I_non_mod_L = self.gCaL * self.dvf
self.I_non_mod_N = self.gCaN * self.dvf
# Ca Extrusion by PMCA
self.pmcaBinding = rxd.Reaction(self.ca + self.pmca, self.pmcaCa, kPMCABindF, kPMCABindR,
region=[self.cyt])
self.pmcaRelease = rxd.Reaction(self.pmcaCa, self.pmca, kPMCARel, region=[self.cyt])
# Na/Ca Exchange
self.ncxBinding = rxd.Reaction(2*self.ca + self.ncx, self.ncx2Ca, kNCXBindF, kNCXBindR, region=[self.cyt])
self.ncxRelease = rxd.Reaction(self.ncx2Ca, self.ncx, kNCXRel, region=[self.cyt])
# SERCA
self.sercaBinding = rxd.Reaction(2*self.ca + self.serca, self.serca2Ca, kSERCABindF, kSERCABindR, region=[self.cyt])
self.sercaRelease = rxd.MultiCompartmentReaction(self.serca2Ca[self.cyt], 2*self.caer + self.serca[self.cyt], kSERCARel, membrane=self.erCytMem, membrane_flux=True)
# Calbindin-D28k
# self.h0l0HBind = rxd.Reaction(self.ca + self.cbdh0l0, self.cbdh1l0)
# self.h0l0LBind
# self.h0l1HBind
# self.h0l1LBind
# self.h0l2HBind
# self.h1l0LBind
# self.h1l0HBind
# self.h1l1LBind
# self.h1l1HBind
# self.h1l2HBind
# self.h2l0LBind
# self.h2l1LBind
# Calreticulin
self.calrBuffering = rxd.Reaction(self.caer + self.calr, self.calrCa, kCalrBindF, kCalrBindR, region=self.er)
# Parvalbumin
self.parvBuffering = rxd.Reaction(self.cai + self.parv, self.parvCa, kParvBindF, kParvBindR, region=self.cyt)
'''
rxd.MultiCompartmentReaction tracks current across multiple areas, in this case, cytoplasm and extracellular fluid
'''
# Current through ion channels
self.INa = rxd.MultiCompartmentReaction(self.nai, self.nao, self.gNa * (v - self.Ena), mass_action=False,
membrane=self.mem, membrane_flux=True)
totalIK = (self.gK + self.gKABar) * (v - self.Ek)
self.IK = rxd.MultiCompartmentReaction(self.ki, self.ko, totalIK, mass_action=False,
membrane=self.mem, membrane_flux=True)
# ICaT = rxd.MultiCompartmentReaction(cai, cao, I_non_mod_T, mass_action=False,
# membrane=mem, membrane_flux=True)
self.ICaN = rxd.MultiCompartmentReaction(self.cai, self.cao, self.I_non_mod_N, mass_action=False,
membrane=self.mem, membrane_flux=True)
self.ICaL = rxd.MultiCompartmentReaction(self.cai, self.cao, self.I_non_mod_L, mass_action=False,
membrane=self.mem, membrane_flux=True)
self.ILeak = rxd.MultiCompartmentReaction(self.xi, self.xo, self.gL * (v - self.el), mass_action=False,
membrane=self.mem, membrane_flux=True)
def stimulate(self, stimNeuron, measureNeuron, location, delay, amp, stim_dur, sim_dur):
'''
:param neuron: Section of neuron to place electrode
:param location: Location within section
:param delay: No current from start of simulation
:param amp: Current amplitude in nA
:param stim_dur: Duration of current injection
:param sim_dur: Duration of simulation
:return: NONE
'''
self.stim = h.IClamp(stimNeuron(location))
self.stim.delay = delay
self.stim.amp = amp
self.stim.dur = stim_dur
'''
Following vectors are necessary in order to create graphs.
'''
self.tvec = h.Vector().record(h._ref_t)
self.vvec = h.Vector().record(measureNeuron(location)._ref_v)
# self.mvec = h.Vector().record(self.mgate[self.cyt].nodes(measureNeuron(location))._ref_value)
# self.nvec = h.Vector().record(self.ngate[self.cyt].nodes(measureNeuron(location))._ref_value)
# self.hvec = h.Vector().record(self.hgate[self.cyt].nodes(measureNeuron(location))._ref_value)
self.caL_mvec = h.Vector().record(self.mgate_L[self.cyt].nodes(measureNeuron(location))._ref_value)
self.caN_mvec = h.Vector().record(self.mgate_N[self.cyt].nodes(measureNeuron(location))._ref_value)
self.caN_hvec = h.Vector().record(self.hgate_N[self.cyt].nodes(measureNeuron(location))._ref_value)
#
# self.kivec = h.Vector().record(measureNeuron(location)._ref_ki)
# self.kovec = h.Vector().record(measureNeuron(location)._ref_ko)
# self.naovec = h.Vector().record(measureNeuron(location)._ref_nao)
# self.naivec = h.Vector().record(measureNeuron(location)._ref_nai)
self.caovec = h.Vector().record(measureNeuron(location)._ref_cao)
self.caivec = h.Vector().record(measureNeuron(location)._ref_cai)
self.caervec = h.Vector().record(self.caer.nodes(measureNeuron(location))._ref_concentration)
# self.icavec = h.Vector().record(measureNeuron(location)._ref_ica)
# self.pmcavec = h.Vector().record(self.pmca.nodes(measureNeuron(location))._ref_value)
# self.pmcaCavec = h.Vector().record(self.pmcaCa.nodes(measureNeuron(location))._ref_value)
#
# self.ncxvec = h.Vector().record(self.ncx.nodes(measureNeuron(location))._ref_value)
# self.ncx2cavec = h.Vector().record(self.ncx2Ca.nodes(measureNeuron(location))._ref_value)
self.calrvec = h.Vector().record(self.calr.nodes(measureNeuron(location))._ref_concentration)
self.calrcavec = h.Vector().record(self.calrCa.nodes(measureNeuron(location))._ref_concentration)
# self.sercavec = h.Vector().record(self.serca[self.cyt].nodes(measureNeuron(location))._ref_concentration)
# self.serca2cavec = h.Vector().record(self.serca2Ca[self.cyt].nodes(measureNeuron(location))._ref_concentration)
# h.cvode.active(1)
#
# h.v_init = -70.0
#
# h.tstop = sim_dur
# h.celsius = 34.0
# This line needed:
# h.finitialize(-70 * mV)
fih = h.FInitializeHandler(callable)
fih.allprint()
print(h.hname())
h.continuerun(sim_dur)
def interp(self, vec):
return vec.as_numpy() - 0.5 * np.diff(vec.as_numpy(), prepend=vec.x[0])
def flux(self, vec):
return np.diff(vec.as_numpy(), prepend=vec.x[0]) / h.dt
def rxd_flux_mols(self, reaction):
self.rate_str = repr(-reaction.f_rate) if reaction.b_rate is None else repr(reaction.b_rate - reaction.f_rate)
self.rxd_to_py = {'v': 'self.interp(self.vvec)',
'log': 'np.log',
'exp': 'np.exp',
'tanh': 'np.tanh',
'fabs': 'np.fabs',
'self.k[self.cyt]': 'self.interp(self.kivec)',
'self.k[self.ecs]': 'self.interp(self.kovec)',
'self.na[self.cyt]': 'self.interp(self.naivec)',
'self.na[self.ecs]': 'self.interp(self.naovec)',
'self.ngate': 'self.interp(self.nvec)',
'self.mgate': 'self.interp(self.mvec)',
'self.hgate': 'self.interp(self.hvec)',
'self.mgate_L': 'self.interp(self.ca_Lvec)',
'self.mgate_N': 'self.interp(self.caN_mvec)',
'self.hgate_N': 'self.interp(self.caN_hvec)',
'self.ca[self.cyt]': 'self.interp(self.caivec)',
'self.ca[self.ecs]': 'self.interp(self.caovec)'}
self.current_rate = functools.reduce(lambda r, kv: r.replace(*kv), self.rxd_to_py.items(), self.rate_str)
self.flux_mols = eval(self.current_rate)
return self.flux_mols
def rxd_flux_mM(self, reaction, node):
self.flux_mols = self.rxd_flux_mols(reaction)
return self.flux_mols * node.surface_area / (node.volume * rxd.constants.molecules_per_mM_um3())
if __name__ == "__main__":
rgc1 = rgcMorph.RGCMorphology()
vgcc = VGCC_2(rgc1, 'rxd_file.json')
vgcc.implement_RXD(rgc1)
vgcc.stimulate(vgcc.somaList[0], vgcc.somaList[0], 0.5, 0, 0.5, 60, 100)
# print(vgcc.somaList[0].diam)
fig = pyplot.figure()
pyplot.plot(vgcc.tvec.as_numpy(), vgcc.vvec.as_numpy(), '-b', label='k')
pyplot.xlabel('t (ms)')
pyplot.ylabel('Voltage (mV)')
pyplot.savefig(os.path.join('MemPotTest.svg'), format="svg")
# fig = pyplot.figure()
# pyplot.plot(vgcc.tvec.as_numpy(), vgcc.caL_mvec.as_numpy(), '-b', label='L m-gate')
# # pyplot.plot(vgcc.tvec.as_numpy(), vgcc.caN_hvec.as_numpy(), '-r', label='N h-gate')
# pyplot.plot(vgcc.tvec.as_numpy(), vgcc.caN_mvec.as_numpy(), '-g', label='N m-gate')
# pyplot.legend()
# pyplot.xlabel('t (ms)')
# pyplot.ylabel('Open Probability')
# pyplot.savefig(os.path.join('OpenProbTest.svg'), format="svg")
#
# fig = pyplot.figure()
# pyplot.plot(vgcc.tvec, vgcc.icavec.as_numpy(), '-g', label='k')
# pyplot.xlabel('t (ms)')
# pyplot.ylabel('Current (mA/cm^2)')
# pyplot.savefig(os.path.join('CaCurrentClass.svg'), format="svg")
fig = pyplot.figure()
pyplot.plot(vgcc.tvec.as_numpy(), vgcc.caivec.as_numpy(), '-g', label='k')
pyplot.xlabel('time (ms)')
pyplot.ylabel('Intracellular Ca Concentration (mM)')
pyplot.savefig(os.path.join('IntraCaTest.svg'), format="svg")
fig = pyplot.figure()
pyplot.plot(vgcc.tvec.as_numpy(), vgcc.caervec.as_numpy(), '-g', label='k')
pyplot.xlabel('time (ms)')
pyplot.ylabel('ER Lumenal Ca Concentration (mM)')
pyplot.savefig(os.path.join('LumenCaTest.svg'), format="svg")
# fig = pyplot.figure()
# pyplot.plot(vgcc.tvec, vgcc.pmcavec.as_numpy(), '-b', label='PMCA')
# pyplot.plot(vgcc.tvec, vgcc.pmcaCavec.as_numpy(), '-r', label='PMCA_CA')
# pyplot.xlabel('t (ms)')
# pyplot.ylabel('Concentration')
# pyplot.legend()
# pyplot.savefig(os.path.join('PMCAClass.svg'), format="svg")
# fig = pyplot.figure()
# pyplot.plot(vgcc.tvec, vgcc.ncxvec.as_numpy(), '-b', label='ncx')
# pyplot.plot(vgcc.tvec, vgcc.ncx2cavec.as_numpy(), '-r', label='ncx_2ca')
# pyplot.xlabel('t (ms)')
# pyplot.ylabel('Concentration')
# pyplot.legend()
# pyplot.savefig(os.path.join('NCXClass.svg'), format="svg")
# fig = pyplot.figure()
# pyplot.plot(vgcc.tvec, vgcc.sercavec.as_numpy(), '-b', label='serca')
# pyplot.plot(vgcc.tvec, vgcc.serca2cavec.as_numpy(), '-r', label='serca_2ca')
# pyplot.xlabel('t (ms)')
# pyplot.ylabel('Concentration')
# pyplot.legend()
# pyplot.savefig(os.path.join('SERCA.svg'), format="svg")
# fig = pyplot.figure()
# pyplot.plot(vgcc.tvec, vgcc.calrvec.as_numpy(), '-b', label='calr')
# pyplot.plot(vgcc.tvec, vgcc.calrcavec.as_numpy(), '-r', label='calr_ca')
# pyplot.xlabel('t (ms)')
# pyplot.ylabel('Concentration')
# pyplot.legend()
# pyplot.savefig(os.path.join('calr.svg'), format="svg")
# print(vgcc.caivec[-1])
rgc1 = vgcc = None
Code: Select all
from neuron import h
h.load_file('stdrun.hoc')
h.load_file("import3d.hoc")
h.load_file("stdlib.hoc")
class RGCMorphology:
def __init__(self):
self.load_morphology()
def load_morphology(self):
cell = h.Import3d_SWC_read()
cell.input("8.swc")
i3d = h.Import3d_GUI(cell, False)
i3d.instantiate(self)