Statistical Estimation of Surface Heat Control and Exchange in Endotherms

Barbara Henning^{1}^{*},
Benilton de Sá Carvalho^{2},
José Luiz Boldrini^{3},
Sérgio Furtado dos Reis^{4},
Denis Otávio Vieira de Andrade^{5}

Show more

1. Introduction

Endotherms (mammals and birds) regulate body temperature primarily by balancing metabolic heat production with heat exchanged with the ambient [1] . Heat is exchanged at the interface between body surface and the ambient [1] , and under standard experimental settings, the differential between surface

temperature ( ${T}_{s}$ ) and ambient temperature ( ${T}_{a}$ ), $\Delta T={T}_{s}-{T}_{a}$ , varies linearly or nonlinearly with ambient temperature [2] . Physiologically, linearity or nonlinearity between $\Delta T$ and ambient temperature derive from the absence or presence of cutaneous vasodilation and vasoconstriction, respectively [3] [4] . At the interface between body surface and the ambient, cutaneous vasodilation and vasoconstriction is a major physiological mechanism of vasomotor control, modulating heat exchange with the ambient [1] . Whereas there is experimental evidence for an association between nonlinearity or linearity and control of heat exchange or lack thereof [2] , a theoretical justification is lacking. Recently, Boldrini et al. [5] demonstrated mathematically, based on physical first principles, that linearity or nonlinerity is expected under realistic assumptions of the thermoregulatory process and relevant thermal feedback control mechanisms. Nonlinearity or linearity thus appears to be signatures of vasomotor control. Furthermore, it is well-established that the body surface of an endotherm is not homogeneous with respect to vasomotor control [1] . Detection of signatures of cutaneous vasoreaction across the body surface is fundamental not only for understanding the dynamics of heat control, but also to quantify the relative importance of different body regions for global heat exchange in endotherms [6] [7] [8] [9] . Body surface temperature is accurately measured using infrared thermal imaging (IRT) either at the large scale of whole organisms [10] or at local, small scales such as that of the cornea [11] . Infrared thermal measurements taken at sensible gradients of ambient temperatures allow to estimate whether body surface temperature varies linearly or nonlinearly with ambient temperature, and, therefore, to determine the existence of heat exchange control [6] . The same thermographic data can be used in connection with biophysical models to estimate how different body regions contribute to total heat balance [6] . Such estimates are instrumental to assess the dynamics of heat control and exchange in endotherms [6] . Infrared thermal data are nevertheless intrinsically prone to within-individual correlations, missing values, time-varying covariates, and irregularly-timed measurements. Additionally, there is a need to discriminate between linear or nonlinear functional relationships between surface body temperature and ambient temperature. Therefore, to infer physiological process from patterns in thermographic data, we need reliable statistical methods. Here, we advance a three-pronged framework to test for signatures of linear or nonlinear relationships between body surface temperature and ambient temperature, and to quantify the relative importance of different surface body regions for the maintenance of total heat balance. First, we use polynomial functions that capture linearity and nonlinearity of the functional relationship between surface body temperature and ambient temperature. Second, we use Generalized Estimating Equations [12] , a modeling formalism particularly suited to handle the idiosyncrasies of thermographic data, to estimate the relevant parameters. Third, we assess model parsimony, that is, the compromise between model complexity and residual variance, using a bootstrap-based strategy [13] . We demonstrate this approach using as a reference system the gracile mouse opossum, a Neotropical small mammal that has become a model for research in population ecology [14] [15] [16] [17] , physiology [18] , and bioenergetics [19] .

2. Methodology

2.1. Animals: Field Collection and Care

Specimens of Gracilinanus microtarsus were live-trapped in savanna-like habitat at the Reserva Biológica de Mogi-Guaçu (RBMG) located in the district of Martinho Prado, Mogi-Guaçu, São Paulo (22˚15S; 47˚08W). Field work was carried out from October 2011 to April 2012. Trapping was done on three consecutive nights. A single 8 ´ 8 trapping grid with 64 trapping stations located 10 m from each other was used to capture individual G. microtarsus. A single Sherman live trap (7.5 cm ´ 9.0 cm ´ 23.5 cm; H. B. Sherman Traps, Inc., Tallahassee, Florida) was set on trees at each trapping station about 1.75 m above ground and baited with banana and peanut butter.

Individuals of G. microtarsus were housed in individual cages in an animal room maintained at approximately 23˚C with a 12h/12h light/dark cycle. Gracile mouse opossums were provided with ad libitum water and an amount of food (dry cat and dog food and mango) designed to keep weight gain (Figure 1) similar to that observed under natural conditions, following data from Martins et al. [15] . Trapping and handling methods followed the guidelines of the American Society of Mammalogists (Animal Care and Use Committee 1998).

2.2. Experimental Procedures

Changes in radiative heat exchange at the surface of different body regions of G. microtarsus were monitored by video thermography as ambient temperature was varied from 8˚C to 38˚C. On the day of an experiment, each individual animal was transferred from its home-cage to an experimental container (7.5 cm ´ 9.0 cm ´ 23.5 cm). This container was framed with wire mesh (1 ´ 10 cm), which allowed the monitoring of body temperature surface without any considerable interference. The animal container was then placed inside a temperature-controlled chamber (FANEM Ltd., São Paulo) and the thermal camera (FLIR SC640, FLIR Systems, Inc.) was positioned under the animal cage, allowing us to track body surface temperature while manipulating ambient temperature. First, the individual was allowed one hour to habituate to the experimental conditions at 23˚C. We set the temperature change protocol for a 1˚C stepwise increment or decrement from the initial temperature (23˚C) up to 38˚C or down to 8˚C, respectively. This range of ${T}_{a}$ covers the range of ${T}_{a}$ the gracile mouse opossum is commonly exposed in its natural habitat (Figure 2). Upon each step change in ${T}_{a}$ , animals were kept for 12 min at the target temperature before the next step change. Previous tests indicated that this time interval assured that animals had reached a steady-state with ${T}_{a}$ . Incrementing and decrementing temperature protocols were separated by an interval of at least

Figure 1. Gracile mouse opossum growth curves. Solid line shows growth curve in captivity. Dashed line shows observed growth under natural conditions according to data from Martins et al. [14] . Grey areas represent confidence intervals.

