Introduction

Many neurons maintain stable patterns of activity over years, despite the fact that the ion channels that give rise to this activity are replaced thousands of times in their membranes, as they turn over on time scales from minutes to weeks. Indeed, a number of computational (Abbott & LeMasson 1993, Alonso et al 2023, Golowasch et al 1999b, LeMasson et al 1993, Liu et al 1998, O’Leary et al 2014) and experimental (Desai et al 1999, Haedo & Golowasch 2006, Mizrahi et al 2001, Santin & Schulz 2019, Thoby-Brisson & Simmers 2002, Turrigiano et al 1994, Turrigiano et al 1995, Viteri & Schulz 2023) studies have argued that neurons accomplish this stability in activity by prespecifying an activity target and achieving it through activity-dependent regulation of ion channel densities.

On the other hand, the gating properties of these ion channels are known to be regulated by second messenger pathways, enabling multiple discharge patterns to arise from the same channel repertoire. The neuron implements these changes by altering the channel’s voltage-dependent properties, viz. their activation curves, inactivation curves, and time constants, through conformational changes induced by post-translational modifications such as phosphorylation or S-nitrosylation (Broillet 1999, Hille 2001, Levitan 1988). These modifications are usually thought to take place on time scales that range over seconds to minutes but can include more long-lasting effects (Levitan 1988). The second messenger modulation of ion channels can have profound effects on excitability and can radically transform the patterns of electrical discharge of single neurons (Levitan 1987).

There is evidence that the neuron can also alter the gating properties of ion channels to achieve a stable activity profile (Gasselin et al 2015, Haedo & Golowasch 2006, Thoby-Brisson & Simmers 2002). This suggests that the second messenger pathways responsible for modifying gating properties may play a role in stabilizing neuronal activity. However, the contribution of changes in channel gating to achieving activity stability remains underexplored, largely due to the challenges associated with measuring such changes. 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).

To evaluate how regulation of channel voltage-dependence affects a neuron’s ability to achieve a stable activity profile amidst the insertion and deletion of ion channels, we developed a computational model. The model was built around a simulated neuron with seven intrinsic currents and a specified target activity profile. The model tracked fluctuations in intracellular Ca²⁺ concentration resulting from the neuron’s activity, with the target profile defined by specific Ca²⁺ concentration targets. The model then evaluated whether the Ca²⁺ fluctuations resulting from the neuron’s discharge pattern met these predefined targets. If they did not, the model iteratively adjusted the density (via alterations in maximal conductance) and voltage-dependence (via alterations in half-(in)activations) of all seven intrinsic currents until the target was achieved.

This approach is justified for two reasons. First, intracellular Ca²⁺ levels reflect neuronal activity (Grienberger & Konnerth 2012). Second, elevated Ca²⁺ concentration levels are known to be associated with the activation of many second messenger pathways that putatively regulate channel density (Northcutt & Schulz 2020) or regulate channel voltage-dependence (Kennedy 1989). However, a limitation of this approach is the assumption that the same Ca²⁺ target deviations driving changes in channel density also drive changes in channel voltage-dependence. We know, also, that not all changes to intrinsic excitability are mediated by activity induced Ca²⁺ concentration fluctuations (Khorkova & Golowasch 2007, MacLean et al 2005, Maclean et al 2003, Zhang et al 2009).

Still, this model has two important features that make it particularly useful. First, it allows independent control over the timescales at which channel voltage-dependence and channel density respond to deviations from the target activity profile. This enabled us to explore the specific role of voltage-dependence regulation in achieving neuronal stability amidst ongoing ion channel turnover by altering the relative timescale of their responses to target deviations. Second, this type of model is known to have a particular biologically realistic feature: it can achieve a given an activity target using different configurations of intrinsic currents. Indeed, experiments show that neurons with putative activity targets, such as constitutive bursters, exhibit varying channel densities while exhibiting similar electrical behaviors (Alonso et al 2023, Golowasch 2014, Golowasch et al 1999a, Liu et al 1998, Marder & Goaillard 2006, O’Leary et al 2014, Schulz et al 2006).

These two aspects of the model converge to highlight the central finding of this paper. The speed (or timescale) of these adjustments can guide the neuron toward subsets of intrinsic current configurations, all consistent with the predefined activity target. This suggests a potential role for second messenger-mediated regulation of channel voltage-dependence in shaping the electrical patterns that emerge as stability is achieved.

Results

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 Box 1, Supplement 1). The neuron model consisted of a fast Na⁺ current, two Ca²⁺ currents, three K⁺ currents, a hypolarization-activated inward current (I_h), 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 Ca²⁺ concentrations as a proxy for the neuron’s activity (Abbott & LeMasson 1993, Alonso et al 2023, Golowasch et al 1999b, Gorur-Shandilya et al 2020, LeMasson et al 1993, O’Leary et al 2014, Siegel et al 1994, Srikanth & 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 Ca²⁺ concentration and three filters (also called sensors) of the Ca²⁺ 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 average value of these filters. When deviations from these Ca²⁺ sensor target 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).

