USEION in mod files for RxD species I created

Extending NEURON to handle reaction-diffusion problems.

Moderators: hines, wwlytton, ramcdougal

Post Reply
bschneiders
Posts: 34
Joined: Thu Feb 02, 2017 11:30 am

USEION in mod files for RxD species I created

Post by bschneiders »

Hi, I am trying to keep track of the maximum value of several concentrations using mod files (so that I don't have to record in a bunch of locations). This worked really well for calcium, but the same syntax doesn't work for species that I created with RxD. Here is what worked for calcium:

Code: Select all

NEURON {
    SUFFIX camax
    USEION ca READ cai WRITE cai
    RANGE cam
}

ASSIGNED {
       ca (mM)
       cam (mM)
}

PARAMETER {
	cai		(mM)
}

INITIAL {
    cam = cai
}

BREAKPOINT { 
   if (cai>cam) { cam=cai }
}
and then "sec.insert('camax')" in my python code for all sections. However, when I try this for species "CG" (and associated "CGi", "cgmax"), where:

Code: Select all

cg = rxd.Species(r, initial=cg0, name='CG', d=dg, charge=0)
it doesn't work. Any idea why or what to change? My guess is it has something to do with calcium already being recognized by Neuron, and my new species not, but I'm not sure how to change that. Thanks!
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: USEION in mod files for RxD species I created

Post by ted »

Here's a couple of minor suggestions for your "max cai" mod file--

It doesn't need to WRITE cai.
It's OK to declare cai in the PARAMETER block, but strictly speaking it's not a parameter--it's a variable whose value comes from NEURON's computational engine, so it would be conceptually more consistent with being declared in the ASSIGNED block.

Unfortunately, neither of these suggestions has anything to do with your question about accessing the values of RxD variables from NMODL.

Now, for your question. One workaround is to use Vector.record to capture the time course of the species of interest, then after the simulation find its maximum with Vector.max. But instead, you might find it better to try this:
In model setup, make sure your "max conc" "mechanism" is inserted before you execute any RxD code. The USEION statement in its NMODL code makes NEURON create a new species with whatever name and charge you specified. When your model setup code finally gets to the point where you're ready to execute RxD statements, use neuron.rxd.Species's name and charge parameters (see https://www.neuron.yale.edu/neuron/stat ... xd.Species) to set the name that RxD will use to refer to the species, and the charge associated with that species. Do this before any statements that instantiate the reaction-diffusion model. This should ensure that the values calculated by RxD will be reflected in the value of the non-RxD variable that your NMODL file created.
ramcdougal
Posts: 267
Joined: Fri Nov 28, 2008 3:38 pm
Location: Yale School of Public Health

Re: USEION in mod files for RxD species I created

Post by ramcdougal »

Ted raises good points.

The underlying issue though is that if you run the straightforward translation of that mod file, you'll get a maximum value of 1 (assuming your concentrations are below 1). Why 1? It has to do with the order of events during initialization: NEURON sets concentrations to 1 unless there's a (in your case) h.CGi0_CG_ion specifying a different initial concentration, then does MOD file INITIAL rules, and finally does rxd initialization rules. In particular, that means that during initialization, the max variable is set to 1 and can only go up from there.

There are at least two reasonable ways around this:
  • Set h.CGi0_CG_ion = cg0
  • Change your initial block to e.g. max = 0
A full working example follows:

cgmax.mod

Code: Select all

NEURON {
    SUFFIX cgmax
    USEION CG READ CGi CHARGE 0
    RANGE max
}

ASSIGNED {
       CGi (mM)
       max (mM)
}

INITIAL {
    max = CGi
}

BREAKPOINT { 
   if (CGi>max) { max=CGi }
}
test.py

Code: Select all

from neuron import h, rxd

cg0 = 0.9
dg = 1

soma = h.Section(name='soma')
soma.insert('cgmax')
r = rxd.Region([soma], nrn_region='i')
h.CGi0_CG_ion = cg0
cg = rxd.Species(r, initial=cg0, name='CG', d=dg, charge=0)
dec = rxd.Rate(cg, 0.1)
cgvec = h.Vector()
cgvec.record(soma(0.5)._ref_CGi)
t = h.Vector()
t.record(h._ref_t)

h.finitialize()
h.fadvance()
h.fadvance()
print(soma(0.5).CGi)
print(max(cgvec))
print(soma(0.5).cgmax.max)
One important thing you'll notice if you run the above: The value of cgmax.max is the maximum from t=0 until the beginning of the last time step. Why? That's when BREAKPOINT is executed, as it has to be, because you can't compute the new e.g. membrane potential without computing the currents first.
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: USEION in mod files for RxD species I created

Post by ted »

The value of cgmax.max is the maximum from t=0 until the beginning of the last time step.
If you don't want that, change

Code: Select all

BREAKPOINT { 
   if (CGi>max) { max=CGi }
}
to

Code: Select all

AFTER BREAKPOINT { 
   if (CGi>max) { max=CGi }
}
and the maximum value will be determined from the values AFTER each time step.
Post Reply