Temperature compensation in a small rhythmic circuit

Abstract
Introduction
Results
Discussion
Materials and methods
Data availability
References
Article and author information
Metrics

Abstract

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.

Introduction

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. Temperature influences all biological processes, to a greater or lesser degree. This poses an inherent difficulty for neuronal signaling: if the currents involved in neuronal and network dynamics are differentially temperature-dependent a system that is well-tuned to work at one temperature may not function at a different temperature (Caplan et al., 2014; O'Leary and Marder, 2016; Tang et al., 2010; Tang et al., 2012). Nonetheless, many ectothermic animals have neurons and circuits that 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 effects of temperature on channel function and neuronal activity have been studied extensively for many years (Taylor and Kerkut, 1958). Increasing temperature generally results in increases of channel maximal conductances and faster activation/inactivation rates (Frankenhaeuser and Moore, 1963). But, importantly ionic channels of different types are affected by temperature to different extents: each of these processes has a different $Q_{10}$ (Schauf, 1973; Kukita, 1982; Ruff, 1999; Tang et al., 2010; Ranjan et al., 2019). The effect of temperature on neuronal intrinsic excitability such as voltage and current thresholds can be diverse (Sjodin and Mullins, 1958; Guttman, 1962; Guttman, 1966; Fitzhugh, 1966). For example, identified neurons in locust showed reversible changes in spike amplitude and duration of spikes as temperature varies from 18 to 35°C (Heitler et al., 1977). Studies in the jump neural circuit of grasshopers showed that neural excitability is differently affected by temperature across neuronal types (Abrams and Pearson, 1982).

There are also examples of neuronal and circuit processes that are relatively temperature-compensated. The frequency/current (f/I) curves in both molluscan and locust neurons are not substantially affected by small temperature changes (6–8°C) (Connor, 1975; Heitler et al., 1977). The f/I curves of auditory sensory neurons in grasshoppers remain largely unaffected by changes in temperature between 21°C and 29°C (Roemschied et al., 2014). Temperature compensation also takes place across 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 their intrinsic excitability properties are differentially modified by temperature.

Here, we explore these issues using computational models of the pyloric network—a subnetwork within the stomatogastric ganglion (STG) of crustaceans (Marder and Bucher, 2007; Maynard, 1972). The pyloric rhythm is a triphasic motor pattern that consists of bursts of action potentials in a specific sequence. This behavior is robust and stable, and the cells and their connections are well-characterized. Robustness to temperature in the pyloric network has been explored experimentally by Tang et al., 2010; Tang et al., 2012; Soofi et al., 2014 and Haddad and Marder, 2018. 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} \approx 2$ ), but their duty cycle (the burst duration in units of the period) stays approximately constant (Rinberg et al., 2013).

