Cellular function given parametric variation in the Hodgkin and Huxley model of excitability.

How is reliable physiological function maintained in cells despite considerable variability in the values of key parameters of multiple interacting processes that govern that function? Here, we use the classic Hodgkin-Huxley formulation of the squid giant axon action potential to propose a possible approach to this problem. Although the full Hodgkin-Huxley model is very sensitive to fluctuations that independently occur in its many parameters, the outcome is in fact determined by simple combinations of these parameters along two physiological dimensions: structural and kinetic (denoted S and K, respectively). Structural parameters describe the properties of the cell, including its capacitance and the densities of its ion channels. Kinetic parameters are those that describe the opening and closing of the voltage-dependent conductances. The impacts of parametric fluctuations on the dynamics of the system-seemingly complex in the high-dimensional representation of the Hodgkin-Huxley model-are tractable when examined within the S-K plane. We demonstrate that slow inactivation, a ubiquitous activity-dependent feature of ionic channels, is a powerful local homeostatic control mechanism that stabilizes excitability amid changes in structural and kinetic parameters.

The canonical Hodgkin–Huxley mathematical model of membrane excitability is embedded in a high-dimensional parameter space. In their original report, Hodgkin and Huxley (1) indicate that the parameters vary substantially between different cells, an observation that is extensively documented by electrophysiologists who have studied excitable membranes over half a century (2, 3). This cell-to-cell intrinsic variability—that is, variability that cannot be attributed to measurement uncertainties (2, 4)—is habitually averaged out in attempts to generalize findings (5, 6).

There are many reasons to assume that cell-to-cell variability is also expressed in a given cell over time (7–13). For instance: kinetic parameters of channel gating change due to continuous modulation and activity-dependent roaming in protein configuration space; the ratio between the number of different channel proteins in the membrane might change due to differential protein expression or turnover; the membrane capacitance and leak conductance change during massive cell growth, movement, or contact of the cell with biological matrices that impact on membrane surface tension.

Randomly and independently pulled from the physiological range indicated by Hodgkin and Huxley, many combinations of parameters give rise to an excitable “solution” (i.e., stimulus-driven excitable membranes that generate single action potentials), but many other combinations lead to either nonexcitable or oscillatory, pacemaking membranes (14). The solutions are sensitive to relatively slight parametric independent variations. This is in contrast to the biological neurons and muscles that maintain relatively invariant patterns of activity that are seemingly more robust compared with the classic Hodgkin–Huxley-type models used to describe them. Cell-to-cell and within-cell parametric variation challenge our understanding of establishment and maintenance of excitability, as well as the methods we use to extract parameters from voltage-clamp data and construct suitable models (2, 15–17). This essential problem goes beyond the regulation of excitability; it belongs to a class of open questions that concern the study of organization in biological systems and the emergence of macroscopic functional order from a large space of potential microscopic “disordered” configurations (18).

Attempts to account for invariant excitability given parametric variation focus on activity-dependent, homeostatic coordination of channel protein expression (11, 19–21). The interpretation is corroborated by correlations between mRNA concentrations of different ionic channel proteins (22). It is also supported by elegant simulations showing how centralized activity-dependent (feedback) regulation that controls protein expression may navigate a cell into one of many functional solutions (11, 23). However, channel protein densities are not the only determinants of membrane excitability status. Even in the relatively simple Hodgkin–Huxley model—a single compartment with two voltage-dependent conductances—more than 10 parameters are involved, possibly varying and impacting each other in a wide range of time scales (subsecond to hours and days).

We acknowledge the difference between the high-dimensional parameter space dictated by the explicit Hodgkin–Huxley model and the dimensionality of the physiological space within which regulation of excitability is embedded. In other words, we ask how many physiologically relevant dimensions are needed to capture the dynamics of excitability and its regulation in the Hodgkin–Huxley formulation, as this may be very different from the number of free parameters in the full Hodgkin–Huxley model. From the physiologist’s perspective, the actual (hopefully not too many) dimensions should be expressed in parameters that can be directly extracted from standard voltage-clamp data. The approach we take below maintains the biophysically measurable parameters and is therefore different from most other reductions previously done.

