Minimal requirements for a neuron to coregulate many properties and the implications for ion channel correlations and robustness

  1. Jane Yang
  2. Husain Shakil
  3. Stéphanie Ratté
  4. Steven A Prescott  Is a corresponding author
  1. Neurosciences and Mental Health, The Hospital for Sick Children, Canada
  2. Institute of Biomedical Engineering, University of Toronto, Canada
  3. Department of Physiology, University of Toronto, Canada

Abstract

Neurons regulate their excitability by adjusting their ion channel levels. Degeneracy – achieving equivalent outcomes (excitability) using different solutions (channel combinations) – facilitates this regulation by enabling a disruptive change in one channel to be offset by compensatory changes in other channels. But neurons must coregulate many properties. Pleiotropy – the impact of one channel on more than one property – complicates regulation because a compensatory ion channel change that restores one property to its target value often disrupts other properties. How then does a neuron simultaneously regulate multiple properties? Here, we demonstrate that of the many 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 multioutput solution set) correspond to the intersection between single-output solution sets. Properties can be effectively coregulated only if the number of adjustable channels (nin) exceeds the number of regulated properties (nout). Ion channel correlations emerge during homeostatic regulation when the dimensionality of solution space (ninnout) is low. Even if each property can be regulated to its target value when considered in isolation, regulation as a whole fails if single-output solution sets do not intersect. Our results also highlight that ion channels must be coadjusted with different ratios to regulate different properties, which suggests that each error signal drives modulatory changes independently, despite those changes ultimately affecting the same ion channels.

Editor's evaluation

Neurons develop and maintain rich electrophysiological properties that enable nervous systems to function. The question of how different neural properties are regulated by internal ion channel expression mechanisms remains unresolved. In this paper, Yang and colleagues address the question from an abstract perspective by asking how multiple constraints on the physiological properties of neurons, such as firing rate curves and energy efficiency, narrow down the available regulatory possibilities. Their results from a mixture of modelling and dynamic clamp experiments point to the existence of multiple parallel internal feedback loops for controlling ion channel expression in neurons, and derive conditions under which co-regulation mechanisms will fail.

https://doi.org/10.7554/eLife.72875.sa0

Introduction

Neurons maintain their average firing rate near a target value, or set point, by adjusting their intrinsic excitability and synaptic weights (Turrigiano et al., 1994; O’Leary et al., 2010; Aizenman et al., 2003; Desai et al., 1999; Hengen et al., 2013; van Welie et al., 2006). Homeostatic regulation of intrinsic excitability is achieved through feedback control of diverse ion channels (Turrigiano et al., 1994; O’Leary et al., 2010; Desai et al., 1999; Joseph and Turrigiano, 2017). Computational models have successfully employed negative feedback to adjust ion channel densities (O’Leary et al., 2014; O’Leary and Marder, 2016; Liu et al., 1998; Olypher and Prinz, 2010; LeMasson et al., 1993) and control theory provides a valuable framework to conceptualize how this occurs (O’Leary and Wyllie, 2011). But most of the mechanistic details remain unclear (Davis, 2006; Turrigiano, 2011) and are not straightforward; for instance, different perturbations can trigger similar changes in excitability via different signaling pathways affecting different ion channels (Kulik et al., 2019), or via different combinations of excitability changes and synaptic scaling (Maffei and Turrigiano, 2008).

Firing rate homeostasis is facilitated by the ability of different ion channel combinations to produce equivalent excitability (Marder, 2011). The ability of distinct elements to produce the same outcome is known as degeneracy (Edelman and Gally, 2001) and has attracted increasing attention in neuroscience (Marder, 2011; Marder and Goaillard, 2006; Cropper et al., 2016; O’Leary, 2018; Ratté and Prescott, 2016; Tononi et al., 1999; Goaillard and Dufour, 2014; Mason et al., 2015; Rathour and Narayanan, 2019). Many aspects of neural function at the genetic (Klassen et al., 2011; Trojanowski et al., 2014), synaptic (Mukunda and Narayanan, 2017; Anirudhan and Narayanan, 2015), cellular (Kim et al., 2017; Taylor et al., 2009; Drion et al., 2015; Ratté et al., 2014b; Mittal and Narayanan, 2018; Jain and Narayanan, 2020; Günay et al., 2008; Migliore et al., 2018), and network (Grashow et al., 2010; Marder and Taylor, 2011; Prinz et al., 2004; Price and Friston, 2002; Onasch and Gjorgjieva, 2020) levels are now recognized as being degenerate. Ion channel degeneracy facilitates robust homeostatic regulation of neuronal excitability by enabling a disruptive change in one ion channel to be offset by compensatory changes in other ion channels (O’Leary, 2018; Drion et al., 2015; Ratté et al., 2014b; Rho and Prescott, 2012; Swensen and Bean, 2005; Zhao and Golowasch, 2012; Achard and De Schutter, 2006; Bhalla and Bower, 1993; Olypher and Calabrese, 2007).

But to maintain good neural coding, neurons must regulate various aspects of their activity beyond average firing rate (Stemmler and Koch, 1999), and must do so while also regulating their energy usage, osmolarity, pH, protein levels, etc. (Frere and Slutsky, 2018; Fisher et al., 2010; Balch et al., 2008). These other cellular properties affect and are affected by neuronal activity, which is to say that properties are not regulated in isolation from one another. And just as each property depends on multiple ion channels, each ion channel can affect multiple properties (Taylor et al., 2009; Günay et al., 2008). This pleiotropy or functional overlap (Goaillard and Marder, 2021) confounds the independent regulation of each property. In response to a perturbation, a compensatory ion channel change that restores one property to its target value may exacerbate rather than mitigate the disturbance of a second property (Figure 1A). In theory, both properties could be regulated by adjusting two or more ion channels (Figure 1B), but this suggests that the number of properties that can be coregulated is limited by the number of adjustable ion channels, even if the channels are pleiotropic.

Simultaneous regulation of >1 property in a theoretical neuron.

(A) Challenge: If two ion channels, gM (pink) and gAHP (purple), both affect firing rate (orange), then the change in firing rate caused by perturbing gM can be offset by a compensatory change in gAHP. However, if gM and gAHP also affect firing pattern (blue), but not in exactly the same way, then the compensatory change in gAHP that restores firing rate to its target value may exacerbate rather than resolve the change in firing pattern. (B) Solution: Adjusting gAHP and at least one additional conductance, gx (green), that affects each property in a different way than gAHP, may enable firing rate and firing pattern to be regulated back to their target values. This suggests that the number of adjustable ion channels (‘dials’) relative to the number of regulated properties is important. Note that adjusting gx to restore the firing pattern necessitates a small extra increase in gAHP (compare pale and dark purple); in other words, conductances must be coadjusted.

To identify the requirements for coregulating many properties, we began with an experiment to confirm that regulating one property by adjusting a single channel is liable to disrupt other properties. We then proceeded with modeling to unravel how multiple channels must be coadjusted to ensure proper regulation of >1 property. We show that the number of adjustable ion channels (nin) must exceed the number of regulated properties (nout). This is consistent with past work (Olypher and Prinz, 2010; Olypher and Calabrese, 2007; Foster et al., 1993) but the implications for homeostatic regulation have not been explored. To that end, we show that the dimensionality of solution space (ninnout) influences the emergence of ion channel correlations in the presence or absence of noise, and that increased correlations presage regulation failure. We also show that regulating different properties requires that ion channels are coadjusted with different ratios, which necessitates separate master regulators.

Results

Adjusting an ion channel to regulate one property risks disrupting other properties – experiments