Figure 2. Average annual air temperature variation at Reserva Biológica de Mogi-Guaçu between 1998 and 2010. Red line represents the average air temperature across months. Green line represents the average maximum air temperature across months. Blue line represents the average minimum air temperature across months. Grey areas represent confidence intervals.

4 days, in which animals were returned to their home-cages and allowed to feed and drink. Whether individual animals were first submitted to the incrementing and decrementing temperature protocols was a random decision.

Thermal images were recorded at 10 frames∙s^{−1} using an IR camera streamed to a computer where data acquisition was managed by the software Thermacam 2.9 (FLIR Systems, Inc). Thermal imaging cameras measure the amount of near infrared radiation (typically wavelengths: 8 - 12 nm) emitted by a surface and then convert this measured radiation to a radiative temperature reading according to the Stefan-Boltzmann equation. This equation states that the energy emitted, R (W∙m^{−2}), is proportional to the fourth power of its absolute temperature, T, in Kelvin degrees,

$R=\u03f5\sigma {T}^{4}$ (1)

where,
$\u03f5$ is the emissivity of the surface and
$\sigma $ is the Stefan-Boltzmann constant (5.67 ´ 10^{-8} W∙m^{−2}∙K^{−4}). We assumed that the body surface of G. microtarsus has an homogenous emissivity of 0.96, which is typical of organic materials [8] . Surface temperatures were visualised as images in grey or colour scales [6] .

We analyzed surface temperature ( ${T}_{s}$ ) separately for different body regions including ears (split into distal and proximal regions), feet, tail, chest, back (evaluated when the animal turned on its back) and ventral areas. Initially, we chose one frame for each ${T}_{a}$ from the last two minutes of the exposure period to any given ${T}_{a}$ . Thus, at the time used for data collection, animals had already been exposed to that particular ${T}_{a}$ for at least 10 min. For each chosen frame, the regions of interest were digitally drawn to obtain the average surface temperature of each body part. Typically, we used the same frame to analyze the ${T}_{s}$ from all body parts. In cases in which the animal positioning did not allow for the analysis of all the body regions on the same frame, we used the temporally closest frame to analyze the missing body region.

2.3. Statistical Analysis for the Relationship between T_{s} and T_{a}

2.3.1. Defining the Model Family

The relationship between body surface and ambient temperature can be one of two possible types depending on whether there is no vasomotor adjustment or whether there is vasomotor adjustment [2] [9] . The first type of relationship is for body regions with no capacity of vasomotor adjustment. Because such regions have a constant blood flow that transports a constant amount of heat to the surface, surface temperature varies linearly with ambient temperature. The second type of relationship is for body regions with capacity of vasomotor adjustment. In this case, because blood flow varies nonlinearly with ambient temperature, surface temperature varies nonlinearly with ambient temperature. A convenient function to describe the functional relationship between body surface temperature ( ${T}_{s}$ ) and ambient temperature ( ${T}_{a}$ ) can be conceived as

