J Neurophysiol 99: 1871-1883, 2008.
First published February 6, 2008; doi:10.1152/jn.00032.2008
0022-3077/08 $8.00
Using Complicated, Wide Dynamic Range Driving to Develop Models of Single Neurons in Single Recording Sessions
Kevin H. Hobbs and
Scott L. Hooper
Neuroscience Program, Department of Biological Sciences, Ohio University, Athens, Ohio
Submitted 10 January 2008;
accepted in final form 31 January 2008
 |
ABSTRACT
|
|---|
Neuron models are typically built by measuring individually, for each membrane conductance, its parameters (e.g., half-maximal voltages) and maximal conductance value (gmax). However, neurons have extended morphologies with nonuniform conductance distributions, whereas models generally contain at most a few compartments. Both the original conductance measurements and the models therefore unavoidably contain error due to the electrical filtering of neurons and the differential placement of conductances on them. Model parameters (typically gmax values) are therefore generally altered by hand or brute force to match model and neuron activity. We propose an alternative method in which complicated, rapidly changing driving input is used to optimize model parameters. This method also ensures that neuron and model dynamics match across a wide dynamic range, a test not performed in most modeling. We tested this concept using leech heartbeat and generic tonically firing models and lobster stomatogastric and generic bursting models as targets and gmax values as optimized parameters. In all four cases optimization solutions excellently matched target activity. Complicated, wide dynamic range driving thus appears to be an excellent method to characterize neuron properties in detail and to build highly accurate models. In these completely defined targets, the method found each target's 8–13 gmax values with high accuracy, and may therefore also provide an alternative, functionally based method of defining neuron gmax values. The method uses only standard experimental and computational techniques, could be easily extended to optimize conductance parameters other than gmax, and should be readily applicable to real neurons.
|
|
INTRODUCTION
|
|---|
Conductance-based neuron modeling plays an increasingly important role in neuroscience. The central goal of such work is to build models that accurately describe neuron activity. However, present model-building techniques do not use neuron activity as their fundamental selection criterion. Models are instead typically constructed by first isolating individual conductances by a combination of pharmacological blockade of other conductances and the use of voltage command protocols that activate only one or a few conductances and various subtractive techniques. The conductances are then characterized on the basis of ion selectivity and the dependence of their opening and closing rates and steady-state values on voltage and, for some conductances, calcium. In these experiments maximum conductance (gmax) values are also typically measured.
The resulting equations are then placed in models typically containing one or a few compartments. In this first step, model and neuron activity often do not match well. One reason for this mismatch is the extended morphology of most neurons and the nonuniform distribution of membrane conductances on this morphology. No matter where a recording electrode is placed, some parts of the neuron are therefore electrically distant. Activity arising at these distant sites thus suffers from slow temporal filtering and amplitude attenuation. Conductance time constant measurements of electrically distant conductances will consequently be slower, and gmax measurements smaller, than they actually are.
These difficulties can theoretically be overcome. With present computational resources models with thousands of compartments are possible. Patch-clamp recordings of individual channels, or sharp electrode recordings from regions where a conductance is known to be present, can accurately measure conductance parameters. Advances in immunology and molecular biology hold the promise of accurately identifying conductance placement across neuron morphology. However, it is presently experimentally impossible to perform all these measurements (full description of a neuron's morphology, accurate description of the parameters defining all its conductances, accurate placement of the conductances across the neuron's processes) on an individual neuron.
Neuron descriptions are therefore typically data averaged from multiple neurons of a given type. For instance, at present neuron conductance properties and amounts are typically measured by pharmacological isolation. These treatments are time consuming and, for some conductances, irreversible. Characterizing all of a neuron's conductances during a single recording session is consequently generally impossible. Membrane conductance equations of a "neuron" are typically instead averages obtained from recordings of multiple neurons of the "same" neuron type, with one set of conductances being obtained from recordings from one group of neurons, another set of conductances from recordings from a second group of neurons, and so forth, until all the conductances of the "neuron" have been included.
The difficulty with this approach, however, is that recent work suggests that neurons of a given type can have very different maximum conductance value sets because maximum conductance values can change in a coordinated fashion that maintains neuron activity (Desai et al. 1999
; Franklin et al. 1992; Goldman et al. 2001; Golowasch et al. 1999a,b, 2002; Haedo and Golowasch 2006; Hong and Lnenicka 1995; Klein et al. 2003; LeMasson et al. 1993; Li et al. 1996; Liu et al. 1998; MacLean et al. 2003, 2005; Marder and Goaillard 2006; Marder and Prinz 2002; Mee et al. 2004; Mizrahi et al. 2001; Schulz et al. 2006; Swensen and Bean 2005; Thoby-Brisson and Simmers 2002; Turrigiano et al. 1995). Averaging across neurons could thus average neurons with very different maximum conductance values, and these average values may be far from those present in any member of the population (Marder and Goaillard 2006). Furthermore, modeling shows that the activity of these "average" neurons can substantially differ from the activity of the neurons from which the averages were derived (Golowasch et al. 2002). Averaging may thus not only fundamentally misrepresent a neuron type's maximum conductance makeup, but also result in models that do not reproduce neuron activity.
These articles deal only with maximum conductance values. But there is no reason to believe that similar compensatory mechanisms may not exist with respect to neuron morphology, conductance parameters (e.g., half-maximal voltage) other than maximum conductance, and placement of conductances on neuron processes. Neurons that are the "same" in that they produce the same activity and serve the same function could thus nonetheless have substantial differences in the details of their construction. Given all these factors, it would be astounding if models using average values of conductance parameters and amounts, placed in one compartment (or a few compartments), satisfactorily reproduced neuron activity. It is therefore typically necessary to tune (generally by altering gmax values) models, either by hand or brute force search, to find models acceptable to the user.
However, this raises another difficulty. The goal presumably is to find the approximation that best reproduces neuron activity. Implicit in this goal is to have in hand the best possible description of the neuron. Clearly, descriptions that include neuron responses to perturbation (e.g., current injection) are more complete than those of only rest activity. Presumably, the more complicated the input, the better the description. However, in most cases only modest alterations (e.g., responses to sustained current injection) are used, far from the most complicated alterations that can be introduced given present technology.
In considering these issues, we wondered whether a reverse, top-down, activity-based approach might result in both more specific neuron definition and more efficient and accurate model optimization. Although all the difficulties associated with imperfect conductance descriptions, and mapping complex neuron morphologies onto geometrically simpler models, still apply, such an activity-based approach would at least ensure that the overriding goal—accurately reproducing neuron activity—was always guiding the process. Such a procedure, if successful, would also allow models to be built for every neuron recorded from.
To test this method under completely defined conditions, we used already developed neuron models as targets. We drove these targets with complicated, rapidly changing input and used these data as input to an evolutionary algorithm optimization procedure. In all cases (four different models) the algorithm found gmax values that matched target activities extremely well when the equations used to describe target and model conductance were identical (e.g., with identical half-maximal voltages and dynamic parameters). The method also improved model performance when 1) noise was added; 2) when noise was added and the equations describing target and model conductances substantially differed; and 3) when noise was added, the equations describing target and model conductances substantially differed, and a conductance was omitted during the optimization. These data suggest that activity-driven optimization algorithms are a promising alternative to classic model-building methods.
The experimental procedures are rapid (minutes) and use only standard techniques. Single runs of the evolutionary algorithm required about 4 days using twenty 2-GHz processors. Several weeks were required for complete optimization, but this time would scale inversely with processor number and speed. Extending the method to optimize conductance properties as well as amounts, and to build multicompartment models, would be straightforward.
|
|
METHODS
|
|---|
Neuron simulation
Neuron simulation was conductance based and standard [see Numeric values and equations for generic and lobster burster (Prinz et al. 2003) models and Nadim et al. (1995) for the leech tonic equations]. All code was written in C using the GNU Scientific library version 1.7. Variable integration steps with Runge–Kutta stepping functions were used. Step size was continuously adjusted until all variable and derivative errors were just <10–3 (millisecond-based models, leech and lobster) or 10–6 (second-based models, both generics). To avoid very small step sizes at current injection or command voltage changes, if the change would have occurred within an otherwise acceptable step, step size was set so that step end occurred at the time of the change.
Neuron perturbation
The models were stimulated with a driving current in current clamp and a command voltage in voltage clamp. The stimulus sequence was generated once, randomly (from a uniform distribution), and was the same for all simulations, except for the data in Fig. 6. Current amplitude (all nA) distributions were –400 to 200 for the leech tonic, –0.4 to 0.2 for the lobster burster, and –20 to 10 for both generic models. Command voltage distributions were –100 to –30 mV for all models. Application of tetrodotoxin, tetraethylammonium, and 4-aminopyridine was simulated by setting gNa, gK1, gK2, and gA (leech tonic), gNa, gKd, and gA (lobster burster), and gNa, gKd, gAfast, and gAslow (both generic) to zero.
View larger version (38K):
[in this window]
[in a new window]
|
FIG. 6. Effect of adding noise and using a different stimulus pattern (A), of adding noise, using a different stimulus pattern, and using imperfect conductance descriptions (B1, B2), and of adding noise, using a different stimulus pattern, using imperfect conductance descriptions, and removing IA (C1, C2) on method performance. In A the method continued to find a model with essentially perfect activity (the second shown is the worst in the entire 30 s). In B and C the method could not exactly match target activity, but optimization still improved model performance. Black traces, target activity; red traces, model activity with target maximal conductance (gmax) values but incorrect or missing conductances; green traces, model activity after gmax optimization.
|
|
Evolutionary algorithm
The evolutionary algorithm was largely standard (Fogel 1994; Yao et al. 1999). The algorithm began on 1,000 model neurons ("individuals") whose maximum conductance values (xij, where i refers to the individual and j to the conductance type) were drawn randomly (within the conductance's ranges; see Numeric values and equations) from uniform distributions, one for each conductance type. Each individual's maximum conductances were associated with mutation size parameters (
ij) that were set in the initial population to one half of the relevant maximum conductance's range. The error scores (see RESULTS) of each individual were then calculated.
Each individual created a child as follows. The children first inherited mutated size parameters according to
where 'ij is the child's new size parameter; ij is the parent's size parameter; N(0, 1) is a number drawn from a standard normal distribution once per generation and used for every maximum conductance and individual in that generation; Nij(0, 1) is a number drawn anew from a standard normal distribution for each individual i and maximum conductance j; and 
' = 
and = (
)–1, where n is the number of maximum conductances. The ij then helped determine the child's mutated maximum conductances according to
where 
ij is a number randomly drawn for each individual i and conductance j from a Cauchy distribution with a scale parameter of 1. Cauchy distributions resemble Gaussian distributions but have longer, flatter tails and are thus more likely to generate numbers far from 0. If a mutation would cause a maximum conductance value to become negative the parameter was set to zero. Each child also inherited its parent's error.
Recombination among the children occurred as follows. A potential mate was chosen randomly from the entire population of children (i.e., a child could mate with itself) for one child. The errors the child and its mate inherited from their parents were then compared. If the child's error was greater than that of its mate, the child's maximum conductance values and mutation size parameters were replaced by a random assortment of both children's values (sexual reproduction). All values were randomly assorted; that is, a maximum conductance value could come from one child and the conductance's mutation size parameter from the other. If the child's error was less than or equal to that of its mate, no recombination occurred. A mate was then selected for the next child and the process repeated for all 1,000 children. Since mates were always chosen from the entire population of children, children that had undergone recombination in prior matings could serve as mates in subsequent ones.
At the end of this process the population consisted of 1,000 parents and 1,000 children. The errors of the mutated and recombined children were then calculated. Each individual was then compared with ten individuals chosen randomly from the entire (2,000) population and assigned a score equal to how many of the ten its error was less than or equal to. The 1,000 individuals with the highest scores (ties being resolved based on error) survived and became parents in the next generation. This procedure allows some individuals with relatively high errors to survive. Evolution continued until the population mean stopped moving, and the population distribution remained confined to a very small range, for about 50 generations.
The mutation size parameter and population error distribution interact as follows. If a parent's error is large relative to the population mean error, then this parent having a large mutation size parameter is generally advantageous because large mutations increase the child's likelihood of successfully competing against the rest of the population. If the parent's error is small relative to the population mean error, then its having a small mutation size parameter is advantageous because it produces children similar to the parent. This interaction results in mutation size parameters being large in early generations and then becoming very small as evolution proceeds and population mean error decreases.
Hardware
Computations were performed on local Beowulf Linux clusters and Beowulf Linux and OSX clusters at the Ohio Supercomputer Center (OSC). Message-passing interface was used to communicate between nodes; one node served as a master node. The master node ran the evolutionary algorithm and assigned individuals of the population one at a time to the other nodes, which performed the simulations and returned errors to the master. The clusters used were the Ohio University Physics (20 Intel Pentium 3 processors) and Biology (48 AMD Athlon MP processors) clusters and the OSC Pentium 4 (112 Intel Pentium 4 processors), Itanium 2 (476 Itanium 2 64-bit processors), and Apple G5 (64 PowerPC G5 processors) clusters. We never used more than 20 processors at a time on any cluster. We used 170,994 CPU hours on the OSC clusters and a similar but unrecorded number on the local clusters, although much of this time was spent developing the method.
Numeric values and equations
Target conductance units were microsiemens (µS) for the two generic models, millisiemens (mS)/cm2 for the lobster burster model, and nanosiemens (nS) for the leech heartbeat model. Target conductance values were gCaF = 16, gCaS = 5, gP = 3, gNa = 350, gK1 = 100, gK2 = 50, gA = 80, gh = 7, gleak = 10 (leech tonic); gleak = 0.05, gNa = 100, gCaT = 1, gCaS = 4, gA = 5, gCa = 15, gKd = 50, gh = 0.02 (lobster burster); gNa = 4,078, gKd = 2.06, gAslow = 6.029, gh = 0.4539, gK = 0.2485, gKCa = 27.41, gAHP = 0.0303, gCaf1 = 1.083, gCaf2 = 0.3573, gCas = 0.0469 (generic tonic); gleak = 0.1, gNa = 1,500, gKd = 0.5, gAfast = 0.07, gAslow = 0.6, gh = 4, gKv3 = 0.06, gK = 0.2, gKCa = 1, gAHP = 0.3, gCaf1 = 0.4, gCaf2 = 0.7, gCas = 0.002 (generic burster).
Initial maximum conductance ranges for the evolutionary algorithm were 0–100 for all leech tonic conductances; 0–0.05 (gleak, gh), 0–500 (gNa), 0–12.5 (gCaT), 0–10 (gCaS), 0–50 (gA), 0–25 (gKCa), 0–125 (gKd) for the lobster burster (as in Prinz et al. 2003); and 500–1,000 (gNa) and 0–10 (all others) for both generic models.
The generic models used the following equations
The lobster burster model equations are identical to those in Prinz et al. (2003). They are included here because this model was the one in which equation parameters (e.g., half-maximal voltages) were altered to test whether the algorithm would continue to find improved solutions even when the equations describing the conductances being used by the algorithm and those describing the target conductances substantially differed. In these cases (Fig. 6, B and C) equation parameters were modified by sequentially going through the equations parameter by parameter and alternately multiplying the parameter by 1.1 or 0.9, except in all conductances but Ih, the first term in the m equations was always multiplied by 1.1 to avoid m becoming negative at certain membrane voltages. In the following equations a parameter was multiplied by 1.1 when followed by a "u" and by 0.9 when followed by a "d."
iNa
iCaT
iCaS
iA
iKCa
iKd
ih
[Ca2+] dynamics
|
|
RESULTS
|
|---|
Neuron electrical activity results from the combined activity of a neuron's membrane conductances. We reasoned that the reverse—that a neuron's activity can be used as the fundamental driver in building models of the neuron—should also be true. Key to this process is altering neuron activity. That is, many different models can reproduce the activity of a neuron whose rest activity is an unvarying –60 mV membrane potential. When current is injected into the neuron, its responses reduce the number of models that reproduce its activity. For instance, if the neuron fires when depolarized, models with too little INa to support action potentials are excluded. More complex perturbation further decreases the number of models that can reproduce neuron responses. The goal is to provide sufficiently complex input that only one model reproduces neuron response.
The most complex input per time is continuously varying random input. We approximated continuously varying input with a series of step changes. Shorter steps give a better approximation. However, very brief steps are difficult to implement experimentally due to membrane and electrode time constants and amplifier settling times. We chose a step short enough (50 ms) to achieve large step numbers in relatively brief recording times, but long enough to be easily implemented with almost all neurons in both current and voltage clamp, and varied only step amplitude.
Target response to input was measured in current and voltage clamp in normal physiological saline and, to provide an additional condition involving changed conductance makeup, saline containing blockers of INa, IK, and the A-type potassium current, IA. Input amplitudes were chosen to induce physiologically relevant and experimentally achievable responses. Current-clamp current injections induced membrane potentials between –125 mV and action potential threshold. To avoid experimental difficulties with action potential control, voltage-clamp commands varied from –125 to approximately 10 mV below spike threshold.
Three error measures were used. Voltage (in current clamp) and current (in voltage clamp) response errors were the areas (mV x s and nA x s) between target and model activity (gray shading, Fig. 1, A–D). Spike time error was the sum of the absolute values of the times between model and nearest target spikes (Fig. 1E1) plus the sum of the absolute values of the times between target and nearest model spikes (Fig. 1E2). Measuring these times both target to model and model to target was required to account for all missing and extra spikes (note that Fig. 1E1 has two arrows and Fig. 1E2 three).
View larger version (20K):
[in this window]
[in a new window]
|
FIG. 1. Error functions measured the area (gray shading) between target and model activity in normal saline in current (A) and voltage (C) clamp and in a drug cocktail (B, D). Bottom panels are current injection (A, B) or command voltage (C, D) trajectories. Times (arrows) from model to target (E1) and target to model (E2) spikes were also measured. All panels: black line, target activity; gray line, activity of a model at an intermediate generation.
|
|
These errors determined survival in an evolutionary algorithm (Fogel 1994; Yao et al. 1999) described in detail in METHODS. An evolutionary algorithm was chosen because it does not require calculating error derivatives, is resistant to being trapped in local minima, and is easily implemented on distributed parallel computers. In brief, it worked on an initial population of 1,000 model neurons with different maximum conductances. Each individual produced one child. Mates were chosen randomly. The pair reproduced sexually if the reproducing individual's error was greater than that of its mate; otherwise, the individual reproduced asexually. Maximum conductances were inherited with mutation. Mutation amplitude was determined by random drawing from a Cauchy distribution and an inherited mutation size parameter. One thousand individuals from the combined parent and child population were then selected to be new parents using a procedure that favored low-error individuals but allowed some high-error individuals to survive. A different error measure and experimental condition (saline current clamp, saline voltage clamp, saline spike time, drug cocktail current clamp, drug cocktail voltage clamp) were used each generation. Note that the algorithm has no knowledge of target gmax values; survival is determined only by activity-based error measures.
It is important to note that the method presented here thus consists of three interacting parts: complicated driving input (Fig. 1), error functions (Fig. 1), and the algorithm used to determine individual survival based on these errors. To make this presentation as clear as possible, we use the term "method" to refer to all three processes operating as a whole and "algorithm" to apply only to the evolutionary algorithm itself.
To test the method against a wide variety of conductance sets, we used multiple neuron models as targets. One (leech tonic firer) had nine conductances whose properties (e.g., half-maximal activation voltage) and maximum conductance values were measured in leech heart interneurons (Nadim et al. 1995). A second (lobster burster) had eight conductances whose properties were measured in lobster stomatogastric neurons and whose maximum conductances produce an endogenous burster (Prinz et al. 2003). A third (generic tonic firer) had ten conductances modified from measurements in a variety of neurons. A fourth (generic burster) had these ten conductances with different maximum conductance values and an additional three conductances. Because the leech, lobster, and generic models were developed from different species, somewhat different half-maximal activation voltages and time constants were used even for conductances of the same type (e.g., Ih). The four models thus differed not only in maximum conductance values, but also in which conductances were present in each model and in each conductance's properties.
Maximum conductance values of the initial population were assigned randomly across initial ranges either arbitrarily chosen or, for the lobster burster, those used in the model's original source (Prinz et al. 2003). These choices were not particularly important because evolutionary algorithm performance depends relatively little on initial conditions (Fogel 1994; Yao et al. 1999) and because we used the results of earlier runs to define the initial ranges of later runs (see following text).
As we subsequently show, the method finds solutions that match target activity extremely well. In practice, the gmax values of the neuron being modeled would be unknown, and thus a user would halt the optimization either 1) because solutions with a predetermined level of error had been obtained (at which point multiple solutions with equally small error may be present) or 2) because further optimization did not decrease solution error (at which point all solutions are presumably very near or at the global minimum of the parameter space, and thus there is little or no further improvement the algorithm can achieve). However, since we were testing method performance and had foreknowledge of target gmax values, we could test the method's ability to find the global minimum by explicitly comparing method to target gmax values. We therefore first present the data in terms of gmax values (Figs. 2–4) and only then (Fig. 5) present activity data.
View larger version (29K):
[in this window]
[in a new window]
|
FIG. 2. Time course of 5 runs from the initial naive range of the evolutionary algorithm for the generic tonic target. Small gray dots show maximum conductance of the delayed rectifier (gKd) for the entire population, and black line the population mean, generation by generation, for one run. Colored lines are population means for the other 4 runs. Runs ended with multiple solutions (inset) grouped around the target (dashed line). The dots in the inset are larger than the range of values present in the population at each run's end.
|
|
View larger version (10K):
[in this window]
[in a new window]
|
FIG. 3. Comparison of target and method final run maximum conductance values. In all panels x is the target, the circle is the average, and the error bars the SD (in all cases too small to be seen) of the final set of runs; the small dots are the final values of the individual runs (in most cases so closely grouped that they appear as single dots). Note logarithmic scale for conductance values.
|
|
View larger version (12K):
[in this window]
[in a new window]
|
FIG. 4. Comparison of target and method final run maximum conductance values normalized to target values: 0% line marks target; other conventions as in Fig. 3.
|
|
View larger version (12K):
[in this window]
[in a new window]
|
FIG. 5. Comparison of 1 s of target and method activity using the runs with the smallest current-clamp saline error. The second shown has the greatest error in the entire simulation. For all 4 targets the 2 activities are indistinguishable on this scale. Method activity is therefore shown with a thick gray line.
|
|
In all cases the populations evolved similarly. Figure 2 shows the evolution of gKd values for the first five runs of the generic tonic target. The gray dots are the maximum conductance values of the entire population for one run. The black line is the population mean. The range of maximum conductance values first expanded (from 0–10 to 0–10,000 µS in the first 50 generations) and then collapsed with very small variation around a single value (runs typically continued for 300–500 generations, and final ranges were thus smaller than that shown here). This variation was so small that all members of the population had essentially equal errors. We therefore used the population's mean maximum conductance values as the run's final value.
The other runs from this initial maximum conductance range (colored lines) collapsed onto final values (inset) whose range (2–3.5 µS) was smaller than the initial range (0–10 µS). This evolution to a set of final values with a smaller range than the initial range occurred simultaneously for all model conductances. These new ranges were used as initial ranges for a second set of runs, which in turn collapsed onto a more tightly grouped set of final values. Each repetition typically resulted in a 106- to 1014-fold decrease in the hypervolume defined by the range of maximum conductance values. We performed two to six such repetitions, each with at least five runs.
This process found maximum conductance values that showed excellent agreement with all four targets (Fig. 3). However, the very large y-axis ranges (up to six orders of magnitude) in this figure make it difficult to see relative differences between method and target values. We therefore normalized the data to target values (Fig. 4). These again showed excellent agreement, with mean method values of 38 of the 40 conductances falling within ±1% of the target and the remaining two within ±5%. Considering the values of the individual runs from which these averages were calculated, the values of all runs were within ±1% of the target for 32 conductances, ±5% for 37, and ±10% for all 40.
Although these data show that the method finds solutions very close to target maximum conductance values, as noted earlier, in practice these values would be unknown. The truly relevant comparison is thus between target activity and the activities of the method solutions. These comparisons show extremely good matches, with mean voltage differences (in current-clamp saline) between method and target activity of <1 mV for all 23 solutions, 0.5 mV for 22, and 0.1 mV for 18. This extremely good match is shown in Fig. 5, which plots the worst second of the best solution for each target; the two activities are indistinguishable for all four targets.
Although the data presented earlier show the method is successful under perfect conditions, real neurons present a number of additional problems, and it was important to test that the method also functioned well under more realistic conditions. One such problem is that other researchers will use different pseudorandom stimulus patterns, and it was therefore important to ensure that the method would work with other pseudorandom stimulus patterns. A second problem is that real neuron recordings inevitably contain extrinsic noise. To test whether these difficulties destroyed the method's ability to find solutions that well reproduce neuron activity, we used a different pseudorandom stimulus pattern and added noise to the lobster burster traces in both current and voltage clamp in control saline and the drug cocktail, and reran the evolutionary algorithm from intermediate initial ranges (those used for the second iteration in the original set of runs). Adding noise only very slightly decreased the ability of the algorithm to find target activity (Fig. 6A shows the 1-s time interval with highest error) and the gmax values found differed from the target by <0.5%.
A potentially greater problem is that (as noted in the INTRODUCTION) it is likely that not only directly measured gmax values, but also other parameters (e.g., half-maximum voltage) in the equations used to describe the neuron conductances, will have errors. It was thus possible that the method might fail to find improved solutions if the algorithm used equations in which these parameters differed from those in the target. For instance, the actual half-maximal voltage of a neuron's INa activation could be –28 mV, but the value measured from the cell body—the value that would be used in the conductance's equations—would be –25.5 mV. To test the effect of such errors on method performance we altered all equation parameters by increasing or decreasing them by 10% (see METHODS) and reran the lobster burster model using target gmax values and the new conductance equations. Model activity of course no longer matched target activity because of the changes in the conductance equations. We then reran the evolutionary algorithm (again with added noise, the new stimulation pattern, and from intermediate initial ranges, as in Fig. 6A), with the algorithm using the altered conductance equations, to test whether it would find new gmax values that produced model activity more closely matching the activity of the target (Fig. 6, B1 and B2). The black traces show target activity with the original gmax values and original conductance equations. The red traces show the activity of a model with the target gmax values and the modified conductance parameters.
The model with target gmax values and modified conductance parameters (red trace, Fig. 6B1) missed some bursts (those at 5, 10.5, and 28 s), had difficulties (four misplaced bursts where there were originally three) from 11 to 14 s, and had a badly timed burst at 23 s. The evolutionary algorithm, working with the conductance equations with modified parameters, optimized gmax values sufficiently to cause the following improvements in model performance (green traces): the missing burst at 10 s was restored, the three bursts between 11 and 14 s were restored, and the timing of the burst at 23 s became closer to the target burst time. These gmax changes also strongly improved spiking activity within the bursts (Fig. 6B2). Using target gmax values and modified conductance parameters (red trace) produced single spike bursts with a pronounced shoulder. The gmax values achieved after the optimization with the modified conductances produced activity (green traces) with spike number and timing much closer to target activity, as were also the slow wave membrane potential variations between the bursts.
When using equations with parameters different from those in the target, the optimization was not able to find a set of gmax values that exactly reproduced target activity (for instance, the bursts at 5 and 28 s are still missing, and intraburst spiking is not perfect). However, this is not surprising given that the optimization was working with incorrect conductances. What is important is that the method still improved model performance, and thus the optimization procedure continues to function even when using equations that are different from those that produced the target activity. To do so with mismatched conductance equations, of course, resulted—required—in the method finding gmax values different from those in the target (see DISCUSSION for further consideration of this issue). However, these are gmax values that, with the mismatched conductances, improved model activity.
An additional potentially serious error that could occur is that the description of the neuron being modeled is missing a conductance, and thus this conductance is not included in the model and hence not available for optimization. To test whether this would cause the method to fail, we removed IA from the conductances the algorithm had at its disposal and reperformed the steps used to produce the data in Fig. 6, B1 and B2, again using the lobster burster model with added noise, the new stimulation pattern, and mismatched conductances (Fig. 6, C1 and C2). Interestingly, removing IA improved nonoptimized (red trace; compare with red trace in Fig. 6B1) burst performance, missing only the bursts at 5 and 28 s and lacking the difficulties between 11 and 14 s. Optimization slightly improved burst performance (e.g., improved burst time at 16 s) (Fig. 6C1). However, it again strongly improved intraburst spiking performance (Fig. 6C2). Thus even when completely lacking conductances present in the target, the optimization method still improved model activity.
|
|
DISCUSSION
|
|---|
The primary goal of this work was to determine whether, using presently available computational resources, neuron response to experimentally feasible perturbation contains sufficient information, using already defined conductance types but allowing maximum conductance amounts to vary, to find models that well reproduce neuron activity. The almost perfect correspondence between method derived and target activity (Fig. 5) shows that this is true for the four targets chosen. Four examples do not prove generality. However, the models chosen include a wide variety of conductances (Ileak, INa, IKd, IP, IA, IKv3, IAfast, IAslow, IAHP, IKCa, Ih, ICaS, ICaF, ICaT) present over a very wide (105) gmax range, which suggests that this method may be useful for building models of a wide variety of neurons.
It is essential to stress the central importance of perturbation of neuron activity and of stringent error measures to the method's success. In early attempts with less complicated perturbation and less stringent error measures the method found widely separated solutions with approximately equally small error. Extensive perturbation is not typically performed when describing neuron properties and testing model adequacy. Our experience suggests that, as is common in analysis of other physiological systems (Marmarelis and Marmarelis 1978), it should be.
Application to real neurons
The ability of the method to continue to improve model activity when faced with noise, imperfect conductance descriptions, and even a missing conductance supports the argument that the method can profitably be applied to real neurons. The fact that the method could not find perfect solutions in the latter two cases suggests that in real neurons it will also not find solutions that perfectly reproduce neuron activity because of the difficulties in conductance characterization and identification mentioned in the INTRODUCTION. However, the central role of neuron activity in the method means that models built using the maximum conductance values it does find should be near the best possible with respect to reproducing neuron activity. Furthermore, the central role of complex perturbation in the method means that these models should also well reproduce neuron activity for all inputs within the perturbation's dynamic and amplitude limits. The ability of the method to determine presumably near-best gmax values for a given set of conductances and conductance equations suggests that, if method solutions do not meet user standards, then model improvement likely requires changes in conductance makeup, equation parameters, or multicompartmentalization.
Model performance can be further enhanced by expanding the parameters it optimizes. For instance, allowing it to adjust conductance time constants and voltage dependencies would presumably further decrease model error. Even more accurate representations could be obtained using multicompartment models. A major difficulty in such models is correctly placing the conductances in the compartments. Since this process simply represents extra parameters, modifying the method described here to identify which placings best reproduce neuron activity is straightforward. In both cases the primary difficulty is the increased parameter space that must be searched. However, experimental data to guide initial ranges (e.g., that spike-generating conductances are located only in distal compartments) are often known, which would greatly reduce computation time.
What does "best" mean with respect to model parameters and neuron description?
As explained in the INTRODUCTION, all measurements of conductance properties and amounts measured from one site (often the neuron cell body), inherently and unavoidably contain errors due to morphological filtering and attenuation. It is to overcome this difficulty that hand tuning or brute force parameter optimization, or the method described here, is required. The method described here, particularly if extended to optimizing conductance parameters and placement as suggested earlier, may seem to some readers insufficiently grounded in "hard" data. That is, it may seem that descriptions of conductance time constants, half-maximal voltages, and amounts derived from direct measurements of individually isolated conductances are more nearly "correct" than those derived from a method that finds those values simply by matching neuron and model activity.
However, it cannot be overstressed that in each case the values obtained are gross oversimplifications. For instance, the very small amounts of INa that standard techniques measure for neurons in which the sodium channels are electrically distant from the cell body are certainly not "correct" for action potential generation in a single-compartment model. The method presented here is specifically designed to find maximum conductance amounts (or, if extended, other conductance properties and compartment placement) that best reproduce neuron activity, and reproducing neuron activity is the ultimate test of whether a model is correct. With respect to building neuron models, we would thus argue that the values found by the optimization method described here or similar methods are more nearly correct.
With respect to describing a neuron's actual membrane complement, which technique gives a truer description depends on the neuron being investigated and investigator goal. For electrically compact neurons standard techniques and the activity-dependent technique described here will give similar (and true, in that little filtering is occurring) results. The ability of the activity-dependent technique to measure all neuron parameters from individual neurons would thus make it the preferred method for such neurons. For neurons with extended anatomies and unequal placement of membrane conductances across the neuron processes, standard techniques will give low-pass-filtered and reduced amplitude (relative to what is actually present where the conductance is located) measurements of distant conductances. The activity-dependent technique will give different values, finding those necessary to reproduce neuron activity. In both cases, more accurate descriptions of neuron anatomy will result in more accurate descriptions of neuron conductances. Given that each is a cluster of approximations, which is better is a matter of opinion, but the ability of the activity-dependent technique to approximate all membrane conductances simultaneously suggests that in many cases it would still be preferred. Membrane conductances are often measured to compare different neuron types (see following text for discussion of comparison within a neuron type), often with the implicit or explicit goal of understanding why the neurons have different activities. In these cases, given its fundamental dependence on activity, the activity-dependent technique would likely again be preferable.
Can multiple conductance sets make the same neuron type?
Whether neurons of the "same" type—that is, neurons with similar intrinsic activity and responses to input—can be constructed with different maximum conductance amounts remains controversial (Nowotny et al. 2007; Szüchs and Selverston 2006). The method described here can help to resolve this issue in three ways. First, the neurons in question have never been subjected to the complicated perturbations used here. It is possible that when the neurons are so perturbed they will have several different types of response (for each of which the method would find different, statistically separate, conductance value sets as solutions). In this case this issue would be resolved in that these neurons actually constitute several separate neuron classes and have appeared to be similar only because they have not been examined in sufficient detail. Second, the neurons may instead show a graded continuum of responses centered on a mean. In this case the method would find a different best solution for each neuron, with the range of solutions found reflecting the variability across the group of neurons. Third, it is possible that all the neurons show similar responses, but the method finds multiple, statistically separate, solutions equally well reproduce the activity of each neuron. This would indicate that truly different combinations of conductances can produce very similar activity even with strong perturbation but, because the method uses only activity to determine error and never directly measures conductance values, would not prove that the neurons actually have different conductance sets. To prove this would require using classical (e.g., voltage-clamp) techniques to measure individual conductance values in the individual neurons. However, having identified which sets of conductances are capable of producing the observed activity would make this last step rapid and easy compared with present approaches.
It is important to make one final comment on this issue. Even in the case in which the neurons fall into separate classes when strongly perturbed (case 1 in the preceding paragraph), this does not necessarily mean the neurons are not functionally the same. It is unclear that real neurons would ever need to respond "correctly" to input this complicated. If real neurons need to perform only a more limited set of tasks, neurons with different conductance compositions, which our method would identify as being different, could nonetheless be functionally equivalent.
|
|
GRANTS
|
|---|
This work was supported by an Ohio Supercomputer Center (OSC) Cluster Ohio Initiative grant for the Ohio University Department of Biology parallel computing facility, support from Ohio University Computer and Network Services, grants totaling 20,000 resource units on the OSC facility, a guest professorship of the University of Köln to S. L. Hooper, National Science Foundation Grant 90250, and National Institute of Mental Health Grant MH-57832.
|
|
DISCLOSURE
|
|---|
The authors have no potential conflicts of interest.
|
|
ACKNOWLEDGMENTS
|
|---|
We thank R. A. DiCaprio for extensive discussion and comments on the manuscript and the Ohio University Physics Department for access to its proprietary cluster.
|
|
FOOTNOTES
|
|---|
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
Address for reprint requests and other correspondence: S. L. Hooper, Department of Biological Sciences, Irvine Hall, Ohio University, Athens, OH 45701 (E-mail: hooper{at}ohio.edu)
|
|
REFERENCES
|
|---|
Desai NS, Rutherford LC, Turrigiano GG. Plasticity in the intrinsic excitability of cortical pyramidal neurons. Nat Neurosci 2: 515–520, 1999.[CrossRef][Web of Science][Medline] Fogel DB. An introduction to simulated evolutionary optimization. IEEE Trans Neural Netw 5: 3–14, 1994.[CrossRef][Web of Science][Medline]
Franklin JL, Fickbohm DJ, Willard AL. Long-term regulation of neuronal calcium currents by prolonged changes of membrane potential. J Neurosci 12: 1726–1735, 1992.[Abstract]
Goldman MS, Golowasch J, Marder E, Abbott LF. Global structure, robustness, and modulation of neuronal models. J Neurosci 21: 5229–5238, 2001.[Abstract/Free Full Text]
Golowasch J, Abbott LF, Marder E. Activity-dependent regulation of potassium currents in an identified neuron of the stomatogastric ganglion of the crab Cancer borealis. J Neurosci 19: RC33 (1–5), 1999a.[Abstract/Free Full Text]
Golowasch J, Casey M, Abbott LF, Marder E. Network stability from activity-dependent regulation of neuronal conductances. Neural Comput 11: 1079–1096, 1999b.[CrossRef][Web of Science][Medline]
Golowasch J, Goldman MS, Abbott LF, Marder E. Failure of averaging in the construction of a conductance-based neuron model. J Neurophysiol 87: 1129–1131, 2002.[Abstract/Free Full Text]
Haedo RJ, Golowasch J. Ionic mechanisms underlying recovery of rhythmic activity in adult isolated neurons. J Neurophysiol 96: 1860–1876, 2006.[Abstract/Free Full Text]
Hong S, Lnenicka G. Activity-dependent reduction in voltage-dependent calcium current in a crayfish motoneuron. J Neurosci 15: 2539–3547, 1995.
Klein JP, Tendi EA, Dib-Hajj SD, Fields RD, Waxman SG. Patterned electrical activity modulates sodium channel expression in sensory neurons. J Neurosci Res 74: 192–198, 2003.[CrossRef][Web of Science][Medline]
LeMasson G, Marder E, Abbott LF. Activity-dependent regulation of conductances in model neurons. Science 259: 1915–1917, 1993.[Abstract/Free Full Text]
Li M, Jia M, Fields RD, Nelson PG. Modulation of calcium currents by electrical activity. J Neurophysiol 76: 2595–2607, 1996.[Abstract/Free Full Text]
Liu Z, Golowasch J, Marder E, Abbott LF. A model neuron with activity-dependent conductances regulated by multiple calcium sensors. J Neurosci 18: 2309–2320, 1998.[Abstract/Free Full Text]
MacLean JN, Zhang Y, Goeritz MI, Casey R, Oliva R, Guckenheimer J, Harris-Warrick RM. Activity-independent co-regulation of IA and Ih in rhythmically active neurons. J Neurophysiol 94: 3601–3617, 2005.[Abstract/Free Full Text]
MacLean JN, Zhang Y, Johnson BR, Harris-Warrick RM. Activity-independent homeostasis in rhythmically active neurons. Neuron 37: 109–120, 2003.[CrossRef][Web of Science][Medline]
Marder E, Goaillard JM. Variability, compensation and homeostasis in neuron and network function. Nat Rev Neurosci 7: 563–574, 2006.[CrossRef][Web of Science][Medline]
Marder E, Prinz AA. Modeling stability in neuron and network function: the role of activity in homeostasis. Bioessays 24: 1145–1154, 2002.[CrossRef][Web of Science][Medline]
Marmarelis PZ, Marmarelis VZ. Analysis of Physiological Systems. New York: Plenum, 1978.
Mee CJ, Pym ECG, Moffat KG, Baines RA. Regulation of neuronal excitability through pumilio-dependent control of a sodium channel gene. J Neurosci 24: 8695–8703, 2004.[Abstract/Free Full Text]
Mizrahi A, Dickinson PS, Kloppenburg P, Fenelon V, Baro DJ, Harris-Warrick RM, Meyrand P, Simmers J. Long-term maintenance of channel distribution in a central pattern generator neuron by neuromodulatory inputs revealed by decentralization in organ culture. J Neurosci 21: 7331–7339, 2001.[Abstract/Free Full Text]
Nadim F, Olsen Ø, de Schutter E, Calabrese RL. Modeling the leech heartbeat elemental oscillator. I. Interactions of intrinsic and synaptic currents. J Comput Neurosci 2: 215–235, 1995.[CrossRef][Web of Science][Medline]
Nowotny T, Szüchs A, Levi R, Selverston AI. Models wagging the dog: are circuits constructed with disparate parameters? Neural Comput 19: 1985–2003, 2007.[CrossRef][Web of Science][Medline]
Prinz AA, Billimoria CP, Marder E. Alternative to hand-tuning conductance-based models: construction and analysis of databases of model neurons. J Neurophysiol 90: 3998–4015, 2003.[Abstract/Free Full Text]
Schulz DJ, Goaillard JM, Marder E. Variable channel expression in identified single and electrically coupled neurons in different animals. Nat Neurosci 9: 356–362, 2006.[CrossRef][Web of Science][Medline]
Swensen AM, Bean BP. Robustness of burst firing in dissociated Purkinje neurons with acute or long-term reductions in sodium conductance. J Neurosci 25: 3509–3520, 2005.[Abstract/Free Full Text]
Szüchs A, Selverston AI. Consistent dynamics suggests tight regulation of biophysical parameters in a small network of bursting neurons. J Neurobiol 56: 1584–1601, 2006.
Thoby-Brisson M, Simmers J. Long-term neuromodulatory regulation of a motor pattern-generating network: maintenance of synaptic efficacy and oscillatory properties. J Neurophysiol 88: 2942–2953, 2002.[Abstract/Free Full Text]
Turrigiano G, LeMasson G, Marder E. Selective regulation of current densities underlies spontaneous changes in the activity of cultured neurons. J Neurosci 15: 3640–3652, 1995.[Abstract]
Yao X, Liu Y, Lin G. Evolutionary programming made faster. IEEE Trans Evol Comput 3: 82–102, 1999.[CrossRef]
This article has been cited by other articles:
|
|
|
|
 
A. L. Taylor, J.-M. Goaillard, and E. Marder
How Multiple Conductances Determine Electrophysiological Properties in a Multicompartment Model
J. Neurosci.,
April 29, 2009;
29(17):
5573 - 5586.
[Abstract]
[Full Text]
[PDF]
|
|
|
Copyright © 2008 by the The American Physiological Society.
TOP
ABSTRACT
INTRODUCTION