Simultaneously regulating many properties requires that neurons adjust diverse ion channels

Different ion channel combinations can yield equivalent neuronal excitability. This degeneracy facilitates robust regulation of excitability by enabling a disruptive change in one channel to be offset by compensatory changes in other channels. But various aspects of excitability plus other cellular properties require regulation. Coordinately regulating multiple properties via control of overlapping sets of ion channels is no easy task. Here we demonstrate that of the many ion channel combinations producing the target value for one property (the single-output solution set), few combinations produce the target value for other properties. Combinations producing the target value for two or more properties (the multi-output solution set) correspond to the intersection between single-output solution sets. We demonstrate that for multi-output solution sets to be degenerate, the number of controllable ion channels (nin) must exceed the number of regulated properties (nout). We further demonstrate how the dimensionality of solution space shapes ion channel correlations that emerge during homeostatic regulation and how homeostatic regulation may fail given the heightened challenge of simultaneously regulating many properties.


INTRODUCTION
Neurons have been shown to adjust their membrane excitability and synaptic weights in order to maintain their average firing rate near a target value, or set point [1][2][3][4][5][6] . This homeostatic regulation is achieved through feedback control of diverse ion channels with intracellular calcium serving as an important feedback signal 1,2,4,7 . Computational models have successfully employed feedback signaling to adjust ion channel densities [8][9][10][11][12] . Control theory provides a valuable framework to conceptualize how this regulation occurs 13 but most of the mechanistic details remain unclear 14,15 and are not straightforward; for instance, different perturbations can trigger similar changes in excitability via different signaling pathways affecting different ion channels 16 , or via different combinations of excitability changes and synaptic scaling 17 .
Just as each property depends on multiple ion channels, each ion channel can affect multiple properties 38 , which confounds the independent regulation of each property. Yet specific aspects of excitability, not just average firing rate, may need adjusting in order to modulate dynamic range and maintain good neural coding as synaptic weights or other conditions change 51 . And like other cells, neurons must also regulate their energy, osmolarity, pH, protein levels, etc. [52][53][54] . These other cellular properties affect and are affected by neuronal activity, and their regulation likely involves feedback signals and proteins (including ion channels, pumps and transporters) affecting excitability; for example, manipulating metabolic processes can change the firing rate set point 55 .
In short, properties are not regulated in isolation from one another, which prompts the question: How does a neuron manage to coordinate the regulation of its many properties?
Though diverse ion channel combinations can produce the target value for a certain property (e.g. firing rate), we hypothesized that only a subset of those ion channel combinations would also produce the target value for a second property (e.g. energy efficiency). In other words, the set of ion channel combinations producing one target output (i.e. the single-output solution set, or SOSS) is potentially quite large thanks to ion channel degeneracy, and though each property might have a large SOSS, the set of ion channel combinations producing multiple target outputs (i.e. the multioutput solution set, or MOSS) is likely much smaller. The relation between single-and multi-output solutions has been previously considered 11,50,56 but is not widely appreciated, nor has its implications been thoroughly explored.
Using computational modeling in which several neuronal properties are homeostatically regulated through feedback control of various ion channel densities, we show that the MOSS corresponds to the intersection between SOSSs and that MOSS dimensionality corresponds to the difference between the number of controllable ion channels (n in ) and the number of regulated properties (n out ). We also show that low-dimensional solutions encourage ion channel correlations and pose an increased risk for regulation to fail. Ion channel correlations and regulation failure are mitigated by having large (high-dimensional) SOSSs made possible by ion channel diversity, but, as our computational results demonstrate, regulating multiple properties in a coordinated manner requires even greater diversity in order for the MOSS to be degenerate and, therefore, robust.