follows. If the rate of change is constant, we must have $\frac{\text{d}{t}_{s}}{\text{d}{t}_{a}}=c$ . Integrating both sides, we have $\int}\frac{\text{d}{T}_{s}}{\text{d}{T}_{a}}\text{d}{t}_{a}={\displaystyle \int}\text{\hspace{0.05em}}c\text{d}{T}_{a$ , then $f\left({T}_{s}\right)=c{T}_{a}+d$ . Therefore, the function that describes the relationship between ${T}_{s}$ and ${T}_{a}$ is linear. If the rate of change is not constant, we must have $\frac{{\text{d}}^{2}{T}_{s}}{\text{d}{T}_{a}^{2}}=c$ . Integrating both sides we have $\int}\frac{{\text{d}}^{2}{T}_{s}}{\text{d}{t}_{a}^{2}}\text{d}{T}_{a}={\displaystyle \int}\text{\hspace{0.05em}}c\text{d}{T}_{a$ , yielding $\frac{\text{d}{T}_{s}}{\text{d}{T}_{a}}=c{T}_{a}+d$ . Integrating both sides again we have $\int}\frac{\text{d}{T}_{s}}{\text{d}{T}_{a}}\text{d}{T}_{a}={\displaystyle \int}\left(c{T}_{a}+d\right)\text{d}{T}_{a$ , yielding $f\left({T}_{s}\right)=c{T}_{a}^{2}+d{T}_{a}+e$ . Therefore, the

function describing the relationship between body surface temperature $\left({T}_{s}\right)$ and ambient temperature $\left({T}_{a}\right)$ is nonlinear, quadratic in ${T}_{a}$ . Continued differentiation and integration leads to higher order polynomials in ${T}_{a}$ . Polynomial functions thus capture the linearity and nonlinearity intrinsic to the relationship between surface temperature and ambient temperature associated with the mechanism of vasomotor adjustment. Therefore, we modeled the relationship between $\left({T}_{s}\right)$ and $\left({T}_{a}\right)$ with a family of polynomial functions.

2.3.2. Model Estimation

A response variable that is repeatedly measured on the same subject at different time points is the key feature of longitudinal data sets. In this setting, the correlation between observations from a given individual must be accounted for. Otherwise, downstream analyses may be affected by a number of factors, including false conclusions due to underestimated variance terms. Generalized estimating equation (GEE) models are an extension of generalized linear models devised to analyze data, which arise commonly in applied sciences [12] . Typically, in such cases, the data are collected in clusters in which observations within a cluster tend to be correlated, whereas observations in separate clusters are independent [12] . Clustered data arise from longitudinal or panel data, family studies or studies with spatial structures [12] . GEE models are used to test hypotheses about the dependent variables that are not necessarily normally distributed. The GEE framework is based on the concept of quasi-likelihood, which requires only an assumed relationship between the expected value (first moment) of the dependent variable and the covariates and between the conditional mean and variance (first and second moments) of the dependent variable. This approach uses the mean and variance of the response variable to derive the quasi-likelihood and its estimating equations [20] [21] . Because only the first two moments are needed, this approach is called semi-parametric. GEE allows for different choices of correlation structures, namely, independent, exchangeable, autoregressive, and unstructured. The method yields asymptotically unbiased and consistent estimates even if the incorrect working correlation structure is chosen. GEE accounts for within-individual correlations, allows for missing data, time-varying covariates, irregularly-timed or mistimed measurements. Fitting models using GEE requires the definition of the link function (which associates the linear predictor to the mean), an assumed distribution for the response variable, and a correlation structure of the response variable, often referred to as the working correlation matrix [12] .

We defined $\Delta T$ as ${T}_{s}-{T}_{a}$ , where ${T}_{s}$ is the surface temperature and ${T}_{a}$ is

the ambient temperature. Note that the $\frac{\text{d}{T}_{s}}{\text{d}{T}_{a}}$ introduced in the previous section

is proportional to $\Delta T$ , representing the rate of change of ${T}_{s}$ as ${T}_{a}$ changes. Our main goal was to describe the relationship between $\Delta T$ and ${T}_{a}$ for each body region. The relationship between $\Delta T$ and ${T}_{a}$ is linear when the body region is not capable of vasomotor adjustment and nonlinear when the body region is capable of vasomotor adjustment. Here, we considered as nonlinear the polynomial functions of degree 2 and 3. We used a two-step strategy to select models to describe variation in $\Delta T$ as a function of ambient temperature ( ${T}_{a}$ ) for each body region. In the first step, we selected a candidate model by minimizing the quasilikelihood under the Independence model Criterion (QIC) statistic over a grid comprised of correlation structures and model complexities (represented by increasing polynomial degrees). The QIC statistic is analogous to the familiar AIC (Akaike’s Information Criterion) statistic used for comparing models fitted with likelihood-based methods [22] . Because the generalized estimating equations (GEE) method is not a likelihood-based method, the AIC statistic is not available [22] . Here we applied the GEE method assuming a Normal distribution (identity link function), considered a third degree polynomial as the global model, and allowed for three possible correlation structures to be tested: independent, exchangeable, and unstructured. We chose as the best candidate the model with the minimum QIC. In the second step, we assumed the correlation structure as fixed (as defined in the first step) and assessed model parsimony (compromise between model complexity and residual variance) using a boostrap-based strategy. We defined the test statistic as ${\Delta}_{QIC}=QI{C}_{0}-QI{C}_{1}$ , where $QI{C}_{0}$ is the $QIC$ value of the simpler model and $QI{C}_{1}$ is the $QIC$ value of the more complex model. This allows us to test whether a complex model leads to a significantly decreased $QIC$ . If there are no gains in using a more complex model, we expect to observe ${\Delta}_{QIC}\le 0$ (null hypothesis). If there is a substantial decrease in $QIC$ when using the more complex model, we expect to observe ${\Delta}_{QIC}>0$ . Therefore, we used the distribution of ${\Delta}_{QIC}$ to decide whether or not we should use a more complex model vis-à-vis a simpler one. Because the theoretical distribution of ${\Delta}_{QIC}$ is unknown, we used the bootstrap to estimate its distribution empirically. Once the empirical bootstrap-based distribution of ${\Delta}_{QIC}$ is available, we used its quantiles to decide whether choosing a simpler model or not. In practical terms, if $0\in {q}_{\mathrm{2.5\%}}\mathrm{,}{q}_{\mathrm{97.5\%}}$ , then we do not have evidence to reject the hypothesis that there are no gains in using a more complex model.

2.3.3. Biophysical Modeling

Biophysical modeling of the heat exchange process allowed us to estimate the amount of heat exchanged (Q) with the ambient by each body region of the gracile mouse opossum. We represented the gracile mouse opossum by a simple geometric model: the ventral body region as a semi-cylinder, the back as a horizontal plate, the chest as a horizontal plate, the tail as a horizontal cylinder, the feet as a flat plate, the fingers as horizontal cylinders (data for fingers were added to that of the feet), and the ears as vertical plates. We measured the characteristic dimension, ${d}_{i}$ , and surface area, ${A}_{i}$ , of each i-th body region from thermal images of the gracile mouse opossum using the shape tools of software ImageJ (Figure 3).

We estimated heat exchange (Q) in Watts for each body surface ( ${Q}_{\text{ears}}$ , ${Q}_{\text{feet}}$ , ${Q}_{\text{tail}}$ , ${Q}_{\text{chest}}$ , ${Q}_{\text{dorsal}}$ , ${Q}_{\text{ventral}}$ ), where ${Q}_{\text{ears}}$ is the sum of the estimated heat exchange by the proximal and distal regions of the ear (multiplied by two to account for the two ears); ${Q}_{\text{feet}}$ is the sum of the estimated heat exchange by the lower and upper soles (multiplied by four to account for the four feet), plus the estimated heat exchange by the fingers (multiplied by 20 to account for the 20 fingers); ${Q}_{\text{tail}}$ is the estimated heat exchange by the tail; ${Q}_{\text{chest}}$ is the estimated heat exchange by the chest; ${Q}_{\text{dorsal}}$ is the estimated heat exchange by the dorsal region; and ${Q}_{\text{ventral}}$ is the estimated heat exchange by the ventral region. We assumed that the gracile mouse opossum was in thermal equilibrium with the ambient. We defined ${Q}_{\text{total}}$ for each gracile mouse opossum by adding the Q values for all body surfaces as,

${Q}_{\text{total}}={Q}_{\text{ears}}+{Q}_{\text{feet}}+{Q}_{\text{tail}}+{Q}_{\text{chest}}+{Q}_{\text{dorsal}}+{Q}_{\text{ventral}}.$ (2)

Figure 3. Infrared images of gracile mouse opossum within the experimental container. Solid lines represent geometric forms used for surface area estimation and dashed lines represent the characteristic dimension (d) of the body region (see values in Table 1).

Percent heat exchange by each i-th body region $\left({Q}_{i}\right)$ was calculated as $\frac{{Q}_{i}}{{Q}_{\text{total}}\times 100}$ . The ${Q}_{i}$ for each body region was estimated using the following equation [2] [7] [23] ,

${Q}_{i}={Q}_{r}+{Q}_{c},$ (3)

where ${Q}_{r}$ is the radiative heat exchange and ${Q}_{c}$ is the free convective heat exchange, such that,

${Q}_{r}=\u03f5\sigma A\left({T}_{s}^{4}-{T}_{a}^{4}\right),$ (4)

and

${Q}_{c}={h}_{c}A\left({T}_{s}-{T}_{a}\right)\mathrm{,}$ (5)

where
$\u03f5$ is the emissivity for biological tissues, assumed to be 0.96 according to Monteith and Unsworth [23] ;
$\sigma $ is the Stefan-Boltzman constant (5.6703 ´ 10^{−8}); A is the surface area (m^{2}) of the i-th body region;
${T}_{s}$ is the surface temperature of the i-th body region;
${T}_{a}$ is ambient temperature (K˚); and
${h}_{c}$ is the heat transfer coefficient given by

${h}_{c}=Nu\frac{k}{d},$ (6)

where k is the thermal conductivity of the air at a particular
${T}_{a}$ (W∙m^{−1}∙˚K). The relationship between k and
${T}_{a}$ was estimated by Tattersall et al. [2] as,

$k=0.00241+7.5907{e}^{-6}\times {T}_{a}.$ (7)

$Nu$ is the dimensionless Nusselt number that is a measure of the ratio of buoyant to viscous forces. In free convection, the Nussel number is a function of the dimensionless Grashof (Gr) and Prandtl (Pr) numbers, written as: $Nu=f\left(Gr\mathrm{,}Pr\right)$ . In forced convection, the Nussel number is a function of the dimensionless Reynolds (Gr) and Prandtl (Pr) numbers: $Nu=f\left(Re\mathrm{,}Pr\right)$ . Grashof, Prandtl, and Reynolds numbers are given by,

$Gr=\frac{ag{d}^{3}\left({T}_{s}-{T}_{a}\right)}{{v}^{2}}\mathrm{,}$ (8)

$Re=\frac{Vd}{v}\mathrm{,}$ (9)

$Pr=\frac{{c}_{p}\mu}{k}\mathrm{,}$ (10)

where a is the coefficient of thermal expansion of air, g is the acceleration of gravity, d is the characteristic dimension of the body region (Table 1), v is the kinematic viscosity of air, ${c}_{p}$ is the specific heat of air, $\mu $ is the dynamic viscosity of air, and k is the thermal conductivity of air. Experimental evidence shows that the free-convection function for heat exchange should be used when $Gr>16R{e}^{2}$ and the forced-convection function when $Gr<0.1R{e}^{2}$ [24] . We obtained values for the Gr number that were under the range of validity available on the literature [23] [24] . Therefore, we modeled heat exchange by forced

Table 1. Biophysical model parameters, coefficients, and dimensionless numbers across the different body regions of the gracile mouse opossum.

convection considering a wind velocity of 0.1 ms^{−1} according to Tattersall et al. [2] , Gates [24] , and Greenberg et al. [9] , to obtain the values for Gr and Re. The relationship between Nu and Gr, Re, and Pr has been determined empirically for a range of geometric shapes and flow regimes (see Table 1 for the relationships for each body region). According to Monteith and Unsworth [23] and Gates [24] we assumed,

$Nu=cR{e}^{n}\mathrm{,}$ (11)

where c and n are constants related to shape.

Data in table are the shape and characteristic dimensions used for heat exchange calculations, calculated d and surface area, and relationship between Grashof (Gr), Reynolds (Re), and Nusselt (Nu) numbers [23] [24] . (Gr), (Re), and (Nu) were calculated according to the equations shown in the Biophysical Modeling subsection.

Assuming that each gracile mouse opossum was in thermal balance during the experimental procedure we should expect that the relationship between ${Q}_{\text{total}}$ and ${T}_{a}$ would follow the relationship between metabolic heat production and ${T}_{a}$ . Based on the rates for oxygen uptake determined by Cooper et al. [18] , we estimated the total metabolic heat production (Watts) for a gracile mouse opossum with a body mass of 40 g at ${T}_{a}=12\u02da\text{C}$ , 20˚C, 28˚C, 30˚C, and 32˚C. Thereafter, we compared these rates with the total heat exchange rate (Watts) estimated by the biophysical model ( ${Q}_{\text{total}}$ ).

2.3.4. Statistical Analysis for the Relationship between Heat Exchange of Each Body Region (Q_{i}) and (T_{a})

To estimate the role of each body region for heat exchange we defined the variable ${Q}_{i}$ as the percent heat exchange of each i-th body region. ${Q}_{i}$ was

calculated as $\frac{{Q}_{i}}{{Q}_{\text{total}}}\times 100$ . The relationship between ${Q}_{i}$ and ${T}_{a}$ was described

by linear regression models. We followed the same framework of Generalized Estimating Equations (GEE) as the estimation method for the model to account for correlated response.

3. Results

Our primary objective is to search for signatures of linearity or non-linearity in the relationship between the differential $\Delta T={T}_{s}-{T}_{a}$ and ambient temperature ( ${T}_{a}$ ). Linearity or nonlinearity between $\Delta T$ and ${T}_{a}$ arise as a consequence of the absence or presence of cutaneous vasodilation and vasoconstriction, respectively, which are major physiological mechanisms of vasomotor control, modulating heat exchange with the ambient [1] [3] [4] . Furthermore, the body surface of an endotherm is not homogeneous with respect to vasomotor control. Therefore, surface temperatures were measured in individuals of the mammalian experimental model at a range of ambient temperatures for different body regions. The data we gathered and the relevant statistical analyses are described below.

Within the range of ${T}_{a}$ from 8˚C to 38˚C, surface temperatures ( ${T}_{s}$ ) for all body regions of the gracile mouse opossum increased with ${T}_{a}$ at increasing rates depending on the body region (Figure 4). For furred body regions such as the back, chest, and ventral region, we observed a linear relationship between $\Delta T$ and ${T}_{a}$ . Conversely, for furless body regions such as the ears, feet, and tail, the relationship between $\Delta T$ and ${T}_{a}$ was nonlinear (Table 2). The nonlinear relationships observed for feet and tail are typical of thermal windows showing a vasodilation at ${T}_{a}>20\u02da\text{C}$ and vasoconstriction at ${T}_{a}<20\u02da\text{C}$ (Figure 4). On the other hand, the nonlinear relationship observed for ears has a close to linear relationship at ${T}_{a}<35\u02da\text{C}$ and a change in this tendency at ${T}_{a}>35\u02da\text{C}$ , probably due to the equalization between ${T}_{s}$ and ${T}_{a}$ .

Figure 4. Variation in surface temperature of the gracile mouse opossum across the different body regions as a function of ${T}_{a}$ . Red lines represent the relationship between $\Delta T$ versus ${T}_{a}$ . Confidence interval bands are represented in grey. Grey dots represent the data for gracile mouse opossum individuals.

Table 2. Polynomial models selected for describing the relationship between $\Delta T$ and ${T}_{a}$ across the different body regions of the gracile mouse opossum.

*indicates which candidate model was selected. If $0\in {q}_{\mathrm{2.5\%}}\mathrm{,}{q}_{\mathrm{97.5\%}}$ , then we do not have evidence to reject the hypothesis that there are no gains in using a more complex model.

The comparison between the rates of heat production estimated from Cooper et al. [18] and the total heat exchange (present study) for the gracile mouse opossum (Figure 5) revealed that our estimates were consistently lower than the rates measured by Cooper et al. [18] . At ${T}_{a}=12\u02da\text{C}$ , the estimated heat exchange ( ${Q}_{\text{total}}$ ) is approximately 97% of the predicted metabolic heat production. This percentage tends to get lower as ${T}_{a}$ increases, 74% at ${T}_{a}=20\u02da\text{C}$ , 83% at ${T}_{a}=28\u02da\text{C}$ , 57% at ${T}_{a}=30\u02da\text{C}$ , and 35% at ${T}_{a}=32\u02da\text{C}$ .

Heat exchanged by each body region varied across the range of ambient temperature as indicated by the significant interactions between ${Q}_{i}$ and ${T}_{a}$ (Table 3 and Figure 6). However, the pattern of variation was considerably different between furred and furless body regions. Furred regions such as the chest, dorsal, and ventral regions exchanged most of the heat, from 50% to 75%, dissipated by the gracile mouse opossum at lower temperatures. The furred ears account for about 15% of the total heat exchange only at mid-ambient temperature. As ambient temperature gets higher, furless regions such as feet and tail become more important and can account for almost 70% of the total heat exchange.

Figure 5. Variation in heat exchange (as the ${Q}_{i}$ percentage) of each i-th body region of the gracile mouse opossum as a function of ${T}_{a}$ . Red lines represent the relationship between percentage of heat exchange versus ${T}_{a}$ . Grey areas represent confidence intervals. Grey dots represent the data of gracile mouse opossum individuals.

Table 3. Coefficient estimation from the model describing the relationship between ${Q}_{i}$ and ${T}_{a}$ across the different body regions of the gracile mouse opossum.

*coefficients statistically significant at p-value £ 0.05; ***coefficients statistically significant at p-value £ 0.001.

Figure 6. Variation in total heat exchange of the gracile mouse opossum ( ${Q}_{\text{total}}$ ) as a function of ${T}_{a}$ . Blue line represents the relationship between total heat exchange due to convection and radiation versus ${T}_{a}$ . Black line represents predicted metabolic rate of gracile mouse opossum based on data from Cooper et al. [18] . Grey area represents confidence intervals.

4. Discussion

Heat exchange in the surface of an endotherm is a process that can be modeled by the equation ${Q}_{\text{total}}=\u03f5\sigma A\left({T}_{s}^{4}-{T}_{a}^{4}\right)+{h}_{c}A\left({T}_{s}-{T}_{a}\right)$ , as described in classic textbooks and review articles of animal physiology and biophysical ecology [10] [24] [25] [26] . From this equation, it becomes evident that heat exchange is determined by the surface temperature, ${T}_{s}$ , the gradient between surface temperature and ambient temperature ${T}_{s}-{T}_{a}$ , the superficial area, A, of the body region, and the intrinsic characteristics of a given body region that define the amount of heat transported to the surface. The surface of an organism, however, is not uniform with respect to all these characteristics [4] [26] . Therefore, in order to understand the heat exchange process as a function of ambient temperature across non-uniform surfaces in endotherms one must: 1) measure ${T}_{s}$ as function of ${T}_{a}$ across the surface of the organism; 2) assess the relationship between ${T}_{s}$ and ${T}_{a}$ across the surface of the organism; 3) understand how ${T}_{s}$ and ${T}_{a}$ interact with each body region to determine the role of each body region for total heat exchange at each ${T}_{a}$ [2] [6] [25] [26] .

