An improved method for detecting torpor entrance and arousal in a mammalian hibernator using heart rate data

Hibernation in mammals, typified by the thirteen-lined ground squirrel [Ictidomys tridecemlineatus (Mitchill 1821)], consists of individual torpor bouts with two extreme metabolic and thermal states (see Staples, 2016 for a comprehensive review). Torpor periods are characterized by low and fairly constant core body temperature (T_b; near 5°C at 5°C ambient) and metabolic rate (MR). Torpor is punctuated approximately every 10 days by periods of interbout euthermia (IBE) with T_b near 37°C and fairly stable MR, similar to those of summer, non-hibernating conspecifics. An arousal between torpor and IBE raises MR as much as 100-fold, followed by a rapid increase in T_b. The decline in MR and T_b during entrance into torpor is somewhat slower than their increase before arousal.

The rapid physiological changes that occur during entrance and arousal have inspired research into the mechanisms underlying the induction and cessation of these periods. Although the hunt for a ‘hibernation induction trigger’ has not yet been successful (see Staples, 2016 for a review), changes in several circulating metabolites during these transitions suggest significant reorganization of tissue metabolism (Epperson et al., 2011). Fully understanding the signals that trigger such changes will require rapid tissue sampling at precise time points before, during and after these transitions.

Most experimental designs define torpor as the period during which T_b either: (1) remains below a certain value [typically 30°C (e.g. Zervanos et al., 2014), but sometimes as low as 25°C (e.g. Bieber et al., 2014)], or (2) declines more than a designated amount below ‘euthermic resting values’ (Ruf and Geiser, 2015). Such designs may be adequate for many field and laboratory studies, especially on small animals, but may not be able to detect key signals underlying the thermal and metabolic transitions of entrance and arousal. Indeed, such criteria would miss the latter stage of arousal and the beginning of the entrance stage, when many physiological changes occur. For example, in thirteen-lined ground squirrels, liver mitochondrial succinate metabolism is substantially suppressed between IBE and torpor (Mathers et al., 2017), but this suppression is fully realized long before T_b falls below 30°C during entrance (Chung et al., 2011). Some recent studies have employed more rigorous sampling regimes, designating ‘early arousal’ as occurring at T_b values between 7 and 12.8°C (Hindle and Martin, 2014). Unfortunately, changes in T_b lag behind those in MR (Ortmann and Heldmaier, 2000), so even experiments that use such rigorous T_b criteria may miss signals that initiate entrance and arousal, and their metabolic effects. More accurate methods for defining the transitions among torpor bout stages are required in order to fully understand the regulation of hibernation.

Improving our identification of these transitions will also advance our understanding of evolutionary and environmental patterns of hibernation. The duration of torpor and IBE periods depend largely on time of year and ambient temperature but, within individual species, these periods may differ among populations depending on latitude (Zervanos et al., 2010) and local environmental factors such as rainfall (Lehmer et al., 2006), as well as food quantity (Humphries et al., 2003) and quality (Munro et al., 2005). Reliance on arbitrary T_b thresholds for quantifying these periods has at least two characteristics that limit its utility. Firstly, to account for variability in T_b, the thresholds are quite conservative (e.g. >5°C below euthermic resting values; Ruf and Geiser, 2015). Secondly, changes in T_b lag behind those in MR (Ortmann and Heldmaier, 2000), so T_b will underestimate the duration of both arousal and entrance periods, and overestimate torpor and IBE. The magnitude of these errors will be greater for larger animals with higher thermal inertia and some large hibernators show only modest drops in T_b despite large changes in MR (Tøien et al., 2011). Improving the accuracy of these durations may be particularly important for intraspecific studies on the subtle effects of body mass on the duration of torpor and IBE periods (e.g. Bieber et al., 2014; Zervanos et al., 2014).

Some researchers have measured MR, using flow-through respirometry, to quantify torpor bout stages (e.g. Ortmann and Heldmaier, 2000). Such measurements found that entrance begins at T_b values significantly higher than 30°C in a variety of hibernators weighing 15–406 g (Willis, 2007). While accurate at detecting transitions between stages, this approach has several disadvantages. Firstly, it requires costly systems for both respirometry and T_b telemetry. Secondly, respirometry measurements require sealing animals in air-tight chambers for extended periods, an unnatural situation that may affect torpor patterns (Willis, 2007). Finally, the most affordable respirometry systems can monitor only eight animals at a time, half that of most telemetry systems.

