Abstract

Neurons can maintain stable activity profiles over their lifetimes despite ion channel turnover over minutes to weeks. Neuronal activity is also influenced by regulating the voltage dependence of ion channels. How do these two forms of plasticity work together to maintain a stable activity profile? We augment a classical model of activity-dependent ion channel density regulation with a mechanism that adjusts channel voltage dependence based on activity. These findings reveal that the timescale of these mechanisms shapes the specific electrical activity patterns that achieve a target activity profile. Moreover, alterations in channel voltage dependence can impact a neuron’s ability to recover from perturbations. These results highlight a potentially distinct role for activity-dependent regulation of channel voltage dependence in maintaining stable neuronal activity profiles.

Introduction

Neurons maintain stable activity patterns over years despite continual turnover in their ion channel composition, which occurs on timescales from days to weeks. Experimental studies have shown that neurons homeostatically regulate both channel number (Desai et al., 1999; Haedo and Golowasch, 2006; Mizrahi et al., 2001; Santin and Schulz, 2019; Thoby-Brisson and Simmers, 2002; Turrigiano et al., 1994; Turrigiano et al., 1995; Viteri and Schulz, 2023) and channel voltage dependence (Gasselin et al., 2015; Haedo and Golowasch, 2006; Thoby-Brisson and Simmers, 2002) to stabilize their activity.

These mechanisms have been widely studied through computational modeling, but most of this work has focused on changes in channel number (Abbott and LeMasson, 1993; Alonso et al., 2023; Golowasch et al., 1999; Gorur-Shandilya et al., 2020; LeMasson et al., 1993; O’Leary et al., 2014; Siegel et al., 1994; Srikanth and Narayanan, 2015). This is largely due to the plethora of work indicating such changes occur (Desai et al., 1999; Haedo and Golowasch, 2006; Mizrahi et al., 2001; Santin and Schulz, 2019; Thoby-Brisson and Simmers, 2002; Turrigiano et al., 1994; Turrigiano et al., 1995; Viteri and Schulz, 2023). However, this emphasis also reflects the challenges associated with precisely measuring changes in channel voltage dependence. Inter-neuronal variability in activation and inactivation curves can obscure clear patterns of activity-mediated modifications, especially when analyzed at the population level (Turrigiano et al., 1995). As a result, our understanding of how modifications to channel voltage dependence contribute to activity stability remains limited.

This gap is important to address given the ubiquity of channel voltage-dependence regulation. Neurons modulate voltage dependence through mechanisms ranging from rapid post-translational modifications (Fraser and Scott, 1999; Johnson et al., 1994) or slower processes like subunit synthesis and degradation (Wu et al., 2023). These processes are known to be associated with calcium-mediated second messengers. Intracellular calcium is also thought to play an important role in maintaining neuronal activity set points (LeMasson et al., 1993; Turrigiano et al., 1994). As such, calcium-mediated homeostasis may operate, at least in part, by tuning channel voltage dependence—yet the consequences of this form of regulation are largely uncharacterized.

To address this, we extended a state-of-the-art computational homeostatic model. This model stabilizes neuronal activity by monitoring intracellular Ca2+ concentrations and inserting and deleting ion channels (Alonso et al., 2023). To this model, we added a mechanism that alters channel voltage dependence. Importantly, this model enabled independent control over the timescales of channel number and channel voltage-dependence alterations. This organization of the model let us simulate fast changes in voltage dependence, such as those driven by phosphorylation, or slower changes, like those resulting from synthesis or degradation of channel subunits. Our central finding is that the timescale of voltage-dependence regulation influences the specific channel configurations a neuron adopts—both while attaining a prespecified activity target and while maintaining it during perturbation.

Results

The model

To explore how modifications in channel voltage-dependence influence a neuron’s capacity to achieve an activity profile as changes in ion channel density occur, we created a computational model. This model was built around a single compartment, Hodgkin–Huxley type conductance-based neuron with seven intrinsic currents (see Turrigiano et al., 1995 or Appendix 1, Neuron model). The neuron model consisted of a fast Na+ current, two Ca2+ currents, three K+ currents, a hyperpolarization-activated inward current (IH), and a leak current. Each current was described by its maximal conductance (number of ion channels) and activation curves. Some currents also possessed inactivation curves as well. The (in)activation curves described the effects of voltage on the opening and closing of the membrane channels.

