Understanding a basic NEURON simulation completely

From Site
Jump to: navigation, search

Contents

My journey into NEURON

I wanted to write my own set of equations, and solve using Matlab or my own solver, because that way I know exactly what is going on. However, to make sure that I was on the right track, I wanted to run an equivalent simulation in NEURON, and compare the results. I needed a simple spatially-compact simulation of membrane dynamics with equations of Hodgkin-Huxley level of detail.

Installation

  • Installed the Windows version as on the website.
  • Got the Windows version to compile from source, following the instructions on the website.
    • One minor problem: ./configure failed, because it couldn't find libncurses. Ran Cygwin's setup.exe, and added that library (the development version). Then ./configure worked.

Goals

Goal I: Want to find the quintessential Hodgkin-Huxley equations somewhere inside NEURON. Are they in HOC files or MOD files? Which ones?
There is a demo that comes included with NEURON which seems to do this. Is this the standard Hodgkin-Huxley system, or a simplified version? Where are the equations to check?


Goal II: Want to run a standard simulation with the Hodgkin-Huxley equations I found above. Basically, to see the response to a current pulse.


Read through the tutorial called Construction and Use of Models: Part 1. Elementary tools, which showed how to get a minimalist HH model working with a current pulse. The script is:

   load_file("nrngui.hoc")
   
   create soma
   access soma
   
   soma {
       diam = 10
       L = 10/PI
       insert hh
   }

This seems like a great place to start with Neuron.

Next question: what happens when we insert hh? Where do the equations come from?

Using the index of the Reference documentation, we find documentation for the insert keyword. Also, there is an entry for insert hh, which points to $NEURONHOME/src/nrnoc/hh.mod.

Next, Kevin E. Martin's tutorial had more information about what's happening in HOC files:

  • A section is a cylinder.
  • The soma above is a section.
    • That's why it has a diameter and a length.
  • A mechanism (like hh) applies to a whole section.
  • A point process (like IClamp) applies to a point in space.

Goal I is now attained: the canonical HH model used by NEURON is in $NEURONHOME/src/nrnoc/hh.mod.

Goal II has also been attained: this is what that Elementary tools tutorial was all about above.


Getting the equations out of a MOD file

Generally, the equations are recognisable in the file, but it's good to know what all the other things are.

TITLE, COMMENT, ENDCOMMENT, and UNITS are straightforward to understand.

According to this description:

  • The NEURON block specifies ions and currents, which variables can be accessed by the user. One type of variable is RANGE.
  • The PARAMETER block has constants which may be set by the user.
    • That's why they have ranges in addition to default values and units.
  • The STATE block lists the state variables. They are usually solved for by code in the BREAKPOINT block.
  • The ASSIGNED block seems to be a list of variables which change, but are somehow different from state variables.
  • The BREAKPOINT block seems to contain some of the equations, and also, in the SOLVE statement, specifies the name of the block where more equations could be found.
  • The INITIAL block initialises variables: state variables and also some assigned variables it seems.

In hh.mod there are also DERIVATIVE, PROCEDURE, and FUNCTION blocks.

Seems to be rather complicated. But then the 124-line file gets turned into a 653-line hh.c, so there must be some advantages in not doing it directly in C.


MOD file reference

Reminder: I have found that the equations are in the hh.mod file. Now I want to get them out — they are not written there in a simple form.

I missed that this page was actually a table of contents for the reference manual. I thought it was an index of some sort — it looks like a bunch of keywords.

There are a number of headings, but the NMODL reference seems to be just two HTML documents: nmodl.html and nmodl2.html.

Unfortunately, this reference says that it only covers the differences between MODL and NMODL and that NMODL is a little extension of the substantial MODL. MODL is described in a manual which is "available from NBSR", but does not seem to be available anywhere on the web, which for us means that it doesn't exist. However, there is a small description in the SCoP Language Summary.

From nmodl.html we learn that:

  • In NMODL, PARAMETER and ASSIGNED variables are practically synonyms.
  • Membrane potential v is never a state variable (I was wondering about that before) because only NEURON is allowed to calculate it.
  • The DERIVATIVE block is for assigning values to the derivatives, which are part of differential equations governing the system.
  • The PROCEDURE block is a function that does not return a value.

From nmodl2.html we learn that:

  • The NEURON block is part of NMODL which has been added on top of those in MODL. It is more or less fully documented in this file.

