## Issue with larger and more equations in Reaction

Extending NEURON to handle reaction-diffusion problems.

Moderators: hines, wwlytton, ramcdougal

jetlim
Posts: 1
Joined: Wed Oct 21, 2020 9:53 am

### Issue with larger and more equations in Reaction

Hello everybody,

I am currently working on a Python rxd simulation where we are simulating multiple mitochondria being transported within an axon, while also keeping track of their molecular outputs (e.g. ATP production and calcium buffering) in both the spatial (1D) and temporal dimension.
I have two questions I need help in.

Question 1: My first question is how I should best implement such a model. Essentially I have one large region (axon), with several moving compartments within it (mitochondria), where it is important to keep track of not only their internal parameters individually, but also their effect on the axon is localized to their position.

My current implementation is to divide the axon Region into Sections of fixed lengths. Then I had defined internal molecular Species (e.g. mitochondrial ATP, Ca) for each individual mitochondria using Species(....,initial = function(node),...),using the function to define the Species as having a value at the mitochondrial position node and 0 everywhere else - essentially my way of saying "all the stuff is here". Then during my simulation, if it detects that a mitochondria's location has resulted in it changing nodes, it shifts the Species value within the array where appropriate.
As for the Reaction Rates (e.g. ATP and Ca transport), I had to define Rates for each individual mitochondria (because they are often dependent on the mitochondria's internal Species concentration). I also used an rxd.Parameter which marks the active node as 1 (where the mitochondria is) and the rest 0, so I can multiply that with the Rate and thus make it only work at that specific node.

While this implementation seems to work fine (for the most part, my 2nd question regards a more specific issue regarding setup), my gut instinct tells me that there could be a better way to do this.

Question 2: My second, more specific question is about setting up the Rates.

Our mitochondria molecular rates are based on papers and we had recently discovered a new paper that showed issues in the simplified equations that we had been using thus far.
However, implementing the new equations resulted in a rather bulky reaction rates. For example:

Code: Select all

# Calcium Reaction Rates, From MK (Bertram), verified,unit: μM/ms
calcrate = []

for i in range(Mito_num):
ratedict = {}
ratedict["u"] = cac/0.2
ratedict["G(u)"] = ratedict["u"]*((1+0.52*ratedict["u"])**2.8)*((1+0.01*ratedict["u"])**3)/(110+((1+0.52*ratedict["u"])**2.8)*(1+0.01*ratedict["u"])**4)
ratedict["w"] = mitolist[i]["PSIm"]/91

ratedict["J_uni_MK"] = ratedict["G(u)"]*0.11*6.73*(ratedict["w"]-1)/(1-rxdmath.exp(-6.73*ratedict["w"]-1))
#new equation from Siqueria2013

ratedict["J_NCX_MK"] = 0.06*(mitolist[i]["cam"]/cac)*rxdmath.exp(0.018*(mitolist[i]["PSIm"]-91))/((1.94**3)*(1+0.375/mitolist[i]["cam"]))

ratedict["Uni_in"] = rxd.Rate(mitolist[i]["cam"], g_ca*ratedict["J_uni_MK"]*mitolist[i]["node_act"])
ratedict["NCX_out"] = rxd.Rate(mitolist[i]["cam"], -g_ca*ratedict["J_NCX_MK"]*mitolist[i]["node_act"])
calcrate.append(ratedict)
When I attempt to run this block of code, it results in a significantly long wait time, such as several minutes for 4-5 mitochondria. Strangely, it does not appear to influence the simulation time, just the setup time.

Based on Spyder profiler, the problem appears to arise from the following check in rate.py, line 63.

Code: Select all

self._voltage_dependent = False if not hasattr(rate,'_voltage_dependent') else rate._voltage_dependent

I cannot seem to find documentation on what this does, but from reading the code, it appears that it checks every single entity within the Rate equation, seeking out for the presence of the membrane potential Vm object. Other than that, I am unsure what this does.
Given the bulk of my equation (probably not helped by having individual rates for individual mitochondria), it results in this check being made a huge number of times. If I change the code in rate.py to:

Code: Select all

self._voltage_dependent = False

The problem is fixed (reducing the setup time to milliseconds), but I am obviously nervous to alter the library code itself without knowing what it does or what alternative is available.

Posts: 43
Joined: Tue Apr 18, 2017 10:05 am

### Re: Issue with larger and more equations in Reaction

1. This is a great question, your workaround sounds reasonable. I don’t think it’s necessary to split the axon into multiple sections, just set nseg as required. An alternative solution is to use node.include_flux, here is a toy example using an nmodl mechanism to allow a leak of ATP between a mitochondria and the axon.

Code: Select all

import numpy
from neuron import h, rxd
from matplotlib import pyplot
cv = h.CVode()

axon = h.Section(name='axon')
axon.L = 100
axon.diam = 2
axon.nseg = 11

mito = h.Section(name='mitochondria')
mito.L = 1
mito.diam = 1
mito.nseg = 1

# create the regions
cyt = rxd.Region([axon], name='cyt', nrn_region='i')
mit = rxd.Region([mito], name='mit', )

# define the species
ATP = rxd.Species([cyt, mit], name='ATP', d=1,
initial=lambda nd: 10 if nd.region == mit else 0)

# insert the mechanism that links the mitochondria and axon
# set the pointer used in the mechanisms to the relevant state/species in the
# mitochondria
axon.insert("atpleak")
for seg in axon:
seg.atpleak._ref_atpm = ATP[mit].nodes[0]._ref_concentration

# place the mitochondria in the first segment
axon(1.0/axon.nseg).atpleak.position = 1

# link the flux in the mechanism to both the concentration of the axon
# and the mitochondria
for nd in ATP[cyt].nodes:
seg = nd.segment
nd.include_flux(seg.atpleak._ref_fluxA, units='mol/ms')
ATP[mit].nodes[0].include_flux(seg.atpleak._ref_fluxM, units='mol/ms')

# a rate for the mitochondria to produce ATP based on it's own state.
ATPsource = rxd.Rate(ATP, 10, regions=[mit])

h.finitialize(-70)

# to move the mitochondria change the 'position' parameter in the used in the
# mechanism.
def mMove():
for seg in axon:
if seg.atpleak.position == 1:
seg.atpleak.position = 0
axon(min(1 - 1.0/axon.nseg, seg.x + 1.0/axon.nseg)).atpleak.position = 1
break

# move the mitochondria every 10ms
for i in range(1,10):
h.CVode().event(i*10, mMove)

# run and plot the results every 25ms
for i in [25,50,75,100]:
h.continuerun(i)
pyplot.plot(numpy.linspace(0,axon.L,axon.nseg),
ATP[cyt].nodes.value, label="%1.0fms"%h.t)

pyplot.legend(frameon=False)
pyplot.xlabel("position ($\mu$m)")
pyplot.ylabel("concentration (mM)")
pyplot.show()


Code: Select all

NEURON {
SUFFIX atpleak
USEION atp READ atpi VALENCE -1
POINTER atpm
RANGE fluxA, fluxM, position
}

UNITS
{
(mol)   = (1)
(molar) = (1/liter)
(mM)    = (millimolar)
(um)    = (micron)
}
ASSIGNED {
atpi (mM)
atpm (mM)
area (um2)
diam (um)
vol  (um3)
fluxA (mol/ms)
fluxM (mol/ms)
}

INITIAL {
vol = area * diam / 4.0
}

PARAMETER {
leak = 1 (1/ms)
position = 0 (1)
}

BREAKPOINT {
fluxA = position * (1e-18) * leak * (atpm - atpi) *  vol
fluxM = -fluxA
}


If speed isn’t an issue for your simulation, you can also use a Python callback for include_flux instead of an nmodl mechanism. While this might be easier to define, it is also slower because execution needs to move from the compiled code back to Python every timestep.

2. You're correct, the method is checking if any of these terms include the rxdmath.v symbol representing the membrane voltage. Looking into this issue it seems the problem was using hasattr was doubling the number of checks at every level of the tree. We’ve fixed this in the current master branch and it is available in the neuron-nightly release (you can “pip install neuron-nightly” on Linux or MacOS). Thank you for bringing this to our attention and let us know if the problem persists.