To test experimentally if regulating one property by adjusting one ion channel is liable to disrupt a second property, we disrupted the firing rate of a CA1 pyramidal neuron by blocking its voltage-gated M-type K+ current (IM) and then we restored firing rate to its baseline value by inserting a virtual calcium-activated AHP-type K+ current (IAHP) using dynamic clamp, while also monitoring the firing pattern. Specifically, the neuron was stimulated by injecting irregularly fluctuating (noisy) current under four conditions (Figure 2A, B): at baseline (blue), after blocking native IM (black), and again after introducing virtual IM (cyan) or virtual IAHP (red). Adding virtual IM demonstrates that we can replace a native current with an equivalent virtual current; compensation was modeled by inserting a distinct current with functional overlap, namely IAHP (see Figure 1). Inserting virtual IM or virtual IAHP reversed the depolarization (Figure 2C), spike amplitude attenuation (Figure 2D), and firing rate increase (Figure 2E) caused by blocking native IM. Replacing native IM with virtual IM did not affect firing pattern, quantified as the coefficient of variation of the interspike interval (CVISI), whereas virtual IAHP reduced CVISI (Figure 2F), often causing the neuron to spike at different times (in response to different stimulus fluctuations) than the neuron with native or virtual IM (Figure 2G). The results confirm predictions from Figure 1, namely that compensatory ‘upregulation’ of IAHP restored firing rate but disrupted firing pattern.

A compensatory change that restores firing rate disrupts firing pattern in a CA1 pyramidal neuron.

(A) Rasters show spiking in a single CA1 pyramidal neuron given the same noisy current injection (gray) on five trials under each of four conditions: baseline (blue), after blocking native IM with 10 μM XE991 (black), and again after inserting virtual IM (cyan) or virtual IAHP (red) using dynamic clamp. See Methods for virtual conductance parameters. (B) Sample voltage traces under each condition. (C) Membrane potential differed across conditions (F3,16 = 298.49, p < 0.001, one-way ANOVA); specifically, it was depolarized by blocking IM (t = 18.71, p < 0.001, Tukey test) but that effect was reversed by inserting virtual IM (t = 24.56, p < 0.001) or virtual IAHP (t = 27.00, p < 0.001). (D) Spike amplitude also differed across conditions (F3,16 = 154.04, p < 0.001); specifically, it was attenuated by blocking IM (t = 20.23, p < 0.001) but that effect was reversed by inserting virtual IM (t = 13.99, p < 0.001) or virtual IAHP (t = 16.09, p < 0.001). (E) Firing rate also differed across conditions (F3,16 = 177.61, p < 0.001); specifically, it was increased by blocking IM (t = 18.38, p < 0.001) but that effect was reversed by inserting virtual IM (t = 20.85, p < 0.001) or virtual IAHP (t = 16.12, p < 0.001). (F) Regularity of spiking, reflected in the coefficient of variation of the interspike interval, also differed across conditions (F3,16 = 188.23, p < 0.001), but whereas blocking IM had a modest effect (t = 2.69, p = 0.048) and inserting virtual IM had no effect (t = 0.39, p = 1.00), inserting virtual IAHP had a large effect (t = 18.23, p < 0.001). Data are summarized as mean ± standard deviation (SD). Each data point represents a different trial from a single neuron. (G) Enlarged view of rasters to highlight spikes that occurred with native or virtual IM but not with virtual IAHP (blue shading) or vice versa (red shading) to illustrate the change in spike pattern caused by correcting the change in firing rate caused by blocking IM with a compensatory change in IAHP.

Figure 2—source data 1

Numerical values for experimental data plotted in Figure 2.

Data include all spikes times (for panel A) as well as trial-averaged values of membrane potential, spike amplitude, firing rate, and coefficient of variation of the interspike interval (for panels C–F, respectively).

https://cdn.elifesciences.org/articles/72875/elife-72875-fig2-data1-v2.xlsx

Few of the ion channel combinations producing the target value for one property also produce the target value for a second property

We proceeded with computational modeling to identify the conditions required to coregulate multiple properties. For our investigation, it was critical to account for all ion channels that are changing. Even for well characterized neurons like pyramidal cells, there is no comprehensive list of which properties are directly regulated and which ion channel are involved in each case. Therefore, we built a simple model neuron in which >1 property can be coregulated by adjusting a small number of ion channels. The relative rates at which ion channel densities are updated using the ‘error’ in each property differ between properties, meaning ion channels are coadjusted according to a certain ratio to regulate one property (O’Leary et al., 2014) but are coadjusted according to a different ratio to regulate a second property, which has some notable implications (see Discussion). We focused on regulation of independent properties like firing rate and energy efficiency per spike – instead of energy consumption rate, for example, which depends on firing rate – but our conclusions do not hinge on which properties are considered. Target values 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, gAHP was included at a fixed density. Densities of all other channels 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, which amounts to a grid search, identifies all density combinations producing a target output (i.e., the solution set). The latter approach finds a subset of those combinations through a homeostatic regulation mechanism. Ion channels with adjustable densities included a sodium conductance (gNa) and potassium conductance (gK) that both activate with kinetics similar to the delayed rectifier channel, a slower-activating M-type potassium conductance (gM), and a leak conductance (gleak).

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 3A shows rheobase for different combinations of g¯Na and g¯K. The contour from point a (g¯Na = 1.54 mS/cm2, g¯K = 0 mS/cm2) to point b (g¯Na = 4 mS/cm2, g¯K = 2.95 mS/cm2) represents all combinations (i.e., the solution set) yielding a rheobase of 30 μA/cm2. Despite yielding the same rheobase, density combinations along the a–b contour did not yield the same minimum sustainable firing rate (fmin) (Figure 3B). The value of fmin reflects spike initiation dynamics: fmin >> 0 spk/s is consistent with class 2 excitability (Hodgkin, 1948) and operation as a coincidence detector (Ratté et al., 2014a), whereas fmin≈ 0 spk/s is consistent with class 1 excitability and operation as an integrator, with several consequences (Figure 3—figure supplement 1). After adding gAHP to expand the stimulus range over which fluctuation-driven spikes occur, we plotted the firing rate driven by irregularly fluctuating (noisy) Istim for different combinations of g¯Na and g¯K (Figure 3C). Density combinations yielding the same rheobase did not yield equivalent stimulation-evoked firing rates (Figure 3D).

Figure 3 with 2 supplements see all
Most ion channel combinations producing the target value for one property produce inconsistent values for other properties.

(A) Color shows the minimum Istim required to evoke repetitive spiking (rheobase) for different combinations of g¯Na and g¯K. Contour linking a and b highlights density combinations yielding a rheobase of 30 μA/cm2. (B) Minimum sustainable firing rate (fmin) varied along the iso-rheobase contour (see Figure 3—figure supplement 1). (C) Color shows firing rate evoked by noisy Istimstim = 5ms, σstim = 10 μA/cm2, µstim = 40 μA/cm2) for different combinations of g¯Na and g¯K. Adding spike-dependent and -independent forms of adaptation (g¯AHP = 1.75 mS/cm2 and g¯M = 0.5 mS/cm2) broadened the dynamic range. Gray curve shows density combinations yielding a firing rate of 40 spk/s (i.e., an iso-firing rate contour; shown in red in Figures 46). Insets show sample responses to equivalent noisy stimulation. (D) Firing rate varied along the iso-rheobase contour from panel A. (E) Color shows energy consumption rate for different combinations of g¯Na and g¯K based on firing rates shown in panel C. Iso-firing rate contour c–d (gray line in panel C) 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. (F) Energy consumption rate varied along the iso-firing rate contour. (G) Color shows energy efficiency per spike for different combinations of g¯Na and g¯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 Figure 3—figure supplement 2). Density combinations along contour i–ii (dashed line) yield energy efficiency of 23.5% (shown in green in Figure 5). (H) Both rheobase and firing rate varied along the iso-energy efficiency contour.

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., fmin or firing rate), we predicted that the same lack of generalization would extend to other cellular properties (e.g., energy efficiency). Spikes are energetically costly (Attwell and Laughlin, 2001) but vary in their energy efficiency based on the temporal overlap in sodium and potassium channel activation (Sengupta et al., 2010; Figure 3—figure supplement 2). Using responses reported in Figure 3C, we measured energy consumption rate for combinations of g¯Na and g¯K (Figure 3E). 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 (Figure 3F). This is due to differences in energy efficiency (Figure 3G) determined as the energy consumed per spike relative to the theoretical minimum (see Methods). Density combinations yielding equally efficient spikes yielded different stimulation-evoked firing rates and rheobase values (Figure 3H). One might presume that spikes are produced as efficiently as possible rather than being regulated to a specific value, but energy efficiency decreases propagation safety factor (Al-Basha and Prescott, 2019) and a target value likely emerges by balancing these competing interests (Speakman et al., 2011).

Adjusting an ion channel to regulate one property risks disrupting other properties – simulations

