Changing parameters in NMODL files from within a Python script?

[espenh]

Hi,

How can you change (e.g., by looping over a range of values) a parameter of an NMODL file (compiled and made available to NEURON) from within a Python script used to run a compartmental model that uses the NMODL mechanism?
One example could be that you want to explore the effect of changing a specific rate constant for a channel implemented as a Markov scheme without having to generate separate NMODL files for each rate constant value.

Thanks a lot for your help,
Espen

[urid]

Here's an example of some code which changes the amplitude of a current to see the change in the mean firing frequency.
It sets up only one object for the soma and current clamp, and loops through different values of the amplitude of injected current.

Code: Select all

from neuron import h
import numpy as np
import matplotlib.pyplot as plt
h.load_file("stdrun.hoc")

def mean_spike_freq(t,v,spike_thresh=-10): # Function to detect spikes and return spike frequency
    crossings = np.where((v[:-1] < spike_thresh) & (v[1:] >= spike_thresh))[0]
    spike_timing = np.diff(t[crossings]) # time between spikes
    if len(spike_timing) == 0: # 0 Hz if no spikes
        return 0
    return 1000 / np.mean(spike_timing) # Mean frequency in Hz

soma = h.Section(name="soma") # Create section
soma.insert("hh") # Insert the Hodgkin-Huxley mechanism

ic = h.IClamp(soma(0.5)) # Attach a current clamp
ic.delay = 200
ic.dur = 600

currents = np.linspace(6,25,num=20) # Rheobase for HH is around 10
frequencies = np.zeros_like(currents)

trec = h.Vector().record(h._ref_t)
vrec = h.Vector().record(soma(0.5)._ref_v)

for i,amp in enumerate(currents):
    ic.amp = amp # set the current amplitude

    h.finitialize(-65) # calls h.frecord_init() and resets dynamical variables
    h.continuerun(1000) # run for 1000 ms

    f = mean_spike_freq(trec.as_numpy(),vrec.as_numpy())
    frequencies[i] = f

plt.plot(currents,frequencies)
plt.xlabel("Injected Current (nA)")
plt.ylabel("Spike frequency (Hz)")

[ted]

Good Python suggestion, urid.

Here are some NEURON-specific comments:

First, just to make sure that everyone who reads this thread knows what we're writing about--

1. NMODL's global variables.

An NMODL-specified mechanism's PARAMETERs are global by default (have the same value in all segments in which the mechanism exists), but this can be overridden by RANGE declarations in the NEURON block. Example: consider a kinetic scheme specification of a potassium channel described by

Code: Select all

NEURON {
        SUFFIX khh
        USEION k READ ek WRITE ik
        RANGE g, gbar
}
PARAMETER {
        gbar = 22 (millimho/cm2)
        d1 = 21 (millivolt)
        k1 = .2 (/millivolt)
        d2 = 43 (millivolt)
        k2 = .036 (/millivolt)

        ta1 = 4.4 (ms)
        tk1 = -.025 (/millivolt)
        ta2 = 2.6 (ms)
        tk2 = -.007 (/millivolt)
 }
ASSIGNED {
        g (millimho/cm2)
 . . .
        a1(/ms) b1(/ms) a2(/ms) b2(/ms)
}
STATE {c1 c2 o}
KINETIC kin {
        rates(v) :
        ~ c1 <-> c2     (a1, b1)
        ~ c2 <-> o      (a2, b2)
        CONSERVE c1 + c2 + o = 1
}

(other details omitted for clarity).

The variable gbar is declared in the PARAMETER block, which is appropriate because gbar is a parameter. By default gbar would be global (would have the same value in every segment that has khh), but the RANGE statement in the NEURON block turns it into a range variable, so the value of gbar may be different in different segments.

A side comment: that same RANGE statement ensures that g is visible to hoc and Python, even though g is declared in the ASSIGNED block (an appropriate declaration because g appears on the left hand side of an assignment statement in one of khh's equation blocks). By default, ASSIGNED variables are not visible outside of the mechanism that declared them, so g would have been inaccessible via hoc or Python.

Now back to the main discussion:

The other PARAMETERs are used in formulas (not shown) from which the forward and backward rate constants (ASSIGNED variables a1, b1, a2, and b2) are calculated, and they do not appear in RANGE statements, so they retain their default status as globals.

2. hoc and Python syntax for accessing NMODL globals.

NMODL globals are not specific to any particular segment, so the syntax for accessing them is different from what you'd use for a range variable. For this particular example, the hoc name for parameter d1 would be d1_khh, and the Python name would be h.d1_khh (no need to specify section or segment).