## Abstract

Neurons in cold-blooded animals remarkably maintain their function over a wide range of temperatures, even though the rates of many cellular processes increase twofold, threefold, or many-fold for each 10°C increase in temperature. Moreover, the kinetics of ion channels, maximal conductances, and *Ca*^{2+} buffering each have independent temperature sensitivities, suggesting that the balance of biological parameters can be disturbed by even modest temperature changes. In stomatogastric ganglia of the crab *Cancer borealis*, the duty cycle of the bursting pacemaker kernel is highly robust between 7 and 23°C (Rinberg et al., 2013). We examined how this might be achieved in a detailed conductance-based model in which exponential temperature sensitivities were given by *Q*_{10} parameters. We assessed the temperature robustness of this model across 125,000 random sets of *Q*_{10} parameters. To examine how robustness might be achieved across a variable population of animals, we repeated this analysis across six sets of maximal conductance parameters that produced similar activity at 11°C. Many permissible combinations of maximal conductance and *Q*_{10} parameters were found over broad regions of parameter space and relatively few correlations among *Q*_{10}s were observed across successful parameter sets. A significant portion of *Q*_{10} sets worked for at least 3 of the 6 maximal conductance sets (∼11.1%). Nonetheless, no *Q*_{10} set produced robust function across all six maximal conductance sets, suggesting that maximal conductance parameters critically contribute to temperature robustness. Overall, these results provide insight into principles of temperature robustness in neuronal oscillators.

- conductance-based models
- electrical coupling
*Q*_{10}- stomatogastric ganglion
- temperature
- voltage-dependent currents

## Introduction

Temperature is an influential environmental factor that perturbs the rate of most biochemical processes in a near-exponential fashion. The maximal conductance and activation/inactivation rate of ionic currents increase by a factor (the *Q*_{10}) of ∼2 or more over a 10°C range (Pena et al., 2006; Peloquin et al., 2008; Tang et al., 2010). These *Q*_{10} coefficients can be highly variable. For example, the rate of potassium current activation increases approximately twice as quickly as the rate of sodium activation in *Xenopus laevis* as temperature increases (Frankenhaeuser and Moore, 1963). Temperature therefore perturbs the magnitude and balance of rate constants in all neurons, which has the potential to drastically alter and potentially disrupt circuit function. Indeed, temperature alters the excitability of cultured neurons and acute brain slices (Watson et al., 1997; Masino and Dunwiddie, 1999) and has been used as a tool to manipulate circuit output and animal behavior *in vivo* (Long and Fee, 2008; Fee and Long, 2011; Goldin et al., 2013). Although these results establish the physiological importance of temperature, our mechanistic understanding of these effects is very incomplete.

Many biological systems in diverse organisms are surprisingly robust to large temperature perturbations (Ruoff et al., 2007). Cold-blooded animals and warm-blooded animals that hibernate have evolved neural circuits that tolerate these changes remarkably well (Thomas et al., 1986; von der Ohe et al., 2006; Robertson and Money, 2012). We used computational modeling to examine temperature robustness in the pyloric central pattern generator (CPG) in the stomatogastric ganglion (STG) of the crab, *Cancer borealis*. Both *in vitro* and *in vivo* experiments have shown that the pyloric motor pattern increases in frequency in response to raised temperature, but maintains its phasing, or relative timing (Tang et al., 2010; Goeritz et al., 2013). This phenomenon, in which certain features of rhythmic behavior are invariant across temperatures, is referred to as “temperature compensation” (Ruoff et al., 2007; Robertson and Money, 2012). Temperature compensation of phase is important for maintaining appropriate motor output, and has been observed in other invertebrate motor networks (Barclay et al., 2002; Katz et al., 2004; Armstrong et al., 2006).

We examined the effects of temperature on a detailed conductance-based model of the pyloric pacemaker kernel (Soto-Treviño et al., 2005). This component of the pyloric CPG produces a rhythmic, single-phased output that tightly maintains its duty cycle over an ∼20°C temperature range (Rinberg et al., 2013). We investigated how tightly tuned model parameters had to be to satisfy this functional constraint and assessed model output over broad regions of *Q*_{10} parameter space to determine whether robust models were confined to specific regions of this space or if these parameters were only loosely constrained. We also investigated whether the maximal conductances of intrinsic ionic currents altered the set of permissive *Q*_{10} parameters and the overall temperature robustness of the pacemaker.

## Materials and Methods

##### Experimental model.

The model used here consists of four compartments representing two cells, each with a somatic and axonal compartment. These two model cells are used to represent the three cells of the pyloric pacemaker circuit. The pair of pyloric dilator (PD) cells, which are strongly electrically coupled, are modeled as a single cell. The single anterior burster (AB) cell is modeled as a separate cell. The AB and PD model cells are electrically coupled by a nonrectifying gap junction.