Temperature robustness was explored in computational models of the pacemaking kernel (Soto-Treviño 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 different 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 reproduce much of the experimental phenomenology previously reported (Tang et al., 2010; Tang et al., 2012; Soofi 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 different models and found that in all cases, the contributions of the currents to the activity can be significantly different across temperatures. The currents are not simply scaled up but instead become reorganized: a current that is important for burst termination at 10°C may no longer play that role at 25°C. Because the contribution of a given current to a neuronal process can be replaced by another at different temperatures, deletion or blockade of a current can produce qualitatively different effects at high and low temperatures. These results provide a plausible hypothesis for why interactions between temperature and a second perturbation can be observed experimentally (Haddad and Marder, 2018; Ratliff et al., 2018).

Results

Duty cycle and phase maintenance in model pyloric networks

The triphasic pyloric rhythm is produced by the periodic sequential activation of the pyloric dilator (PD) neurons, which are electrically coupled to 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. The model network consists of three cells, each modeled by a single compartment with eight currents as in previous studies (Golowasch and Marder, 1992; 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 $A B - P D$ (here we refer to this compartment as $P D$ ). Following Liu et al., 1998, each neuron has a sodium current, $I_{N a}$ ; transient and slow calcium currents, $I_{C a T}$ and $I_{C a S}$ ; a transient potassium current, $I_{A}$ ; a calcium-dependent potassium current, $I_{K C a}$ ; a delayed rectifier potassium current, $I_{K d}$ ; a hyperpolarization-activated inward current, $I_{H}$ ; and a leak current $I_{l e a k}$ . 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$ , $L P$ and $P Y$ cells. This pattern remains approximately periodic so we can measure the phases of the cycle at which each burst begins and terminates using the burst start of the $P D$ cell as reference (indicated as $B_{R E F}$ in magenta in Figure 1A), as indicated by the color labels.

Figure 1 with 3 supplements see all

Download asset Open asset

Effects of temperature in a model pyloric network.

(A) Schematic diagram of the model pyloric network in Prinz et al., 2004. The three groups interact via seven inhibitory chemical synapses. The red synapses are cholinergic from the PD neurons and all others are glutamatergic. 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$ - $L P$ - $P Y$ . (B) Activity of a temperature robust model network at 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) Top: average burst frequency of each cell over temperature (values are nearly identical so dots overlap). Middle: average duty cycle of each cell (cell type indicated in colors). Bottom: average phases of the cycle at which bursts begin and terminate using the start of the PD burst as reference (indicated by label $B_{R E F}$ in magenta). As temperature increases these phases remain approximately constant. The panels show average values over 30 s for 16 values of temperature between 10°C - 25°C.

Temperature sensitivity was modeled for each conductance as in Caplan et al., 2014 and O'Leary and Marder, 2016 (see Materials and methods). We increased each conductance by an Arrhenius type factor $R (T) = Q_{10}^{\frac{T - T_{r e f}}{10}}$ (with $Q_{10} \geq 1$ and $T_{r e f}$ = 10°C) and also increased the rates of the channel kinetics $τ$ in similar fashion. We assume that $Q_{10_{G}} \in [1, 2]$ for the maximal conductances and $Q_{10_{τ}} \in [1, 4]$ for the timescales as these ranges are consistent with experimental measurements of these quantities (Tang et al., 2010). The $Q_{10}$ values are a property of the channel proteins, so in a given model we use the same $Q_{10}$ for a given current (and synapse) type in all cells. Different values of the maximal conductances $G$ can produce an innumerable variety of activities. The space of solutions of the model—the different 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 (31 in total, 8 intrinsic conductances ×3 cells +7 synaptic conductances). The picture is more complicated as we consider the effect of temperature because it affects both the conductances and the time scales of the different processes (14 in total: 12 intrinsic, 2 synaptic). Therefore in the models, changes in temperature correspond to a coordinated change or a path in a 45-dimensional parameter space. The complexity of this problem led us to confine this 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 two- to three-fold factor over a 15°C range (Tang et al., 2010). For this to happen in the model, there should be many regions in the parameter space for which the activity is triphasic at different frequencies and paths in the 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 31 maximal conductances $G$ and 24 $Q_{10}$ values (see Materials and methods). 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 Materials and methods). 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 $1 H z$ , duty cycle ≈ 20% for $P D$ , and duty cycle ≈ 20% for $L P$ 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 robust pyloric rhythm.

Figure 1B shows the activity of one example network at several temperatures, similar to 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 simulated the model at several temperatures over the working range (10°C to 25°C) and computed the average duty cycle, frequency and phases of the bursts. While the burst frequency increases with $Q_{10} \approx 2$ , the duty cycle of the cells, and the phases of the onsets/offsets of activity, remain approximately constant. Figure 1—figure supplement 1A shows the same analysis performed on three different models. In all cases, the duty cycle and the phases stay approximately constant while the frequency increases. Figure 1—figure supplement 1B shows the values of the maximal conductances $G$ and the $Q_{10}$ values for each of these models.

While the duty cycle remains approximately constant on average, the spiking patterns of the cells show discernible differences across temperatures. We subjected the models to temperature ramps from 10°C to 25°C over 60 min and recorded the spike times (Figure 1—figure supplement 2). As temperature is increased the spiking patterns change in complicated ways and the largest $I S I$ decreases monotonically, consistent with the overall increase in bursting frequency of the cells. Different values of temperature result in slightly different spiking patterns which in turn result in different 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 (Figure 1—figure supplement 2D, blue box in PD cell). There are temperatures for which the duty cycle takes two values preferentially (Figure 1—figure supplement 2D, pink box in $P D$ cell) and there are temperatures for which the duty cycle differs noticeably from burst to burst. Although in all models, the average duty cycle and phases stay approximately constant, as temperature is changed the precise way in which the spiking patterns change show marked differences across models: different maximal conductances and different $Q_{10}$ values produce different patterns. Figure 1—figure supplement 3 shows the temperature dependence of the duty cycle in the $L P$ cell over the working range for all 36 models and provides a notion of how variable the bursting patterns of the cells can be across temperatures and individuals.

Changes in membrane potential over temperature

Experimental studies on other systems reported measurable effects 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 also 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 101 values of temperature. Temperature changes features of the voltage waveforms, resulting in changes in these distributions. Figure 2 shows example traces and the result of this analysis for the $L P$ cell in one model. The membrane potential changes with temperature in several ways: the spike amplitude decreases and the spike threshold depolarizes as temperature increases. 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 distribution on the right. The distributions permit visualizing changes in the waveform as the control parameter is changed (Alonso and Marder, 2019). The effect of temperature on membrane potential is not consistent across models and depends on the precise values of conductances and temperature sensitivities. Figure 2—figure supplement 1 shows the membrane potential distributions of each cell for three example models. In all cases, the waveforms show visible differences across temperatures and there is considerable variability in the amplitudes of the oscillations across models. In some models, the peak voltage of spikes, as indicated by the upper envelope of the distributions, decreases (left), increases (middle), or remains relatively constant (right). In addition, we found that cells can either become more depolarized or more hyperpolarized as temperature increases.

Figure 2 with 1 supplement see all

Download asset Open asset

Changes in membrane potential over temperature.

(left) Representative traces of the $L P$ cell in one model at different temperatures. (right) Membrane potential distribution at each temperature. The gray scale indicates the proportion of time the cell spends at that potential. The distribution facilitates inspecting how features such as the total amplitude of the oscillations and spiking thresholds change with temperature.

Dynamics of the currents at different temperatures

Neuronal activity is governed by currents that result from the precise activation and inactivation of ion channels with different kinetics. Experimental access to these quantities is limited because it is hard to measure currents individually without blocking all other currents, and thus changing the activity. Here, we employ models to explore how these currents are differentially 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 compared the dynamics of the currents in two model $P D$ neurons with different maximal conductances $G$ and temperature sensitivities $Q_{10}$ by direct inspection of the currents’ time series. Figure 3A shows the time series of each intrinsic current in $P D$ at 10°C and 25°C for one model. The corresponding membrane potential activity $V$ is shown in the blue traces on top. The red horizontal dashed lines indicate the peak amplitude of the currents at control temperature (10°C) to visualize whether it increases or decreases at 25°C. This is indicated by the arrows in each row (in colors: ⇑ red increases, ⇓ green decreases). In this example, the $N a$ current increases its amplitude by more than two fold over the 15°C range. Despite the fact that $g C a T$ increases with temperature, the amplitude of the current decreases to about half its control value at 10°C. Similar observations can be made for other currents (Figure 3A). Note that all $Q_{10}$ values in these models are greater than 1. Therefore, when the current decreases it is because the relative timing of the currents’ activation is changed by the neurons’ firing properties. The way currents are modified by temperature is different across models. Figure 3B shows the currents in the PD cell at two temperatures for another model. In this case, the $N a$ current decreases by a small amount despite the fact that the maximal conductance and timescales of this current increase with temperature. The $K C a$ current which in Figure 3A grows by a factor of six, is moderately decreased by temperature in Figure 3B.

Figure 3

Download asset Open asset

Dynamics of the currents at two temperatures.

Currents of the $P D$ cell for two different models (**A-B**). Temperature increases the peak amplitudes of some currents and it decreases it for others, as indicated by the ⇑ and ⇓ arrows in red and green. The way these currents are modified can vary substantially across models.

Because the total number of currents in the circuits is 31, it becomes cumbersome to compare them separately 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 the currents in a compact fashion. Figure 4 shows the currentscapes of each cell in one model at 10°C and 25°C. 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 on a logarithmic scale at the top and bottom. The currentscapes are noticeably different at 10°C and 25°C indicating that the currents contribute differentially across temperatures. For example, in the $P D$ cell the $L e a k$ 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 evenly distributed during the burst at 10°C but at 25°C its contribution is visibly larger at the beginning of the burst. In all cells, but most notably in the $L P$ and $P Y$ neurons in this example model, temperature substantially increases the contribution of the $K C a$ current toward the end of the bursts.

Figure 4

Download asset Open asset

Currentscapes at different temperatures.

The currentscapes use colors to display the percent contribution of each current type to the total inward and outward currents over time. The filled lines at the top and bottom indicate the total inward/outward currents in logarithmic scale. The panel shows the currentscapes for the three cells in the model at two temperatures. The panels show 2 periods $τ$ of the oscillation.

The currentscapes show that in all models and at all temperatures, there are periods of time when the activity is dominated by one or two currents, and periods of time when several currents contribute to it by similar amounts. During spike production, the activity is first driven by $I_{N a}$ and later by $I_{K d}$ . During these periods, the total inward and outward currents differ by orders of magnitude and are entirely dominated by one current type. Similarly, when the cells are most hyperpolarized the dominant currents are $I_{H}$ and the synaptic currents $I_{S y n}$ . There are periods of time—such as the beginnings and ends of bursts—when the total inward and outward currents 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.

Figure 5 shows 20 ms after the burst in a $P D$ neuron over temperatures between 10°C and 25°C. Almost all of the currents change over this temperature range, most notably there is almost a tradeoff between the contributions of the $A$ and $K C a$ currents, with $A$ decreasing smoothly and $K C a$ increasing smoothly. The specifics of how the contributions change is different across models, but in all cases there is a smooth transition between different relative fractions of currents. This allows the neuron to slide smoothly through changes in burst termination mechanisms.

Figure 5

Download asset Open asset

Current contributions at the end of bursts across temperature.

The figure shows the currentscapes of the PD cell for the 20 ms following burst termination, and how these change over temperature (same model as in Figure 7A).

Different compensation mechanisms fail in different ways. Figure 6 shows full currentscapes for two network models as they ‘crash’ at high temperature. Note that the complete loss of activity in the $P D$ neurons results from different currents being active in the two examples. Further examples of high temperature ‘crashes’ are shown in Figure 6—figure supplement 1. These examples show that the networks are crashing by different mechanisms in each case.

Figure 6 with 1 supplement see all

Download asset Open asset

Currentscapes reveal that ‘crashes’ occur by different mechanisms.

Currentscapes of two networks at 10°C and 25°C showing phase compensation, and ‘crashes’ at 30°C.

High temperatures produce disordered circuits or ‘crashes’

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). At temperatures higher than 25°C some models do not display irregular behavior but just become quiescent, or produce subthreshold oscillations. Other models display irregular and seemingly aperiodic regimes. As expected, irregular regimes typically take place at values of the temperature near a transition between qualitatively different patterns of activity. Figure 7 shows the activity of five models at temperatures higher than their working regimes. The left panels show 30 s of data and the right panels show 5 s. In all cases, there are timescales in the behavior that are much longer than the period of the rhythm at 10°C (≈ 1 s). The dynamics of the models in these regimes is daunting and their 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).