Different body regions of the gracile mouse opossum exhibited distinct heat exchange properties in response to variation in ambient temperature. As a consequence, the relative contribution of each body region to total heat balance also changed with ambient temperature. Whereas furred body regions such as the chest, dorsal, and ventral region, exhibited a linear relationship between $\Delta T$ and ${T}_{a}$ , furless body regions such as the ears, feet, and tail exhibited a nonlinear relationship between $\Delta T$ and ${T}_{a}$ . The linear relationship between $\Delta T$ and ${T}_{a}$ for the furred regions of G. microtarsus indicates that mechanisms other than vasomotor adjustment dominated the dynamics of heat exchange at these regions. Indeed, this is readily observed by the decay of $\Delta T$ as the ambient temperature approaches ${T}_{b}$ (around 30.6˚C and 34.7˚C, [18] [27] ), or even surpasses it. Although quite similar in their general behavior, furred regions also varied in the relationship between $\Delta T$ and ${T}_{a}$ . The chest and ventral regions always had $\Delta T$ values considerably higher than those recorded for the dorsal region for most ${T}_{a}$ considered. Also, the $\Delta T$ for the chest and ventral regions exhibited more variation as a function of ${T}_{a}$ (i.e., a higher slope value) compared to the dorsal region. Together, these differences indicate a better insulation of the dorsal region compared to the chest and ventral regions and, possibly, differences in the capacity to modulate body insulation by piloerection. These interpretations also supported by the results on heat exchange.