The activity-dependent regulation of channel number and (in)activation curves employed a strategy that uses intracellular Ca2+ concentrations as a proxy for the neuron’s activity (Abbott and LeMasson, 1993; Alonso et al., 2023; Golowasch et al., 1999; Gorur-Shandilya et al., 2020; LeMasson et al., 1993; O’Leary et al., 2014; Siegel et al., 1994; Srikanth and Narayanan, 2015). Specifically, we followed the approach used in Liu et al., 1998 and Alonso et al., 2023, for altering channel density to guide our approach to altering (in)activation curves. In those models, the intracellular Ca2+ concentration and three filters (also called sensors) of the Ca2+ concentration were computed. The three sensors were employed to better differentiate between the activity of tonic spiking and bursting neurons (Alonso et al., 2023; Liu et al., 1998). A target activity profile was chosen by prespecifying the target values of these sensors. These values were selected through trial and error to consistently produce regular bursting across a specific range of initial conditions (see Methods). When deviations from these Ca2+ sensor targets were detected, channel densities were altered to up- or down-regulate the number of channels in the membrane (Alonso et al., 2023; Liu et al., 1998). Details on how the sensors were computed and how the target activity—defined by the target sensor values—was specified are provided in the Methods.

We expanded on these models by incorporating a mechanism that shifts the (in)activation curves in response to deviations from the Ca2+ sensor targets. This was accomplished by changing the half-(in)activation voltages of the (in)activation curves, mirroring strategies employed in Liu et al., 1998 and Alonso et al., 2023. The details of the Model, with their equations and assumptions, are described in the Methods section. A brief overview is provided in Appendix 1.

Importantly, the response times of the mechanisms responsible for ion channel insertion and deletion, as well as those for shifts in (in)activation curves, were independently controlled. In the model, these response times were initially set two orders of magnitude apart. This reflected the slower processes of protein synthesis and degradation associated with channel insertion and deletion. This contrasted with faster timescales of post-translational modifications such as channel phosphorylation. Specifically, the timescale for changes in maximal conductance was set at τg=600 s, and the timescale for half-(in)activation shifts at τhalf=6 s. Later, we adjusted the response times of the (in)activation curve shifts. In prior computational studies of activity-dependent regulation, τg was set to be less than 10 s, substantially faster than the timescales used in our model (Alonso et al., 2023; Liu et al., 1998). This study lengthens this timescale and explores the impact of varying the relative timescales of maximal conductance changes (τg) and half-(in)activation shifts (τhalf).

Bursters assembled by activity-dependent alterations in ion channel density and channel voltage dependence

Figure 1 shows the full model, with the two regulation mechanisms in operation. From a set of starting randomly chosen parameters (maximal conductances and half-(in)activations), referred to as Starting Parameters 1 (SP1), the activity-dependent regulation mechanism altered the ion channel density and voltage dependences to assemble an electrical activity pattern that satisfied the Ca2+ targets. The targets were chosen so that a burster self-assembled (see Methods). The final ion channel density and voltage dependences are shown in the bar graphs (Figure 1A). The first two rows in Figure 1A show the electrical activity of the neuron as the model seeks to satisfy the Ca2+ targets. As it does, the model assembled a burster. The second row shows the activity over the entire trace, and the first row shows blows up at: the start (a), about 100 min (b), and at the end of the simulation (c). At the time shown in (b), the neuron was bursting, but the Ca2+ targets had not been achieved. There were two indications that this burster did not satisfy the specified Ca2+ targets: (1) the model still continued to adjust the properties of the ion channel repertoire here and (2) α was not yet zero.

Activity-dependent regulation mechanism (τhalf = 6 s and τg = 600 s) was used to assemble a burster that meets the Ca2+ targets from ‘random’ starting parameters (i.e. SP).

Both this set of starting parameters and its resultant model are called SP1. (A) Row 1: Blow-ups of specific points a, b, and c in the voltage trace in Row 2. Row 2: The neuron’s electrical activity as its channel density and (in)activations curves were being adjusted. α is a measure of how close the model is to satisfying Ca2+ targets (see Methods). It took values between 0 and 1, with 0 meaning the targets were satisfied. Row 3: Plots the adjustments made to the half-activations as the model attempted to satisfy the calcium targets. They are color-coded to match the intrinsic currents shown in the left plot of Row 5. Row 4: Same as Row 3, except for (in)activation curves. Note that not all currents in this model have inactivation curves. Row 5: Same as Row 3, except for maximal conductances (i.e. channel density). Row 6: The final levels of maximal conductances are shown in the left plot. The final half-(in)activations are plotted on the right. The final half-activation of an intrinsic current, Vmi,1/2 , was the sum of the initially posited level of the that current’s half-activation Vm^6,1/20 and the shifts the model made to the activation curve, Vsim (shown in Row 3). The final half-inactivation of an intrinsic current, Vhi,1/2 , was the sum of the initially posited level of the that current’s half-(in)activation Vh^0i,1/2 and the shifts the model made to the inactivation curve, Vsih (shown in Row 4). [The process is illustrated in Figure 1—figure supplement 1 for a different set of starting parameters (i.e. initial half-(in)activation shifts and maximal conductance levels).] (B) This displays the neuron’s electrical activity if the model was restricted to varying only half-(in)activations (τhalf = 6 s), with maximal conductances fixed (τg = ∞ s). The same starting parameters as in Panel A were used for this simulation.

