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.

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.

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

  • 1,220
    views
  • 208
    downloads
  • 34
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

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

Share this article

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

Further reading

    1. Neuroscience
    Yiting Li, Wenqu Yin ... Baoming Li
    Research Article

    Time estimation is an essential prerequisite underlying various cognitive functions. Previous studies identified ‘sequential firing’ and ‘activity ramps’ as the primary neuron activity patterns in the medial frontal cortex (mPFC) that could convey information regarding time. However, the relationship between these patterns and the timing behavior has not been fully understood. In this study, we utilized in vivo calcium imaging of mPFC in rats performing a timing task. We observed cells that showed selective activation at trial start, end, or during the timing interval. By aligning long-term time-lapse datasets, we discovered that sequential patterns of time coding were stable over weeks, while cells coding for trial start or end showed constant dynamism. Furthermore, with a novel behavior design that allowed the animal to determine individual trial interval, we were able to demonstrate that real-time adjustment in the sequence procession speed closely tracked the trial-to-trial interval variations. And errors in the rats’ timing behavior can be primarily attributed to the premature ending of the time sequence. Together, our data suggest that sequential activity maybe a stable neural substrate that represents time under physiological conditions. Furthermore, our results imply the existence of a unique cell type in the mPFC that participates in the time-related sequences. Future characterization of this cell type could provide important insights in the neural mechanism of timing and related cognitive functions.

    1. Neuroscience
    Bhanu Shrestha, Jiun Sang ... Youngseok Lee
    Research Article

    Sour taste, which is elicited by low pH, may serve to help animals distinguish appetitive from potentially harmful food sources. In all species studied to date, the attractiveness of oral acids is contingent on concentration. Many carboxylic acids are attractive at ecologically relevant concentrations but become aversive beyond some maximal concentration. Recent work found that Drosophila ionotropic receptors IR25a and IR76b expressed by sweet-responsive gustatory receptor neurons (GRNs) in the labellum, a peripheral gustatory organ, mediate appetitive feeding behaviors toward dilute carboxylic acids. Here, we disclose the existence of pharyngeal sensors in Drosophila melanogaster that detect ingested carboxylic acids and are also involved in the appetitive responses to carboxylic acids. These pharyngeal sensors rely on IR51b, IR94a, and IR94h, together with IR25a and IR76b, to drive responses to carboxylic acids. We then demonstrate that optogenetic activation of either Ir94a+ or Ir94h+ GRNs promotes an appetitive feeding response, confirming their contributions to appetitive feeding behavior. Our discovery of internal pharyngeal sour taste receptors opens up new avenues for investigating the internal sensation of tastants in insects.