For the range of temperatures tested, temperature equalization (i.e., $\Delta T=0$ ) for the furred body regions occurred only at the highest temperatures tested (differently from what happened with the furless body regions). At these temperatures (>35˚C), ambient temperature most likely had already surpassed internal body temperature. Therefore, we suspect that, in this case, fur insulation of the animals reduced heat gain from the environment at greater rates, shifting the point of temperature equalization to a temperature a few degrees higher than their normal body temperature. As at these temperatures energy expensive heat loss mechanisms are likely to intervene, preventing extra heat gain from the ambient might be highly relevant for the total heat balance and energetics of G. microtarsus. Finally, we should acknowledge the possibility that the atainement of equalization at higher than expected temperatures (normal ${T}_{b}$ ) could be attributable to hyperthermia. We did not monitored internal body temperature of our animals during the experiments. Nevertheless, the results from temperature equalization obtained for the furless regions, particularly the ear, suggest that the animals did not experience hyperthermia. Thus, we believe that the linear relationship between $\Delta T$ and ${T}_{a}$ for the furred regions of G. microtarsus reveals not only the insulative properties of the fur in preventing excessive heat loss at low temperatures, but also the less commonly acknowledged role of preventing excessive heat gain at higher temperatures.

Furless body regions such as the ears, feet, and tail exhibited a nonlinear relationship between $\Delta T$ and ${T}_{a}$ . However, there were differences among different furless regions. The relationship found for the ears was quite interesting. For ${T}_{a}$ higher than 35˚C, $\Delta T$ tended to zero and remained unchanged. Our interpretation of this result is that around 35˚C ${T}_{a}$ equalized with body temperature [18] [27] , which induced a vasomotor response to limit of blood flow to the ears, thus preventing more heat to be gained through this region. Interestingly, below 35˚C and even at the lowest temperature tested, we found no indication of vasoconstriction for the ears, which could have been expected as a heat conservation strategy under cooler conditions. The dynamics of heat exchange variation in response to changes in ${T}_{a}$ for the ears was quite different from that observed for the other furless regions, such as the feet and tail. Feet and tail both exhibited clear vasomotor responses (Figure 4). For feet, the threshold for vasodilation was detected at 20˚C and increased until around 30˚C, declining thereafter. For temperatures below 20˚C, interestingly, there was also a vasodilation response for the feet indicating that the temperature of this body region was, to some extense, defended from external cooling. The tail had a vasodilation threshold a few degrees higher than the feet, around 24˚C, that increased to ambient temperatures around 32˚C, declining thereafter. Differently from the feet, however, it seems that the vasoconstriction response of the tail was maintained as the animals were cooled below 24˚C. In this regard, the tail behaved as a typical thermal window, promoting non-evaporative heat dissipation at higher temperatures and heat conservation as ambient temperature is lowered [2] [28] [29] .

