What is a NEURON equivalent of a point-neuron?

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

Moderator: hines

Post Reply
mi2n
Posts: 5
Joined: Wed Apr 12, 2023 6:36 am

What is a NEURON equivalent of a point-neuron?

Post by mi2n »

It is possible to define a dimensionless point neuron without using NEURON simulation environment as below:

Code: Select all


import numpy as np
from scipy.integrate import solve_ivp

# Define the Hodgkin-Huxley model parameters
Cm = 1.0   # membrane capacitance (uF/cm^2)
gNa = 120.0   # maximum sodium conductance (mS/cm^2)
ENa = 50.0    # sodium reversal potential (mV)
gK = 36.0     # maximum potassium conductance (mS/cm^2)
EK = -77.0    # potassium reversal potential (mV)
gL = 0.3      # maximum leak conductance (mS/cm^2)
EL = -54.387  # leak reversal potential (mV)

# Define the Hodgkin-Huxley model equations
def hh_model(t, y, I):
    V, m, h, n = y

    # Calculate the gating variable time constants
    alpha_m = 0.1 * (V + 40.0) / (1.0 - np.exp(-0.1 * (V + 40.0)))
    beta_m = 4.0 * np.exp(-0.0556 * (V + 65.0))
    alpha_h = 0.07 * np.exp(-0.05 * (V + 65.0))
    beta_h = 1.0 / (1.0 + np.exp(-0.1 * (V + 35.0)))
    alpha_n = 0.01 * (V + 55.0) / (1.0 - np.exp(-0.1 * (V + 55.0)))
    beta_n = 0.125 * np.exp(-0.0125 * (V + 65.0))

    # Calculate the membrane currents
    INa = gNa * m**3 * h * (V - ENa)
    IK = gK * n**4 * (V - EK)
    IL = gL * (V - EL)
    dVdt = (I - INa - IK - IL) / Cm

    # Calculate the gating variable changes
    dmdt = alpha_m * (1.0 - m) - beta_m * m
    dhdt = alpha_h * (1.0 - h) - beta_h * h
    dndt = alpha_n * (1.0 - n) - beta_n * n

    return [dVdt, dmdt, dhdt, dndt]

# Set the initial conditions and input current
y0 = [-65.0, 0.05, 0.6, 0.32]
I = 10.0   # input current (uA/cm^2)

# Set the time range for the simulation
t_start = 0.0
t_end = 100.0
dt = 0.01
t_span = np.arange(t_start, t_end+dt, dt)

# Solve the differential equations using solve_ivp
sol = solve_ivp(fun=lambda t, y: hh_model(t, y, I),
                t_span=[t_start, t_end],
		y0=y0,
                method='RK45',
                t_eval=t_span)

# Plot the results
import matplotlib.pyplot as plt
plt.plot(sol.t, sol.y[0], color='k')
plt.show()

The same neuron, when implemented using NEURON, look like this:

Code: Select all


from neuron import h, gui

# Create a new point neuron
soma = h.Section()
soma.L = 20.0    # length (um)
soma.diam = 20.0 # diameter (um)

# Define the Hodgkin-Huxley model parameters
soma.insert("hh") # insert the HH mechanism
soma(0.5).hh.gnabar = 0.12  # sodium conductance (S/cm^2)
soma(0.5).hh.gkbar = 0.036  # potassium conductance (S/cm^2)
soma(0.5).hh.gl = 0.0003  # leak conductance (S/cm^2)
soma(0.5).hh.el = -54.3  # reversal potential for the leak current (mV)

# Create a current clamp to stimulate the neuron
stim = h.IClamp(soma(0.5))
stim.delay = 100.0  # delay before the stimulus (ms)
stim.dur = 500.0    # duration of the stimulus (ms)
stim.amp = 0.1885      # amplitude of the stimulus (nA)

# Create a vector to record the membrane potential
time = h.Vector()
voltage = h.Vector()
time.record(h._ref_t)
voltage.record(soma(0.5)._ref_v)

# Run the simulation
h.finitialize(-65.0)
h.continuerun(700.0)

# Plot the results
import matplotlib.pyplot as plt
plt.plot(time, voltage)
plt.xlabel('Time (ms)')
plt.ylabel('Voltage (mV)')
plt.show()