Like the experiment in Figure 2, we tested if restoring firing rate to its target value via compensatory changes in either of two ion channels (Figure 4A) after perturbing a third channel would disrupt energy efficiency in our model neuron. Gray circles on Figure 4B show randomly chosen combinations of g¯Na and g¯leak that yield a firing rate of 40 spk/s when g¯K = 2 mS/cm2. When g¯K was ‘blocked’ (abruptly reset to 0 mS/cm2), firing rate jumped to ~93 spk/s before being restored to 40 spk/s via feedback control of either g¯Na (pink) or g¯leak (cyan) (Figure 4C). As g¯Na or g¯leak converged on new, compensated densities, firing rate returned to its target value but energy efficiency was affected in opposite ways (Figure 4D). If energy efficiency is also regulated, then g¯Na and g¯leak must be coadjusted in a way that restores firing rate without disrupting energy efficiency.

Ion channel changes mediating the same effect on firing rate can oppositely affect energy efficiency.

(A) Schematic shows how a difference in firing rate from its target value creates an error that is reduced by updating g¯Na or g¯leak. (B) Iso-firing rate contours for 40 spk/s are shown for different combinations of g¯Na and g¯leak with g¯K at baseline (2 mS/cm2, dashed red curve) and after g¯K was ‘knocked out’ (0 mS/cm2, solid red curve). When g¯K was abruptly reduced, starting models (gray dots) spiked rapidly (~93 spk/s) before firing rate was regulated back to its target value by compensatory changes in either g¯Na (pink) or g¯leak (cyan). Models evolved in different directions and settled at different positions along the solid curve. (C) Trajectories show evolution of g¯Na and g¯leak. Trajectories are terminated once target firing rate is reached. Sample traces show responses before (gray) and immediately after (black) g¯K was reduced, and again after compensatory changes in g¯Na (pink) or g¯leak (cyan). (D) Distributions of energy efficiency are shown before (gray) and after firing rate regulation via control of g¯Na (pink) or g¯leak (cyan).

Geometrical explanation for the relationship between single- and multioutput solutions

Next, we sought a geometrical explanation for how to coadjust >1 ion channel to coregulate >1 property. The top panel of Figure 5A shows all combinations of g¯Na and g¯K (nin = 2) yielding the target value for firing rate (red) or energy efficiency (green) (nout = 1). Like in Figures 3 and 4, the solution set for each property corresponds to a curve. The multioutput solution set for firing rate and energy efficiency (nout = 2) corresponds to where the two single-output solution sets intersect, which occurs at a point – only one density combination yields the target values for both properties, meaning the multioutput solution is unique. But if another conductance like gleak is adjustable (nin = 3), the curves in 2D parameter space (top panel) transform into surfaces in 3D parameter space (bottom panel) and those surfaces intersect along a curve – many density combinations yield the target values for both properties, meaning the multioutput solution is degenerate. The same patterns are evident for firing rate and input resistance (Figure 5B) or energy efficiency and input resistance (Figure 5C). Figure 5D shows that if all three properties – firing rate, energy efficiency and input resistance – are regulated (nout = 3), the multioutput solution set is empty for nin = 2 (i.e., the three curves do not intersect at a common point [top panel] unless g¯leak is reset to 1.95 mS/cm2 [inset]) and is unique for nin = 3 (bottom panel). From this, we conclude that robust regulation of n properties requires >n adjustable ion channels even if each ion channel contributes to regulation of >1 property. This is like a system of linear equations, which is said to be underdetermined if unknowns (inputs) outnumber equations (outputs) (see Discussion).

Figure 5 with 1 supplement see all
A degenerate solution for n properties requires at least n + 1 adjustable ion channels.

Curves in top panels depict single-output solutions sets based on control of g¯Na and g¯K (nin = 2); g¯leak = 2 mS/cm2. Surfaces in bottom panels depict single-output solution sets based on control of g¯Na, g¯K, and g¯leak (nin = 3). Intersection (yellow) of single-output solutions at a point constitutes a unique multioutput solution, whereas intersection along a curve (or higher-dimensional manifold) constitutes a degenerate multioutput solution. (A) Curves for firing rate (40 spk/s) and energy efficiency (23.5%) intersect at a point whereas the corresponding surfaces intersect along a curve. Solutions for firing rate and input resistance (0.65 kΩ cm2) (B) and for energy efficiency and input resistance (C) follow the same pattern as in panel A. (D) For nin = 2 (top), curves for firing rate, energy efficiency and input resistance do not intersect at a common point unless g¯leak is reset to 1.95 mS/cm2 (inset). For nin = 3 (bottom), the three surfaces intersect at the same point as in the inset. See Figure 5—figure supplement 1 for the effects of tolerance on solution sets.

In all simulations involving homeostatic regulation, properties were regulated to within certain bounds of the target rather than to a precise value, despite how solution sets are depicted. When tolerances are depicted, single-output solution sets correspond to thin strips (rather than curves) or shallow volumes (rather than surfaces) in 2D and 3D parameter space, respectively (Figure 5—figure supplement 1). The width, depth, etc. are proportional to the tolerance. How does this impact multioutput solutions? Thin 2D strips intersect at a small patch (unlike curves intersecting at a point) but that patch in 2D parameter space (Figure 5—figure supplement 1A) is unlike the long 1D curve formed by broad 2D surfaces intersecting in 3D parameter space (see bottom panels of Figure 5). Likewise, shallow 3D volumes intersect at a narrow tube (unlike surfaces intersecting at a curve) but that tube in 3D parameter space (see Figure 5—figure supplement 1B) is unlike the broad 2D surface formed by deep 3D volumes intersecting in 4D parameter space. In other words, increasing tolerance does not increase solution space dimensionality the same way as increasing nin.

The need for ≥n adjustable channels to regulate n properties has been shown before (Olypher and Prinz, 2010; Olypher and Calabrese, 2007; Foster et al., 1993) and may be obvious to some for mathematical reasons, but the relationship has some important implications (e.g., for parameter estimation Sarkar and Sobie, 2010). The implications for homeostatic regulation have not been thoroughly explored, thus prompting the next steps of our study.

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 single- and multioutput solution sets found through homeostatic regulation. Figure 6A shows all combinations of g¯Na and g¯K (nin = 2) producing the target firing rate (nout = 1). Correlation between g¯Na and g¯K 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 g¯M is also allowed to vary (nin = 3), homeostatically determined solutions spread across a surface and pairwise correlations become predictably weaker (Figure 6B) since variation in one channel can be offset by variation in two other channels, and covariance is diluted across >1 pairwise relationship. If additional channels are allowed to vary (nin ≥ 4), solutions distribute over higher-dimensional manifolds and pairwise correlations further weaken (Figure 6C). However, if firing rate and energy efficiency are both regulated (nout = 2; Figure 6D), homeostatically determined solutions are once again constrained to fall along a curve when nin = 3, and pairwise correlations strengthen (Figure 6E). Increasing nin to four while keeping nout at two caused correlations to weaken (Figure 6F). These results confirm the predicted impact of solution space dimensionality on ion channel correlations.

Dimensionality of solution space affects ion channel correlations.