The differences in the heat exchange properties observed among the different body regions of G. microtarsus had a direct impact on the contribution of each body region to the total heat balance of the gracile mouse opossum. Furred regions were the largest in surface area and collectively responded for more than half of the total heat dissipated, especially at temperatures below 20˚C - 24˚C. However, as these regions lacked the capacity for vasomotor adjustment, as temperatures were elevated above that limit, the contribution of the furless regions of the feet and tail to the total heat balance of the animals gained prominence. Indeed, from almost negligible at low temperatures, the contribution of the feet and tail was elevated to more than half of the total heat exchanged at the highest temperature tested. This effect can be directly attributed to the vasomotor response exhibited by these body regions. Also, it is conceivable that some evaporative cooling mechanism concurrently improved heat dissipation through tail and feet at temperatures above the vasomotor threshold. At this point, it remains uncertain which evaporative heat dissipation mechanisms are employed by G. microtarsus.

Whereas examining heterothermy across different body regions can be straightforward with IRT, addressing the question of how ${T}_{s}$ can inform about ${T}_{b}$ regulation and the existence of thermal windows lacks a sound statistical approach. Here, we have shown that combining the use of polynomial functions and estimation by GEE can be a robust and informative strategy for addressing questions on the dynamics of either ${T}_{s}$ or surface heat exchange across animal body regions. Polynomial functions capture the linearity and nonlinearity in the relationship between ${T}_{s}$ and ${T}_{a}$ , making the task of detecting regions capable or not of vasomotor adjustments objective and unbiased. In addition, the GEE estimation method used with generalized linear models incorporates structures typical of IRT data, such as correlated responses, within-individual correlations, allows for missing data, time-varying covariates, irregularly-timed or mistimed measurements [12] .

Tail and feet comprise less than 25% of G. microtarsus body surface although their contribution to total heat exchange at temperatures above 20˚C - 24˚C was considerable. A similar response is found in rabbits and elephant ears [28] [29] [30] [31] , guanaco fur [32] , penguin feet [33] , bird bills [2] [34] , digits [35] , combs and wattles of several fowl [36] , and ungulate horns and antlers [37] [38] . The tail, in particular, is commonly acknowledged as an important organ to modulate radiative heat dissipation via vasomotor adjustments in rats [39] [40] [41] . However, the potential for controlling heat exchange rate through the tail is probably higher in rodents than in G. microtarsus, which has semi-arboreal habits and uses its tail for climbing [14] [42] . We suspect that conflicting functional demands between thermoregulation and climbing may explain the difference in the potential to control heat exchange through the tail in rodents and the gracile mouse opossum. A similar situation seems to occur with young toco-toucans, in which the vasomotor control of heat exchange through the bill is thought to conflict with a demand associated with bill growth [2] . We are not aware of other studies focusing on the heat exchange properties of other small marsupials and/or bearers of prehensile tails. Therefore, it remains uncertain whether a functional trade-off in tail use may occur in G. microtarsus or whether it reflects a more general difference between different mammalian groups.

Our estimates of total heat exchange for the gracile mouse opossum were consistently lower than the rates of heat production estimated by Cooper et al. [18] using metabolic measurements (Figure 6). Assuming that the gracile mouse opossum reached a thermal steady state at each ${T}_{a}$ tested, this comparison shows that our estimates of heat dissipation of gracile mouse opossum using biophysical modeling are good, given that they balance its internal heat production. The difference observed between the two curves might be attributable to mechanisms of heat exchange, such as evaporative heat exchange, which were not considered in our model.

5. Conclusion