We show that—congruent with low-dimensional models of excitability (24–28)—the phenomenon of excitability may be reduced to two dimensions. We identify these dimensions as cellular-level structural (denoted S) and protein-intrinsic kinetic (denoted K) dimensions and express them as combinations of actual Hodgkin–Huxley parameters that may be extracted from voltage-clamp data. Structural parameters refer to membrane surface area, number of different voltage-dependent conductances, and ionic concentrations. Kinetic parameters are those that describe the rate constants associated with channel opening and closing. Reduced to the S–K space, the manifold of functional solutions is simpler to understand amidst parametric variations, enabling regulation of excitability by one activity-dependent principle that is tightly related to a ubiquitous physiological process: slow inactivation of ionic channels.

Results and Discussion

Hodgkin–Huxley Model, Multiple Solutions, and Their Sensitivity to Parametric Variations.

The dispersion of parameters measured by Hodgkin and Huxley is summarized in Fig. 1—a collage of data and images from the original report, adapted to modern notation. The table (Fig. 1) indicates ranges of cellular-level structural parameters. The term “structural” is used, as these parameters are fully determined by physical measures that are characteristic of the cell. These include membrane surface area, the number of voltage-dependent conductances, and relevant equilibrium potentials. The graphs of Fig. 1 depict protein-level kinetic parameters, expressed as six transition rate functions superposed with data points. The mathematical expressions of these six rate functions—each of which describes the change of transition rate with membrane voltage ()—involve >10 different “hidden” parameters. Note the dispersion of points around the fitted rate functions, depicting repeated measurements in different axons.

Fig. 1.

A collage of slightly modified images extracted from the 1952 original report of Hodgkin and Huxley (1). (A) The table indicates ranges of cellular-level structural parameters. (B) The graphs depict protein-level kinetic parameters, expressed as six transition rate functions superposed with data points. The shaded areas added to the graphs suggest that a linear scaling of transition rate functions captures most of the underlying variance: <> = [0.85,1.15], <> = [0.7,1.3], <> = [0.8,1.3], <> = [0.7,1.35], <> = [0.85,1.15], <> = [0.8,1.5]. B is modified with permission from ref. 1.

To simplify matters, the following analysis ignores variations in equilibrium potentials and focuses on 10 parameters: membrane capacitance (); maximal sodium, potassium, and leak conductance (, , and , respectively); and the six transition rates underlying the opening and closing of “gates”—, , , , , and —as explained below.

We begin creating a straw man, considering marginal (i.e., independent) uniformly distributed variations for all 10 parameters. Values of parameters are expressed in terms of their scaling relative to the values chosen by Hodgkin and Huxley (1952). Thus, for instance, = 1.2 stands for = 144 mS/ (i.e., ×1.2 the value chosen by Hodgkin and Huxley; Fig. 1). Transition rates are similarly scaled by multiplication. Thus, for instance, the expression = 0.75 stands for . The shaded areas added to the graphs of Fig. 1 suggest that such a linear scaling of transition rate functions is justified, as it captures most of the underlying variance.

Hence, a realization of a Hodgkin–Huxley model is defined by a list of 10 scaling parameters: , , , , , , , , , and . Assuming independence of the 10 parameters, we randomly generated 10,000 such lists of scaling parameters with values range [0.75, 1.25] and numerically instantiated, each one of them, in a full Hodgkin–Huxley model. The resulting behaviors may be classified as nonexcitable, excitable (i.e., a membrane that generates one spike in response to a short above-threshold stimulus), and oscillatory (i.e., pace-making). Resting membrane potentials of the excitable and nonexcitable outcomes did not differ much over the 25% deviation from the original Hodgkin and Huxley parameters, being −64.5 1.4 mV and −66.2 1.5 mV (respectively).

To promote effective visualization of the 10-dimensional space, results are presented in the form of polar plots (Fig. 2): Given a vector of 10 scaling parameters, the value of each parameter is depicted along its own (angular) axis, and the entire vector is depicted as one line that connects the 10 scaling parameters. The standard Hodgkin–Huxley model would be a line passing through 1 for all scaling parameters. Three separate polar plots (Fig. 2, A–C) show that practically all three classes of excitability status (depicted by three different colors) are distributed throughout the 10-dimensional parameter space. Mean vectors in each of these cases are depicted by black dashed lines. For comparison, these three mean vectors and their corresponding standard deviations are plotted together in the polar plot of Fig. 2. The Euclidean distance between the mean vector of excitable and the mean vectors of the other two solutions (nonexcitable or oscillatory) is 0.15, similar to the SD of distances within each of them (Fig. 2 , Inset). In other words, assuming complete independence of the parameters within the 25% range of parametric variation, almost any randomly chosen vector of Hodgkin–Huxley parameters, regardless of its outcome (nonexcitable, excitable, or oscillatory), may be “pushed” to display any other excitability status by a minor manipulation of parameters.

