Chapter 4. Development of simulation tools-- culminating in NEURON.

You may jump to these later topics:

Inhomogeneous Morphologies

Other Problems

Mike Hines Develops NEURON


Having benefitted from long exposure to the HH equations and their power, I arrived at Duke in 1961 fired by the possibility of using them to gain insight into the functioning of whole nerve cells. Of course neither the software nor laboratory hardware required for such computer simulations was available at that time. Nevertheless, I began a 3+ decade effort of accumulate information on how bits and pieces of nerve cells function; this necessarily was in a step-wise fashion as both software and computer power developed. The key to the success of this process was in the choosing nad being chosen by students and postdocs who were capable mathematicians and who were excellent computer programmers. I learned almost I know from their instruction. Over these years many were generous in their efforts to teach me; they include Ron Joyner (as an undergraduate student, graduate student and later as a postdoc), Ed Cox, Fidel Ramon, Monte Westerfield, Nelson Arispe, Stewart Jaslove, and Norman Stockbridge.

Because my goal was always to gain more insight into the functioning of nerves, these developments were closely linked to experiments in our own laboratory. Early digital laboratory computers, Digital Equipment Corp. (Dec) A Linc8 and a PDP8, were used to run computer controlled experiments and analyze data (Joyner & Moore, 1973). Students and Postdocs also programmed and ran simulations to compliment experiments such as the spread of potential and current in squid giant axons. Later we turned to branches and changes in diameter as our programs and computing power expanded. We also used simulations to evaluate the quality of our instrumentation and understand their limits. Only later, with Mike Hines, did we begin to develop much more general purpose simulations.

My first simulations of HH equations were done on an Electronic Associates TR48 analog computer. While these simulations were very fast (0.1sec or less) and could be displayed on an oscilloscope they were limited to a uniform patch of membrane. These simulations were helpful in evaluating expressions for ionic channel kinetics as alternatives to the H& H equations. One of my first investigations was to see if a voltage sweep (as used in my study the tri-phasic K current-voltage characteristics of a squid axon in high potassium) could be used to quidkly display the sodium current-voltage characteristics. At first glance I was excited find that the current records obtained with this method had the characteristic shape of the plot of the peak sodium currents as a function of voltage obtained from a family of voltage steps. Further runs made it this similarity was only apparent and misleading. The outward arm of the "sodium" current at voltages above the sodium equilibrium potential was in fact carried by potassium ions. Thus I was surprised later to find that the so-called "ramp clamp" was being proposed as a fast way to obtain the current-voltage relation for the sodium channel.

IMPROVING METHODS AND INCREASING COMPLEXITY. Fidel Ramon and I (Moore & Ramon, 1974) were interested in optimizing computational speed without compromising accuracy. We used the early interpretive language FOCAL to evaluate a variety of methods of numerical integration of the Hodgkin and Huxley equations for a membrane action potential. We developed novel hybrid "m- corrected" Euler integration scheme which was 2-3 times faster than a Runga-Kutta integration of the same accuracy. This method is presently available in the simulation program NODUS (Eric De Schutter).

This work was followed by a group of four papers in the January 1975 issue of the Biophysical Journal on simulations of realistic voltage-clamp circuits connected to axons in a variety of experimental situations. In the first (Moore et. al., 1975a) we developed and tested our methods. The second paper (Moore et. al., 1975b) examined the quality of the double sucrose-gap voltage-clamp in use in our lab on lobster and squid giant axons. We found that the "artificial node" length which we used (equal to or less than one-half the axon diameter) provided a very good voltage clamp. Only when this length was increased dramatically and longitudinal voltage gradients were large did we see infamous "notches" in the voltage clamp currents - notches which indicated space clamp failure.