We expanded on these models by incorporating a mechanism that shifts the (in)activation curves in response to deviations from the Ca²⁺ 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 are described in the Methods section and in Box 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, could be 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 time scale 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. We note that in prior computational studies of activity-dependent regulation, τ_g is set to be less than 10 seconds, substantially faster than the timescales used in our model (Alonso et al 2023, Liu et al 1998). In contrast, our study lengthens this timescale 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 that 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 Ca²⁺ 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 graph (Figure 1A). The first two rows in Figure 1A show the electrical activity of the neuron as model seeks to satisfy the Ca²⁺ target. 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 minutes (b) and at the end of the simulation (c). At the time shown in (b), the neuron was bursting, but the Ca²⁺ targets had not been achieved. There were two indications that this burster did not satisfy the specified Ca²⁺ targets: (1) the model still continued to adjust the properties of the ion channel repertoire here and (2) ⍺ was not yet zero (see Methods).

Figure 1

Figure 1 (Supplement)

The third and fourth traces in Figure 1A show how the model altered half-activation and half-inactivation values to meet the Ca²⁺ 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-activation curves followed suit, adapting to the new repertoire of ion channels to close the distance to Ca²⁺ target. Figure 1 Supplement 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.

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. Five sets of Starting Parameters (SPs 1-5) were selected based on their periods being within 2%, with an average period of 416 ms. This group was then expanded to include an additional 15 model instances (SP6-SP20), whose periods are within 20% of the mean.

Figure 1B shows that without the activity-dependent regulation of the maximal conductances, the model failed to assemble SP1 into a burster. This was seen in the other 19 Starting Parameters (SP2-SP20), as well. We note, however, that for specific, suitably organized sets of maximal conductances, the half-activation regulation mechanism was able to assemble a burster. For instance, when we took the randomly initialized starting half-(in)activations and the finally assembled maximal conductances for SP1 the resulting activity pattern did not satisfy the activity target (Figure 1C, Blow-Up 1). However, by only altering the half-(in)activations, the model assembled a burster from the previous activity pattern (Figure 1C, Blow-Up 3).

Statistics of the Population

The burst characteristics of these 20 bursters were measured (Figure 2A). The blue circles represent the wider Group of 20 while the black dots denote the Group of 5. The stringency of period selection is seen in how the black dots are more concentrated within the more dispersed blue circles in the plot of Period in Figure 2A. In contrast, the black dots are more dispersed amongst the blue dots in plots of the Burst Duration, Spike Height, Maximum Hyperpolarization and Slow Wave Amplitude. This suggests that the narrow constraints on period did not constrain other characteristics of the burster. Indeed, bootstrap statistical analysis revealed that the relative standard deviation measurements between the Group of 5 and Group of 20 were not significant for any of the burst characteristics (see Figure 2 legend and Methods).

Figure 2

Figure 2B assesses how the stringency of period selection affected intrinsic currents. The left panel shows the range of conductances, and the right panel illustrates the variation in half-(in)activation curves. The asterisks represent, within 0.5 mV, where all the half-(in)activations were initialized (for more see Methods). It initially appeared that the impact of period stringency on parameter variability might reveal itself in some of these properties. In all intrinsic current properties, the black dots had a smaller relative standard deviation than the circles (Figure 2 legend). For example, g_Na for Group of 5 and Group of 20 had a relative standard deviation of 11.4% and 129.2% (Figure 2-B1), and the relative standard deviation of the half-activations for the sodium current were 2.5% vs 10.0%, respectively (Figure 2-B2). However, bootstrap analysis revealed that these differences were statistically insignificant. Evidently, the constraints on the Group of 5’s period were not reflected in the intrinsic current’s properties. Together, with the findings on burst characteristics, this suggested that the Group of 5 is quite diverse despite the strict period selection.

Half-(in)activation Mechanism Timescale Impacts Activity Pattern

In previous the 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). To start, we present an example of what occurred when SP1 was assembled under various timescales of half-(in)activation alterations: fast, with the same time scale of ion channel density alterations, and infinitely slow (i.e. off). Evidently, the timescale of half-(in)activation alterations did not impact whether a burster assembled. Rather, it impacted the burster characteristics of the assembled electrical pattern. This can be seen by comparing the blow-ups of the SP1 activity profiles in the three conditions in Figure 3. Notably, burst characteristics such as maximum hyperpolarization and period varied between the different homeostatic mechanisms originating from SP1. In fact, these differences were observed the entire SP1-SP20 population. Figure 4A 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.

Figure 3

Figure 4

The timescale of half-(in)activation alterations also impacted some intrinsic properties of the population. Figure 4B 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 & 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 4B). 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 higher maximal conductance for I_KCa than the burster constructed with half-(in)activation alterations (Figure 4C, Mauve dot). Moreover, the latter burster had more depolarized half-activations (Figure 4C, Mauve dot). In fact, this observation held true for SP2-SP20 as well (Figure 4C).

