hh.mod + Python + membrane potential at several points of the soma

When Python is the interpreter, what is a good
design for the interface to the basic NEURON
concepts.

Moderator: hines

Post Reply
smetens
Posts: 2
Joined: Wed Jun 09, 2021 4:48 am

hh.mod + Python + membrane potential at several points of the soma

Post by smetens »

Hello,
I am new with NEURON (8.0) and have some difficulties in the understanding of some structures and designs.
I use the Python (3.8) interface through the import of neuron libraries.
In particular, I begin with a soma with a given geometry and I add the Hodgkin–Huxley model for the dynamical part with the command
soma.insert('hh')

My questions are :

Question 1
: When we analyse what is behind hh.mod in the following link
https://raw.githubusercontent.com/neuro ... noc/hh.mod

I do not understand the structure of the lines :

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)
}


A-How do we introduce the rate equation for v in these lines of command ?
Do we add all currents on the right part of the equation ? And where is the left part of the equation C dv/dt ?

B- What is precisely the numerical method use to solve the 4 ODE’s of the HH model. Implicit, explicit, semi-implicit ? I expect that cnexp is the method but to what cnexp refers exactly ?
C- Can we change the method used ? e.g by an option somewhere ?
D- is there an hh.py version somewhere that we can use instead of hh.mod in the insert ("hh") command
E- If I want to change the model not using ModelDB but directly with Python is it possible ? is it possible to do it with PyNN http://neuralensemble.org/docs/PyNN/0.7/index.html

Question 2 :

In Neuron we can add some current clamp at a point of the soma by
stim = h.IClamp(soma(0.5))

A-How this current is added in the rate equation for v ? I do not see the way to introduce this current in hh.mod

B-In Neuron, soma has some extension, we can evaluate the membrane potential at a point of the membrane by soma_v.record(soma(0.5)._ref_v) (here in the middle of the soma). If I choice another point I have another profile for the membrane potential. Thus we have no space clamp, but hh.mod seems to solve ODE’s and not the PDE’s for the cable equation. Thus from where come the space here? Which equation is solved to have a space dependence of the membrane potential along the soma ?

Thanks a lot for those who have time to answer those beginners questions.
ramcdougal
Posts: 267
Joined: Fri Nov 28, 2008 3:38 pm
Location: Yale School of Public Health

Re: hh.mod + Python + membrane potential at several points of the soma

Post by ramcdougal »

I suggest in Python using the object-oriented form (h.hh) rather than the string ("hh") because it gives you an object whose properties you can explore (although your exploration options depend on whether or not you have the nmodl Python module installed...); i.e.

Code: Select all

soma.insert(h.hh)
You can then see the full code directly from Python via

Code: Select all

print(h.hh.code)
You'll notice from the PARAMETER block there are four parameters that you can change for this model: gnabar, gkbar, gl, and el. e.g. if you've ran the above, then you could always e.g.:

Code: Select all

for seg in soma:
    seg.hh.gnabar *= 2
From the ASSIGNED block, we see this model also depends on voltage (obviously), temperature (celsius), and the ek and ena reversal potentials. It also computes a bunch of things which you may or may not be interested in. (Speaking of temperature, note that the default is for h.celsius = 6.3, which is what the Hodgkin-Huxley channels expect as squid live in a cold environment, but this is not appropriate for mammals.)