In conclusion, the process of heat exchange across an endotherm surface as a function of ambient temperature is complex and depends on as many dimensions of the organism and its surrounding ambient. Assessing this process under simplified conditions within the unifying principles of thermoregulation unveils and clarifies the rich diversity of mechanisms underlying the process. Under this approach, one is able to better describe how the dynamics of heat exchange mechanisms vary between body regions, allowing endotherms to thermoregulate across ambient temperatures. This approach is based on detailed understanding of the relative importance of each body region for the entire budget across the range of ambient temperature.

Acknowledgements

We are very grateful to Marcus A. M. de Aguiar, Eduardo G. Martins, Luiza C. Duarte, and Mathias M. Pires for carefully reading the first drafts and improving the clarity of the manuscript. João Del Giudice Neto and Marcos Mecca Pinto from Mogi-Guaçu Biological Reserve for logistical support. Research supported by Fundação de Amparo à Pesquisa do Estado de São Paulo, Brazil (FAPESP) [Denis O. V. de Andrade - 2007/05080-1 and 2013/04190-9; Sérgio F. dos Reis - 2005/51353-4]; Conselho Nacional de Desenvolvimento Científico e Tecnológico, Brazil (CNPq) [Barbara Henning - 140773/2013-4; Denis O. V. de Andrade - 302045/2012-0 and 306811/2015-4; Sérgio F. dos Reis - 303544/2011-2]; and CAPES [Barbara Henning - 99999.004965/2014-00].

References

[1] Schmidt-Nielsen, K. (1997) Animal Physiology: Adaptation and Environment. Cambridge University Press, Cambridge.

[2] Tattersall, G.J., Andrade, D.V. and Abe, A.S. (2009) Heat Exchange from the Toucan Bill Reveals a Controllable Vascular Thermal Radiator. Science, 325, 468-470.

https://doi.org/10.1126/science.1175553

[3] Kingma, B.R., Frijns, A.J., Saris, W.H., van Steenhoven, A.A. and van Marken Lichtenbelt, W.D. (2012) Mathematical Modeling of Human Thermoregulation: A Neurophysiological Approach to Vasoconstriction. In: Madani, K., Dourado Correia, A., Rosa, A. and Filipe, J., Eds., Computational Intelligence, Studies in Computational Intelligence, Springer Verlag, Berlin.

[4] Wang, Y., Kimura, K., Inokuma, K.I., Saito, M., Kontani, Y., Kobayashi, Y., Mori, N. and Yamashita, H. (2006) Potential Contribution of Vasoconstriction to Suppression of Heat Loss and Homeothermic Regulation in UCP1-Deficient Mice. Pflügers Archiv, 452, 363-369.

https://doi.org/10.1007/s00424-005-0036-3

[5] Boldrini, J.L., Viana, M.P., Reis, S.F. and Henning, B. (2018) A Mathematical Model for Thermoregulation in Endotherms including Heat Transport by Blood Flow and Four Thermal Feedback Control Mechanisms: Changes in Coat, Metabolic Rate, Blood Flux, and Ventilation Rate. (In Preparation)

[6] McCafferty, D.J., Gilbert, C., Paterson, W., Pomeroy, P., Thompson, D., Currie, J. and Ancel, A. (2011) Estimating Metabolic Heat Loss in Birds and Mammals by Combining Infrared Thermography with Biophysical modeling. Comparative Biochemistry and Physiology A, 158, 337-345.

https://doi.org/10.1016/j.cbpa.2010.09.012

[7] McCafferty, D.J., Gilbert, C., Thierry, A.M., Currie, J., Le Maho, Y. and Ancel, A. (2013) Emperor Penguin Body Surfaces Cool below Air Temperature. Biology Letters, 9, 20121192.

https://doi.org/10.1098/rsbl.2012.1192

[8] Phillips, P.K. and Sanborn, A.F. (1994) An Infrared, Thermographic Study of Surface Temperature in Three Ratites: Ostrich, Emu and Double-Wattled Cassowary. Journal of Thermal Biology, 19, 423-430.

https://doi.org/10.1016/0306-4565(94)90042-6

[9] Greenberg, R., Cadena, V., Danner, R.M. and Tattersall, G. (2012) Heat Loss May Explain Bill Size Differences between Birds Occupying Different Habitats. PLoS ONE, 7, e40933.

https://doi.org/10.1371/journal.pone.0040933

[10] Tattersall, G.J. (2016) Infrared Thermography: A Non-Invasive Window into Thermal Physiology. Comparative Biochemistry and Physiology A, 202, 78-98.

https://doi.org/10.1016/j.cbpa.2016.02.022

[11] Sniegowski, M.C., Erlanger, M. and Olson, J. (2017) Thermal Imaging of Corneal Transplant Rejection. International Ophthalmology, 2017, 1-5.

https://doi.org/10.1007/s10792-017-0731-z

[12] Liang, K.Y. and Zeger, S.L. (1986) Longitudinal Data Analysis Using Generalized Linear Models. Biometrika, 73, 13-22.

https://doi.org/10.1093/biomet/73.1.13

[13] Lubke, G.H., Campbell, I., McArtor, D., Miller, P., Luningham, J. and van den Berg, S.M. (2016) Assessing Model Selection Uncertainty Using a Bootstrap Approach: An Update. Structural Equation Modeling: A Multidisciplinary Journal, 24, 230-245.

https://doi.org/10.1080/10705511.2016.1252265

[14] Martins, E.G., Bonato, V., Pinheiro, H. and Reis, S.F. (2006) Diet of the Gracile Mouse Opossum (Gracilinanus microtarsus) (Didelphimorphia: Didelphidae) in a Brazilian Cerrado: Patterns of Food Consumption and Intrapopulation Variation. Journal of Zoology, 269, 21-28.

https://doi.org/10.1111/j.1469-7998.2006.00052.x

[15] Martins, E.G., Bonato, V., Da-Silva, C. and Reis, S.F. (2006) Partial Semelparity in the Neotropical Didelphid Marsupial Gracilinanus microtarsus. Journal of Mammalogy, 87, 915-920.

[16] Martins, E.G., Bonato, V., Da-Silva, C. and Reis, S.F. (2006) Seasonality in Reproduction, Age Structure and Density of the Gracile Mouse Opossum Gracilinanus microtarsus (Marsupialia: Didelphidae) in a Brazilian Cerrado. Journal of Tropical Ecology, 22, 461-468.

