Changing RxD species concentration from a point process using USEION

[wthun]

I've been using RxD with a point process changing the concentration of a RxD species which is then read by a density mechanism.

Small example:

example.py

from neuron import h
import neuron.crxd as rxd

import numpy as np

import plotly.graph_objects as go
from plotly.subplots import make_subplots
import plotly.io as pio

pio.templates.default = "seaborn"  # big deal


if __name__ == "__main__":
    # basic setup
    h.load_file("stdrun.hoc")
    soma = h.Section(name="soma")
    soma.L = soma.diam = 10
    soma.Ra = 100
    soma.cm = 1

    soma.insert("kirDA")
    soma.insert("pas")

    # iclamp
    iclamp = h.IClamp(0, soma(0.5))
    iclamp.dur = 2e3
    iclamp.delay = 500
    iclamp.amp = 0.01

    for seg in soma:
        seg.pas.g = 1e-4
        seg.pas.e = -65
        seg.kirDA.gbar = 1e-2

    # rxd
    soma_region = rxd.Region([soma], nrn_region="i")

    DA = rxd.Species(soma_region, charge=0, d=1, name="DA", initial=0)
    DA_decay = rxd.Rate(DA, -0.01 * DA)

    # insert point process
    da_syn = h.DASyn(soma(0.5))
    da_syn.quanta = 1e-1

    net_stim = h.NetStim()
    net_stim.number = 10
    net_stim.start = 1000
    net_stim.interval = 100

    nc = h.NetCon(net_stim, da_syn)
    nc.weight[0] = 1

    # this is the part I'd like to replace using USEION
    h.setpointer(DA.nodes[0]._ref_concentration, "DA_conc", da_syn)

    # run and record
    recordings = {}
    recordings["time"] = h.Vector().record(h._ref_t)
    recordings["voltage"] = h.Vector().record(soma(0.5)._ref_v)
    recordings["kirDA-levelDA"] = h.Vector().record(soma(0.5).kirDA._ref_levelDA)
    recordings["kirDA-p_open"] = h.Vector().record(soma(0.5).kirDA._ref_p_open)
    recordings["DA-conc"] = h.Vector().record(DA.nodes[0]._ref_concentration)
    recordings["DA-conc"] = h.Vector().record(DA.nodes[0]._ref_concentration)
    spike_times = h.Vector()
    nc.record(spike_times)
    recordings["spike_times"] = spike_times

    h.tstop = 3e3
    h.run()

    for key, val in recordings.items():
        recordings[key] = val.as_numpy().copy()

    # plots
    fig = make_subplots(
        rows=3,
        cols=1,
        row_heights=[0.2, 0.2, 0.6],
        shared_xaxes=True,
        vertical_spacing=0.02,
    )

    fig.add_trace(
        go.Scatter(
            x=recordings["time"], y=recordings["voltage"], name="membrane voltage"
        ),
        row=3,
        col=1,
    )
    fig.add_trace(
        go.Scatter(
            x=recordings["time"],
            y=recordings["kirDA-levelDA"],
            name="kirDA-levelDA",
            opacity=0.75,
        ),
        row=2,
        col=1,
    )
    fig.add_trace(
        go.Scatter(
            x=recordings["time"],
            y=recordings["DA-conc"],
            name="DA-conc",
            line=dict(dash="dash"),
            opacity=0.75,
        ),
        row=2,
        col=1,
    )
    fig.add_trace(
        go.Scatter(
            x=recordings["spike_times"],
            y=np.ones(len(recordings["spike_times"])),
            name="spike times",
            mode="markers",
        ),
        row=1,
        col=1,
    )
    fig.show()

    import pdb

    pdb.set_trace()

kirDA.mod

NEURON {
    SUFFIX kirDA
    USEION k READ ek WRITE ik
    USEION DA READ DAi VALENCE 0
    RANGE levelDA
    RANGE gbar, gk, ik, p_open
}

UNITS {
    (S) = (siemens)
    (mV) = (millivolt)
    (mA) = (milliamp)
}

PARAMETER {
    gbar = 0.0 	(S/cm2) 
}

ASSIGNED {
    v (mV)
    ek (mV)
    ik (mA/cm2)
    gk (S/cm2)
    minf
    mtau (ms)
    p_open (1)

    : modulation example
    DAi (mM)
    levelDA (mM)
}

STATE { m }

BREAKPOINT {
    SOLVE states METHOD cnexp
    levelDA = DAi
    p_open = m*levelDA
    gk = gbar*p_open
    ik = gk*(v-ek)
}

DERIVATIVE states {
    rates()
    m' = (minf-m)/mtau
}

INITIAL {
    rates()
    m = minf
}

PROCEDURE rates() {
    UNITSOFF
    minf = 1/(1+exp((v-(-82))/13))
    mtau = 1/(exp((v-(-103))/(-14.5))+0.125/(1+exp((v-(-35))/(-19))))
    UNITSON
}

DASyn.mod

NEURON {
    POINT_PROCESS DASyn
    USEION DA READ DA i: WRITE DAi VALENCE 0
    RANGE quanta, tau, open
    POINTER DA_conc
}

PARAMETER {
    quanta = 1e-4 (mM)
    tau = 10 (ms)
}

INITIAL {
    open = 0
}

ASSIGNED {
    DAi (mM)
    DA_conc 
}

STATE { 
    : DAi (mM) 
    open (1)
}

BREAKPOINT {
    SOLVE state METHOD cnexp
}

DERIVATIVE state { 
    open' = -open/tau
    : DAi' = open*quanta
}

NET_RECEIVE(weight) {
    open = open + weight
    DA_conc = DA_conc + quanta
}

Instead of tracking the concentration with a pointer, I'd like to use USEION to write the change in concentration from the synapse (since I think it'd be easier to track interaction between several NT systems, pathways, and sources that way).

If I try to write DAi from the point process directly instead of using a pointer, the RxD concentration doesn't seem to react. When printed, it seems to reset between each stimulation event. Is there some way I could use USEION still?

[Hjorthmedh]

@ramcdougal

[wthun]

I've found out how to do it now using include_flux: DA.nodes[0].include_flux(da_syn._ref_open), a large NetCon weight tis needed to see effects in the example (nc.weight[0] = 10_000_000.0). (So USEION is not necessary).

[ramcdougal]

With the include_flux, the default units are in molecules/ms, which is probably part of why you needed a large weight.

The docstring is probably more helpful than the documentation on this:

nrn/share/lib/python/neuron/rxd/node.py

Line 268 in 0e781ec

def include_flux(self, *args, **kwargs):