Interpreting hh.mod

Now I'm going to go through each line in hh.mod and attack anything that I don't understand.

First, I don't understand what the following line means:

? interface

It's not a comment. A comment line is indicated by a leading colon. (Now I've seen everything.)
There doesn't seem to be an explanation of this anywhere. Since it seems that it doesn't do anything, ignore it.

In the NEURON block:

  • A SUFFIX is not a functional item, just a namespace control item.
  • The USEION block specifies which ions should be tracked. Each ion has a number of variables created to describe it: iion the current, ioni the internal concentration, iono the external concentration, and eion the reversal potential. The READ part tells us which variables will be read from, and the WRITE part tell us which variables will be updated. Each MOD file is a mechanism. One mechanism may be able to read ena and write ina, while another one could perhaps write ena. A mechanism could work with variables of a few different species (like a calcium-activated potassium channel). The combination of what mechanisms read and write which variables in a given section, should be kept sensible.
  • The NONSPECIFIC_CURRENT is not associated with any ion.
  • The RANGE part: these are variables which are declared in PARAMETER or ASSIGNED blocks which become range variables (which are?). I guess they are not global variables. Is there anything else we should know about "range" variables? Yes, see next section.
  • The GLOBAL part: these are PARAMETER or ASSIGNED variables which become global.
  • Some of the ASSIGNED variables were not mentioned in either RANGE or GLOBAL parts of the NEURON block, which means that they will become hidden from the user.

In the BREAKPOINT block we specify how to calculate the state of the model. First we SOLVE states and states is the name of the DERIVATIVE block that appears later. Then we assign currents ina, ik, and il.

In the DERIVATE states block, we first call PROCEDURE rates, which calculates right-hand-sides of the differential equations, based on the membrane potential. There is a TABLE statement in rates, which is a caching mechanism to speed up calculations.

Also, from the SCoP Language Summary, it seems that the mathematical statements are as one could expect.

Range variables

According to The NEURON Simulation Environment, range variables are actually continuous functions of normalised distance along a segment: from 0 at one end of the segment to 1 at the other. NEURON hides compartmentalisation from the user. There are actually a number of compartments underneath, so that a range variable will actually be an array, with one value per compartment.

What about global variables? Are they global to the model or to the segment (i.e. cylinder)? In the HH model, we have:

GLOBAL minf, hinf, ninf, mtau, htau, ntau

All of these are functions of local voltage. Although it's not explicitly clear whether or not the membrane voltage is a range variable or not, it must be. This suggests that this is a different global than conventional global variables in programming languages.

ODE integration methods used

The DERIVATIVE states block is solved by the METHOD cnexp. What exactly is that method? This line appears in hh.mod:

m' = (minf-m)/mtau

It gets translated into the following in hh.c:

m = m + (1 - exp(-dt  / mtau))
   *(minf - m) ;

After grepping $NEURONHOME, and finding cnexp mentioned in src/nmodl/solve.c and src/nmodl/deriv.c, there was no explanation there.

The document Expanding NEURON's Repertoire of Mechanisms with NMODL asserts that cnexp is a second-order efficient method which is suitable for stiff systems. It also asserts that cnexp is described elsewhere in the document, but it actually isn't.

However, in this tutorial for GENESIS, we find that the default numerical method in GENESIS is Exponential Euler, and the formula they give is, after some manipulation, exactly equivalent to the formula above.

Therefore cnexp is Exponential Euler.

GENESIS is a simulation tool similar to NEURON, but apparently developed independently. Here is a document showing that NEURON, GENESIS, and MOOSE (yet another simulator) all produce nearly identical results when simulating a series of cells.

The Exponential Euler method

Looks like it is based on the observation that (in Matlab):

dsolve('Dy = A - B*y')

ans =

   (A - C/exp(B*t))/B

Now, if we make this equal to y at t = 0, then we get C = A - B*y, and we get out the Exponential Euler method exactly. So it seems like a good method. This document gives hints about its applicability to stiff systems.

Conclusion

Finding these things out was enough to understand the first order idea of NEURON.

Further reading

To understand the details, I've read the following:

  • Kevin E. Martin's tutorial #1. Note that on Windows I had to use load_file("nrngui.hoc") rather than load_proc("nrnmainmenu"), which didn't work.