Chapter 3. Professional Style for Simulations

3.1 Fitting Data
3.2 Iterative Manual Searches
3.3 Sensitivity Testing & "Sensible" fitting
3.4 Numerical Accuracy

3.1 Fitting Data

At the heart of simulations is the attempt and goal to "fit" the available data with curves calculated from a model. However there are traps which modeller should be aware of. This chapter illustrates some of the pitfalls of fitting data which the informed modeler avoids.

A truly professional approach to simulations involves

A Nobel prize-winning early example of elegant fitting of experimental data

In their classic 1952 paper, Hodgkin & Huxley (HH) developed their equations for currents in the squid axon membrane, providing a stellar example of how to fit experimental data with mathematical expressions. They discussed the difficulties in finding a suitable form to fit the very unusual kinetic properties of the ionic conductances: High order differential equations, chosen to mimic the delay in onset, would be unsatisfactory because of a similar delay upon repolarization. They solved the problem by using the property that the state variable of a first-order differential equation raised to a suitable power provided both a delayed onset and rapid reset with repolarization.

This formulation was so powerful and accurate that it is still - some four decades later - the most widely used set of equations (expressions) for describing the ionic channels in nerve cells. That this is true in spite of the fact that other models have been developed which give better fits to other data observed later is truly remarkable and enviable.

Fitting multiple types of data