(A) Starting from a normally distributed cluster of g¯Na and g¯K (white dots) yielding an average firing rate of 73 spk/s, g¯Na and g¯K (nin = 2) were homeostatically adjusted to regulate firing rate (nout = 1) to its target value of 40 spk/s. Gray lines show trajectories. Because solutions (gray dots) converge on a curve, the pairwise correlation between g¯Na and g¯K is predictably strong. Scatter plots 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 panel A (nout = 1) but via control of g¯Na, g¯K, and g¯M (nin = 3). Homeostatically found solutions converge on a surface and ion channel correlations are thus weaker. (C) If g¯leak is also controlled (nin = 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 panel B (nin = 3), but for regulation of firing rate and energy efficiency (nout = 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 panel A. (F) If nin is increased to four while nout remains at 2, solutions converge on a surface (not shown) and ion channel correlations weaken.

O’Leary et al., 2013 demonstrated how the relative rates at which different ion channel densities are controlled impact their correlation. This is reproduced in Figure 7A, where, from the same initial conditions (density combinations), homeostatic control with different relative rates (shown in pink and cyan) produces solutions with different correlations. Relative regulation rates can affect not only the strength of pairwise correlations, but also the sign (compare correlation between g¯Na and g¯M). However, if firing rate and energy efficiency are both homeostatically regulated, correlations strengthen (consistent with results from Figure 6) and become independent of the relative regulation rates because solutions are limited to a lower-dimensional solution space (Figure 7B). Recall that dimensionality of the multioutput solution set corresponds to ninnout. For relative regulation rates to influence ion channel correlations, homeostatically determined solutions must fall on a solution manifold with dimensionality >1, but not >>1 lest pairwise correlations be diluted.

Effect of relative regulation rates on ion channel correlations depends on the dimensionality of solution space.

(A) Homeostatic regulation of firing rate via control of g¯Na, g¯K, and g¯M. Same as Figure 6B, but for two new sets of regulation rates (see Supplementary file 1). 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., 2013. (B) Same as panel A (nin = 3), but for homeostatic regulation of firing rate and energy efficiency (nout = 2). Solutions converge on a curve (yellow), giving rise to virtually identical ion channel correlations regardless of regulation rates.

However, correlations produced through the homeostatic regulation mechanism described by O’Leary et al., 2013, which is essentially the same mechanism used here, were recently shown to be sensitive to noise (Franci et al., 2020). This occurs because this regulation mechanism brings conductance densities to the solution manifold (by minimizing the error signal) but cannot control how solutions drift on the manifold (since the error signal is 0 everywhere on the manifold). Solutions do not drift in the absence of noise but, of course, noise is ubiquitous in biological systems and its effects must be considered. Accordingly, we reran simulations from Figures 6 and 7 with noise added to the conductance densities. After a few hundred noise-update iterations, solutions spread across the available solution space, which corresponds to a surface when only firing rate is regulated (Figure 8A, top). The distribution of solutions can produce correlations (Figure 8A, bottom). Lower and upper bounds on conductance densities shape the solution space and influence correlations (Figure 8—figure supplement 1A). Correlations also depend on how solutions distribute across the solution space, which depends on regulation rates (Figure 8—figure supplement 1B, C), which means correlations can still depend on regulation rates under noisy conditions, but correlations are not the same as in the absence of noise (Figure 8B). When firing rate and energy efficiency were regulated, noise caused solutions to spread across the intersection (Figure 8C) but correlations were relatively unaffected since they are limited by the low dimensionality of this solution space. Spread was bounded by one or another conductance density reaching 0 mS/cm2. Franci et al., 2020 proposed another regulatory scheme in which an attractive subspace emerges through cooperative molecular interactions; correlations created by fluctuations along that subspace are robust to noise. That scheme should preclude correlations from arising as in Figure 8A because solutions are prevented from spreading across the solution space; conversely, a low-dimensional solution space like in Figure 8C may prevent correlations unless the attractive subspace arising from molecular interactions aligns with the solution space. Further work is required to explore these regulatory mechanisms and how they interact.

Figure 8 with 1 supplement see all
Noise can affect ion channel correlations depending on the dimensionality of solution space.

(A) Starting from solution sets in Figure 6B, g¯Na, g¯K, and g¯M had noise added to them and were then updated by homeostatic regulation (using regulation rates from Figure 6) to correct the noise-induced disruption of firing rate. Conductance density combinations are depicted before and after 50, 200, 1000, or 3000 noise-update iterations (left to right). Noise caused solutions to spread across the surface (top). Some solutions (number in italics) drifted beyond the illustrated region but this is prevented by imposing an upper bound (Figure 8—figure supplement 1A). Ion channel correlations reflect the available solution space; for example, changes in g¯Na and g¯K are limited by g¯M remaining positive (dashed green line), while changes in g¯Na and g¯M are limited by g¯K remaining positive (dashed purple line). How solutions distribute within that space depends on regulation rates, and also affects correlations (Figure 8—figure supplement 1B, C). Conductance density distributions and pairwise correlations (bottom) were not centered on the mean or normalized by standard deviation (unlike in other figures) in order to visualize how distributions evolve over iterations. (B) Comparison of regression lines from Figure 6 (gray) and 7 (pink, cyan) without noise (top) and with noise (bottom). Correlations are affected by regulation rates in both conditions, but not in the same way. (C) Same as panel A but for regulation of firing rate and energy efficiency. Solutions spread along the intersection of the two surfaces, but the spread is bounded by one or another conductance density reaching 0 mS/cm2. Correlations slightly increase under noisy conditions (compare with Figure 6).

The dimensionality of solution space affects 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 9A shows an example in which coadjusting g¯Na, g¯K, and g¯M successfully regulates firing rate. Figure 9B shows successful regulation of energy efficiency by coadjusting the same three channels. Using the same initial conditions and relative regulation rates as above, the system failed to coregulate firing rate and energy efficiency (Figure 9C). Notably, coordinated regulation of both properties was achieved using other relative regulation rates (see Figures 6C and 7B) or using the same relative rates but starting from different initial conditions (not shown), suggesting that low-dimensional solutions are less accessible (i.e., a smaller set of relative regulation rates succeed in finding the solution space). This may be a challenge for the Franci et al. regulation model if molecular interactions create an attractive subspace that does not align with a low-dimensional solution space.

Low-dimensional solutions can be hard for homeostatic regulation to ‘find’.

(A) Homeostatic regulation of firing rate (nout = 1) via control of g¯Na, g¯K, and g¯M (nin = 3), like in Figures 6B and 7A but using a different set of regulation rates (see Supplementary file 1). 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 (nin = 3) but for homeostatic regulation of energy efficiency (nout = 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 (nout = 2) via control of the same ion channels (nin = 3) using the same relative rates and initial values as in panels A and B. Neither firing rate (red trajectories) nor energy efficiency (green trajectories) reached its target value. Conductance densities were capped at 4 mS/cm2 but this does not account for trajectories not reaching their target. Noise was not included in these simulations.

Figure 10 shows additional examples of regulating firing rate and energy efficiency by coadjusting g¯Na, g¯K, and g¯M. In these simulations, firing rate is regulated to a precise target value whereas energy efficiency is maintained above a lower bound; the single-output solution set for energy efficiency thus corresponds to a volume rather than a surface for nin = 3. For energy efficiency ≥22%, homeostatically determined solutions converge on the iso-firing rate surface without regulation of energy efficiency having much effect (Figure 10A). 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 being met, solutions then move across the surface until reaching the intersection with the iso-energy efficiency volume (Figure 10B). By converging on the intersection, which is a curve, ion channel correlations are stronger than in Figure 10A. In other words, ion channel correlations increase as single-output solution sets start to disconnect, meaning increased correlations presage regulation failure; in general, this is consistent with the impact of increasing constraints (see Figure 8), which go hand in hand with decreasing degrees of freedom. 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 single-output solution sets constrained by g¯K and g¯M reaching 0 mS/cm2 (Figure 10C). In this last example, neither property is optimally regulated but the outcome is a reasonable compromise. Regulation might have failed outright if, in the absence of an upper bound, ion channel densities increased (wound-up) without properties ever reaching their target values. This highlights that coregulation of multiple properties might fail or find suboptimal solutions not because single-output solutions do not exist, but because even large single-output solution sets might not intersect, meaning no multioutput solution exists.

The outcome of homeostatic regulation depends on if and how single-output solution sets intersect.

Homeostatic regulation of firing rate and energy efficiency (nout = 2) via control of g¯Na, g¯K, and g¯M (nin = 3). For these simulations, energy efficiency was maintained above a lower bound rather than being regulated to a specific target value; accordingly, the single-output solution set 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 single-output solution sets (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 g¯Na since g¯K and g¯M cannot become negative (bottom panel). Noise was not included in these simulations.

Discussion

This study has identified conditions required for a neuron to coregulate >1 property. Being able to produce the same output using diverse ion channel combinations ensures that single-output solution sets are large. This is critical because of the many ion channel combinations that produce the desired output for one property, few also produce the desired output for other properties (Figure 3), which means ion channel adjustments that regulate one property are liable to disrupt other properties due to ion channel pleiotropy (Figures 2 and 4). Indeed, the multioutput solution set corresponds to the intersection between single-output solution sets and the dimensionality of the multioutput solution corresponds to the difference between the number of adjustable ion channels (nin) and the number of regulated outputs (nout) (Figure 5). Coregulation of n properties requires at least n adjustable ion channels for a unique solution, and at least n + 1 channels for a degenerate solution. This constraint is not alleviated by pleiotropy, but nin would need to exceed nout by an even wider margin if ion channels were not pleiotropic. Moreover, channels must be coadjusted with different ratios to regulate different properties, which constrains how feedback loops are organized (see below). These important issues get overlooked if regulation of each property is considered in isolation.

With respect to the number of channels required to coregulate a certain number of properties, 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. Degeneracy is synonymous with underdetermination, which highlights that degeneracy depends not only on ion channel diversity, but also on those channels not having to coregulate ‘too many’ properties. Extra degrees of freedom broaden the range of available solutions, meaning a good solution is liable to exist over a broader range of conditions. On the other hand, if constraints outnumber the degrees of freedom, the system becomes overdetermined and solutions disappear (e.g., Figure 5D), which can cause regulation to fail (see below). One might reasonably speculate that ion channel diversity has been selected for because it enables the addition of new functionality (e.g., excitability) without the cell becoming overdetermined, lest pre-existing functions (e.g., osmoregulation) become compromised.

There are notable similarities and differences between a neuron adjusting its ion channel densities to regulate properties 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 (Foster et al., 1993), just as regulating more properties makes the solution set smaller. But how good are those parameter estimates? Assuming there is no measurement noise, 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 (Sarkar and Sobie, 2010; Aster et al., 2013). The neuron solves the forward problem, producing a firing rate, input resistance, etc. based on its channel density combination. That combination is determined by negative feedback but the precise densities are unimportant so long as the combination produces the target values for all regulated properties (which happens for every density combination in the multioutput solution set). As conditions change, the multioutput solution will evolve and negative feedback will adjust densities accordingly. A neuron needs to regulate its properties but does not do so by regulating its conductance densities to particular values; in that respect, a neuron does not solve an inverse problem. (Convergence of conductance densities to the same ‘attractive’ set regardless of initial conditions, or after a perturbation, may suggest otherwise [Franci et al., 2020], but may reflect other, unaccounted for constraints like regulation of another property.) Neurons have evolved under selective pressure to be degenerate (see above), not parsimonious, contrary to how most models are constructed.

Previous studies using grid searches to explore degeneracy have tended to apply selection criteria simultaneously, finding the ‘multioutput’ solution set in one fell swoop. In contrast, we considered one criterion (property) at a time in order to find each single-output solution set, which we then combined to find the multioutput solution. The former approach is akin to aggregating several error functions to create a single-objective problem, whereas the latter resembles multiobjective optimization. Nevertheless, past studies have observed that certain parameter changes (in a particular direction through parameter space) dramatically affect some properties but not others (Drion et al., 2015; Goldman et al., 2001), consistent with manifolds representing the single-output solution sets for sensitive and insensitive properties lying orthogonal to one another in parameter space. The effects of solution dimensionality on ion channel correlations and regulation failure (see below) highlight the value of a multiobjective perspective, which implicitly recognizes that there are multiple properties, each with its own feedback loop, even if those feedback loops intersect. Our results show that coregulating many properties using a common set of ion channels is feasible only if enough different channels are involved. Indeed, for coadjustments to work, channels must differ from each other in how they affect each property (see Figure 1), which emphasizes the distinction between degeneracy and redundancy.

Our results also highlight the need for >1 ‘master regulator’ because regulation of each property requires that ion channels are coadjusted with different ratios (Figure 11). O’Leary et al., 2014 argued in favor of a single master regulator, but that was predicated on using a single factor, namely global calcium, to encode the error signal; because a single factor cannot reach two different targets (i.e., simultaneously minimize two errors), it cannot support two master regulators. But local (Bootman et al., 2001) or kinetically distinct (Liu et al., 1998) variations in calcium could encode >1 error signal if the calcium sensors coupling intracellular calcium to gene expression or channel modulation are spatially segregated or operate with different filter properties. Beyond calcium, evidence points to firing rate homeostasis being activity independent (MacLean et al., 2003), using membrane potential as feedback (Santin and Schulz, 2019), or involving dual feedback loops (Kulik et al., 2019). Furthermore, the AMP:ATP ratio is used to track energy level (Hardie, 2014; Garcia and Shaw, 2017) and free amino acids can also be sensed (Efeyan et al., 2012). This is just the tip of the iceberg. Rather than error signals being encoded by a single factor like calcium, we favor the notion of a multi-input/multioutput system (Papin and Palsson, 2004) in which factors combine to encode multiple error signals. Take mTOR (mechanistic target of rapamycin) for example; it is modulated by diverse factors and, in turn, modulates many processes including transcription, translation and protein degradation (Switon et al., 2017). The regulation model proposed by Franci et al., 2020 provides additional insights. These issues require more investigation but, for now, our results argue that coregulation of multiple properties is incompatible with the amalgamation of error signals or certain control steps. The analogy with single- and multiobjective optimization (see above) seems apropos. Critically, nin exceeding nout will not guarantee solutions if degrees of freedom are reduced by bottlenecks elsewhere in the feedback loops.

Each error signal must be able to coadjust channels with different ratios.

Cartoons depict regulation of two properties (A and B) via control of three ion channels (X, Y, and Z). Channel levels are homeostatically adjusted to minimize the difference (error) between each property and its target value. The error signal divided by each channel’s regulation time constant (τ) dictates the magnitude of the change in that channel (represented by length of curved arrows in panel B). (A) The set of regulation time constants define a ‘master regulator’ (black dial) that coadjusts ion channels according to a certain ratio. Convergence of error signals at or before the master regulator limits how ion channels are coadjusted. (B) An additional master regulator (red dial) is needed to enable error B to coadjust channels with a different ratio than error A. Conditions like in panel B were required for all coregulation simulations (see Figure 6C). Different coadjustment ratios are required to regulate each property because of the different way each ion channel affects each property (see Figure 1). Note that master regulators may be more distributed than depicted here, with each regulation time constant τi,j likely reflecting the net effect of regulating multiple (transcriptional, translational, etc.) processes. The positive molecular regulatory network proposed by Franci et al., 2020 involves additional feedback mechanisms not easily depicted in this sort of cartoon.

Ion channel correlations have been studied in simulations (O’Leary et al., 2014; Mukunda and Narayanan, 2017; Taylor et al., 2009; Jain and Narayanan, 2020; O’Leary et al., 2013; Franci et al., 2020; Soofi et al., 2012; Hudson and Prinz, 2010; Ball et al., 2010) and experiments (Zhao and Golowasch, 2012; Temporal et al., 2014; Schulz et al., 2006; Tobin et al., 2009; Schulz et al., 2007; Khorkova and Golowasch, 2007). These correlations have been ascribed to the coregulation of ion channels. 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 (O’Leary et al., 2013). Our results are consistent with that explanation (Figure 7), notwithstanding noise effects (see below), but they also highlight the importance of manifold dimensionality: Correlations are stronger but are insensitive to relative regulation rates if trajectories reach a 1D manifold (curve) rather than a higher-dimensional manifold (surface, volume, etc.) (Figure 6). Under noisy conditions, correlations reflect the available solution space and how solutions drift across it, which depends on regulation rates (Figure 8). The cooperative molecular interactions proposed by Franci et al., 2020 can produce correlations in a high-dimensional solution space, but may be prevented from doing so by a low-dimensional solution space. Notwithstanding molecular interactions, the existence of pairwise correlations suggests that the dimensionality of solution space (=ninnout) is relatively low. This may be surprising since nin is high, but makes sense if nout is also high; in other words, there are many channels but they are responsible for regulating many properties. Though nout is constrained by nin (i.e., the number of regulated properties cannot exceed the number of adjustable ion channels), neurons might take advantage of regulating as many properties as their ion channel diversity safely allows, leading to a relatively low-dimensional solution space (as nout approaches nin).

Homeostatic regulation can fail for different reasons. If there are no solutions (i.e., the solution set is empty), negative feedback cannot regulate a property to its target values. A multioutput solution set can be empty because single-output solution sets do not intersect (Figure 10C), meaning regulation of different properties to their respective target values is incompatible. Regulation can also fail because negative feedback fails to converge on available solutions (Figure 9C). Notably, a multioutput solution set may become less accessible (lower dimensional) as single-output solution sets start to separate (Figure 10B), and may foreshadow the eventual disjunction of those solution sets (Figure 10C). 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 regulation rates are normally flexible but might fail to reach those solutions if regulation rate flexibility is reduced. Cooperative molecular interactions (Franci et al., 2020) are notable in this regard, insofar as they would tend to constrain regulation and thus influence failure modes. 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 (Ramocki and Zoghbi, 2008), but subtle differences in exactly how the failure transpires may provide important clues to help pinpoint the underlying mechanism(s).

In conclusion, neurons can coregulate multiple properties by coadjusting ion channels. Despite reusing channels to regulate more than one property, many channels are required to coregulate many properties and feedback loops must allow for coadjustments with different ratios. This may account for why evolution has yielded such diverse channels and why their transcriptional regulation has become so complicated (Gabel et al., 2015). Given how difficult it is to study regulation of any one property, studying the coregulation of multiple properties seems truly daunting, especially if feedback loops intermingle and error signals involve combinatorial codes. But if that is indeed the case, then accounting for the coregulation of multiple properties may be crucial for appreciating design features that might otherwise elude us.

Methods

Slice electrophysiology

All procedures were approved by the Hospital for Sick Children Animal Care Committee. Using the same procedures and equipment previously described (Khubieh et al., 2016), coronal slices of hippocampus were prepared from an adult mouse and a CA1 pyramidal neuron was recorded using whole-cell patch clamp. A junction potential correction of −9 mV was applied to the recorded membrane potential. A noisy stimulus was generated through an Ornstein–Uhlenbeck process (see Equation 7) with τstim = 5 ms, μstim = 60 pA, and σstim = 10 pA; the same stimulus was replayed on each trial. Native IM was blocked by bath application of 10 µM XE991 (Tocris). The block was continued throughout dynamic clamp experiments. Virtual IM and IAHP were modeled as per Equation 6 (with τz = 100 ms, βz = –35 mV, and γz = 4 mV for IM, and τz = 600 ms, βz = 0 mV, and γz = 1 mV for IAHP) and were applied using the dynamic clamp capabilities of Signal v6 (Cambridge Electronic Design). The density of each virtual conductance was adjusted to produce the desired firing rates (g¯M = 10 nS, g¯AHP = 120–150 nS).

Neuron model

Our base model includes a fast-activating sodium conductance (gfast) and a slower-activating potassium conductance (gslow) which together are sufficient to produce spikes. Other conductances were added to modulate excitability.

(1) CdVdt=Istimg¯fastm(V)(VENa)g¯sloww(VEK)g¯Nan(VENa)g¯Kn(VEK)g¯aHPp(VEK)g¯Mq(VEK)g¯leak(VEl),

where V is voltage and m changes instantaneously with V whereas other gating variables change more slowly. The spike-generating conductances gfast and gslow were modeled using a Morris–Lecar formalism (Rho and Prescott, 2012), with C = 2 µF/cm2, ENa = 50 mV, EK= −100 mV, Eleak = −70 mV, φw= 0.15, g¯fast = 20 mS/cm2, g¯slow = 20 mS/cm2, g¯leak = 2 mS/cm2, βm = −1.2 mV, γm = 18 mV, βw = −10 mV, and γw = 10 mV. A generic sodium conductance (gNa) and potassium conductance (gK) were modeled using Hodgkin–Huxley formalism as described by Ratté et al., 2014b

(2) INa,K= g¯Na,Kn(VENa,K),
(3) dndt=α(1n)βn,
(4) α=kα(VVα)e(VVαsα)1 ,
(5) β=kβe(VVαsα),

where Vα,β = −24 mV, sα,β = −17 mV, and kα,β = 1 ms−1. Note that gNa and gK differ only in their reversal potentials. Channels approximating a calcium-activated potassium conductance (gAHP) and an M-type potassium conductance (gM) were modeled as described by Prescott and Sejnowski, 2008

(6) dzdt={11+e(βzVγz)z}/τz,

where τz= 100 ms, γz = 4 mV, and βz = 0 mV or −35 mV for gAHP and gM, respectively. Maximal conductance densities for gleak, gNa, gK, and gM were systematically varied or adjusted by a homeostatic feedback mechanism (see below). Any other parameters that differ from the sources cited above are reported in the relevant figure legends. Injected current Istim was applied as either a constant step or as noisy fluctuations modeled with an Ornstein–Uhlenbeck process

(7) dIstimdt=Istimμstimτstim+SσstimdtN(t),

where τstim is a time constant that controls the rate at which Istim drifts back toward the mean μstim, and N (t) is a random number drawn from a normal distribution with 0 mean and unit variance that is scaled by Sσstim, where a scaling factor S = √2/τstim makes the standard deviation σstim independent of τstim. All simulations were conducted in MATLAB using the forward Euler integration method and a time step of 0.05–0.1 ms.

Energy calculations

Energy consumption rate was calculated as in Hasenstaub et al., 2010 Briefly, models were stimulated with a fast fluctuating stimulus (µstim = 40 µA/cm2, σstim = 10 μA/cm2). Sodium and potassium current through all channels was integrated for 1 s to determine the charge for each ion species, which was then converted to ion flux based on the elementary charge 1.602 × 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/cm2 s). Energy efficiency was calculated as the ratio between capacitive minimum and total Na+ flux during an action potential (Sengupta et al., 2010). 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.

Grid search

Models were tested with conductance density combinations chosen from a 100 × 100 or 30 × 30 × 30 grid for 2D and 3D plots, respectively. To depict each single-output solution set, 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Ω cm2 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 5—figure supplement 1 but were not shown in other figures for sake of clarity.

Feedback control

We used a homeostatic learning rule similar to O’Leary et al., 2014; O’Leary et al., 2013 with two notable differences: (1) we did not use intracellular calcium or other biological signals (e.g., AMP:ATP ratio) as intermediaries for our error signals and (2) the error for each output was determined as the difference between the current and target values at the end of each iteration, and conductance densities were adjusted before the start of the next iteration, rather than updating and feeding back error signals in ‘real time’. The separation between fast and slow (regulatory) timescales justifies the former approach. For each conductance g¯i, error was divided by the regulation time constant (τi) and added to the conductance (see Figure 4A). For coregulating >1 property, we used the sum of scaled errors to update conductance densities (see Figure 6D). A single run consists of a maximum of 200 iterations, during which a model must reach and maintain its regulated property within the tolerance for five consecutive iterations. All models that reached the target output(s) did so in well under 100 iterations; regulation was deemed to have failed for models not reaching their target output(s) within 200 iterations. Conductance densties during the last five iterations were averaged and reported as the final value. See Supplementary file 1 for the initial conductance densities and regulation time constants used for each figure.

Conductance noise

For simulations in Figure 8, noisy variations in each conductance density, g¯Na, g¯K, and g¯M, was applied by adding a random number drawn from a Gaussian distribution with a mean of 0 and a standard deviation of 0.05 mS/cm2. Noise was independent for each conductance. After applying noise, the simulation was run, errors were calculated, and conductance densities were updated according to the feedback control described above; this qualifies as one noise-update iteration. Up to 3000 noise-update iterations were run. If the addition of noise caused a conductance density to become negative, the density was reset to 0 mS/cm2 before applying feedback control. Where indicated, conductance density was likewise reset to 4 mS/cm2 if noise caused it to increase above that value.

Code availability

All computer code is available at http://modeldb.yale.edu/267309 and at http://prescottlab.ca/code-for-models.

Data availability

All computer code is available at http://modeldb.yale.edu/267309 and at http://prescottlab.ca/code-for-models. Key parameter values are provided in Supplementary file 1. Other parameter values are identified in the Methods. Source data are provided for Figure 2.

The following data sets were generated
    1. Yang J
    (2022) ModelDB
    ID 267309. A modified Morris-Lecar model with gM and gAHP.

References

    1. Goldman MS
    2. Golowasch J
    3. Marder E
    4. Abbott LF
    (2001)
    Global structure, robustness, and modulation of neuronal models
    The Journal of Neuroscience 21:5229–5238.

Decision letter

  1. Timothy O'Leary
    Reviewing Editor; University of Cambridge, United Kingdom
  2. John R Huguenard
    Senior Editor; Stanford University School of Medicine, United States
  3. Alessio Franci
    Reviewer; UNAM, Mexico

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Minimal requirements for a neuron to co-regulate many properties and the implications for ion channel correlations and robustness" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and John Huguenard as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Alessio Franci (Reviewer #1).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1) Return to the analysis and modeling to address whether the results may substantially alter depending on the regulation model assumed. In particular, take account of the updated mathematical model described in Franci et al., (2020) IEEE Control Systems Letters 4 (4), 946-951.

2) Address concerns about the experimental data. In particular, stronger statistical support is needed regarding reported changes in spike height and resting membrane potential and the reversal of these changes with dynamic clamp. In addition, please provide a better example neuron for Figure 2, which seems to show a very depolarized resting potential of approximately -30 mV after blocking I-m. This level of depolarization is inconsistent with any spiking at all.

3) Address the specific comments detailed in the reviewer assessments below.

Reviewer #1 (Recommendations for the authors):

I really appreciated the effort of this paper but the modeling assumptions underlying it hamper its biological relevance. My suggestion would be to redo the analysis when the regulation loops are not independent (in particular, when the target for one loop depends on the value of all the other variables) and when the revised model in Franci et al., 2020 is used for each loop.

Reviewer #2 (Recommendations for the authors):

I strongly suggest changing the way the solution surfaces are displayed in Figures 6 – 9. The hatching patterns make it difficult to see other elements of the figure, such as the dots and regulation trajectories. The hatching also does not add anything to the perception how the surface is situated in parameter space. For example, in Figure 5 supplement 1 the surfaces are easy to see without hatching.

https://doi.org/10.7554/eLife.72875.sa1

Author response

Essential revisions:

1) Return to the analysis and modeling to address whether the results may substantially alter depending on the regulation model assumed. In particular, take account of the updated mathematical model described in Franci et al., (2020) IEEE Control Systems Letters 4 (4), 946-951.