Figure 7

Download asset Open asset

Disordered circuits at high temperatures.

The models become dysfunctional—or crash—at temperatures higher than 25°C and do so in different ways. (**A-E**) Membrane potential of five models. The traces on the *left* show 30 s of data and the traces on the *right* show an expanded trace of 5 s. Note the emergence of time scales much longer than that of the original rhythm. Temperatures in the simulations are (**A–E**): 34°C, 34°C, 33°C, 29°C, 30°C.

Hysteresis and multistability

It is natural to ask whether these models show hysteresis when the temperature is returned to control levels. Moreover, multistability is also often observed in conductance-based models of neuronal activity and biological neurons (Cymbalyuk and Shilnikov, 2005; Malashchenko et al., 2011; Marin et al., 2013; Lechner et al., 1996; Paydarfar et al., 2006), and the network studied here also exhibits such features. We explored hysteresis in the models by increasing temperature linearly from 10°C to 35°C in 30 min and then decreasing it back to 10°C at the same rate. Figure 8 shows the spiking patterns for two $P D$ neurons during such temperature ramps (schematized on top). In Figure 8A the ISI distributions show bursting activity up to a critical temperature at which the model becomes quiescent. The quiescent state remains stable as temperature continues to increase and eventually recovers activity during the decreasing ramp, but at a lower temperature than when it stopped spiking during the increasing ramp. This implies that there is a range of temperatures where the model can be either quiescent or spiking, and thus is multistable. In most models (34 of 36), the spiking patterns over the working temperature range are similar across the up and down ramps but are not identical, suggesting there are multiple stable attractors that correspond to the triphasic rhythm but differ slightly in their spiking patterns (Figure 1—figure supplement 2E). This is consistent with the findings in Tang et al., 2010 where temperature was ramped slowly and no significant differences were found when ramping up and down in temperature.

Figure 8

Download asset Open asset

Hysteresis and multistability.

Temperature was increased from 10°C to 35°C and then decreased back to 10°C symmetrically. (Top) Schematic representation of the temperature ramp. Panels A and B show the ISI distributions of the $P D$ cell over time for two different models. (A) The cell ceases to produce spikes before 600 s, which is also before reaching 35°C, and remains quiescent until temperature is ramped down. The temperature at which the cell resumes spiking is ≈ 2°C lower than the temperature at which it became quiescent, indicating that the system is multi-stable over that temperature range. (B) Hysteresis is more evident in this model. The spiking patterns during the down ramp differ visibly from those during the up ramp over a wide temperature range.

We also found models with pathological stable attractors that coexist with the pyloric rhythm at working temperatures. Figure 8B shows a model in which hysteresis is salient. At high temperatures (times between 25 and 35 min), this model produces $I S I$ values as long as ≈ 1 s and remains spiking at all temperatures. The spiking patterns differ noticeably on the up and down ramps. The spiking pattern of the $P D$ cell changes from regular bursting to a pathological state that remains stable until the model is ramped back down to ≈ 20°C. In our simulations, hysteresis mainly occurs at high temperatures indicating that states such as those in Figure 7 can be multistable. The extent such hysteresis can be revealed in the biological experiments depends on the details of the perturbation protocol, and its currently under study.

Response to perturbations at different temperatures

In all models, the dynamics of the currents are different at different temperatures and their contributions become reorganized in complex ways. The combined activity of all the variables in the model—the dynamical attractor—is therefore expected to have different stability properties at different temperatures. This means that an extreme perturbation can have qualitatively different effects at different temperatures. The responses of different models to extreme perturbations such as partially or completely removing a current can be diverse at any temperature (Alonso and Marder, 2019). While this is expected due to the complexity of these models, the fact that these responses can also be qualitatively different across temperatures indicates that changes in the current contributions are meaningful and affect the stability of the network. To shed light on this issue, we performed several simple perturbations.

Figure 9 shows one perturbation in two models. In Figure 9A, we explored the interaction between temperature and the effect of removing the $A$ -type K⁺ channel at two temperatures: the models were simulated for 20 s in control conditions and then $I_{A}$ was completely removed for the next 20 s. Figure 9A (top) shows that removing the $A$ current at 10°C has catastrophic consequences for the network activity which becomes 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 9B shows a similar pertubation in which we remove the $K C a$ current. This model displayed normal triphasic behavior at 10°C (top) and when the $K C a$ current was removed its activity slowed but the triphasic rhythm was maintained. The same perturbation at 25°C results in quiescence. In this model, the $K C a$ current is necessary for the activity at 25°C but not at 10°C.

Figure 9

Download asset Open asset

Response to perturbations at different temperatures.

The response of the models to extreme perturbations can be different at different temperatures. The figure shows the responses of models to the same perturbation at two different temperatures. (A) (top) Membrane potential over time at 10°C. The first 4 s correspond to the control condition with all currents intact. At $T = 4 s e c$ we removed the $A$ current. (bottom) Same perturbation performed at 25°C (same model as in Figure 4). (B) Same protocol as in A but removing the $K C a$ current in a different model (Figure 7A).

