Connecting rxd to custom mechanism where custom mechanism models production/secretion

Extending NEURON to handle reaction-diffusion problems.

Moderators: hines, wwlytton, ramcdougal

Post Reply
tiki
Posts: 5
Joined: Sun Jun 28, 2020 3:07 pm

Connecting rxd to custom mechanism where custom mechanism models production/secretion

Post by tiki »

Hello,

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
}
The relevant parts of the mod file above are Sst and JSS where Sst is the internal concentration of somatostatin and JSS is the amount of somatostatin secreted at a given time step.

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 produces the following result
Image
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
or add a rxd.Rate statement, that uptake and release occur. I am not sure however that this uptake and release is according to the mechanism I have in the mod file however. Would the best approach be to write the whole mechanism in the python script (like here)or is there a way to connect release rate and internal concentration to a custom mechanism?
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Connecting rxd to custom mechanism where custom mechanism models production/secretion

Post by ted »

adamjhn or ramcdougal will have to answer your rxd-specific questions.

I have some questions about the NMODL code:

All those variable names and parameter values, but no units. You're sure there are no units inconsistencies? that no scale factors are necessary?

Why
USEION sst READ ssti, ssto WRITE isst
and why bother with sstin and sstout if a value is never computed or assigned for ssti or ssto?
tiki
Posts: 5
Joined: Sun Jun 28, 2020 3:07 pm

Re: Connecting rxd to custom mechanism where custom mechanism models production/secretion

Post by tiki »

Hello Ted,

Should I repost my above post in another section of the forum for adamjhn or ramcdougal?

As for the code, much of it was translated from the xpp files at this link: https://lbm.niddk.nih.gov/sherman/gallery/bad/ and we plan to add units in the near future, but I don't see a reason for inconsistencies outside of a typo.

I used sstin and sstout because ssti and ssto were not output when I ran

Code: Select all

delta.psection()['density_mechs']
It was just a way to store the values of ssti and ssto at each time step in the csv generated by the script.

I am also unsure whether writing isst would changes the value of ssti within and outside of the cell appropriately when there is no production (rxd.Rate) of sst.

Just for reference, here are the other files I used to produce the above figure:
constants.ini

Code: Select all

[Alpha]
eCa
gCaL
gcapq
vcapqm
scapqm
vcapqh
scapqh
tcapqh1
tcapqh2
gCaT
eGIRK
gGIRKko 
sombara2
ssom2
alpha 
Ba 
fcyta 
fVpqa 
tmsb 
vcella
vmdpq 
kpmcaa 
ksercaa 
pleaka 
fera 
sigmava 
fmda 
k1a 
km1a 
r1a 
rm1a 
r20a 
r30a 
rm3a 
u1a 
u2a 
u3a 
kpa 
kp2a 
GlucFacta 
knockoutda 
ra 
sombara 
rako 
ssom 
fa 
vc 
eK
gKa
tau_mKa
knockoutba
ka1
gKATPbara
gkdr
vkdrm
skdrm
eL
gL
eNa
gNa

[Beta]
vmd 
vcell 
B 
kpmca 
cbas 
kserca2b 
kserca3 
per 
phigk 
KGPDH 
kappa 
Jgk 
fcyt 
delta 
p23 
p24 
psim 
p21 
p22 
fer 
sigmav 
fmd 
lambda 
VmaxPFK 
weight1 
topa1 
bottom1 
atot 
k4 
k3 
f43 
k2 
f42 
f23 
amp 
k1 
f41 
f13 
gamma 
p19 
Amtot 
p20 
FRT
p16 
p13 
p14 
p15 
fmito 
p8 
p9 
p10 
p11 
p17 
p18 
Cmito 
p1 
p2 
p3 
NADmtot 
Jgpdh_bas 
p4 
p5
p6 
p7 
vm 
sm 
gca 
nca 
vCa 
raL 
khyd 
JhydSS 
taua 
u1 
u2 
u3 
gthresh 
Kp2 
Kp 
r30 
factor
amplify 
bas_r3 
rm3 
rb 
knockoutdb 
sombarb
ssomb 
r1 
rm1 
km1 
bas_cmd 
max_cmd 
cmdp  
kcmd 
exo_k1 
fb 

