Flickering body temperature anticipates criticality in hibernation dynamics

Hibernation has been selected for increasing survival in harsh climatic environments. Seasonal variability in temperature may push the body temperatures of hibernating animals across boundaries of alternative states between euthermic temperature and torpor temperature, typical of either hibernation or summer dormancy. Nowadays, wearable electronics present a promising avenue to assess the occurrence of criticality in physiological systems, such as body temperature fluctuating between attractors of activity and hibernation. For this purpose, we deployed temperature loggers on two hibernating edible dormice for an entire year and under Mediterranean climate conditions. Highly stochastic body temperatures with sudden switches over time allowed us to assess the reliability of statistical leading indicators to anticipate tipping points when approaching a critical transition. Hibernation dynamics showed flickering, a phenomenon occurring when a system rapidly moves back and forth between two alternative attractors preceding the upcoming major regime shift. Flickering of body temperature increased when the system approached bifurcations, which were also anticipated by several metric- and model-based statistical indicators. Nevertheless, some indicators did not show a pattern in their response, which suggests that their performance varies with the dynamics of the biological system studied. Gradual changes in air temperature drove transient between states of hibernation and activity, and also drove hysteresis. For hibernating animals, hysteresis may increase resilience when ending hibernation earlier than the optimal time, which may occur in regions where temperatures are sharply rising, especially during winter. Temporal changes in early indicators of critical transitions in hibernation dynamics may help to understand the effects of climate on evolutionary life histories and the plasticity of hibernating organisms to cope with shortened hibernation due to global warming.

We log-transformed data (using log(T+1)) because of the presence of very fluctuating values between active BT (around 38º) and dormancy BT (close to 0º).
We filtered extreme data by fitting linear AR models, and we also generated partial autocorrelation functions for setting the best choice of time lags for the linear models fitted to the TS Conditional heteroskedasticity (CH) It measures whether variance at one time step is positively associated with variance at one or more previous steps; CH should increase when approaching a critical transition.
We used two lengths of rolling windows (10% and 30% of the total, the latter because of the slow timescale of varying temperatures in dormice). The significance threshold was set to 0.1 Data was standardized because of the presence of extreme values of BT in the dataset.
We filtered by fitting linear AR models Nonparametric Drift-Diffusion-Jump (DDJ) Models yield: conditional variance or drift, diffusion, jumps, and total variance BT was log-transformed (using log(T+1)) Generic early warning signals EWS include: temporal autocorrelation of BT at lag-1, standard deviation SD and skewness We used an overlapping moving window since indicators should be estimated as data become available BT was log-transformed (using log(T+1)) and results were compared with raw BT Data were filtered using Gaussian smoothing 3

Table S4
Coefficients of the parameters of regression AR(p) models of body temperature T, in which air temperature was added (A), compared to AR(p) models not including this covariate. All temperatures were standardized for each of the three periods considered: pre-hibernation, hibernation and post-hibernation (a, b and c respectively, see figure; note that periods a and b correspond to transition to hibernation (TH) and to activity (TA)). Models have the form: , in which p = 1. Standard errors of the mean (SE) are shown in parentheses; AIC = difference in AIC value with the best model in each of the comparisons. All models including air temperature as a covariate showed lower AIC values.

Figure S1
Partial autocorrelation functions for the body temperature time series partitioned for the period before the transition to hibernation (panels a and b for the female and the male, respectively) and before the transition to activity (panels c and d for the female and the male, respectively). Body temperature when dormice were active and prior to hibernation (panels a and b for the female and the male respectively) showed a circadian cyclicity, but not during hibernation.

Figure S2
Results from conditional heteroskedasticity (expressed as the R 2 values of the regressed residuals) for the transitions to hibernation (panels (a) and (b) for female and male respectively) and transitions to activity (panels (c) and (d) for female and male respectively), considering two sizes of rolling windows (lengths 10% and 30% of the total). The red lines show the limit above which R 2 values are statistically significant (P = 0.1). CH was sensitive to rolling window size.

Figure S3
Nonparametric drift-diffusion-jump model-based indicators for body temperature dynamics of the female (panels a to e) and the male (panels f to j) edible dormouse. Left panels show (from top to bottom) how body temperature, conditional variance, diffusion, jump intensity and total variance vary with time, whereas right panels show how the same indicators vary with body temperature. Dashed lines show the two breaking points signalling transition to hibernation and activity. When plotted against body temperature, conditional variance peaked around value 3 (20.1ºC), which was the temperature that separates the two states of hibernation and activity, and jump intensity peaked at value 2.5 (12.2ºC), which measures the average jump in temperature from dormancy to T e and vice versa.  Figure S5 Metric-based rolling window indicators estimated for the body temperature time series for the two dormice for the whole time series (a and b for female and male, respectively). Red lines show the Gaussian filtering of the time series. Panels 1, 2 and 3 show autocorrelation at lag-1 (AR(1)), standard deviation and skewness respectively, estimated within sliding windows of 20% the size of the time series. The Kendall  indicates the strength of the trend in the indicators. Vertical dashed lines show the breaking points separating the two stable states.

Figure S6
Metric-based rolling window indicators estimated for the body temperature time series and air temperature time series (yellow lines) for the two dormice and separated for each alternative state (activity and hibernation). (a and b): activity state and hibernation state for the female; (c and d): activity state and hibernation state for the male. Red lines show the Gaussian filtering of the time series. Panels 1, 2 and 3 show autocorrelation at lag-1 (AR(1)), standard deviation and skewness respectively, estimated within sliding windows of 20% the size of the time series. The Kendall  indicates the strength of the trend in the indicators. Vertical dashed lines show the breaking points separating the two stable states.

Figure S7
Sensitivities of generic EWS (autocorrelation at lag-1 AR(1), standard deviation SD and skewness Sk) for testing the robustness of the rolling window metrics. For each of the two transitions and each dormouse, left panels are the contour plots showing the influence of the width of the rolling window and Gaussian filtering on the observed trend in the metrics as measured by Kendall's ; right panels are the histograms with the frequency distribution of the trend for . Small triangles show the parameter choice used in each of our analyses. Windows sizes ranged from 20% to 80% in increments of 10 points and bandwidths ranged from 5 to 200 in an increment of 20.

Female
Transition to hibernation Male AR (1) AR (1) SD SD Sk Sk