In many endotherms (reviewed in Green, 2011), including marsupial hibernators (Swoap et al., 2017), heart rate (f_H) changes follow MR closely, with little lag, and thus f_H is a good candidate for monitoring hibernation states. Indeed, f_H is a good correlate of MR during torpor in bats (Currie et al., 2014). Electrocardiogram (ECG) telemeters have been used to characterize f_H during torpor in Glis glis (Elvert and Heldmaier, 2005) and a marsupial hibernator (Swoap et al., 2017), both facing a changing thermal regime. To our knowledge, however, no studies have used telemetric ECG data to identify transitions between torpor bout stages in a hibernator, perhaps because f_H can vary greatly within and among individuals.

Mathematical transformations of ECG data have been used for several years in clinical medicine to analyze f_H variability caused by cardiac malfunctions (e.g. John et al., 2017). Furthermore, spectral analysis of interbeat intervals has been used extensively to analyze regulation of cardiac function, for example in ectotherms under conditions that change metabolism using thermal (Seebacher and Franklin, 2004) or hormonal (Little and Seebacher, 2014) treatments. In this paper, we describe methods for analyzing telemetric ECG data from hibernators that accurately determine the beginning of the arousal and entrance periods. Our data show that this technique offers superior temporal resolution than using T_b alone. These data also permitted us to examine potential changes to cardiac regulation during these transitions.

Data were collected from thirteen-lined ground squirrels, one of the most widely used experimental models for laboratory-based hibernation studies. All procedures were approved by the University of Western Ontario Animal Care Committee (protocol 2012-016) and followed Canadian Council on Animal Care guidelines. These animals were also used in our studies of mitochondrial metabolism (Mathers et al., 2017), so no additional animals were required for the current study. Details of ground squirrel trapping and husbandry can be found in our recent publications (e.g. Mathers et al., 2017). In summer, radio telemeters (ETA-F10 Data Sciences International, St Paul, MN, USA) were implanted under isoflurane anesthesia (Cesarovic et al., 2011), but were not switched on until hibernation began in late October. The transmitter body was placed within the abdominal cavity with the two leads stitched to the muscle of the thorax and abdomen, respectively. The surgery is only slightly more complicated than implantation of simple body temperature telemeters, and each procedure took no more than 30 min total. These small transmitters (1.6 g and 1.1 cm³ volume) represented 0.8–1% of animal mass (range 152–208 g). The transmitters were designed for mice (approximately 25 g) so would be suitable even for small hibernators but, since the time of data collection (winter 2014), new transmitter models have been introduced with similar capabilities. Transmitters were calibrated by the manufacturer for temperature (resolution 0.05°C, accuracy 0.25°C).

Telemetry data were collected using Data Sciences International data exchange matrix receivers placed directly under cages, allowing for continuous monitoring within the animals' home cages without further disturbances. This system can monitor up to 16 animals at once. ECG and T_b signals were collected in 1-min intervals every 3 min using the Acquisition application of Dataquest Advance Research Technology software (V4.3, Data Sciences International). From eight different animals, data from one randomly selected torpor bout between early January and early March were analyzed. Five of these animals were male and three female. Masses did not differ significantly between sexes [t-test, P=0.93; males 174.4±17.3 g (mean±s.e.m.), females 176.1±10.2 g]. For one additional animal (male, 151.9 g), we analyzed data from three consecutive torpor bouts spanning 20 January to 11 February.

f_H was calculated from the R–R intervals of ECGs using Ponemah software (Data Sciences International; V5.20-SP2). The R peaks were first identified using the automated ECG Analysis Attributes tool with the ‘Rats’ setting. We used a QRS detection threshold of 40% and a minimum R deflection of 0.25 mV. This automated process did not accurately detect R peaks at f_H values below 80 beats min^–1, so we resorted to manual identification of R peaks. All data were analyzed post hoc, i.e. not in real time.

During arousal and IBE, f_H is highly variable over short time frames, with differences between consecutive sampling points as high as 150 beats min^–1. Such variability makes it impractical to use absolute or relative thresholds of f_H to define entrance into torpor. We chose, therefore, to digitally filter the f_H data. We then used discrete differential analysis to define relative thresholds as criteria to determine when an animal was entering into a torpor bout.

