The formulas for rxd.Rate, etc need to be continuous, but you can approximate a threshold function using a sigmoidal shaped curve. For example, you can use rxdmath.tanh to remove calcium at a constant rate when it exceeds a threshold with a narrow transition period:
Code: Select all
from neuron import h, rxd
from neuron.rxd import rxdmath
from matplotlib import pyplot
h.load_file('stdrun.hoc')
# NEURON model
dend = h.Section('dend')
# Where -- define the regions
cyt = rxd.Region(h.allsec(), name = 'cyt', nrn_region = 'i')
# Who -- define the species
ca = rxd.Species(cyt, name = 'ca', charge = 2, initial = 1.0)
# What -- define the reactions
# remove calcium at 0.01mM/ms when it is above 0.1mM using a sigmoid function
threshold = 0.1
steepness = 100
ca_removal = rxd.Rate(ca, -0.01 * (1+rxdmath.tanh(steepness*(ca-threshold)))/2)
# record & run the simulation
t_vec = h.Vector()
t_vec.record(h._ref_t)
ca_vec = h.Vector()
ca_vec.record(ca.nodes[0]._ref_concentration)
h.finitialize(-65)
h.continuerun(200)
# plot the result
pyplot.plot(t_vec, ca_vec)
pyplot.xlabel('t (ms)')
pyplot.ylabel('Ca$^{2+}$ (mM)')
pyplot.show()
You may also find the rxd tutorial useful;
https://www.neuron.yale.edu/neuron/stat ... index.html