Integral Control from a Simple Model of Ion Channel Expression

Experiments have shown that the processes responsible for regulating intrinsic neuronal properties are slow relative to fluctuations in electrical activity (Desai et al., 1999; O’Leary et al., 2010; Thoby-Brisson and Simmers 2000). These processes behave as a feedback control mechanism that monitors average activity and adjusts membrane conductances to achieve some kind of target activity. An important readout signal appears to be intracellular [Ca²⁺], which correlates with electrical activity due to voltage-dependent Ca²⁺channels and buffering mechanisms that average out fluctuations in time and space (Berridge 1998; Wheeler et al., 2012). Moreover, long-term changes in [Ca²⁺] are found to regulate many ion channel types (Barish 1998; Mermelstein et al., 2000; O’Leary et al., 2010; Turrigiano et al., 1994; Wheeler et al., 2012).

Hence we developed a model of activity-dependent conductance regulation using intracellular [Ca²⁺] as a feedback control signal. There are many ways to implement feedback control of membrane conductances (Günay and Prinz 2010; LeMasson et al., 1993; Liu et al., 1998; Olypher and Prinz 2010; Stemmler and Koch 1999). We wanted to focus on a rule that captures essential biological principles and has experimentally testable properties. Ion channels are proteins, and their expression depends on the level of channel mRNA in the cell. A simple way of capturing this leads to a canonical model of regulation (also known as the “central dogma” of molecular biology):

Here, m is the concentration of mRNA for channel protein g and α_x and β_x are synthesis and degradation rates; dots denote time derivatives. The biochemical scheme underlying this model is shown in Figure 1A. In spite of its simplicity, this model has proven useful for understanding gene expression dynamics in systems biology (Alon 2007). Neurons possess a rich repertoire of other regulatory mechanisms, including alternative splicing, alternative promoter usage, RNA interference, regulated protein trafficking, and posttranslational modifications to channel proteins. Therefore, the simplified scheme we use is a first approximation that can be refined to take into account more intricate aspects of regulation when and where sufficient experimental data are available.

Where does activity dependence enter this model? mRNA expression rates depend on transcription factor activation. Many important transcription factors such as CREB are known to be Ca²⁺ dependent or dependent on other Ca²⁺-sensing enzymes (Finkbeiner and Greenberg 1998; Mermelstein et al., 2000; Mihalas et al., 2013; Wheeler et al., 2012). Furthermore, transcriptional changes in ion channel genes occur in response to activity perturbations (Kim et al., 2010) and may underlie homeostatic regulation of network activity (Thoby-Brisson and Simmers 2000). We therefore assume that mRNA production depends on some Ca²⁺-sensitive enzyme, or enzyme complex, T, whose production rate is Ca²⁺ dependent and whose rate of degradation is saturated (Drengstig et al., 2008). Incorporating this into the model, we have

We do not know in general how the forward rate, α_T depends on [Ca²⁺]. The assumption that the degradation rate of T is saturated means the equilibrium of the system occurs at a unique value of [Ca²⁺]. Thus average [Ca²⁺] will be maintained at a specific “target” value, Ca_tgt, given by solving the steady-state, = 0. We do not know in general how the forward rate, α_T depends on [Ca²⁺]. If this rate is determined by a single reaction involving Ca²⁺ binding, then it will typically have a monotonic dependence in the form of a Hill equation (Supplemental Experimental Procedures). Thus, for simplicity, we assume a linear approximation, α_T([Ca^{2 +}]) = α_T·[Ca²⁺]. In this case the target [Ca²⁺] is simply the ratio of two rate constants: Ca_tgt = β_T/α_T.

We can now show how mRNA and conductances are regulated to keep the system at Ca_tgt. Inspection of Equation 2 reveals that T explicitly integrates the difference between [Ca²⁺] and Ca_tgt overtime:

This integrated Ca²⁺ “error” signal is then fed into the synthesis term of the channel mRNA (α_mT, Equation 2). Similarly, m directly controls the expression rate of g via the term α_gm. Finally, g controls the membrane potential and Ca²⁺ dynamics. The scheme in Equation 2 therefore constitutes a feedback loop that maintains average [Ca²⁺] by continually modifying the expression rates of channels in the membrane. This is illustrated in Figure 1B, where the shaded area shows the accumulated error signal over time. If [Ca²⁺] is different from the target, error will accumulate and drive changes in the expression of mRNA and membrane conductances until the system reaches equilibrium at [Ca²⁺] = Ca_tgt.

So far, we assumed that a global regulator (T) controls downstream precursors of membrane conductances. What if these conductances are controlled by independent pathways that have the same integrating mechanism acting on the same error signal? An immediate problem arises if the set points for each controller are not tuned so that they all agree precisely. In Figure 2A, an inward and outward conductance are under the control of a single Ca²⁺-integrating regulator, T₁, with target t₁. Now suppose (Figure 2B) that each conductance is controlled by separate regulators, T1, T2, with different targets, t₁ ≠ t₂. Two possibilities exist: either one target will become satisfied, in which case error will accumulate without bound in the other controller, or, as will be the case more generally, neither target will be satisfied and both controllers will accumulate error without bound. This is shown in Figure 2B, where the two conductances upregulate without bound. In control theory, this accumulation of error is known as “windup.” In biological terms, windup would result in unbounded (eventually saturating) production of mRNAs and channels and loss of regulatory control. This could be avoided in this scenario if the reaction rates determining the independent targets are precisely matched; however, precise tuning seems unlikely in biological systems. We therefore conclude that for this model to work across a set of conductances, a single master regulator pathway is preferable. In more complex schemes with several distinct regulatory signals, it is possible to have separate targets for each signal (Liu et al., 1998). However, windup can still occur if these multiple signals cannot be satisfied simultaneously.

Specifying and Maintaining Cell Types with Multiple Regulated Conductances

The above analysis shows how activity in neurons can be maintained using a simple model of gene regulation. Overall, the regulation scheme can be written in simplified form for a neuron with multiple conductances, g_i:

In this simplified form, the equations for and m are lumped together (Experimental Procedures), concentrations are scaled, and reaction rates are replaced by single time constants. τ_g represents the characteristic time constant of channel expression and τ_i represents the coupling of channel gene expression to [Ca²⁺]. We use this simplified system in what follows.

The model achieves target [Ca²⁺], but what combination of conductances will the neuron express at this target? There are typically many ways to reach the same average activity using different maximal conductances (Bhalla and Bower 1993; Golowasch et al., 2002; Olypher and Calabrese 2007; Prinz et al., 2003; Sobie 2009; Swensen and Bean 2005; Taylor et al., 2006, 2009). Recent mathematical work has made the relationships between conductances and excitability clearer and more precise and can be understood in terms of ratios of conductances that act on different timescales (Drion et al., 2012; Franci et al., 2012; Franci et al., 2013). We see next that the integral control rule produces “nice” conductance distributions with constant ratios that are qualitatively similar to those observed biologically and can thus generate “cell types.”

Previous work (O’Leary et al., 2013) showed how regulation time constants determine correlations in conductance expression at steady state. What is the relation between conductance ratios and the regulation time constants τ_i in the present model? From the simplified scheme (Equation 3), each m_i converges to a value that depends on the time integral of average [Ca²⁺], scaled by the inverse expression time constant τ_i. Thus, we can estimate steady-state g_i for positive time constants and small initial conductances:

When taking ratios, the integrals cancel, so that:

In summary, different ratios of the τ_i ’s specify correlations between each conductance. Correlations in conductance expression, in turn, promote defined electrophysiological characteristics, because ratios of different kinds of voltage-gated conductances largely determine single neuron dynamics (Drion et al., 2012; Franci et al., 2013; Hudson and Prinz 2010).