Another measure of true professional style in modeling lies in testing the fit of simulations to as many different types of experimental data as possible. In their pioneering paper in 1952, Hodgkin & Huxley demonstrated monumental and exhaustive thoroughness in testing their model. The credibility of this work was immediately established because they were able to show, in amazing detail, an almost precise match of their model to their experimental observations:
  • Fig 12. Membrane action potentials generated with a range of initial voltage displacements
  • Fig 13. Hyperpolarization of the membrane following an action potential
  • Fig 14. Temperature induced changes in the membrane action potential and after- hyperpolarization
  • Fig 15. Form and velocity of a propagating action potential
  • Fig 16. Time course of the conductance change during an action potential - matching the observations of Cole and Curtis (1939)
  • Table 5. Net movement of sodium and potassium ions during an impulse
  • Fig 20 Refractory period following an impulse
  • Fig 21 Near threshold membrane action potentials
  • Fig 22. Reproduction of anode break excitation
  • Fig 23. Oscillations in membrane potentials in response to prolonged subthreshold current pulses

    This work must be the supreme example of craftsmanship in the field and, while every effort should be made to emulate it with each new system simulated, it may well be that no one else will be able to approach it's completeness.

    Another example of where multiple data sets were employed to achieve a "reasonable" representation of neuronal function comes from the work of Fred Dodge on a model of the spinal motoneuron of the cat. He took a Rallian electro-morphological description for soma and dendrites of this cell. From the work of Coombs, Curtis, & Eccles (1957) he postulated that the threshold in the axon was lower than that in the soma and adjusted the rate constants in the HH equations to model this. Furthermore he took voltage clamp data on this motoneuron to set the densities of the HH channels and sought credibility for this model by demonstrating that several independent experimental observations of antidromic and orthodromic propagation were reproduced by the simulations.

    Traub (1977) also has employed a variety of types of experimental data to set the density and types of channels in a model of the motorneuron: (a) firing frequency vs current injected into the soma, (b) the shape of the interspike voltage trajectory, (c) afterhyperpolarizations, antidromic spike invasion properties, etc.

    In unpublished simulations of a stylized motorneuron, I have adjusted the location and density of HH channels to fit experimental observations: (a) current observed during voltage clamp of the soma, and (b) the shapes and amplitudes of orthodromically generated spikes. In order to match the data on antidromic impulse invasion I had to search for the appropriate first internodal distance. However changes in the assignment of this parameter degraded the quality of fit to (a) or (b) or both. In order to achieve a reasonable fit to all experimental observations I had to iterate the parameter search cycle several times.

    Fitting algorithms

    There are many algorithms used to fit curves to data points. They all seek a point of minimum sum of the errors in the curve at each point (squared) but the search path differs. The quality of an algorithm consists in its ability to find the lowest global minimum in a terrain of small valleys, or "local" minima, which may give a good but not the best fit. Because some algorithms touted to be excellent sometimes stray far afield even when started at the known global minimum, testing is essential. One should test one's method of choice with data generated from known parameters and over a wide range of values for those parameters. When I have done such testing to see how well an algorithm could find the parameters underlying a given HH ionic conductance curve, I have been amazed at the poor quality of the values returned by the searches, even when I have given the correct values as the initial starting point!

    Apparent fitting of data can be deceptive and give unwarranted confidence. Examples:

    1. MISMATCH between the number of states in model and biological system

    True confidence in one's understanding of a system requires not only

    The praxis fitting algorithm has been used to demonstrate this problem. The accuracy of this algorithm is demonstrated by using it in a NEURON application program to find the parameters for a curve generated by a pair of known exponentials. Here it finds both the amplitudes and time constants precisely.

    However when the number of states in a model differs from that in the biological system, one can be deluded very easily into thinking that one has a good fit to data. In fact, there may be NO relation between the model parameters and those in the system. This is demonstrated in the figure below where data points taken from a decay curve generated by the sum of 3 exponentials is overlaid by the curve resulting from forcing a model having only 2 exponentials to fit those data. Here the three components generating the data points have equal amplitudes of 1 and decay time constants which increase in the ratios of 1: 3: 9. The "fitted" components have amplitudes (Amp1=1.23 & Amp2=1.56) and time constants (tau1=6.77 & tau2=0.37) as seen in the output panel below. These "fitted" parameters have no relation to the data values. Nevertheless the fact that the model curve overlays the data points so well often provokes an inappropriate confidence in one's knowledge of the properties of the data.

    Only when the ratio between the time constants is increased to 10-fold does the lack of quality in fitting become apparent. This 10-fold difference in the taus strains the algorithm and reveals the poor quality of the parameters found.

    A copy of this program is available and one may run it with one's own choices to examine the effect of) changes in relative amplitude and time constants, etc., of the generating function.

    2. The "Ramp Clamp": Facination with shapes can be Misleading

    For certain restricted rates of voltage sweep, the ramp clamp current-voltage relationship resembles the form of conventional I/V plots. For example, a very slow ramp has been used to plot the potassium current in a squid giant axon (Moore, 1959) (and in Cole & Moore, 1960). Because the ramp clamp was useful for obtaining a "quick and dirty" record under these restricted conditions, it was tempting to hope that a fast ramp may yield a plot resembling the sodium current.

    Simulations were obviously required and showed up the flaws and possible misinterpretations immediately. Early in the 1960's I ran analog computer simulations on a patch of HH membrane in order to examine the quality of the ramp clamp to generate realistic sodium current patterns. For other rates however, the shape of the curves deviated sharply because both the sdium and potassium conductances vary with time as well as voltage and BOTH contributed to the current observed at any time. I found the distortions so great and the confusing mixture of contributions of potassium to the "sodium" current that I did not judge it worth while to write a paper on a technique which I felt was completely useless.

    This is an example of where the verisimilitude of a hoped for match to data could give unwarrented confidence in a completely erroneous interpretation. Nonetheless the ramp clamp was later rediscovered and proposed as a useful way to record data very quickly because the output traces resembled the conventional peak and steady state plots of the step voltage clamp. Because of this I deem it necessary to show figures from ramp clamp simulations.Because my previous analog simulations were unpublished, a colleague, Kevin Martin, used NEURON to recreate them for me and are are shown here. A particular rate of voltage sweep (2mSec for -70mV to +70mV) was chosen to yield a plot resembling the sodium current. Of course for other rates, the shape of the curves deviated sharply because both the sodium and potassium conductances vary with time as well as voltage and BOTH contributed to the current observed at any time.

    The most obvious failing with the ramp clamp showed up when the sodium conductance had been completely blocked by application of tetrodotoxin (TTX). While the inward current disappeared for a fast ramp, an outward current persisted at voltages above the sodium equilibrium potential! Having done the original voltage clamp experiments showing that TTX blocked the sodium current in both directions, I knew that this outward current masquerading as Na current was specious. Simulations of this condition showed similar plots and revealed that the apparent "outward sodium current" was actually carried by potassium!

    Back to: Top of 3.1 Fitting
    Forward to: 3.2 Manual Iterative Searches
    Forward to: 3.3 Sensitivity Testing
    Forward to: 3.4 Numerical Accuracy