The third paper (Joyner et al. 1975) considered a two microelectrode voltage clamp of the postsynaptic region of the squid giant synapse. It included the addition of a drop of oil injected to isolate the postsynaptic region from the rest of the giant axon. Dramatic improvement in the quality of the clamp was shown to occur when the intracellular region was infused with tetraethylammonium chloride made the membrane potential more uniform. The last paper (Ramon et al 1975) dealt with the then very serious problem of the quality of voltage clamps of bundles of nerve or muscle in sucrose-gaps which were frequently being used at that time. Our simulations of these multicellular preparations showed that, in spite of an apparently good clamp as judged by the fact that the recorded voltage matched the command potential and no infamous "notches" were seen, propagated impulses were observed in much of the preparation. Nevertheless, we also found that the measured currents in the "artificial node" were not too badly distorted because the longitudinal currents associated with the propagating spikes were relatively small compared to the radial membrane currents.

These simulations were carried out using FORTRAN on a DEC PDP8 computer. In spite of its very limited memory (only 28Kbytes of 12 bit memory) creative and masterful programming on the part of my students allowed us to find solutions of the partial differential equations describing these time-voltage profiles in axons under a variety of conditions .


Over the next few years we examined impulse propagation under a variety of inhomogeneous conditions:

MYELINATED FIBERS: In collaboration with Steve Waxman's lab, Ron Joyner and I developed a simulation program (Brill et al, 1977). Later we (Moore et al,1978) also used it to examine the sensitivity of the conduction velocity to a wide variety of parameters and found that it was quite insensitive to the nodal area or to the exact description of its excitable processes. The normal internodal length of 2 mm seems optimal because conduction velocity is almost independent of a 50% change. On the other hand the velocity of impulse propagation varies almost linearly with diameter of the fiber and is relatively sensitive to the axoplasmic resistivity and the myelin capacitance. This program was useful also in gaining detailed insight into the problems of propagation failure in multiple sclerosis (Waxman & Wood, 1984).

CHANGES IN AXON DIAMETER (non- myelinated). We found that impulses propagating in squid giant axons underwent striking alterations of shape and amplitude as they progress from a normal region into one containing a low resistance axial wire. We drew attention to the fact that these shapes were very similar to those seen in antidromic impulse propagation into the soma of cat spinal motoneurons by using the title "Squid giant axons: a model for the neuron soma?" (Ramon et al 1976).