Figure 3 shows an example of a complex model neuron with seven voltage-dependent conductances, all regulated by the integral control rule. The time evolution of the membrane conductances for multiple, randomly initialized runs of the model is shown (Figure 3A) with the membrane potential of an example neuron at four different time points. Each model starts in a nonexcitable nascent state (left-most trace, Figure 3A), which is ensured by choosing the random initial distribution of conductances to be small relative to the leak conductance that is itself randomly distributed (Experimental Procedures). The steady-state behavior of the model neuron is stereotyped and develops a rhythmic bursting activity (right-most trace, Figure 3A). In spite of varying initial conditions, the models at steady state all have similar membrane potential activity, as can be seen in the example traces of six model neurons that developed from different initial conditions (Figure 3B).

Figure 3C (lower panel) shows steady-state conductance distributions and intrinsic properties of 20 independent runs of the model. The conductance densities vary several-fold over the population, but all of the neurons have a similar bursting phenotype. Both the conductances and the firing properties show clear pairwise correlations that are reminiscent of experimental data in identified crustacean as well as mammalian neurons (Amendola et al., 2012; Liss et al., 2001; Schulz et al., 2007; Tobin et al., 2009). In contrast to previous modeling work that used a less biologically realistic regulation rule (Liu et al., 1998; O’Leary et al., 2013), the conductance correlations are approximately linear. Second, the use of a single activity sensor in the present model differs from this previous work, which used three sensors with different timescales to promote bursting behavior. While we do not rule out the possibility that more than one activity sensor is used biologically, the current model shows that this is not necessary.

Generating Cell Types

The simple relationship between regulation time constants and conductance ratios in the model means we can determine arbitrary correlations between conductances and thus construct self-regulating cells with specific intrinsic properties using only a single activity sensor. These intrinsic properties can encompass any excitability type, provided we have a sufficiently rich set of available conductances (i.e., a sufficiently rich “genome” in the model). Figure 4 shows five distinct neuron types that are achieved using the same underlying model and random initial conditions but with appropriately chosen sets of regulation time constants and [Ca²⁺] targets. We can thus specify cells that establish and maintain specific input-output relations, as quantified by the Type I/II FI curves in the first two examples. Type I excitability is characterized by the existence of arbitrarily low firing frequencies at spiking threshold (Figure 4A, example 1); in contrast, Type II excitability does not support firing below a fixed nonzero rate (example 2, Figure 4A) (Rinzel and Ermentrout, 1989). Notably, in spite of having different firing properties, both of these models have the same [Ca²⁺] target. Similarly, we can specify cells that respond reliably to input, as exemplified by an excitable rebound bursting cell that generates action potentials coupled to slow membrane potential oscillations in response to both depolarizing and hyperpolarizing current (example 3, Figure 4A). Finally, we can specify cells that are active in either a tonic spiking mode or bursting mode (examples 4 and 5, Figure 4A). Each of these cell types has a unique correlation structure in its steady-state conductance distribution (Figure 4B) following multiple runs from the same initial conditions. Furthermore, the straight line calculated from the ratio of expression time constants in Equation 4 predicts the pairwise conductance distributions in each case (Figure 4B).

The conductance distributions produced by the model may explain cell-specific linear correlations that are found biologically. Figure 4C reproduces data from Schulz et al. (2007) in which the expression levels of multiple ion channel genes were measured in single stomatogastric ganglion (STG) neurons using quantitative PCR. This revealed linear correlations in the expression that are specific the cell type (Figure 4C). The model predicts that the slopes of the correlations in the data should equal the (time-averaged) expression rates of the respective mRNAs in each plot.

There are two important biological assumptions in this instantiation of the model. First, the leak conductance, which can be thought of as an aggregate of multiple conductances, is static but varies between cells. This serves as a model of conductances that are not regulated by the integral control rule and that may vary across a population of neurons. Second, we have assumed that expression rates are fixed. Biologically, this corresponds to steady values in average promoter activity, binding affinities of signaling enzymes, translation rates, protein trafficking, and degradation rates. While this may be a reasonable assumption at any given stage in nervous system development, it is entirely plausible that these relationships change over time. Neurons may thus cycle through several physiological “types” as they develop, and this process will be highly specific to the species and brain area in question. We have not explicitly attempted to model these transient stages as our goal is more general; they can however be incorporated by switching the rates in the regulatory rule–an idea we examine in the final part of this work.