RESULTS
In even the best characterized neuron types, it is not fully understood which properties are directly regulated and which ion channels are involved. We therefore studied basic principles using a generic model neuron in which multiple properties are regulated via control of the same ion channels. The relative rates at which ion channels are updated using the "error" in each property differ across properties, consistent with different feedback signals being used. We focused on independent properties like firing rate and energy efficiency per spike instead of, for example, energy consumption rate, which depends on firing rate. Our conclusions do not hinge on which properties are considered. Target values, or set points, were chosen arbitrarily.
Simulations were conducted in a single-compartment model neuron whose spikes are generated by a fast sodium conductance and a delayed rectifier potassium conductance with fixed densities. A spike-dependent adaptation mechanism was included at a fixed density. Densities of all other conductances were either systematically varied (and values producing the target output were selected) or they were adjusted via negative feedback to produce the target output (see Methods). The former approach, or brute force grid search, identifies all density combinations producing a target output (i.e. the level set or SOSS). The latter approach finds a subset of those combinations through a biologically plausible mechanism. Ion channels with controllable densities included a generic sodium conductance ( "# ) and potassium conductance ( $ ) that both activate with kinetics similar to the delayed rectifier channel, a slower-activating M-type potassium conductance ( % ), and a leak conductance ( &'#( ).

Most of the ion channel combinations producing a consistent value for one property produce
inconsistent values for other properties. We began by testing whether ion channel density combinations that produce the target value for one property also produce a consistent value for other properties. Figure 1A shows rheobase, defined here as the minimum injected current (I stim ) required to evoke repetitive spiking, for different combinations of "# and $ . Combinations yielding a rheobase of 30 µA/cm 2 are highlighted with a contour extending from point a ( "# = 1.54 mS/cm 2 , $ = 0 mS/cm 2 ) to point b ( "# = 4 mS/cm 2 , $ = 2.95 mS/cm 2 ); that contour represents the SOSS. Despite yielding the same rheobase, density combinations along the a-b contour did not yield the same minimum sustainable firing rate (f min ) (Fig 1B). The value of f min reflects spike initiation dynamics: f min >>0 spk/s is consistent with class 2 excitability 57 and operation as a coincidence detector 58 , whereas f min » 0 spk/s is consistent with class 1 excitability and operation as an integrator. Figure 1C shows a family of input-output curves for different models along the a-b contour, including models with density combinations a (pink) and b (cyan); these models exhibited notable differences in their operating mode, membrane potential oscillations, and spiking pattern (top, middle and bottom insets, respectively). After adding )*+ to expand the stimulus range over which fluctuation-driven spikes occur, we plotted the firing rate driven by irregularly fluctuating (noisy) I stim for different combinations of "# and $ (Fig.   1D). Density combinations yielding the same rheobase did not yield equivalent stimulation-evoked firing rates (Fig. 1E).
Having demonstrated that ion channel density combinations yielding the target value for one aspect of excitability (e.g. rheobase) yield differing values for other aspects of excitability (e.g. f min or firing rate), we predicted that the same lack of generalization would extend to other cellular properties (e.g. energy efficiency). Ion concentration gradients are depleted by the transmembrane ion flux associated with spikes and must be replenished directly or indirectly via energy-dependent pumps. Indeed, spikes are energetically costly 59 but vary in their energetic efficiency based on the temporal overlap in sodium and potassium channel activation 60 (Fig. S1).
Using the responses reported in Figure 1D, we measured energy consumption rate for combinations of "# and $ ( Fig. 2A). Energy consumption rate increased with firing rate, as expected, but did not increase equivalently across all density combinations, as evident from the variation in energy consumption rate along the iso-firing rate contour (Fig. 2B). This is due to differences in energy efficiency (Fig. 2C), which was determined as the energy consumed per spike relative to the theoretical minimum (see Methods). Density combinations yielding equally efficient spikes yielded very different stimulation-evoked firing rates and rheobase values (Fig. 2D).
One might presume that spikes are produced as efficiently as possible and, therefore, that energy efficiency is maintained above some lower bound rather than being regulated to a specific target value, and indeed this can be modeled (see Fig. 8). But one must also consider that energy efficiency trades off with functional considerations like propagation safety factor 61 . Target values might emerge from this competition 62 but such mechanisms were beyond the scope of the current study.

Regulation of one property risks disturbing other properties.
Based on the results described above, we predicted that compensatory ion channel changes that restore one property to its target value following a perturbation are likely to disturb other properties. We tested this by regulating firing rate via control of either of two ion channels (Fig. 3A) after removal of a third channel.
Grey circles on Figure 3B show randomly chosen combinations of "# and &'#( that yield a firing rate of 40 spk/s when $ = 2 mS/cm 2 . When $ was abruptly reset to 0 mS/cm 2 , mimicking blockade of the channel, firing rate jumped to ~93 spk/s before being restored to 40 spk/s via feedback control of either "# (pink) or &'#( (cyan) (Fig. 3C). As "# or &'#( converged on their new (compensated) densities, firing rate returned to its target value but energy efficiency was affected in opposite ways (Fig. 3D). If energy efficiency is also regulated, then the compensatory ion channel changes capable of restoring firing rate without disrupting energy efficiency is far more limited than the changes capable of restoring firing rate without regard for energy efficiency.

Geometrical explanation for the relationship between single-and multi-output solutions.
The disruption of one property by changes serving to regulate another property prompted us to ask how two or more properties are co-regulated. The left panel of Figure 4A shows all combinations of "# and $ (n in = 2) yielding the target value for firing rate (red) or energy efficiency (green) (n out = 1). Like in Figures 1-3, the SOSS for each property corresponds to a contour, or curve. The MOSS for firing rate and energy efficiency (n out = 2) corresponds to where the two curves intersect, which occurs at a point, indicating that only one density combination yields the target values for both properties -the MOSS is unique. But if another conductance like g leak is controllable (n in = 3), the curves in 2D parameter space (left panel) transform into surfaces in 3D parameter space (right panel) and those surfaces intersect along a curve, meaning many density combinations yield the target values for both properties -the MOSS becomes degenerate. The same patterns are evident for the MOSS for firing rate and input resistance (Fig. 4B) or energy efficiency and input resistance (Fig. 4C). Figure 4D shows that if all three properties -firing rate, energy efficiency and input resistance -are regulated (n out = 3), the MOSS is empty for n in = 2 (i.e. the three curves do not intersect at a common point [left panel] unless &'#( is reset to 1.95 mS/cm 2 [inset]) and the MOSS is unique for n in = 3 (right panel). We conclude that coordinated regulation of n properties requires at least n controllable ion channels, but n+1 controllable ion channels are required for the MOSS to be degenerate. This is like a system of linear equations, which is said to be underdetermined if there are fewer equations (outputs) than unknowns (inputs) (see Discussion).
If properties are regulated to within certain bounds rather than to a precise target value, SOSSs correspond to thin strips (rather than curves) or shallow volumes (rather than surfaces) in 2D and 3D parameter space, respectively (Fig. S2). In other words, less precise regulation allows SOSSs to expand in proportion to the tolerance. How does this impact the MOSS? Thin 2D strips intersect at a small patch (unlike curves intersecting at a point) but that small 2D patch in 2D parameter space (Fig. S2A) is unlike the long 1D curve formed by broad 2D surfaces intersecting in 3D parameter space. Likewise, shallow 3D volumes intersect at a narrow tube (unlike surfaces intersecting at a curve) but that 3D tube in 3D parameter space (Fig S2B) is unlike the broad 2D surface formed by deep 3D volumes intersecting in 4D parameter space. To summarize, increased tolerance does not expand MOSS dimensionality the same way as increased n in . On the other hand, if a property is maintained above or below some bound (see Fig. 8), the solution space is higher dimensional than for regulation to a target value under otherwise equivalent conditions. 8 The dimensionality of solution space affects ion channel correlations. We predicted that the dimensionality of solution space -point (0D), curve (1D), surface (2D), volume (3D), etc.affects ion channel correlations by limiting the degrees of freedom. To explore this, we measured ion channel correlations within SOSSs and MOSSs found through homeostatic feedback. Figure   5A shows all combinations of "# and $ (n in = 2) producing the target firing rate (n out = 1).
Correlation between "# and $ is high because homeostatically determined solutions are constrained to fall along a curve; under those conditions, variation in one channel is offset solely by variation in the other channel, and all covariance is thus captured in a single pairwise relationship. If % is also allowed to vary (n in = 3), homeostatically determined solutions spread across a surface and pairwise correlations become predictably weaker (Fig. 5B) since variation in one channel can be offset by variation in two other channels, and covariance is thus diluted across >1 pairwise relationship. If additional channels are allowed to vary (n in ³ 4), solutions distribute over higher-dimensional manifolds and pairwise correlations further weaken ( Fig. 5C). However, if firing rate and energy efficiency are both regulated (n out = 2; Fig. 5D), homeostatically determined solutions are once again constrained to fall along a curve when n in = 3, and pairwise correlations strengthen (Fig. 5E). Increasing n in to 4 while keeping n out at 2 caused correlations to weaken (Fig. 5F). These results confirm the predicted impact of solution space dimensionality on ion channel correlations.
O'Leary et al. 63 demonstrated how the relative rates at which different ion channel densities are controlled impact their correlation. This is reproduced in Figure 6A, where, from the same initial conditions (density combinations), homeostatic control with different relative rates (shown in pink and cyan) produce solutions with very different correlations. Relative control rates can affect not only the strength of pairwise correlations, but also the sign (compare correlation between "# and % ). However, if firing rate and energy efficiency are both homeostatically regulated, correlations strengthen (consistent with results from Fig. 5) and become independent of the relative control rates because solutions are limited to a lower-dimensional solution space ( Fig. 6B). Recall that MOSS dimensionality corresponds to n in -n out . For relative control rates to influence ion channel correlations, homeostatically determined solutions must fall on a solution space with dimensionality >1, but not >>1 lest pairwise correlations be diluted.

The dimensionality of solution space can also affect the success of homeostatic regulation.
Beyond affecting the ion channel correlations that emerge through homeostatic regulation, we predicted that the dimensionality of solution space affects whether multiple properties can be successfully regulated to their target values. Figure 7A shows an example in which control of "# , $ , and % successfully regulates firing rate. Figure 7B shows successful regulation of energy efficiency via control of the same three channels. Using the same initial conditions and relative control rates as in Figures 7A and B, the system failed to regulate firing rate and energy efficiency (Fig. 7C). Notably, coordinated regulation of both outputs was achieved using other relative control rates (see Figs. 5C and 6B) or using the same relative rates but starting from other initial conditions (not shown), suggesting that low-dimensional solutions are less 'accessible' (i.e. a smaller set of relative control rates succeed in finding the solution space). Figure 8 shows additional examples in which control of "# , $ , and % regulates firing rate and energy efficiency. In these simulations, firing rate is regulated to a precise target value whereas energy efficiency is maintained above a lower bound; the SOSS for energy efficiency thus corresponds to a volume rather than a surface for n in = 3. For energy efficiency ≥22%, homeostatically determined solutions converge on the iso-firing rate surface without regulation of energy efficiency having much effect (Fig. 8A). For energy efficiency ≥27%, solutions initially converge onto part of the iso-firing rate surface that sits outside the targeted energy efficiency volume, but, without energy efficiency requirements having been met, the solutions then move across the surface until reaching the intersection with the energy efficiency volume (Fig. 8B). By converging on the intersection, which is a curve, ion channel correlations are stronger than in Figure 8A. For energy efficiency ≥30%, the iso-firing rate surface and targeted energy efficiency volume do not intersect and solutions thus proceed to a point between the two SOSSs, constrained by $ and % reaching 0 mS/cm 2 (Fig. 8C). In this last example, neither property is optimally regulated but the outcome is a reasonable compromise. In other scenarios, regulation might have failed outright, e.g. if, in the absence of an upper bound, ion channel densities continued to increase without properties ever reaching their target values. This highlights that coordinated regulation of multiple properties might fail (or find sub-optimal solutions) not because singleoutput solutions do not exist, but because even large SOSSs might not (all) intersect, thus producing an empty MOSS.

DISCUSSION
This study has explored how neurons manage to homeostatically regulate multiple properties.
Being able to produce the same output using diverse ion channel combinations is key insofar as degeneracy allows solution sets to be large. But of the many ion channel combinations that produce the desired output for one property, only a subset also produce the desired output for other properties ( Figs. 1 and 2). Specifically, the multi-output solution set (MOSS) corresponds to the intersection between component single-output solution sets (SOSSs) (Fig. 4). The dimensionality of solution space therefore corresponds to the difference between the number of controllable ion channels (n in ) and the number of regulated outputs (n out ). Coordinated regulation of n properties thus requires at least n controllable ion channels for a unique solution, and at least n+1 channels for a degenerate solution.
A direct analogy can be made with a system of linear equations. Each unknown constitutes a degree of freedom and each equation constitutes a constraint that reduces the degrees of freedom by one. The system is said to be overdetermined if equations outnumber unknowns, and underdetermined if unknowns outnumber equations. An underdetermined system can have infinite solutions, and while frustrating mathematically, a biological system benefits from having multiple (degenerate) solutions available to it. 'Extra' degrees of freedom offer a broader range of solutions, which means a good solution is liable to exist over a broader range of conditions. If constraints outnumber the degrees of freedom, the system becomes overdetermined and solutions disappear (Fig. 4D), which can cause regulation to fail (see below).
There are notable similarities and differences between a neuron adjusting its ion channel densities to regulate properties to their target values and a scientist trying to infer conductance densities based on the measured values of those properties (i.e. fitting a model to experimental data). Fitting a model with n parameters to many outputs (e.g. firing rate and input resistance and rheobase and spike height) is more difficult but yields better parameter estimates than fitting the same model to just one output 56 , just as regulating more properties leads to a smaller MOSS. But how good are those parameter estimates? Can one confidently infer the true conductance densities in a particular neuron by measuring and fitting enough properties? Degeneracy makes solving this inverse problem difficult, if not impossible 64 . The neuron solves the forward problem, producing a firing rate, input resistance, etc. based on its ion channel density combination. That density combination is determined by negative feedback, but the precise density combination is not important so long as it falls within the MOSS (i.e. produces the target values for all regulated properties). As conditions change, the MOSS will change and negative feedback will adjust the density combination accordingly. A neuron needs to regulate its properties robustly but does not do so by regulating conductance densities to target values; in that respect, a neuron does not solve an inverse problem. By this logic, neurons have evolved under selective pressure to be degenerate, not parsimonious, contrary to how scientists typically construct their models.
Previous studies using brute force searches to explore degeneracy have tended to apply selection criteria simultaneously, thus finding the MOSS in one fell swoop. In contrast, we considered one criterion (property) at a time in order to find each SOSS, which we then combined to find the MOSS. The former approach is akin to aggregating several error functions to create a single-objective problem, whereas the latter resembles multi-objective optimization, which can reveal trade-offs between objectives 65 . Nevertheless, past studies have observed that certain parameter changes (in a particular direction through parameter space) dramatically affect some properties but not others 34,66 , which is consistent with manifolds representing the SOSS for 'sensitive' and 'insensitive' properties lying orthogonal to one another in parameter space. The impact of solution dimensionality on ion channel correlations and regulation failure (see below) highlight the value of a 'multi-objective' perspective, which implicitly recognizes that different properties are regulated independently but not in isolation.
The relative rates with which different conductance densities are controlled dictate the direction in which trajectories move through parameter space, which in turn dictates how trajectories approach and distribute across the solution manifold 63 . Our results are consistent with that explanation (Fig. 6) but also highlight the importance of manifold dimensionality (Fig. 5). If trajectories reach a 1D manifold (curve), correlations will be stronger (but less sensitive to relative control rates) than if they reach a 2D manifold (surface). Pairwise correlations will be even weaker if the manifold is higher-dimensional (volume, etc.). The existence of pairwise correlations thus suggests that the dimensionality of solution space (= n in -n out ) is relatively low. This may be surprising since n in is presumably high, but makes sense if n out is also high; in other words, there are many channels, but they are responsible for regulating many properties.
Homeostatic regulation can fail for different reasons. If there are no solutions (i.e. the solution set is empty), negative feedback cannot regulate properties to their target values. A MOSS can be empty because SOSSs do not intersect (Fig. 8C), meaning regulation of different properties to their respective target values is incompatible. This is distinct from a SOSS being empty, which means regulation of a particular property is impossible. Regulation can also fail because negative feedback fails to converge on available solutions (Fig. 7C). Notably, a MOSS may become less 'accessible' (lower dimensional) as SOSSs start to separate (Fig. 8B), and may foreshadow the eventual disjunction of SOSSs (Fig. 8C). Regulation can also fail because feedback signaling is compromised. These failure modes are not mutually exclusive; for example, a system might reach less 'accessible' solutions if control rates are normally flexible, but might fail to reach those solutions if control rate flexibility is reduced. Failed homeostatic regulation may have similar consequence regardless of how the failure occurs, consistent with the emergence of common disease phenotypes despite vastly different underlying pathologies 75 . That said, subtle differences in exactly how the failure transpires may provide clues to help pinpoint the underlying mechanism(s).
In conclusion, neurons manage to homeostatically regulate multiple properties thanks to ion channel degeneracy. This may account for why ion channels have evolved to a point where their diversity seems superfluous -'excess' degeneracy accommodates complex regulation of diverse properties. Tantalizing clues support this possibility 76  Energy calculations. Energy consumption rate was calculated as in Hasenstaub et al. 79 Briefly, models were stimulated with a fast fluctuating stimulus (µ stim = 40 µA/cm 2 , σ stim = 10 µA/cm 2 ).
Sodium and potassium current through all channels was integrated for 1 second to determine the charge for each ion species, which was then converted to ion flux based on the elementary charge 1.602 x 10 -19 C. Based on the 3:2 stoichiometry of the Na + /K + pump, we divided the number of sodium and potassium ions by 3 and 2, respectively, and used the maximum of those two values as the energy consumption rate (in ATP/cm 2 ·s). Energy efficiency was calculated as the ratio between capacitive minimum and total Na + flux during an action potential 60 . Briefly, voltage was reset to -40 mV to evoke a single spike in all models. Models were treated as pure capacitors to calculate the minimum capacitive current as CΔV, where C is the capacitance and ΔV is the difference between the resting membrane potential and the spike peak.
Brute force grid search. Models were tested with conductance density combinations chosen from a 100x100 or 30x30x30 grid for 2D and 3D plots, respectively. To depict the SOSS, all models with outputs within a range (tolerance) of the target value (± 3 spk/s for firing rate, ± 0.25 % for energy efficiency, and ± 0.003 kΩ·cm 2 for input resistance) were selected and a curve or surface was fit to those successful models. The same tolerances were implemented in the homeostatic learning rule (see below) to minimize ringing. Tolerances are illustrated in Figure S2 but were not shown in other figures for sake of clarity.   Table S1 for the initial conductance values and control rates used for each figure.
All code is available at http://prescottlab.ca/code-for-models.  A Color shows the minimum I stim required to evoke repetitive spiking (rheobase) for different combinations of ḡ Na and ḡ K . Contour linking a and b highlights density combinations yielding a rheobase of 30 μA/cm 2 . B Minimum sustainable firing rate (f min ) varied along the iso-rheobase contour. C Firing rate is shown as a function of I stim for combination a (pink), b (cyan), and other values along the iso-rheobase contour a -b tested at 0.5 mS/cm 2 increments in ḡ Na (black). Grey arrows point to responses to noisy stimulation modeled as an Ornstein-Uhlenbeck (OU) process with τ stim = 5 ms and σ stim = 5, 1 or 2 µA/cm 2 and μ stim = 25, 29.5 or 35 μA/cm 2 for sub-, peri-and suprathreshold conditions (top, middle, and bottom insets, respectively). Top: Spike-triggered averages (STAs) show that model a operates as an integrator with a longer integration time than model b, which operates as a coincidence detector. Middle: Noise-induced membrane potential oscillations were slower for model a than for model b, consistent with differences in f min , and average membrane potential differed by >10 mV. Bottom: spike-dependent adaptation (ḡ AHP = 2 mS/cm 2 ) caused bursting in model b but not in model a. Voltage threshold also differed by 10 mV. D Color shows firing rate evoked by noisy I stim (τ stim = 5 ms, σ stim = 10 μA/cm 2 , μ stim = 40 μA/cm 2 ) for different combinations of ḡ Na and ḡ K . Adding spike-dependent and -independent forms of adaptation (ḡ AHP = 1.75 mS/cm 2 and ḡ M = 0.5 mS/cm 2 ) broadened the dynamic range. Grey curve shows density combinations yielding a firing rate of 40 spk/s (i.e. an iso-firing rate contour; shown in red in Energy Efficiency (%) Figure 2. Different ion channel combinations yielding equivalent firing rate do not produce energetically equivalent spikes. A Color shows energy consumption rate for different combinations of ḡ Na and ḡ K based on firing rates shown in Figure  1D. Iso-firing rate contour c-d (grey line in Fig. 1D) does not align with energy contours. ATP consumed by the Na + /K + pump was calculated from the total Na + influx and K + efflux determined from the corresponding currents. B Energy consumption rate varied along the iso-firing rate contour. C Color shows energy efficiency per spike for different combinations of ḡ Na and ḡ K . Energy efficiency was calculated as the capacitive minimum C·ΔV divided by total Na + influx, where C is capacitance and ΔV is spike amplitude (see Fig. S1). Density combinations along contour i-ii (dashed line) yield energy efficiency of 23.5% (shown in green in Fig. 4). D Both rheobase and firing rate varied along the iso-energy efficiency contour. A Schematic shows how a difference in firing rate from its target value creates an error that is reduced by updating ḡ Na or ḡ leak . B Iso-firing rate contours for 40 spk/s are shown for different combinations of ḡ Na and ḡ leak with ḡ K at baseline (2 mS/cm 2 , dashed red curve) and after ḡ K was "knocked out" (0 mS/cm 2 , solid red curve). When ḡ K was abruptly reduced, starting models (grey dots) spiked rapidly (~93 spk/s) before firing rate was regulated back to its target value by compensatory changes in either ḡ Na (pink) or ḡ leak (cyan). Models evolved in different directions and settled at different positions along the solid curve. C Trajectories show evolution of ḡ Na and ḡ leak . Trajectories are terminated once target firing rate is reached. Sample traces show responses before (grey) and immediately after (black) ḡ K was reduced, and again after compensatory changes in ḡ Na (pink) or ḡ leak (cyan). D Distributions of energy efficiency are shown before (grey) and after firing rate regulation via control of ḡ Na (pink) or ḡ leak (cyan).  A Starting from a normally distributed cluster of ḡ Na and ḡ K (white dots) yielding an average firing rate of 73 spk/s, ḡ Na and ḡ K (n in = 2) were homeostatically adjusted to regulate firing rate (n out = 1) to its target value of 40 spk/s. Grey lines show trajectories. Because solutions (grey dots) converge on a curve, the pairwise correlation between ḡ Na and ḡ K is predictably strong. Scatterplots show solutions centered on the mean and normalized by the standard deviation (z-scores). Correlation coefficient (R) is shown in the bottom right corner of each scatter plot. B Same as a (n out = 1) but via control of ḡ Na , ḡ K and ḡ M (n in = 3). Homeostatically found solutions converge on a surface and ion channel correlations are thus weaker. C If ḡ leak is also controlled (n in = 4), solutions converge on a hard-to-visualize volume (not shown) and pairwise correlations are further weakened. D Schematic shows how errors for two regulated properties are combined: the error for each property is calculated separately and is scaled by its respective control rate τ to calculate updates, and all updates for a given ion channel (i.e. originating from each error signal) are summed. E Same as B (n in = 3), but for regulation of firing rate and energy efficiency (n out = 2). Homeostatically found solutions once again converge on a curve (yellow), which now corresponds to the intersection of two surfaces; ion channel correlations are thus strong, like in A. F If n in is increased to 4 while n out remains at 2, solutions converge on a surface (not shown) and ion channel correlations weaken.  Figure 6. The dimensionality of solution space and relative control rates shape ion channel correlations. A Homeostatic regulation of firing rate via control of ḡ Na , ḡ K and ḡ M . Same as Figure 5B, but for two new sets of control rates (see Table S1).
Solutions found for each set of rates (pink and cyan) approached the surface from different angles and converged on the surface with different patterns, thus producing distinct ion channel correlations (consistent with O'leary et al. 63 ). B Same as A (n in = 3), but for homeostatic regulation of firing rate and energy efficiency (n out = 2). Solutions converge on a curve (yellow), giving rise to virtually identical ion channel correlations regardless of control rates. ḡ Na 4 2 0 ḡ K ḡ Na ḡ K ḡ M Figure 7. Low-dimensional solutions can be hard for homeostatic regulation to "find". A Homeostatic regulation of firing rate (n out = 1) via control of ḡ Na , ḡ K and ḡ M (n in = 3), like in Figures 5B and 6A but using a different set of control rates (see Table  S1). Solutions converged onto the iso-firing rate surface (top panel) and firing rate was regulated to its target value in <30 iterations (bottom panel). B Same as A (n in = 3) but for homeostatic regulation of energy efficiency (n out = 1). Solutions converged on the iso-energy efficiency surface (top panel) and energy efficiency was regulated to its target value in ~10 iterations (bottom panel). C Homeostatic regulation of firing rate and energy efficiency (n out = 2) via control of the same ion channels (n in = 3) using the same relative rates and initial values as in A and B. Neither firing rate (red trajectories) nor energy efficiency (green trajectories) reached its target value. Conductance densities were capped at 4 mS/cm 2 but this does not account for trajectories not reaching their target. Homeostatic regulation of firing rate and energy efficiency (n out = 2) via control of ḡ Na , ḡ K and ḡ M (n in = 3). For these simulations, energy efficiency was maintained above a lower bound rather than being regulated to a specific target value; accordingly, the SOSS for energy efficiency corresponds to a volume (green) rather than a surface. A For energy efficiency ≥22%, homeostatically determined solutions converge on the iso-firing rate surface (red) in a region sitting within the green volume (top panel). The rate of convergence and resulting ion channel correlations are shown in the middle and bottom panels, respectively. B For energy efficiency ≥27%, solutions initially converge on the red surface in a region outside the green volume, but trajectories then bend and proceed across the red surface until that surface reaches the green volume. Because solutions converge on a curve, ion channel correlations are stronger than in A, where solutions distributed across a surface. C For energy efficiency ≥30%, the red surface and green volume do not intersect. Consequently, solutions settle between the two SOSSs (top panel) without either property reaching its target (middle panel). The outcome represents the balance achieved by the opposing pull of control mechanisms regulating different properties, and depends entirely on ḡ Na since ḡ K and ḡ M cannot become negative (bottom panel). K + Na + Figure S1. Overlap between Na + and K + currents dictates energy efficiency. Voltage, total Na + current and total K + current during an action potential are shown for the most efficient (left) and least efficient (right) models. Simultaneous activation of Na + and K + channels creates "waste" current because Na + influx and K + efflux cancel each other out. Overlap between Na + and K + currents is larger in the model on the right, meaning more current is used relative to the theoretical minimum required to spike, which is calculated as CΔV, where C is the membrane capacitance and ΔV is the difference between the resting membrane potential and the spike peak.  Figure S2. Increasing tolerance does not increase MOSS dimensionality the same way as increasing n in . A Curves for firing rate (40 spk/s) and energy efficiency (23.5%) expand into strips when tolerance (±3 spk/s and ±0.25%, respectively) is depicted. The strips intersect as a patch (yellow highlighting), unlike curves which intersect at a point (see Fig. 4A, left panel). Note that 2D strips are limited to a 2D parameter space. The broad 2D surfaces that exist in 3D parameter space intersect along a long curve (see Fig. 4A, right panel), unlike the small patch. B If ḡ Leak is also controllable (n in = 3), strips in 2D parameter space expand into shallow volumes in 3D parameter space. Those shallow volumes intersect along a narrow tube (yellow highlighting), which is unlike the broad surface formed by the intersection of deep volumes in 4D parameter space (not illustrated). Table S1. Initial conductance values and regulation rates for each figure.