Half-(in)activation Alterations Contribute Activity Homeostasis

This model can also simulate the homeostatic recovery of activity in the assembled bursters following perturbations. In Figure 5A, SP1 received a depolarizing perturbation of extracellular potassium concentration increase for 100 minutes (in olive green) after which the model neuron was returned to normal conditions. The initially bursting neuron (Figure 5A, Blow-Up 1) received a perturbation that results in the cessation of activity (Figure 5A, Blow-Up 2). However, during the perturbation, the neuron altered its channel density and voltage-dependences to reattain bursting activity pattern (Figure 5A, Blow-Up 3). This bursting pattern, in fact, satisfied the calcium target (alpha is near zero in Figure 5A, Blue Arrow). Following the perturbation, the neuron briefly ceased activity (Figure 5A, Blow-Up 4), but returned to a bursting pattern that satisfied the calcium target (Figure 5A, Blow-Ups 5 and 6).

Figure 5

The model might offer some perspectives on the role that alterations in channel voltage-dependence have during homeostasis. First, these alterations may play a unique role in the recovery of activity following a perturbation. While the channel densities were substantially altered during the perturbation, after the perturbation they did not change much (Figure 5B, Black Arrow). In contrast, channel voltage-dependences appeared to change greatly post-perturbation (Figure 5C and 5D, Black Arrow). The channel densities appeared to have been altered in a way so that only changes in channel voltage-dependences were necessary to recover bursting following the perturbation. Second, the changes resulting from shifting activation and inactivation curves may, perhaps, be reconceptualized as alterations in window currents. To see this, observe the relative positions of the activation curves and inactivation curves for I_CaT, I_CaS, I_A, and I_Na (the other currents in the model do not possess both curves). At the time the neuron recovered bursting during the perturbation (Figure 5A, Blow-Up 3), the model shifted the activation curves to more hyperpolarized potentials and inactivation curves to more depolarized potentials (Figure 5C and 5D, Blue Arrow). In effect, the excitability of the neuron increased following the perturbation due to changes in window currents.

Discussion

Timescale of Channel Voltage-Dependence Alterations

Post-translational modifications alter the voltage-dependent gating mechanisms of ion channels, affecting their open probability and thereby influencing neuronal excitability (Hille 2001, Levitan 1988). Increases in intracellular calcium are associated with signaling pathways that trigger these modifications (Kandel et al 2000, Levitan 1988). Intracellular calcium is also thought to play a crucial role in maintaining neuronal activity set points (LeMasson et al 1993, Turrigiano et al 1994). Together, these observations suggest that calcium-mediated post-translational modifications impact a neuron’s ability to maintain a target activity profile. Indeed, experimental evidence supports the idea that neurons can modify their channel voltage-dependence properties to achieve activity targets (Gasselin et al 2015, Haedo & Golowasch 2006, Thoby-Brisson & Simmers 2002). While alterations in channel voltage-dependence are known to plastically modulate a neuron’s excitability (Levitan 1988), their role in stabilizing a neuron around an activity target has not been extensively explored. As such, we designed a study to explore how modifications to channel gating properties can influence the attainment of an activity target.

We demonstrated that altering the timescale of these modifications is important in determining the neuron’s electrical activity pattern and underlying intrinsic properties. A way to understand this work is outlined in Figure 6. The space of activity characteristic (Figure 6A) or underlying intrinsic parameters (Figure 6B) contains a region representing all electrical patterns consistent with the neuron’s activity target. When channel gating properties are altered quickly in response to deviations from the target activity, the resulting electrical patterns are shown in Figure 6 as the orange bubble labeled “τ_half = 6 s”. As the timescale of alterations slow, the activity characteristics and model parameters of the electrical patterns shift, as illustrated by the other orange subregions in Figure 6. This 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.

Figure 6

Figure 6 tacitly suggests the existence of a continuous region within the space of activity measurements or model parameters that aligns with the calcium target. If this continuity is indeed present, neurons could, in response to perturbations, seamlessly transition to new patterns that meet the activity target. Such adaptability may be important in scenarios like temperature perturbations, where it is postulated that neurons continuously adjust their intrinsic currents to manage temperature changes (Alonso & Marder 2020). In certain rhythmic neuronal networks, the relative phase of neuronal firing is preserved even as the overall network frequency increases (Hooper 1997, Tang et al 2010). To target intrinsic properties that may be amenable to maintaining phase in higher temperatures, a neuron may control how quickly it adjusts its (in)activation curve in responses to deviations from its target activity as temperature changes.