Here we have to provide a dimension for the soma. I gave 20 um length and 20 um diameter - just a set of dummy values. But the outputs of these models don't match. The NEURON version doesn't produce any action potential. I calculated the neuron current (0.1885 nA) as the product of point neuron current (10 uA/cm2) and the surface area of soma (18.85e-6 cm2). Was there something more to this? How to get the stimulus amplitude that is equivalent to that of the point-neuron version?

Any help would be greatly appreciated.
ted
Site Admin
Posts: 6286
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: What is a NEURON equivalent of a point-neuron?

Post by ted »

Nice, concise python.
The same neuron, when implemented using NEURON
"The same neuron"? Well, your own results demonstrate that the two implementations are not the same. And many of the differences are easy to discover.

If you haven't already done this, get NEURON's source code from github and find the file hh.mod (I see it in nrn/src/nrnoc/). Compare the rate constants' parameters, and the reversal potentials of the various ionic currents, with the corresponding values in your pure python implementation?

Another easily discoverable difference is in the initialization of the gating states m, h, and n. In your implementation of the NEURON model, right after
finitialize()
print out the values of soma(0.5).hh.m, soma(0.5).hh.h, and soma(0.5).hh.n
The difference is small, but it would matter if the stimulus current started before about 20 ms into the simulation. Hmm, when does the stimulus start in your pure python implementation?

What about the stimulus amplitude? FYI an HH model with a surface area of 100 um2 generates robust repetitive spiking in response to a sustained stimulus current of 0.01 nA. A 0.01 nA current applied to a model with 100 um2 surface area corresponds to a transmembrane current density of 0.01 mA/cm2. Given a specific membrane capacitance of 1 uf/cm2, the ratio of stimulus current to total membrane capacitance is 0.01 mA/uf.

The surface area of the model you implemented for NEURON is about 1257 um2, so it needs 0.1257 nA to achieve the same stimulus current/membrane capacitance ratio. What is the actual ratio that you implemented?

Given the parameters in your pure python implementation, what is the ratio of its stimulus current to its total membrane capacitance? Is this the same as in your NEURON implementation?
mi2n
Posts: 5
Joined: Wed Apr 12, 2023 6:36 am

Re: What is a NEURON equivalent of a point-neuron?

Post by mi2n »

Thanks a lot for your clarification. I wasn't very thorough with the experiment. There was a mistake in evaluating the area as well. I included the flat surfaces of the soma for the computation. It seems only the curved surface area of the cylinder is involved in the ion conductance - is that so?

Earlier, I didn't wait 20 ms to apply the current stimulation. After adding the delay of 100 ms, the outputs seem to match.

Image

The small difference might be caused by the small difference in the m,n,h, and Vm values at the time of the application of the stimulus. The parameter values at 100 ms (at the onset of stimulation) are given below:

Code: Select all

-------------------------------------------
for Pure Python implementation
-------------------------------------------
Vm	=-64.97337788385995 mV

m	=0.05305467054932449
n	=0.31807373890240004
h	=0.5952146975685233

I/Cm	=0.01 mA/uF

-------------------------------------------------
for Python-NEURON implementation
-------------------------------------------------
Vm	=-64.97405244991657 mV

m	=0.05309465179576391
n	=0.31807461731315
h	=0.595213019322433

I/Cm	=0.009999999999999997 mA/uF
Once again, thank you so much for helping me with this.
ted
Site Admin
Posts: 6286
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: What is a NEURON equivalent of a point-neuron?

Post by ted »

That's impressive. Thanks for following up on this!
There was a mistake in evaluating the area as well. I included the flat surfaces of the soma for the computation. It seems only the curved surface area of the cylinder is involved in the ion conductance - is that so?
Yes. That should really be stated quite explicitly in the documentation. Actually it is, but somehow very few people stumble across it.

Go to
nrn.readthedocs.io
and search for
section

Scroll down the results until you find
Conceptual Overview of Sections

Be sure to read that whole page. It might contain other helpful information.

I wonder how many current NEURON users have ever read that entire page . . . and then re-read it to catch the stuff they missed the first time.
Earlier, I didn't wait 20 ms to apply the current stimulation. After adding the delay of 100 ms, the outputs seem to match.
Yep, initialization is important.
The small difference might be caused by the small difference in the m,n,h, and Vm values at the time of the application of the stimulus.
Those differences are pretty small. Could easily be due to a combination of < differences in execution sequence > and < accumulation of roundoff error > over the course of 4000 time steps.
Post Reply