combining currents from mod files with rxd¶
A version of this notebook may be run online via Google Colab at https://tinyurl.com/neuron-rxd-and-mod (make a copy or open in playground mode).
Overview¶
NEURON's reaction-diffusion infrastructure can be used to readily allow intracellular concentrations to respond to currents generated in MOD files, as long as:
nrn_region='i'
is specified for therxd.Region
(so that it knows it corresponds to the electrophysiology region of the inside of the cell), AND the name and charge of the ion/etc are given in therxd.Species
declaration. Satisfying the above two rules also allows MOD files to see intracellular concentrations.3D extracellular concentrations also interoperate with electrophysiology automatically as long as name and charge are specified.
As a simple example, we consider a model with just a single point soma, of length and diameter 10 microns, with Hodgkin-Huxley kinetics (which uses the built in mod file hh.mod
), and dynamic sodium (declared using rxd but without any additional kinetics).
Setup NEURON library and imports¶
Let's import our usual NEURON libraries and definitions. Remember you can use either um
or µm
for micron.
from neuron import h, rxd
from neuron.units import mV, ms, um, mM
# needed for standard run system
h.load_file('stdrun.hoc')
1.0
Now import plotly
, a graphics library. (You could easily modify this code to use other graphics libraries like matplotlib
, plotnine
, or bokeh
.)
import plotly.graph_objects as go
from plotly.subplots import make_subplots
Setup the model¶
# define morphology
soma = h.Section(name='soma')
soma.L = soma.diam = 10 * um
# add ion channels (h.hh is built in, so always available)
h.hh.insert(soma)
# define cytosol. MUST specify nrn_region for concentrations to update
cyt = rxd.Region([soma], name='cyt', nrn_region='i')
# define sodium. MUST specify name and charge for concentrations to update
na = rxd.Species(cyt, name='na', charge=1)
Alternatively, we could have written h.hh.insert(soma.wholetree())
to put Hodgkin-Huxley channels everywhere in the cell that the soma is part of, but since there is only one section, this is not required.
Let's also add an excitatory synapse to receive events (these will trigger the cell to spike).
syn = h.ExpSyn(soma(0.5))
syn.tau = 1 * ms
syn.e = 0 * mV
Add a stimulus¶
The spike events themselves (two events, 15 ms apart starting at h.t=10*ms):
stim = h.NetStim()
stim.interval = 15 * ms
stim.number = 2
stim.start = 10 * ms
Send those events to our synapse:
nc = h.NetCon(stim, syn)
nc.weight[0] = 0.001
Setup recording variables¶
t = h.Vector().record(h._ref_t)
v = h.Vector().record(soma(0.5)._ref_v)
na_vec = h.Vector().record(soma(0.5)._ref_nai)
ina = h.Vector().record(soma(0.5)._ref_ina)
Initialize and run the simulation¶
h.finitialize(-65 * mV)
h.continuerun(50 * ms)
0.0
Plot it¶
fig = make_subplots(rows=3, cols=1)
fig.add_trace(go.Scatter(x=t, y=v, name="v"), row=1, col=1)
fig.update_yaxes(title_text="v (mV)", row=1, col=1)
fig.add_trace(go.Scatter(x=t, y=ina, name="ina"), row=2, col=1)
fig.update_yaxes(title_text="ina (mA/cm^2)", row=2, col=1)
fig.add_trace(go.Scatter(x=t, y=na_vec, name="[Na+]"), row=3, col=1)
fig.update_xaxes(title_text="t (ms)", row=3, col=1)
fig.update_yaxes(title_text="[Na+] (mM)", row=3, col=1)
fig.show()
Note: On homeostatic mechanisms or the lack thereof¶
Without any additional homeostatic mechanisms (Hodgkin and Huxley did not model sodium concentration so they did not need to include homeostatic mechanisms for it), intracellular sodium concentration will not return to baseline, and each spike will move intracellular sodium concentration closer to the extracellular concentration. Potassium concentration in this model is constant as we did not declare a potassium rxd.Species, but if we did it would also approach its extracellular concentration with each spike, until eventually the cell is not able to fire action potentials anymore.
Finally: The mod file¶
For any distributed mechanism, including the built-in h.hh
, one can always get the source code for the corresponding mod file using the .code
property; e.g.
print(h.hh.code)
TITLE hh.mod squid sodium, potassium, and leak channels COMMENT This is the original Hodgkin-Huxley treatment for the set of sodium, potassium, and leakage channels found in the squid giant axon membrane. ("A quantitative description of membrane current and its application conduction and excitation in nerve" J.Physiol. (Lond.) 117:500-544 (1952).) Membrane voltage is in absolute mV and has been reversed in polarity from the original HH convention and shifted to reflect a resting potential of -65 mV. Remember to set a squid-appropriate temperature (e.g. in HOC: "celsius=6.3" or in Python: "h.celsius=6.3"). See squid.hoc for an example of a simulation using this model. SW Jaslove 6 March, 1992 ENDCOMMENT UNITS { (mA) = (milliamp) (mV) = (millivolt) (S) = (siemens) } ? interface NEURON { SUFFIX hh REPRESENTS NCIT:C17145 : sodium channel REPRESENTS NCIT:C17008 : potassium channel USEION na READ ena WRITE ina REPRESENTS CHEBI:29101 USEION k READ ek WRITE ik REPRESENTS CHEBI:29103 NONSPECIFIC_CURRENT il RANGE gnabar, gkbar, gl, el, gna, gk : `GLOBAL minf` will be replaced with `RANGE minf` if CoreNEURON enabled GLOBAL minf, hinf, ninf, mtau, htau, ntau THREADSAFE : assigned GLOBALs will be per thread } PARAMETER { gnabar = .12 (S/cm2) <0,1e9> gkbar = .036 (S/cm2) <0,1e9> gl = .0003 (S/cm2) <0,1e9> el = -54.3 (mV) } STATE { m h n } ASSIGNED { v (mV) celsius (degC) ena (mV) ek (mV) gna (S/cm2) gk (S/cm2) ina (mA/cm2) ik (mA/cm2) il (mA/cm2) minf hinf ninf mtau (ms) htau (ms) ntau (ms) } ? currents BREAKPOINT { SOLVE states METHOD cnexp gna = gnabar*m*m*m*h ina = gna*(v - ena) gk = gkbar*n*n*n*n ik = gk*(v - ek) il = gl*(v - el) } INITIAL { rates(v) m = minf h = hinf n = ninf } ? states DERIVATIVE states { rates(v) m' = (minf-m)/mtau h' = (hinf-h)/htau n' = (ninf-n)/ntau } :LOCAL q10 ? rates PROCEDURE rates(v(mV)) { :Computes rate and other constants at current v. :Call once from HOC to initialize inf at resting v. LOCAL alpha, beta, sum, q10 : `TABLE minf` will be replaced with `:TABLE minf` if CoreNEURON enabled) TABLE minf, mtau, hinf, htau, ninf, ntau DEPEND celsius FROM -100 TO 100 WITH 200 UNITSOFF q10 = 3^((celsius - 6.3)/10) :"m" sodium activation system alpha = .1 * vtrap(-(v+40),10) beta = 4 * exp(-(v+65)/18) sum = alpha + beta mtau = 1/(q10*sum) minf = alpha/sum :"h" sodium inactivation system alpha = .07 * exp(-(v+65)/20) beta = 1 / (exp(-(v+35)/10) + 1) sum = alpha + beta htau = 1/(q10*sum) hinf = alpha/sum :"n" potassium activation system alpha = .01*vtrap(-(v+55),10) beta = .125*exp(-(v+65)/80) sum = alpha + beta ntau = 1/(q10*sum) ninf = alpha/sum } FUNCTION vtrap(x,y) { :Traps for 0 in denominator of rate eqns. if (fabs(x/y) < 1e-6) { vtrap = y*(1 - x/y/2) }else{ vtrap = x/(exp(x/y) - 1) } } UNITSON