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).