Overriding _w_currents() in rxd.py

Extending NEURON to handle reaction-diffusion problems.

Moderators: hines, wwlytton, ramcdougal

Post Reply
davidsterratt
Posts: 26
Joined: Tue Jan 17, 2006 11:36 am

Overriding _w_currents() in rxd.py

Post by davidsterratt » Tue Jul 08, 2014 11:43 am

I'm trying to override the _w_currents() callback in rxd.py so that I can add a contribution to the membrane current generated by a rule-based model that I've embedded using a subclass of a GeneralizedReaction.

The following snippet works for me - almost:

Code: Select all

import neuron.rxd.rxd as nrr
import random

def _kn_currents(rhs):
    nrr._currents(rhs)
    # global nrr._rxd_induced_currents
    print "adding some noise"
    sign = 1
    cur = random.random()/10
    ## This line alters ica, but does not seem to have effect on voltage
    print nrr._curr_ptrs[0][0]
    nrr._curr_ptrs[0][0] += -sign * cur
    print nrr._curr_ptrs[0][0]
    # print nrr._rxd_induced_currents
    # nrr._rxd_induced_currents[0] += 0.0
    nrr._rxd_induced_currents[0] += sign * cur
    print nrr._rxd_induced_currents

nrr._callbacks[2] = _kn_currents
When I run my simulation the above code changes ica (as evidenced by a recording of _ref_ica, in which the noise is evident), but the noisy current doesn't seem to feed back into noise in the voltage (_ref_v). There is a mod file in my simulation which writes ica, and this I know that the ica that is definitely feeding into the simulation, but the noisy component doesn't seem to be. It's as though the changes in _currents are being ignored by the voltage integrator altogether.

Any thoughts on this would be welcome. I realise that I'm probably missing a trick by not using the apparatus in multiCompartmentalReaction.py, but this is my attempt to understand how the code is working without losing myself completely in all the pointer maps!

davidsterratt
Posts: 26
Joined: Tue Jan 17, 2006 11:36 am

Re: Overriding _w_currents() in rxd.py

Post by davidsterratt » Tue Jul 08, 2014 6:36 pm

I've now got a minimal complete (if nonsensical) example that shows my problem more clearly. There's also a possible solution, which I got to by examining the NEURON source code, in particular by doing some printf debugging in setup_tree_matrix() in treeset.c. It seems that the rhs argument feeds into the solver directly, but any species-specific currents (e.g. ica) are not incorporated as components of the total current before the solver is called. So setting the relevant element of the current pointer neuron.rxd.rxd._curr_ptrs and the element of rhs corresponding to the same segment is necessary to see current injected by the callback routine manifest in the current, voltage and concentration. However, I'm not sure if this is the intended mode of operation

I should add that I'm using python2.7 and the version of NEURON from the default mercurial branch.

Code: Select all

from neuron import *
from neuron import rxd
import neuron.rxd.rxd as nrr
import random
import matplotlib.pyplot as plt

def _kn_currents(rhs):
    nrr._currents(rhs)
    cur = -random.random()/10
    ## This line alters ica, but on its own does not seem to affect
    ## the voltage
    nrr._curr_ptrs[0][0] += cur
    ## This line is necessary to change the voltage
    rhs[1] += -cur

nrr._callbacks[2] = _kn_currents

sh = h.Section()
sh.insert("pas")
r = rxd.Region([sh], nrn_region='i')
ca = rxd.Species(r, name='ca', charge=2, initial=0.0)

rec_t = h.Vector()
rec_t.record(h._ref_t)
rec_v = h.Vector()
rec_v.record(sh(0.5)._ref_v)
rec_ica = h.Vector()
rec_ica.record(sh(0.5)._ref_ica)
rec_cai = h.Vector()
rec_cai.record(sh(0.5)._ref_cai)

init()
run(15)

fig, ax = plt.subplots(nrows=3, ncols=1, figsize=(2.25*2, 2.5*2))
plt.subplots_adjust(left=0.2, top=0.95, bottom=0.1)

ax[0].plot(rec_t, rec_v)
ax[0].set_xlabel("")
ax[0].set_ylabel("V (mV)")
ax[0].axis(ymin=-80, ymax=50)

ax[1].plot(rec_t, rec_ica)
ax[1].set_xlabel("")
ax[1].set_ylabel("ICa (mA/cm2)")
ax[1].axis(ymin=-0.1, ymax=0.1)

ax[2].plot(rec_t, rec_cai)
ax[2].set_xlabel("")
ax[2].set_ylabel("[Ca] (mM)")
ax[2].axis(ymin=0, ymax=0.001)

fig.show()

ramcdougal
Posts: 122
Joined: Fri Nov 28, 2008 3:38 pm
Location: Yale School of Medicine

Re: Overriding _w_currents() in rxd.py

Post by ramcdougal » Wed Jul 09, 2014 5:32 pm

You are correct; when _currents(rhs) is invoked, rhs already contains the collected currents (ica + ina + ik + etc...), so the current and the voltage derivative (contained in rhs) both must be changed separately. This is because nrn_nonvint_block_current is invoked after nrn_rhs in src/nrnoc/treeset.c's setup_tree_matrix and after rhs_memb in src/nrncvode/cvtrset.cpp's Cvode::rhs. I'm not sure if those lines can be safely swapped (with the necessary changes in rxd), but even if they can that would cause your code to behave differently in 7.3 vs 7.4. Do we think fixing this awkwardness is worth the resulting incompatibilities?

The other challenge is to identify what index in rhs corresponds to your segment of interest. In NEURON 7.3+, the index in rhs for the voltage derivative for a segment seg is seg.node_index().

davidsterratt
Posts: 26
Joined: Tue Jan 17, 2006 11:36 am

Re: Overriding _w_currents() in rxd.py

Post by davidsterratt » Sun Jul 13, 2014 4:58 pm

Thanks for the clarification and the tip. Re the suggestion about changing the order of calling the blocks, I'm not sure. Intuitively I expected the callback to be like WRITE ica I'm a mod file, but it may actually suit my purposes as well if it doesn't.

davidsterratt
Posts: 26
Joined: Tue Jan 17, 2006 11:36 am

Re: Overriding _w_currents() in rxd.py

Post by davidsterratt » Sun Jul 13, 2014 4:59 pm

Thanks for the clarification and the tip. Re the suggestion about changing the order of calling the blocks, I'm not sure. Intuitively I expected the callback to be like WRITE ica I'm a mod file, but it may actually suit my purposes as well if it doesn't. I'll get back to this once I'm back in the office.

Post Reply