According to our understanding of Franci et al., (2020), that study was prompted in part by the observation that ion channel correlations emerging through the regulation mechanism used by O’Leary et al., (PNAS 2013, Neuron 2014) are disrupted by noise. Noise is ubiquitous and we totally agree that noise effects must be considered (see below). Other considerations like whether conductance densities converge on the same values regardless of initial conditions (i.e. whether the solution set is attractive) may influence the choice of regulation models in certain contexts (e.g. emergence of cell types during development; O’Leary et al., 2014) but we do not think this is critical for the regulation we studied; for instance, we are not aware of strong experimental evidence that conductance densities necessarily return to their original values after a perturbation. Thus, the most important question, it seems to us, is whether our conclusions are invalidated by noise.

To investigate this, we added noise to our model and explored the effects on ion channel correlations under different conditions (i.e. for solution spaces with different dimensionality). By adding noise to each conductance density and allowing feedback regulation to restore output(s) to target value(s), we show that solutions spread out across the solution manifold, consistent with Franci et al., (2020). This is now shown in Figure 8, which is an entirely new figure on the effects of noise. When the solution space is low-dimensional (when a submanifold is formed by intersecting manifolds), none of the conductance densities can rise very high without one of the other conductance densities hitting zero (Figure 8C) and the spread is thus bounded. As solutions spread across the submanifold, correlations strengthen but are not qualitatively altered. For higher-dimensional solution spaces, solutions drift across the manifold and sometimes proceed outside the illustrated region (Figure 8A). But saturation of rate-limiting production steps almost certainly occurs, and would limit how high conductance densities can rise (lines 657-659); the effect of including an upper bound is shown in Figure 8—figure supplement 1. These bounds prevent infinite spread. Noise-induced spreading of solutions on a higher-dimensional solution manifold destroys the correlations formed according to the O’Leary et al., (2013) mechanism, but other correlations emerge, as explained in Figure 8 and its supplement. These results are exciting insofar as they affirm our main conclusion about the importance of solution space dimensionality: effects of noise on ion channel correlations depend on the dimensionality of the solution space.