Expression Rates versus Activity Targets as Determinants of Electrophysiological Properties

We have shown that different sets of regulation rates/time constants determine cell types and that distinct cell types can have the same [Ca²⁺] target in principle. What happens when the [Ca²⁺] target is scaled within a cell type? While the regulation time constants determine the direction in which the cell moves in conductance space, the target determines how far it travels along a trajectory before reaching equilibrium. Thus, targets can determine the location of the conductance distribution as well as scaling activity. Figure 5 shows the steady-state activity of three example neuron types as target [Ca²⁺] is scaled. Below the traces are plots showing the regulation time constants for each cell type. Typically, as the target is raised, spiking activity elevates because this corresponds to greater average Ca²⁺ influx, as can be seen in the first two examples. In some cases, moving the target can also cause a qualitative change in activity as seen in the third example, which transitions from bursting to spiking as [Ca²⁺] target is increased. Thus, the combined contributions of ion channel expression dynamics can be dissociated from activity set points in neurons, but both have a role in determining physiological properties.

Open in a separate window

Figure 5

Changing Targets within Cell Types

Each column shows 500 ms segments of steady-state membrane potential activity in a different self-regulating model at steady state with the [Ca²⁺] target (= 4 µM) scaled. The regulation time constants for each conductance are shown below, normalized to τ_g.

A Self-Assembling Motor Circuit

If we can reliably specify cell types using this model, it should be possible to construct a self-assembling, homeostatically regulated network whose activity depends on specific properties of the component cells (Golowasch et al., 1999b). The pyloric central pattern generating network of the crustacean STG consists of three distinct cell type modules: a pacemaker complex and two follower cell types that fire in successive phases. Activity in this network consists of a triphasic pattern of firing starting with the AB/PD complex, followed by the LP cell and then PY (Marder and Bucher 2007). The synaptic connectivity is known (Figure 6A) and consists of slow inhibitory cholinergic synapses as well as fast inhibitory glutamatergic synapses (Marder and Eisen 1984). We reasoned that by finding steady-state conductances and synaptic strengths that produced a triphasic pattern, we could then find regulation time constants for intrinsic conductances and synapses that would dynamically specify and maintain a characteristic network activity pattern.

Open in a separate window

Figure 6

A Self-Assembling, Self-Regulating Central Pattern Generating Network

(A) Connectivity diagram of the model CPG, based on the synaptic connectivity of the pyloric network in the crustacean STG (PD/AB, pyloric dilator/ anterior burster; LP, lateral pyloric cell; PY, pyloric cell). The PD/AB pacemaker kernel is modeled as a single cell. All synapses are inhibitory and graded; glutamate (Glu) synapses are instantaneous, acetylcholine (ACh) synapses are slow (activation time constant = 50 ms).

(B) (Top) Example membrane potential traces for random initial conductances. (Second from top) Example steady-state behavior of the model. The triphasic order (PD, LP, PY) is highlighted with shaded boxes. (Third from top) Perturbation of network activity by addition of hyperpolarizing (reversal potential = −80 mV) conductance to PD. (Bottom) steady-state recovery of the network with hyperpolarizing conductance still present. All traces = 1 s.

(C) Example time evolution of intrinsic and synaptic conductances in a self-regulating pyloric network model for a single run. Onset of the PD/AB perturbation is indicated by the vertical line. Insets show detail of the conductance dynamics on a linear timescale.

After randomly searching conductance values to find combinations that produced bursting pacemaker cell types, we hand tuned an unregulated network to produce a qualitatively realistic triphasic rhythm. We then converted maximal conductances and synaptic conductances in the hand-tuned network to expression time constants using Equation 4. We next searched around this initial set of time constants using a log-normal distribution to find those that reliably produced triphasic networks from random initial conditions (Experimental Procedures).