Briefly, α measures whether the model has successfully tuned the neuron’s ion channel densities and half-(in)activations to create an activity pattern that matches the preset sensor targets; values near zero indicate successful tuning (more in Methods). Even after the model assembles a channel repertoire that satisfies the targets, brief excursions of α away from zero can occur. These excursions exist because the neuron’s electrical activity is quasistable—meaning it may lose/add a spike or experience a slight drop in slow wave amplitude, temporarily causing the sensors to be unsatisfied. When this happens, the model quickly makes minor adjustments to the channel densities and half-(in)activations to restore alignment with the targets (Figure 1—figure supplement 1, Row 2, between minutes 300 and 400). These adjustments are minimal (Figure 1, Rows 3–5).

The third and fourth traces in Figure 1A show how the model altered half-(in)activation and half-(in)activation values to meet the Ca2+ target. Very early on, there were large changes in the activation curves, which then slowly moved to a steady-state value. The fifth trace shows the evolution of maximal conductances. These started to change early on but did so slowly. As the maximal conductances slowly changed, the half-(in)activation curves followed suit, adapting to the new repertoire of ion channels to close the distance to the Ca2+ target. Figure 1—figure supplement 1 shows the assembly of a burster from a different set of starting parameters (maximal conductances and half-(in)activations (SP4)), with qualitatively similar results to those obtained with SP1 but quantitatively different final parameters.

The results shown in Figure 1A require activity-dependent regulation of the maximal conductances. When activity-dependent regulation of the maximal conductances is turned off, the model failed to assemble SP1 into a burster (Figure 1B). This was seen in the other 19 starting parameters (SP2–SP20), as well.

In general, when the model was initialized with different starting parameters, it produced bursters with different final intrinsic properties. Twenty degenerate bursters (SP1–SP20), conforming to standards specified in Methods, were selected from over 100 bursters produced with different starting parameters. Figure 2A shows the activity characteristics of these 20 bursters (period, burst duration, spike height, maximum hyperpolarization, and slow wave amplitude). The left panel of Figure 2B shows the range of conductances in this group. The right panel illustrates the variation in half-(in)activation curves. The asterisks indicate the baseline half-(in)activation values; all neuron models were assembled from starting parameters that lay within ±0.5 mV of these baseline values (see Methods for details).

All 20 SPs were analyzed, including those illustrated in Figure 1.

(A) The top left panel displays the periods for all SPs. Burst duration, interburst interval, spike height, maximum hyperpolarization, and slow wave amplitude measurements for these SPs are provided in subsequent plots. (B) B1 shows maximal conductances, and B2 displays half-(in)activation voltages for each of the 20 SPs. An asterisk on the right indicates the initially specified level of half-(in)activation from which random deviations were made.

Half-(in)activation mechanism timescale impacts activity pattern

In the previous section, we studied the types of bursters that arose from slowly changing the density of the channel repertoire and quickly altering the half-(in)activations of the resulting repertoire. We next changed the timescale over which we altered half-(in)activations and assessed how the burst characteristics of the population changed (the timescale of channel density alterations remains the same). Figure 3A shows that alterations in the timescale of half-(in)activation alterations significantly impacted the period, interburst interval, spike height, and maximum hyperpolarization of the entire bursting population.

The speed at which the model changed half-(in)activations in response to deviations from Ca2+ target impacted the type of bursters the model assembled.

In all simulations for this figure, τg = 600. (A) These plots extend the measurements from Figure 2 to the same 20 SPs, now with half-(in)activations regulated at a different speed. In blue, τhalf = 6 s; recapitulating results in Figure 2. In green and orange, the SPs were reassessed with slower half-(in)activation adjustments, τhalf = 600 s, and τhalf = ∞ s (i.e. halted), respectively. Brackets indicate significant differences in medians (p < 0.05). These pairwise differences were assessed using a Kruskal–Wallis test (to assess whether any differences in medians existed between all groups) and post hoc Dunn tests with Bonferroni corrections to assess which pairs of groups differed significantly in their medians. (B) The model was used to stabilize SP1 to the same Ca2+ target but using different timescales of half-(in)activation alterations. τhalf = 6 s is shown on the left and τhalf = ∞ s (representing the slowest possible response) on the right. Rows 1–4 below each voltage trace show: the total outward current, percentage contribution of each outward intrinsic current to the total outward current, percentage contribution of each inward intrinsic current to the total inward current, and the total inward current. A dashed vertical line across all rows marks the point of maximum hyperpolarization in each activity pattern to guide focused comparison. (C) When the model stabilized around a burster, the speed at which the model changed half-(in)activations impacted maximal conductance and half-activation location of the calcium-activated potassium current. Shown here in blue, green, and orange are τhalf = 6 s, τhalf = 600 s, and τhalf = ∞ s (i.e. halted), respectively. SP1–SP20 are in the colored circles. Trends are illustrated with a black line, with large circles marking median positions. The mauve dot highlights the levels to which maximal conductance and half-(in)activation of SP1 evolved. A Kruskal–Wallis test followed by post hoc Dunn tests (with Bonferroni corrections) was conducted to assess significant (p < 0.05) pairwise differences in medians (brackets).