[Delta]
sCaLm
sCaLh
vCaLm 
vCaLh
tCaLh1
tCaLh2
gCaPQ 
gCaL
gGABAbar
gKa
vCaPQm 
sCaPQm 
vCaPQh 
sCaPQh 
tCaPQh1 
tCaPQh2 
knockoutbd
vGABA
vKam 
sKam 
vKah 
sKah 
tauKam 
tKah1 
tKah2 
gKATPbar 
vKdrm 
sKdrm 
gKdr
vNam
vNah 
sNam 
sNah 
tNah1 
tNah2 
con
vmdl 
vmdPQ 
fVl 
fVPQ   
kserca 
pleak  
tausom 
bas 
f 
Delta_CaL.mod

Code: Select all

NEURON{
SUFFIX D_CaL
USEION CaL WRITE iCaL VALENCE 2
USEION Ca WRITE eCa VALENCE 2
RANGE sCaLm, sCaLh, vCaLm, vCaLh, tCaLh1, tCaLh2, gCaL
RANGE mCaL_inf, hCaL_inf, tauCaLm, tauCaLh, iCaL
}

PARAMETER{   
sCaLm
sCaLh
vCaLm 
vCaLh
tCaLh1
tCaLh2
gCaL
eCa
v : This is the voltage when I run h.initial.....

}

ASSIGNED{
mCaL_inf
hCaL_inf
tauCaLm
tauCaLh
iCaL
}

STATE{
mCaL
hCaL
}

INITIAL{
eCa = 65 : mV
gCaL = 0.7 : nS
mCaL = 0.8218051702003508
hCaL = 0.6672499701175263 
sCaLm = 10
sCaLh = -5
vCaLm = -30
vCaLh = -33
tCaLh1 = 60
tCaLh2 = 51
}

BREAKPOINT{
mCaL_inf = 1/(1 + exp(-(v - vCaLm)/sCaLm))
hCaL_inf = 1/(1 + exp(-(v - vCaLh)/sCaLh))    
tauCaLm = (1/(exp(-(v + 23)/20) + exp((v + 23)/20))) + 0.05
tauCaLh = (tCaLh1/(exp(-(v + 0)/20) + exp((v + 0)/20)))+tCaLh2
SOLVE states METHOD cnexp :were putting this before current variable    
iCaL = gCaL*pow(mCaL,2)*hCaL*(v-eCa)
}

DERIVATIVE states{
mCaL'= (mCaL_inf - mCaL)/tauCaLm
hCaL'=(hCaL_inf - hCaL)/tauCaLh
}
Delta_CaPQ.mod

Code: Select all

NEURON{
SUFFIX D_CaPQ
USEION CaPQ WRITE iCaPQ VALENCE 2
USEION Ca READ eCa
RANGE gCaPQ, mCaPQ, hCaPQ, vCaPQm, sCaPQm, vCaPQh, sCaPQh, tCaPQh1, tCaPQh2  
RANGE mCaPQ_inf, hCaPQ_inf, tauCaPQm, tauCaPQh, iCaPQ
}

PARAMETER{
gCaPQ  
vCaPQm 
sCaPQm 
vCaPQh 
sCaPQh 
tCaPQh1 
tCaPQh2   
eCa
v : This is the voltage when I run h.initial.....

}

ASSIGNED{
mCaPQ_inf
hCaPQ_inf
tauCaPQm
tauCaPQh
iCaPQ
}

STATE{
mCaPQ
hCaPQ
}

INITIAL{
eCa = 65
gCaPQ = 0.7
mCaPQ = 0.5089033026581684  
hCaPQ = 0.6672499701175263  
vCaPQm = -15
sCaPQm = 10
vCaPQh = -33
sCaPQh = -5
tCaPQh1 = 60
tCaPQh2 = 51
}

BREAKPOINT{
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
iCaPQ = gCaPQ * mCaPQ * hCaPQ * (v - eCa)
}

DERIVATIVE states{
mCaPQ'= (mCaPQ_inf - mCaPQ)/tauCaPQm
hCaPQ' = (hCaPQ_inf - hCaPQ)/tauCaPQh
}
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Connecting rxd to custom mechanism where custom mechanism models production/secretion

Post by ted »