From the NEURON block, we see that three of the values calculated, namely ina, ik, and il are to be used as currents. These show up in the right hand side of the Cm v' = ... equation. NEURON is by no means limited to this, but to a large extent, you can think of NEURON as a tool for solving the cable equation, and so it's your responsibility to give it the parameters, the morphology, and the currents, and it knows to put everything into the right place in that PDE. (For your second question, this is why hh.mod doesn't say anything about current clamp currents... it only provides the currents that it generates, which NEURON will combine with other currents as needed.)

Okay, so what about Cm? Well that's the parameter cm and can be set either for an entire section (e.g. soma.cm = 1) or for a specific segment (e.g. seg.cm = 1).

There's of course one other major parameter in cable theory, the axial resistivity or Ra. This is set on a per-Section basis via e.g. soma.Ra = 100. Note that the default is 35.4, which is appropriate for squid but too low for mammals.

The cnexp METHOD here is a statement about how the states should be advanced. This is explained in more detail elsewhere in the forum, but in brief it indicates that the state ODEs (for m, n, and h) are linear in v and thus can be advanced using analytical solutions.

By default, NEURON uses a first-order fixed step implicit solver with timestep 0.025 ms. (The latter can be changed via, e.g. h.dt = 0.0125). You can request a second-order crank-nicolson solution via h.secondorder = 2 (see the docs for more and why we use 2 here).

NEURON also supports implicit variable step, variable order solvers from the SUNDIALS suite, namely IDA/DASPK and CVode. Switching to variable step mode can be enabled via h.load_file("stdrun.hoc"); h.cvode_active(True) (there's a shorter version without loading the library, but this has some nice properties for the built-in graphics). The default is CVode unless you're simulating extracellular potentials or using linear mechanisms (e.g. simulating the experimental rig beyond a simple current injection), at which point it switches to IDA/DASPK. You can explicitly request IDA/DASPK from that point with h.cvode.use_daspk(True), but this is in general slower because it is not able to take advantage of the tree-structure of the matrix which was proven in 1984 to allow solving the matrix problem necessary for an implicit solver in O(n) time instead of O(n**3) time. Obviously, when you're talking about adaptive solvers, there are questions about tolerances and such (I generally use an absolute tolerance and a per-state variable scaling factor rather than relative tolerances), but you can learn all about that in the CVode docs.

MOD files are only one way of specifying kinetics that run at machine speed. Here's an example showing an equivalent specification of Hodgkin-Huxley in pure Python. Behind the scenes, starting to integrate will trigger the creation of a C file which is compiled and loaded in automatically. NEURON also supports a tool known as the ChannelBuilder which generates KSChan objects for channels meeting a technically-restricted but practically-reasonable set of kinetics.
smetens
Posts: 2
Joined: Wed Jun 09, 2021 4:48 am

Re: hh.mod + Python + membrane potential at several points of the soma

Post by smetens »

Dear ramcdougal,

thank you very much for this very precise answer to my questions.
I have to read very carefuly all the references you send us.
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: hh.mod + Python + membrane potential at several points of the soma

Post by ted »

Here https://neuron.yale.edu/neuron/publicat ... isms-nmodl you will find a link to a much expanded preprint of a paper that explains how NMODL is used to implement ion channels, transport mechanisms (pumps), and accumulation and diffusion of solutes. Very useful if you are trying to understand the contents of a .mod file.
If I want to change the model not using ModelDB but directly with Python is it possible ? is it possible to do it with PyNN
I don't understand this question. ModelDB is a searchable database of published computational neuroscience models. What does it mean to "use ModelDB to change a model"? With regard to what can and cannot be done with PyNN, you'll have to ask whoever supports PyNN. Use Python to change a model? Of course, if the model was implemented for NEURON. All user-specifiable aspects of a NEURON model's structure and parameters are changeable by hoc or Python statements.
How this current is added in the rate equation for v ? I do not see the way to introduce this current in hh.mod
You don't introduce it. And it isn't "introduced" in hh.mod. It is "introduced" in NEURON's computational engine when the statement
stim = h.IClamp(soma(0.5))
is executed.
we have no space clamp, but hh.mod seems to solve ODE’s and not the PDE’s for the cable equation. Thus from where come the space here? Which equation is solved to have a space dependence of the membrane potential along the soma ?
Every section is a cable of finite length. Every section has a discretization parameter nseg that specifies the number of "internal nodes" (spatial locations) at which the discretized cable equation will be solved during a simulation. Now go read chapter 5 of The NEURON Book, and you'll have a much clearer understanding of these and related issues. If you don't have the book, get https://www.neuron.yale.edu/ftp/ted/boo ... xedref.pdf and read it.
yuzhang
Posts: 2
Joined: Fri Oct 08, 2021 3:00 am

Re: hh.mod + Python + membrane potential at several points of the soma

Post by yuzhang »

Hello, excuse me, I find that the neuron engine sloves the cable hh equations, which have a second order derivation of v about the distance x, and thus are different from the point hh equations. In the cable hh equations, there are some parameters time the second order derivation, such as the diameter d, Rm, etc., how to choose these parameters reasonablely? In the attachment named "Ball_and_ Stick.py", If I choose the different values of "soma.L" and "soma.diam " , the results are not the same as that in the calssical HH model.

"Ball_and_ Stick.py":

from __future__ import division

from neuron import h, gui
import matplotlib.pyplot as plt

plt.ion()
###################
# # Build Model # #
###################

soma = h.Section(name="soma")
soma.L = 10 # um how to choose these parameters
soma.diam = 10 # um how to choose these parameters
soma.Ra = 100
soma.insert('pas')
soma.g_pas = 1 / 10000

dend = h.Section(name="dend")
dend.L = 500 # um
dend.diam = 1 # um
dend.Ra = 100 # ohm*cm
dend.insert('pas')
dend.g_pas = 1 / 10000

dend.connect(soma, 1, 0) # connect the end of the soma to the start of the dendrite

h("forall { nseg = int((L/(0.1*lambda_f(100))+0.9)/2)*2 + 1 }")
#########################
# # Set up experiment # #
#########################

stim = h.IClamp(soma(0.5)) # add a current clamp the the middle of the soma
stim.delay = 10 # ms
stim.dur = 100 # ms
stim.amp = 0.1 # nA

soma_v = h.Vector() # set up a recording vector
soma_v.record(soma(0.5)._ref_v) # record voltage at the middle of the soma

# Record voltage from all segments in the dendrite
dend_vs = []
for seg in dend:
dend_vs.append(h.Vector())
dend_vs[-1].record(seg._ref_v)

t = h.Vector()
t.record(h._ref_t) # record time.
h.v_init = -70 # set starting voltage
h.tstop = 200 # set simulation time
h.run() # run simulation
################
# Plot Results #
################
plt.figure(figsize=(8, 5))
plt.plot(t, soma_v, color='k', label='soma(0.5)')
for i, v in list(enumerate(dend_vs))[::3]:
plt.plot(t, v, color=(0, 0, 0.5 + 0.5 * i / len(dend_vs)),
label='dend({:.2})'.format(i / len(dend_vs)))

plt.xlim(0, 200)
plt.xlabel('Time (ms)', fontsize=15)
plt.ylabel('Voltage (mV)', fontsize=15)
plt.legend(fontsize=14, frameon=False)
plt.tight_layout()

###########################
# Add HH channels to soma #
###########################

h.tstop = 25
stim.dur = 5
soma.insert('kv') # add potassium channel
soma.gbar_kv = 2000 # set the potassium conductance

soma.insert('na') # add sodium channel
soma.gbar_na = 10000 # set the sodium conductance
h.celsius = 30

h.run()
plt.figure(figsize=(8, 5))
plt.plot(t, soma_v, color='k', label='soma(0.5)')
for i, v in list(enumerate(dend_vs))[::3]:
plt.plot(t, v, color=(0, 0, 0.5 + 0.5 * i / len(dend_vs)),
label='dend({:.2})'.format(i / len(dend_vs)))

plt.xlabel('Time (ms)', fontsize=15)
plt.ylabel('Voltage (mV)', fontsize=15)
plt.legend(fontsize=14, frameon=False)
plt.tight_layout()
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: hh.mod + Python + membrane potential at several points of the soma

Post by ted »

If I choose the different values of "soma.L" and "soma.diam " , the results are not the same as that in the calssical HH model.
Are you asking a question, or just making a statement? The model created by the python code in your post doesn't have the same anatomical structure as the classical HH model, nor does it have the same biophysical properties as the classical HH model (your code (1) uses a different value for cytoplasmic resistivity than the classical HH model does, (2) uses density mechanisms called na and kv whose source code is not available to me (so how do I know if they're the same as the classical HH sodium and potassium channels?), (3) even if they are the same, their channel densities are orders of magnitude larger than those used in the classical HH model, (4) uses pas instead of the HH mechanism's leak channels but specifies different channel density and reversal potential than HH's leak does, and (4) sets celsius to 30, which is way too hot for the classical HH mechanism (which expects celsius to be 6.3). And, unlike the classical HH model, it neither initializes to steady state, nor does it wait for the model to settle to steady state before the stimulus is applied.

Why should it produce the same results as the classical HH model does?
Post Reply