At the reference temperature of 11°C, the model is directly adapted from Soto-Treviño et al. (2005), with two exceptions. First, the value for maximal conductance for *I*_{A} in the AB cell was 21.6 μS, not 200 μS. Second, the activation exponents of *I*_{A} in the two cells should be: AB = *m*^{4} and PD = *m*^{3}. These discrepancies were discovered through observation of model behavior and verified by checking the source code used in the original study.

The membrane potential, *V*, evolves as follows:
where *C* is the compartment's capacitance, *I*_{ionic} is the current due to channels in a compartment's membrane, *I*_{axial} is the current flowing between the somatic and axial compartments of a cell, and *I*_{gap} is current flowing through the gap junction between cells. The voltage-dependent currents included in the model are as follows: sodium (*I*_{Na}), transient calcium (*I*_{CaT}), slowly inactivating calcium (*I*_{CaS}), persistant sodium (*I*_{NaP}), hyperpolarization-induced cation current (*I*_{H}), delayed-rectifier potassium (*I*_{Kd}), calcium-dependent potassium (*I*_{KCa}), A-type potassium (*I*_{A}), and the mixed modulatory inward current (*I*_{MI}).

The current flowing through each conductance is given as follows:
where *ḡ* is the maximal conductance, *m*_{i} is the activation variable, *h*_{i} is the inactivation variable, *p*_{i} is an integer between 1 and 4, *q*_{i} is either 0 or 1, and *E*_{i} is the reversal potential of the conductance.

The activation and inactivation variables approach their steady-state values as follows:
where *m*_{∞} and *h*_{∞} are the steady-state values of the activation and inactivation variables, respectively, and τ_{m} and τ_{h} are their respective voltage-dependent time constants. These functions are listed in Table 1, along with values for *p*_{i} and *q*_{i}.

Calcium buffering is implemented as an exponential decay to steady-state as follows:
where τ_{Ca} is the time constant for *Ca*^{2+} buffering, [*Ca*] is the internal calcium concentration, and [*Ca*]_{min} = 0.5 μm is the minimal intracellular *Ca*^{2+} concentration. *I*_{Ca} is the total *Ca*^{2+} current into the cell through *Ca*^{2+} channels in the membrane and *S* is a constant that converts a current into concentration and is related to the ratio of the surface area of the cell to the volume in which *Ca*^{2+} concentration is considered. Table 2 gives the values for necessary parameters involving the calculation of [*Ca*] and coupling conductances and capacitances for each compartment.

##### Temperature dependence in the model.

The *Q*_{10} temperature coefficient represents the multiplicative change in a parameter over a 10 degree increase in temperature. This can be expressed mathematically as follows:
where ρ_{1} denotes the parameter value at temperature *T*_{1} and ρ_{2} denotes the parameter value at temperature *T*_{2}. We assigned *Q*_{10} values for the gating kinetics of ion channels, maximal conductances, compartmental coupling conductances, and the calcium buffering time constant for the pacemaker model described above. The value of a parameter any given temperature *T* can be calculated given its value at a reference temperature, *T*_{ref}, by rearranging the equation as follows:
where ρ(*T*) is the function describing the temperature dependence of a parameter ρ, and ρ_{ref} is the value of ρ at the reference temperature. For the second equality, we have made the substitution *Q*(*T*) = *Q*_{10}^{(Tref−T)/10}; we refer to *Q* as the temperature scaling factor. Because the model developed by Soto-Treviño et al. (2005) was fit to experimental data at 11°C, we chose *T*_{ref} = 11°*C*.