Our simulations show that altering the timescale of channels’ voltage-dependent gating properties impact the targeted 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, for instance Protein Kinase A (Fraser & Scott 1999) or cAMP dependent protein kinase (Johnson et al 1994). Over slightly longer timescales, these modifications may proceed more slowly as second messengers diffuse across the cytosol and potentially degrade (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 time scales, such as S-nitrosylation. In fact, all the channels used in this study are known change their voltage gating properties in response to phosphorylation or S-nitrosylation: I_Na (Ashpole et al 2012, Bielefeldt et al 1999), I_Kd (Bai et al 2004, Dilly et al 2004, Thompson et al 2017), I_H (Difrancesco & Tortora 1991, Ingram & Williams 1996), I_CaS (Bai et al 2004, Peterson et al 1999), I_KCa (Lang et al 2000, Weiger et al 2002), I_A (Hoffman & Johnston 1998, Varga et al 2004), and I_CasT (Welsby et al 2003).

On intermediate timescales, the channel’s voltage-dependence properties can be modulated by either the association of calcium activated subunits, for example via Potassium Channel-Interacting Proteins (Wu et al 2023), or cleavage of subunits by calcium mediated proteins, such as calpain (Friedrich et al 2004, Sandoval et al 2006). On even slower timescales, the neuron can modulate the transcription and translation of these auxiliary subunits based on intracellular calcium levels. For instance, Potassium Channel-Interacting Proteins are known to be controlled by CREB (Wu et al 2023). Additionally, the composition of the lipid rafts that house these ion channels may impact their gating properties (Dart 2010, Levitan et al 2010). Documented cross talk between the pathways that modify the lipid composition of these rafts and calcium signaling pathways may further influence a channel’s gating properties during slow changes in intracellular calcium levels (Pulli et al 2018).

Model Assumptions

This model of channel voltage-dependence alterations is based on certain assumptions. The concept outlined in Figure 6 is not expected to depend on these choices. However, examining these assumptions reveal additional insights.

The notion of ion channel gating alteration speed is simplified in two ways in this model. First, the model assumes a single, uniform time constant for modifying all intrinsic currents. In reality, each intrinsic current is likely influenced by a combination of mechanisms. For example, the relative contributions of anchored second messengers and freely diffusing kinases in phosphorylating key residues in L-Type calcium channels are not well understood and may well involve both processes (Weiss et al 2013). Consequently, the timescale of modulation for I_CaS—the current to which L-Type calcium channels putatively contributes—likely represents a weighted average of these mechanisms. Similarly, I_Kr, is modulated by PIP2 (a metabolite of calcium-mediated second messenger pathways) and a beta subunit (Li et al 2011). While the presence or absence of PIP2 can change on a fast timescale, the synthesis of the beta subunit occurs over much slower timescales, resulting in mixed temporal dynamics. These examples illustrate that the balance between fast and slow regulatory mechanisms likely varies among intrinsic currents, complicating the assumption of a single time constant. Second, our framework assumes post-translational modifications are added as quickly (or slowly) as they are removed. This simplification fails to account for processes where rapid calcium increases lead to prolonged downstream effects. For instance, the calcium-mediated second messenger CaMKII undergoes autophosphorylation, enabling a brief calcium signal to induce long-lasting changes in protein activity (Hudmon & Schulman 2002, Kennedy 1989). This autophosphorylation is known to influence intrinsic excitability neurons (Sametsky et al 2009) and may modulate the voltage-dependence of Ca_V2 channels (Jiang et al 2008), which putatively contribute to I_CaS. Accounting for these processes would require a more complex model capable of capturing asymmetries in modification rates.

Another assumption we make is that there are limits on the extent to which (in)activation curves can shift. This limitation is implemented through the cubic term in Equations 16 and 17. These constraints reflect the fact that shifts in activation curves are not unlimited. For example, while cAMP can shift the half-activation of I_h, this effect saturates at high cAMP concentrations (Difrancesco & Tortora 1991). In our study, however, a uniform constraint is applied across all intrinsic currents. In reality, different intrinsic currents likely have unique limits on how much they can shift, requiring distinct functional forms to model their unique restrictions.

Furthermore, we assumed that activity-dependent alterations in channel voltage-dependence are implemented by changing half-(in)activations. We have not, however, considered the (in)activation curve’s slope or the gating variable’s time constant. In this study, the gating variable’s voltage-dependent time constant shift by the same amount as the half-(in)activation. This preserves how quickly the gating variable respond to changes in voltage – regardless of where the activation curve is centered. However, this coupling between the time constant and the activation curve is not necessarily required. For instance, activity-dependent shifts in voltage-dependent time constants have been observed without concomitant shifts in activation curves (Thoby-Brisson & Simmers 2002). The voltage dependence of the time constant can, in principle, change in other ways as well. Also, we presume that the slopes of the (in)activation curves do not respond to changes in activity. The slope of (in)activation curves defines the steepness of the voltage dependence for ion channel gating, thereby influencing the voltage sensitivity of the transition between states (e.g., from closed to open or from inactivated to deinactivated). In our model, the slopes of these curves are held constant. Consequently, the voltage ranges remain fixed regardless of shifts in the centering of the activation or inactivation curves. However, this assumption may not be true. For example, with certain calcium-activated potassium currents, the center and slope of the activation curve depends on calcium concentration (Wang & Brenner 2006). To the extent that calcium entry during activity initiates downstream signaling processes that alter ion channel function, this observation supports the idea that activity could influence the slope of the I_KCa activation curve.

In the model, the same activity targets are used to drive changes in both channel density and voltage dependence. This assumes extensive crosstalk between the calcium pathways posited to be responsible for ion channel insertion, deletion, and post-translational modifications. However, the extent to which this crosstalk occurs is not known. However, it is known that cells can detect the temporal patterns and specific pathways of calcium entry, leading to specific activity-dependent changes (Bito et al 1997, Fields et al 1997, Gallin & Greenberg 1995). Therefore, it is likely that activity-induced calcium influxes trigger changes along distinct pathways – with some impacting only channel gating properties, some impacting channel density, and some impacting both.

It is not necessarily true that all activity-dependent changes in the properties of intrinsic currents are directly controlled by calcium influx. For example, neurons coregulate I_A and I_H expression, so that a change in one current automatically induces a change in the other (MacLean et al 2005). Neuromodulatory environment may also play a role in coordinating ion channel levels (Khorkova & Golowasch 2007). In such cases, one may use this model for the channels that are controlled by activity and augment it by explicitly modelling the remaining mechanisms that are independent of activity.

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 no electrical activity to possessing it (Baccaglini & Spitzer 1977). As shown in Figure 6, the timescale of (in)activation curve alterations define 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 & Spitzer 1977, Gao & Ziskind-Conhaim 1998, Goldberg et al 2011, Hunsberger & Mynlieff 2020, McCormick & Prince 1987, Moody & Bosma 2005, O’Leary et al 2014, Picken Bahrey & Moody 2003).