The timescale of half-(in)activation alterations also impacted some intrinsic properties of the population. Figure 3B plots the burster assembled with the timescale of half-(in)activation alterations at two extremes: fast and infinitely slow (i.e. off). This visualization is known as a currentscape (Alonso and Marder, 2019). In both cases, the resulting burster appeared to terminate using the same mechanism: hyperpolarization via the calcium-activated potassium current. This is seen by the relative contribution the calcium-activated potassium current makes at the point of maximum hyperpolarization (see at the dotted line in Figure 3B). However, the manner in which this current came to dominate the burster termination differs. Note that the burster constructed without half-(in)activation alterations possessed a higher maximal conductance for IKCa than the burster constructed with half-(in)activation alterations (Figure 3C, mauve dot). Moreover, the latter burster had more depolarized half-(in)activations (Figure 3C, mauve dot). In fact, this observation held true for SP2–SP20 as well (Figure 3C).

Potential mechanism

To explore why the properties of the resulting bursters depend on the timescale of half-(in)activation adjustments, we examined what happens when SP1 is assembled under different half-(in)activation timescales: (1) fast, (2) intermediate (matching the timescale of ion channel density changes), and (3) infinitely slow (i.e. effectively turned off). The effects of these timescales can be seen by comparing the zoomed-in views of the SP1 activity profiles under each condition (Figure 4).

Changes to the timescale of half-(in)activation alterations steer the model through different trajectories in intrinsic parameter space.

Each row shows the behavior of Starting Parameter 1 (SP1) as the model adjusts toward a configuration that satisfies the calcium activity target, using different timescales for half-(in)activation regulation. (A) τhalf = 6 s, (B) τhalf = 600 s, and (C) τhalf = ∞ s (i.e. half-(in)activations do not change). In all panels, the speed at which the model changed channel density is fixed in all rows: τg = 600 s. Superimposed on these traces is the α parameter, which indicates whether the activity patterns were meeting the Ca2+ targets. Blow-ups of the activity pattern when the calcium sensors were satisfied are displayed to the right.

When half-(in)activations are fast, the time evolution of α—which tracks how far the activity pattern is from its targets (see Methods)—shows an abrupt jump as it searches for a voltage-dependence configuration that meets calcium targets (Figure 4A). As this happens, the channel densities are slightly altered, and this process continues again. Slowing the half-(in)activations alterations reduces these abrupt fluctuations (Figure 4B). Making the alterations infinitely slow effectively removes half-(in)activation changes altogether, leaving the system reliant solely on slower alterations in maximal conductances (Figure 4C). Because each timescale of half-(in)activation produces a different channel repertoire at each time step, different timescales of half-(in)activation alteration led the model through a different path in the space of activity profiles and intrinsic properties. Ultimately, this resulted in distinct final activity patterns, all of which were consistent with the Ca2+ targets.

Half-(in)activation alterations contribute to activity homeostasis

This model also simulates how the assembled bursters recover their activity following perturbations. In particular, we modeled an increase in extracellular potassium concentration by altering both the potassium reversal potential and the leak reversal potential (see Methods). For these simulations, we used the previously tuned models SP1–SP20. They begin with the homeostatic mechanism engaged but idling because these models already met the prespecified activity targets. However, once extracellular potassium concentration increases, the shifts in reversal potentials cause all 20 models to experience depolarization block. In response, the homeostatic mechanism becomes operative and adjusts the neuron’s maximal conductances and half-(in)activations.

We anticipated that the perturbation would depolarize the membrane. To accommodate this shift, we increased the allowable deviation for the DC calcium sensor, ΔD (see Methods), from 0.015 to 0.035. Figure 5A shows SP1’s activity before perturbation (blue) and after 265 min of recovery (red), under three conditions: fast half-(in)activation regulation, intermediate (same timescale as conductance regulation), and no half-(in)activation regulation. The depolarized membrane potential is evident in the more positive maximum hyperpolarizations during perturbation. Notably, although all three cases eventually re-stabilized to activity patterns that satisfied the calcium targets, the resulting burst patterns were different.

The speed at which the model changed half-(in)activations in response to deviations from Ca2+ targets impacted the type of bursters the model assembled following perturbation.

(A) Voltage traces from SP1 before (blue) and 265 min after (red) a simulated increase in extracellular potassium concentration. Each panel shows recovery under a different timescale for half-(in)activation regulation: fast (τhalf = 6 s), intermediate (τhalf = 600 s), and infinitely slow (i.e. off, τhalf = ∞ s). In all conditions, alteration of maximal conductances was held fixed (τg = 600 s). (B) Summary of activity features for all SP1–SP20 models that recovered to a stable burst pattern. The number of models that successfully recovered varied with the homeostatic mechanism: fast (n = 17), intermediate (n = 18), and no half-(in)activation modulation (n = 20). Gray points indicate pre-perturbation measurements of SP1–SP20 (from Figure 3A). Significant pairwise differences in median values were assessed using a Kruskal–Wallis test followed by post hoc Dunn tests with Bonferroni correction. Significant comparisons (p < 0.05) are marked with horizontal brackets. (C) The speed of half-(in)activation regulation influenced both the maximal conductance and half-activation voltage of KCa. Data are grouped by as they are in B. Statistical comparisons were performed as in B.

