Chapter 9 How to expand NEURON's library of mechanisms Neuronal function involves the interaction of electrical and chemical signals that are distributed in time and space. The mechanisms that generate these signals and regulate their interactions are marked by a wide diversity of properties, differing across neuronal cell class, developmental stage, and species (e.g. chapter 7 in [Johnston, 1995 #218]; also see [McCormick, 1998 #288]). To be useful in research, a simulation environment must provide a flexible and powerful means for incorporating new biophysical mechanisms in models. It must also help the user remain focused on the model instead of programming. Such a means is provided to NEURON by NMODL, a high-level language that was originally implemented for NEURON by Michael Hines and later extended by him and Upinder Bhalla to generate code suitable for linking with GENESIS [Wilson, 1989 #289]. This chapter shows how to use NMODL to represent biophysical mechanisms by presenting a sequence of increasingly complex examples. Overview of NMODL A brief overview of how NMODL is used will clarify its rationale. First one writes a text file (a "mod file") that describes a mechanism as a set of nonlinear algebraic equations, differential equations, or kinetic reaction schemes. The description employs a syntax that closely resembles familiar mathematical and chemical notation. This text is passed to a translator that converts each statement into many statements in C, automatically generating code that handles details such as mass balance for each ionic species and producing code suitable for each of NEURON's integration methods. The output of the translator is then compiled for computational efficiency. This achieves tremendous conceptual leverage and savings of effort, not only because the high-level mechanism specification is much easier to understand and far more compact than the equivalent C code, but also because it spares the user from having to bother with low-level programming issues like how to "interface" the code with other mechanisms and with NEURON itself. NMODL is a descendant of the MOdel Description Language (MODL [Kohn, 1994 #290]), which was developed at Duke University by the National Biomedical Simulation Resource project for the purpose of building models that would be exercised by the Simulation Control Program (SCoP [Kootsey, 1986 #286]). NMODL has the same basic syntax and style of organizing model source code into named blocks as MODL. Variable declaration blocks, such as PARAMETER, STATE, and ASSIGNED, specify names and attributes of variables that are used in the model. Other blocks are directly involved in setting initial conditions or generating solutions at each time step (the equation definition blocks, e.g. INITIAL, BREAKPOINT, DERIVATIVE, KINETIC, FUNCTION, PROCEDURE). Furthermore, C code can be inserted inside the model source code to accomplish implementation-specific goals. NMODL recognizes all the keywords of MODL, but we will address only those that are relevant to NEURON simulations. We will also examine the changes and extensions that were necessary to endow NMODL with NEURON-specific features. To give these ideas real meaning, we will consider them in the context of models of the following kinds of mechanisms: a passive "leak" current and a localized transmembrane shunt (density mechanisms vs. point processes) an electrode stimulus (discontinuous parameter changes with variable time step methods) voltage-gated channels (differential equations vs. kinetic schemes) ion accumulation in a restricted space (extracellular K+) buffering, diffusion, and active transport (Ca2+ pump) Features of NMODL that are used in models of synaptic transmission and networks are examined in Chapter 10. Example 9.1: a passive "leak" current A passive "leak" current is one of the simplest biophysical mechanisms. Because it is distributed over the surface of a cell, it is described in terms of conductance per unit area and current per unit area, and therefore belongs to the class of "density" or "distributed mechanisms" (see Distributed mechanisms in Chapter 5). Figure 9.1 Figure 9.1 illustrates a branch of a neuron with a distributed leak current (left) and the equivalent circuit of a model of this passive current mechanism (right): a distributed constant conductance gleak in series with a voltage source Eleak that represents the equilibrium potential for the ionic current. The leak current density is given by ileak = gleak (Vm - Eleak), where Vm is membrane potential. Since this is a model of a physical system that is distributed in space, the variables ileak and Vm and the parameters gleak and Eleak are all functions of position. Listing 9.1 presents an implementation of this model with NMODL. Inline comments start with a colon and terminate at the end of the line. NMODL also allows comment blocks, which are demarcated by the keywords COMMENT . . . ENDCOMMENT. In passing, we mention that a similar syntax can be used to embed C code in a mod file, e.g. VERBATIM /* c statements */ ENDVERBATIM The statements between VERBATIM and ENDVERBATIM will appear without change in the output file that is written by the NMODL translator. Although this should be done only with great care, VERBATIM can be a convenient and effective way for individual users to add new features or even to employ NEURON as a "poor man's C compiler." : A passive leak current NEURON { SUFFIX leak NONSPECIFIC_CURRENT i RANGE i, e, g } PARAMETER { g = 0.001 (siemens/cm2) < 0, 1e9 > e = -65 (millivolt) } ASSIGNED { i (milliamp/cm2) v (millivolt) } BREAKPOINT { i = g*(v - e) } Listing 9.1. leak.mod Named blocks have the general form KEYWORD { statements }, and keywords are all upper case. User-defined variable names in NMODL can be up to 20 characters long. Each variable must be defined before it is used. The variable names chosen for this example were i, g, and e for the leak current, its specific conductance, and its equilibrium potential, respectively. Some variables are not "owned" by any mechanism but are available to all mechanisms; these include v, celsius, t, diam, and area. Another variable that is available to all mechanisms is dt. However, using dt in NMODL is neither necessary nor good practice. Before variable time step methods were added to NEURON, analytic expressions involving dt were often used for efficient modeling of voltage-sensitive channel states. This idiom is now built-in and employed automatically when such models are described in their underlying derivative form. The NEURON block The principal extension that differentiates NMODL from its MODL origins is that there are separate instances of mechanism data, with different values of states and parameters, in each segment (compartment) of a model cell. The NEURON block was introduced to make this possible by defining what the model of the mechanism looks like from the "outside" when many instances of it are sprinkled at different locations on the cell. The specifications entered in this block are independent of any particular simulator, but the detailed "interface code" requirements of a particular simulator determine whether the output C file is suitable for NEURON (NMODL) or GENESIS (GMODL). For this paper, we assume the translator is NMODL and that it produces code accepted by NEURON. The actual name of the current NMODL translator is nocmodl (nocmodl.exe on the PC). This translator is consistent with the object-oriented extensions that were introduced with version 3 of NEURON. However, the older translator which predated these extensions was called nmodl, and we will use the generic name NMODL to refer to NEURON-compatible translators. The SUFFIX keyword has two consequences. First, it identifies this to be a density mechanism, which can be incorporated into a NEURON cable section by an insert statement (see Usage below). Second, it tells the NEURON interpreter that the names for variables and parameters that belong to this mechanism will include the suffix _leak, so there will be no conflict with similar names in other mechanisms. The stipulation that i is a NONSPECIFIC_CURRENT also has two consequences. First, the value of i will be reckoned in charge balance equations. Second, this current will make no direct contribution to mass balance equations (it will have no direct effect on ionic concentrations). In later examples, we will see how to model mechanisms with specific ionic currents that can change concentrations. The RANGE keyword asserts that the values of i, e, and g are functions of position. That is, each of these variables can have a different value in each of the segments that make up a section. These variables can be accessed by the hoc interpreter using RANGE variable syntax (see Range and range variables in Chapter 5). The alternative to RANGE is GLOBAL, which is discussed below in The PARAMETER block. The membrane potential v is not mentioned in the NEURON block because it is one of the variables that are available to all mechanisms, and because it is a RANGE variable by default. However, for model completeness in non-NEURON contexts, and to enable units checking, v should be declared in the ASSIGNED block (see below). Variable declaration blocks As noted above, each user-defined variable must be declared before it is used. Even if it is named in the NEURON block, it still has to appear in a variable declaration block. Mechanisms frequently involve expressions that mix constants and variables whose units belong to different scales of investigation and which may themselves be defined in terms of other, more "fundamental" units. This can easily produce arithmetic errors that are difficult to isolate and rectify. Therefore NMODL has special provisions for establishing and maintaining consistency of units. To facilitate unit checking, each variable declaration includes a specification of its units in parentheses. The names used for these specifications are based on the UNIX units database, and are defined in a file called nrnunits.lib (located in nrn-5.4/share/lib/ under UNIX/Linux, and C:\nrn54\lib\ under MSWindows). A variable whose units are not specified is taken to be dimensionless. The user may specify whatever units are appropriate except for variables that are defined by NEURON itself. These include v (millivolts), t (milliseconds), celsius (şC), diam (µm), and area (µm2). Currents, concentrations, and equilibrium potentials created by the USEION statement also have Listing 9.1 specific units (see The NEURON block in Example 9.6: extracellular potassium accumulation below). In this particular density mechanism, i and g are given units of current per unit area (milliamperes/cm2) and conductance per unit area (siemens/cm2), respectively. The PARAMETER block Variables whose values are normally specified by the user are parameters, and are declared in a PARAMETER block. In the NEURON graphical user interface, a parameter is viewed with a special field editor which is designed to facilitate the entry of new values (see Usage below). PARAMETERs generally remain constant during a simulation, but they can be changed in mid-run if necessary to emulate some external influence on the characteristic properties of a model. To avoid confusion, such changes should only be performed through the hoc interpreter or the GUI, and not by statements in the mod file. The PARAMETER block in this example gives default values of 0.001 siemens/cm2 and -65 mV to g and e, respectively. The pair of values in angle brackets specifies the default minimum and maximum values for g that can be entered into the field editor of the GUI. In this case, we are trying to prevent negative values of conductance g. This protection, however, only holds for field editors and does not prevent a hoc statement from giving g a negative value. Because g and e are PARAMETERs, their values are visible at the hoc level and can be overridden by hoc commands or altered through the GUI. PARAMETERs ordinarily have global scope, which means that changing the value of a PARAMETER affects every instance of that mechanism throughout an entire model. However, the NEURON block for this particular mechanism stipulates that g and e are RANGE variables, so they can be given different values in every segment that has this leak current. The ASSIGNED block The ASSIGNED block is used to declare two kinds of variables: those that are given values outside the mod file, and those that appear on the left hand side of assignment statements within the mod file. The first group includes variables that are potentially available to every mechanism, such as v, celsius, t, and ionic variables (ionic variables are discussed in connection with The NEURON block in Example 9.6: extracellular potassium accumulation below). The second group specifically omits variables that are unknowns in a set of simultaneous linear or nonlinear algebraic equations, or that are dependent variables in differential equations or kinetic reaction schemes, which are handled differently (see Example 9.4: a voltage-gated current below for a discussion of the STATE block). Mechanism-specific ASSIGNED variables are RANGE variables by default. For a mechanism-specific ASSIGNED variable to be visible outside of the mod file, it must be declared as RANGE or GLOBAL in the NEURON block. ASSIGNED variables that are not "owned" by any mechanism (v, celsius, t, dt, diam, and area) are not mentioned in the NEURON block. The current i is not a state variable because the model of the leak current mechanism does not define it in terms of a differential equation or kinetic reaction scheme; that is to say, i has no dynamics of its own. Furthermore it is not an unknown in a set of equations. Instead, it is calculated by direct assignment. Therefore it is declared in the ASSIGNED block. For similar reasons membrane potential v is also declared in the ASSIGNED block. Although membrane potential is unquestionably a state variable in a model of a cell, to the leak current mechanism it is a driving force rather than a state variable. Equation definition blocks One equation suffices to describe this simple leak current model. This equation is defined in the BREAKPOINT block. As we shall see later, more complicated models may require invoking NMODL's built-in routines to solve families of simultaneous algebraic equations or perform numeric integration. The BREAKPOINT block The BREAKPOINT block is the main computation block in NMODL. Its name derives from SCoP, in which simulations are executed by incrementing an independent variable through a sequence of steps or "breakpoints" at which the dependent variables of the model are computed and displayed [Kohn, 1994 #290]. At exit from the BREAKPOINT block, all variables should be consistent with the independent variable. The independent variable in NEURON is always time t, and neither t nor the time step dt should be changed in NMODL. Usage The following hoc code illustrates how this mechanism might be used. Note the use of RANGE syntax to examine the value of i_leak near one end of cable. cable { nseg = 5 insert leak // override defaults g_leak = 0.002 // S/cm2 e_leak = -70 // mV } // show leak current density near 0 end of cable print cable.i_leak(0.1) Because of the interface code generated as a consequence of the definitions in the NEURON block, the leak mechanism will appear with the other density mechanisms in the Distributed Mechanism Manager and Viewer windows. This is illustrated in Fig. 9.2, which shows the Distributed Mechanism Inserter. The check mark signifies that the leak mechanism has been inserted into the section named cable. Figure 9.2. Compiling the leak mechanism automatically makes it available through NEURON's graphical user interface, as in this Distributed Mechanism Inserter (brought up by NEURON Main Menu / Tools / Distributed Mechanisms / Managers / Inserter). Example 9.2: a localized shunt At the opposite end of the spatial scale from a distributed passive current is a localized shunt induced by microelectrode impalement [Durand, 1984 #291][Staley, 1992 #151]. A shunt is restricted to a small enough region that it can be described in terms of a net conductance (or resistance) and total current, i.e. it is a point process (see Point processes in Chapter 5). Most synapses are also best represented by point processes. Figure 9.3 The localized nature of the shunt is emphasized in a cartoon of the neurite (Fig. 9.3 left). The equivalent circuit of the shunt (right) is similar to the equivalent circuit of the distributed leak current (Fig. 9.1 right), but here the resistance and current are understood to be concentrated in a single, circumscribed part of the cell. We will focus on how the NMODL code for this model (Listing 9.2) differs from the density mechanism presented earlier. The NEURON block The NEURON block identifies this mechanism as a point process, so it will be managed in hoc using an object-oriented syntax (see Usage below). Making i, e, and r RANGE variables means that each instance of this point process can have separate values for these variables. If a variable were asserted to be GLOBAL, then its value would be shared among all instances of the mechanism. Variable declaration blocks These are nearly identical to the PARAMETER and ASSIGNED blocks of the leak mechanism. However, Shunt is a point process so all of its current flows at one site instead of being distributed over an area. Therefore its i and r are in units of nanoamperes (total current) and gigaohms (0.001 ? total conductance in microsiemens), respectively. This code specifies default values for the PARAMETERs r and e. Allowing a minimum value of 10-9 for r prevents an inadvertent divide by 0 error (infinite conductance) by ensuring that a user cannot set r to 0 in its GUI field editor. However, as we noted in the leak model, the syntax does not prevent a hoc statement from assigning r a value outside of the desired range. : A shunt current NEURON { POINT_PROCESS Shunt NONSPECIFIC_CURRENT i RANGE i, e, r } PARAMETER { r = 1 (gigaohm) < 1e-9, 1e9 > e = 0 (millivolt) } ASSIGNED { i (nanoamp) v (millivolt) } BREAKPOINT { i = (0.001)*(v - e)/r } Listing 9.2. shunt.mod Equation definition blocks Like the leak current mechanism, the shunt mechanism is extremely simple and involves no state variables. The single equation is defined in the BREAKPOINT block. The BREAKPOINT block The sole "complication" here is that computation of i includes a factor of 0.001 to reconcile the units on the left and right hand sides of this assignment (nanoamperes vs. millivolts divided by gigaohms). The parentheses surrounding this conversion factor are a convention that is necessary for units checking: they disambiguate it from mere multiplication by a number. When NEURON's unit checking utility modlunit is used to check the NMODL code in Listing 9.2, it will find no errors and will exit without an error message. [ted@fantom dshunt]$ modlunit shunt.mod model 1.1.1.1 1994/10/12 17:22:51 Checking units of shunt.mod [ted@fantom dshunt]$ However if the parentheses were omitted, an error message would be emitted that reports inconsistent unit factors. [ted@fantom dshunt]$ modlunit shunt.mod model 1.1.1.1 1994/10/12 17:22:51 Checking units of shunt.mod The previous primary expression with units: 1-12 coul/sec is missing a conversion factor and should read: (0.001)*() at line 20 in file shunt.mod i = 0.001*(v - e)/r<> [ted@fantom dshunt]$ An error message would also result if parentheses surrounded a number which the user intended to be a quantity, since the unit factors would be inconsistent. The convention of using single numbers enclosed in parentheses to signify unit conversion factors is simple and minimizes the possibility of mistakes either by the user or by the software. It is important to note that expressions that involve more than one number, such as "(1 + 1)", will not be interpreted as conversion factors. Usage This hoc code illustrates how the shunt mechanism might be applied to a section called cable; note the object syntax for specifying the shunt resistance and current (see Point processes in Chapter 5). objref s // put near 0 end of cable cable s = new Shunt(0.1) // not bad for a sharp electrode s.r = 0.2 // show shunt current print s.i The definitions in the NEURON block of this particular model enable NEURON's graphical tools to include the Shunt object in the menus of its PointProcessManager and Viewer windows (Fig. 9.4). The check mark on the button adjacent to the numeric field for r indicates that the shunt resistance has been changed from its default value (0.2 gigaohm when the shunt was created by the hoc code immediately above) to 0.1 gigaohm. Figure 9.4. The properties of a specific instance of the Shunt mechanism are displayed in this Point Process Viewer (brought up by NEURON Main Menu / Tools / Point Processes / Viewers / PointProcesses / Shunt and then selecting Shunt[0] from the displayed list). Example 9.3: an intracellular stimulating electrode An intracellular stimulating electrode is similar to a shunt in the sense that both are localized sources of current that are modeled as point processes. However, the current from a stimulating electrode is not generated by an opening in the cell membrane but instead is injected directly into the cell. This particular model of a stimulating electrode (Listing 9.3) has the additional difference that the current changes discontinuously, i.e. it is a pulse with distinct start and stop times. The NEURON block This mechanism is identical to the built-in IClamp. Calling it IClamp1 allows the reader to test and modify it without conflict with the existing IClamp point process. This model of a current clamp generates a rectangular current pulse whose amplitude amp in nanoamperes, start time del in milliseconds, and duration dur in milliseconds are all adjustable by the user. Furthermore, these parameters are individually adjustable for each separate instance of this mechanism. Therefore they are declared as RANGE variables. The current i delivered by IClamp1 is declared in the NEURON block to make it available for examination. The ELECTRODE_CURRENT statement has two important consequences: positive values of i will depolarize the cell (in contrast to the hyperpolarizing effect of positive transmembrane currents), and when the extracellular mechanism is present there will be a change in the extracellular potential vext. Extracellular fields are discussed further in Chapter X. Equation definition blocks The BREAKPOINT block The logic for deciding whether i = 0 or i = amp is straightforward, but the at_time() calls need explanation. To work properly with variable time step methods, e.g. CVODE, models that change parameters discontinuously during a simulation must notify NEURON when such events take place. With fixed time step methods, users implicitly assume that events take place on time step boundaries (integer multiples of dt), and they would never consider defining a pulse duration narrower than dt. Neither eventuality can be left to chance with variable time step methods. During a variable time step simulation, the first at_time() call guarantees that a time step boundary will be at del - ?, where ? is on the order of 10-9 ms. Integration will restart from its new initial condition at del + ?. Here we must point out that at_time() has become a "deprecated" function, i.e. it still works but it should not be used in future model development. It is important to understand what at_time() does because it has been used in legacy code, and we discuss it in further detail later in this chapter under the heading Models with discontinuities. However, NEURON's event delivery system (see Chapter 10) provides a far better way to implement discontinuities. : Current clamp NEURON { POINT_PROCESS IClamp1 RANGE del, dur, amp, i ELECTRODE_CURRENT i } UNITS { (nA) = (nanoamp) } PARAMETER { del (ms) dur (ms) < 0, 1e9 > amp (nA) } ASSIGNED { i (nA) } INITIAL { i = 0 } BREAKPOINT { at_time(del) at_time(del+dur) if (t < del + dur && t > del) { i = amp } else { i = 0 } } Listing 9.3. iclamp1.mod The INITIAL block The code in the INITIAL block is executed when the hoc function finitialize() is called. The initialization here consists of making sure that IClamp1.i is 0 when t = 0. Initialization of more complex mechanisms is discussed below in Example 9.4: a voltage-gated current and Example 9.6: extracellular potassium accumulation, and Chapter 8 considers the topic of initialization from a broader perspective. Usage Regardless of whether a fixed or variable time step integrator is chosen, IClamp1 looks the same to the user. In either case, a current stimulus of 0.01 nA amplitude that starts at t = 1 ms and lasts for 2 ms would be created by this hoc code or through the GUI (Fig. 9.5). objref ccl // put at middle of soma soma ccl = new IClamp1(0.5) ccl.del = 1 ccl.dur = 2 ccl.amp = 0.01 Figure 9.5. A PointProcessManager configured as an Iclamp1 object. Example 9.4: a voltage-gated current One of the particular strengths of NMODL is its flexibility in dealing with ion channels whose conductances are not constant but instead are regulated by factors such as membrane potential and/or ligand concentrations on one or both sides of the membrane. Here we use the well-known Hodgkin-Huxley (HH) delayed rectifier to show how a voltage-gated current can be implemented, and later we will examine a model of a potassium current that depends on both voltage and intracellular calcium concentration. The delayed rectifier and all other voltage-gated channels that are distributed over the cell surface are distributed mechanisms. Therefore their NMODL representations and hoc usage be similar to those of the passive leak current presented in Example 9.1. The following discussion focuses on the significant differences between the models of the delayed rectifier and the passive leak current. In this example, membrane potential is in absolute millivolts, i.e. reversed in polarity from the original Hodgkin-Huxley convention and shifted to reflect a resting potential of -65 mV. : HH voltage-gated potassium current NEURON { SUFFIX kd USEION k READ ek WRITE ik RANGE gbar, g, i } UNITS { (S) = (siemens) (mV) = (millivolt) (mA) = (milliamp) } PARAMETER { gbar = 0.036 (S/cm2) } ASSIGNED { v (mV) ek (mV) : typically ~ -77.5 ik (mA/cm2) i (mA/cm2) g (S/cm2) } STATE { n } BREAKPOINT { SOLVE states METHOD cnexp g = gbar * n^4 i = g * (v - ek) ik = i } INITIAL { : Assume v has been constant for a long time n = alpha(v)/(alpha(v) + beta(v)) } DERIVATIVE states { : Computes state variable n at present v & t n' = (1-n)*alpha(v) - n*beta(v) } FUNCTION alpha(Vm (mV)) (/ms) { LOCAL x UNITSOFF x = (Vm+55)/10 if (fabs(x) > 1e-6) { alpha = 0.1*x/(1 - exp(-x)) } else { alpha = 0.1/(1 - 0.5*x) } UNITSON } FUNCTION beta(Vm (mV)) (/ms) { UNITSOFF beta = 0.125*exp(-(Vm+65)/80) UNITSON } Listing 9.4. kd.mod The NEURON block As with the passive leak model, SUFFIX marks this as a density mechanism, whose variables and parameters are identified in hoc by a particular suffix. Three RANGE variables are declared in this block: the peak conductance density gbar (the product of channel density and "open" conductance per channel), the macroscopic conductance g (the product of gbar and the fraction of channels that are open at any moment), and the current i that passes through g. At the level of hoc, these will be available as gbar_kd, g_kd, and i_kd. This model also has a fourth RANGE variable: the gating variable n, which is declared in the STATE block (see The STATE block below). STATE variables are automatically RANGE variables and do not need to be declared in the NEURON block. A mechanism needs a separate USEION statement for each of the ions that it affects or is affected by. This example has one USEION statement, which includes READ ek because the potential gradient that drives i_kd depends on the equilibrium potential for K+ (potassium). Since the resulting ionic flux may change local [K+], this example also includes WRITE ik. The WRITE ionname syntax enables NEURON to keep track of the total outward current that is carried by an ionic species, its internal and external concentrations, and its equilibrium potential. We will return to this point in the context of a model with extracellular K+ accumulation. The UNITS block The statements in the UNITS block define new names for units in terms of existing names in the UNIX units database. This can increase legibility and convenience, and is helpful both as a reminder to the user and as a means for automating the process of checking for consistency of units. Variable declaration blocks The ASSIGNED block This is analogous to the ASSIGNED block of the leak mechanism. For the sake of clarity, variables whose values are computed outside this mod file are listed first. Note that ek is listed as an ASSIGNED variable, unlike e of the leak mechanism which was a PARAMETER. The reason for this difference is that mechanisms that produce K+ fluxes may cause the equilibrium potential ek to change in the course of a simulation. However, the equilibrium potential for the leak current was not linked to a specific ionic species and therefore remains fixed unless explicitly altered by hoc statements or the GUI. The STATE block If a model involves differential equations, families of algebraic equations, or kinetic reaction schemes, their dependent variables or unknowns are to be listed in the STATE block. Therefore gating variables such as the delayed rectifier's n are declared here. In NMODL, variables that are declared in the STATE block are called STATE variables, or simply STATEs. This NMODL-specific terminology should not be confused with the physics or engineering concept of a "state variable" as a variable that describes the state of a system. While membrane potential is a "state variable" in the engineering sense, it would never be a STATE because its value is calculated only by NEURON and never by NMODL code. Likewise, the unknowns in a set of simultaneous equations (e.g. specified in a LINEAR or NONLINEAR block) would not be state variables in an engineering sense, yet they would all be STATEs. All STATEs are automatically RANGE variables. This is appropriate, since channel gating can vary with position along a neurite. Equation definition blocks In addition to the BREAKPOINT block, this model also has INITIAL, DERIVATIVE, and FUNCTION blocks. The BREAKPOINT block This is the main computation block of the mechanism. By the end of the BREAKPOINT block, all variables are consistent with the new time. If a mechanism has STATEs, this block must contain one SOLVE statement that tell how the values of the STATEs will be computed over each time step. The SOLVE statement specifies a block of code that defines the simultaneous equations that govern the STATEs. Currents are set with assignment statements at the end of the BREAKPOINT block. There are two major reasons why variables that depend on the number of executions, such as counts or flags or random numbers, should generally not be calculated in a BREAKPOINT block. First, the assignment statements in a BREAKPOINT block are usually called twice per time step. Second, with variable time step methods the value of t may not even be monotonically increasing. The right way to think about this is to remember that the BREAKPOINT block is responsible for making all variables consistent at time t. Thus assignment statements in this block are responsible for trivially specifying the values of variables that depend only on the values of STATEs, t, and v, while the SOLVE statements perform the magic required to make the STATEs consistent at time t. It is not belaboring the point to reiterate that the assignment statements should produce the same result regardless of how many times BREAKPOINT is called with the same STATEs, t, and v. All too often, errors have resulted from an attempt to explicitly compute what is conceptually a STATE in a BREAKPOINT block. Computations that must be performed only once per time step should be placed in a PROCEDURE, which in turn would be invoked by a SOLVE statement in a BREAKPOINT block. We must also emphasize that the SOLVE statement is not a function call, and that the body of the DERIVATIVE block (or any other block specified in a SOLVE statement) will be executed asynchronously with respect to BREAKPOINT assignment statements. Therefore it is incorrect to invoke rate functions from the BREAKPOINT block; instead these must be called from the block that is specified by the SOLVE statement (in this example, from within the DERIVATIVE block). Models of active currents such as i_kd are generally formulated in terms of ionic conductances that are functions of voltage- and time-dependent gating variables. The SOLVE statements at the beginning of the BREAKPOINT block specify the differential equations or kinetic schemes that govern the kinetics of the gating variables. The algebraic equations that compute the ionic conductances and currents follow the SOLVE statements. For mechanisms whose STATEs are described by differential equations, it is often most convenient and efficient to use one of NEURON's built-in numerical integrators. A good choice for this particular mechanism is cnexp, which is described below in connection with the DERIVATIVE block. The INITIAL block The INITIAL block may contain any instructions that should be executed when the hoc function finitialize() is called. Though often overlooked, proper initialization of all STATEs is as important as correctly computing their temporal evolution. This is accomplished for the common case by finitialize(), which executes the initialization strategy defined in the INITIAL block for each mechanism (see also INITIAL blocks in NMODL in Chapter 8). Prior to executing the INITIAL block, STATE values are set to their values in the STATE declaration block (or set to 0 if it was not given a specific value in the STATE declaration block). For this delayed rectifier mechanism, n is set to its steady-state value for the membrane potential that exists in the compartment. This potential itself can be "left over" from a previous simulation run, or it can be specified by the user, e.g. on a compartment by compartment basis using statements such as dend.v(0.2) = -48 before calling finitialize(), or uniformly over the entire cell with a statement like finitialize(-55). The INITIAL block should be used to initialize STATEs with respect to the initial values of membrane potential and ionic concentrations. There are several other ways to prepare STATEs for a simulation run, (see Chapter 8 for details, especially Examples of custom initializations). Among these, the most direct is simply to assign values explicitly using hoc statements such as cable.n_kd(0.3) = 0.9, but this can create arbitrary initial conditions that would be quite unnatural. A more "physiological" approach, which may be appropriate for models of oscillating or chaotic systems or whose mechanisms show other complex interactions, is to perform an initialization run during which the model converges toward its limit cycle or attractor. A practical alternative for systems that settle to a stable equilibrium point when left undisturbed is to assign t a large negative value and then advance the simulation over several large time steps (keeping t < 0 prevents the initialization steps from triggering scheduled events such as stimulus currents or synaptic inputs); this takes advantage of the strong stability properties of NEURON's implicit integration methods. The DERIVATIVE block This is used to assign values to the derivatives of those STATEs that are described by differential equations. The statements in this block are of the form y? = expr, where a series of apostrophes can be used to signify higher-order derivatives. For NEURON's fixed time step integration method, these equations are integrated using the numerical method specified by the SOLVE statement in the BREAKPOINT block. The SOLVE statement should explicitly invoke one of the integration methods that is appropriate for systems in which state variables can vary widely during a time step (stiff systems). The cnexp method used in this example combines second order accuracy with computational efficiency. It is appropriate when the right hand side of y? = f(v,y) is linear in y and involves no other states, so it is well-suited to models with HH-style ionic currents. This method calculates the STATEs analytically under the assumption that all other variables are constant throughout the time step. If the variables change but are second order correct at the midpoint of the time step, then the calculation of STATEs is also second order correct. If f(v,y) is not linear in y, then the SOLVE statement in the BREAKPOINT block should specify the implicit integration method derivimplicit. This provides first-order accuracy and is usable with general ODEs regardless of stiffness or nonlinearity. With variable time step methods, no variable is assumed to be constant. These methods not only change the time step, but adaptively choose a numerical integration formula with local error that ranges from first-order up to O(?t6). The present implementation of NMODL creates a diagonal Jacobian approximation for the block of STATEs. This is done analytically if yi? = fi(v,y) is polynomial in yi; otherwise, the Jacobian is approximated by numerical differencing. In the rare case where this is inadequate, the user may supply an explicit Jacobian. Future versions of NMODL may attempt to deal with Jacobian evaluation in a more sophisticated manner. This illustrates a particularly important benefit of the NMODL approach: improvements in methods do not affect the high level description of the membrane mechanism. The FUNCTION block The functions defined by FUNCTION blocks are available at the hoc level and in other mechanisms by adding the suffix of the mechanism in which they are defined, e.g. alpha_kd() and beta_kd(). Functions or procedures can be simply called from hoc if they do not reference RANGE variables (references to GLOBAL variables are allowed). If a function or procedure does reference a RANGE variable, then prior to calling the function from hoc it is necessary to specify the proper instance of the mechanism (its location on the cell). This is done by a setdata_ function that has the syntax section_name setdata_suffix(x) where section_name is the name of the section that contains the mechanism in question, suffix is the mechanism suffix, and x is the normalized distance along the section where the particular instance of the mechanism exists. The functions in our kd example do not use RANGE variables, so a specific instance is not needed. The differential equation that describes the kinetics of n involves two voltage-dependent rate constants whose values are computed by the functions alpha() and beta(). The original algebraic form of the equations that define these rates is and Eq. 9.1 The denominator for goes to 0 when v = -55 mV, which could cause numeric overflow. The code used in alpha() avoids this by switching, when v is very close to -55, to an alternative expression that is based on the first three terms of the infinite series expansion of ex. As noted elsewhere in this paper, NMODL has features that facilitate establishing and maintaining consistency of units. Therefore the rate functions alpha() and beta() are introduced with the syntax FUNCTION f_name(arg1 (units1), arg2 (units2), . . . )(returned_units) to declare that their arguments are in units of millivolts and that their returned values are in units of inverse milliseconds ("/ms"). This allows automatic units checking on entry to and return from these functions. For the sake of legibility the UNITSOFF . . . UNITSON directives disable units checking just within the body of these functions. This is acceptable because the terms in the affected statements are mutually consistent. Otherwise the statements would have to be rewritten in a way that makes unit consistency explicit at the cost of legibility, e.g. x = (Vm + 55 (millivolt))/(10 (millivolt)) Certain variables exist solely for the sake of computational convenience. These typically serve as scale factors, flags, or temporary storage for intermediate results, and are not of primary importance to the mechanism. Such variables are often declared as LOCAL variables within an equation block, e.g. x in this mechanism. LOCAL variables that are declared in an equation block are not "visible" outside the block and they do not retain their values between invocations of the block. LOCAL variables that are declared outside an equation block have very different properties and are discussed under Variable declaration blocks in Example 9.8: calcium diffusion with buffering. Usage The hoc code and graphical interface for using this distributed mechanism are similar to those for the leak mechanism (Fig. 9.2). However, the kd mechanism involves more RANGE variables, and this is reflected in the choices available in the RANGE variable menu of NEURON's Plot what? tool for graph windows. Since kd uses potassium, the variables ek and ik (total K+ current) appear in this list along with the variables that are explicitly declared as RANGE and STATE in kd.mod (see Fig. 9.6). The total K+ current ik will differ from i_kd only if another mechanism that WRITEs ik is present in this section. Figure 9.6. A Plot what? tool from a graph created after the kd mechanism was inserted into a section called cable. Note the hoc names of variables associated with the kd mechanism. Example 9.5: a calcium-activated voltage-gated current This model of a potassium current that depends on both voltage and intracellular calcium concentration [Ca2+]i. is based on the work of Moczydlowski and Latorre [Moczydlowski, 1983 #306]. It is basically an elaboration of the HH mechanism in which the forward and backward rates depend jointly on membrane potential and [Ca2+]i. Here we point out the salient implementational differences between this and the previous model. : Calcium activated K channel NEURON { SUFFIX cagk USEION ca READ cai USEION k READ ek WRITE ik RANGE gkbar GLOBAL oinf, tau } UNITS { (mV) = (millivolt) (mA) = (milliamp) (S) = (siemens) (molar) = (1/liter) (mM) = (millimolar) FARADAY = (faraday) (kilocoulombs) R = (k-mole) (joule/degC) } PARAMETER { gkbar = 0.01 (S/cm2) d1 = 0.84 d2 = 1.0 k1 = 0.18 (mM) k2 = 0.011 (mM) bbar = 0.28 (/ms) abar = 0.48 (/ms) } ASSIGNED { cai (mM) : typically 0.001 celsius (degC) : typically 20 v (mV) ek (mV) ik (mA/cm2) oinf tau (ms) } STATE { o } : fraction of channels that are open BREAKPOINT { SOLVE state METHOD cnexp ik = gkbar*o*(v - ek) } DERIVATIVE state { rate(v, cai) o' = (oinf - o)/tau } INITIAL { rate(v, cai) o = oinf } : the following are all callable from hoc FUNCTION alp(v (mV), ca (mM)) (/ms) { alp = abar/(1 + exp1(k1,d1,v)/ca) } FUNCTION bet(v (mV), ca (mM)) (/ms) { bet = bbar/(1 + ca/exp1(k2,d2,v)) } FUNCTION exp1(k (mM), d, v (mV)) (mM) { : numeric constants in an addition or subtraction : expression automatically take on the unit values : of the other term exp1 = k*exp(-2*d*FARADAY*v/R/(273.15 + celsius)) } PROCEDURE rate(v (mV), ca (mM)) { LOCAL a : LOCAL variable takes on units of right hand side a = alp(v,ca) tau = 1/(a + bet(v, ca)) oinf = a*tau } Listing 9.5. cagk.mod The NEURON block This potassium conductance depends on [Ca2+]i, so two USEION statements are required. The RANGE statement declares the peak conductance density gkbar, but not g, so this mechanism's ionic conductance will not be visible from hoc (in fact, this model doesn't even calculate the activated ionic conductance density). Likewise, there is no i_cagk to report this particular current component separately, even though it will be added to the total K+ current ik because of WRITE ik. The variables oinf and tau, which govern the gating variable o, should be accessible in hoc for the purpose of seeing how they vary with membrane potential and [Ca2+]i. At the same time, the storage and syntax overhead required for a RANGE variable does not seem warranted because it appears unlikely to be necessary or useful to plot either oinf or tau as a function of space. Therefore they have been declared GLOBAL rather than RANGE. On first examination, this might seem to pose a problem. The gating of this K+ current depends on membrane potential and [Ca2+]i, both of which may vary with location, so how can it be correct to use GLOBALs for oinf and tau? And if some reason did arise to examine the values of these variables at a particular location, how could this be done? The answers to these questions lie in the DERIVATIVE and PROCEDURE blocks, as we shall see below. The UNITS block The last two statements in this block require some clarification. The first parenthesized item on the right hand side of the equal sign is the numeric value of a standard entry in the UNIX units database, which may be expressed on a scale appropriate for physics rather than membrane biophysics. The second parenthesized item acts like a scale factor that converts it to the specific units chosen for this model. Thus (faraday) appears in the units database in terms of coulombs/mole and has a numeric value of 96,485.309, but for this particular mechanism we prefer to use a constant whose units are kilocoulombs/mole. The statement FARADAY = (faraday) (kilocoulombs) results in FARADAY having units of kilocoulombs and a numeric value of 96.485309. The item (k-mole) in the statement R = (k-mole) (joule/degC) is not kilomoles but instead is a specific entry in the units database equal to the product of Boltzmann's constant and Avogadro's number. The end result of this statement is that R has units of joules/°C and a numeric value of 8.313424. These special definitions of FARADAY and R pertain to this mechanism only; a different mechanism could assign different units and numeric values to these labels. Another possible source of confusion is the interpretation of the symbol "e". This is always the electronic charge (~ 1.6 · 10-19 coulombs), except outside the UNITS block where a single number in parentheses is treated as a conversion factor, e.g. the expression (2e4) is treated as a conversion factor of 2 · 104. Errors involving "e" in a units expression are easy to make, but they are always caught by modlunit. Variable declaration blocks The ASSIGNED block Comments in this block can be helpful to the user as reminders of "typical" values or usual conditions under which a mechanism operates. For example, the cagk mechanism is intended for use in the context of [Ca2+]i on the order of 0.001 mM. Similarly, the temperature sensitivity of this mechanism is accommodated by including the global variable celsius. NEURON's default value for celsius is 6.3şC, but as the comment in this mod file points out, the parameter values for this particular mechanism were intended for an "operating temperature" of 20şC. Therefore the user may need to change celsius through hoc or the GUI. The variables oinf and tau, which were made accessible to NEURON by the GLOBAL statement in the NEURON block, are given values by the procedure rate and are declared as ASSIGNED. The STATE block This mechanism needs a STATE block because o, the fraction of channels that are open, is described by a differential equation. Equation definition blocks The BREAKPOINT block This mechanism does not make its ionic conductance available to hoc, so the BREAKPOINT block just calculates the ionic current passing through these channels and doesn't bother with separate computation of a conductance. The DERIVATIVE block The gating variable o is governed by a first-order differential equation. The procedure rate() assigns values to the voltage-sensitive parameters of this equation: the steady-state value oinf, and the time constant tau. This answers the first question that was raised above in the discussion of the NEURON block. The procedure rate will be executed individually for each segment in the model that has the cagk mechanism. Each time rate is called, its arguments will equal the membrane potential and [Ca2+]i of the segment that is being processed, since v and cai are both RANGE variables. Therefore oinf and tau can be GLOBAL without destroying the spatial variation of the gating variable o. The FUNCTION and PROCEDURE blocks The functions alp(), bet(), exp1(), and the procedure rate() implement the mathematical expressions that describe oinf and tau. To facilitate units checking, their arguments are tagged with the units that they use. The rate() procedure achieves some efficiency by calling alp() once and using the returned value twice; calculating oinf and tau separately would have required two calls to alp(). The procedure rate() helps answer the second question that was raised in the discussion of the NEURON block: how to examine the variation of oinf and tau over space. This is easily done in hoc with nested loops like this: forall { // iterate over all sections for (x) { // iterate over each segment rate(v(x), cai(x)) // here put statements to plot // or save oinf and tau } } Usage This mechanism involves both K+ and Ca2+, so the list of RANGE variables displayed by Plot what? has more entries than it did for the kd mechanism (Fig. 9.7; compare this with Fig. 9.6). However, cai, cao, and eca will remain constant unless the section in which this mechanism has been inserted also includes something that can affect calcium concentration (e.g. a pump or buffer). Figure 9.7. A Plot what? tool from a graph created after the cagk mechanism was inserted into a section called cable. Note the hoc names of variables associated with the cagk mechanism. Example 9.6: extracellular potassium accumulation Because mechanisms can generate transmembrane fluxes that are attributed to specific ionic species by the USEION x WRITE ix syntax, modeling the effects of restricted diffusion is straightforward. The kext mechanism described here emulates potassium accumulation in the extracellular space adjacent to squid axon (Fig. 9.8). The experiments of Frankenhaeuser and Hodgkin [Frankenhaeuser, 1956 #307] indicated that satellite cells and other extracellular structures act as a diffusion barrier that prevents free communication between this space and the bath. When there is a large efflux of K+ ions from the axon, e.g. during the repolarizing phase of an action potential or in response to injected depolarizing current, K+ builds up in this "Frankenhaeuser-Hodgkin space" (F-H space). This elevation of [K+]o shifts EK in a depolarized direction, which has two important consequences. First, it reduces the driving force for K+ efflux and causes a decline of the outward IK. Second, when the action potential terminates or the injected depolarizing current is stopped, the residual elevation of [K+]o causes an inward current that decays gradually as [K+]o equilibrates with [K+]bath. Figure 9.8. Restricted diffusion may cause extracellular potassium accumulation adjacent to the cell membrane. From Fig. 1 in [Hines, 2000 #323]. : Extracellular potassium ion accumulation NEURON { SUFFIX kext USEION k READ ik WRITE ko GLOBAL kbath RANGE fhspace, txfer } UNITS { (mV) = (millivolt) (mA) = (milliamp) FARADAY = (faraday) (coulombs) (molar) = (1/liter) (mM) = (millimolar) } PARAMETER { kbath = 10 (mM) : seawater (squid axon!) fhspace = 300 (angstrom) : effective thickness of F-H space txfer = 50 (ms) : tau for F-H space <-> bath exchange = 30-100 } ASSIGNED { ik (mA/cm2) } STATE { ko (mM) } BREAKPOINT { SOLVE state METHOD cnexp } DERIVATIVE state { ko' = (1e8)*ik/(fhspace*FARADAY) + (kbath - ko)/txfer } Listing 9.6. kext.mod The NEURON block A compartment may contain several mechanisms that have direct interactions with ionic concentrations (e.g. diffusion, buffers, pumps). Therefore NEURON must be able to compute the total currents and concentrations consistently. The USEION statement sets up the necessary "bookkeeping" by automatically creating a separate mechanism that keeps track of four essential variables: the total outward current carried by an ion, the internal and external concentrations of the ion, and its equilibrium potential. In this case the name of the ion is "k" and the automatically-created mechanism is called "k_ion" in the hoc interpreter. The k_ion mechanism has variables ik, ki, ko, and ek, which represent IK, [K+]i, [K+]o, and EK, respectively. These do not have suffixes; furthermore, they are RANGE variables so they can have different values in every segment of each section in which they exist. In other words, the K+ current through Hodgkin-Huxley potassium channels near one end of the section cable is cable.ik_hh(0.1), but the total K+ current generated by all sources, including other ionic conductances and pumps, would be cable.ik(0.1). This mechanism computes [K+]o from the outward potassium current, so it READs ik and WRITEs ko. When a mechanism WRITEs a particular ionic concentration, this means that it sets the value for that concentration at all locations in every section into which it has been inserted. This has an important consequence: in any given section, no ionic concentration should be "written" by more than one mechanism. The bath is assumed to be a large, well-stirred compartment that envelops the entire "experimental preparation." Therefore kbath is a GLOBAL variable so that all sections that contain the kext mechanism will have the same numeric value for [K+]bath. Since this would be one of the controlled variables in an experiment, the value of kbath is specified by the user and remains constant during a simulation. The thickness of the F-H space is fhspace, the time constant for equilibration with the bath is txfer, and both are RANGE variables so they can vary along the length of each section. Variable declaration blocks The PARAMETER block The default value of kbath is set to 10 mM, consistent with the composition of seawater [Frankenhaeuser, 1956 #307]. Since kbath is GLOBAL, a single hoc statement can change this to a new value that will affect all occurrences of the kext mechanism, e.g. kbath_kext = 8 would change it to 8 mM everywhere. The STATE block Ionic concentration is a STATE of a mechanism only if that mechanism calculates the concentration in a DERIVATIVE or KINETIC block. This model computes ko, the potassium concentration in the F-H space, according to the dynamics specified by an ordinary differential equation. Equation definition blocks The BREAKPOINT block This mechanism involves a single differential equation that tells the rate of change of ko, the K+ concentration in the F-H space. The choice of integration method in NMODL is based on the fact that the equation is linear in ko. The total K+ current ik might also vary during a time step (see the DERIVATIVE block) if membrane potential, some K+ conductance, or ko itself is changing rapidly. In a simulation where such rapid changes are likely to occur, proper modeling practice would lead one either to use NEURON with CVODE, or to use a fixed time step that is short compared to the rate of change of ik. The INITIAL block How to provide for initialization of variables is a recurring question in model implementation, and here it comes again. The answer is important because it bears directly on how the model will be used. The only STATE in this mechanism is the ionic concentration ko, which we could initialize in several different ways. The simplest might be with the INITIAL block INITIAL { ko = kbath } but this seems too limiting. One alternative is to declare a new RANGE variable ko0 in the NEURON block, specify its value in the PARAMETER block PARAMETER { . . . ko0 = 10 (mM) = 10 (mM) } and use this INITIAL block INITIAL { ko = ko0 } This would be a very flexible implementation, allowing ko0 to vary with location wherever kext has been inserted. But some care is needed in its use, because ion concentration assignment in an INITIAL block can result in an inconsistent initialization on return from finitialize() (see Ion concentrations and equilibrium potentials in Chapter 8). So for this example we have decided to let the initial value of ko be controlled by the built-in hoc variable ko0_k_ion (see Initializing concentrations in hoc in Chapter 8). To make our mechanism rely on ko0_k_ion for the initial value of ko, we merely omit any ko = . . . assignment statement from the INITIAL block. Since ko is kext's only STATE, we don't need an INITIAL block at all. This might seem a less flexible approach than using our own ko0 RANGE variable, because ko0_k_ion is a global variable (has the same value wherever ko is defined), but Initializing concentrations in hoc in Chapter 8 shows how to work around this apparent limitation. The DERIVATIVE block At the core of this mechanism is a single differential equation that relates d[K+]o/dt to the sum of two terms. The first term describes the contribution of ik to [K+]o, subject to the assumption that the thickness F-H space is much smaller than the diameter of the section. The unit conversion factor of 108 is required because fhspace is given in Ċngstroms. The second term describes the exchange of K+ between the bath and the F-H space. Usage If this mechanism is present in a section, the following RANGE variables will be accessible through hoc: [K+] inside the cell and within the F-H space (ki and ko); equilibrium potential and total current for K (ek and ik); thickness of the F-H space and the rate of equilibration between it and the bath (fhspace_kext and txfer_kext). The bath [K+] will also be available as the global variable kbath_kext. General comments about kinetic schemes Kinetic schemes provide a high level framework that is perfectly suited for compact and intuitively clear specification of models that involve discrete states in which "material" is conserved. The basic notion is that flow out of one state equals flow into another (also see Chemical reactions in Chapter 3). Almost all models of membrane channels, chemical reactions, macroscopic Markov processes, and ionic diffusion are elegantly expressed through kinetic schemes. It will be helpful to review some fundamentals before proceeding to specific examples of mechanisms implemented with kinetic schemes. The unknowns in a kinetic scheme, which are usually concentrations of individual reactants, are declared in the STATE block. The user expresses the kinetic scheme with a notation that is very similar to a list of simultaneous chemical reactions. The NMODL translator converts the kinetic scheme into a family of ODEs whose unknowns are the STATEs. Hence STATE { mc m } KINETIC scheme1 { ~ mc <-> m (a(v), b(v)) } is equivalent to DERIVATIVE scheme1 { mc' = -a(v)*mc + b(v)*m m' = a(v)*mc - b(v)*m } The first character of a reaction statement is the tilde "~", which is used to immediately distinguish this kind of statement from other sequences of tokens that could be interpreted as an expression. The expression to the left of the three character reaction indicator "<->" specifies the reactants, and the expression immediately to the right specifies the products. The two expressions in parentheses are the forward and reverse reaction rates (here the rate functions a(v) and b(v)). Immediately after each reaction, the variables f_flux and b_flux are assigned the values of the forward and reverse fluxes respectively. These can be used in assignment statements such as ~ cai + pump <-> capump (k1,k2) ~ capump <-> pump + cao (k3,k4) ica = (f_flux - b_flux)*2*Faraday/area In this case, the forward flux is k3*capump, the reverse flux is k4*pump*cao, and the positive-outward current convention is consistent with the sign of the expression for ica (in the second reaction, forward flux means positive ions move from the inside to the outside). More complicated reaction sequences such as the wholly imaginary KINETIC scheme2 { ~ 2A + B <-> C (k1,k2) ~ C + D <-> A + 2B (k3,k4) } begin to show the clarity of expression and suggest the comparative ease of modification of the kinetic representation over the equivalent but stoichiometrically confusing DERIVATIVE scheme2 { A' = -2*k1*A^2*B + 2*k2*C + k3*C*D - k4*A*B^2 B' = -k1*A^2*B + k2*C + 2*k3*C*D - 2*k4*A*B^2 C' = k1*A^2*B - k2*C - k3*C*D + k4*A*B^2 D' = - k3*C*D + k4*A*B^2 } Clearly a statement such as ~ calmodulin + 3Ca <-> active (k1, k2) would be easier to modify (e.g. so it requires combination with 4 calcium ions) than the relevant term in the three differential equations for the STATEs that this reaction affects. The kinetic representation is easy to debug because it closely resembles familiar notations and is much closer to the conceptualization of what is happening than the differential equations would be. Another benefit of kinetic schemes is the simple polynomial nature of the flux terms, which allows the translator to easily perform a great deal of preprocessing that makes implicit numerical integration more efficient. Specifically, the nonzero elements (partial derivatives of with respect to ) of the sparse matrix are calculated analytically in NMODL and collected into a C function that is called by solvers to calculate the Jacobian. Furthermore, the form of the reaction statements determines if the scheme is linear, obviating an iterative computation of the solution. Voltage-sensitive rates are allowed, but to guarantee numerical stability the rate constants should not be functions of STATEs. Thus writing the calmodulin example as ~ calmodulin <-> active (k3*Ca^3, k2) will work but is potentially unstable if Ca is a STATE in other simultaneous reactions in the same mod file. Variable time step methods such as CVODE will compensate by reducing dt, but this will make the simulation run more slowly. Kinetic scheme representations provide a great deal of leverage because a single compact expression is equivalent to a large amount of C code. One special advantage from the programmer's point of view is the fact that these expressions are independent of the solution method. Different solution methods require different code, but the NMODL translator generates this code automatically. This saves the user's time and effort and ensures that all code expresses the same mechanism. Another advantage is that the NMODL translator handles the task of interfacing the mechanism to the remainder of the program. This is a tedious exercise that would require the user to have special knowledge that is not relevant to neurophysiology and which may change from version to version. Special issues are raised by mechanisms that involve fluxes between compartments of different size, or whose reactants have different units. The first of the following examples has none of these complications, which are addressed later in models of diffusion and active transport. Example 9.7: kinetic scheme for a voltage-gated current This illustration of NMODL's facility for handling kinetic schemes implements a simple three-state model for the conductance state transitions of a voltage-gated potassium current Eq. 9.2 The closed states are C1 and C2, the open state is O, and the rates of the forward and backward state transitions are calculated in terms of the equilibrium constants and time constants of the isolated reactions through the familiar expressions and . The equilibrium constants are the Boltzmann factors and , where the energies of states C1, C2, and O are 0, , and respectively. The typical sequence of analysis is to determine the constants k1, d1, k2, and d2 by fitting the steady-state voltage clamp data, and then to find the voltage-sensitive transition time constants and from the temporal properties of the clamp current at each voltage pulse level. In this example the steady-state information has been incorporated in the NMODL code, and the time constants are conveyed by tables (arrays) that are created within the interpreter. : Three state kinetic scheme for HH-like potassium channel : Steady-state v-dependent state transitions have been fit : Needs v-dependent time constants : from tables created under hoc NEURON { SUFFIX k3st USEION k READ ek WRITE ik RANGE g, gbar } UNITS { (mV) = (millivolt) } PARAMETER { gbar = 33 (millimho/cm2) d1 = -38 (mV) k1 = 0.151 (/mV) d2 = -25 (mV) k2 = 0.044 (/mV) } ASSIGNED { v (mV) ek (mV) g (millimho/cm2) ik (milliamp/cm2) kf1 (/ms) kb1 (/ms) kf2 (/ms) kb2 (/ms) } STATE { c1 c2 o } BREAKPOINT { SOLVE kin METHOD sparse g = gbar*o ik = g*(v - ek)*(1e-3) } INITIAL { SOLVE kin STEADYSTATE sparse } KINETIC kin { rates(v) ~ c1 <-> c2 (kf1, kb1) ~ c2 <-> o (kf2, kb2) CONSERVE c1 + c2 + o = 1 } FUNCTION_TABLE tau1(v(mV)) (ms) FUNCTION_TABLE tau2(v(mV)) (ms) PROCEDURE rates(v(millivolt)) { LOCAL K1, K2 K1 = exp(k2*(d2 - v) - k1*(d1 - v)) kf1 = K1/(tau1(v)*(1+K1)) kb1 = 1/(tau1(v)*(1+K1)) K2 = exp(-k2*(d2 - v)) kf2 = K2/(tau2(v)*(1+K2)) kb2 = 1/(tau2(v)*(1+K2)) } Listing 9.7. k3st.mod The NEURON block With one exception, the NEURON block of this model is essentially the same as for the delayed rectifier presented in Example 9.4: a voltage-gated current. The difference is that, even though this model contributes to the total K+ current ik, its own current is not available separately (i.e. there will be no ik_k3st at the hoc level) because ik is not declared as a RANGE variable. Variable declaration blocks The STATE block The STATEs in this mechanism are the fractions of channels that are in closed states 1 or 2 or in the open state. Since the total number of channels in all states is conserved, the sum of the STATEs must be unity, i.e. c1 + c2 + o = 1. This conservation rule means that the k3st mechanism really has only two independent state variables, a fact that underscores the difference between a STATE in NMODL and the general concept of a state variable. It also affects how NMODL sets up the equations that are to be solved, as we will see in the discussion of the KINETIC block below. Not all reactants or products need to be STATEs. If the reactant is an ASSIGNED or PARAMETER variable, then a differential equation is not generated for it, and it is treated as constant for the purposes of calculating the declared STATEs. Statements such as PARAMETER {kbath (mM)} STATE {ko (mM)} KINETIC scheme3 { ~ ko <-> kbath (r, r) } are translated to the single ODE equivalent ko' = r*(kbath - ko) i.e. ko tends exponentially to the steady state value of kbath. Equation definition blocks The BREAKPOINT block The recommended idiom for integrating a kinetic scheme is BREAKPOINT { SOLVE scheme METHOD sparse . . . } which integrates the STATEs in the scheme one dt step per call to fadvance() in NEURON. The sparse method is generally faster than computing the full Jacobian matrix, though both use Newton iterations to advance the STATEs with a fully implicit method (first-order correct). Additionally, the sparse method separates the Jacobian evaluation from the calculation of the STATE derivatives, thus allowing variable time step methods, such as CVODE, to efficiently compute only what is needed to advance the STATEs. Non-implicit methods, such as Runge-Kutta or Euler, should be avoided since kinetic schemes commonly have very wide ranging rate constants that make these methods numerically unstable with reasonable dt steps. In fact, it is not unusual to specify equilibrium reactions such as ~ A <-> B (1e6*sqrt(K), 1e6/sqrt(K)) which can only be solved by implicit methods. The INITIAL block Initialization of a kinetic scheme to its steady state is accomplished with INITIAL { SOLVE scheme STEADYSTATE sparse } Appropriate CONSERVE statements should be part of the scheme (see the following discussion of the KINETIC block) so that the equivalent system of ODEs is linearly independent. It should be kept in mind that source fluxes (constant for infinite time) have a strong effect on the steady state. Finally, it is crucial to test the scheme in NEURON under conditions in which the correct behavior is known. The KINETIC block The voltage-dependent rate constants are computed in the separate procedure rates(). That procedure computes the equilibrium constants K1 and K2 from the constants k1, d1, k2, and d2, whose empirically-determined default values are given in the PARAMETER block, and membrane potential v. The time constants tau1 and tau2, however, are found from tables created under hoc (see The FUNCTION_TABLEs below). The other noteworthy item in this block is the CONSERVE statement. As mentioned above in General comments about kinetic schemes, the basic idea is to systematically account for conservation of material. If there is neither a source nor a sink reaction for a STATE, the differential equations are not linearly independent when steady states are calculated (dt approaches infinity). For example, in scheme1 above the steady state condition m' = mc' = 0 yields two identical equations. Steady states can be approximated by integrating for several steps from any initial condition with large dt, but roundoff error can be a problem if the Jacobian matrix is nearly singular. To help solve the equations while maintaining strict numerical conservation throughout the simulation (no accumulation of roundoff error), the user is allowed to explicitly specify conservation equations with the CONSERVE statement. The conservation law for scheme1 is specified in NMODL by CONSERVE m + mc = 1 The CONSERVE statement does not add to the information content of a kinetic scheme and should be considered only as a hint to the translator. The NMODL translator uses this algebraic equation to replace the ODE for the last STATE on the left side of the equal sign. If one of the STATE names is an array, the conservation equation will contain an implicit sum over the array. If the last STATE is an array, then the ODE for the last STATE array element will be replaced by the algebraic equation. The choice of which STATE ODE is replaced by the algebraic equation is implementation-dependent and does not affect the solution (to within roundoff error). If a CONSERVEd STATE is relative to a compartment size, then compartment size is implicitly taken into account for the STATEs on the left hand side of the CONSERVE equation (see Example 9.8 for discussion of the COMPARTMENT statement). The right hand side is merely an expression, in which any necessary compartment sizes must be included explicitly. Thus in a calcium pump model Eq. 9.3 the pump is conserved and one could write CONSERVE pump + pumpca = total_pump * pumparea The FUNCTION_TABLEs As noted above, the steady-state clamp data define the voltage dependence of K1 and K2, but a complete description of the K+ current requires analysis of the temporal properties of the clamp current to determine the rate factors at each of the command potentials. The result would be a list or table of membrane potentials with associated time constants. One way to handle these numeric values would be to fit them with a pair of approximating functions, but the tactic used in this example is to leave them in tabular form for NMODL's FUNCTION_TABLE to deal with. This is done by placing the numeric values in three Vectors, say v_vec, tau1_vec, and tau2_vec, where the first is the list of voltages and the other two are the corresponding time constants. These Vectors would be attached to the FUNCTION_TABLEs of this model with the hoc commands table_tau1_k3st(tau1_vec, v_vec) table_tau2_k3st(tau2_vec, v_vec) Then whenever tau1(x) is called in the NMODL file, or tau1_k3st(x) is called from hoc, the returned value is interpolated from the array. A useful feature of FUNCTION_TABLEs is that, prior to developing the Vector database, they can be attached to a scalar value, as in table_tau1_k3st(100) effectively becoming constant functions. Also FUNCTION_TABLEs can be declared with two arguments and attached to doubly-dimensioned hoc arrays. In this case the table is linearly interpolated in both dimensions. This is useful with rates that depend both on voltage- and calcium. Usage Inserting this mechanism into a section makes the STATEs c1_k3st, c2_k3st, and o_k3st available at the hoc level, as well as the conductances gbar_k3st and g_k3st. Example 9.8: calcium diffusion with buffering This mechanism illustrates how to use kinetic schemes to model intracellular Ca2+ diffusion and buffering. It differs from the prior example in several important aspects: Ca2+ is not conserved but instead enters as a consequence of the transmembrane Ca2+ current; diffusion involves the exchange of Ca2+ between compartments of unequal size; Ca2+ is buffered. Only free Ca2+ is assumed to be mobile, whereas bound Ca2+ and free buffer are stationary. Buffer concentration and rate constants are based on the bullfrog sympathetic ganglion cell model described by Yamada et al. [Yamada, 1998 #580]. For a thorough treatment of numeric solution of the diffusion equations the reader is referred to Oran and Boris [Oran, 1987 #93]. Modeling diffusion with kinetic schemes Diffusion is modeled as the exchange of Ca2+ between adjacent compartments. For radial diffusion, the compartments are a series of concentric shells around a cylindrical core, as shown in Fig. 9.9 for Nannuli = 4. The index of the outermost shell is 0 and the index of the core is Nannuli - 1. The outermost shell is half as thick as the others so that [Ca2+] will be second order correct with respect to space at the surface of the segment. Concentration is also second order correct midway through the thickness of the other shells and at the center of the core. These depths are indicated by "x" in Fig. 9.9. The radius of the cylindrical core equals the thickness of the outermost shell, and the intervening Nannuli - 2 shells each have thickness ?r = diam / 2 (Nannuli - 1), where diam is the diameter of the segment. Figure 9.9. Diagram of the concentric shells used to model radial diffusion. The × mark the radial distances at which concentration will be second order correct in space. Because segment diameter and the number of shells affect the dimensions of the shells, they also affect the time course of diffusion. The flux between adjacent shells is ?[Ca2+] DCa A / ?r, where ?[Ca2+] is the concentration difference between the shell centers, DCa is the diffusion coefficient for Ca2+, A is the area of the boundary between shells, and ?r is the distance between their centers. This suggests that diffusion can be described by the basic kinetic scheme FROM i = 0 TO Nannuli-2 { ~ ca[i] <-> ca[i+1] (f[i+1], f[i+1]) } where Nannuli is the number of shells, ca[i] is the concentration midway through the thickness of shell i (except for ca[0] which is the concentration at the outer surface of shell 0), and the rate constants f[i+1] equal DCa Ai+1 / ?r. For each adjacent pair of shells, both Ai+1 and ?r are directly proportional to segment diameter. Therefore the ratios Ai+1 / ?r depend only on shell index, i.e. once they have been computed for one segment, they can be used for all segments that have the same number of radial compartments regardless of segment diameter. As it stands, this kinetic scheme is dimensionally incorrect. Dimensional consistency requires that the product of STATEs and rates be in units of STATE per time (also see Compartment size in the section on Chemical reactions in Chapter 3). In the present example the STATEs ca[] are intensive variables (concentration, or mass/volume), so the product of f[] and ca[] must be in units of concentration per time. However, the rates have units of volume per time, so this product is in units of mass per time, i.e. a flux that signifies the rate at which Ca2+ is entering or leaving a compartment. This flux is the time derivative of an extensive variable. This disparity is corrected by specifying STATE volumes with the COMPARTMENT statement, as in COMPARTMENT volume {state1 state2 . . . } where the STATEs named in the braces have the same compartment volume given by the volume expression after the COMPARTMENT keyword. The volume merely multiplies the dSTATE/dt left hand side of the equivalent differential equations, converting it to an extensive quantity and making it consistent with flux terms in units of absolute quantity per time. The volume of each cylindrical shell depends on its index and the total number of shells, and is proportional to the square of segment diameter. Consequently the volumes can be computed once for a segment with unit diameter and then scaled by diam^2 for use in each segment that has the same Nannuli. The equations that describe the radial movement of Ca2+ are independent of segment length. Therefore it is convenient to express shell volumes and surface areas in units of µm2 (volume/length) and µm (area/length), respectively. : Calcium ion accumulation with radial and longitudinal diffusion NEURON { SUFFIX cadifus USEION ca READ cai, ica WRITE cai GLOBAL vrat, TotalBuffer : vrat must be GLOBAL--see INITIAL block : however TotalBuffer may be RANGE } DEFINE Nannuli 4 : must be >=2 (i.e. at least shell and core) UNITS { (molar) = (1/liter) (mM) = (millimolar) (um) = (micron) (mA) = (milliamp) FARADAY = (faraday) (10000 coulomb) PI = (pi) (1) } PARAMETER { DCa = 0.6 (um2/ms) k1buf = 100 (/mM-ms) : Yamada et al. 1989 k2buf = 0.1 (/ms) TotalBuffer = 0.003 (mM) } ASSIGNED { diam (um) ica (mA/cm2) cai (mM) vrat[Nannuli] (1) : dimensionless : numeric value of vrat[i] equals the volume : of annulus iof a 1um diameter cylinder : multiply by diam^2 to get volume per um length Kd (/mM) B0 (mM) } STATE { : ca[0] is equivalent to cai : ca[] are very small, so specify absolute tolerance ca[Nannuli] (mM) <1e-10> CaBuffer[Nannuli] (mM) Buffer[Nannuli] (mM) } BREAKPOINT { SOLVE state METHOD sparse } LOCAL factors_done INITIAL { if (factors_done == 0) { : flag becomes 1 in the first segment factors_done = 1 : all subsequent segments will have factors() : vrat = 0 unless vrat is GLOBAL } Kd = k1buf/k2buf B0 = TotalBuffer/(1 + Kd*cai) FROM i=0 TO Nannuli-1 { ca[i] = cai Buffer[i] = B0 CaBuffer[i] = TotalBuffer - B0 } } LOCAL frat[Nannuli] : scales the rate constants for model geometry PROCEDURE factors() { LOCAL r, dr2 r = 1/2 : starts at edge (half diam) dr2 = r/(Nannuli-1)/2 : full thickness of outermost annulus, : half thickness of all other annuli vrat[0] = 0 frat[0] = 2*r FROM i=0 TO Nannuli-2 { vrat[i] = vrat[i] + PI*(r-dr2/2)*2*dr2 : interior half r = r - dr2 frat[i+1] = 2*PI*r/(2*dr2) : outer radius of annulus : div by distance between centers r = r - dr2 vrat[i+1] = PI*(r+dr2/2)*2*dr2 : outer half of annulus } } LOCAL dsq, dsqvol : can't define local variable in KINETIC block : or use in COMPARTMENT statement KINETIC state { COMPARTMENT i, diam*diam*vrat[i] {ca CaBuffer Buffer} LONGITUDINAL_DIFFUSION i, DCa*diam*diam*vrat[i] {ca} ~ ca[0] << (-ica*PI*diam/(2*FARADAY)) : ica is Ca efflux FROM i=0 TO Nannuli-2 { ~ ca[i] <-> ca[i+1] (DCa*frat[i+1], DCa*frat[i+1]) } dsq = diam*diam FROM i=0 TO Nannuli-1 { dsqvol = dsq*vrat[i] ~ ca[i] + Buffer[i] <-> CaBuffer[i] (k1buf*dsqvol, k2buf*dsqvol) } cai = ca[0] } Listing 9.8. cadif.mod The NEURON block This model READs cai to initialize the buffer (see The INITIAL block), and it WRITEs cai because it computes [Ca2+] in the outermost shell during a simulation run. It also READs ica, which is the Ca2+ influx into the outermost shell. There are two GLOBALs. One is the total buffer concentration TotalBuffer, which is assumed to be uniform throughout the cell. The other is vrat, an array whose elements will be the numeric values of the (volume/length) of the shells for a segment with unit diameter. These values are computed by PROCEDURE factors() near the end of Listing 9.8. As noted above, a segment with diameter diam has shells with volume/length equal to diam^2 * vrat[i]. Because each instance of this mechanism has the same number of shells, the same vrat[i] can be used to find the shell volumes at each location in the model cell where the mechanism exists. The DEFINE statement sets the number of shells to 4. Many of the variables in this model are arrays, and NMODL arrays are not dynamic. Instead, their lengths must be specified when the NMODL code is translated to C. The UNITS block Faraday's constant is scaled here in order to avoid having to include this scale factor as a separate term in the statement in the KINETIC block where transmembrane current ica is reckoned as the efflux of Ca2+ from the outermost shell. Since each statement in a UNITS block must include an explicit assertion of the units that are involved, the statement that assigns the value 3.141 . . . to PI includes a (1) which signifies that this is a dimensionless constant. Variable declaration blocks The ASSIGNED block The variable vrat is declared to be an array with Nannuli elements. As with C, array indices run from 0 to Nannuli - 1. The variables Kd and B0 are the dissociation constant for the buffer and the initial value of free buffer, which are computed in the INITIAL block (see below). Both the total buffer and the initial concentration of Ca2+ are assumed to be uniform throughout all shells, so a scalar is used for B0. The STATE block In addition to diffusion, this mechanism involves Ca2+ buffering Eq. 9.4 This happens in each of the shells, so ca, Buffer and CaBuffer are all arrays. The declaration of ca[] uses the syntax state (units) to specify the absolute tolerance that will be employed by CVODE. The solver tries to use a step size for which the local error ?i of each statei satisfies at least one of these two inequalities: ?i < relative_tolerance · |statei| or ?i < absolute_tolerance The default values for these tolerances are 0 and 10-3, respectively, so only a STATE that is extremely small (such as intracellular [Ca2+]) needs to have its absolute tolerance specified. As an alternative to specifying a smaller absolute tolerance, ca[] could have been defined in terms of units such as micromolar or nanomolar, which would have increased the numeric value of these variables. This would necessitate different scale factors in many of the statements that involve ca[]. For example, the assignment for cai (which is required to be in mM) would be cai = (1e-6)*ca[0]. LOCAL variables declared outside of equation definition blocks A LOCAL variable that is declared outside of an equation definition block is equivalent to a static variable in C. That is, it is visible throughout the mechanism (but not at the hoc level), it retains its value, and it is shared between all instances of a given mechanism. The initial value of such a variable is 0. This particular mechanism employs four variables of this type: factors_done, frat[], dsq, and dsqvol. The meaning of each of these is discussed below. Equation definition blocks The INITIAL block Initialization of this mechanism is a two step process. The first step is to use PROCEDURE factors() (see below) to set up the geometry of the model by computing the scale factor arrays vrat[] and frat[] that are applied to the shell volumes and rate constants. This only has to be done once because the same scale factors are used for all segments that have the same number of shells, as noted above in Modeling diffusion with kinetic schemes. The variable factors_done is a flag that indicates whether vrat[] and frat[] have been computed. The NMODL keyword LOCAL means that the value of factors_done will be the same in all instances of this mechanism, but that it will not be visible at the hoc level. Therefore factors() will be executed only once, regardless of how many segments contain the cadifus mechanism. The second step is to initialize the mechanism's STATEs. This mechanism assumes that the total buffer concentration and the initial free calcium concentration are uniform in all shells, and that buffering has reached its steady-state. Therefore the initial concentration of free buffer is computed from the initial [Ca2+] and the buffer's dissociation constant. It should be noted that the value of cai will be set to cai0_ca_ion just prior to executing the code in the INITIAL block (see Ion concentrations and equilibrium potentials in Chapter 8). It may be instructive to compare this initialization strategy with the approach that was used for the voltage-gated current of Listing 9.7 (k3st.mod). That previous example initialized the STATE through numeric solution of a kinetic scheme, so its KINETIC block required a CONSERVE statement to ensure that the equivalent system of ODEs would be linearly independent. Here, however, the STATEs are initialized by explicit algebraic assignment, so no CONSERVE statement is necessary. PROCEDURE factors() The arrays vrat[] and frat[], which are used to scale the shell volumes and rate constants to ensure consistency of units, are computed here. Their values depend only on the number of shells, so they do not have to be recomputed if diam or DFree is changed. The elements of vrat[] are the volumes of a set of concentric cylindrical shells, whose total volume equals the volume of a cylinder with diameter and length of 1 µm. These values are computed in two stages by the FROM i=0 TO Nannuli-2 { } loop. The first stage finds the volume of the outer half and the second finds the volume of the inner half of the shell. The frat array is declared to be LOCAL because it applies to all segments that have the cadifus mechanism, but it is unlikely to be of interest to the user and therefore does not need to be visible at the hoc level. This contrasts with vrat, which is declared as GLOBAL within the NEURON block so that the user can see its values. The values frat[i+1] equal Ai+1 / ?r, where Ai+1 is the surface area between shells i and i+1 for 0 ? i < Nannuli, and ?r is the distance between shell centers (radius / (Nannuli - 1)). The KINETIC block The first statement in this block specifies the shell volumes for the STATEs ca, CaBuffer, and Buffer. As noted above in Modeling diffusion with kinetic schemes, these volumes equal the elements of vrat[] multiplied by the square of the segment diameter. Because this mechanism involves many compartments whose relative volumes are specified by the elements of an array, this example takes care of all compartments with a single statement of the form COMPARTMENT index, volume[index] { state1 state2 . . . } where the STATEs that are diffusing are listed inside the braces. Next in this block is a LONGITUDINAL_DIFFUSION statement, which specifies that this mechanism includes nonlocal diffusion, i.e. longitudinal diffusion along a section and into connecting sections. The syntax for scalar STATEs is LONGITUDINAL_DIFFUSION flux_expr { state1 state2 . . . } where flux_expr is the product of the diffusion constant and the cross-sectional area between adjacent compartments. Units of the flux_expr must be (micron4/ms), i.e. the diffusion constant has units of (micron2/ms) and the cross-sectional area has units of (micron2). For cylindrical shell compartments, the cross-sectional area is just the volume per unit length. If the states are arrays then all elements are assumed to diffuse between corresponding volumes in adjacent segments and the iteration variable must be specified as in LONGITUDINAL_DIFFUSION index, flux_expr(index) { state1 state2 . . . } A COMPARTMENT statement is also required for the diffusing STATEs and the units must be (micron2), i.e. (micron3/micron). The compactness of LONGITUDINAL_DIFFUSION specification contrasts nicely with the great deal of trouble imposed on the computational methods used to solve the equations. The standard fixed time step implicit method, historically the default method used by NEURON (see Chapter 4), can no longer find steady states with extremely large (e.g. 109 ms) steps since not every Jacobian element for both flux and current with respect to voltage and concentration is presently accurately computed. The CVODE method works well for these problems since it does not allow dt to grow beyond the point of numerical instability. Despite these occasional limitations on numerical efficiency, it is satisfying that, as methods evolve to handle these problems more robustly, the specification of the models does not change. The third statement in this block is equivalent to a differential equation that describes the contribution of transmembrane calcium current to Ca2+ in the outermost shell. The << signifies an explicit flux. Because of the COMPARTMENT statement, the left hand side of the differential equation is not d[Ca2+]0/dt but d(total Ca2+ in the outermost shell)/dt. This is consistent with the right hand side of the equation, which is in units of mass per time. Next is the kinetic scheme for radial diffusion. The rate constants in this scheme equal the product of DCa and the factor frat[] for reasons that were explained earlier in Modeling diffusion with kinetic schemes. It may not be immediately clear why the rate constants in the kinetic scheme for Ca2+ buffering are scaled by the compartment volume dsqvol; however, the reason will become obvious when one recalls that the COMPARTMENT statement at the beginning of the KINETIC block has converted the units of the dSTATE/dt on the left hand side of the equivalent differential equations from concentration per time to mass per time. If the reaction rate constants were left unchanged, the right hand side of the differential equations for buffering would have units of concentration per time, which is inconsistent. Multiplying the rate constants by compartment volume removes this inconsistency by changing the units of the right hand side to mass per time. The last statement in the KINETIC block updates the value of cai from ca[0]. This is necessary because intracellular [Ca2+] is known elsewhere in NEURON as cai, e.g. to other mechanisms and to NEURON's internal routine that computes ECa. When developing a new mechanism or making substantive changes to an existing mechanism, it is generally advisable to check for consistency of units with modlunit. Given the dimensional complexity of this model, such testing is absolutely indispensable. Usage If this mechanism is inserted in a section, the concentrations of Ca2+ and the free and bound buffer in all compartments will be available through hoc as ca_cadifus[], Buffer_cadifus[], and CaBuffer_cadifus[]. These STATEs will also be available for plotting and analysis through the GUI. The PARAMETERs DCa, k1buf, k2buf, and TotalBuffer will also be available for inspection and modification through both the graphical interface and hoc statements (with the _cadifus suffix). All PARAMETERs are GLOBALs by default, i.e. they will have the same values in each location where the cadifus mechanism has been inserted. Therefore in a sense it is gratuitous to declare in the NEURON block that TotalBuffer is GLOBAL. However, this declaration does serve the purpose of underscoring the nature of this important variable which is likely to be changed by the user. In some cases it might be useful for one or more of the PARAMETERs to be RANGE variables. For example, TotalBuffer and even DCa or the buffer rate constants might not be uniform throughout the cell. To make TotalBuffer and DCa RANGE variables only requires replacing the line GLOBAL vrat, TotalBuffer in the NEURON block with GLOBAL vrat RANGE TotalBuffer, DCa The GLOBAL volume factors vrat[] are available through hoc for inspection, but it is inadvisable to change their values because they would likely be inconsistent with the frat[] values and thereby cause errors in the simulation. All occurrences of this mechanism will have the same number of shells, regardless of the physical diameter of the segments in which the mechanism has been inserted. With Nannuli = 4, the thickness of the outermost shell will be ? 1 µm in segments with diam ? 6 µm. If this spatial resolution is inadequate, or if the model has segments with larger diameters, then Nannuli may have to be increased. NMODL does not have dynamic arrays, so in order to change the number of shells one must recompile the mechanism after assigning a new value to Nannuli by editing the NMODL source code. Example 9.9: a calcium pump This mechanism involves a calcium pump based on the reaction scheme outlined in the description of the KINETIC block of Example 9.7: kinetic scheme for a voltage-gated current. It is a direct extension of the model of calcium diffusion with buffering in Example 9.8: calcium diffusion with buffering, the principal difference being that a calcium pump is present in the cell membrane. The following discussion focuses on the requisite changes in Listing 9.8, and the operation and use of this new mechanism. For all other details the reader should refer to Example 9.8. The NEURON block Changes in the NEURON block are marked in bold. The first nontrivial difference from the prior example is that this mechanism READs the value of cao, which is used in the pump reaction scheme. NEURON { SUFFIX cdp USEION ca READ cao, cai, ica WRITE cai, ica RANGE ica_pmp GLOBAL vrat, TotalBuffer, TotalPump } The mechanism WRITEs a pump current that is attributed to ica so that its transmembrane Ca2+ flux will be factored into NEURON's calculations of [Ca2+]i. This current, which is a RANGE variable known as ica_pmp_cdp to the hoc interpreter, constitutes a net movement of positive charge across the cell membrane, and it follows the usual sign convention (outward current is "positive"). The pump current has a direct effect on membrane potential, which, because of the rapid activation of the pump, is manifest by a distinct delay of the spike peak and a slight increase of the postspike hyperpolarization. This mechanism could be made electrically "silent" by having it WRITE an equal but opposite NONSPECIFIC current or perhaps a current that involves some other ionic species, e.g. Na+, K+, or Cl-. The variable TotalPump is the total density of pump sites on the cell membrane, whether free or occupied by Ca2+. Making it GLOBAL means that it is user adjustable, and that the pump is assumed to have uniform density wherever the mechanism has been inserted. If local variation is required, this should be a RANGE variable. The UNITS block This mechanism includes the statement (mol) = (1) because the density of pump sites will be specified in units of (mol/cm2). The term mole cannot be used here because it is already defined in NEURON's units database as 6.022169 · 1023. Variable declaration blocks The PARAMETER block Five new statements have been added because this mechanism uses the rate constants of the pump reactions and the density of pump sites on the cell membrane. k1 = 1 (/mM-ms) k2 = 0.005 (/ms) k3 = 1 (/ms) k4 = 0.005 (/mM-ms) : to eliminate pump, set TotalPump to 0 in hoc TotalPump = 1e-11 (mol/cm2) These particular rate constant values were chosen to satisfy two criteria: the pump influx and efflux should be equal at [Ca2+] = 50 nM, and the rate of transport should be slow enough to allow a slight delay in accelerated transport following an action potential that included a voltage-gated Ca2+ current. The density TotalPump is sufficient for the pump to have a marked damping effect on [Ca2+]i transients; lower values reduce the ability of the pump to regulate [Ca2+]i. The ASSIGNED block These three additions have been made. cao (mM) ica_pmp (mA/cm2) parea (um) This mechanism treats [Ca2+]o as a constant. The pump current and the surface area over which the pump is distributed are also clearly necessary. The CONSTANT block Consistency of units requires explicit mention of an extracellular volume in the kinetic scheme for the pump. CONSTANT { volo = 1e10 (um2) } The value used here is equivalent to 1 liter of extracellular space per micron length of the cell, but the actual value is irrelevant to this mechanism because cao is treated as a constant. Since the value of volo is not important for this mechanism, there is no need for it to be accessible through hoc commands or the GUI, so it is not a PARAMETER. On the other hand, there is a sense in which it is an integral part of the pump mechanism, so it would not be appropriate to make volo be a LOCAL variable since LOCALs are intended for temporary storage of "throwaway" values. Finally, the value of volo would never be changed in the course of a simulation. Therefore volo is declared in a CONSTANT block. The STATE block The densities of pump sites that are free or have bound Ca2+, respectively, are represented by the two new STATEs pump (mol/cm2) pumpca (mol/cm2) Equation definition blocks The BREAKPOINT block This block has one additional statement BREAKPOINT { SOLVE state METHOD sparse ica = ica_pmp } The assignment ica = ica_pmp is needed to ensure that the pump current is reckoned in NEURON's calculation of [Ca2+]i. The INITIAL block The statement parea = PI*diam must be included to specify the area per unit length over which the pump is distributed. If it is correct to assume that [Ca2+]i has been equal to cai0_ca_ion (default = 50 nM) for a long time, the initial levels of pump and pumpca can be set by using the steady-state formula pump = TotalPump/(1 + (cai*k1/k2)) pumpca = TotalPump - pump An alternative initialization strategy is to place ica = 0 SOLVE state STEADYSTATE sparse at the end of the INITIAL block, where the ica = 0 statement is needed because the kinetic scheme interprets transmembrane Ca2+ currents as a source of Ca2+ flux. This idiom makes NEURON compute the initial values of STATEs, which can be particularly convenient for mechanisms whose steady state solutions are difficult or impossible to express in analytical form. This would require adding a CONSERVE statement to the KINETIC block to insure that the equations that describe the free and bound buffer are independent (see also The INITIAL block in Example 9.7: kinetic scheme for a voltage-gated current). Both of these initializations explicitly assume that the net Ca2+ current generated by other sources equals 0, so the net pump current following initialization is also 0. If this assumption is incorrect, as is almost certainly the case if one or more voltage-gated Ca2+ currents are included in the model, then [Ca2+]i will start to change immediately when a simulation is started. Most often this is not the desired outcome. The proper initialization of a model that contains mechanisms with complex interactions may involve performing an "initialization run" and using SaveState objects (see Examples of custom initializations in Chapter 8). The KINETIC block Changes in this block are marked in bold. The new COMPARTMENT statements and the scale factor (1e10) are required for dimensional consistency in the pump scheme. KINETIC state { COMPARTMENT i, diam*diam*vrat[i] {ca CaBuffer Buffer} COMPARTMENT (1e10)*parea {pump pumpca} COMPARTMENT volo {cao} LONGITUDINAL_DIFFUSION i, DCa*diam*diam*vrat[i] {ca} :pump ~ ca[0] + pump <-> pumpca (k1*parea*(1e10), k2*parea*(1e10)) ~ pumpca <-> pump + cao (k3*parea*(1e10), k4*parea*(1e10)) CONSERVE pump + pumpca = TotalPump * parea * (1e10) ica_pmp = 2*FARADAY*(f_flux - b_flux)/parea : all currents except pump ~ ca[0] << (-(ica - ica_pmp)*PI*diam/(2*FARADAY)) FROM i=0 TO Nannuli-2 { ~ ca[i] <-> ca[i+1] (DCa*frat[i+1], DCa*frat[i+1]) } dsq = diam*diam FROM i=0 TO Nannuli-1 { dsqvol = dsq*vrat[i] ~ ca[i] + Buffer[i] <-> CaBuffer[i] (k1buf*dsqvol, k2buf*dsqvol) } cai = ca[0] } The pump reaction statements implement the scheme outlined in the description of the KINETIC block of Example 9.7: kinetic scheme for a voltage-gated current. Also as described in that section, the CONSERVE statement ensures strict numerical conservation, which is helpful for convergence and accuracy. In the steady state, the net forward flux in the first and second reactions must be equal. Even during physiologically-relevant transients, these fluxes track each other effectively instantaneously. Therefore the transmembrane Ca2+ flux generated by the pump is taken to be the net forward flux in the second reaction. This mechanism WRITEs ica in order to affect [Ca2+]i. The total transmembrane Ca2+ flux is the sum of the pump flux and the flux from all other sources. Thus to make sure that ica_pmp is not counted twice, it is subtracted from total Ca2+ current ica in the expression that relates Ca2+ current to Ca2+ flux. Usage The STATEs and PARAMETERs that are available through hoc and the GUI are directly analogous to those of the cadifus mechanism, but they will have the suffix _cdp rather than _cadifus. The additional pump variables pump_cdp, pumpca_cdp, ica_pmp_cdp, and TotalPump_cdp will also be available and are subject to similar concerns and constraints as their counterparts in the diffusion reactions (see Usage in Example 9.8: calcium diffusion with buffering). Models with discontinuities The incorporation of variable time step integration methods in NEURON made it necessary to provide a way to ensure proper handling of abrupt changes in PARAMETERs, ASSIGNED variables, and STATEs. This was accomplished by adding two new statements to NMODL: at_time() and state_discontinuity(), which we discuss here. The later addition of an event delivery system to NEURON has made at_time() a deprecated function, and has almost obviated the need for state_discontinuity(). However, these statements have been used in many mechanisms of recent vintage, e.g. models of pulse generators and synaptic transmission, so it is important to understand why they were used and what they do. Having said that, we strongly recommend against using these functions in any new model development. Discontinuities in PARAMETERs Before CVODE was added to NEURON, sudden changes in PARAMETERs and ASSIGNED variables, such as the sudden change in current injection during a current pulse, had been implicitly assumed to take place on a time step boundary. This is inadequate with variable time step methods because it is unlikely that a time step boundary will correspond to the onset or offset of the pulse. Worse, the time step may be longer than the pulse itself, which may thus be entirely ignored. The at_time() function was added so that a model description could explicitly notify NEURON of the times at which any discontinuities occur. This function has no effect on fixed time step integration. However, during variable time step integration, the statement at_time(event_time) guarantees that the integrator reduces the step size so that it completes at t = event_time - ?, where ? ~ 10-9 ms. The next step resets the integrator to first order, thereby discarding any previous solution history, and immediately returns after computing all the at t = event_time + ?. Note that at_time() returns a value of 1 ("true") only during the "infinitesimal" step that ends at t = event_time + ?; otherwise it returns 0. This is how the built-in current clamp model IClamp notifies NEURON of the time of onset of the pulse and its offset (see the BREAKPOINT block of Example 9.3: an intracellular stimulating electrode). More recently, however, the preferred way to implement abrupt changes in PARAMETERs and ASSIGNED variables is to take advantage of NEURON's event delivery system (specifically, self-events) because of improved computational efficiency and greater conceptual clarity (see Chapter 10). During a variable time step simulation, a missing at_time() call may cause one of two symptoms. If a PARAMETER changes but returns to its original value within the same interval, the pulse may be entirely missed. More often, a single discontinuity will take place within a time step interval, causing the integrator to start what seems like a binary search for the location of the discontinuity in order to satisfy the error tolerance on the step; of course, this is very inefficient. Time dependent PARAMETER changes Time dependent PARAMETER changes can be specified explicitly in a model description, as in BREAKPOINT { i = imax*sin(w*t) } which works with both fixed and variable time step integration. Time-dependent changes can also be specified at the hoc interpreter level, but care is needed to ensure they are properly computed in the context of variable time steps. For instance, it might seem convenient to change PARAMETERs prior to fadvance() calls, e.g. proc advance() { IClamp[0].amp = imax*sin(w*t) fadvance() } This does work with fixed dt but is discouraged because it produces unacceptable results with variable dt methods. An alternative that works well with fixed and variable time step integration is to use the Vector class's play() method with linear interpolation, which became available in NEURON 5.4. This is invoked with vec.play(&rangevar, tvec, 1) in which vec and tvec are a pair of Vectors that define a piecewise linear function of time y = f(t), i.e. tvec contains a monotonically increasing sequence of times, and vec holds the corresponding y values. The rangevar is the variable that is to be driven by f(). In the future, Vector.play will be extended to cubic spline interpolation and will allow "continuous" play of a smooth function defined by a Vector. Another strategy would be to use the NetCon class's event() method, which exploits NEURON's event delivery system (see Chapter 10). The handler() procedure in netcon.event(te, "handler()") can contain statements that change anything discontinuously, as long as the last statement in handler is cvode.re_init(). Discontinuities in STATEs Some kinds of synaptic models involve a discontinuity in one or more STATE variables. For example, a synapse whose conductance follows the time course of an alpha function (for more detail about the alpha function itself see Rall [Rall, 1977 #108] and Jack et al. [Jack, 1983 #90]) can be implemented as a kinetic scheme in the two state model KINETIC state { ~ a <-> g (k, 0) ~ g -> (k) } where a discrete synaptic event results in an abrupt increase of STATE a. This formulation has the attractive property that it can handle multiple streams of events with different weights, so that g will be the sum of the individual alpha functions with their appropriate onsets. Abrupt changes in STATEs require particularly careful treatment because of the special nature of states in variable time step ODE solvers. Before the advent of an event delivery system in NEURON, this required not only an at_time() call to notify NEURON about the time of the discontinuity, but also a state_discontinuity() statement to specify how the affected STATE would change. Furthermore, state_discontinuity() could only be used in an if (at_time()){} block. Thus if onset is the time of the synaptic event and gmax is the desired maximum conductance change, the BREAKPOINT block would look like this BREAKPOINT { if (at_time(onset)) { : scale factor exp(1) = 2.718... ensures : that peak conductance will be gmax state_discontinuity(a, a + gmax*exp(1)) } SOLVE state METHOD sparse i = g*(v - e) } The first argument to state_discontinuity() is interpreted as a reference to the STATE, and the second argument is an expression for its new value. The first argument will be assigned the value of its second argument just once for any time step. This is important because, for several integration methods, BREAKPOINT assignment statements are often executed twice to calculate the di/dv terms of the Jacobian matrix. This synaptic model works well with deterministic stimulus trains, but it is difficult for the user to supply the administrative hoc code for managing the onset and gmax variables to take advantage of the promise of "multiple streams of input events with different weights." The most important problem is how to save events that have significant delay between their generation and their handling at time onset. As is, an event can be passed to this model by assigning values to onset and gmax only after the previous onset event has been handled. These complexities have been eliminated by the event delivery system. Instead of handling the state discontinuity in the BREAKPOINT block, the synaptic model should now be written in the form BREAKPOINT { SOLVE state METHOD sparse i = g*(v - e) } NET_RECEIVE(weight (microsiemens)) { a = a + weight*exp(1) } in which event distribution is handled internally from a specification of network connectivity (see next section). Note that there is no need to use either at_time or state_discontinuity. Also, the BREAKPOINT block should not have any if statements. All discontinuities should be handled in a NET_RECEIVE block. For further details of how to deal with streams of synaptic events with arbitrary delays and weights, see Chapter 10.