Our new results also demonstrate that our simple regulation mechanism can successfully co-regulate multiple properties under noisy conditions, and can do so without feedback loops being coupled. This is all that we asked of the regulation mechanism insofar as we did not set out to explain the emergence of ion channel correlations or other such phenomena; instead, we focused on whether such outcomes differed depending on the number of regulated properties and the number of adjustable ion channels. That said, outcomes might differ depending on the regulation mechanism (see below). But we do not think it is fair to impose additional requirements on the regulation mechanism that are not essential for our focus (e.g. that the solution set is attractive). Whether a neuron regulates its properties or its conductance densities is addressed at some length in our Discussion (lines 296-315).

We strove to implement the cooperative molecular interactions proposed by Franci et al., but ran into difficulties because we do not fully understand the mathematical explanations in the 2020 paper. Apart from these technical difficulties, we also have reservations about certain suggestions. For example, coupling the feedback loops in the proposed manner requires use of a common feedback signal, which is something we believe does not happen (see Lines 342-344) and is unnecessary from our perspective (see above). Our own unpublished experimental data make us skeptical about the premise for cooperative interactions; specifically, we have found that different sodium channel types are negatively correlated and that this is not reflected in transcript levels because it depends on the regulation of translation, not transcription. There is clearly much more work to do in this area, but insofar as our model provides valuable insights using a feedback mechanism that does not have obvious shortcoming (including under noisy conditions), exploring alternative feedback mechanisms arguably exceeds the scope of the current study.