To further characterize the differential response to perturbations at different temperatures, we performed a second assay in which we gradually decreased a current from control (100%) to complete deletion (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 different states depending on the temperature at which this perturbation is performed. The responses are so diverse that a coarse classification scheme is sufficient 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 sufficient to tease apart the different regimes we observe at a coarse level (Materials and methods).

Figure 10A shows the result of classifying the responses of one model to gradually decreasing a conductance in five steps, at two different temperatures, for four conductances. We employed 5 × 4 response grids to summarize the effects of partially or completely removing a conductance at 10°C and 25°C. The rows in the grids correspond to different values of 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 color coded for condition as indicated by the labels in the figure. The top left panel in Figure 10A shows the responses of one model to removing the sodium current $I_{N a}$ . At 10°C, when $g N a \to 75 % g N a$ the $L P$ 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 triphasic. Removing the $K C a$ current in this model has little effect 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 $I_{K C a}$ . The same observation that the responses to perturbation are different at each temperature holds regardless of the perturbation type. In the cases of the $C a T$ and $C a S$ currents complete removal results in different non-triphasic activity at each temperature. The responses of the model are in general qualitatively different at different temperatures, and they are also diverse across models (not shown). The fact that the models ‘break-down’ in different ways is consistent with the observation that these network mechanisms are diverse. In all models, the responses to extreme perturbation are temperature dependent for most perturbation types. In Figure 10B we quantify how many models respond differently to complete removal of a current at 10°C and 25°C. 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 different ways at each temperature. The exceptions are the $C a T$ current for which most models respond in different ways, and the $H$ current where complete removal results in equivalent states (quiescence in this case) at both temperatures.

Figure 10

Download asset Open asset

Classification of responses at different temperatures.

The models’ responses to perturbations are diverse and are in general different at different temperatures. We classified the responses of the models at two temperatures for five values of gradual decrements of each current. (A) Response to removal of $N a$ , $K C a$ , $C a T$ and $C a S$ currents for one model. (B) The pie charts show the number of models for which the responses to complete removal (0%) are classified as the same (green) or different (blue) at 10°C and 25°C.

Discussion

Enzymes and ion channels show altered behavior when temperature is changed, as temperature influences their conformation. Yet many animals live successfully over wide temperature ranges. For example, North Atlantic lobsters and crabs routinely experience more than 25°C temperature swings during a year, and may see 8 - 10°C changes within a short time period (Haefner, 1977; Stehlik et al., 1991). This raises a series of fascinating questions of how behavior and circuit performance can be maintained while all biological molecules are, to a greater or lesser degree, differentially temperature dependent (Robertson and Money, 2012).

It is important to distinguish between altered activity that occurs rapidly in response to acute, or immediate changes in temperature and those changes that occur during slow adaptation and acclimation, processes that can take days, weeks or months. In this study, we have modeled the fast changes that we assume would occur ‘automatically’ as proteins change conformation and properties rapidly, and as are usually measured in their $Q_{10}$ ’s. This is to be contrasted with the slower, long-term effects that result in acclimation in which a variety of molecular processes are triggered that may in turn change the proteins present in the membrane. These mechanisms, and much adaptation to long-term changes in temperature depend on the activation of heat shock proteins, RNA editing and a whole array of molecular mechanisms (Garrett and Rosenthal, 2012; Rosenthal, 2015; Martin Anduaga et al., 2019). We will first discuss the issues that arise from the examination of robustness to acute temperature changes, as modeled here, and then offer some speculations for additional processes that may be important to enhance robustness in response to long-term temperature changes.

Unfortunately, there are relatively few data relevant to whether the $Q_{10}$ of a channel in a given animal changes in response to acclimation or seasonal conditions. Therefore, although it is often assumed that a $Q_{10}$ for a given protein is invariant over the lifetime of the animal this may not necessarily be the case. Another difficulty is that there are few, if any, studies of many ion channel $Q_{10}$ s in the same animal in the same temperature range. Computational modeling will not answer what the ‘real’ conductance and $Q_{10}$ combinations are.

The assumption that these quantities are stationary throughout the life of an animal (and/or the same across species) may not be warranted. For example, long-term acclimation over months resulting in RNA editing, could alter some $Q_{10}$ s. Moreover, it is possible that minor genetic variations in channel proteins exist in the population so that different animals in the population could show modestly different $Q_{10}$ s. Therefore, the point of this study is to ask generically how compensation occurs acutely assuming different temperature sensitivities for the processes, within the experimental range. For these reasons, in the absence of data, it is valuable to explore the $Q_{10}$ space to have on overview of how compensation takes place acutely at any point during the lifetime of an animal.

One of the most important and non-intuitive results from this work is seen in Figure 5, which demonstrates that a neuron can smoothly slide through mechanisms that govern a physiological process as temperature is changed. In the case of Figure 5, the repolarization of the last spike in the burst shows a large contribution from $I_{A}$ at 10°C but becomes small at 25°C while $I_{K C a}$ contributes little at 10°C and becomes dominant at 25°C. This shows that while both of these $K$ currents can contribute to spike and burst repolarization, their shares in these processes change, and the neuron is robust as it smoothly slides through these mechanisms. This smooth transition can occur precisely because these three K⁺ currents depend on voltage, time, and temperature differently.

We speculate that one advantage conferred by the expression of many different ion channels that differ in their properties is to provide many sets of such ‘sliding stability mechanisms’. Individual neurons in all species have large numbers of ion channel genes that encode proteins that give rise to the membrane currents that are usually measured and characterized in voltage-clamp. The model neurons used are loosely based on voltage-clamp measurements from crab and lobster stomatogastric ganglion (STG) neurons (Golowasch et al., 1992; Turrigiano et al., 1995) and have eight currents. That said, we now know that there are at least twenty to thirty and probably more voltage-gated ion channel genes in individual STG neurons (Northcutt et al., 2016). In some cases, several of these genes may contribute to one of the currents that we model here. For example, the delayed rectifier current modeled here may have contributions from at least three different genes, and $I_{A}$ , modeled here, also probably has contributions from at least two different genes. Given that these genes encode similar but not identical protein subunits, it is not unreasonable to expect that these proteins may respond differently to temperature, and that the actual composition of the activated $K$ currents may reflect a change in these contributions as a function of temperature. This again, could provide still another ‘sliding stability mechanism’, and it would be interesting in the future to create models using data from currents that have been measured in isolation, such as in Ranjan et al., 2019.

The same principle of sliding circuit mechanisms 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 conjoint action of a number of different 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 different cellular mechanisms.

Long-term mechanisms for temperature adaptation and acclimation may extend robustness by changing the functional properties of ion channel proteins. In biological systems, sustained exposure to temperature often activates heat shock proteins and has effects on many cellular and molecular functions (Sharp et al., 1999). Interestingly, changes in ion channel composition can result from RNA editing, or changes in splice variants (Johnson et al., 2011; Lin and Baines, 2015). It is easy to imagine that these kinds of molecular mechanisms can result in shifts in the temperatures at which neurons and circuits may function. Acclimation to high temperatures usually requires at least 3–4 weeks of sustained alteration of the environmental temperature, and this is enough time to replace a large subset of ion channel proteins in the membrane. This raises a series of fascinating questions for how homeostatic regulation of excitability may be preserved if the molecular mechanisms that are engaged in homeostasis are also influenced by temperature.

Hysteresis in response to perturbation can result from a variety of time-dependent mechanisms. Figure 8 shows examples of hysteresis that can occur in response to temperature ramps, even in models in which the effects of temperature on the properties of each ion channel are viewed as essentially instantaneous. This hysteresis can occur because at high temperatures (> 25°C) the network can produce several different patterns of activity or dynamical states (multistability). The activation and inactivation rates that characterize ion channel function are different for each of these states and therefore, as temperature is lowered, they lose stability in different ways. Of course, any long-term molecular mechanisms (not modeled here) that alter channel expression or splicing or phosphorylation, could produce long-lasting hysteresis.

In this work, we used a genetic algorithm (Alonso and Marder, 2019) to find temperature robust networks. It is important to reiterate that temperature robust neurons and networks are difficult to find doing random searches (Caplan et al., 2014). Therefore, despite the significant animal-to-animal differences in conductance densities seen across the population, randomly sampling conductance values is unlikely to result in successful solutions for individual animals. This argues that evolution has arrived at a series of biological mechanisms that give rise to successful temperature compensation. These are likely to include 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., 1999a; 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; Golowasch et al., 1999b; Golowasch et al., 1999a; Santin and Schulz, 2019; Golowasch, 2019). The temperature sensitivities of ion channels are dictated by their molecular structures, and for this reason it is tempting to assume that the temperature sensitivity of a given channel is similar throughout the lifetime of the animal. Nonetheless, as all channels are subject to splice variants and other mechanisms of molecular regulation, it is possible that environmental temperature results in changes in $Q_{10}$ that are consistent with acute robustness, even if there have been changes in conductance densities. This suggests that conductance densities and $Q_{10}$ regulation may be coordinately controlled in biological systems.

One of the most significant results of this paper is that temperature robustness of network function can take place in multiple different ways. In these models, both the membrane potential and the spiking patterns are affected by temperature changes and these changes appear to be different in every model we inspected. Across the models, different values of the maximal conductances G and temperature sensitivities $Q_{10}$ result in consistent differences 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., 2006; Schulz et al., 2007; Temporal et al., 2014; Northcutt et al., 2016) 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 different 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; Temporal 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 different compositions of channel densities—which produce similar pyloric rhythms—are differentially resilient to any given perturbation. By the same token, as currents change as a function of temperature we expect that a second global perturbation of an individual may have qualitatively different effects at different temperatures, as is illustrated with the models in Figure 10. The interaction between temperature and a second perturbation is the subject of recent experimental studies (Haddad and Marder, 2018; Ratliff et al., 2018). These studies are consistent with the interpretation that different current configurations have different stability properties and that temperature changes these configurations.