Should I repost
Not at all. They routinely monitor this part of the Forum. I was just declaring the limited scope of my own comments.
I don't see a reason for inconsistencies outside of a typo
I have just seen too much user-written NMODL code. One of the advantages of the rxd module is the convenience with which different mixes of units can be handled. For NMODL code there is a tool called modlunit that can detect and report inconsistencies, but it's not quite as easy to work with.

Now something I just saw:
the USEION statement includes
WRITE isst
and indeed the BREAKPOINT block calculates
isst = tmsb * som * con
As far as I know, isst will contribute to the charge balance equation for any compartment in which it is nonzero, regardless of the VALENCE of sst. This could become a problem if your mechanism ever exists in a model cell that that has a voltage-sensitive process (e.g. an ion channel).

Second, NONSPECIFIC_CURRENT i is declared, but never assigned a value. Did you perhaps intend to follow the calculation of isst with the statement
i = -isst
in order to cancel out the polarizing effect of isst on membrane potential?
tiki
Posts: 5
Joined: Sun Jun 28, 2020 3:07 pm

Re: Connecting rxd to custom mechanism where custom mechanism models production/secretion

Post by tiki »

Thank you, I did not know that a current for an ion with valence 0 would affect the charge balance equation. Do you know what the role of this current is on the concentration of ssti and ssto. The ion concentration seem to change only when the urrate and rrate that are input into the multicompartment reaction are changed.
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Connecting rxd to custom mechanism where custom mechanism models production/secretion

Post by ted »

I did not know that a current for an ion with valence 0 would affect the charge balance equation.
I don't know for sure, would have to run a test to verify, but it seems reasonable. As far as I know, valence is used only by NEURON's computational engine when the ion_style settings tell NEURON to calculate reversal potential from the Nernst equation.
Do you know what the role of this current is on the concentration of ssti and ssto.
It would do nothing unless there is some NMODL-specified mechanism that READs isst and WRITEs ssti or ssto (the NMODL code you've shown here doesn't do that), or some bit of Python that makes the rxd module calculate the values of ssti or ssto from relevant fluxes and volumes of distribution. I don't know enough about rxd to tell you if the latter is the case.
tiki
Posts: 5
Joined: Sun Jun 28, 2020 3:07 pm

Re: Connecting rxd to custom mechanism where custom mechanism models production/secretion

Post by tiki »

Right now isst is set to the flux of sst, which I felt would be more appropriately represented by ssti' in the derivative block. Unfortunately, writing ssti in this block did not change the above result.
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Connecting rxd to custom mechanism where custom mechanism models production/secretion

Post by ted »

tiki wrote:now isst is set to the flux of sst . . . writing ssti in this block did not change the above result.
Another couple of misunderstandings appear to have emerged. First, the trivial one. NMODL expects currents to be in units of absolute current (nA) or current density (mA/cm2), depending on whether the mechanism in question is a point process or a density mechanism. Fluxes are generally in units of millimoles per ms or millimoles per (ms * cm2), so conversion between flux and current involves Faraday's constant.

Now for the nontrivial misunderstanding. I should have explicitly said I was only explaining why calculating a value for isst in the NMODL file had no effect on ssti or ssto--not actually suggesting that you do anything about that in the NMODL file. Instead, use the rxd module's methods, as you were originally trying to do in the Python file you first posted in this thread.

At this point I must pass the baton to someone else because I am not finding documentation of many of the parameters used in your Python cod.
adamjhn
Posts: 54
Joined: Tue Apr 18, 2017 10:05 am

Re: Connecting rxd to custom mechanism where custom mechanism models production/secretion

Post by adamjhn »

Sorry for the delay, you're right about not using charge=0. Rxd reads the currents and then (if they have a charge or valence) divides it by the charge and scale it appropriately to determine the change in concentration. Unfortunately this does not work with charge=0.

As you have found, to work around this constraint, change “charge=0” in your python file and “VALENCE 0” in the mod file to 1. Then, as Ted suggested, you can subtract the current isst in the mod file from the NONSPECIFIC_CURRENT i to prevent the current from altering the voltage.

If you decide to use a point process rather than a range process, you can use a pointer to link values in the mod file with those in rxd, as shown in the mGlur example.
Post Reply