This document describes how to use the NBSR model description language to add membrane mechanisms to NEURON. The kinds of mechanisms that can be added are:
Mechanisms are normally local. That is they do not depend on what is happening at other places on the neuron. However, a method exists for writing mechanisms that depend on variables of mechanisms at other locations. For example the calcium concentration at a presynaptic mechanism can be used to calculate the conductance change at a postsynaptic mechanism. (See, Importing variables from other mechanisms.) Also, FUNCTION's written in a model are global and may be used in other models if they do not involve range variables.
MODL (model description language) was originally developed at the NBSR (National Biomedical Simulation Resource) to specify models for simulation with SCoP (Simulation Control Program). With MODL one specifies a physical model in terms of simultaneous nonlinear algebraic equations, differential equations, or kinetic schemes. MODL translates the specification into the C language which is then compiled and linked with the SCoP program. It turned out that only modest extensions to the MODL syntax were necessary to allow it to translate model descriptions into a form suitable for compiling and linking with NEURON V2. The extended version was called NMODL. In NEURON V3 the advent of the object oriented interpreter, OC, allowed Point Processes to be treated as objects instead of parallel arrays of variables. The model description translator that emits code suitable for NEURON V3 is called NOCMODL. NMODL and NOCMODL handle identical input model descriptions, they differ merely in the output interface code. A prototype model description translator has been written to generate code suitable for linking with GENESIS.
This document discusses only the differences between NMODL and MODL. A complete user manual for the standard model description language is available from NBSR. A brief description of MODL is in the document entitled, ``SCoP Language Summary''.
The easiest way to write membrane mechanisms is by analogy with the examples. The example files come in pairs with a .mod and .hoc extension. Models (membrane mechanisms) are linked into neuron with the command:
nrnivmodl file1 file2 file3 ...In the list of .mod files you do not type the extension, only the prefix. If there are no files specified then nrnivmodl will use all the mod files in the current working directory. When nrnivmodl is finished, there will exist a version of NEURON called ``special'' which contains those membrane mechanisms described in the files. ``Special'' should be renamed to something more suitable. The associated .hoc files can be executed by ``special'' to test various aspects of the models.
It is extremely important that mechanisms have consistent units. To ensure this use the command:
modlunit fileleaving off the file extension. For more information about units.
Our first nerve simulation program, CABLE, contained several built-in membrane mechanisms, including radial calcium diffusion, calcium channel, calcium activated potassium channel, calcium pump, etc. However, in practice, only the Hodgkin-Huxley squid channels were enough of a standard to be used ``as is'' across more than one series of simulations. The other channels all required some type of modification to be useful as new situations arose. Sometimes the modifications were minor, such as changing the coordinate system for radial calcium diffusion so that there were more compartments near the membrane, but often we were forced to add an entirely new mechanism from scratch such as Frankenhaeuser-Huxley channels for Xenopus node. The problem was greatly compounded for other users of CABLE who needed to add new channels but were not familiar with the numerical issues or the detailed interface requirements. NMODL with NEURON is a significant improvement over CABLE with regard to adding new membrane mechanisms:
Membrane mechanisms deal with currents, concentrations, potentials, and state variables and it is helpful to know how NEURON treats these variables in order to correctly write a new membrane mechanism.
NEURON integrates its equations using the function fadvance(). During a call to this function the value of the global time variable, t, is increased by the value of dt (t = t + dt), and all the voltages, currents, concentrations, etc. are changed to new values appropriate to the new value of time. The default numerical method used by NEURON produces values which have an error proportional to dt. That is, it makes no sense to ask at what time in the interval are the values most accurate. However, by setting the global variable secondorder equal to 2, the values produced by fadvance have errors proportional to dt^2 and it is important to realize that
Fadvance() goes about its business by first setting up the current conservation matrix equation to be used in the calculation of membrane potential. To do this it calls the current functions for each mechanism in each segment which compute conductance using the old values of states and current using the old values of states and membrane potential. The value of time when the BREAKPOINT block is called is t+dt/2 so models which depend explicitly on time will be second order correct it they use the value of t}. Fadvance then solves the matrix equation for the new value of the membrane potential. Depending on the value of secondorder it then may re-call these current functions with the average of the new and old membrane potentials to get an accurate final value of the current. It then calls the state integrator functions using the new value of the membrane potential and the second order correct currents to calculate the new values of the states. The details of this method can be gleaned from the file nrn/src/nrnoc/fadvance.c.
It is therefore necessary for NMODL to divide up the statements properly into a current function and a state function. It also has to create the interface between model variables and NEURON and create a memory allocation function so that segments have separate copies of each variable. Finally, it has to make sure that local model currents get added to the correct global ionic currents.
Note: This simulation method is very effective and highly efficient when currents depend on membrane potential and ionic concentrations do not change on the same time scale as the membrane potential. When these conditions are not met, however, such as in a calcium pump mechanism in which the current depends on the concentrations of calcium next to the membrane, one must be careful to use a dt small enough to prevent the occurrence of numerical instabilities. (Or else using a single model to describe both the pump current and that current's effect on concentration so that the concentrations and pump states may be computed simultaneously. An example of such a model is in nrn/demo/release/cabpump.mod) A future version of NEURON will have the option (slightly less efficient) of calculating all state variables simultaneously so that numerical stability is guaranteed. Further discussion of the numerical method's used by NEURON are found here.
Only a small part of the full model description language is relevant to neuron mechanisms. The important concepts held in common are the declaration of all variables as
INCLUDE
The INCLUDE statement replaces itself with the contents of the indicated file.
eg.
INCLUDE "units.inc"If the full path to the file is not given, the file is first looked for in the current working directory, then in the directory where the original .mod file was located, and then in the directories specified by the colon separated list in the environment variable MODL_INCLUDES. Notice that the INCLUDE filename explicitly requires a complete file name --- don't leave off the suffix, if any.
Other blocks which play similar roles in NMODL and MODL are
SOLVE deriv METHOD euler or SOLVE deriv METHOD runge or SOLVE deriv METHOD derivimplicitbecause the default integration method is a variable time step runge-kutta method which cannot work in a NEURON context. The first two methods above are computationally cheap but are numerically unstable when equations are stiff (states vary a lot within a dt step). HH type mechanisms have state equations which are particularly simple and extra efficiency and accuracy is easily obtained by integrating the states analytically. The hh2.mod example shows how to do this.
~ expr = exprWhen this block is called by the SOLVE statement, the values of the states are computed so that the equations are true. The default method used is newton's method. These kinds of equations can also appear within a DERIVATIVE block.
However if a procedure is called by the user, and it makes use of any range variables, then the user is responsible for telling the mechanism from what location it should get its range variable data. This is done with the hoc function:
setdata_mechname(x)where mechname is the mechanism name. For range variables there must of course be a currently accessed section. In the case of Point processes, one calls procedures using the object notation
pp_objref.procname()In this case procname uses the instance data of the point process referenced by pp_objref.
Sometimes, state equations are so simple, e.g. HH states, that significant efficiency gains and extra accuracy are obtainable by a special integration procedure. In this case the procedure can be called by a SOLVE statement and actually integrates the states (but don't call it directly at the user level!). If a PROCEDURE is solved by a SOLVE statement it may return an error code (By default it returns an error code of 0 which denotes success.) To return a non-zero error code use the idiom
VERBATIM return ...; ENDVERBATIM
TABLE
The TABLE statement is very useful in a NEURON context because of
the potentially great increase in speed of simulation. Often rate functions are
complicated functions of the voltage and is is very expensive to calculate
their values over and over at every segment. By using tables of rate
coefficients, it is not uncommon to improve simulation speed by a factor
of 5.
In the context of a procedure taking one argument, TABLE has the syntax
TABLE variables DEPEND dependencies FROM lowest TO highest WITH tablesize
where: variables is a list of variable names each of which will have its own table. dependencies is a list of parameters that, when any of them changes their value, cause the tables to be recalculated. lowest is the least arg value for the first table entry. highest is the greatest arg value for the last table entry. table size is the number of elements in each table.
Each model that has a table also has a flag associated with it that can be changed by the user called usetable_suffix which specifies that the tables are to be used (1, default) or not used (0).
With usetable_suffix = 0, when the procedure is called it ignores the tables and just computes the values using the assignment statements as any normal procedure.
With usetable_suffix = 1, when the procedure is called, the arg value is used to assign values to the "variables" by looking them up in the tables and the time normally spent executing the assignment statements is saved. If the tables are out of date (a "dependency" has a different value from its value the last time the tables were constructed) or had never been created, the tables are created.
Note that updating tables with tablesize=200 is equivalent to calling the procedure 200 times with different values of the argument. This investment is only repaid if the tables remain valid for many more than 200 subsequent calls to the procedure and if the calculation takes more time than an interpolated table lookup.
INITIAL
The INITIAL block is called when the user executes the finitialize() function
from hoc. Just prior to executing the user code in the INITIAL block (and
if an INITIAL block does not exist) all
states are set to the values in the state0 variables (default 0). It may
be useful to declare some state0 variables as GLOBAL or RANGE in the NEURON
block in order to give some user control over the default initialization
of states. In the INITIAL block these or any other variables may be set
to calculated values. Note that
states can also be initialized explicitly by the user at the hoc level.
The case where an ionic variable is also a STATE requires some care to deal with properly in the INITIALIZE block. The problem is that the ionic variable, eg. cai, is actually the value of a local copy of the ionic variable which is located in the variable named _ion_cai. Because of the order of copying and default initialization, cai is always initialized to 0 regardless of the global value of cai and on exit the global value of cai is then set to 0 as well. The way to avoid this is either to make sure the state0 variable, cai0, is set properly or (I believe more preferably), set the local cai variable explicitly using the global cai variable with a VERBATIM statement within the INITIAL block. The idiom is:
VERBATIM cai = _ion_cai; ENDVERBATIM
Many other features of the model description language, such as DISCRETE blocks, sensitivity analysis, optimization are not relevant in the NEURON context and may or may not produce meaningful translations. Since NMODL produces a c file, it is possible for the highly motivated to modify that file in order to do something implementation dependent. In this regard, the VERBATIM block can be used to place c code within the model description file.
This section describes the special NEURON block that has been added to the standard model description language in order to allow translation of a model into a form suitable for linking with NEURON.
The keyword NEURON introduces a special block which contains statements that tell NMODL how to organize the variables for access at the NEURON user level. It declares:
NEURON {
SUFFIX ...
RANGE ...
GLOBAL ...
NONSPECIFIC_CURRENT ...
USEION ... READ ... WRITE ... VALENCE real
POINT_PROCESS ...
POINTER ...
EXTERNAL ...
}
The suffix, ``_name'' is appended to all variables, functions, and procedures that are accessible from the user level of NEURON. If the SUFFIX statement is absent, the file name is used as the suffix (with the addition of an underscore character. If there is a POINT_PROCESS statement, that name is used as the suffix. Suffixes prevent overloading of names at the user level of NEURON. At some point in the future I may add something similar to the access statement which will allow the omission of the suffix for a specified mechanism. Note that suffixes are not used within the model description file itself. If the SUFFIX name is the word, "nothing", then no suffix is used for variables, functions, and procedures explicitly declared in the .mod file. However, the mechanism name will be the base file name. This is useful if you know that no conflict of names will exist or if the .mod file is primarily used to create functions callable from NEURON by the user and you want to specify those function names exactly.
These names will be become range variables. Do not add suffixes here. The names should also be declared in the normal PARAMETER or ASSIGNED statement outside of the NEURON block. Parameters that do not appear in a NEURON RANGE statement will become global variables. Assigned variables that do not appear in this statement or in the NEURON GLOBAL statement will be hidden from the user. When a mechanism is inserted in a section, the values of these range variables are set to the values specified in the normal PARAMETER statement outside the NEURON block.
These names, which should be declared elsewhere as ASSIGNED or PARAMETER variables, become global variables instead of range variables. Notice here that the default for a PARAMETER variable is to become a global variable whereas the default for an ASSIGNED variable is to become hidden at the user level.
This signifies that we are calculating local currents which get added to the total membrane current but won't contribute to any particular ionic concentration. This current should be assigned a value after any SOLVE statement but before the end of the BREAKPOINT block. This name will be hidden at the user level unless it appears in a NEURON RANGE statement.
This statement declares that a specific ionic species will be used within this model. The built-in HH channel uses the ions na and k. Different models which deal with the same ionic species should use the same names so that total concentrations and currents can be computed consistently. The ion, Na, is different from na. The example models using calcium call it, ca. If an ion is declared, suppose it is called, ion, then a separate mechanism is internally created within NEURON, denoted by ion, and automatically inserted whenever the ``using'' mechanism is inserted. The variables of the mechanism called ion are outward total current carried by this ion, iion; internal and external concentrations of this ion, ioni and iono; and reversal potential of this ion, eion. These ion range variables do NOT have suffixes. Prior to 9/94 the reversal potential was not automatically calculated from the Nernst equation but, if it was used it had to be set by the user or by an assignment in some mechanism (normally the Nernst equation). The usage of ionic concentrations and reversal potential has been changed to more naturally reflect their physiological meaning while remaining reasonably efficient computationally.
The new method governs the behaviour of the reversal potential and concentrations with respect to their treatment by the GUI (whether they appear in PARAMETER, ASSIGNED, or STATE panels; indeed, whether they appear at all in these panels) and when the reversal potential is automatically computed from the concentrations using the Nernst equation. The decision about what style to use happens on a per section basis and is determined by the set of mechanisms inserted within the section. The rules are defined in the reference to the function ion_style(). Three cases are noteworthy.
Assume only one model is inserted in a section.
USEION ca READ ecaThen eca will be treated as a PARAMETER and cai/cao will not appear in the parameter panels created by the gui.
Now insert another model at the same section that has
USEION ca READ cai, caoThen 1) eca will be "promoted" to an ASSIGNED variable, 2) cai/cao will be treated as constant PARAMETER's, and 3) eca will be computed from the Nernst equation when finitialize() is called.
Lastly, insert a final model at the same location in addition to the first two.
USEION ca WRITE cai, caoThen eca will still be treated as an ASSIGNED variable but will be computed not only by finitialize but on every call to fadvance(). Also cai/cao will be initialized to the global variables cai0_ca_ion and cao0_ca_ion respectively and treated as STATE's by the graphical interface.
The idea is for the system to automatically choose a style which is sensible in terms of dependence of reversal potential on concentration and remains efficient.
Since the nernst equation is now automatically used as needed it is necessary to supply the valence (charge carried by the ion) except for the priviledged ions: na, k, ca which have the VALENCE 1, 1, 2 respectively.
Only the ion names na, k, and ca are initialized to a physiologically meaningful value --- and those may not be right for your purposes. Concentrations and reversal potentials should be considered parameters unless explicitly calculated by some mechanism.
The READ list of a USEION specifies those ionic variables which will be used to calculate other values but is not calculated itself. The WRITE list of a USEION specifies those ionic variables which will be calculated within this mechanism. Normally, a channel will read the concentration or reversal potential variables and write a current. A mechanism that calculates concentrations will normally read a current and write the intracellular and/or extracellular; it is no longer necessary to ever write the reversal potential as that will be automatically computed via the nernst equation. It usually doesn't make sense to both read and write the same ionic concentrations. It is possible to READ and WRITE currents. One can imagine, a large calcium model which would WRITE all the ion variables (including current) and READ the ion current. And one can imagine models which READ some ion variables and don't WRITE any. It would be an error if more than one mechanism at the same location tried to WRITE the same concentration.
A bit of implementation specific discussion may be in order here. All the statements after the SOLVE statement in the BREAKPOINT block are collected to form a function which is called during the construction of the charge conservation matrix equation. This function is called several times in order to compute the current and conductance to be added into the matrix equation. This function is never called if you are not writing any current. The SOLVE statement is executed after the new voltages have been computed in order to integrate the states over the time step, dt. Local static variables get appropriate copies of the proper ion variables for use in the mechanism. Ion variables get updated on exit from these functions such that WRITE currents are added to ion currents.
Point models are used for synapses, electrode stimuli, etc. They are distinguished from standard mechanisms in that instead of ``insert''ing the mechanism into a section and accessing parameters via range variables, point mechanisms are created as interpreter objects, eg.
objref stim stim = new IClamp(x)This differs from the previous release of NEURON V2 where point processes were created as a vector of mechanisms each of which is given a particular location within a specified section. Now, instead of the vector of mechanisms being represented as a set of parallel vectors containing the values of STATE's, PARAMETER's, etc., values are accessed via the standard object syntax, eg. stim.amp = 2. Since standard mechanisms are considered in terms of density, the appropriate current units for standard mechanisms are mA/cm2 and conductance units are mho/cm2. However, point process current units are nA and conductance units are umho. These conventions ensure that the simulation is independent of the number of segments in a section (assuming the number of segments is large enough so spatial discretization error is small.
At the NEURON user level, all variables and functions associated with a POINT_PROCESS are accessed via the normal object syntax. A point process, call it pp is inserted into (or moved to) the currently specified section at location, 0 < x < 1, with the function, pp.loc(x). If a point process is created with no argument then it is not located anywhere. If an argument is present and there is a currently accessed section then the point process is placed there. At this time point process's are placed at the center of the nearest segment and may be placed explicitly at 0 or 1 (0 area segments at the ends of sections).
pp.has_loc() returns 1 if the point process is located in a section and returns 0 if not located. If a point process has no location then attempts to access its variables or get its location will produce an error message.
One finds the location of a point process via the function, x = pp.get_loc(). The function returns the x location at the center of the segment where the process was placed and pushes the section name onto the stack so that it becomes the currently accessed section. The stack must be popped with pop_section() at a subsequent time. BE SURE TO POP THE SECTION STACK! This can be a dangerous function in the sense that if the stack is not popped, then section access is completely screwed up.
These names are pointer references to variables outside the model. They should be declared in the body of the description as normal variables with units and are used exactly like normal variables. The user is responsible for setting these pointer variables to actual variables at the hoc interpreter level. Actual variables are normal variables in other mechanisms, membrane potential, or any hoc variable. See below for how this connection is made. If a POINTER variable is ever used without being set to the address of an actual variable, NEURON may crash with a memory reference error, or worse, produce wrong results. Unfortunately the errors that arise can be quite subtle. For example, if you set a POINTER correctly to a mechanism variable in section a. And then change the number of segments in section a, the POINTER will be invalid because the memory used by section a is freed and might be used for a totally different purpose. It is up to the user to reconnect the POINTER to a valid actual variable.
These names, which should be declared elsewhere as ASSIGNED or PARAMETER variables allow global variables in other models or NEURON c files to be used in this model. That is, the definition of this variable must appear in some other file. Note that if the definition appeared in another mod file this name should explicitly contain the proper suffix of that model. You may also call functions from other models (but don't ignore the warning; make sure you declare them as
extern double fname_othermodelsuffix();in a VERBATIM block and use them with the proper suffix.
Occasionally mechanisms need information from other mechanisms which may be located elsewhere in the neuron. Connecting pre and post synaptic point mechanisms is an obvious example. In the same vein, it may be useful to call a function from hoc which modifies some mechanism variables at a specific location. (Normally, mechanism functions callable from HOC should not modify range variables since the function doesn't know where the mechanism data for a segment is located. Normally, the pointers are set when NEURON calls the BREAKPOINT block and the associated SOLVE blocks.)
One kind of connection between mechanisms at the same point is through ionic mechanisms invoked with the USEION statement. In fact this is entirely adequate for local communication although treating an arbitrary variable as an ionic concentration may be conceptually strained. However it doesn't solve the problem of communication between mechanisms at different points.
Basically what is needed is a way to implement the hoc statement
section1.var1_mech1(x1) = section2.var2_mech2(x2)efficiently from within a mechanism without having to explicitly connect them through assignment at the HOC level everytime the var2 might change.
First of all, the variables which point to the values in some other mechanism are declared within the NEURON block via
NEURON {
POINTER var1, var2, ...
}
These variables are used exactly like normal variables in the sense that
they can be used on the left or right hand side of assignment statements
and used as arguments in function calls. They can also be accessed from HOC
just like normal variables.
It is essential that the user set up the pointers to point to the correct
variables. This is done by first making sure that the proper mechanisms
are inserted into the sections and the proper point processes are actually
``located'' in a section. Then, at the hoc level each POINTER variable
that exists should be set up via the command:
setpointer pointer, variablewhere pointer and variable have enough implicit/explicit information to determine their exact segment and mechanism location. For a continuous mechanism, this means the section and location information. For a point process it means the object. The variable may also be any hoc variable or voltage, v.
For example, consider a synapse which requires a presynaptic potential in order to calculate the amount of transmitter release. Assume the declaration in the presynaptic model
NEURON { POINTPROCESS Syn POINTER vpre }
Then
objref syn
somedendrite {syn = new Syn(.8)}
setpointer syn.vpre, axon.v(1)
will allow the syn object to know the voltage at the distal end of the axon
section. As a variation on that example, if one supposed that the synapse
needed the presynaptic transmitter concentration (call it tpre) calculated
from a point process model called ``release'' (with object reference
rel, say) then the
statement would be
setpointer syn.tpre, rel.AcH_release
The caveat is that tight coupling between states in different models may cause numerical instability. When this happens, merging models into one larger model may eliminate the instability.