Temperature compensation in a small rhythmic circuit

Temperature affects the conductances and kinetics of the ionic channels that underlie neuronal activity. Each membrane conductance has a different characteristic temperature sensitivity, which raises the question of how neurons and neuronal circuits can operate robustly over wide temperature ranges. To address this, we employed computational models of the pyloric network of crabs and lobsters. We produced multiple different models that exhibit a triphasic pyloric rhythm over a range of temperatures and explored the dynamics of their currents and how they change with temperature. Temperature can produce smooth changes in the relative contributions of the currents to neural activity so that neurons and networks undergo graceful transitions in the mechanisms that give rise to their activity patterns. Moreover, responses of the models to deletions of a current can be different at high and low temperatures, indicating that even a well-defined genetic or pharmacological manipulation may produce qualitatively distinct effects depending on the temperature.


I. Introduction
Dynamic biological systems depend on many interacting nonlinear processes that together produce complex outputs. In the nervous system, neuronal activity requires the coordinated activation and inactivation of many inward and outward currents, mediated by a variety of di↵erent ion channels. Temperature influences all biological processes, to a greater or lesser degree. This poses an inherent di culty for neuronal signaling: if the currents involved in neuronal and network dynamics are di↵erentially temperature-dependent a system that is well-tuned to work at one temperature may not function at a di↵erent temperature (Caplan et al., 2014;O'Leary and Marder, 2016;Tang et al., 2010). Nonetheless, many ectothermic animals have neurons and circuits that do function well over an extended temperature range (Robertson and Money, 2012). It then becomes important to understand how and to what extent this can occur.
The e↵ects of temperature on channel function and neuronal activity have been studied extensively for many years (Taylor and Kerkut, 1958). Increasing temperature a) Electronic mail: lalonso@brandeis.edu generally results in increases of channel maximal conductances and faster activation/inactivation rates (Frankenhaeuser and Moore, 1963). But, importantly ionic channels of di↵erent types are a↵ected by temperature to different extents: each of these processes has a di↵erent Q 10 (Schauf, 1973;Kukita, 1982;Ru↵, 1999;Tang et al., 2010). The e↵ect of temperature on neuronal intrinsic excitability such as voltage and current thresholds can be diverse (Sjodin and Mullins, 1958;Guttman et al., 1962Guttman et al., , 1966Fitzhugh, 1966). For example, temperature experiments (18 35°C) in identified neurons in locust showed reversible changes in spike amplitude and duration of spikes (Heitler et al., 1977), and studies in the jump neural circuit of grasshopers showed that neural excitability is di↵erently a↵ected by temperature across neuronal types (Abrams and Pearson, 1982).
There are also examples of neuronal and circuit processes that are relatively temperature-compensated. For example, the f /I curves in both in molluscan and locust neurons are not substantially a↵ected by small temperature changes (6 8°C) (Connor, 1975;Heitler et al., 1977). A recent study showed that the f /I curves of auditory sensory neurons in grasshoppers remains largely una↵ected by temperature changes over a range (21 29°C) (Roemschied et al., 2014). Temperature compensation also takes place the behavioral level. For example, a recent study in hunting archerfish showed that the duration of the two major phases of the C-start-a fast escape reflex which follows prey release-were temperature compensated (Krupczynski and Schuster, 2013). The question then arises of how these behaviors are preserved when the currents in the cells and the cells' intrinsic excitability properties are di↵erentially modified by temperature.
well-characterized. Temperature compensation in the pyloric network has been explored experimentally by Tang et al. (2010Tang et al. ( , 2012 and Soofi et al. (2014). These studies show that as temperature increases the frequency of the pyloric rhythm increases but the phases of the cycle at which each cell is active remain approximately constant. Additional work on the pacemaker kernel of the pyloric network (three cells connected by gap junctions that burst synchronously) showed that temperature increases the frequency of these bursts (Q 10 ⇡ 2), but their duty cycle (the burst duration in units of the period) stays approximately constant (Rinberg et al., 2013).
Temperature compensation was explored in computational models of the pacemaking kernel (Soto-Trevino et al., 2005) by Caplan et al. (2014) and O'Leary and Marder (2016). They showed that it was possible to find multiple sets of Q 10 values for di↵erent membrane conductances processes, so that the duty cycle of the cells remained constant as their bursting frequency increased. In this work we build on the results in Caplan et al. (2014) and implemented temperature sensitivity in a model of the pyloric network (Prinz et al., 2004). We show that in these models there are multiple sets of maximal conductances and temperature sensitivities that give rise to temperature compensated rhythms, and that they can reproduce much of the phenomenology previously reported (Tang et al., 2010(Tang et al., , 2012Soofi et al., 2014;Haddad and Marder, 2018). In addition, we explored how the dynamics of the currents are modified to sustain the correct activity at each temperature. We performed this study for 36 di↵erent models and found that in all cases, the contributions of the currents to the activity can be significantly di↵erent across temperatures. The currents are not simply scaled up but instead become reorganized in a complex way: for instance, a current that is important for burst termination at 10°C may no longer play that role at 25°C. The contribution of a given current to a neuronal process can be partially replaced by another. As a consequence of this, the stability properties of the activity can be di↵erent at di↵erent temperatures. Therefore, a second perturbation for which the system is not compensated, may trigger qualitatively di↵erent responses in individual circuits depending on the temperature. These results provide a plausible hypothesis for why interactions between temperature and a second perturbation are often observed experimentally (Haddad and Marder, 2018;Ratli↵ et al., 2018).

A. Duty cycle and phase compensation in model pyloric networks
The pyloric rhythm is produced by the periodic sequential activation of three neuronal types: the pyloric dilator (PD) neurons which are electrically coupled to the the anterior burster (AB) neuron forming a pacemaking kernel, the lateral pyloric (LP) neuron and five to eight pyloric (PY) neurons. We modified the model of the pyloric circuit in Prinz et al. (2004) to include temperature sensitivity. Figure 1A shows a schematic representation of the model pyloric network studied here (diagram adapted from Prinz et al. (2004)). The model network consists of three cells, each modeled by a single compartment with eight voltage-gated currents as in previous studies Buchholtz et al., 1992;Goldman et al., 2001). The synaptic connections are given by seven chemical synapses of two types as in Prinz et al. (2004). In this model the pacemaking kernel is aggregated into a single compartment AB P D (here we refer to this compartment as P D). Following Liu et al. (1998), each neuron has a sodium current, I Na ; transient and slow calcium currents, I CaT and I CaS ; a transient potassium current, I A ; a calciumdependent potassium current, I KCa ; a delayed rectifier potassium current, I Kd ; a hyperpolarization-activated inward current, I H ; and a leak current I leak . The traces in Figure 1A show a solution of this model for one set of maximal conductances G. The traces exhibit a triphasic pyloric rhythm that consists of the sequential bursting of the P D, LP and P Y cells. This pattern remains approximately periodic so we can measure the phases of the cycle at which each cell is active, as indicated by the color labels.
Temperature sensitivity was modeled for each conductance as in Caplan et al. (2014) and O'Leary and Marder (2016) (see Methods). We increased each conductance by an Arrhenius type factor R(T ) = Q T T ref 10

10
(with Q 10 1 and T ref = 10°C) and also increased the rates of the channel kinetics ⌧ in similar fashion. We assume that Q 10 G 2 [1, 2] for the maximal conductances and Q 10⌧ 2 [1, 4] for the timescales as these ranges are consistent with experimental measurements of these quantities. Because the properties of each channel type are expected to be similar across cells we assume the Q 10 values for any given current are the same across cell and synapse types. Di↵erent values of the maximal conductances G can produce an innumerable variety of activities. The space of solutions of the model-the di↵erent patterns produced for each set of maximal conductances G-is formidably complex. This is in part due to the nonlinear nature of neuronal dynamics but also because the number of conductances that need to be specified is large (dim(G) = 8 ⇥ 3 + 7 = 31). The picture is more complicated as we consider the e↵ect of temperature because it a↵ects both the conductances and the time scales of the di↵erent processes (14 in total). Therefore in the models, changes in temperature correspond to a coordinated change-a path-in a 45-dimensional parameter space. The complexity of this problem led us to confine our study to specific questions suggested by biological observations over a permissible temperature range.
In the crab, as temperature is increased pyloric activity remains triphasic-the cells fire in the same order and at the same relative phases-while the pyloric network frequency increases by a 2 to 3 fold factor over a  Prinz et al. (2004). The three cells interact via seven chemical synapses of two types. The traces below show a representative solution that exhibits a triphasic rhythm: the activity is approximately periodic and the cells burst in a specific sequence: P D-LP -P Y . B Activity of a temperature compensated model network at di↵erent temperatures 10°C 25°C. As temperature increases the frequency of the rhythm increases but the duty cycle of the cells-the burst duration in units of the period-remains approximately constant. C The panels show average values over 30 seconds for 16 values of temperature between 10°C 25°C. Top: average burst frequency of each cell. Middle: average duty cycle of each cell (cell type indicated in colors). Bottom: average phases of the cycle at which bursts begin and terminate. As temperature increases these phases remain approximately constant.
15°C range (Tang et al., 2010). For this to happen in the model, there should be many regions in parameter space for which the activity is triphasic at di↵erent frequencies and paths in parameter space-sequences of values of G and ⌧ -that interpolate these regions. One of the main results in this work is that such paths do exist and that temperature compensation is possible in this model. Finding these temperature compensated networks is nontrivial: it requires the specification of 8 ⇥ 3 + 7 = 31 maximal conductances G and 24 Q 10 values. Because it is virtually impossible to find these sets of parameters at random we employed a landscape optimization scheme described previously (Alonso and Marder, 2019) and adapted it to this particular scenario (see Methods). An alternative approach was described in O' Leary and Marder (2016). We first searched for sets of maximal conductances G that displayed triphasic rhythms at control temperature (10°C). Inspired by physiological recordings of the pyloric rhythm of crabs and lobsters we targeted control activities with a network frequency of 1Hz, duty cycle ⇡ 20% for P D, and duty cycle ⇡ 25% for LP and P Y (Bucher et al., 2006;Hamood et al., 2015). We then selected 36 of these models and for each of them we searched for Q 10 values so that the activity is preserved over a temperature range. We found that, regardless of the maximal conductances G, it was always possible to find multiple sets of Q 10 values that produced a temperature compensated pyloric rhythm. Figure 1B shows the activity of one example network at several temperatures. As temperature increases the frequency of the rhythm increases but the phases of the cycle at which each cell is active remains approximately constant. This is consistent with the experimental results reported by Tang et al. (2010). Figure 1C shows the same analysis as in Tang et al. (2010) performed on this particular model. We simulate the model at several temperatures over the working range and compute the average duty cycle, frequency and phases of the bursts. While the burst frequency increases with Q 10 ⇡ 2, the duty cycle of the cells, and the phases of the onsets/o↵sets of activity, remain approximately constant. Figure  1SA supplemental shows the same analysis performed on three di↵erent models. In all cases the duty cycle and the phases stay approximately constant while the frequency increases. Figure 1SB shows the values of the maximal conductances G and Q 10 for each of these models.

B. E↵ect of temperature on membrane potential
Experimental studies on other systems reported e↵ects of temperature on resting membrane potentials (Klee et al., 1974), spiking thresholds, and amplitudes of spikes (Heitler et al., 1977). In the models studied here, some temporal properties of the activity remain approximately constant over a temperature range but the precise shape of the waveforms can change. To inspect how the membrane potential changes over the working temperature range (10°C 25°C) we computed the distribution of V of each cell for R = 101 values of temperature. Temperature changes features of the voltage waveforms that result in changes in these distributions. Figure 2A shows example traces and the result of this analysis for the LP cell in one model. The membrane potential changes with temperature in several ways: the spike amplitude decreases and the spike threshold depolarizes. At each temperature, we computed the distribution of V and used a gray scale to indicate the fraction of time that the cell spends at each voltage value (y-axis). The color lines match features of the voltage waveforms at each temperature, with features in the distributions on the right. The distributions permit visualizing changes in the waveform as the control parameter is changed (Alonso and Marder, 2019). Figure 2B shows the membrane potential distributions of each cell for 6 models. In all cases the waveforms show visible di↵erences across temperatures. In some cases the peak voltages of the spikes-indicated by the upper envelope of the distributions-remains relatively constant while in other cases it changes. There is considerable variability in the amplitudes of the oscillations across models. The distributions also show how some features of the sub-threshold activity change with temperature. In particular, in the models on the top, we observe an increase in the spiking thresholds while in the models in the bottom the spiking threshold remains relatively constant or decreases. In addition, note that cells can either become more depolarized or more hyperpolarized as temperature increases.

C. Spiking patterns during temperature ramps
In the models, as temperature is increased the values of the maximal conductances increase and the time scales of the channels kinetics decrease (the processes become faster). While we showed previously that the duty cycle remains-on average-approximately constant, the spiking patterns of the cells do show discernible di↵erences across temperatures. To show this we subjected the models to linear temperature ramps from 10 C to 25 C over 60 minutes, as depicted in Figure 3A, and recorded the spike times. Figure 3B shows the inter-spike interval (ISI) distributions for each cell over the working temperature range (10 25 C). Notice that at each temperature, there is one large ISI that corresponds to the inter-burst interval and many small ISI values, which correspond to the spike intervals within a burst. As temperature is increased the spiking pattern changes in complicated ways and the largest ISI decreases monotonically, consistent with the overall increase in bursting frequency of the cells.
We define a spike to be a part of a burst if it occurs within 200msec of another spike. Using this criterion we can compute the duty cycle and instantaneous frequency of each burst. Figure 3C shows the instantaneous burst frequency of every burst in the P D cell. As temperature increases the frequency of the bursts increases and there is little variability from burst to burst. Figure 3D shows the distribution of duty cycles for each cell over the working temperature range. Notice that di↵erent values of temperature result in slightly di↵erent spiking patterns which in turn result in di↵erent values of the duty cycle. For some temperatures, the duty cycle is nearly identical in every burst resulting in a single point in the y-axis. For example, this is the case for the temperatures indicated by the pink box in the PY cell. There are temperatures for which the duty cycle takes two values (blue box in PY cell) and there are temperatures for which the duty cycle di↵ers noticeably from burst to burst, and produce a filled region in the diagrams. In most if not all cases, this di↵erence corresponds to a last spike that may be absent in some bursts. While this is unsurprising it is interesting that it happens in some temperature ranges and not in others.
Although in all models the average duty cycle stays approximately constant as temperature is changed, the duty cycle distributions can show marked di↵erences across models: di↵erent maximal conductances and different Q 10 values produce di↵erent patterns. We focus on these spiking patterns because they can be compared directly to both intracellular and extracellular experimental recordings. Figure 4 shows the duty cycle distributions of the LP cell over the working range for all 36 models. We chose to show this cell because the biological network contains 2 P D cells and 5 P Y cells but only one LP cell which makes comparisson to experiments simpler. In all cases there are temperatures for which the duty cycle takes almost exactly one value, ranges where it takes two values, and ranges where multiple values are produced. Overall Figure 4 shows that these patterns can be diverse across models and provides a notion of how variable the bursting patterns of the cells can be across temperatures and individuals.

D. Hysteresis and multistability
Hysteresis and multistability are often observed in conductance based models of neuronal activity (Malashchenko et al., 2011;Marin et al., 2013), and the network studied here also exhibits such features. We explored hysteretic behavior in the models by increasing temperature linearly from 10°C to 35°C in 30 minutes and then decreasing it back to 10°C at the same pace. Figure 5 shows the spiking patterns for two models during one such temperature ramp (schematized on top). The ISI distributions show normal bursting activity up to a critical temperature where the model becomes quiescent. The quiescent state remains stable as tempera- FIG. 3. Spiking patterns during temperature ramps. The spiking patterns of the cells exhibit interesting dependences with temperature. A Schematic representation of the temperature ramp used in these simulations. Temperature was increased linearly between 10°C and 25°C over 60 minutes. B ISI of each cell over temperature (y-axis is logarithmic). The spiking patterns consist of one large ISI corresponding to the inter-burst interval ( 200msec) and several smaller values produced by spikes within bursts. There are ranges of temperature over which the ISI takes almost the same values and the patterns result in nearly identical sets of points, and there are ranges where there is more variability. C Instantaneous frequency of each burst of the P D cell. D Duty cycle of each burst in each cell. While the duty cycle remains approximately constant on average, there are ranges of temperature for which the duty cycle alternates between two or more values (pink and blue shaded boxes). ture continues to increase and eventually recovers activity during the decreasing ramp, at a lower temperature than when it stopped spiking during the increasing ramp. The temperature at which the model switches from spiking to quiescence is di↵erent in each side of the ramp. This implies that there is a range of temperatures where the model can be either quiescent or spiking, and thus is multistable.
Multistability is ubiquitous in these models. The bottom panel in Fig. 6 shows a di↵erent model in which these features are more salient. At high temperatures (for times between 25 and 35 minutes) this model produces ISI values as long as ⇡ 1sec and it remains producing spikes at all temperatures. The spiking patterns di↵er noticeably on the up and down ramps. These observations suggest that at many temperatures, there are multiple stable attractors that correspond to the triphasic rhythm but di↵er slightly in their spiking patterns. In many cases these di↵erences come down to the precise timing of the last spike of a burst.

E. Dynamics of the circuits at high temperatures
For temperatures higher than a critical temperature, the biological network displays anomalous regimes characterized by the emergence of slow time scales and intermittency between what appear to be metastable states (Tang et al., 2012;Haddad and Marder, 2018). We explored how such regimes may arise in the models by direct inspection of the traces at di↵erent temperatures. We found that some models do not display irregular behavior but just become quiescent, or produce subthreshold oscillations, after a critical temperature. Many models do display irregular and seemingly aperiodic regimes. As expected, irregular regimes typically take place for values of the temperature near a transition between qualitatively di↵erent patterns of activity. Figure 6 shows the activity of 5 models at one of their critical temperatures. The left panels show 30 seconds of data and the right panels show 5 seconds. Notice that in FIG. 4. Duty cycle distributions of all models. There is significant variability across models in the way the duty cycle distributions change with temperature. The figure shows the distributions of duty cycle of the LP cell for all 36 models during the same temperature ramp as in Fig. 3. The precise shape of these patterns depends on the values of the maximal conductances and temperature sensitivities.
all cases there are timescales that are much larger than the control period of 1sec. The dynamics of the models in these regimes is daunting and its characterization is beyond the scope of this work. However, qualitative features of these states-such as the appearance of slow timescales or absence/presence of activity in one or more cells-are captured by these models. This, together with the multistable regimes shown before plus a source of noise, may provide a reasonable model to account for the irregular regimes in the biological network (Tang et al., 2012;Haddad and Marder, 2018).

F. Dynamics of the currents at di↵erent temperatures
Neuronal activity is governed by currents that result from the precise activation and inactivation of ion chan-nels with di↵erent kinetics. Experimental access to these quantities is limited since it is hard to measure currents individually without blocking all other currents, and thus changing the activity. Here, we employ the models to explore how these currents are di↵erentially altered by temperature, and yet are able to remain balanced in such a way that the network activity speeds up while preserving its behavior and phase relationships.
We first compare the dynamics of the currents in two model P D neurons with di↵erent maximal conductances G and temperature sensitivities Q 10 by direct inspection of the currents' time series. Figure 7A shows the time series of each intrinsic current in P D at 10°C and 25°C for one example model. The corresponding membrane po-FIG. 5. Hysteresis and multistability. We ramped the temperature in the models from 10°C to 35°C and then back to 10°C in a symmetric fashion. The figure shows a schematic representation of the temperature ramp on the top. The bottom part shows the ISI distributions of the P D cell over time for two di↵erent models. (top) The cell ceases to produce spikes before 600 secs, which is also before reaching 35°C, and remains quiescent until temperature is ramped back down. The temperature at which the cell resumes spiking is about ⇡ 2°C lower than the temperature at which it became quiescent. This indicates that the system is multi-stable over that temperature range. (bottom) Hysteresis is more evident in this model. The spiking patterns during the down ramp di↵er visibly from those during the up ramp over a wide temperature range. tential activity V is shown in the blue traces on top. The red dashed lines indicate the peak amplitude of the currents at control temperature (10°C) to visualize whether it grows or decreases. This is indicated by the boxes in each row (in colors: red grows, green decreases, white stays the same). In this example, the Na current increases its amplitude by more than two fold over the 15°C range. The CaT current also shows a strong dependence with temperature: despite the fact that gCaT increases with temperature, the amplitude of this current decreases to about half its control value at 10°C. Similar observations can be made for other currents (such as I A in Fig.  7A). The general case is that any current can increase or decrease depending on the dynamics of neuronal firing, the maximal conductances G and temperature sensitivities Q 10 . In addition, we found that whether any given current increases or decreases with temperature can be di↵erent across cells in the same model.
The way currents are modified by temperature is different across models. Figure 7B shows the currents in the PD cell at two temperatures for another model. In this case, the Na current decreases by a small amount despite the fact that the maximal conductance and timescales of this current increase with temperature. Notice how the trends are di↵erent across models ( Fig. 7A and 7B). The KCa current which in Fig. 7A grows by a factor six, is moderately decreased by temperature in Fig. 7B. Temperature not only a↵ects the amplitudes of the currents but it also a↵ects their waveforms in interesting ways. Because the total number of currents in the circuits is 31 it becomes cumbersome to compare them across temperatures (and models). For this reason, we computed and inspected their currentscapes (Alonso and Marder, 2019). The currentscapes use colors to show the percent contribution of each current to the total inward (or outward) current and are useful to display the dynamics of FIG. 6. Dynamics of the circuits at high temperatures. The models become dysfunctional -or crash -at temperatures higher than 25°C and the do so in very di↵erent ways. The figure shows the membrane potential of five models A-E. The traces on the left show 30 seconds of data and the traces on the right show an expanded trace of 5 seconds. One hallmark of these regimes is the emergence of time scales much longer than that of the original rhythm.
the currents in a compact fashion. Figure 8 shows the currentscapes of each cell in one model at two di↵erent temperatures. The shares of each current of the total inward and outward currents are displayed in colors over time. The total inward and outward currents are represented by the filled black curves in logarithmic scale in the top and bottom. The currentscapes are noticeably di↵erent indicating that the currents con-tribute di↵erentially across temperatures. For example, in the P D cell the Leak current contributes a visible share of the inward current both at the beginning and the end of the burst at 10°C, but at 25°C its contribution is mainly confined to the beginning of the burst. The A current is more evenly distributed during the burst at 10°C and at 25°C its contribution is visibly larger at the beginning of the burst. In all cells in this example FIG. 7. Dynamics of the currents at di↵erent temperatures. Currents of the P D cell at two temperatures for two di↵erent models A-B. Temperature increases the peak amplitudes of some currents and it decreases it for others, as indicated by the + andsymbols in red and green. The way these currents are modified can vary substantially across models. model, temperature massively increases the contribution of the KCa current towards the end of the bursts. Figure 9 shows the same analysis for a second example. In this case also, the KCa currents increase their share at 25°C but is almost absent in the P D cell at 10°C. The CaS current increases with temperature in the PD cell, so temperature has the opposite e↵ect on this current as it did in the example in Fig. 8. The Leak current contributes a large percent of the current both at the beginning and end of bursts at 10°C. At 25°C its contribution at the end of the burst becomes negligible while it still contributes significantly at the beginning of bursts. The contribution of the CaT current increases visibly in P Y . In the third example in Figure 10 the contributions of the KCa current in the PD and LP cells are smaller and stay approximately constant. The contribution of the A current decreases in these cells as in previous examples, but because there is not an increment of KCa, the contribution of the synaptic currents-represented by the gray areas in the currentscapes-increases instead of decreasing. Together, these examples illustrate how a current can play di↵erent roles at di↵erent temperatures, and how diverse these mechanisms can be across individual solutions.
The currentscapes show that in all models and at all temperatures, there are periods of time where the activity is dominated by one or two currents, and periods of time where several currents contribute to it by similar shares. During spike production the activity is first driven by I Na and later by I Kd . During these periods the total inward and outward currents di↵er by orders of magnitude, and are entirely dominated by one current type. Similarly, when the cells are most hyper-polarized the dominant currents are I H and the synaptic currents I Syn . In contrast, there are periods of time-such as the beginings and ends of bursts-where the total inward and outward currents are are composed of comparable contributions of several current types. It is during these periods-when multiple currents act together to change the membrane potential-where we observe the most dramatic changes with temperature. To show this we computed the currentscapes for the 20 milliseconds following the end of bursts in the PD cells of the models in Figures 8, 9 and 10. Figure 11A shows how the currentscapes at burst termination change with temperature in the model in Fig.  8. In this case the inward currents show more changes than the outward currents: the contribution of the CaT current decreases with increasing temperatures and the contribution of CaS increases. The most salient changes in the outward currents are that I A decreases its contribution and the contribution of I KCa increases (orange areas). Di↵erences across temperature are more visible for the model in Fig. 11B (also Fig. 9). In this model the A current also decreases its share as temperature increases but this e↵ect is more pronounced than in Fig.   11A. The decrement in the share of I A is accompanied by a massive increase of the KCa current contribution. The contributions of the Leak and synaptic (SY N S) currents decrease for increasing temperatures. The inward currents also show visible di↵erences. The contribution of the Na current increases and the contribution of the calcium currents I CaT and I CaS change visibly (their overall contribution appears to decrease, but their timing is also a↵ected). The average contribution of the H current appears similar across temperatures but its dynamics is changed. It tends to increase its contribution after the end of the burst, and this increment becomes faster at higher temperatures. The third example in Fig. 11C (same model as in Figure 10) shows similarities and differences with the previous cases. The A current decreases its share but the contribution of the KCa current does not increase massively as in the previous cases. Instead, the contribution of the Leak current shows a large in- crement (pink areas). The inward currents also change with increasing temperatures: the contributions of the calcium currents CaT and CaS decrease, and the contribution of the H current increases. In total, Figure 11 shows that the ionic mechanisms that accomplish burst termination can be di↵erent across temperatures, and the way in which they change can be dramatically di↵erent across models. But importantly, in all cases there is a smooth transition between di↵erent relative fractions of currents. While all currents change their contributions with temperature, the most dramatic changes involve the relative roles of the K + currents. Figure 11 shows how the currentscapes at a specific instant change with temperature. To explore how the full currentscapes change over temperature we computed the distributions of current shares (Alonso and Marder, 2019). We computed the currentscapes at each value of temperature over the working range and inspected the time series of the currents' shares, or equivalently, the widths of the colored regions in the currentscapes. We then proceeded in the same way as in Fig. 2 and computed the distributions of these time series at each temperature. Supplemental Figure S11A shows the time series of the P D cell in one model at two temperatures (10°C and 25°C) and their corresponding membrane potential distributions (left panels). The distributions are normalized so the y-axis indicates the probability of V (p(V )) in logarithmic scale. The quantity p(V ) is proportional to the time the system spends at each voltage bin. We colored the traces using this probability as a color map. In this way, di↵erent portions of the cycle correspond to di↵erent colors: spikes appear in blue, sub-threshold activity in yellow, and brown and pink correspond to near-threshold activity. Notice that when the activity is faster this color code still parses the activity in a similar fashion. Supplemental Figure S11A and S11B shows the distributions of the current shares of all currents for 101 values of temperature between 10°C and 25°C (x-axis), for the models in Figures 8 and 9. In each panel, the y-axis corresponds to the share of each current to the total (inward or outward) and the colors indicate the portions of the cycle during which the current contributes by each amount. The distributions feature far more details than we are able to describe so we will only describe a few. For example in Fig. 11A, the H current in the PD cell contributes with less that 75% most of the time (yellow, hyper-polarized) at 10°C but its contribution increases almost linearly as temperature increases. The Na and Kd currents contributions span the full range (0% 100%) but they contribute large amounts only during spikes (in blue), while most of the time (yellow, hyper-polarized) their contributions are small (⇡ 5%). The Leak current contributes between 0% and 25% and while the distribution changes in a complicated manner, its average contribution does not appear to change noticeably. As a last example, notice that this is not the case in the P Y where a clear decreasing trend can be observed for the Leak current. Figure 11B shows the same analysis for the model in Figure 10. These examples make the case that the way the currents contribute to activity can change quite dramatically with temperature and that the details of how these changes take place can be di↵erent across models, as well as across cells.

G. Multiple mechanisms for temperature compensation
In these models, temperature increases the timescales and maximal conductances of all channels but its effect on the currents is far more diverse. As a way to characterize these changes we performed a series of simple measures on the time series of each current, and monitored these quantities over the working temperature range (10 25°C). For each current, we measured the peak amplitude (see Fig. 7 red dashed lines), the mean value over 30 seconds, and the mean share of the total (average along the y-axis in Fig. S11 for each temperature, or areas of the colored regions in the currentscapes) also over 30 seconds. Figure 12 shows examples of how FIG. 11. Currents' contributions at the end of bursts across temperature. The figure shows the currentscapes of the PD cell for the 20 milliseconds following burst termination, and how these change over temperature. Panels A-B-C correspond to the models in Figs. 8, 9, and 10.
FIG. 12. Summarizing changes in currents contributions with slopes. We quantified the e↵ect of temperature on the currents by measuring their peak amplitudes, their average values and the average percent contributions to the total over 20 seconds. A The panels show each of these quantities for the CaT current in the LP cell, for 15 models indicated in di↵erent colors. The left sub-panels show the actual values and the right sub-panels show a linear fit to these trends. B Slopes of the linear fits in colors for each of these measures for the models in Figs.8-10 (top to bottom). The average share (right panels) shows more clearly that the way in which the currents change their contributions is diverse across models. The colors indicate whether a current increases its share (red) or decreases it (blue).C Principal component analysis of the average share trends. Color scale is the same as in B. The first three components capture ⇡ 50% of the variance and show some correlated changes that occur frequently in the models studied here (N = 36).
the CaT currents in LP change over temperature, for many models indicated in di↵erent colors. The left subpanels show the peak amplitude, average, and average share of the CaT current over temperature, and the right sub-panels show a linear fit to each of those curves using the same color scheme. Note that in most cases the linear fit is a good approximation to the curve so this permits using the slope of the fits to compare trends across cells and models.
These measures provide di↵erent information. The peak currents show the largest variations. The peak currents in one cell can either increase or decrease with temperature, and this trend can be di↵erent across cells despite their having the same Q 10 values. The peak amplitude is not necessarily the most informative measure and it can be sensitive to large infrequent peaks such as those that occur in the synaptic currents when the spikes of two cells occur close in time (not shown). The average current already reveals di↵erent trends. The same currents whose peak amplitude decreased over temperature increase in the average. The third way in which we characterized these changes was by computing the average value of the shared contribution to the total. This corresponds to the width of each colored region in the currentscapes averaged over time. It is proportional to the total area of each colored region and therefore provides a measure of how the currentscapes change with temperature. This quantity behaves di↵erently to the average current and is normalized to the interval [0, 1]. Figure 12B shows the slopes of the fits of all current types in the three example models of Figs. 8, 9 and 10, for each of the quantities we monitored. In the first model (top row), the Kd current increases its peak amplitude by about ⇡ 40 nA°C in all cells (less in PY). The Na peak current also increases in PD and PY but less than Kd (⇡ 30 nA°C ) and it increases maximally in LP (⇡ 50 nA°C ). The peak amplitude of the calcium currents increases a small amount and the rest of the currents do not change their peak amplitudes. The second column shows the same analysis performed for the average value of the currents. This measure indicates similar trends as before showing growth of the Na and Kd currents, and it highlights that the average CaS currents increase more than the rest of the currents. Finally, the third column shows the trends of the average shares. In this case we used red for currents that increase their share and blue for currents that decrease it. The average contributions of currents can change by about ⇡ 7% over the working temperature range and can increase or decrease. The second and third rows in Fig. 12B show the same analyses for the models in Figs. 9 and 10. The way the peak and average currents change with temperature appears to be similar across models, but if we consider the average contributions we see that each model achieves compensation by balancing its currents in di↵erent ways.
The mean current contributions or shares reveal more di↵erences on how the currents become reorganized as temperature changes. The analyses in the third column of Fig. 12B can be compared directly with the currentscapes in Figs. 8, 9, and 10. The first row in Figure 12B, third column, shows in color whether a current increases or decreases its average contribution, for the model in Figure 8. In this model KCa grows in all cells. The contribution of the synaptic currents stays approximately constant in PD and LP but decreases in PY. The contributions of the CaT and A currents decrease in all cells by di↵erent amounts. The contribution of the CaS current changes di↵erentially with temperature across cells: it decreases in P D, it increases in P Y , and stays approximately constant in LP. All these trends can be seen by comparing the currentscapes at 10°C and 25°C as shown in Figure 8. In all cells in the second model (middle row, third column) temperature increases the average contribution of the KCa current and decreases the contribution of the synaptic currents. This is visible in the currentscapes as the KCa current takes more area at 25°C than at 10°C and reduces the average contribution of the synaptic current (area of the gray region in the currentscapes in Fig. 9). The third row in Fig. 12B shows this analysis for the model in Fig. 10. In this model the contribution of the A currents decreases in all cells. The contribution of KCa stays approximately constant in P D and LP , but it grows in PY. The contribution of the synaptic currents increases in P D and LP , but not in PY where it decreases. Overall, these three models achieve temperature compensation by completely di↵erent mechanisms.
There are multiple ways in which the currents can reorganize to produce a faster version of the activity while preserving its morphology. To gain intuition on how diverse these arrangements can be, we performed a PCA analysis across models on the trends described in Fig.  12A and 12B. We found that in the case of the peak amplitude, the first 2 components captured most of the variance (84%), which correspond to large changes in the Na and Kd currents. The same analysis performed with the average current also yields 2 components which capture most of the variance and are also very similar to the peak amplitude case (87%). The mean shares on the other hand require 6 components to explain a similar portion of the variance (84%), suggesting that this measure better highlights the di↵erent compensation mechanisms across models. Figure 12C shows the first three Principal Components of the average shares' slopes for all 36 models (right column 13B). The analysis suggests that there are certain currents that are more likely to change together. One example of this is indicated by the first component. If the KCa current increases (or decreases) its share, the synaptic current decreases (or increases) its share by a similar amount, and the H and Leak currents behave similarly to the synaptic current. The second component suggests that large increments in the KCa current in the LP neuron are usually correlated with changes in the H and synaptic currents in P Y and anti-correlated with changes in the Na current. The third component shows that in both LP and P Y an overall increment or decrement of the Na and Kd shares is anti-correlated with the H and synaptic currents.

H. Response to perturbations at di↵erent temperatures
In all models the dynamics of the currents are di↵erent at di↵erent temperatures and their contributions become reorganized in complex ways. Therefore the role any given current plays at 10 C can be di↵erent at a different temperature (25 C). The combined activity of all the variables in the model-the dynamical attractor-is therefore expected to have di↵erent stability properties at di↵erent temperatures. This means that in the models, an extreme perturbation can have qualitatively di↵erent e↵ects at di↵erent temperatures. The responses of di↵erent models to extreme perturbations such as partially or completely removing a current can be remarkably diverse This manuscript has not been peer-reviewed L.M. Alonso and E. Marder  (Alonso and Marder, 2019). While this is expected due to the complexity of these models, the fact that these responses can be qualitatively di↵erent across temperatures is perhaps less clear. To shed light on this issue we performed simple perturbation protocols. Figure 13 shows one perturbation protocol in two models. In this case we explored the interaction between temperature and removing the A-type K + channel because Tang et al. (2010) measured the e↵ects of temperature on I A . We performed the same perturbation at two temperatures: the models were simulated for 20 seconds in control conditions and then I A was completely removed for the next 20 seconds. Figure 13A (top) shows that removing the A current at 10°C has catastrophic consequences for the rhythm and the activity becomes completely irregular. When the same perturbation is performed at 25°C (bottom) the triphasic rhythm is almost normal except for the irregular bursting of P Y . Figure 13B shows the same pertubation protocol in a di↵erent model. This model displayed normal triphasic behavior at 10°(top) and when the A current was removed its activity slowed down but the triphasic rhythm was mantained. The same perturbation at 25°C results in quiescence. In this model, the A current is necessary for the activity at 25°C but not at 10°C.
To further characterize the di↵erential response to perturbations at di↵erent temperatures, we performed a second assay where we gradually decreased a current from control (100%) to completely removed (0%) in five steps. We performed these simulations for all 36 models and all current types at 10°C and 25°C and compared their responses. In general, gradually removing a current can result in qualitatively di↵erent states depending on the temperature at which this perturbation is performed. The responses are so diverse that a coarse classification scheme is su cient to highlight this observation. We classified the network activity based on the spiking patterns and waveforms of each cell, and also by their relative activities. The activities of each cell were classified as: regular bursting, irregular bursting, tonic spiking, irregular spiking, single spike bursting, quiescent and other, while the network activity was classified as triphasic or not triphasic. The classification scheme consists of a simple decision tree based on the statistics of spiking and is su cient to tease apart the di↵erent regimes we observe at a coarse level (Methods). Figure 14A shows the result of classifying the responses of one model to gradually decreasing a conductance in 5 steps, at two di↵erent temperatures (10°C and 25°C), for four conductances. We employed 5 ⇥ 4 response grids to summarize the e↵ects of gradually removing a conductance at any given temperature. The rows in the grids correspond to di↵erent conductance removals, with the control condition (100%) on top, and the completely removed condition (0%) at the bottom. The first three columns correspond to each of the cells (PD-LP-PY) and the fourth column corresponds to the network (NET). The colors indicate the type of activity of each cell and the network at each condition as indicated by the labels in the figure. The top left panel in Figure 15A shows the responses of one model to removing the sodium current I Na at two temperatures. At 10°C, when gN a ! 75%gN a the LP cell bursts irregularly and the activity of the network is not triphasic, but if the same perturbation is performed at 25°C the activity remains normal triphasic. Removing KCa in this model has little e↵ect at 10°C as the activity remains triphasic. However, this current plays an important role at 25°C and the model becomes quiescent upon complete removal of KCa. The same observation that the responses to perturbation are di↵erent at each temperature holds regardless of the perturbation type. In the cases of CaT and CaS complete removal results in di↵erent non-triphasic activity at each temperature. Figure 14B shows the results of this analysis for all 36 models for the case of gradually removing Kd (labels are as in Fig. 14A and were dropped for clarity). The figure shows that the responses of the model are in general qualitatively di↵erent at di↵erent temperatures, and also that they are diverse across models. The fact that the models "break-down" in di↵erent ways is consistent with the observation that these compensation mechanisms are indeed diverse. In all models, the responses to extreme perturbation are temperature dependent for most perturbation types. In Figure 14C we quantify how many models respond di↵erently to complete removal of a current at each temperature. This corresponds to comparing the bottom row of the response grids at each temperature, for each conductance type. In most cases, about half of the models respond to complete removal in di↵erent ways at each temperature. The exceptions are CaT for which most models respond in di↵erent ways, and H where complete removal results in equivalent states (quiescence in this case) at both temperatures. Overall our results indicate that the responses to perturbations can be qualitatively di↵erent at di↵erent temperatures.

III. Discussion
Neurons characteristically have multiple ion channels that di↵er from each other in terms of their selectivity, voltage-dependence and time courses of activation and inactivation. For example, it is not uncommon that neurons might display 12 or more di↵erent voltage and timedependent currents, including multiple classes of potassium and calcium currents. Undoubtedly, many of these channels are targeted to specific sites on neurons, but nonetheless they all will contribute to the overall electrical activity of neurons and the networks in which they are found. One of the potential advantages of this multiplicity of channels is that they may confer resilience to perturbations of many kinds (Drion et al., 2015). In this paper, we explored some ways in which multiple ionic currents help confer temperature robustness to neurons and a small network. FIG. 14. Classification of responses at di↵erent temperatures. The models' responses to perturbations are extremely diverse and are in general di↵erent at di↵erent temperatures. We classified the responses of the models at two temperatures for 5 values of gradual decrements of each current. A Response to removal of Na, A, CaT and CaS for one model. B Response to removal of Kd at two temperatures for all models. C The pie charts show the number of models for which the responses to complete removal (0%) are classified as the same (green) or di↵erent (blue) at 10°C and 25°C.
Temperature a↵ects the properties of ion channels differentially and yet many animals live successfully over wide temperature ranges. The North Atlantic lobsters and crabs routinely experience more than 25°C temperature swings, and may see 8 10°C changes within a short time period.
The temperature sensitivities of ion channels are dictated by their molecular structures, and for this reason it may be reasonable to assume that the temperature sensitivity of a given channel is similar across cells and individuals, and throughout the lifetime of the animal. This raises the question of whether it is possible for one set of Q 10 values to function appropriately over wide ranges of conductance densities. Alternatively one may consider molecular processes which could alter channel function. For example, RNA editing can change the temperature sensitivities of specific channels to facilitate adaptation (Garrett and Rosenthal, 2012;Rosenthal, 2015). A recent study in Drosophila showed that acute temperature perturbations result in changes in RNA editing levels (Rieder et al., 2015). Therefore, it is also reasonable to explore the possibility that the Q 10 values can be di↵er-ent across individuals and change over time.
Nervous systems must have mechanisms that allow them to rapidly and reliably deal with acute changes in temperature. We provide simulations that suggest that some sets of ion channels are able to deal with changes in temperature of this sort by smoothly changing the mechanisms by which neuronal activity is produced. This can be seen most dramatically in Figure 11 which shows that the role of the three di↵erent K + currents in spike repolarization gradually changes as the temperature increases. While this is true in all three examples shown in Figure 11, the details by which the currents trade-o↵ their role in spike repolarization are di↵erent, and depend on the di↵erences in their starting conductance values. This smooth transition can occur precisely because these three K + depend on voltage and time di↵erently and also are a↵ected di↵erently by temperature. Our speculation is that this kind of smooth transition is one of the reasons that neurons have so many ion channels with overlapping functions.
Transitions in the relative importance of di↵erent currents are also seen as a function of changes in membrane potential  and as a function of spike broadening in a repetitively firing neuron (Ma and Koester, 1996). So it is likely that the smooth transition in the relative roles of multiple conductances in neuronal dynamics is a general feature of all neurons with complex dynamics and many kinds of channels. This illustrates the benefit for a neuron to have multiple channels with overlapping functions.
The same principle is likely to hold when we approach network dynamics. The phase relationships of the pyloric rhythm are maintained over a considerable temperature range (Tang et al., 2010;Soofi et al., 2014). Likewise, the pyloric phase relationships are also maintained over a considerable frequency range (Hamood et al., 2015). In this latter case, it has been argued that the phase compensation depends on the conjoined action of a number of di↵erent cellular mechanisms, including synaptic depression (Manor et al., 1997;Hooper, 1997), and the activation and inactivation of I A and I H (Nadim et al., 1999;Nadim and Manor, 2000;Tang et al., 2010;Harris-Warrick et al., 1995). So, here the principle also holds: resilience and robust function may require smoothly moving between a variety of di↵erent cellular mechanisms.
In this work we employed biophysically detailed models of the pyloric network to explore how temperature compensation can be achieved in small neural circuits. We extended a model of the pyloric network developed previously to include temperature sensitivity and produced multiple temperature compensated models that capture much of the phenomenology reported in experiments (Tang et al., 2010(Tang et al., , 2012Soofi et al., 2014). One of the most significant results of this paper is that temperature compensation can take place in multiple di↵erent ways. In these models, both the membrane potential and the spiking patterns are a↵ected by temperature changes and these changes appear to be di↵erent in every model we inspected. Across the models, di↵erent values of the maximal conductances G and temperature sensitivities Q 10 result in consistent di↵erences in the duty cycle distributions. Measurements of the values of the maximal conductances in the STG show large variability across individual animals (Goaillard et al., 2009;Schulz et al., 2006Schulz et al., , 2007 so our expectation is that the duty cycle distributions of the biological cells will also display intricate dependencies with temperature, and that these distributions will be di↵erent across individuals. Temperature is not the only perturbation that crabs and lobsters experience. As with temperature, the responses of the STG to changes in pH are diverse across individuals (Haley et al., 2018), again consistent with the large amount of animal-to-animal variability in the expression of ion channels (Schulz et al., 2006(Schulz et al., , 2007Temporal et al., 2014;Tran et al., 2019). It is therefore reasonable to assume that the responses to a global perturbation are diverse across individuals because di↵erent compositions of channel densities-which produce similar pyloric rhythms-are di↵erentially resilient to any given perturbation. By the same token, as currents change as a function of temperatures we expect that a global perturbation of the same individual may have qualitatively di↵erent e↵ects at di↵erent temperatures, as is illustrated in Figure 13. The interaction between temperature and a second perturbation are the subject of recent experimental studies (Haddad and Marder, 2018;Ratli↵ et al., 2018). These studies are consistent with the interpretation that di↵erent current configurations have di↵erent stability properties and that temperature changes these configurations.
In this work we used a genetic algorithm (Alonso and Marder, 2019) to find temperature-compensated networks. It is important to reiterate that temperaturecompensated neurons and networks are di cult to find doing random searches (Caplan et al., 2014). Therefore, despite the significant animal-to-animal di↵erences in conductance densities seen across the population, randomly sampling conductance values is unlikely to result in successful solutions for individual animals. This argues that the biological mechanisms that give rise to successful temperature compensation are likely to result from rules by which correlated values of conductances are produced, such as those found in homeostatic models (O'Leary and Marder, 2016;O'Leary et al., 2013).
Taken at face value, models suggest that conductance densities drift throughout the life of an animal as a result of homeostatic processes (LeMasson et al., 1993;Liu et al., 1998;Golowasch et al., 1999b;O'Leary et al., 2014). Conductance densities can also change in an activity dependent manner as a result of perturbations (Turrigiano et al., 1994;Golowasch et al., 1999a;Santin and Schulz, 2019;Golowasch, 2019). Understanding how temperature compensation can be compatible with homeostatic and adaptation mechanisms remains an open theoretical and experimental question.
Although the models used here are semi-realistic they are far from incorporating the complexity of actual biological neurons and the pyloric circuit. Therefore, we do not expect that specific observations derived from our cohort of models-such as the co-regulation of some current pairs being more likely than other pairs-have predictive value in detail. Nonetheless, we do expect that the general principles highlighted here will hold in biological systems. In particular, the observation that the dynamics of the currents are reorganized in non-trivial ways in response to changes in temperature is likely to be a common finding in biological systems.

A. The model
The activity of the cells was modeled using singlecompartment models similar to those described previously (Turrigiano et al., 1995;Liu et al., 1998;Goldman et al., 2001;Alonso and Marder, 2019). Each neuron has a sodium current, I Na ; transient and slow calcium currents, I CaT and I CaS ; a transient potassium current, I A ; a calcium-dependent potassium current, I KCa ; a delayed rectifier potassium current, I Kd ; a hyperpolarizationactivated inward current, I H ; and a leak current I leak . The number of state variables per cell is 13. The units for voltage are mV , the conductances are expressed in µS and currents in nA. The pyloric network consists of 3 compartments with the same ion channel types but di↵erent conductance densities. The interactions in the network consist of 7 chemical synapses and are similar to Prinz et al. (2004). The synaptic current is given by I s = g s s(V post -E s ), where g s is the synapse strength, V post is the membrane potential of the postsynaptic neuron and E s is the reversal potential of the synapse. The activation of a synapse s(t) is given by with, and These equations are identical to Prinz et al. (2004) except for the inclusion of a bound for the timescale of activation ⌧ r = 20msec (we want to avoid the case that as s 1 ! 1 then ⌧ s ! 0 andṡ ! 1). All other parameters (except g s ) are identical to Prinz et al. (2004). Following Prinz et al. (2004) we set E s = -70mV and k = 1 40 ms for glutamatergic synapses, and E s = -80mV and k = 1 100 msec for cholinergic synapses. We set V th = -35mV and = 5mV for both synapse types.
Temperature e↵ects were included in this model as done previously by Caplan et al. (2014) and O'Leary and Marder (2016). Temperature dependence was introduced in the time constants of the channel-gating variables ⌧ mi and ⌧ hi , the maximal conductances g i , the time constants of calcium bu↵ering ⌧ Ca , and the maximal conductances and time constants of the synapses. This was done by replacing all conductances g i by R i (T )g i and all time scales Here T is the temperature, T ref = 10 is the reference temperature and Q 10i is defined as the fold change per 10 C from the reference temperature (i indicates the process type). Finally, the calcium reversal potential E Ca depends on temperature through the Nernst equation.
The total dimension of the model is 46 = 3 ⇥ 13 + 7 with 13 state variables per cell, and 7 variables for the state of the synapses. In this work some parameters were considered fixed while others were allowed to take values in a range. The number of parameters we varied per cell is 27 = 9 ⇥ 3: the 8 maximal conductances plus the calcium time constant ⌧ Ca . The temperature sensitivities were assumed to be the same for all cells and amount to 24 additional parameters. We need to specify 8 Q 10 values for the intrinsic maximal conductances (one per channel type), plus 12 Q 10 values for the time scales of each intrinsic process (not all currents inactivate). The number of Q 10 parameters for the synapses is 4: 2 for the maximal conductances (of each synapse type), and 2 for the timescales of activation. The models were simulated using a Runge-Kutta order 4 (RK4) method with a time step of dt = 0.05msec (Press et al., 1988). We used the same set of initial conditions for all simulations in this work V = 51mV , m i , h i = 0 and [Ca +2 ] = 5µM .

B. Finding parameters
Each model is specified by the maximal conductances G and calcium time constants of each cell (9 ⇥ 3 parameters), and the temperature sensitivities Q 10 of each process (24 parameters). We recently introduced a function that upon minimization, results in values of the maximal conductances for which the activity of a single compartment corresponds to periodic bursting with a target frequency (f tg ) and duty cycle (dc tg ). The function uses thresholds to obtain temporal information of the waveform, such as spike times, and then uses them to assign a score or error that measures how close the solutions are to a target activity. This function is described in detail in Alonso and Marder (2019) and was used here to find temperature compensated networks.
For any given set of conductances we simulated the model for 20 seconds and dropped the first 10 seconds to minimize the e↵ects of transient activity. We then computed the average (<>) burst frequency < f b >, the average duty cycle < dc >, the number of crossings with a slow wave threshold # sw = 50mV , the num-ber of bursts # b , and the average lags between bursts < P D LP > and < LP P Y >. To discard unstable solutions we checked if the standard deviation of the burst frequency and duty cycle distributions was small; a solution was discarded if std({f b }) < f b > ⇥0.1 or std({dc}) < dc > ⇥0.2. If a solution is not discarded we can use these quantities to measure how close it is to a target behavior, Here E f measures the mismatch of the bursting frequency of each cell with a target frequency f tg and E dc accounts for the duty cycle. E sw measures the di↵erence between the number of bursts and the number of crossings with the slow wave threshold t sw = 50mV (we ask that # sw = # b ). E ph compares the lags between bursts (in units of the bursting period ⌧ b = 1 <f b > ) to a target lag tg . These measures are discussed in more detail in Alonso and Marder (2019).
Temperature compensation is achieved in our models by searching sets of Q 10 values that produce the target pyloric rhythm at each temperature. As Caplan et al. (2014) and O'Leary and Marder (2016) we searched for values of Q 10i between 1 and 2 for conductances and between 1 and 4 for activation/inactivation time scales. We built a new objective function by evaluating the previous objective function (6) Here, E(G, T i ) is the score of the solutions of the set of conductances G at temperature T i and E crash is the total number of spikes in all three cells. This last term is introduced to enforce that when the temperature is close to 35°C the network stops working as in experiments. Evaluation of the objective functions (6) and (7) requires that the model is simulated for a number of seconds and this is the part of the procedure that requires most computing power. Longer simulations will provide better estimations for the burst frequency and duty cycle of the solutions, but it will linearly increase the time it takes to evaluate them. If the simulations are shorter, evaluations of the objective function are faster but its minimization may be more di cult due to transient behaviors and its minima may not correspond to stable pyloric rhythms.
Because we only target the activity at 3 temperature values, it is not guaranteed that the models will behave properly for temperatures in between. We found that if the system meets the target activity at the control temperatures, in the vast majority of the cases, the models also displayed the correct activity for temperatures in between, but there were exceptions. This means that even if a solution achieves low scores, we are still required to screen the temperature values in between to make sure it is indeed temperature compensated over the full range. Increasing the number of control temperatures makes the estimation procedure slower and we met a good balance in 4 control temperatures plus screening as described.

C. Membrane potential and current shares distributions
The membrane potential distributions and current shares distributions were computed for 101 values of temperature between 10°C 25°C. For each temperature, the models were simulated from identical initial conditions for 30 seconds with high numerical resolution (dt = 0.001). The distributions were computed using the last 10 seconds of each simulation. The number of sam-ples of the distributions at each temperature is 5 ⇥ 10 5 .

D. Classification scheme
We classified the activity of each cell into several categories by inspecting their spiking patterns and features of their waveforms. For this we first compute the spike times and their ISI distributions. If the coe cient of variation CV ISI < 0.1 we declare the activity as spiking. We then measure the lag between the spike onset and the crossing times with a slow wave threshold at 50mV . If this < 20msec we label the trace as tonic spiking and else, we label it as single-spike bursting. If CV ISI > 0.1 then we ask if the cell is bursting. For this we group spikes into bursts using a temporal threshold of 200msec and compute the distributions of duty cycle dc and instantaneous burst frequency f b . If CV dc < 1, CV f b < 0.1 and mean(dc) > 0 we label the trace as regular bursting. If mean(dc) > 0 and either CV dc > 1 or CV f b > 0.1 we label the trace as irregular bursting. If mean(dc) = 0 we label the trace as irregular spiker. If none of these conditions were met we label the trace as other and finally, if the total number of ISI values is not greater than 1 we label the trace as quiescent. The network state was classified as triphasic or not triphasic. We declared the network (NET) activity as triphasic if the frequency of the cells was similar and they fired in the right order, and not triphasic elsewhere.

E. Parameters used in this study
Each model was assigned a 6 character name. Here we list which model was used in each figure. Code to simulate the network and the parameters of each model are supplemental to this work. V. Supplemental Figures  FIG. 1 Supplemental. Equivalent phenotype for di↵erent sets of maximal conductances and temperature sensitivities. For each set of maximal conductances we studied, it was possible to find multiple sets of Q10 that yielded temperature compensated solutions. A Burst frequency (top), duty cycle (middle), and phases (bottom) over the working temperature range for three di↵erent models (indicated in colors). B Values of the maximal conductances (in nS) and Q10 for each model .  FIG. 11 Supplemental. Distributions of current shares over temperature. A Traces and membrane potential distribution at two temperatures (top: 10°C, bottom: 25°C). Traces were colored using log(p(V )) as a color scale. In this way, di↵erent portions of the cycle appear in di↵erent colors: blue for spikes, yellow for sub-threshold activity and brown and pink for nearthreshold activity. B-C The panels show the percent contribution of each current to the total inward (and outward) currents, over the working temperature range (x-axis, 10°C 25°C)). The colors indicate the portions of the cycle at which each current contributes by each amount (y-axis).