Though we did not get cooperative molecular interactions working in our model, we have added discussion on how anticipated results would differ from the results obtained with our simple regulation mechanism. In a high-dimensional solution space, we expect that cooperative interactions would preclude correlations arising from the mechanism shown in Figure 8A since the Franci mechanism prevents solution from drifting (explained on lines 226-227) and would, instead, produce other correlations. On the other hand, in a low-dimensional solution space, the constraints of the solution space (Figure 8C) may prevent correlations unless the attractive subspace emerging from molecular interactions align with the solution space (lines 227-229). This is addressed also in the Discussion (lines 362-364). We have also touched upon how other issues addressed by Franci et al., (2020) apply to our results, including attractiveness of the solution set (lines 312-314) and how cooperative molecular interactions might influence success of regulation (lines 242-243; 381-383).

2) Address concerns about the experimental data. In particular, stronger statistical support is needed regarding reported changes in spike height and resting membrane potential and the reversal of these changes with dynamic clamp. In addition, please provide a better example neuron for Figure 2, which seems to show a very depolarized resting potential of approximately -30 mV after blocking I-m. This level of depolarization is inconsistent with any spiking at all.

We have analyzed the experimental data as requested and have incorporated this analysis into Figure 2. Statistical analysis has been reported in detail in the figure legend (lines 511-524). With respect to membrane potential, a junction potential correction had not been applied to data reported in the original submission; applying the -9 mV correction (lines 400-401) renders the values more accurate. Abrupt blockade of Im is expected to cause strong depolarization which will in turn cause some sodium channel inactivation and attenuated spike height. That inactivation is relieved by the hyperpolarization produced by either of the two virtual conductances, despite continued presence of the drug, 4-AP. In short, the chosen traces exemplify what we are trying to illustrate.

We want to clarify here that this experiment was conducted in a single neuron. Certain comments suggest that this was not clear, and we have revised the text to avoid potential confusion. We did not repeat the experiment in multiple neurons because our goal was not to quantify the impact of the compensatory mechanism (which depends on the strength of the virtual conductance we insert); rather, we sought simply to show that different cell properties are differently affected by each compensatory mechanism (as proposed in Figure 1), and that such differences are not an artifact of our model neuron, which is why we conducted this initial test experimentally rather than with simulations. Insofar as we would apply the same virtual conductance if we were to repeat the experiment in additional neurons, we would not expect notable differences. Trials were repeated to assess spike timing across trials rather than with statistical analysis in mind, but effects are evidently very reproducible across trials.

3) Address the specific comments detailed in the reviewer assessments below.

We have responded to all specific comments in the reviewer assessments below.

Reviewer #1 (Recommendations for the authors):