Additionally, our results suggest that that activity-dependent regulation of channel voltage dependence can help to homeostatically stabilize a neuron’s activity in response to perturbations (Figure 5). Although homeostatic activity-dependent alterations in channel voltage-dependence have been reported in a few studies of activity-dependent homeostasis (Gasselin et al 2015, Haedo & Golowasch 2006, Thoby-Brisson & Simmers 2002), many studies report these changes are small, variable, and ultimately statistically insignificant (Desai et al 1999, Mizrahi et al 2001, Turrigiano et al 1995). However, even small, coordinated shifts in the half-(in)activations of multiple ion channels may collectively have a significant impact. In our simulations, the half-(in)activation changes are relatively minor for some currents. Despite this, these currents act in concert to shape the observed electrical activity profiles, and their combined effects are evident in the overall activity characteristics (Figure 4A). Thus, averaging small shifts in (in)activation curves across a population may obscure the important role of these small, coordinated changes in multiple (in)activation curves have on an individual neuron’s activity. As such, the small changes in these channel voltage dependences observed in studies of activity-dependent homeostasis may be significant. A similar issue regarding averaging maximal conductances has been reported (Golowasch et al 2002). As such, our findings may have implications for conditions such as addiction (Kourrich et al 2015) and Alzheimer’s disease (Styr & Slutsky 2018), where disruptions in homeostatic processes are thought to contribute to pathogenesis.

In conclusion, our findings suggest that when the 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 are insufficient to attain bursting (Figure 2). 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

Box 1 presents the equations used for the model: a Hodgkin-Huxley type equation with 7 intrinsic currents, Ca²⁺ sensors, and an activity-dependent regulation mechanism.

The Hodgkin-Huxley type neuron consisted of 7 intrinsic currents (Box 1, Supplement 1, Equations 1-3) and an intracellular Ca²⁺ concentration calculation (Box 1, Supplement 1, Equation 4). The neuron’s capacitance (C) was set to 1 nF. The currents were labeled by the index i and were in order: fast sodium (I_Na), transient calcium (I_CaT), slow calcium (I_CaS), hyperpolarization-activated inward cation current (I_H), potassium rectifier (I_Kd), calcium-dependent potassium (I_KCa), and fast transient potassium (I_A). Additionally, there was a leak current (I_L). The maximal conductances were given by g̅_i. Activation and inactivation dynamics were given by m_i and h_i, respectively. The exponents q_i 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 mV, −20 mV, −80 mV and −50 mV respectively. The calcium equilibrium potential was computed via the Nernst Equation (Box 1, Supplement 1, Equation 5) using the computed internal Ca²⁺ concentration ([Ca²⁺]_i) and an external Ca²⁺ concentration of 3 mM (temperature was 10 °C).

The activation curves and inactivation curves were given by m_i,∞ and h_i,∞. Their associated time constants were and , respectively. The (in)activation curves and their associated time constants were functions of the membrane voltage, V (see tables in Box 1, Supplement 1). Note that I_KCa had an activation curve that depends on internal Ca²⁺ concentration. The internal calcium concentration was computed using I_CaT and I_CaS (Box 1, Supplement 1, Equation 4).

