I was looking through the past posts on this forum and I felt that my case might not neatly fit into the categories covered by previous posts. I am working on a model of endocrine cell and currently have a custom mechanisms which models the secretion of insulin as follows:
Code: Select all
NEURON{
SUFFIX D_Somatostatin
USEION sst READ ssti, ssto WRITE isst VALENCE 0
NONSPECIFIC_CURRENT i
USEION CaL READ iCaL
USEION CaPQ READ iCaPQ
RANGE isst, ssti, ssto, iCaL, iCaPQ, tmsb, con, alpha, vmdl, vmdPQ, fVl, B, fVPQ, kpmca, kserca, pleak, vCaPQm, sCaPQm, vCaPQh, sCaPQh, tCaPQh1, tCaPQh2, tausom, bas, fcyt, fmd, fer, sigmav, vc, f, sst_init
RANGE JL, JPQ, Jserca, Jer, mCaPQ_inf, hCaPQ_inf, tauCaPQm, tauCaPQh, Jmem, y, Jleak, som
RANGE sstin, sstout
}
PARAMETER{
tmsb
con
alpha
vmdl
vmdPQ
fVl
B
fVPQ
kpmca
kserca
pleak
vCaPQm
sCaPQm
vCaPQh
sCaPQh
tCaPQh1
tCaPQh2
tausom
bas
fcyt
fmd
fer
sigmav
vc
f
}
ASSIGNED{
sstin
sstout
ssti
ssto
isst
i
iCaL
iCaPQ
JL
JPQ
Jserca
Jer
mCaPQ_inf
hCaPQ_inf
tauCaPQm
tauCaPQh
Jmem
y
Jleak
som
JSS
v : This is the voltage when I run h.initial.....
}
STATE{
mCaPQ
hCaPQ
c
cmdl
cmdPQ
cer
Sst
}
INITIAL{
tmsb = 0.001
con = 0.00000000594
alpha = 5.18e-15
vmdl = 2.12e-15
vmdPQ = 1.41e-15
fVl = 0.00340
B = 1
cmdl = 19.82903375122306
c = 0.3051356309084454
fVPQ = 0.00226
cmdPQ = 27.93917378868966
kpmca = 0.3
kserca = 0.4
pleak = 0.0003
cer = 413.8135591677398
mCaPQ = 0.5089033026581684
hCaPQ = 0.6672499701175263
vCaPQm = -15
sCaPQm = 10https://imgur.com/a/b9seZWU
vCaPQh = -33
sCaPQh = -5
tCaPQh1 = 60
tCaPQh2 = 51
tausom = 2000
bas = 0.0009
fcyt = 0.01
fmd = 0.01
fer = 0.01
sigmav = 31
vc = 1e-13
f = 0.003
Sst = 18.71318922819339
}
BREAKPOINT{
sstin = ssti
sstout = ssto
JL = -alpha * iCaL/vmdl
JPQ = -alpha * iCaPQ/vmdPQ
Jserca = kserca * c
Jer = (Jleak - Jserca)
mCaPQ_inf = 1/(1 + exp(-(v - vCaPQm)/sCaPQm))
hCaPQ_inf = 1/(1 + exp(-(v - vCaPQh)/sCaPQh))
tauCaPQm = (1/(exp(-(v + 23)/20) + exp((v + 23)/20))) + 0.05
tauCaPQh = (tCaPQh1/(exp(-(v + 0)/20) + exp((v + 0)/20))) + tCaPQh2
SOLVE states METHOD cnexp : Put current equation below this
Jmem = fVl * B * (cmdl - c) + fVPQ * B * (cmdPQ - c) - kpmca * c
y = pow(cmdPQ/200,4)/(pow(0.2,4) + pow(cmdPQ/200,4))
Jleak = pleak * (cer - c)
som = (200 * mCaPQ * hCaPQ * y/tausom) + bas
JSS = tmsb * som * con
isst = tmsb * som * con
}
DERIVATIVE states{
mCaPQ' = (mCaPQ_inf - mCaPQ)/tauCaPQm
hCaPQ' = (hCaPQ_inf - hCaPQ)/tauCaPQh
c' = fcyt * (Jmem + Jer)
cmdl' = fmd * JL - fmd * B * (cmdl - c)
cmdPQ' = fmd * JPQ - fmd * B * (cmdPQ-c)
cer' = -fer * sigmav * Jer
Sst' = JSS/vc - f * Sst
}
I am inserting the above mechanism into a section with the following script:
Code: Select all
from neuron import h, rxd
import sys
import csv
import numpy as np
from matplotlib import pyplot as plt
import configparser
h.load_file('stdrun.hoc')
# exclude constants from output
configc = configparser.ConfigParser(allow_no_value= True)
configc.optionxform = str
configc.read('constants.ini')
# create a cell
delta = h.Section('delta')
delta.pt3dclear()
delta.pt3dadd(-9,0,0,10)
delta.pt3dadd(1,0,0,10)
delta.insert('D_Somatostatin')
delta.insert('D_CaL')
delta.insert('D_CaPQ')
# Where?
# the intracellular spaces
cyt = rxd.Region(h.allsec(), name='cyt', nrn_region='i')
# plasma membrane
mem = rxd.Region(h.allsec(), name='mem', geometry=rxd.membrane)
# the extracellular space
ecs = rxd.Extracellular(-35, -15, -15, 35, 15, 15, dx=10, volume_fraction=0.2, tortuosity=1.6)
Sst = rxd.Species([cyt, ecs], name='sst', charge=0, d=1.0, initial= 18.71318922819339)
Sst_cyt = Sst[cyt]
Sst_ecs = Sst[ecs]
# interaction between intracellular and extracellular glucagon
R = 1e4
U = 1e4
rrate = R*Sst_cyt # release rate [molecules per square micron per ms]
urate = U*Sst_ecs # uptake rate [molecules per square micron per ms]
sst_release = rxd.MultiCompartmentReaction(Sst_cyt, Sst_ecs, rrate, urate,
membrane=mem,
custom_dynamics=True)
# record the concentrations in the cells
t_vec = h.Vector()
t_vec.record(h._ref_t)
delta_Sst = h.Vector()
delta_Sst.record(delta(0.5)._ref_ssti)
delta_Sst_ecs = h.Vector().record(Sst_ecs.node_by_location(0,0,0)._ref_concentration)
# variables to store data for delta cell
t = []
v = []
rec = {}
header = []
# record mechanisms delta
for i in delta.psection()['density_mechs']:
for j in delta.psection()['density_mechs'][i]:
if j not in configc['Delta']:
header.append(j+'_'+i)
rec[str(j+'_'+i)] = []
# record variables of every mechanism in every segment
for k in delta:
v.append(h.Vector().record(k._ref_v))
mechRecord = getattr(k, '_ref_'+j+'_'+i)
rec[str(j+'_'+i)].append(h.Vector().record(mechRecord))
# fix header / record voltage of every segment
head = ['Time']
count = 0
for i in delta:
temp = 'VC'+str(count)
# ease writing to csv by keeping same format even though it is not necessary
rec[temp] = []
rec[temp].append(h.Vector().record(i._ref_v))
count += 1
head.extend(header)
head.append(temp)
t = h.Vector().record(h._ref_t)
h.finitialize(-65)
h.continuerun(1000)
# write data cell for delta cell
with open('data/'+sys.argv[1],'w') as file:
writer = csv.writer(file,quoting = csv.QUOTE_NONE,escapechar=' ')
writer.writerow(head)
for i in range(len(t)):
out = [t[i]]
for q in rec:
out.append(rec[q][0][i])
writer.writerow(out)
plt.plot(t_vec, delta_Sst, label='delta')
plt.plot(t_vec, delta_Sst_ecs, label='ecs')
plt.legend(frameon=False)
plt.xlabel('t (ms)')
plt.ylabel('sst (mM)')
plt.show()
Which is not what I want since there is no uptake or release of sst. I have found when I change the 'R' and 'U' variables
Code: Select all
R = 1e4
U = 1e4