One of the general take-home lessons from this work is that models can be resistant to a pharmacological perturbation at one temperature but the same model can be sensitive to the same perturbation at another temperature (Figure 10), and that different individuals may differ in their responses to the same perturbation. The potential implications of this for therapeutics are evident. But also evident, is that the search for therapeutic agents can be fraught if the assay itself is differentially sensitive to temperature or another perturbation. Moreover, the ‘sliding mechanisms’ provided by multiple ion channel genes and their molecular regulation suggests that over time therapeutic sensitivity to a given pharmacological intervention could be quite variable. These findings may contribute to our understanding of individual variable therapeutic efficacies within the population and across time within an individual.

Materials and methods

The model

Request a detailed protocol

The activity of the cells was modeled using single-compartment 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_{N a}$ ; transient and slow calcium currents, $I_{C a T}$ and $I_{C a S}$ ; a transient potassium current, $I_{A}$ ; a calcium-dependent potassium current, $I_{K C a}$ ; a delayed rectifier potassium current, $I_{K d}$ ; a hyperpolarization-activated inward current, $I_{H}$ ; and a leak current $I_{l e a k}$ . The number of state variables per cell is 13. The units for voltage are $m V$ , the conductances are expressed in $μ S$ and currents in $n A$ . The pyloric network consists of three cells with the same ion channel types but different conductance densities. The interactions in the network consist of seven chemical synapses and are similar to Prinz et al., 2004. The synaptic current is given by $I_{s} = g_{s} s (V_{p o s t} - E_{s})$ , where $g_{s}$ is the synapse strength, $V_{p o s t}$ 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

\frac{d s}{d t} = \frac{s_{\infty} (V_{p r e}) - s}{τ_{r} + τ_{s}}

with,

s_{\infty} (V_{p r e}) = \frac{1}{1 + e x p ((V_{t h} - V_{p r e}) / Δ)},

and

τ_{s} = \frac{1 - s_{\infty} (V_{p r e})}{k_{-}} .

These equations are identical to Prinz et al., 2004 except for the inclusion of a bound for the timescale of activation $τ_{r} = 20 m s e c$ (we want to avoid the case that as $s_{\infty} \to 1$ then $τ_{s} \to 0$ and $\dot{s} \to \infty$ ). All other parameters (except $g_{s}$ ) are identical to Prinz et al., 2004. Following Prinz et al., 2004, we set $E_{s} = - 70 m V$ and $k_{-} = \frac{1}{40} m s$ for glutamatergic synapses, and $E_{s} = - 80 m V$ and $k_{-} = \frac{1}{100} m s e c$ for cholinergic synapses. We set $V_{t h} = - 35 m V$ and $Δ = 5 m V$ for both synapse types.

Temperature effects 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 $τ_{m_{i}}$ and $τ_{h_{i}}$ , the maximal conductances $g_{i}$ , the time constants of calcium buffering $τ_{C a}$ , 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 $τ_{i}$ by $R_{i} {(T)}^{- 1} τ_{i}$ where

R_{i} (T) = Q_{10_{i}}^{\frac{T - T_{r e f}}{10}} .

Here $T$ is the temperature, $T_{r e f} = 10^{\circ}$ is the reference temperature and $Q_{10_{i}}$ is defined as the fold change per 10°C from the reference temperature ( $i$ indicates the process type). Finally, the calcium reversal potential $E_{C a}$ 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 9: the 8 maximal conductances plus the calcium time constant $τ_{C a}$ . 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 11 $Q_{10}$ values for the time scales of each intrinsic process (not all currents inactivate), plus the $Q_{10}$ of the calcium time constant. 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 $d t = 0.05 m s e c$ (Press et al., 1988). We used the same set of initial conditions for all simulations in this work $V = - 51 m V$ , $m_{i}, h_{i} = 0$ and $[C a^{+ 2}] = 5 μ M$ .

Finding parameters

Request a detailed protocol

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_{t g}$ ) and duty cycle ( $d c_{t g}$ ). The function uses thresholds to obtain temporal information of the waveform, such as spike times, and then uses it 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 robust networks.

For any given set of conductances we simulated the model for 20 s and dropped the first 10 s to minimize the effects of transient activity. We then computed the average ( $< >$ ) burst frequency $< f_{b} >$ , the average duty cycle $< d c >$ , the number of crossings with a slow wave threshold $#_{s w} = - 50 m V$ , the number of bursts $#_{b}$ , and the average lags between bursts $< Δ_{P D - L P} >$ and $< Δ_{P D - 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 $s t d ({f_{b}}) \geq< f_{b} > \times 0.1$ or $s t d ({d c}) \geq< d c > \times 0.2$ . If a solution is not discarded we can use these quantities to measure how close it is to a target behavior,