There were three Ca²⁺ sensors that responded to different measures associated with fluctuations of [Ca²⁺]_i during a burst. They were used to determine how close these measures are to their targets. These measures were the amount of [Ca²⁺]_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 6-8 (Box 1, Supplement 2). S̅, F̅, D̅ were prespecified targets of these sensors that encoded the target activity pattern. Box 1, Supplement 2, Equations 9-11 created moving averages of the mismatches between the calcium sensors and their targets. These moving averages were taken over two seconds (τ_S = 2000 ms).

These mismatches between the moving averages and their targets were combined in Equation 12 (Box 1, Supplement 2) to create a “match score”, S_F (Equation 12). The L⁸ norm, defined as , is used in Equation 12 to combine the errors. This expression was slightly modified from Alonso et al (2023). This formulation is particularly sensitive to the largest error values. By raising each error to the 8th power, large deviations became even more prominent in the combined score, allowing us to examine the errors in each sensor independently.

Equations 13 and 14 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 14 (Box 1, Supplement 2), 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 13 (Box 1, Supplement 2) then effectively computed a 2 second time average of Equation 14. Equation 13 was an activity sensor, ⍺, that was plotted in Figure 1. It took the value zero when all three Ca²⁺ sensors were close to the Ca²⁺ 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 Ca²⁺ sensors reached their targets. It also served as a measure for how well the model satisfied the Ca²⁺ 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, S_f, was greater than ∼0.32, the model considered the calcium targets satisfied.

Using the information regarding how well the Ca²⁺ 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, , and inactivation curve shifts, , from their specified locations ( for activation curves and for inactivation curves), were altered by the model in response to deviations from the calcium targets (Box 1, Supplement 3, Equations 15-17). The latter were adjusted on a time scale given by τ_half and the former by τ_g. The index i ran from 1 through 7 for the maximal conductances (Equation 15) and half-(in)activations (Equation 16). In Equation 17, 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 16 and 17 were new equations added to the model. They were built in analogy to Equation 15. In Equation 15, the maximal conductances were altered when the calcium sensors didn’t match their targets. There were two parts to Equation 15: a part that converted calcium sensor mismatches into changes in maximal conductances and a cubic term. The cubic term prevented unbounded growth in the maximal conductances. We set the constant associated to the cubic term, γ, to 10⁻⁷. The other term adjusted the maximal conductance based on fluctuations of calcium sensors from their target levels. The parameters A_i, B_i and C_i informed how the sensor fluctuations contributed 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. This prevented the maximal conductance from becoming negative and scaled the speed of evolution so that smaller values evolved slower and larger values faster. As such, multiplication by g̅_i effectively induced exponential changes in the maximal conductance, allowing g̅_i to vary over orders of magnitude.

Equations 16 and 17, introduced in this paper, enabled the adjustment of the ion-channel voltage dependence. Equations 16 and 17 also had two components: a component that converted calcium sensor mismatches into changes in half-(in)activations and a cubic term. The cubic term prevented unbounded translations. The cubic term’s constant, γ̂, was set to 10⁻⁵. The other term adjusted the half-activations (Equation 16) and half-inactivation (Equation 17) based on fluctuations of calcium sensors from their target levels.

Finally, we discuss how the coefficients A_i, B_i and C_i in Equation 15, and the corresponding parameters in Equation 16 and 17 were chosen. The same coefficients used to adjust maximal conductances were used for the inactivation curve shifts, and . The coefficients of the activation curve shift are their negatives, Âi, B̂i, and Ĉ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, I_Na, harder to activate, by shifting the activation curve of I_Na to a more depolarizing potential. This implies 0. Therefore, the coefficients Âi, B̂i, and Ĉi must be non-positive. The coefficients for I_Na in Liu et al (1998) were given as A₁ = 1, B₁ = 0 and C₁ = 0 and were chosen so that , thereby decreasing the maximal conductance of sodium in the same scenario. So, we flipped the signs to give and . In this way, only F Ca²⁺ sensor contributed to the modulation of I_Na, but the activation curves were now adjusted in a homeostatic way. Furthermore, the model makes it harder to de-inactivate I_Na by shifting its activation curve to more hyperpolarizing values . In the scenario where the activity sensors are greater than their activity targets, this corresponds to setting , and ,

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

Note, the coefficients in Liu et al (1998) don’t perfectly correspond with this reasoning. In Liu et al (1998), some coefficients were assigned values to ensure convergence of the model. Nevertheless, we followed the prescription above and kept this caveat in mind (Box 1, Supplement 3, Table 1).

Starting Parameters