I really appreciated the effort of this paper but the modeling assumptions underlying it hamper its biological relevance. My suggestion would be to redo the analysis when the regulation loops are not independent (in particular, when the target for one loop depends on the value of all the other variables) and when the revised model in Franci et al., 2020 is used for each loop.

Thank you for the detailed suggestions. We hope the responses we provided above have addressed many of the concerns. We will address additional issues here, including our attempts to introduce cooperative molecular interactions and coupled feedback loops into our model.

Without solid experimental data for many aspects of our model, we tried to minimize assumptions by keeping our model simple (see lines 79-90). This extends to using minimalist feedback loops (i.e. not implicating specific biophysical signals or specific transcriptional/ translational/ post-translational processes). This limits the scope of our conclusions but, on the other hand, this makes us more confident in the conclusions we have drawn and provides a solid foundation for future work. We think it is reasonable to start by assuming that feedback loops are independent and we believe that experimental data support their independence (see above), and so it is valuable to know whether such an arrangement can produce ion channel correlations (with or without noise). This is now addressed in Figure 8.

As explained under point 1 of the essential revisions, we cannot follow all the mathematical arguments presented by Franci et al., (2020). This complicated our efforts to implement cooperative molecular interactions in our model. Combined with our reservations about the impetus for implementing certain changes (as expressed above), we stopped short of testing a fundamentally different regulation mechanism as part of the current study. As evident from a number of elegant studies motivated by and published since the 2014 O’Leary paper in Neuron, we believe that the simple regulation mechanism they proposed can lead to many valuable insights. Our study continues in that vein. We are intrigued by the regulation mechanism proposed by Franci et al., and we think it requires further investigation in connection with some of the concepts raised in our study, but we also think certain concerns about unrobustness must be tempered based on the focus of our study. We welcome deeper discussions with the reviewers, unfettered the peer review process, to push this work forward collaboratively. We also think there is an opportunity for Franci to write a piece to accompany our article, especially if he can implement his regulation mechanism in our model and show simulation results that challenge or support our conclusions.

Reviewer #2 (Recommendations for the authors):

I strongly suggest changing the way the solution surfaces are displayed in Figures 6 – 9. The hatching patterns make it difficult to see other elements of the figure, such as the dots and regulation trajectories. The hatching also does not add anything to the perception how the surface is situated in parameter space. For example, in Figure 5 supplement 1 the surfaces are easy to see without hatching.

We have lightened the hatching on all solution surfaces. When we removed the hatching, we realized that the hatching helped give some perspective. That perspective can be achieved with light hatching that does not obscure the solutions or their trajectories. We think our changes serve to optimize visualization.

https://doi.org/10.7554/eLife.72875.sa2

Article and author information

Author details

  1. Jane Yang

    1. Neurosciences and Mental Health, The Hospital for Sick Children, Toronto, Canada
    2. Institute of Biomedical Engineering, University of Toronto, Toronto, Canada
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Visualization, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0114-5503
  2. Husain Shakil

    1. Neurosciences and Mental Health, The Hospital for Sick Children, Toronto, Canada
    2. Department of Physiology, University of Toronto, Toronto, Canada
    Present address
    Division of Neurosurgery, Department of Surgery, University of Toronto, Toronto, Canada
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3995-6811
  3. Stéphanie Ratté

    Neurosciences and Mental Health, The Hospital for Sick Children, Toronto, Canada
    Contribution
    Data curation, Formal analysis, Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7005-6726
  4. Steven A Prescott

    1. Neurosciences and Mental Health, The Hospital for Sick Children, Toronto, Canada
    2. Department of Physiology, University of Toronto, Toronto, Canada
    Contribution
    Conceptualization, Funding acquisition, Methodology, Project administration, Supervision, Writing – original draft, Writing – review and editing
    For correspondence
    steve.prescott@sickkids.ca
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3827-4512

Funding

Canadian Institutes of Health Research (Foundation Grant 167276)

  • Steven Alec Prescott

Natural Sciences and Engineering Research Council of Canada (Discovery Grant RGPIN 436168)

  • Steven Alec Prescott

The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.

Acknowledgements

This research was funded by a Discovery Grant from the Natural Sciences and Engineering Research Council (NSERC) and by a Foundation Grant from the Canadian Institutes of Health Research (CIHR). H.S. was supported by Ontario Graduate Scholarship, the University of Toronto Centre for the Study of Pain Scholarship, and James F Crothers Family Scholarship. J.Y. was supported by a SickKids Restracomp Studentship and Ontario Graduate Scholarship. We thank Eve Marder, Etay Hay, Shreejoy Tripathy, and Simon Hardy for constructive feedback on the manuscript.

Ethics

All experimental procedures were approved by The Hospital for Sick Children Animal Care Committee (protocol #53451) and were conducted in accordance with guidelines from the Canadian Council on Animal Care.

Senior Editor

  1. John R Huguenard, Stanford University School of Medicine, United States

Reviewing Editor

  1. Timothy O'Leary, University of Cambridge, United Kingdom

Reviewer

  1. Alessio Franci, UNAM, Mexico

Publication history

  1. Preprint posted: December 4, 2020 (view preprint)
  2. Received: August 6, 2021
  3. Accepted: March 3, 2022
  4. Accepted Manuscript published: March 16, 2022 (version 1)
  5. Version of Record published: April 6, 2022 (version 2)

Copyright

© 2022, Yang et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 685
    Page views
  • 114
    Downloads
  • 2
    Citations

Article citation count generated by polling the highest count across the following sources: PubMed Central, Crossref, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Jane Yang
  2. Husain Shakil
  3. Stéphanie Ratté
  4. Steven A Prescott
(2022)
Minimal requirements for a neuron to coregulate many properties and the implications for ion channel correlations and robustness
eLife 11:e72875.
https://doi.org/10.7554/eLife.72875
  1. Further reading

Further reading

    1. Developmental Biology
    2. Neuroscience
    Anadika R Prasad, Inês Lago-Baldaia ... Vilaiwan M Fernandes
    Research Article Updated

    Neural circuit formation and function require that diverse neurons are specified in appropriate numbers. Known strategies for controlling neuronal numbers involve regulating either cell proliferation or survival. We used the Drosophila visual system to probe how neuronal numbers are set. Photoreceptors from the eye-disc induce their target field, the lamina, such that for every unit eye there is a corresponding lamina unit (column). Although each column initially contains ~6 post-mitotic lamina precursors, only 5 differentiate into neurons, called L1-L5; the ‘extra’ precursor, which is invariantly positioned above the L5 neuron in each column, undergoes apoptosis. Here, we showed that a glial population called the outer chiasm giant glia (xgO), which resides below the lamina, secretes multiple ligands to induce L5 differentiation in response to epidermal growth factor (EGF) from photoreceptors. By forcing neuronal differentiation in the lamina, we uncovered that though fated to die, the ‘extra’ precursor is specified as an L5. Therefore, two precursors are specified as L5s but only one differentiates during normal development. We found that the row of precursors nearest to xgO differentiate into L5s and, in turn, antagonise differentiation signalling to prevent the ‘extra’ precursors from differentiating, resulting in their death. Thus, an intricate interplay of glial signals and feedback from differentiating neurons defines an invariant and stereotyped pattern of neuronal differentiation and programmed cell death to ensure that lamina columns each contain exactly one L5 neuron.

    1. Neuroscience
    Brittany J Bush, Caroline Donnay ... J Christopher Ehlen
    Research Article

    Resilience, the ability to overcome stressful conditions, is found in most mammals and varies significantly among individuals. A lack of resilience can lead to the development of neuropsychiatric and sleep disorders, often within the same individual. Despite extensive research into the brain mechanisms causing maladaptive behavioral-responses to stress, it is not clear why some individuals exhibit resilience. To examine if sleep has a determinative role in maladaptive behavioral-response to social stress, we investigated individual variations in resilience using a social-defeat model for male mice. Our results reveal a direct, causal relationship between sleep amount and resilience-demonstrating that sleep increases after social-defeat stress only occur in resilient mice. Further, we found that within the prefrontal cortex, a regulator of maladaptive responses to stress, pre-existing differences in sleep regulation predict resilience. Overall, these results demonstrate that increased NREM sleep, mediated cortically, is an active response to social-defeat stress that plays a determinative role in promoting resilience. They also show that differences in resilience are strongly correlated with inter-individual variability in sleep regulation.