SP1 was one of eleven models that successfully recovered under all three homeostatic conditions. Eight models required some form of half-(in)activation modulation to stabilize; in these cases, conductance regulation alone was not sufficient for recovery. Among these eight, two only recovered when half-(in)activation changes were fast, and another two only recovered when they were slow. Evidently, for some models, the ability to recover bursting activity depended on both the presence and speed of half-(in)activation modulation.

Figure 5B, C quantifies the variability in post-recovery activity patterns and intrinsic parameters across different half-(in)activation regulation timescales. Comparing the fast (τhalf=6 s, blue) and intermediate (τhalf=600 s, green) conditions reveals significant differences in burst period, interburst interval, and spike height (Figure 5B). Notably, these same features also showed significant variation in Figure 3A, where timescale-dependent effects were observed during the initial assembly of bursters. Similarly, the half-activation voltage of the KCa current differs significantly between these two timescales (Figure 5C), paralleling the differences observed in Figure 3C. In models that successfully recovered from perturbation, the timescale of voltage-dependence modulation influenced both the resulting activity patterns and the intrinsic parameters of the stabilized bursters.

Discussion

Timescale of channel voltage-dependence alterations

While alterations in channel voltage dependence are known to modulate a neuron’s excitability (Levitan, 1988), their role in stabilizing a neuron around an activity target has not been extensively explored. This study explores how modifying the timescale of voltage-dependence regulation influences both the initial attainment of a target activity pattern and its maintenance during perturbation. While all models met the activity target regardless of whether voltage-dependence modulation was present, recovery after perturbation depended on both the presence and timescale of this mechanism. Furthermore, in both scenarios, the timescale shaped the stabilized neuron’s specific activity characteristics and underlying intrinsic parameters. This suggests a general role for the timescale of half-(in)activation modulation in sculpting stable activity states.

Figure 6 provides a conceptual summary of this idea. The space of possible intrinsic parameter configurations (Figure 6A) and electrical activity patterns (Figure 6B) each contain regions that satisfy a given activity target—depicted in blue for activity characteristics and green for intrinsic parameters. The perturbations administered in this study simply shift these target-satisfying regions. When channel gating properties are altered quickly in response to deviations from the target activity, the neuron ultimately settles into a stable activity pattern. The resulting electrical patterns are shown in Figure 6 as the orange bubble labeled ‘τhalf=6 s’. At a slower timescale of alterations, the activity characteristics and model parameters of the electrical patterns shift, as illustrated by the other orange subregions in Figure 6. This shift conceptually highlights how changing the timescale of alterations to channel voltage dependence can move the neuron through different regions within the space of activity characteristics or intrinsic parameters while still maintaining its activity target.

The impact of adjusting the timescale of ion channel (in)activation curve alterations can be illustrated conceptually using the space of activity characteristics (right) and underlying intrinsic parameters (left).

The regions in green and blue are the activity patterns that are consistent with the Ca2+ targets in the space of activity characteristics and underlying intrinsic parameters, respectively. As the rate of ion channel half-(in)activation adjustments in response to deviations from Ca2+ targets is changed, the regions targeted by the model also shift (shown in orange). Regardless of the timescale, all targeted regions reside within a larger region encompassing all activity profile measurements and intrinsic parameters consistent with the Ca2+ targets.

Our simulations show that altering the timescale of the channels’ voltage-dependent gating properties impacts its activity characteristics. However, what biological processes can implement these alterations? Fast timescale post-translational alterations can be implemented by second messengers that are anchored next to ion channels (Fraser and Scott, 1999; Johnson et al., 1994) or second messengers that diffuse across the cytosol (Agarwal et al., 2016; Zaccolo et al., 2006). Phosphorylation by kinases is a commonly studied post-translational modification; however, there are other modifications that also occur on these timescales, such as S-nitrosylation (Broillet, 1999). On slower timescales, the neuron can modulate the transcription and translation of auxiliary subunits based on intracellular calcium levels. For instance, Potassium Channel-Interacting Proteins are known to be controlled by CREB (Wu et al., 2023).

Biological relevance