All conductances (including maximal and compartmental coupling conductances) were directly scaled by the above equation with a *Q*_{10} of 1.6, which is representative of published experimental data on the temperature dependence of maximal conductances (Frankenhaeuser and Moore, 1963; Bennetts et al., 2001). The time constants for the activation and inactivation variables (τ_{m} and τ_{h}) and the time constant of calcium buffering (τ_{Ca}) were modified by dividing by *Q*: τ(*T*) = τ_{ref}/*Q*. For the activation and inactivation variables, this is equivalent to scaling the rate constants for the opening and closing of the gate by the same scaling factor *Q*. The activation and inactivation gating variables follow the standard Hodgkin-Huxley channel model as follows:
where α and β are voltage-dependent functions that represent, respectively, the rates of the forward and backwards chemical reactions. The time constants and steady-state functions for *m* and *h* can be shown to equal (Koch, 1999):
We made the simplifying assumption that the *Q*_{10} for α and β were always equal because this greatly reduces the number of free parameters in the model and is somewhat supported by experimental data (Frankenhaeuser and Moore, 1963; but see, Voets et al., 2004; Tang et al., 2010). Scaling α and β in the above equations by *Q*(*T*) yields the following:
Therefore, given our assumptions, the time constant parameters are divided by the temperature scaling factor *Q*, whereas the steady-state values are insensitive to temperature changes.

Finally, the reversal potentials of all ionic currents were temperature dependent, as specified by the Nernst equation. These changes were relatively small (e.g., *E*_{Na} increases only 1.76 mV over a 10°C range) and did not appear cause substantial effects on model behavior, but were included for realism. To calculate the reversal potential for the mixed cation channels, *I*_{H} and *I*_{MI}, the relative permeabilities of the channel to *K*^{+} and *Na*^{+} were calculated as follows:
where *P* is the relative permeability of the ion, and *E* is the reversal potential of the ion or mixed cation channel. These relative permeabilities, along with the reversal potential of *Na*^{+} and *K*^{+} at the chosen temperature, are used to calculate the reversal potential of these mixed cation channels at that temperature as follows:

##### Model implementation and data analysis.

The above differential equations were numerically integrated using an exponential Euler method (Dayan and Abbott, 2001). This algorithm was implemented in a custom-written C++ program. Each set of model parameters was simulated for 30 s of model time. We characterized each set of parameters based on the final 15 s of simulation, at which time the models displayed their steady-state behavior.

Duty cycle was calculated by examining the activation of the A current. Burst onset was defined as the point when the activation of the A current passed a set threshold of 0.15 and burst offset was defined as the point then activation of the A current fell below 0.10. This definition was chosen because it closely follows the slow-wave oscillation that drives transmitter release in this system and proved more consistent than definition based on spikes per burst. The duty cycle of tonically firing neurons was set equal to one; the duty cycle of silent neurons was set equal to zero. If the interburst interval was not regular (coefficient of variation <0.05), then the simulation was classified as an irregular burster and was excluded from further analysis. This behavior was quite rare (<0.1% of total).

Data analysis was done in the MATLAB computing environment. The source code for the model has been posted to ModelDB (http://senselab.med.yale.edu/ModelDB/ShowModel.asp?model=152636) and may be used with attribution. MATLAB data analysis scripts are available upon request.

##### Mathematical analysis of ‘perfect’ temperature scaling.

In the context of the present study, a “perfectly robust” bursting pacemaker is one that precisely maintains its waveform when normalized to cycle period. By definition, such a model would exactly maintain the same duty cycle at any temperature. Mathematically, this means that the limit cycle is invariant to temperature changes, even though the speed of the system along the limit cycle is multiplied by *Q*. This occurs if the first time derivatives for all state variables (membrane potentials, gating variables, and intracellular calcium), *s*_{i}, are all directly proportional to a common temperature-scaling factor *Q* as follows:
Here, *s*_{i} at the reference temperature (*T* = 11°*C* in our study). The above condition implies that the stability and the shape of the limit cycle are invariant to temperature changes; if a vector field is multiplied by a scalar, all trajectories within that field are preserved (though the system moves quicker or slower along them) and the stability of fixed points and limit cycles are unchanged. This condition is not perfectly satisfied in the model, but approximately holds when all *Q*_{10} parameters are set equal. We are unaware of any simple mathematical analysis to determine whether a system with unequal *Q*_{10} coefficients is temperature compensated. We address this question extensively in the Results section using numerical simulations. If all *Q*_{10} parameters are equal, then the above condition holds exactly for all activation gating variables as follows:
An equivalent argument can be trivially made for all inactivation variables, *h*_{j}, because their dynamics follow the same form as the activation variables.

If we approximate all ionic reversal potentials as being temperature independent, then the membrane potential for each compartment at the reference temperature evolves according to the following:
where *g*_{k,a} represents the junctional conductance between the *k*^{th} and *a*^{th} electrical compartment and *C*_{k} is the capacitance of the compartment. After multiplying all maximal and junctional conductances by *Q*, we get the following:
The above analysis assumes that all reversal potentials, *E*_{j}, do not vary as a function of temperature; however, this is known to be false—the reversal potential of any current is inversely proportional to temperature, as specified by the Nernst equation. Although changes in reversal potentials due to temperature are small (on the order of a few millivolts over a 10°C change), this can potentially have important effects on neuronal activity (Fröhlich et al., 2008). In our model, we did not notice substantial temperature-dependent effects that could be attributed to reversal potential changes (data not shown).

More importantly, the above analysis does not necessarily hold if intracellular calcium or other secondary processes are included in the model. In the pacemaker model analyzed in this paper, the steady-state intracellular calcium, *Ca*_{∞}, increases in proportion to the total calcium current, *I*_{Ca}, as follows:
Here, Ω represents the set of calcium conductances in an electrical compartment. In our model of temperature, all maximal conductances are multiplied by the temperature coefficient *Q* and the time constant of *Ca*^{2+} buffering is divided by *Q*. Implementing these changes we get the following:
Therefore, in our model, the time derivative of [*Ca*] is not simply changed by a scaling factor, and the model's trajectory is therefore temperature dependent. The effects of intracellular calcium are particularly influential in this model because the calcium-dependent potassium current (*I*_{KCa}) plays a crucial role in burst termination. Specifically, the temperature-dependent increase in *Ca*_{∞} may increase the rate of activation of *I*_{KCa} and promote early burst termination. If *I*_{KCa} is activated too quickly, the fast/slow separation of dynamics required for bursting is lost altogether. This analysis is specific to our particular model of intracellular calcium and the details will change depending on the level of biophysical detail that is used (see Discussion).

## Results

### Temperature robustness across a population of pacemaker models

Figure 1*A*, left, shows the connectivity of the pyloric CPG circuit. Spontaneous oscillations are dependent on the pacemaker kernel, which consists of two PD neurons and the AB neuron. The lateral pyloric (LP) neuron and pyloric neurons sequentially burst on rebound from the inhibitory synaptic drive delivered by the pacemaker kernel. Therefore, the duty cycle of the pacemaker plays a critical role in determining the phase relationships of the full motor rhythm. The only synaptic feedback onto the pacemaker kernel from the pyloric rhythm is from the LP neuron and it can be blocked pharmacologically by bath application of picrotoxin (Fig. 1*A*, right; Eisen and Marder, 1982).

Rinberg et al. (2013) demonstrated that the pharmacologically isolated pacemaker kernel precisely maintained its duty cycle as the bath temperature increased despite increasing in frequency (Fig. 1*B,C*). We developed a temperature-sensitive model of this system (Fig. 2*A*) and searched uniformly across *Q*_{10} space for models that produced qualitatively similar temperature stability. In addition to varying *Q*_{10} parameters, we also varied the maximal conductance densities. It is now well established that models with disparate maximal conductance densities can produce similar behaviors (Goldman et al., 2001; Prinz et al., 2003, 2004; Achard and De Schutter, 2006; Sobie, 2009; Marder and Taylor, 2011; Rathour and Narayanan, 2012; Lamb and Calabrese, 2013). Similar parameter variability is observed in experimental measurements of conductance density (Golowasch et al., 1999; Swensen and Bean, 2005; Schulz et al., 2006, 2007; Goaillard et al., 2009; Tobin et al., 2009; Amendola et al., 2012).

We assessed the temperature robustness of the model for six different maximal conductance sets, all of which produced similar activity at the reference temperature of 11°C (Fig. 2*B*, Tables 3, 4). The first maximal conductance set (set #0) was identical to the one developed by Soto-Treviño et al. (2005), which was hand tuned to replicate biological behavior both under standard conditions and in response to current injections. We found the other five maximal conductance sets by uniformly searching over maximal conductance space; we selected five maximal conductance sets with simulated activity that was visually similar to that of set #0. Despite producing nearly identical behaviors (Fig. 2*B*, left, traces), the maximal conductance values were substantially different from each other (Fig. 2*B*, right, bar plots).

We then randomly generated 125,000 distinct *Q*_{10} sets for the time constants of channel gating and calcium buffering; each *Q*_{10} was randomly drawn from a uniform distribution from 1 to 4. A *Q*_{10} coefficent of 1.6 was assigned to all maximal conductances, gap junction conductances, and axial conductances (see Materials and Methods). We simulated each of the six maximal conductance sets with all 125,000 *Q*_{10} sets. The same *Q*_{10} sets were applied to each maximal conductance set to make individual cases directly comparable. This also removed any possible inconsistencies that using different *Q*_{10}s for different maximal conductance sets might introduce. We also simulated maximal conductance set #0 (same as Soto-Treviño et al., 2005) with an additional 125,000 randomly selected *Q*_{10} sets, bringing the total to 825,000 parameter sets (all combinations of *Q*_{10} sets and maximal conductance sets). We simulated the additional 125,000 *Q*_{10} sets for maximal conductance set #0 to ensure that we had sampled *Q*_{10} space densely enough for the remaining maximal conductance sets.

For each combination of *Q*_{10} and maximal conductance parameters, we simulated model output at five different physiological temperatures (7°C, 11°C, 15°C, 19°C, and 23°C). Models with irregular bursting were quite rare (<0.1% of total) and were excluded from further analysis. We observed a range of model behaviors (Fig. 3*A*) that were consistent with preliminary modeling work from a previous study (Tang et al., 2010). Although some individual models maintained their duty cycle robustly over this temperature range (Fig. 3*A*, “pass”), other individual models were fragile to temperature perturbations (Fig. 3*A*, “fail”).

We quantified the variability of pacemaker duty cycle as the total sum of squares difference (*SST*) from the duty cycle at the reference temperature of 11°C as follows:
Lower *SST _{DC}* values correspond to greater robustness (less variability in duty cycle). Figure 3

*B*displays the distribution of this metric over all

*Q*

_{10}sets for each maximal conductance set. Importantly, the chances of finding a

*Q*

_{10}set that produced a robust oscillator varied across maximal conductance sets. The distributions for maximal conductance sets #1 and #3 tended toward lower (more robust)

*SST*

_{DC}values than the distributions of sets #4 and #5.

The black histograms in Figure 3*B* are those individual models that burst over the entire temperature range. Those to the left of the dashed red line maintained approximately constant duty cycle, whereas those to the right of the red line failed to do so. The blue histograms include those models that failed to burst at one temperature, the green models failed to burst at two temperatures, and the red models failed at three temperatures. In this figure, all simulated *Q*_{10} sets are shown for maximal conductance set #0, which is twice as many compared with the other maximal conductance sets.

### Analysis of robust parameter sets

We performed further analysis on all individual models that produced an *SST*_{DC} <0.01 (Fig. 3*B*, models to the left of the dashed red lines in each maximal conductance set). Simulations within this subpopulation were fairly robust to temperature changes and were considered to be “successful” models. Although this threshold was somewhat arbitrarily chosen, changing it to be more or less restrictive did not qualitatively change the results. Likewise, using an alternative measure of oscillator robustness that was based on the number of spikes per burst rather than duty cycle did not substantially alter our results (Caplan, 2013).

To understand which *Q*_{10}s were most influential in conferring temperature robustness to the models, we plotted histograms of each *Q*_{10} parameter value within the subset of successful models (Fig. 4). All *Q*_{10} sets were originally drawn from a uniform distribution; therefore, flat histograms in Figure 4 correspond to *Q*_{10} parameters that were equally likely to produce robust oscillators across the entire sampled interval, whereas peaks (or troughs) of nonflat distributions show the values of *Q*_{10}s that were likely (or, respectively, unlikely) to be found in the successful models. For example, robust oscillations were more likely to be seen if the *Q*_{10} for *m*_{KCa} was small, because the distribution for this parameter is shifted to the left across all six maximal conductance sets.

Interestingly, when viewed on an individual basis, most *Q*_{10} parameters had qualitatively similar effects across many or all maximal conductance sets. The *Q*_{10}s for *m*_{Na}, *h*_{CaT}, *m*_{NaP}, *h*_{NaP}, *m*_{H}, *m*_{A}, *h*_{A}, and *m*_{MI} displayed flat distributions across all maximal conductance sets (shown in white)—that is, they did not noticeably influence the model's ability to maintain a constant duty cycle (looking down the six rows, the results are similar). In yellow are the *Q*_{10}s that were influential in at least one maximal conductance set. The *Q*_{10}s for *m*_{KCa} and the time constant for calcium buffering had typically lower values in robust oscillators for all conductance sets, whereas the *Q*_{10} for *m*_{CaT} was typically lower for conductance sets #1–4, but noninfluential (flat) in set #5. The *Q*_{10}s for *h*_{Na}, *m*_{Kd}, and *m*_{CaS} showed qualitatively opposite effects across some of the maximal conductance sets. For example, the *Q*_{10} distribution for *m*_{Kd} in Figure 4 is shifted toward higher values in all maximal conductance sets but #4, in which it is shifted toward lower values.

We next examined the pairwise interactions between the six particularly influential *Q*_{10} parameters (highlighted in Fig. 4) in the subpopulation of successful models. In Figure 5, the histograms along the diagonal of each panel show the individual distributions for these *Q*_{10}s (same as the histograms in Fig. 4). The blue-to-red heat maps above the diagonal of each panel show the joint distribution of successful models for each pair of *Q*_{10} parameters; these are 2D histograms for each pair of *Q*_{10}s, with red representing areas of relatively high density of robust models and blue representing areas with relatively few successful models. It is important to note that viable *Q*_{10} sets are spread over large regions of parameter space; they are not tightly confined to specific regions and few strong correlations are observed.

The blue-to-gray-to-pink heat maps below the diagonal on each were calculated by subtracting the product of the two individual distributions (histograms on diagonal) from the joint distribution (plots above diagonal). If the two *Q*_{10}s were completely independently distributed, this would produce a value of zero (gray) in each bin. Magenta areas denote positive values, bins that hold more successful models than is predicted by assuming the two individual distributions are independent. Blue areas denote negative values, bins that hold fewer successful models than expected. One can also interpret these plots as a representation of the information gained from considering the joint distribution of two *Q*_{10}s rather than their separate individual distributions; little additional information is gained in uniformly gray plots, whereas more is gained in plots that have substantial pink and blue patches. These plots are useful to visualize correlations that are not easily perceptible in the joint distribution plots.

As in Figure 4, many similar features can be seen in Figure 5 across the six maximal conductance sets. Most differences arise from the fact that maximal conductance set #4 and #5 have fewer successful models than sets #0–3 and therefore have sparser distributions. A particularly striking relationship is present between the *Q*_{10}s for *m*_{KCa} and calcium buffering; a downward sloping pink curve in the plots below the diagonal for all conductance sets. Other patterns are similarly visible across the maximal conductance sets (e.g., the positive correlation between the *Q*_{10}s for *m*_{CaS} and calcium buffering). An exception is the positive correlation between *m*_{CaS} and *m*_{Kd} in set #4, which is not easily seen in any of the other maximal conductance sets.

### No *Q*_{10} sets produced robust oscillators across all maximal conductance sets

To understand whether the *Q*_{10} combinations that work for individual maximal conductance sets are unique to that maximal conductance set or if they provide more general solutions that can work for many maximal conductance sets, we investigated whether there are *Q*_{10} sets that satisfied *SST*_{DC} < 0.01 for multiple maximal conductance sets.

No *Q*_{10} set satisfied this criterion for all maximal conductance sets. When we required that a model work for 5 of the 6 conductance sets, only ∼1.2% of *Q*_{10}s passed this criterion (Fig. 6*A*). Given the dense sampling of *Q*_{10} parameter space, this result suggests that it is difficult (and perhaps biologically implausible) to find a set of *Q*_{10} values that supports temperature compensation for all viable maximal conductance sets. Therefore, tuning maximal conductances appears to be important for designing neuronal systems that are robust to temperature perturbations.

Next, we relaxed the test by requiring that each *Q*_{10} set only pass in 3 of 6 conductance sets (Fig. 6*B*). Approximately 11.1% of the *Q*_{10} sets passed this criterion. These *Q*_{10} sets were distributed similarly to the *Q*_{10} sets that produced temperature compensation for each maximal conductance set (cf. Fig. 5). We then examined *Q*_{10} sets that worked for any one of the six maximal conductance sets (Fig. 6*C*). Only ∼36.2% of *Q*_{10} sets satisfied this relaxed criterion. Overall, these results suggest that many *Q*_{10}sets are incompatible with temperature compensation of duty cycle regardless of underlying maximal conductance profile of the neuron. This suggests that there exists a subset of nonpermissible *Q*_{10} parameters that cannot be compensated for by tuning maximal conductances. Overall, Figure 6 suggests that both maximal conductances and *Q*_{10} parameters must fall within certain ranges to produce temperature robustness across a variable population of models. This suggests that both of these parameters types may be tuned by evolutionary or other processes (see Discussion).

### Do *Q*_{10}s have to be similar to produce robust temperature scaling?

In principle, if all of the *Q*_{10}s in a complex model were identical, one would get close to perfect temperature scaling (see Materials and Methods; Robertson and Money, 2012). This raises the question of whether models with *Q*_{10}s more similar to each other are, on average, more robust than models with dissimilar *Q*_{10}s.

We calculated the coefficient of variation (CV) for the *Q*_{10} parameters to quantify the level of similarity within each *Q*_{10} set. The CV is the corrected sample SD normalized to the sample mean; large values correspond to dissimilar *Q*_{10} sets and small values correspond to sets in which *Q*_{10} coefficients are all similar. Surprisingly, this variable held extremely little predictive power for the robustness of the oscillator as measured by *SST*_{DC}; the coefficient of determination (*R*^{2}) between these two variables was <0.01 for all maximal conductance sets. This was true even if we calculated the CV only across the six “influential” *Q*_{10}s highlighted in Figure 5 and further analyzed in Figure 6 (data not shown), suggesting that this result does not simply arise from variability in the deconstrained *Q*_{10} parameters (i.e., the ones not highlighted in Fig. 5). We also did not find notable relationships between oscillator robustness and alternative measures of *Q*_{10} variability such as the interquartile range, which is a robust measure of variability in the presence of outliers (data not shown).

Therefore, we conclude that one cannot determine the robustness of a *Q*_{10} set simply by inspecting the similarity of the *Q*_{10} coefficients. This result is illustrated in Figure 7, which shows the distribution of the CV for all *Q*_{10} sets in each maximal conductance set. Note that the CV for the robust *Q*_{10} sets (blue histograms; *SST*_{DC} < 0.01) were similarly distributed to the CV of all simulated *Q*_{10} sets (red histograms). If robust models were more likely to have similar *Q*_{10} values, the blue histograms would have been shifted to the left with respect to the red histograms.

## Discussion

Acute temperature changes are global perturbations that simultaneously alter many cellular processes and can therefore disrupt neuronal function (Robertson and Money, 2012; Tang et al., 2012). Animals have evolved to flourish despite significant daily and seasonal changes in temperature using compensatory mechanisms that work well, even if they remain mysterious to us. Cold-blooded animals also regulate their body temperature through behavioral preferences (Lagerspetz and Vainio, 2006; Garrity et al., 2010; Sengupta and Garrity, 2013) and acclimate physiologically to long-term temperature changes over days, weeks, or months (Camacho et al., 2006; Garrett and Rosenthal, 2012b; Tang et al., 2012).

In contrast to long-term adaptation or active tuning, this study examines the intrinsic robustness of neural systems to acute temperature perturbations. We specifically focused on temperature compensation of duty cycle in the pyloric pacemaker kernel of the crustacean STG (Rinberg et al., 2013). A number of other invertebrate circuits also show temperature compensation of phase, including the circuits controlling *Tritonia* swimming (Katz et al., 2004), *Drosophila* larval crawling (Barclay et al., 2002), and locust ventilation (Armstrong et al., 2006).

We examined a population of temperature-sensitive computational models by varying *Q*_{10} and maximal conductance parameters. Our results address two major issues: (1) how a rhythmic dynamical system can maintain its phasing (or duty cycle if single-phase) despite having many components with highly variable temperature sensitivities, and (2) how temperature robustness can be achieved reliably across an intrinsically variable population of neurons (Swensen and Bean, 2005; Schulz et al., 2006, 2007; Goaillard et al., 2009; Tobin et al., 2009; Amendola et al., 2012; Roffman et al., 2012).

An important result is that the same sets of *Q*_{10}s do not uniformly produce temperature compensation across the six maximal conductance sets we examined. No *Q*_{10} set worked for all maximal conductance sets; however, a reasonable fraction of *Q*_{10} sets (∼11%) worked for at least three maximal conductance sets. This strongly suggests that both maximal conductances and *Q*_{10} parameters must be tuned in biological systems that display robust temperature compensation. One possibility is that *Q*_{10}s and conductance densities are coregulated during development and the lifetime of the animal to ensure that they are appropriately matched. Alternatively, *Q*_{10} coefficients may be genetically determined and fixed over an animal's lifetime. In this case, developmental mechanisms must establish a maximal conductance profile that both produces appropriate activity and is well matched to the *Q*_{10} parameters. In turn, the *Q*_{10}s may have been tuned over evolutionary timescales to accommodate variable conductance expression.

Ion channel protein structure determines the *Q*_{10} of channel activation and inactivation (Zheng, 2013). Mutagenic analyses of the canonical heat-sensing channel TRPV1 have demonstrated that temperature sensitivity varies dramatically as protein structure and composition are manipulated (Vlachová et al., 2003; Brauchi et al., 2006; Yao et al., 2011). Therefore, one might imagine that minor allelic fluctuations in ion channel proteins found in the wild-type population produce animal-to-animal variations in the *Q*_{10}s for some channel proteins. This would allow these parameters to be tuned over many generations by natural selection. To our knowledge, the extent of animal-to-animal variability in *Q*_{10} parameters has not been characterized experimentally in detail.

Neurons may also have the ability to actively regulate *Q*_{10} parameters to foster temperature stability. Indeed, crabs acclimated at 19°C are less susceptible to circuit “crashes” at high temperatures than crabs acclimated at 7 or 11°C (Tang et al., 2012). Our results suggest this could be achieved by tuning maximal conductance parameters, perhaps through homeostatic regulatory pathways (Liu et al., 1998; Turrigiano, 2011; Davis, 2013; Williams et al., 2013). It is also plausible that these neurons tune their channel *Q*_{10}s to achieve greater robustness; this could potentially be achieved by altering ion channel composition through posttranslational modifications, RNA editing, or altering the relative expression of ion channel subtypes. If this were the case, then it is also possible that maximal conductances are coregulated with channel *Q*_{10} parameters.

### Other possible mechanisms for biological temperature robustness

The robustness of a neuron to temperature perturbations is partially dependent on its *Q*_{10} parameters. We also show that maximal conductance parameters strongly influence temperature robustness. In addition to tuning maximal conductances, neurons might also manipulate the voltage and time dependence of ionic currents to support robust function. Indeed, this strategy appears to be used by polar octopus species that use RNA editing to express *K*^{+} channel isoforms with accelerated gating kinetics relative to tropical octopus species (Garrett and Rosenthal, 2012a). Intriguingly, RNA-editing mechanisms have been proposed to promote long-term temperature robustness within individual animals (temperature acclimation) in addition to explaining these cross-species differences (Garrett and Rosenthal, 2012b).

Tuning intracellular *Ca*^{2+} regulation is another potential mechanism to improve neural robustness. Intracellular *Ca*^{2+} is controlled by a rich set of temperature-dependent processes, such as sequestration into internal compartments or by chemical chelators (Berridge et al., 2000). We did not attempt to model these details because they are still poorly characterized within crustacean stomatogastric neurons. Nevertheless, given the influential nature of intracellular *Ca*^{2+} in our results, tuning these regulatory pathways is another potential mechanism for achieving temperature robustness.

### Importance of intracellular calcium dynamics to temperature robustness

The *Q*_{10} coefficients for *m*_{KCa} and the time constant for calcium buffering were particularly important for establishing temperature stability in these models. This result likely stems from the fact that both of these parameters strongly influence *I*_{KCa}, which plays a critical role in terminating each burst (Soto-Treviño et al., 2005). Because the steady-state intracellular calcium level (*Ca*_{∞}) increases with temperature (see Materials and Methods), *I*_{KCa} is likely to activate earlier in the burst phase at higher temperatures and disrupt the burst duty cycle. Models with low *Q*_{10} coefficients for *m*_{KCa} and τ_{Ca} were overrepresented in the successful subset of simulated models, suggesting at least two possible compensations for this disruption. These results are consistent with a simple conceptual model, in which several degenerate model parameters contribute to the temperature dependence of *I*_{KCa} and, by extension, duty cycle constancy.

A previous study (Rinberg et al., 2013) examined temperature stability in a Morris-Lecar oscillator as a reduced model of the pyloric pacemaker kernel. This simple model also predicted that disparate maximal conductance sets would produce variable responses to temperature perturbations, even when consistent *Q*_{10} sets are used. Rinberg et al.'s (2013) results suggested that the *Q*_{10} coefficients for the time constants of outward currents should be generally larger than the *Q*_{10}s for inward current gates in robust models; however, we did not observe this. This discrepancy may result from the intracellular *Ca*^{2+} dynamics that are absent in the Morris-Lecar model, but were present in the detailed model analyzed herein.

### General lessons for parameter fitting

Finding appropriate parameter values for neuron and network models is a longstanding problem in computational neuroscience. A successful approach for moderately sized models has been to randomly search parameter space for models that match biological data by some set of physiological measures (Goldman et al., 2001; Prinz et al., 2003, 2004; Achard and De Schutter, 2006; Taylor et al., 2009; Ball et al., 2010; Rathour and Narayanan, 2012; Lamb and Calabrese, 2013). Models that satisfy these constraints are considered “successful.” This population-based approach prevents investigators from studying the idiosyncrasies of individual parameter sets (Marder and Taylor, 2011).

One concern with this approach has been that biological neurons are subject to more constraints than are typically considered by investigators; therefore, these modeling studies may be too lenient and overestimate the number of feasible parameter sets (Nowotny et al., 2007). In this study, we considered effective temperature compensation as an additional biological constraint. This additional constraint suggests that certain maximal conductance sets may not be feasible solutions, even if they appear “successful” at 11°C. However, we still uncovered a large and widely spread population of successful parameter sets (∼12.5% of all simulated sets) by searching over both maximal conductance and *Q*_{10} parameter space.

Incorporating further constraints into a model (e.g., temperature) can introduce more free parameters. Here, we added many *Q*_{10} parameters to make the model temperature dependent and found many feasible models within this expanded parameter space. Although biological systems probably do not explore this entire space of possible solutions (O'Leary et al., 2013), computational modeling can provide insight by exploring this larger space. Biological systems ultimately find solutions in high-dimensional parameter spaces because they must simultaneously satisfy many complex functional constraints (Brezina and Weiss, 1997). In addition to temperature compensation, other constraints could include reliable neuromodulation (Grashow et al., 2009, 2010), energy efficiency (Laughlin and Sejnowski, 2003), and gain control (Abbott et al., 1997). Our results emphasize that there are multiple mechanisms for neurons to achieve temperature compensation (and presumably other constraints). This contrasts with the intuition one might gain from studying reduced models, which typically identify a handful of “critical” variables. This illustrates the utility of comparing results from models with varying degrees of complexity and detail.

## Footnotes

This work was supported by the NIH (Grant NS 081013). We thank Tilman Kispersky and Timothy O'Leary for providing helpful feedback on an early draft of this manuscript.

The authors declare no competing financial interests.

- Correspondence should be addressed to Dr. Eve Marder, Volen Center for Complex Systems, Brandeis University, 415 South Street, Waltham, MA 02454. marder{at}brandeis.edu