SOURCE OF INFLECTIONS IN RISING PHASE OF ACTION POTENTIALS. We refined our question further to compare changes in impulse shapes in regions of alterations in morphological and membrane channel properties. Monte Westerfield used experiments with squid giant axons employing:

  • (a) isolated compartments in which agents were added to change membrane characteristics and
  • (b) axial wires to effect an electrical change in diameter. These experiments and parallel simulations (Moore & Westerfield, 1983) showed that reduced excitability reduced the rate of rise and peak amplitude of the impulses but they were clearly inflection-free. On the other hand, impulses propagating into a region of diameter increase readily showed marked inflection points and, furthermore, a second impulse following too closely could be blocked at such a transition. In addition these studies showed a rapid biphasic change in threshold for impulse initiation, with a reduced threshold in the smaller diameter region adjacent to the larger.

    SIMPLE NEURON MODEL. These results encouraged us to pursue the more general question of the roles of morphology and channel densities in the initiation of impulses in neurons.

    The proposal in 1957 by two groups of investigators (Coombs , Curtis and Eccles; Frank, Fourtes and Becker) that orthodromic impulses in motoneurons are initiated in the initial segment (i.s.) was generally accepted. The soma membrane potential follows the i.s., usually with an inflection point in the rising phase at a depolarization of about 30 mV (Coombs el al, 1957). Because of this both groups concluded that the threshold in the soma was much higher than in the axon. Fred Dodge (Dodge, 1979; Dodge & Cooley,1973) developed a computer simulation of a motoneuron based on this conclusion and altered the threshold in these regions by using different sets for the voltage sensitive rate constants (in the HH equations) for the soma and axon. These assumptions were troubling because it required markedly different voltage sensitivity of Na channels in adjacent regions of a cell. This necessitated both the manufacture of differing protein structures and a means tolocate and keep them separated.

    In our experiments on squid axons, we had consistently seen inflections at 10's of millivolts in the rising phases where propagating spikes approached a transition to a larger diameter axon. Here it was known that the channel densities were uniform throughout. Because of this, we raised the possibility of whether, with identical membrane excitability in the soma and axon, an impulse could be initiated in the initial segment of the axon. Another possibility we considered was that the channel densities were higher in the initial segment.

    We undertook simulations whose purpose was to gain insight into the individual contributions of several parameters. Because few of the parameters values necessary for modeling were known, we made what we thought were reasonable assignments and then examined the sensitivity of the model neuron's response to variations in these parameter values. We carried out extensive simulations with a simple model of a neuron with the dendrites cast into an electrically equivalent Rall cylinder. We investigated both orthodromic and antidromic initiation of impulses, examining the effect of HH ion channel densities and locations, temperature changes, and the abruptness of the transition from axon to soma.

    We found that the most important factor in determining the shape, size and site of the initiation of the spike was the magnitude of the electrical (capacitive) load relative to the conductance of the source driving it. In the case of antidromic invasion, the action potential is affected by

  • (a) the abruptness of the junction between axon and soma,
  • (b) the ionic channel rate constants (temperature dependent),
  • (c) the length of the active soma, and
  • (d) the ionic channel density. Our simulations showed that, for passive dendrites and uniform density of HH Na channels in the soma, axon hillock and axon (non-myelinated), action potentials were generated in the axon hillock region. This initiation of spikes in the region of reduced threshold noted above (Moore et al. 1983) in spite of the fact that here the channels and their densities were the same throughout held true for a wide range of parameter values. An increased channel density in the axon and hillock area simply reduced the threshold there. Since then anatomical and immunological evidence has increasingly indicated that Na channel density was indeed higher in the hillock and initial segment area than in the soma.

    Recent models for more realistic motoneurons use NEURON to incorporate myelination of the axon, starting with short internodal lengths gradually increasing to the normal length. We find that myelination of the axon also plays an important part in the location of the initiation of the spike. In fact, the spike is initiated several nodes beyond the initial segment and back propagate through the axon into the soma. At synaptic input strengths just below that at which an action potential appears in the soma, it is also possible to find an area of parameter space where spikes initiated in the myelinated region - but are unable to invade the axon-soma complex.

    TEMPERATURE CHANGES. It is well known that cooling and warming axons have dramatic effects on the duration of the impulse and its conduction velocity. We expected to see significant changes in the shape of action potentials in the region of transition in diameter but nevertheless we were surprised at the sharpness of the sensitivity at axon branch points to temperature changes (Joyner et al, 1980). We found not only a significant increase in the peak sodium current but also a marked change in the shape of its time course; The shape of the potassium m current was markedly altered but its peak amplitude was about the same. The time integral of both currents also changes significantly in this transition region.


    NONLINEAR CABLES, Internal current injection. Steady-state distribution of current and voltage as a function of distance. Arispe & I (1979) extended the Hodgkin & Rushton (1946) treatment of axons in the linear region to cover large voltage excursions where the ion sensitive channels introduced dramatic nonlinearities. Both experiments on squid axons and simulations with the HH equations showed significant nonlinearities when the depolarization exceeded 1-2 mV. We were also able to validate the accuracy and power of Cole's equation for extracting the nonlinear membrane characteristics simply from the measurement of the input resistance.

    NONLINEAR CABLES, External current electrodes. Nelson Arispe & I also carried out experiments and simulations to determine the nonlinear steady-state current and voltage distributions for an axon in a close-fitting channel (Moore & Arispe, 1979). Our experimental data could be well fitted by the simulations when the leakage conductance was reduced to one-fifth the standard HH value. This both validated another application of the HH equations and further indicated that our axons were probably in better physiological condition than theirs. This is a particularly interesting and important problem for prosthetic devices for stimulation of nerves with external electrodes. In this complex situation, the entry and exit of current from the axon has sharp discontinuities and varies nonlinearly with the membrane depolarization. Because of the complexity of the problem, we used a hybrid method of simulation in which an analog computer solved the cable equations and a digital computer calculated the membrane current densities.

    EPHAPTIC TRANSMISSION. Fidel Ramon proposed that it might be useful in the study of transmission between two axons in close apposition to each other and where the free space for extracellular solutions was very limited. We found that it was relatively easy to produce such conditions and transmission by putting the two axons in close proximity and submerging this overlap region in oil (Ramon & Moore, 1978). In fact, noting that the squid axon could be dissected to such a length that the time for an impulse to propagate from one end to the other was relatively long compared to the duration of the action potential, Fidel was able to set up the situation were the impulse continued to propagate in a circular loop through a long axon doubled back on itself!

    PROPAGATION IN A TWO DIMENSIONAL SHEET. Joyner and Ramon were also interested in impulse propagation in cardiac tissue. They extended these methods to solve for spread of action potentials in a two dimensional inhomogenous sheet representing such muscle (Joyner, Ramon and Moore, 1975).

    KINETIC MODEL for the SODIUM CHANNEL interacting with CALCIUM. The striking and strong effect of the calcium ion concentration on the amplitude and time courses of both sodium and potassium currents in a voltage clamp was demonstrated by Frankenhauser and Hodgkin (1957). Although the changes in the different kinetic and amplitudes parameters of these conductances varied significantly, they chose to simplify the situation and summarize their observations by saying that a five-fold reduction in the calcium concentration was equivalent to a 10-15mV depolarization. Because simulations of the HH equations with a depolarization of this magnitude gave poor representations of the effects of reductions in [Ca], Ed Cox and I decided to try to develop a model for the sodium conductance in the form of kinetic reactions between states and included calcium as a specific reactant. We (Moore & Cox, 1976) investigated a variety of kinetic sequences and were able to find one which gave (1) a good fit to the standard HH conductances, (2) gave an excellent fit to the data for reduced [Ca], and (3) when substituted for the sodium conductance in the HH equations, produced a more realistic falling phase of the action potential - missing the aberrant or "gratuitous" hump in the decline of an HH spike.

    KINETIC MODEL for the POTASSIUM CHANNEL. Earlier, Kacy Cole and I had tested the validity of the Hodgkin and Huxley expression the n4 for the squid axon's K channel conductance (in response to depolarization in a voltage clamp). We preceeded these depolarizing voltage steps with a wide range of levels of conditioning voltages and found that the K currents could be superimposed by translation along the time axis. These records validated the HH single-state variable description. In addition, we found a most striking temporal delay in the onset of caused by a preceeding strong hyperpolarization (Cole & Moore, 1960). Although the onset of the K conductance upon depolarization from rest could be described better by n6 than the n4 used by Hodgkin and Huxley, we found that, in order to express it in the HH format, we had to raise n to the 25th power to match the data! Subsequently this delay has been dubbed the "Cole- Moore effect". Steve Young and I (Moore and Young, 1981a, 1981b) reexamined and extended these measurements in squid and crayfish giant axons by using tetrodotoxin to block the Na current, leaving uncontaminated K currents over a much wider voltage range. We showed that a 4-state kinetic model was sufficient to describe the crayfish axon K channel over the range of conditioning and test voltages used.

    NEW NUMERICAL METHOD. Ron Joyner lead the group in a significant leap by extending the Crank-Nicholson method for implicit integration of cable equations (describing axons with excitable membranes) to simulate propagation in complex morphologies, including branching (Joyner et al, 1978). This method was general and powerful enough to handle radius changes as a function of distance along an axis, any number of branched processes (any or all may be nonuniform in diameter), and with no restriction on the branching pattern. With this method, a systematic investigation of impulse propagation through a variety of transitions in diameter were simulated.

    Back to top of Chapter
    On to NEURON's Development by Mike Hines.