https://doi.org/10.1017/S0266467406003269

[17] Martins, E.G., Araújo, M.S., Bonato, V. and Reis, S.F. (2008) Sex and Season Affect Individual-Level Diet Variation in the Neotropical Marsupial Gracilinanus microtarsus (Didelphidae). Biotropica, 40, 132-135.

[18] Cooper, C., Withers, P. and Cruz-Neto, A. (2009) Metabolic, Ventilatory, and Hygric Physiology of the Gracile Mouse Opossum (Gracilinanus agilis). Physiological and Biochemical Zoology, 82, 153-162.

https://doi.org/10.1086/595967

[19] Ronchi, J.A., Henning, B., Ravagnani, F.G., Figueira, T.R., Castilho, R.F., Reis, S.F. and Vercesi, A.E. (2015) Increased Susceptibility of Gracilinanus microtarsus Liver Mitochondria to Ca2+-Induced Permeability Transition Is Associated with a More Oxidized State of NAD(P). Oxidative Medicine and Cellular Longevity, Article ID: 940627.

[20] Wedderburn, R.W. (1974) Quasi-Likelihood Functions, Generalized Linear Models, and the Gauss-Newton Method. Biometrika, 61, 439-447.

[21] McCullagh, P. (1983) Quasi-Likelihood Functions. Annals of Statistics, 11, 59-67.

https://doi.org/10.1214/aos/1176346056

[22] Pan, W. (2001) Akaike’s Information Criterion in Generalized Estimating Equations. Biometrics, 57, 120-125.

https://doi.org/10.1111/j.0006-341X.2001.00120.x

[23] Monteith, J. and Unsworth, M. (2008) Principles of Environmental Physics. Academic Press, New York.

[24] Gates, D.M. (2012) Biophysical Ecology. Springer Verlag, New York.

[25] Randall, D., Burggren, W.W. and French, K. (2002) Eckert Animal Physiology. Macmillan, New York.

[26] Hill, R.W., Wyse, G. and Anderson, M. (2012) Animal Physiology. Sinauer Associates, Sunderland.

[27] Morrison, P. and McNab, B. (1962) Daily Torpor in a Brazilian Murine Opossum (Marmosa). Comparative Biochemistry and Physiology A, 6, 57-68.

https://doi.org/10.1016/0010-406X(62)90043-9

[28] Hill, R.W. and Veghte, J.H. (1976) Jackrabbit Ears: Surface Temperatures and Vascular Responses. Science, 194, 436-438.

https://doi.org/10.1126/science.982027

[29] Phillips, P.K. and Heath, J.E. (1992) Heat Exchange by the Pinna of the African Elephant (Loxodonta africana). Comparative Biochemistry and Physiology A, 101, 693-699.

https://doi.org/10.1016/0300-9629(92)90345-Q

[30] Hill, R.W., Christian, D.P. and Veghte, J.H. (1980) Pinna Temperature in Exercising Jackrabbits, Lepus californicus. Journal of Mammalogy, 61, 30-38.

https://doi.org/10.2307/1379954

[31] Mohler, F.S. and Heath, J.E. (1988) Comparison of IR Thermography and Thermocouple Measurement of Heat Loss from Rabbit Pinna. American Journal of Physiology, Regulatory, Integrative and Comparative Physiology, 254, R389-R395.

https://doi.org/10.1152/ajpregu.1988.254.2.R389

[32] Morrison, P. (1966) Insulative Flexibility in the Guanaco. Journal of Mammalogy, 47, 18-23.

https://doi.org/10.2307/1378061

[33] Wilson, R.P., Adelung, D. and Latorre, L. (1998) Radiative Heat Loss in Gentoo Penguin (Pygoscelis papua) Adults and Chicks and the Importance of Warm Feet. Physiological Zoology, 71, 524-533.

https://doi.org/10.1086/515955

[34] Hagan, A.A. and Heath, J.E. (1980) Regulation of Heat Loss in the Duck by Vasomotion in the Bill. Journal of Thermal Biology, 5, 95-101.

https://doi.org/10.1016/0306-4565(80)90006-6

[35] Moritz, G.L. and Dominy, N.J. (2012) Thermal Imaging of Aye-Ayes (Daubentonia madagascariensis) Reveals a Dynamic Vascular Supply during Haptic Sensation. International Journal of Primatology, 33, 588-597.

https://doi.org/10.1007/s10764-011-9575-y

[36] Buchholz, R. (1996) Thermoregulatory Role of the Unfeathered Head and Neck in Male Wild Turkeys. The Auk, 113, 310-318.

https://doi.org/10.2307/4088897

[37] Picard, K., Thomas, D.W., Festa-Bianchet, M. and Lanthier, C. (1994) Bovid Horns: An Important Site for Heat Loss during Winter? Journal of Mammalogy, 75, 710-713.

https://doi.org/10.2307/1382520

[38] Taylor, C.R. (1966) The Vascularity and Possible Thermoregulatory Function of the Horns in Goats. Physiological Zoology, 39, 127-139.

https://doi.org/10.1086/physzool.39.2.30152426

[39] Rand, R.P., Burton, A.C. and Ing, T. (1965) The Tail of the Rat, in Temperature Regulation and Acclimatization. Canadian Journal of Physiology and Pharmacology, 43, 257-267.

https://doi.org/10.1139/y65-025

[40] Young, A.A. and Dawson, N.J. (1982) Evidence for On-Off Control of Heat Dissipation from the Tail of the Rat. Canadian Journal Physiology and Pharmacology, 60, 392-398.

https://doi.org/10.1139/y82-057

[41] Dawson, N.J. and Keber, A.W. (1979) Physiology of Heat Loss from an Extremity: The Tail of the Rat. Clinical Experimental Pharmacology and Physiology, 6, 69-80.

https://doi.org/10.1111/j.1440-1681.1979.tb00009.x

[42] Costa, L.P., Leite, Y.L.R. and Patton, J.L. (2003) Phylogeography and Systematic Notes on Two Species of Gracile Mouse Opossums, Genus Gracilinanus (Marsupialia: Didelphidae) from Brazil. Proceedings of the Biological Society of Washington, 116, 275-292.