The network always starts in a nonfunctional state (Figure 6B, top). The membrane potential activity of the cells after the network has reached steady state (Figure 6B, second from top) shows a regular triphasic rhythm. Furthermore, the network activity is robust to perturbations at steady state, as exemplified by recovery from the addition of a hyperpolarizing leak conductance (0.02 µS, E_rev = −80 mV) that silences the PD/AB pacemaker cell (Figure 6B, third and fourth panels from top). Over multiple runs (n = 507) of this model, 99.6% produced stable triphasic rhythms. Of these, 93.5% recovered after the perturbation (which abolished rhythms in 99.2% of networks). Figure 6C shows the evolution of intrinsic and synaptic conductances in the example network of Figure 6B. Notably, the synaptic and intrinsic conductances in all cells respond to the perturbation in the PD/AB cell.

Activity-Dependent Regulation Can Be Compensatory or Pathological

Under what conditions does activity-dependent regulation compensate for mutation or pharmacological blockade? The conductances in the model neurons, as in biological neurons, overlap in some of their properties. Thus, if certain conductances are lost, others can be upregulated or downregulated to compensate. Figure 7A shows the steady-state behavior of a self-regulating bursting pacemaker neuron. Upon deletion of the I_h conductance, the models become silent, leading to a decrease in average [Ca²⁺]. Following deletion, the integral control rule restores bursting activity by altering conductance expression to achieve target [Ca²⁺]. Similar outcomes are possible when the deletion has variable effects owing to variability in the cells produced by the model. Figure 7B shows the effect of deleting a slow Ca²⁺ conductance in two different examples (that have converged to different maximal conductances) of the same cell type. In one example, the model increased in frequency; in the other, the model became silent. Again, average [Ca²⁺] encodes this increase or decrease in activity, and the resulting conductance regulation restores bursting.

Open in a separate window

Figure 7

Outcome of Homeostatic Compensation after Channel Deletion Depends on Cell and Channel Type

Membrane potential activity for a self-regulating bursting ([A]–[C]) and tonic ([D] and [E]) pacemaker models in which specific conductances are deleted. The first column (“wild-type”) shows model behavior at steady state with all conductances present. Acute deletion of the indicated conductance produces the behavior shown in the middle column (“acute KO”). Following conductance deletion, each model is allowed to reach steady state (third column, “compensated KO”).

Can compensation lead to loss of function? The model assumes that neurons sense a gross physiological variable, [Ca²⁺], which cannot always distinguish specific activity patterns. Previous work identified this as a potential problem for regulation (Liu et al., 1998), but as we have shown here, it is nonetheless possible to use [Ca²⁺]to generate and maintain specific electrical properties using differential ion channel expression rates. However, this model will fail to preserve neuronal properties if the relationship between electrophysiological properties and [Ca²⁺] activity changes drastically. Such a change occurs in the example neuron in Figure 7C. Deletion of the transient Ca²⁺ conductance, g_CaT, silences the neuron, but following compensation to target [Ca²⁺], the neuron no longer bursts and instead fires tonically. In this case, the deletion of g_CaT resulted in changes in Ca²⁺ dynamics so that target [Ca²⁺] occurs for a fundamentally different pattern of membrane potential activity.

In a different tonic spiking model, deletion of one of the two Ca²⁺ currents has distinct effects. In Figure 7D, deletion of the slow Ca²⁺ conductance, g_CaS, slightly alters the spiking frequency, and this is compensated by regulation. However, deletion of g_CaT (Figure 7E) results in faster spiking, and compensation to the [Ca²⁺] target instead renders the cell silent.

Model Interpretation