\begin{array}{ll} E_{f} & = \sum_{i = cell} (f_{t g} - < f_{b} >_{i})^{2} \\ E_{d c} & = \sum_{i = cell} (d c_{t g} - < d c >_{i})^{2} \\ E_{s w} & = \sum_{i = cell} (#_{s w} - #_{b})^{2} \\ E_{p h} & = (Δ_{P D - L P t g} - \frac{< Δ_{P D - L P} >}{τ_{b}})^{2} + (Δ_{P D - P Y t g} - \frac{< Δ_{P D - P Y} >}{τ_{b}})^{2} . \end{array}

Here $E_{f}$ measures the mismatch of the bursting frequency of each cell with a target frequency $f_{t g}$ and $E_{d c}$ accounts for the duty cycle. $E_{s w}$ measures the difference between the number of bursts and the number of crossings with the slow wave threshold $t_{s w} = - 50 m V$ (we ask that $#_{s w} = #_{b}$ ). $E_{p h}$ compares the lags between bursts (in units of the bursting period $τ_{b} = \frac{1}{< f_{b} >}$ ) to a target lag $Δ_{t g}$ . These measures are discussed in more detail in Alonso and Marder, 2019.

Let $G$ denote a set of maximal conductances (and $τ_{C a}$ ), we can then define an objective function

E (G) = α E_{f} + β E_{d c} + γ E_{s w} + η E_{p h},

where the coefficients $(α, β, γ, η)$ define how the different sources of penalties are weighted. In this work we used $α = 10$ , $β = 1000$ , $γ = 1$ , $η = 10$ . These weights were found by trial and error and kept fixed. The penalties $E_{i}$ were calculated using $T = 10$ secs with $d t = 0.05$ msecs. The target control behavior was defined as having all cells bursting with 20% duty cycle for $P D$ ( $d c_{t g_{P D}} = 0.2$ ) and 25% for the LP and PY cells ( $d c_{t g_{L P, P Y}} = 0.25$ ). The lag between bursts was targeted to be $Δ_{P D - L P} t g = 0.5$ and $Δ_{P D - P Y} t g = 0.75$ . The target burst frequency of all cells was set to $f_{t g} = 1 H z$ (Bucher et al., 2006; Tang et al., 2010; Hamood et al., 2015).

Minimization of the objective function produces sets of maximal conductances $G$ for which the resulting circuit activity mimics the pyloric rhythm. The minimization was performed over a search space of allowed values listed here: for each cell we searched for $g_{N a} \in [0, 10^{3}]$ ([ $μ S$ ]), $g_{C a T} \in [0, 10^{2}]$ , $g_{C a S} \in [0, 10^{2}]$ , $g_{A} \in [0, 10^{2}]$ , $g_{K C a} \in [0, 10^{3}]$ , $g_{K d} \in [0, 10^{2}]$ , $g_{H} \in [0, 10^{2}]$ , $g_{L} \in [0, 10]$ , $τ_{C a} \in [0, 2 \times 10^{3}]$ (msec). All synaptic conductances were searched in the range $g_{s y n} \in [0, 5 \times 10^{- 2}]$ ([ $μ S$ ]). The choice of search space contains our knowledge that some conductances like $g N a$ and $g K d$ are significantly larger than $g_{l e a k}$ or the synaptic conductances. This choice makes optimizations faster and produces diverse solutions. We discretized our search space by taking 10³ equally spaced values for each parameter. We minimized the objective function (Equation 6) using a custom genetic algorithm (Holland, 1992) on a desktop computer with 32 cores. For each network, we ran the algorithm using 1000 individuals taken at random in the search space for 10,000 generations. Finding a successful population takes several hours.

Temperature robustness is achieved in our models by searching sets of $Q_{10}$ values that produce the target pyloric rhythm at each temperature. As in Caplan et al., 2014 and O'Leary and Marder, 2016 we searched for values of $Q_{10_{i}}$ 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 (Equation 6) over a set of control temperatures $T_{i}$ taken at 15°C, 20°C, 25°C and 35°C,

E_{C} (Q_{10}) = E (G, T = 15^{\circ}) + E (G, T = 20^{\circ}) + E (G, T = 25^{\circ}) + E_{c r a s h} (G, T = 35^{\circ}) .

Here, $E (G, T_{i})$ is the score of the solutions of the set of conductances $G$ at temperature $T_{i}$ and $E_{c r a s h}$ 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 (Equation 6) and (Equation 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 difficult due to transient behaviors, and its minima may not correspond to stable pyloric rhythms.

Because we only target the activity at three 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 (2 in 40 cases). 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 with 4 control temperatures. All models were subject to the same temperature perturbations and temperature compensation was confirmed in all cases by recording spikes and phases of bursting (Figures 1 and 3), membrane potential (Figure 2), and currents’ dynamics (Figures 7, 8 and 9).

Membrane potential and current shares distributions

Request a detailed protocol

The membrane potential distributions were computed for 101 values of temperature between 10°C and 25°C. For each temperature, the models were simulated from identical initial conditions for 30 s with high numerical resolution ( $d t = 0.001$ ). The distributions were computed using the last 10 s of each simulation. The number of samples of the distributions at each temperature is 5 × 10⁵.

Classification scheme

Request a detailed protocol

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 coefficient of variation $C V_{I S I} < 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 $- 50 m V$ . If this $δ < 20 m s e c$ we label the trace as tonic spiking and else, we label it as single-spike bursting. If $C V_{I S I} > 0.1$ then we ask if the cell is bursting. For this we group spikes into bursts using a temporal threshold of 200 ms and compute the distributions of duty cycle $d c$ and instantaneous burst frequency $f_{b}$ . If $C V_{d c} < 1$ , $C V_{f_{b}} < 0.1$ and $m e a n (d c) > 0$ we label the trace as regular bursting. If $m e a n (d c) > 0$ and either $C V_{d c} > 1$ or $C V_{f_{b}} > 0.1$ we label the trace as irregular bursting. If $m e a n (d c) = 0$ we label the trace as irregular spiker. If none of these conditions is 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 (within 5%) and they fired in the right order, and not triphasic elsewhere.

Parameters used in this study

Request a detailed protocol

Each model was assigned a six-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.

Figure 1 4 × 6W66. Figure 1—figure supplement 1: from left to right, WWZ3CN, FGFKQS, 71G6LA. Figure 1—figure supplement 2: PLKTEM. Figure 1—figure supplement 3: all models. Figure 2: G7P2WE. Figure 2—figure supplement 1: from left to right, G7P2WE, 4 × 6W66, MAXTTP. Figure 3: 4 × 6W66, SJR46Y.Figure 4: WWZ3CN. Figure 5: FGFKQS.Figure 6: FGFKQS, 71G6LA. Figure 6—figure supplement 1: in lexicographical order, B7SDOL, J9U8SQ, 4R2719, R3VIEW, Q8FBN1, RMZ8KZ, EF84RN, SJR46Y. Figure 7: (top to bottom), J9U8SQ, G81PUL, G7EOZ8, SAXT8Y, KNMLBC. Figure 8: OSHNMD, H5AIZR. Figure 9: WWZ3CN, FGFKQS. Figure 10: 4 × 6W66.

Data availability

Code and parameters to reproduce the results in this work was uploaded to Dryad, http://doi.org/10.5061/dryad.m31b1h7.

The following data sets were generated

1. Alonso LM
2. Marder E
(2020) Dryad Digital Repository
Data from: Temperature compensation in a small rhythmic circuit.

https://doi.org/10.5061/dryad.m31b1h7

References

1. Abrams TW
2. Pearson KG
(1982) Effects of temperature on identified central neurons that control jumping in the grasshopper
The Journal of Neuroscience 2:1538–1553.

https://doi.org/10.1523/JNEUROSCI.02-11-01538.1982
- PubMed
- Google Scholar
1. Alonso LM
2. Marder E
(2019) Visualization of currents in neural models with similar behavior and different conductance densities
eLife 8:e42722.

https://doi.org/10.7554/eLife.42722
- PubMed
- Google Scholar
(2006) Central pattern generating neurons simultaneously express fast and slow rhythmic activities in the stomatogastric ganglion
Journal of Neurophysiology 95:3617–3632.

https://doi.org/10.1152/jn.00004.2006
- PubMed
- Google Scholar
(1992) Mathematical model of an identified stomatogastric ganglion neuron
Journal of Neurophysiology 67:332–340.

https://doi.org/10.1152/jn.1992.67.2.332
- PubMed
- Google Scholar
(2014) Many parameter sets in a multicompartment model oscillator are robust to temperature perturbations
The Journal of Neuroscience 34:4963–4975.

https://doi.org/10.1523/JNEUROSCI.0280-14.2014
- PubMed
- Google Scholar
1. Connor JA
(1975) Neural repetitive firing: a comparative study of membrane properties of crustacean walking leg axons
Journal of Neurophysiology 38:922–932.

https://doi.org/10.1152/jn.1975.38.4.922
- PubMed
- Google Scholar
1. Cymbalyuk G
2. Shilnikov A
(2005) Coexistence of tonic spiking oscillations in a leech neuron model
Journal of Computational Neuroscience 18:255–263.

https://doi.org/10.1007/s10827-005-0354-7
- Google Scholar
1. Fitzhugh R
(1966) Theoretical effect of temperature on threshold in the Hodgkin-Huxley nerve model
The Journal of General Physiology 49:989–1005.

https://doi.org/10.1085/jgp.49.5.989
- PubMed
- Google Scholar
1. Frankenhaeuser B
2. Moore LE
(1963) The effect of temperature on the sodium and potassium permeability changes in myelinated nerve fibres of Xenopus laevis
The Journal of Physiology 169:431–437.

https://doi.org/10.1113/jphysiol.1963.sp007269
- PubMed
- Google Scholar
1. Garrett S
2. Rosenthal JJ
(2012) RNA editing underlies temperature adaptation in K+ channels from polar octopuses
Science 335:848–851.

https://doi.org/10.1126/science.1212795
- PubMed
- Google Scholar
(2009) Functional consequences of animal-to-animal variation in circuit parameters
Nature Neuroscience 12:1424–1430.

https://doi.org/10.1038/nn.2404
- PubMed
- Google Scholar
(2001) Global structure, robustness, and modulation of neuronal models
The Journal of Neuroscience 21:5229–5238.

https://doi.org/10.1523/JNEUROSCI.21-14-05229.2001
- Google Scholar
(1992) Contribution of individual ionic currents to activity of a model stomatogastric ganglion neuron
Journal of Neurophysiology 67:341–349.

https://doi.org/10.1152/jn.1992.67.2.341
- Google Scholar
(1999a) Activity-dependent regulation of potassium currents in an identified neuron of the stomatogastric ganglion of the crab Cancer borealis
Journal of Neuroscience 19:RC33.

https://doi.org/10.1523/JNEUROSCI.19-20-j0004.1999
- Google Scholar
(1999b) Network stability from activity-dependent regulation of neuronal conductances
Neural Computation 11:1079–1096.

https://doi.org/10.1162/089976699300016359
- PubMed
- Google Scholar
1. Golowasch J
(2019) Neuronal homeostasis: voltage brings it all together
Current Biology 29:R641–R644.

https://doi.org/10.1016/j.cub.2019.05.029
- PubMed
- Google Scholar
1. Golowasch J
2. Marder E
(1992) Ionic currents of the lateral pyloric neuron of the stomatogastric ganglion of the crab
Journal of Neurophysiology 67:318–331.

https://doi.org/10.1152/jn.1992.67.2.318
- PubMed
- Google Scholar
1. Guttman R
(1962) Effect of temperature on the potential and current thresholds of axon membrane
The Journal of General Physiology 46:257–266.

https://doi.org/10.1085/jgp.46.2.257
- PubMed
- Google Scholar
1. Guttman R
(1966) Temperature characteristics of excitation in space-clamped squid axons
The Journal of General Physiology 49:1007–1018.

https://doi.org/10.1085/jgp.49.5.1007
- PubMed
- Google Scholar
1. Haddad SA
2. Marder E
(2018) Circuit robustness to temperature perturbation is altered by neuromodulators
Neuron 100:609–623.

https://doi.org/10.1016/j.neuron.2018.08.035
- PubMed
- Google Scholar
1. Haefner PA
(1977) Aspects of the biology of the jonah crab, Cancer borealis Stimpson, 1859 in the mid-Atlantic Bight
Journal of Natural History 11:303–320.

https://doi.org/10.1080/00222937700770221
- Google Scholar
(2018) Two central pattern generators from the crab, Cancer borealis, respond robustly and differentially to extreme extracellular pH
eLife 7:e41877.

https://doi.org/10.7554/eLife.41877
- PubMed
- Google Scholar
(2015) Quantitative Reevaluation of the Effects of Short- and Long-Term Removal of Descending Modulatory Inputs on the Pyloric Rhythm of the Crab, Cancer borealis
Eneuro 2:ENEURO.0058-14.2015.

https://doi.org/10.1523/ENEURO.0058-14.2015
- Google Scholar
(1995) Dopamine modulation of transient potassium current evokes phase shifts in a central pattern generator network
The Journal of Neuroscience 15:342–358.

https://doi.org/10.1523/JNEUROSCI.15-01-00342.1995
- PubMed
- Google Scholar
(1977) The effects of temperature on the threshold of identified neurons in the locust
Journal of Comparative Physiology ? A 117:163–182.

https://doi.org/10.1007/BF00612785
- Google Scholar
1. Holland JH
(1992) Genetic algorithms
Scientific American 267:66–72.

https://doi.org/10.1038/scientificamerican0792-66
- Google Scholar
1. Hooper SL
(1997) Phase maintenance in the pyloric pattern of the lobster (Panulirus interruptus) stomatogastric ganglion
Journal of Computational Neuroscience 4:191–205.

https://doi.org/10.1023/a:1008822218061
- PubMed
- Google Scholar
(2011) Alternatively spliced domains interact to regulate BK potassium channel gating
PNAS 108:20784–20789.

https://doi.org/10.1073/pnas.1116795108
- Google Scholar
(1974) Temperature effects on resting potential and spike parameters of cat motoneurons
Experimental Brain Research 19:478–492.

https://doi.org/10.1007/BF00236112
- PubMed
- Google Scholar
1. Krupczynski P
2. Schuster S
(2013) Precision of archerfish C-starts is fully temperature compensated
Journal of Experimental Biology 216:3450–3460.

https://doi.org/10.1242/jeb.088856
- PubMed
- Google Scholar
1. Kukita F
(1982) Properties of sodium and potassium channels of the squid giant axon far below 0 degrees C
The Journal of Membrane Biology 68:151–160.

https://doi.org/10.1007/BF01872261
- PubMed
- Google Scholar
(1996) Bistability and its regulation by serotonin in the endogenously bursting neuron R15 in Aplysia
Journal of Neurophysiology 75:957–962.

https://doi.org/10.1152/jn.1996.75.2.957
- PubMed
- Google Scholar
(1993) Activity-dependent regulation of conductances in model neurons
Science 259:1915–1917.

https://doi.org/10.1126/science.8456317
- PubMed
- Google Scholar
1. Lin WH
2. Baines RA
(2015) Regulation of membrane excitability: a convergence on voltage-gated sodium conductance
Molecular Neurobiology 51:57–67.

https://doi.org/10.1007/s12035-014-8674-0
- PubMed
- Google Scholar
1. Liu Z
2. Golowasch J
3. Marder E
4. Abbott LF
(1998) A model neuron with activity-dependent conductances regulated by multiple calcium sensors
The Journal of Neuroscience 18:2309–2320.

https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
- PubMed
- Google Scholar
(2011) Six types of multistability in a neuronal model based on slow calcium current
PLOS ONE 6:e21782.

https://doi.org/10.1371/journal.pone.0021782
- PubMed
- Google Scholar
1. Manor Y
2. Nadim F
3. Abbott LF
4. Marder E
(1997) Temporal dynamics of graded synaptic transmission in the lobster stomatogastric ganglion
The Journal of Neuroscience 17:5610–5621.

https://doi.org/10.1523/JNEUROSCI.17-14-05610.1997
- PubMed
- Google Scholar
1. Marder E
2. Bucher D
(2007) Understanding circuit dynamics using the stomatogastric nervous system of lobsters and crabs
Annual Review of Physiology 69:291–316.

https://doi.org/10.1146/annurev.physiol.69.031905.161516
- PubMed
- Google Scholar
(2013) High prevalence of multistability of rest states and bursting in a database of a model neuron
PLOS Computational Biology 9:e1002930.

https://doi.org/10.1371/journal.pcbi.1002930
- PubMed
- Google Scholar
(2019) Thermosensitive alternative splicing senses and mediates temperature adaptation in Drosophila
eLife 8:e44642.

https://doi.org/10.7554/eLife.44642
- PubMed
- Google Scholar
1. Maynard DM
(1972) Simpler networks
Annals of the New York Academy of Sciences 193:59–72.

https://doi.org/10.1111/j.1749-6632.1972.tb27823.x
- PubMed
- Google Scholar
1. Nadim F
2. Manor Y
3. Kopell N
4. Marder E
(1999) Synaptic depression creates a switch that controls the frequency of an oscillatory circuit
PNAS 96:8206–8211.

https://doi.org/10.1073/pnas.96.14.8206
- PubMed
- Google Scholar
1. Nadim F
2. Manor Y
(2000) The role of short-term synaptic dynamics in motor control
Current Opinion in Neurobiology 10:683–690.

https://doi.org/10.1016/S0959-4388(00)00159-8
- PubMed
- Google Scholar
1. Northcutt AJ
2. Lett KM
3. Garcia VB
4. Diester CM
5. Lane BJ
6. Marder E
7. Schulz DJ
(2016) Deep sequencing of transcriptomes from the nervous systems of two decapod crustaceans to characterize genes important for neural circuit function and modulation
BMC Genomics 17:868.

https://doi.org/10.1186/s12864-016-3215-z
- PubMed
- Google Scholar
(2013) Correlations in ion channel expression emerge from homeostatic tuning rules
PNAS 110:E2645–E2654.

https://doi.org/10.1073/pnas.1309966110
- PubMed
- Google Scholar
(2014) Cell types, network homeostasis, and pathological compensation from a biologically plausible ion channel expression model
Neuron 82:809–821.

https://doi.org/10.1016/j.neuron.2014.04.002
- PubMed
- Google Scholar
1. O'Leary T
2. Marder E
(2016) Temperature-Robust neural function from Activity-Dependent ion channel regulation
Current Biology 26:2935–2941.

https://doi.org/10.1016/j.cub.2016.08.061
- PubMed
- Google Scholar
(2006) Noisy inputs and the induction of on-off switching behavior in a neuronal pacemaker
Journal of Neurophysiology 96:3338–3348.

https://doi.org/10.1152/jn.00486.2006
- PubMed
- Google Scholar
Book
(1988)
Numerical Recipes in C

Cambridge University Press.
- Google Scholar
(2004) Similar network activity from disparate circuit parameters
Nature Neuroscience 7:1345–1352.

https://doi.org/10.1038/nn1352
- PubMed
- Google Scholar
1. Ranjan R
2. Logette E
3. Marani M
4. Herzog M
5. Tâche V
6. Scantamburlo E
7. Buchillier V
8. Markram H
(2019) A kinetic map of the homomeric Voltage-Gated potassium channel (Kv) Family
Frontiers in Cellular Neuroscience 13:358.

https://doi.org/10.3389/fncel.2019.00358
- PubMed
- Google Scholar
Preprint
(2018) Neural circuit robustness to acute, global physiological perturbations
bioRxiv.

https://doi.org/10.1101/480830
- Google Scholar
(2013) The effects of temperature on the stability of a neuronal oscillator
PLOS Computational Biology 9:e1002857.

https://doi.org/10.1371/journal.pcbi.1002857
- Google Scholar
1. Robertson RM
2. Money TG
(2012) Temperature and neuronal circuit function: compensation, tuning and tolerance
Current Opinion in Neurobiology 22:724–734.

https://doi.org/10.1016/j.conb.2012.01.008
- PubMed
- Google Scholar
(2014) Cell-intrinsic mechanisms of temperature compensation in a grasshopper sensory receptor neuron
eLife 3:e02078.

https://doi.org/10.7554/eLife.02078
- PubMed
- Google Scholar
1. Rosenthal JJ
(2015) The emerging role of RNA editing in plasticity
Journal of Experimental Biology 218:1812–1821.

https://doi.org/10.1242/jeb.119065
- PubMed
- Google Scholar
1. Ruff RL
(1999) Effects of temperature on slow and fast inactivation of rat skeletal muscle na ⁺ channels
American Journal of Physiology-Cell Physiology 277:C937–C947.

https://doi.org/10.1152/ajpcell.1999.277.5.C937
- Google Scholar
1. Santin JM
2. Schulz DJ
(2019) Membrane voltage is a direct feedback signal that influences correlated ion channel expression in neurons
Current Biology 29:1683–1688.

https://doi.org/10.1016/j.cub.2019.04.008
- PubMed
- Google Scholar
1. Schauf CL
(1973) Temperature dependence of the ionic current kinetics of Myxicola giant axons
The Journal of Physiology 235:197–205.

https://doi.org/10.1113/jphysiol.1973.sp010384
- PubMed
- Google Scholar
(2006) Variable channel expression in identified single and electrically coupled neurons in different animals
Nature Neuroscience 9:356–362.

https://doi.org/10.1038/nn1639
- PubMed
- Google Scholar
(2007) Quantitative expression profiling of identified neurons reveals cell-specific constraints on highly variable levels of gene expression
PNAS 104:13187–13191.

https://doi.org/10.1073/pnas.0705827104
- PubMed
- Google Scholar
(1999) Heat-shock protein protection
Trends in Neurosciences 22:97–99.

https://doi.org/10.1016/S0166-2236(98)01392-7
- PubMed
- Google Scholar
1. Sjodin RA
2. Mullins LJ
(1958) Oscillatory behavior of the squid axon membrane potential
The Journal of General Physiology 42:39–47.

https://doi.org/10.1085/jgp.42.1.39
- PubMed
- Google Scholar
1. Soofi W
2. Goeritz ML
3. Kispersky TJ
4. Prinz AA
5. Marder E
6. Stein W
(2014) Phase maintenance in a rhythmic motor pattern during temperature changes in vivo
Journal of Neurophysiology 111:2603–2613.

https://doi.org/10.1152/jn.00906.2013
- PubMed
- Google Scholar
(2005) Computational model of electrically coupled, intrinsically distinct pacemaker neurons
Journal of Neurophysiology 94:590–604.

https://doi.org/10.1152/jn.00013.2005
- PubMed
- Google Scholar
(1991)
Distribution and abundance of four brachyuran crabs on the northwest Atlantic shelf

Fishery Bulletin 89:473–492.
- Google Scholar
1. Tang LS
2. Goeritz ML
3. Caplan JS
4. Taylor AL
5. Fisek M
6. Marder E
(2010) Precise temperature compensation of phase in a rhythmic motor pattern
PLOS Biology 8:e1000469.

https://doi.org/10.1371/journal.pbio.1000469
- PubMed
- Google Scholar
1. Tang LS
2. Taylor AL
3. Rinberg A
4. Marder E
(2012) Robustness of a rhythmic circuit to short- and long-term temperature changes
Journal of Neuroscience 32:10075–10085.

https://doi.org/10.1523/JNEUROSCI.1443-12.2012
- PubMed
- Google Scholar
1. Taylor BJR
2. Kerkut GA
(1958) The effect of temperature changes on the activity of poikilotherms
Behaviour 13:259–279.

https://doi.org/10.1163/156853957X00047
- Google Scholar
(2014) Activity-dependent feedback regulates correlated ion channel mRNA levels in single identified motor neurons
Current Biology 24:1899–1904.

https://doi.org/10.1016/j.cub.2014.06.067
- PubMed
- Google Scholar
1. Tran T
2. Unal CT
3. Severin D
4. Zaborszky L
5. Rotstein HG
6. Kirkwood A
7. Golowasch J
(2019) Ionic current correlations are ubiquitous across phyla
Scientific Reports 9:1687.

https://doi.org/10.1038/s41598-018-38405-6
- Google Scholar
(1994) Activity-dependent changes in the intrinsic properties of cultured neurons
Science 264:974–977.

https://doi.org/10.1126/science.8178157
- PubMed
- Google Scholar
(1995) Selective regulation of current densities underlies spontaneous changes in the activity of cultured neurons
The Journal of Neuroscience 15:3640–3652.

https://doi.org/10.1523/JNEUROSCI.15-05-03640.1995
- PubMed
- Google Scholar

Article and author information

Author details

Leandro M Alonso

Volen Center and Biology Department, Brandeis University, Waltham, United States

Contribution
Conceptualization, Software, Formal analysis, Investigation, Visualization, Writing - original draft, Writing - review and editing

For correspondence
lalonso@brandeis.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-8211-2855
Eve Marder

Volen Center and Biology Department, Brandeis University, Waltham, United States

Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Writing - original draft, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-9632-5448

Funding

National Institutes of Health (T32 NS07292)

Leandro M Alonso

National Institutes of Health (MH046742)

Eve Marder

National Institutes of Health (R35 NS097343)

Eve Marder

Swartz Foundation (Swartz Foundation 2017)

Leandro M Alonso

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

LMA acknowledges useful conversations with Hillary Rodgers.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.