To determine the transition between IBE and torpor, raw f_H data were transformed using a forward–backward Butterworth low-pass filter with a parametric critical frequency of 0.03 half cycles per f_H data point. Much more than a simple moving average, this technique strips high-frequency content from the signal above the parametric critical frequency, thereby minimizing sensitivity to volatility over short time frames. For most of the periods under study, the transformed f_H data (f_H-LP) tracks raw f_H smoothly (Fig. 1), allowing for small-window discrete differential analysis. Low-pass and related frequency domain filters have been used for f_H analysis since the early 1970s (Saykrs, 1973). Our algorithm can be adapted to different sampling protocols by changing the parametric critical frequency in inverse proportion to the sampling rate, i.e. if the f_H is determined every 1.5 min (instead of 3 min) the critical frequency must be divided by a factor of 2. Because the filter removes only high frequency content, a filtered signal from animals with lower f_H variability, such as hibernating bears (Laske et al., 2010), would resemble the raw signal more closely than it would from animals with more variable f_H, such as our ground squirrels. It is important, however, that users verify the efficacy of the filtering and threshold criteria for different animal models. All analyses were performed using custom software comprised of signal processing functions available in SciPy, an open-source library for scientific computing, and standard Python packages. The software information can be found in the Appendix. The use of forward–backward, first-order Butterworth filtering together with small-window discrete differential analysis is computationally inexpensive and suitable for real-time implementation.

From the f_H-LP data, we determined entrance into torpor using two criteria. Firstly, f_H-LP had to fall below a threshold proportion of the maximum values. For each data set, we compared three such thresholds: 75, 70 and 65% of the maximum f_H-LP value observed during IBE. Secondly, when these thresholds were reached, f_H-LP had to decrease for three consecutive data collection points. For each data set, the accuracy of the determinations was confirmed using T_b. For example, if the f_H-LP criteria determined that the animal was entering torpor, but the T_b showed that the animal subsequently continued in IBE, this result was deemed a false positive.

We also compiled T_b data from each animal as they reached our f_H-LP criteria for entrance into torpor, a similar approach as used by Willis for MR (Willis, 2007). We then explored the possibility that these T_b values could serve as a proxy to detect entrance. We calculated the minimum T_b and mean T_b–95% confidence interval for both the 70 and 65% f_H-LP criteria. We then analyzed T_b data from individual bouts of 30 different I. tridecemlineatus hibernating under the same conditions and over the same time periods. We determined in how many of these bouts animals entered torpor directly if T_b fell below these values. False positive predictions were identified as bouts during which T_b fell below these levels but subsequently rose above them for longer than 30 min, indicating that animal was still in IBE.

To detect the beginning of arousal, we used a similar HR threshold technique but, because HR during torpor is remarkably invariable within individuals, we were able to use raw f_H data without filtering. We set criteria for the beginning of arousal as the time at which f_H increased beyond 5 beats min^–1 and had increased for three consecutive data collection points.

We also assessed the potential for differential autonomic control of cardiac function between entrance and arousal using the Data Parser tool of the Ponemah software to assess HR variability in the raw data. This analysis converts R–R intervals from successive heart beats to frequency domains. Higher frequency domains (e.g. 0.193–0.700 Hz) indicate parasympathetic regulation, whereas sympathetic (adrenergic) regulation is indicated by lower frequencies (0.022–0.07 Hz) (Altimiras, 1999). The data were compared by two-way ANOVA.

Examples of ECG signals collected from one animal in different stages of the same torpor bout are presented in Fig. S1. Although somewhat noisy in arousal, a clear R signal was always visible. Around the midpoint of a torpor bout (3–5 days after entrance), f_H was remarkably stable at 2.2 beats min^–1 (Table 1). During IBE, f_H was quite variable but, following our data filtering algorithm, f_H-LP tracked raw f_H closely, at least at high f_H values (Fig. 1). At lower f_H values, e.g. early during arousal and late during entrance, the filtering consistently overestimated f_H, but our analysis did not consider f_H-LP data for these periods. In all stages, the low-pass filter reduced data variability substantially (Table 1). Around the IBE midpoint, f_H was approximately 140-fold higher than during steady-state torpor. For eight individuals analyzed, maximum f_H occurred late in arousal when T_b had risen to approximately 33°C. These patterns agree closely with those from the same individual over three consecutive bouts (Fig. 1).