A technical message of this work is that a canonical model of channel expression can be interpreted as a well-known control law: the integral controller. Integral control has been suggested as a mechanism of neuronal homeostasis based on the available molecular machinery for integrating Ca²⁺ signals in neurons (Davis 2006; O’Leary and Wyllie, 2011). We showed in this work how activity-dependent transcription can be an instantiation of integral control. The essential component of integral control is a variable whose rate of change depends on error. In the model presented here, error is deviation in [Ca²⁺] from a specific value, resulting in a change in the equilibrium of a putative regulator enzyme. The rate of change of ion channel mRNA is proportional to this error; consequently, ion channel mRNA concentration can be interpreted as the “accumulated error signal.” Biologically, the regulator enzyme could be a Ca²⁺-dependent transcription factor complex, or a Ca²⁺-binding enzyme upstream of a set of transcription factors. The biological counterpart of the [Ca²⁺] signal we consider is therefore a somatic or nuclear [Ca²⁺].

The form of the model placed a strong constraint on its implementation. If multiple, parallel integral control pathways using the same error signal exist within a cell, the targets for each pathway need to agree, otherwise the continual (and deleterious) accumulation of the molecules that encode error (such as mRNAs) will occur. While in principle it is possible that multiple parallel controllers are tuned so that their set points are equal, in biological reality, slight deviations are unavoidable. Thus for this model to work as a means of jointly regulating conductances in a neuron, a “master regulator” may be required. However, this does not rule out the possibility that other controllers using different error signals may coexist, provided there are conductance combinations that can simultaneously satisfy all controllers. For example, some currents could have targets specified by [Ca²⁺] transients or concentrations of other biological molecules provided these signals are sufficiently independent.

Integral control exists as a regulatory mechanism in simple organisms such as bacteria, where it permits sensitivity to environmental chemical cues and robust chemotaxis (Alon et al., 1999; Yi et al., 2000). It is thus a plausible and testable hypothesis that neurons have developed integral control pathways to regulate membrane conductances. Integral control implies perfect compensation in the control variable (average [Ca²⁺] in our case). Conversely, in systems that can be locally linearly approximated, perfect compensation implies integral control (Yi et al., 2000). Therefore, an experimental test of one assumption of this model is whether a relevant physiological variable such as average [Ca²⁺] is perfectly compensated over a range of perturbations. It is important that the perturbations do not exceed the capacity of the system to compensate, so a carefully controlled range of perturbations may be required along with precise monitoring of [Ca²⁺] to do this test.

The biochemical framework also allows a straightforward interpretation of an “activity target.” The nature of this target has been a source of speculation and even controversy since homeostatic regulation was first proposed (Maffei and Fontanini 2009; Marder and Prinz 2002).Weshowed that target [Ca²⁺] can be encoded by the rates of the underlying molecular mechanisms. Because these rates ultimately depend on chemical properties of enzymes, such as substrate binding affinity, the target can be reliably defined in a given cell or cell type.

How literally should one interpret this model? The mechanisms involved in regulating neuronal conductances are the focus of ongoing research and have many intricate components that we have omitted. Transcriptional control is involved in ion channel regulation (Weston and Baines 2007), and transcript editing, alternative splicing, and RNA interference can occur at the early stages of the process (Lin et al., 2012; Seeburg and Hartner 2003; Wang 2013). Similarly, at the stage when functional channels are expressed in the plasma membrane, phosphorylation and auxiliary subunit interactions can alter the biophysical properties of channels (Lipscombe et al., 2013). We did not attempt to model the effects of all such processes; instead we focus on the major events underlying channel expression that are encapsulated in the canonical model of gene expression: channel genes are activated, channel mRNA is transcribed, and channel protein is produced from mRNA. This simplification can be thought of as averaging out the contribution of more intricate processes, or as forming a backbone onto which the additional processes can be added. The task of refining the model will not be trivial; while we would expect the canonical model to hold across species and cell types, it would be surprising if more detailed models generalize without incorporating data that are specific to each experimental preparation.

Model Predictions