One application for the simulations involving the self-assembly of activity may be to model the initial phases of neural development, when a neuron transitions from having little or no electrical activity to possessing it (Baccaglini and Spitzer, 1977). As shown in Figure 6, the timescale of (in)activation curve alterations defines a neuron’s activity characteristics and intrinsic properties. As such, neurons may actively adjust these timescales to achieve a specific electrical activity aligned with a developmental phase’s activity targets. Indeed, developmental phases are marked by changes in ion channel density and voltage dependence, leading to distinct electrical activity at each stage (Baccaglini and Spitzer, 1977; Gao and Ziskind-Conhaim, 1998; Goldberg et al., 2011; Hunsberger and Mynlieff, 2020; McCormick and Prince, 1987; Moody and Bosma, 2005; O’Leary et al., 2014; Picken Bahrey and Moody, 2003).

Additionally, our results show that activity-dependent regulation of channel voltage dependence can play a critical role in restoring neuronal activity during perturbations (Figure 5). Specifically, the presence and timing of half-(in)activation modulation influenced whether the model neuron could successfully return to its target activity pattern. Many model neurons only achieved recovery when a half-(in)activation mechanism was present. Moreover, the speed of this modulation shaped recovery outcomes in nuanced ways: some model neurons reached their targets only when voltage dependence was adjusted rapidly, while others did so only when these changes occurred slowly. These observations all suggest that impairments in a neuron’s ability to modulate the voltage dependence of its channels may lead to disruptions in activity-dependent homeostasis. This may have implications for conditions such as addiction (Kourrich et al., 2015) and Alzheimer’s disease (Styr and Slutsky, 2018), where disruptions in homeostatic processes are thought to contribute to pathogenesis.

In conclusion, our findings suggest that when a neuron begins from a state of no activity and small randomly specified channel densities, changes in channel density are essential for attaining activity. In this context, activity-dependent changes in channel voltage dependence alone did not assemble bursting from these low-conductance initial states (Figure 1B). Based on this evidence alone, one might conclude that shifts in (in)activation curves have negligible impact, effectively being overshadowed by the larger changes in excitability driven by alterations in channel density. This perspective aligns with a common assumption in systems analysis: that short- and long-timescale processes are sufficiently distinct and can be studied independently without significant interaction. However, our results demonstrate that fast shifts in half-(in)activation constrain the types of activity characteristics and intrinsic properties a neuron ultimately expresses (Figure 6). These fast alterations, while subtle, manifest over long timescales required to assemble activity. This study exemplifies how interactions between processes operating on different timescales can influence one another. Such principles are particularly relevant for understanding multiscale interactions in neuroscience, such as the connections between subcellular, cellular, and network-level properties (Bhalla, 2014; Marom, 2010).

Methods

Details of the model and its implementation

Appendix 1 presents a scheme for the equations used in the model: a Hodgkin–Huxley type equation with seven intrinsic currents, Ca2+ sensors, and an activity-dependent regulation mechanism.

The Hodgkin–Huxley type neuron consisted of seven intrinsic currents (Equations A1–A3) and an intracellular Ca2+ concentration calculation (Equation A4). The neuron’s capacitance (C) was set to 1 nF. The currents were labeled by the index i and were in order: fast sodium (INa), transient calcium (ICaT), slow calcium (ICaS), hyperpolarization-activated inward cation current (IH), potassium rectifier (IKd), calcium-dependent potassium (IKCa), and fast transient potassium (IA). Additionally, there was a leak current (IL). The maximal conductances were given by g¯i. Activation and inactivation dynamics were given by mi and hi, respectively. The exponents qi were given by 3,3,3,1,4,4,3, respectively. The equilibrium/reversal potentials of the sodium current, hyperpolarization-activated cation current, potassium currents, and leak current were +30, –20, –80, and –50 mV, respectively. The calcium equilibrium potential was computed via the Nernst equation (Equation A5) using the computed internal Ca2+ concentration ([Ca2+]i) and an external Ca2+ concentration of 3 mM (temperature was 10°C).

We modeled perturbations in neuronal activity by a 2.5-fold increase in extracellular potassium concentration. Using the Nernst equation, this shift corresponds to a change in the potassium reversal potential from −80 to −55 mV. We assume that leak current in our model includes a potassium component, so that the perturbation also impacts the leak reversal potential. Additionally, assuming that the remaining leak current comes from sodium channels, we can compute how to perturb the leak reversal potential. The leak maximal conductance for our model is fixed at 0.01 μS. The new leak reversal potential was calculated using standard computation involving the conductance-weighted average of potassium and sodium reversal potentials, yielding a value of approximately –32 mV.

The activation curves and inactivation curves were given by mi,∞ and hi,∞. Their associated time constants were τmi and τhi, respectively. The (in)activation curves and their associated time constants were functions of the membrane voltage, V (see tables in Appendix 1, Neuron model). Note that IKCa had an activation curve that depends on internal Ca2+ concentration. The internal calcium concentration was computed using ICaT and ICaS (Equation A4).

There were three Ca2+ sensors that responded to different measures associated with fluctuations of [Ca2+]i during a burst. They were used to determine how close these measures are to their targets. These measures were the amount of [Ca2+]i that entered the neuron: as a result of the slow wave of a burst (S), as a result of the spiking activity of a burst (F), and on average during a burst (D). They were computed using Equations A6–A8 (Appendix 1, Ca2+ Sensors). Each sensor had a corresponding target value—S¯,F¯,D¯—which represented the target activity pattern.