In total, the model was given by 40 ordinary differential equations: 13 equations described the neuron’s electrical activity, 9 equations described calcium sensor components, and 18 equations that altered the maximal conductances and half-(in)activations (Box 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 nS and 0.9 nS and half-(in)activation shifts are randomly offset between −0.5 mV and +0.5 mV from their specified half-(in)activation values ( and in Box 1, Supplement 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 [Ca²⁺]_i = 0.4 μM. Equations 9-11 (Sensor Errors) and Equation 13 (Activity sensor, ⍺) were initialized with the value 1. H_X (Equation 7) and M_X (Equation 6) 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)activations.

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 milliseconds regardless of step size (step size of the integrator was adaptive and left unspecified). Maximum number of iterations was set to 10¹⁰ and relative tolerance was set to 0.0001.

Burster Measurements

We defined regular bursters as activity patterns that exhibited groups of spikes separated by regular quiescent gaps in activity. We developed methods to analyze the activity patterns and identify regular bursters. To this end, we measured the following activity characteristics: period, burst duration, interburst interval, spike height, maximal hyperpolarization, and slow wave amplitude. Peaks were identified as local maxima in the activity pattern. We used the ‘findpeaks’ function in MATLAB on the last 90 secs of activity to identify local maxima in the trace, setting the minimum ‘prominence’ at 2. The latter parameter rejected certain local maxima from counting as peaks. This parameter is one we arrived at by trial and error. Figures 7-B1 and 7-B2 show red arrows where this function identified peaks. Note in Figure 7-B2, panels A-C have local maxima that are not counted as peaks (there is one local maxima following the sequence of peaks in the burst). Our choice of ‘prominence’ in the ‘findpeaks’ function was intentionally set to exclude these. To find troughs, we took the negative of the activity profile and found local maxima using the ‘findpeaks’ function. For this, the same ‘prominence’ parameter was used over the same time window. The arrows in Figure 7C show an example of what was identified by the algorithm as a peak in this case. Troughs were defined as the negative of the identified points. The red arrow in Figure 7C shows the maximum hyperpolarization of a cycle obtained using the method described below. The green arrows are the troughs associated with each spike in the burst (Figure 7C).

Figure 7

Maximum hyperpolarization was found as follows. First all troughs were identified in the manner described above. Then, starting at −85 mV, we checked whether there were any troughs below this value. If there were none, we incremented by 1 mV to −84 mV and checked again if there were any troughs below this value. This process was repeated until troughs are detected. The troughs identified this way demarcated the ends of a cycle in the activity pattern. They were averaged to obtain the maximum hyperpolarization measurement for the activity pattern.

The time interval between successive maximum hyperpolarizations in the activity was used to measure the cycle period (Figure 7C). We averaged these cycle period measurements to give a measurement for the activity pattern’s period. To obtain the burst duration, we segmented the activity pattern into cycles using the time of maximum hyperpolarization (Figure 7-B1). The burst duration of a cycle was measured as the time between the first and last spike in the cycle. The interburst interval of the cycle was computed by subtracting the cycle’s burst duration from the cycle’s period measurement.

To measure the slow wave amplitude, we used the segmentation described above to identify the troughs of each cycle. The slow wave amplitude of a cycle is the difference between the maximum hyperpolarization (Figure 7C, red arrow) during a cycle and the average of the membrane potential at the troughs associated with each spike in a burst (Figure 7C, green arrows). The slow wave amplitudes across all cycles were averaged to obtain a measurement for the activity pattern’s slow wave amplitude.

The spike height was obtained by finding the peaks of the membrane potential during each cycle. A measure for the cycle’s spike height was obtained by averaging over all peaks, subtracting the cycle’s maximum hyperpolarization, then subtracting the cycle’s slow wave amplitude. This value was averaged over all cycles to obtain a measure of the activity pattern’s spike height.

We distinguished regular bursting from tonic spiking or irregular bursting using the following procedure. We classified an activity pattern as a regular burster when the following three criteria were met: 1) there was more than 1 spike per cycle, 2) the maximum hyperpolarization measurements for all cycles were within 3 mV of one another, and 3) the spikes per cycle varied by no more than 1 spike. The first condition eliminated tonic spiking, the second condition is a crude measure to check the average voltage wasn’t increasing, and the third was a condition that checked for irregular bursting. In addition, we verified that the cycle-to-cycle measurements exhibited low variability to ensure the appropriateness of averaging across all cycles. The relative standard deviation for all cycle-to-cycle measurements was less than 3% for the assembled bursters.

Activity Target Selection

We needed to identify a set of sensor targets that consistently produced regular bursters, regardless of the starting parameters selected. To determine whether a set of calcium targets consistently produced a burster, we monitored the activity sensor ⍺. ⍺ approached 0 when the calcium sensors were close to their specified targets. There were small oscillations in ⍺ that resulted from calcium fluctuations associated with the neuron’s activity. These fluctuations were smoothed out by taking a running average over the 8 previous seconds. When a simulation fell below ⍺ = 0.02 for the last time and remained there for 100 minutes, we considered the simulation to have “assembled” a burster that satisfied the calcium targets.

To ensure the activity pattern obtained did not depend on feedback from the activity regulation mechanism, the mechanism was then turned off (τ_half = ∞, τ_g= ∞) and the simulation run for an additional 9000 secs to ensure ⍺ < .02 for the duration. This step was taken in the initial generation of bursters (111 in total, described below). This step was not taken for the analysis where the time constants of the half-activation alterations were changed from τ_half = 6 s (Figures 3 and 4).