Fig. 2.

Minor change in parameters significantly affects membrane excitability. Realizations (10,000) of a full Hodgkin–Huxley model; each realization is uniquely defined by a vector of 10 parameters, expressed in terms of their scaling relative to the values chosen by Hodgkin and Huxley. Responses are classified to three excitability statuses (different colors): excitable (2,225), nonexcitable (4,884), and oscillatory (2,891). Subsets of the results (200 for each excitability status) are presented in polar plots (A–C): Given a list of 10 scaling parameters, the value of each parameter is depicted along its own (angular) axis, and the entire vector is depicted as a line that connects the 10 scaling parameters. The standard Hodgkin–Huxley model would be a line passing through 1 for all scaling parameters. Mean vectors are depicted by dashed lines. The histograms in D depict Euclidean distance between vectors of scaled parameters within each of the three excitability classes.

Lower Dimension Hodgkin–Huxley Parameter Space.

We focus on the conditions for transition between excitable and nonexcitable statuses. Several schematic momentary current–voltage relations of excitable membranes, during an action potential, are plotted in Fig. 3. Grossly speaking, the lower curve depicts current–voltage relations sampled by voltage-clamp steps from deeply hyperpolarized holding potential. The upper curve depicts current–voltage relations sampled by voltage-clamp steps from a relatively depolarized holding potential. During an action potential, where membrane voltage is a dynamical free variable, current–voltage relations slowly shift between these two extremes due to an evolving voltage-dependent restoring force, mediated by the opening of potassium channels and inactivation of sodium channels. The slow change in restoring force gives rise to a current–voltage closed trajectory depicted in Fig. 3 (black continuous line).

Fig. 3.

Theory inspired dimensionality reduction. (A) Idealized momentary current–voltage relations at different ratios of available sodium and potassium conductances [see chapter 11 in Jack et al. (27)]. In different phases of the action potential, different momentary current–voltage relations determine the dynamics. The black continuous line depicts the resulting current–voltage trajectory during an action potential. (A, Inset) A magnified version of the area at threshold (indicated in the main image by a circle), about which the system is linearized. (B) Histograms of the three excitability statuses, constructed from the data of Fig. 2 (10,000 Hodgkin–Huxley realizations), for each of the scaling parameters. Note that all parameters are freely fluctuating, simultaneously, over 25%.

The general differential equation of the system is . As pointed out by Jack et al. (27) (chapter 11), linearization about the threshold potential (Fig. 3 , Inset) leads to an expression of in terms of the momentary conductance () at threshold and the time constant () for evolving restoring force. The condition for instability near threshold, one of the solutions to this equation (expressed in conductance units), is . Note that these lumped entities (, , and ) may naturally be classified into the above groups of structural and kinetic parameters: is obviously structural; likewise, the fast conductance () that depends on the relative contribution of maximal sodium conductance. In contrast, the time scale for introduction of restoring force () is a kinetic parameter because it depends on the actual transition rate functions governing the gating of sodium and potassium channels.