The model computed moving averages of the mismatch between each sensor and its target, denoted, ES, EF, and ED (Equations A9–A11). These averages were taken over a 2-s window (τS=2000 ms), allowing the model to smooth out short-term fluctuations and capture sustained deviations from the target. Ideally, ES, EF, and ED, are zero; however, ongoing membrane fluctuations, such as bursting or spiking, still created small variations in these measures. Allowable tolerances are specified in the model as, Δs, ΔF, and ΔD, respectively. The moving averages of the mismatches between the calcium sensors and their targets are compared against the allowable tolerances and combined in Equation A12 (Appendix 1, Ca2+ Sensors) to create a ‘match score’, SF (Equation A12).

To compute this match score, we adapted a formulation from Alonso et al., 2023, who originally used a root-mean-square (RMS) or L2 norm to combine the sensor mismatches. In that approach, each error (ES, EF, and ED) is divided by its allowable tolerance (Δs, ΔF, and ΔD) to produce a normalized error. These normalized errors are then squared, summed, and square-rooted to produce a single scalar score that reflects how well the model matches the target activity pattern.

In our version, we instead used an L8 norm, which raises each normalized error to the 8th power before summing and taking the 1/8th root. This formulation emphasizes large deviations in any one sensor, making it easier to pinpoint which feature of the activity is limiting convergence. By amplifying outlier mismatches, this approach provided a clearer view of which sensor was driving model mismatch, helping us both interpret failure modes and tune the model’s sensitivity by adjusting the tolerances for individual sensor errors.

Although the L8 norm emphasizes large deviations more strongly than the L2 norm, the choice of norm does not fundamentally alter which models can converge—a model that performs well under one norm can also be made to perform well under another by adjusting the allowable tolerances. The biophysical mechanisms by which neurons detect deviations from target activity and convert them into changes in ion channel properties are still not well understood. Given this uncertainty, and the fact that using different norms ultimately should not affect the convergence of a given model, the use of different norms to combine sensor errors is consistent with the broader basic premise of the model: that intrinsic homeostatic regulation is calcium mediated.

Equations A13 and A14 were used to evaluate whether the match score was sufficient for the model to cease altering maximal conductances and (in)activation curves, indicating that the target had been achieved. In Equation A14 (Appendix 1, Ca2+ sensors), the ‘match score’ was compared to a preset threshold value the model would accept, ρ. This equation returned 1 if the match score is below ρ, a 0 if the match score is above ρ, and some number in between 0 and 1 if the match score was close to the threshold, ρ. Equation A13 (Appendix 1, Ca2+ sensors) then effectively computed a 2-s time average of Equation A14. Equation A13 was an activity sensor, α, that was plotted in Figure 1. It took the value zero when all three Ca2+ sensors were close to the Ca2+ targets but remained close to one otherwise. The value of this activity sensor was used to slow the evolution of maximal conductances and half-(in)activations as the Ca2+ sensors reached their targets. It also served as a measure for how well the model satisfied the Ca2+ sensors. In this model, ρ=0.3. The ‘sharpness’ of the threshold at ρ=0.3 was set by Δα. We set Δα=0.01 indicating that when the match score, Sf, was greater than ~0.32, the model considered the calcium targets satisfied.

Using the information regarding how well the Ca2+ sensors were satisfied, the model regulated conductance densities and (in)activation curves. The maximal conductances g¯i, and the ion channel voltage sensitives, as measured by activation curve shifts, Vsim, and inactivation curve shifts, Vsih, from their specified locations (Vm^0i,1/2 for activation curves and Vh^0i,1/2 for inactivation curves), were altered by the model in response to deviations from the calcium targets (Equations A15–A17). The latter were adjusted on a timescale given by τhalf and the former by τg. The index i ran from 1 through 7 for the maximal conductances (Equation A15) and half-(in)activations (Equation A16). In Equation A17, the index only took on a subset of values: i=1,2,3,7—because only some intrinsic currents in this model had inactivation (the associated currents being listed above).

Equations A16 and A17 were new equations added to the model. They were built similarly to Equation A15. In Equation A15, the maximal conductances were altered when the calcium sensors didn’t match their targets. There are two parts to Equation A15: a part that converts calcium sensor mismatches into changes in maximal conductances and a cubic term. The cubic term prevents unbounded growth in the maximal conductances. We set the constant associated with the cubic term, γ, to 10-7. The other term adjusted the maximal conductance based on fluctuations of calcium sensors from their target levels. The parameters Ai,Bi, and Ci reveal how the sensor fluctuations contribute to the changes of each maximal conductance. These values were provided in Liu et al., 1998. Importantly, note that this term is multiplied by g¯i, which prevents the maximal conductance from becoming negative and scales the speed of evolution so that smaller values evolve more slowly and larger values evolve more quickly. As such, multiplication by g¯i effectively induced exponential changes in the maximal conductance, allowing g¯i to vary over orders of magnitude.