‘Consistently’, in this context, meant that the model assembled bursters from over 95% of the randomly chosen starting parameters. The remaining simulations did not complete assembly within the given simulation period. For the calcium targets we found (S̅ = 0.03, F̅ = 0.25, and D̅ = 0.02), our model ‘consistently’ produced bursters. To obtain these targets, we tried different combinations of calcium targets until a suitable set of targets consistently produced regular bursters.

Groups of Bursters

Using the obtained targets, we assembled 111 regular bursters by running the model from many ‘random’ starting parameters. From this collection, we created a Group of 5 that were very similar to each other but representative of the population. We expanded to a Group of 20 bursters by relaxing the restrictions used for the Group of 5.

To obtain this Group of 5, from this original 111 bursters, we went through three steps. First, we selected all activity patterns that had slow wave amplitude less than 25 mV. This choice had the effect of only including activity patterns that looked like those given in Figure 7-B1 and 7-B2. This procedure reduced the number of bursters to 80. Second, we restricted this group by period. The distribution of these bursters’ periods are shown in Figure 7-A1. The peak of this distribution is between 380 and 500 ms (Figure 7A1, highlighted in Red). We further restricted the range of bursters, resulting in 33 bursters (the “selection group”). This range contained about 30% of the distribution.

Through a third round of restriction, we obtained a Group of 5 from this “selection group”. This was done in two steps. In the first step, we attempted to find a group of 5 bursters that were even more closely matched in period. In the second step, we selected activity patterns that looked very similar based on two other criteria noted below. To complete the first step, we ordered the bursters from shortest to longest period. For each burster from the 3rd to the 31st position, we examined the two bursters immediately preceding and the two bursters immediately succeeding in this sorted list and computed the range of burst periods across these five bursters (Figure 7-A2). This method provided a rough estimate of the variability in periods among neighboring bursters. The burster exhibiting the smallest inter-instance period variability was selected as the first member of the group of 5 (Figure 7-A2, Red Arrow). This first burster, referred to as SP1 in Figure 1 (also shown in Figure 7-B1). This burster had a period of about 416 ms. In the second step, we symmetrically expanded the inter-instance period range around SP1’s period until four additional models with similar activity patterns were found. Similarity of activity pattern to SP1 was judged on the following criteria: spikes per burst and the slow wave only contained spikes and no other features (like minor bumps). The bursters that were deemed similar to SP1 by this procedure are shown in Figure 7B-1. These bursters constituted the Group of 5. Those that were rejected are shown in Figure 7B-2. The Group of 5 had periods within 2% of each other. We then extended this range to 20% and randomly selected an additional 15 models from the selection group to create the Group of 20.

Bootstrap Analysis

We used a bootstrap method to statistically compare the variability of burster characteristics (e.g., Period, Burst Duration, Interburst Interval, Spike Height, Maximum Hyperpolarization, and Slow Wave Amplitude) and intrinsic properties between two groups of bursters: Group of 5 (SP1-SP5) and Group of 20 (SP1-SP20). This variability is reported as relative standard deviation (Figure 2). To compare variability between groups, relative standard deviation was converted to coefficient of variation (CV). For each activity measurement, we generated 10,000 bootstrap samples to estimate the coefficient of variation (CV) for the Group of 5 and the Group of 20, separately, then calculated the difference in CVs between the resampled groups. We then computed the 95% confidence interval for the differences in CVs by identifying the 2.5th and 97.5th percentiles of the bootstrap distribution. The observed difference was computed by finding the actual difference between CVs between groups. If the observed difference fell inside the 95% confidence interval, it was considered statistically insignificant.

Data availability

The code and parameters required to reproduce the results in this work are accessible and at Zenodo. Repository: [https://github.com/YugarshiM/fastTimescaleHomeostasis] [DOI: 10.5281/zenodo.12585617]

Acknowledgements

The authors acknowledge the support of NIMH (R01MH046742) and The Swartz Center for Theoretical Neuroscience at Brandeis University. YM acknowledges valuable discussions with lab members. Also, YM acknowledges assistance with code preparation from Gwendolyn Harris.

Box 1

Box 1, Supplement 1

I_CaT and I_CaS are the i = 2 and i = 3 parts of the sum in Equation 1.

Table 1

Table 2

Table 3

Table 4

Box 1, Supplement 2

G_F, G_S and G_D are set to 53, 3, and 1, respectively. Also, Δ_F = .1, Δ_S = .008, and Δ_D = .015.

Table 1

Table 2

Box 1, Supplement 3

Table 1:

In this paper, the values Â_i, B̂_i, and Ĉ_i are the negative of A_i, B_i, and C_i, respectively. The values , and are equal to A_i, B_i, and C_i, respectively. The manner at which we arrive at these values are outlined in Methods. As outlined in the discussion, these coefficients do not need to take these values.