Our criteria detected the beginning of arousal accurately in all cases: after f_H reached 5 beats min^–1, with increases over three consecutive sampling points, all animals continued to arouse into IBE. At the time points where arousal criteria were met, T_b was statistically indistinguishable from torpor (P=0.195; Fig. 2, Table 1). An average of 51.9±3.1 min elapsed between the time at which our f_H criteria indicated arousal and the time it took T_b to increase from 5 to 7°C, an indicator of ‘early arousal’ previously used in this species (Hindle and Martin, 2014). This time lag did not depend on animal mass (P>0.05; Fig. S2). By the time T_b reached 7°C, f_H had increased more than 100-fold (P<0.001, paired t-test; Table 1). We suggest that studies that use T_b to detect arousal, even with stringent criteria, risk missing signals that underlie substantial increases in HR and, presumably, many other physiological parameters. Although the HR criteria detected arousal almost an hour before T_b, this temporal difference represents <1% of a typical torpor bout duration in this species (Russell et al., 2010), although in larger hibernators this proportion will likely be higher.

The data-filtering algorithm decreased variability, allowing for discrete differential analysis of f_H-LP data during entrance. Applying two of our criteria (f_H-LP falling to 70 or 65% maximal f_H-LP, and decreasing over three consecutive sampling points) predicted entrance into torpor bouts with no false positives in the 11 data sets analyzed. We rejected the 75% maximal f_H-LP criteria because it produced one false positive: after the criteria were met, the animal did not proceed directly into torpor, and both T_b and f_H-LP increased temporarily. We conclude that this technique can accurately detect entrance if appropriate criteria are employed.

Even using conservative criteria, our method provides considerably better temporal resolution for detecting entrance than previous methods. When our 70 and 65% f_H-LP were reached in eight individual ground squirrels, T_b was significantly lower (P<0.001, repeated measures ANOVA) than T_b in mid-IBE, but by only 1.3–2.0°C (Table 1). These T_b values were close to those determined in one individual over three consecutive bouts (35.7±0.4°C and 35.3±0.2°C, respectively) but considerably higher than the 30°C used previously to identify entrance in this species (Chung et al., 2011) and other hibernators (Zervanos et al., 2014). These T_b values agree very closely with MR data used to identify entrance in various hibernators, including ground squirrels of a similar size and held at similar temperatures as the animals in the present study (Willis, 2007), demonstrating that f_H serves as an excellent proxy for MR in these systems. By the time T_b fell to 30°C, f_H fell by close to 60% compared with values at mid-IBE. The average time that elapsed between f_H-LP reaching the 65% criteria and T_b falling to 30°C was 54.9±10.7 min. Using 70% criteria significantly increased (one-tailed paired t-test, P=0.000178) this elapsed time to an average of 76.0±14.9 min. There was no significant effect of sex on these lag times (t-test, P>0.05), nor any significant correlation (P>0.05; Fig. S2) with animal mass, but the mass range was quite small (151.9–208.3 g). Larger hibernating species, such as marmots (mass ∼4000 g) have greater thermal inertia, so these lag times will presumably be greater. Our method will improve detection of both entrance and arousal, and provide more accurate determinations of torpor and IBE durations, especially in large hibernators. Indeed, we would be interested to see this method applied to the largest hibernators, such as black bears, from which detailed data sets on ECG, T_b and even MR exist (Tøien et al., 2011).

We recognize that it is more challenging to collect and analyze f_H data than simple T_b data, so we explored the possibility that the T_b at which our f_H-LP occurred could accurately detect entrance into torpor. The minimum T_b at which the 70% criterion was reached was 33.4°C and the mean T_b–95% confidence interval was 34.3°C. For the 65% criterion these values were 32.1 and 33.3°C, respectively. Using data from 30 separate animals we found that these T_b values were poor predictors of entrance. A T_b of 34.3°C accurately predicted entrance in only 36.7% of bouts. Reducing this threshold to the more conservative level of 33.3°C still accurately predicted entrance in only 60% of bouts. Even the most conservative T_b derived from our f_H-LP analysis, 32.1°C, produced false positive predictions of entrance in 20% of the bouts. The T_b threshold of 30°C, frequently used in hibernation studies, accurately predicted entrance in all but one of 30 bouts (96.7% accuracy). Although T_b may be easier to measure than HR, it is still an unreliable predictor of entrance into torpor unless very conservative criteria are applied, by which time the average animal is already approximately an hour into entrance.