This work makes three general predictions. First, it predicts linear correlated variability in ion channel expression and that the slopes of the pairwise correlations between two ion channel expression measures should correspond to the ratio of their expression rates. For example, if one were to measure the average mRNA expression rates of two ion channels that are known to show a positive linear correlation in single-cell quantitative PCR measurements, then the ratio of the expression rates should equal the measured correlation slope (for example, the K⁺ channel genes shaw and shab in LP cells of the crab STG–see Figure 4B and 4C). Measuring mRNA expression dynamics is challenging and has not, to our knowledge, been performed in single neurons, although tools that may permit such measurements are being developed. On the other hand, single-cell quantification of steady-state ion channel gene expression does indeed show cell-type-specific correlations that are close to linear (Liss et al., 2001; Schulz et al., 2006, 2007; Temporal et al., 2012; Tobin et al., 2009).

Second, the model predicts that neurons do not necessarily perfectly compensate their electrical properties when perturbed, or when an ion channel type is knocked out, even if average [Ca²⁺] (or the relevant activity signal) is perfectly compensated over long timescales. This is illustrated in Figure 7, where we see that average [Ca²⁺] is always compensated, while the physiological behavior of the neuron can be compensated, partially compensated, or can even show pathological changes in behavior caused by the regulatory mechanism. In nonpathological cases, the model works because the regulation signal, [Ca²⁺], distinguishes different regions of conductance space. However, this sensing mechanism is dependent on certain combinations of conductances being present together. For example, if the ratio of delayed-rectifier K⁺ to fast Na⁺ conductances is within a certain (possibly large) range, then (ignoring other conductances and assuming a source of Ca²⁺ influx) low average [Ca²⁺] will correspond to silent cells, while high average [Ca²⁺] will only be achievable if the cell is firing tonically. Removal of one or more conductances can drastically alter the relationship between firing properties and [Ca²⁺], causing aberrant compensation. On the other hand, if the conductances overlap in their properties with other conductances, then removal may have only a subtle effect, or a substantial acute effect that can be compensated by the regulation mechanism.

This disconnect between nominally homeostatic behavior in one variable and nonhomeostatic behavior in the larger system has been suggested previously (O’Leary and Wyllie 2011) and illustrates the need for a careful definition of what homeo-stasis means. The safest definition is that homeostasis is an emergent phenomenon and occurs because the components in biological systems (such as ion channels) are often regulated using feedback. In general, the feedback signals can be a subset of those available and may act on a subset of the systems parameters (i.e., the system may be underactuated). Thus, neurons can exhibit firing rate set points (Hengen et al., 2013) or even maintenance of a coordinated motor pattern (Figure 6), but this does not necessarily mean the system directly measures and maintains these specific properties. This point is perhaps underappreciated but important because it is difficult to assess experimentally which are the controlled features of a homeostatic process.

The third broad prediction of this work is that changes in the regulatory rule itself may be part of nervous system development. The sets of regulation rates that define cell types in the model bring each cell to a steady state from random initial conditions with low conductance densities. This is a reasonable model of the early stage of differentiation from a nonneuronal cell with a nonexcitable (or weakly excitable) membrane (Moody and Bosma 2005; Spitzer et al., 2002). To reach mature levels of conductance expression, all conductances need to increase initially. This coordinated increase is inherent in the model of cell types (Figure 4) and predicts that the appearance of each ion channel type above some detection threshold will show a cell-type-specific ordering, in agreement with experimental observations in developing nervous systems (Baccaglini and Spitzer 1977; Moody and Bosma 2005; Spitzer 1991).

However, other experiments in nominally mature systems have shown that conductances can change their expression in opposite directions in response to perturbations in activity (Desai et al., 1999; O’Leary et al., 2010). Furthermore, while one rule may be sufficient for establishing a broad phenotype, changes to the rule could fine-tune conductances so the cell can preserve specific properties more effectively. We explored this idea speculatively in Figure 8, where we showed that tighter control of specific properties entails a switch in the regulation rates. Moving to a biological interpretation, this idea incorporates the observation that molecular switching events alter the expression rates of different genes (including ion channels) early in development and that some developmental changes have strictly sequential critical periods.

The authors thank Michael Wright and Paul Miller for comments. Funding was provided by NIH 1P01NS079419 and the Charles A. King Trust.