Inspired by the above and related mathematical reductions of excitability (24–26, 28), we turned to the data of Fig. 2 in search for these two dimensions, expressed in terms of Hodgkin–Huxley parameters. To this aim, we constructed count histograms of the three excitability statuses for each of the scaling parameters (Fig. 3). These histograms show that the most critical determinants of excitability status are the rates of opening and closure of sodium and potassium channels (, , , and ) and the maximal conductance of the membrane to the two ions (, ). (Similar results are obtained using the two “far apart” excitability statuses—nonexcitable and oscillatory—to extract principal components; see https://www.wolframcloud.com/objects/marom/SM.nb.) Distribution of excitability status is significantly less sensitive to transition rates involved in sodium conductance inactivation, as well as leak conductance. Membrane capacitance seems to have some effect. A possible interpretation of the histograms of Fig. 3 is that, at least as a first approximation, may be assumed to be proportional to a ratio of structural parameters S = , whereas may be assumed to be proportional to the ratio of kinetic parameters K = .

The murky enmeshment of solutions shown in Fig. 2 (a sample of which is replotted in Fig. 4, superposed) is significantly clarified when the data are arranged according to the values of S and K (Fig. 4): The three different excitability regimes are nicely clustered in three clouds. The oscillatory and nonexcitable phases are well separated in the S–K plane (Fig. 4 , Inset), whereas the borders separating the excitable phase from these other two are “soft” rather than sharp. The nonexcitable cloud in the upper left corner is due to excessive sodium conductance that stabilizes the membrane at a depolarized potential.

Fig. 4.

Hodgkin–Huxley model in a 2D S–K plane. (A) Subsets of 100 realizations from each of the three excitability statuses of Fig. 2 are plotted together. (B) Realizations (30,000) of a full Hodgkin–Huxley model, covering parametric variations over the entire range indicated by Hodgkin and Huxley (Fig. 1), classified (different colors) to three excitability statuses: excitable (4,660), not excitable (12,271), and oscillatory (13,069). Linear regression through the excitable status cloud (blue) is depicted by a line, the equation of which is S = 4.4K – 1.6. (B, Inset) Same plot with excitable points omitted.

Note that in Fig. 4, four parameters (two sodium inactivation rates, capacitance, and leak conductance) are not taken into account, even though they are allowed to freely fluctuate. And, yet, when examined in the 2D S–K space, the Hodgkin–Huxley model reveals order that is literally impossible to detect in the more explicit, higher dimensional representation of Fig. 4. The effectiveness of the dimensionality reduction is further supported in Fig. 5, where the outcomes of multiple realizations of three different S;K pairs are shown. Each of these three predefined S;K pairs (0.60;0.45, 0.50;0.50, and 0.40;0.55) was realized 30 times by adjusting and to the other four, randomly generated, Hodgkin–Huxley parameters (, , , and ). Five of 30 are shown for each value; the rest are comparable. Clearly, the values of the lumped S and K dimensions are better predictions of the outcome than the individual Hodgkin–Huxley parameters.

Fig. 5.

Response features are sensitive to the values of and . (A) Five different instantiations for each of three S;K pairs; stimulation amplitude 14 A. (B) The integral of voltage response emitted during a simulated trace is sensitive to the position within the S–K plain. The integral, calculated by summing the voltage values of all data points along the trace, relative to −65 mV, is presented in arbitrary units: B, Left, gradual change of at ; B, Right, gradual change of at . Point color depicts excitability class.

Several points deserve attention in relation to the numerically calculated S–K plane of Fig. 4. First, the borders between the three phases are fairly steep (note different ranges of S and K axes). The immediate implication of this steepness is that in its two extreme statuses (nonexcitable and oscillatory), the system is relatively immune to variations in maximal conductances. Stated differently, it is sufficient to use ionic channels that set the K dimension below 0.45 to obtain a pacemaker that is insensitive to fluctuations in density of channel proteins; the system maintains its pace-making nature over factor 3 in the value of the structural () dimension. The same can be said about nonexcitable membranes: setting the K dimension above 0.55 to obtain a nonexcitable system that is insensitive to fluctuations in density of channel proteins. Second, the K dimension is a rational function of first degree, a simple combination of concrete Hodgkin–Huxley kinetic parameters; as such it buffers the effect of changes in individual rates. This would also be the case for the S dimension when more than two voltage-dependent conductances are involved; viewed from this angle, excitability status becomes more robust when more conductances are involved. Third, moving within the S–K plane has an interpretable effect on the response shape (Fig. 5): The integral of voltage response during a simulated trace is sensitive to the position within the S–K plane, revealing intermediate excitability states within phases (e.g., lower amplitude or broader spikes). Naturally, difference between points in the nonexcitable phase is very small, if at all. Fourth is a seemingly technical point but of potential interest: Given scaled Hodgkin–Huxley parameters, one can calculate the resulting excitability status without resorting to simulation. This means that the nonlinearity of the model does not change the behavior actually predicted from a low-dimensional representation of the system (17). And, fifth, admittedly, our theoretically inspired choice of the rational functions that express S and K is one of many possible interpretations to relations between scaling parameters. To further justify this choice, we submitted the whole dataset to a linear support vector machine (SVM) algorithm. The results are presented in Fig. 6, where the test set and the probability distribution of each class are plotted as a function of S and K. The accuracy of classifying the outcome of a full Hodgkin–Huxley model, based on the values of S and K is 0.89, suggests that our theoretically-inspired reduction of the Hodgkin–Huxley model to an S–K plane is judicious.

Fig. 6.

Machine learning classification of excitability status. Classification of excitability statuses (data of Fig. 2) using a linear kernel SVM; 80% training set. Each surface (oscillatory, excitable, and nonexcitable; color coded) represents the probability (z axis) of a given combination of and to give rise to its corresponding excitability status. Thus, for instance, the probability (Prob.) of a point and to yield an oscillatory (depicted green) excitability status is practically 1. The colored points at the top of the image are the actual data points, similar to the form of presentation used in Fig. 4.

The roots of our analysis are the same as those motivating well-established simplifications of Hodgkin–Huxley, with linked to the cubic polynomial expression that provides fast positive feedback and to the recovery variable that introduces slow negative feedback. Another class of simplification involves reduction of the multistate gating diagrams of sodium and potassium channels while maintaining spiking features (e.g., ref. 29); this class of models may be considered as modifications to the dimension. One might contemplate the S–K representation as allowing those other simplified models to be seen coherently in a geometric representation, a point demanding systematic analysis that is beyond the scope of the present study.

We opened with a question about the possibility of maintaining reliable excitability in cells despite considerable variability in the values of key underlying parameters. Below, a possible solution is offered: when changes of excitability state are considered as dynamics within the S–K plane, slow inactivation—a ubiquitous activity-dependent feature of ionic channels—may act as a powerful local homeostatic control mechanism over a wide range of time scales.

Closed-Loop Control of Excitability in the S–K Plane: The Case of Sodium Conductance Slow Inactivation.

Being embedded in S–K, the seemingly complicated and parameter-sensitive system becomes tractable, enabling regulation by an activity-dependent rule acting on one physiological entity. A most straightforward regulation rule would involve inverse relations between electrical activity (say, integral of membrane potential depolarizations) and the effective or actual value of the structural dimension (S). Many physiological processes that modulate membrane ion channels may realize such adaptation, covering a wide range of spatial and temporal scales (30). For instance: (i) Slow inactivation of sodium conductance, which is a local modulatory mechanism that operates over seconds to many minute time scales (14, 31–37); (ii) calcium-dependent activation of potassium conductance, a local and relatively fast time scale mechanism (38–40); or (iii) regulation of sodium and/or potassium channel protein expression, an arguably global but definitely slow time scale mechanism (22, 41). Each of these in itself naturally constrains the system to hover about the excitable phase by pushing the value of S downward when above the diagonal, or upward when below. Other regulatory mechanisms may be envisioned, implementing (for instance) temporal integration of subthreshold activity by slow inactivation of potassium channels (42, 43) or maintenance of pace-making activity by regulation of conductance (44).

Of the above-mentioned spectrum of physiological modulatory mechanisms, slow inactivation of sodium channels is especially interesting. While acknowledged from the early days of membrane electrophysiology (reviewed in ref. 35), the concept of slow inactivation of sodium channels as a means to maintain excitability status amidst parametric variations has remained relatively neglected. What makes slow inactivation a powerful regulatory mechanism is its impacts on the effective value of , covering a range of time scales (33–36, 45, 46). Furthermore, inactivation is “local” in the sense that it does not require central control; it occurs automatically as a consequence of activity. Thus, one might picture it acting as a distributed normalizing force in extended excitable tissues (e.g., long axons or electrically coupled excitable cells).

To demonstrate the potential impacts of slow inactivation on dynamics within the S–K plane, we focus on channel gating beyond the time scale of a single action potential. Slow inactivation is represented as a macroscopic system, where channels move—in an activity-dependent manner—between two states: available and not available (Fig. 7, Left), depicted . The first (A) is the set of states that includes, besides the open state itself, all of the states from which the channel may arrive to the open state within the time scale of a single action potential—that is, the closed states and the very first inactive states that are treated in standard Hodgkin and Huxley formalism for the action potential generation. In other words, a channel in A is available for conducting ions within the time scale of a single action potential. In contrast, is a large, interconnected pool of slow inactive states (depicted Ij’s in the scheme of Fig. 7) from which transition to the open state within the time scale of an action potential is impossible.

Fig. 7.

Slow inactivation of sodium channels stabilizes excitability in S–K plane. (Left) A schematic representation of sodium channel states, with many slow inactivation states. (Right) Demonstration of maintenance of excitability in S–K plane given parametric variation, controlled by activity-dependent transitions of sodium channels between available and not-available sets of states. The simulation describes a 200,000-step random walk process, beginning at 0.5;0.5 (). The gray trajectory (squiggles-like) depicts a random walk where and . The blue line depicts a walk where , and . The slope corresponds to the fitted function of Fig. 4.

Recent structural analyses suggest that the large space of slow inactive states might reflect the many distorted versions of the functional protein under conditions where the organization, otherwise enforced by hyperpolarized membrane potential, is compromised upon extensive depolarizing activity (31). Theory and experiments (14, 30, 34, 45, 47–49) show that, in such a scheme, the multiplicity of slow inactive states entails a power-law scaling of recovery from to as a function of time spent in . This implies a potential to become dormant in an activity-dependent manner for a duration ranging from tens of milliseconds to many minutes and possibly hours. Thus, unlike standard Hodgkin–Huxley gates, the rate of recovery from slow inactivation does not have a uniquely defined characteristic time scale. Rather, the time scale is determined by the distribution of channels in the space of inactive states, which, in turn, is dictated by the history of activation. The kinetics of may be qualitatively described by an adaptive rate model (30, 50–53), a logistic-like function of the form: , where is a function of some general activity measure , and is a monotonically increasing function of the system state . The model gives rise to a wide range of time scales of recovery from inactivation and assures, when simplified, a nonzero stable point at , on the edge between excitable and nonexcitable (50).

Mapping the above picture to the terms used in the present work, it is instructive to think of as (i.e., a scaling parameter of maximal sodium conductance). In the original Hodgkin and Huxley formalism, is a structural constant that sets limits on the instantaneous (at the scale of milliseconds) input–output relations of the membrane. But when long-term effects are sought, might be treated as a dynamic variable that modulates residence in a reservoir of slow-inactivation configurations, pulling channels away from the system as a function of activity. Note that where is constant and where = , the adaptive rate model qualitatively captures the dynamics of /( + ). Indeed, Fig. 7, Right demonstrates that application of a simple adaptive rate model , where is substituted by , reveals the potential of sodium conductance slow inactivation to maintain excitability amid parametric variations. Slow inactivation restrains the system to a diagonal in the S–K plane: The blue line depicts a case where the kinetic dimension () walks randomly while follows an adaptive rate formalism; a gray line that connects between points (looks like squiggles) depicts the path of excitability status in a control condition, where both the kinetic dimension () and the structural dimension () walk randomly.

Concluding Remarks.

In an era marked by capacity to collect data at ever-increasing resolution, it is important to identify proper scales in analyses of cellular phenomena, scales that matter to the system (54), scales where the phenomenon of interest is low-dimensional, regulatable by simple physiological processes and explainable in simple physiological terms. What makes membrane excitability—a fundamental physiological phenomenon—particularly attractive to study in this context is the existence of a relatively sound theory in general, its application in the Hodgkin–Huxley formalism in particular, and its amenability to experimental manipulations at both microscopic (channel protein) and macroscopic (membrane potential) levels.

Practically all of the homeostatic-based models of excitability regulation used in the past have kept kinetics constant and only looked at channel densities (11, 19, 21). The present study suggests that when the problem is examined in a lower dimension, a simple control rule that relies on slow inactivation—a ubiquitous protein intrinsic process—can deal with fluctuations in both structural and kinetic parameters. This homeostatic mechanism is local, independent of protein synthesis, and operates over a wide range of time scales (milliseconds to many minutes).

We speculate that activity dependence of protein kinetics at relatively slow time scales, entailed by multiplicity of protein states, is a general “automatic” and local means for stabilization of cellular function.

Materials and Methods

All of the simulations and analyses were performed within the Mathematica (Wolfram Research, Inc.) environment. The Mathematica notebook is available at https://www.wolframcloud.com/objects/marom/SM.nb. Data of Figs. 2–5 were generated using Hodgkin–Huxley equations, scaled to 6 °C, as they appear in the original manuscript (1952). The duration of each simulation epoch was 90 ms, including an initial 50-ms relaxation phase. Stimulus (7 A, 1 ms) was delivered 70 ms into the epoch. A sorting algorithm for excitability status (excitable, nonexcitable, or oscillatory) was constructed, which is based on the time and the number of spikes following the relaxation phase; spike detection threshold was set to 50 mV above resting potential. The sorting algorithm was validated by eye-inspection of multiple sets of 300 sorted epochs. To generate Fig. 6, we implemented an SVM algorithm with linear kernel within Mathematica, using an 80% training set.