Characteristics of the ECGs (e.g. shape, magnitude, phase duration) appear to differ among torpor bout stages (Fig. S1), suggesting differential regulation, as reported recently in a marsupial hibernator (Swoap et al., 2017). Our data suggest that both sympathetic and parasympathetic regulation is significantly lower in the early stages of arousal than during entrance (Fig. S3). During arousal there appears to be a trend towards greater sympathetic than parasympathetic regulation, but the data are quite variable, precluding statistical significance.

Our technique uses fairly simple data filtering that is readily available and widely used in computer sciences and clinical cardiology. We believe that it could be easily adapted to current data acquisition software to provide real-time detection of entrance and arousal. It can also be applied post hoc to improve the accuracy of existing data sets. This method significantly advances the time at which the beginning to entrance and arousal can be identified, so it will improve the ability to detect signals that underlie metabolic suppression and other physiological changes associated with torpor. Moreover, our method will improve quantification of the duration of torpor bout stages and, therefore, our understanding of how they are regulated within and among hibernating species.

The technique is currently best suited to laboratory experiments, and applying it to field studies may be challenging. The telemetry receivers used here must be very close to the instrumented animals and are not designed for outdoor deployment. Some systems (e.g. TSE Systems, Chesterfield, MO, USA) have implants that can transmit up to 5 m, potentially allowing use in semi-natural outdoor enclosures (e.g. Bieber et al., 2014), but, to our knowledge, these transmitters have not been tested at temperatures below 15°C and it is unclear whether the receivers are weatherproof. ECG signals can produce very large data files, a challenge for on-board data logging. Recently, small data storage tags have been used to monitor f_H and T_b in free-moving crustaceans (McGaw et al., 2018), but their storage capacity is limited. Custom-built data loggers have been used to record ECGs in free-ranging black bears (Laske et al., 2010), but their size precludes use in all but the largest hibernators. Hopefully, advances in transmitter and storage technology will allow for broader use of this technique in the future.

Signal filtering and f_H-based threshold detection was performed using a custom code written in Python 3.6 with NumPy for basic mathematical functions and SciPy for signal processing functions. NumPy and SciPy are both free, open-source software libraries included in most distributions of Python for scientific applications.

To use the software, we recommend the Anaconda distribution of Python 3.6, which includes all necessary packages. Anaconda can be obtained freely (https://www.anaconda.com/download/) and is compatible with all major operating systems.

The software assumes that there is a directory called Data/ in the same directory as Analysis.py. Data/ should contain separate .csv (comma separated values) input files for each data set being analyzed, e.g. one torpor–IBE event. The columns of each .csv file should be arranged in the following order: absolute time (any format), T_b (°C) and f_H (beats min^–1). In our data, the sample interval for both T_b and f_H was 4 min. This parameter must be changed on line 204 of Analysis.py if a different sample interval is used. Other parameters, such as the critical frequency for low-pass filtering or f_H-LP values used to define thresholds, are modifiable according to comments in the code.

For each .csv input file processed, two output .csv files are created in a subdirectory Analysis/. The first contains both the original and filtered time series used to compute the various thresholds discussed in this paper. This file is named Data-X.csv, where X is the prefix of the corresponding input file. In this file the columns are arranged as follows: time point index (starting at 0), time (same format as input), HR (beats min^–1), f_H-LP (beats min^–1) and T_b (°C). The second output .csv file, named Analysis-X.csv, contains a summary of T_b (°C), f_H (beats min^–1), f_H-LP (beats min^–1) and absolute time (same format as input) as computed according to various thresholds or conditions.

The Python code, further information and an example of input and output are available at https://github.com/ethancjackson/HR-Analysis.

We thank Alvin Iverson and the staff of the Carman and Area Research Centre (University of Manitoba) for their help in obtaining animals. Sarah McFarlane and Dr Jason Brown helped with surgeries.