Reaction-Diffusion in HOC

Extending NEURON to handle reaction-diffusion problems.

Moderators: hines, wwlytton, ramcdougal

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

Reaction-Diffusion in HOC

Post by ramcdougal » Fri Oct 27, 2017 10:09 am

The reaction-diffusion module can be used in HOC almost identically to how it is used in Python if you begin your HOC code with:

Code: Select all

objref pyobj, h, rxd

{
    // load reaction-diffusion support and get convenient handles
    nrnpython("from neuron import h, rxd")
    pyobj = new PythonObject()
    rxd = pyobj.rxd
    h = pyobj.h
}
This provides easy access to NEURON's rxd and h Python modules. You might not need to use h since you're working in HOC, but it does provide certain convenient functions like h.allsec(), which returns an iterable of all sections usable with rxd without you having to explicitly construct a SectionList.

The main gotchas of using rxd in HOC is that (1) rxd in Python uses operator overloading to specify reactants and products; in HOC, you have to use __add__, etc instead. (2) rxd in Python is usually used with keyword arguments; in HOC, these must be specified using positional notation.

Here's a full working example that simulates a calcium buffering reaction: Ca + Buf <> CaBuf:

Code: Select all

objref pyobj, h, rxd, cyt, ca, buf, cabuf, buffering, g

{
    // load reaction-diffusion support and get convenient handles
    nrnpython("from neuron import h, rxd")
    pyobj = new PythonObject()
    rxd = pyobj.rxd
    h = pyobj.h
}

{
    // define the domain and the dynamics
    create soma
    
    cyt = rxd.Region(h.allsec(), "i")
    ca = rxd.Species(cyt, 0, "ca", 2, 1)
    buf = rxd.Species(cyt, 0, "buf", 0, 1)
    cabuf = rxd.Species(cyt, 0, "cabuf", 0, 0)

    buffering = rxd.Reaction(ca.__add__(buf), cabuf, 1, 0.1)
}

{
    // if launched with nrniv, we need this to get graph to update automatically
    // and to use run()
    load_file("stdrun.hoc")
}

{
    // define the graph
    g = new Graph()
    g.addvar("ca", &soma.cai(0.5), 1, 1)
    g.addvar("cabuf", &soma.cabufi(0.5), 2, 1)
    g.size(0, 10, 0, 1)
    graphList[0].append(g)
}

{
    // run the simulation
    tstop = 20
    run()
}

Post Reply