Equations A16 and A17, introduced in this paper, enable the adjustment of the ion channel voltage dependence. Equations A16 and A17 also each have two components: a component that converts calcium sensor mismatches into changes in half-(in)activations and a cubic term. The cubic term prevents unbounded translations. The cubic term’s constant, γ^, was set to 10-5. The other term adjusts the half-activations (Equation A16) and half-inactivation (Equation A17) based on deviation of calcium sensors from their target levels.

Finally, we discuss how the coefficients Ai,Bi, and Ci in Equation A15, and the corresponding parameters in Equations A16 and A17 were chosen. The same coefficients used to adjust maximal conductances were used for the inactivation curve shifts, A^^i, B^^i, and C^^i. The coefficients of the activation curve shift are their negatives, A^i, B^i, and C^i. To understand this choice, consider a scenario where the neuron is bursting more rapidly than desired. In that case, then F>F¯, S>S¯, and D>D¯. To reduce the excitability of the neuron, we made the depolarizing current, for example, INa, harder to activate, by shifting the activation curve of INa to a more depolarizing potential. This implies dVs1mdt>0. Therefore, the coefficients A^i, B^i, and C^i must be non-positive. The coefficients for INa in Liu et al., 1998 were given as A1=1,B1=0, and C1=0 and were chosen so that dg¯1dt<0, thereby decreasing the maximal conductance of sodium in the same scenario. So, we flipped the signs to give A^1=−1,B^1=0, and C^1=0. In this way, only the F Ca2+ sensor contributed to the modulation of INa, but the activation curves were now adjusted to promote homeostasis. Furthermore, the model makes it harder to de-inactivate INa by shifting its inactivation curve to more hyperpolarized values (dVs1hdt<0). In the scenario in which the activity sensors are greater than their activity targets, this corresponds to setting A^^1=0,B^^1=0, and C^^1=0.

A similar logic can be used to understand how the activation curves of the outward currents were adjusted. In the scenario described above, the neuron was made less excitable by increasing the amount of an outward current, for example IA. This was done by increasing its maximal conductance or moving its activation curve to more hyperpolarized values. Liu et al., 1998 assigned the coefficients A7=0,B7=−1, and C7=-1 to adjust the maximal conductance of IA. In this scenario where all sensors are above their targets, this implies dg¯7dt>0. Following the prescription we outlined above, flipping the signs gives A^7=0,B^7=1, and C^7=1. This, in turn, implied dVs7mdt<0—creating hyperpolarizing shifts in IA’s activation curve when the activity is greater than its target, as we expect. In a similar spirit, we assign the coefficients A^^7=0, B^^7=−1, and C^^7=−1, so that IA’s inactivation curve shifts to more depolarizing potentials (dVs7hdt>0) to make it easier to de-inactivate.

Note, the coefficients in Liu et al., 1998 do not perfectly correspond with this reasoning. In Liu et al., 1998, some coefficients were assigned values to ensure convergence of the model. Still, we followed the prescription outlined above (Appendix 1—table 7).

Starting parameters

The model has 40 ordinary differential equations: 13 equations describe the neuron’s electrical activity, 9 equations describe calcium sensor components, and 18 equations alter the maximal conductances and half-(in)activations (Appendix 1). This means we needed to specify 40 starting values for the model. We defined starting parameters (SPs) as the bursters the model assembles from a particular set of randomly chosen starting maximal conductances and half-(in)activations. The maximal conductances were randomly selected between 0.3 and 0.9 nS and half-(in)-activation shifts are randomly offset between –0.5 and +0.5 mV from their specified half-(in)activation values (Vm^i,1/20 and Vh^i,1/20 in Appendix 1—table 2). Other starting values were V=-50 mV and the gating variables of activation and inactivation were randomly selected between 0.2 and 0.3. The starting value of intracellular calcium concentration was set as [Ca2+]i=0.4 µM. Equations A9–A11 (sensor errors) and Equation A13 (activity sensor, α) were initialized with the value 1. HX (Equation A7) and MX (Equation A6) were initialized to 0 and 1, respectively.

The random initialization of gating variables was not expected to significantly impact the formation of an activity profile that satisfied the activity target. Bursters may exhibit bistability if they are initialized with different gating variables. However, this is for a fixed set of maximal conductances and half-(in)activations. Any bistability in the initially specified profile was likely lost as the model altered the maximal conductances and half-(in)activationss.

The model was integrated using the Tsitouras 5 implementation of Runga–Kutta 4 (5) in Julia’s DifferentialEquations.jl package (Tsitouras, 2011). The simulation was saved every 0.5 ms regardless of step size (step size of the integrator was adaptive and left unspecified). The maximum number of iterations was set to 1010 and relative tolerance was set to 0.0001.

Burster measurements

We defined regular bursters as activity patterns that exhibite