A mathematical model for thermoregulation in endotherms including heat transport by blood flow and thermal feedback control mechanisms: changes in coat, metabolic rate, blood fluxes, ventilation and sweating rates

ABSTRACT Thermoregulation in endotherms allows the maintenance of the body temperature independent of ambient temperature. Experimental data have revealed complex interactions between the physiological mechanisms of thermoregulation and environmental conditions. We derive a nonlinear partial integro-differential dynamical model based on physical first principles and fundamental physiological mechanisms to understand the role of some thermal control mechanisms in the thermoregulation process of endotherms. The model is composed of four layers representing different tissues and it incorporates six thermal feedback control mechanisms. These mechanisms are heat production due to metabolic rate and heat exchange within the body given its internal structure, and the model considers heat exchange due to conduction, heat transport by blood flow, heat exchange with the ambient through convection, radiation, and evaporation from the respiratory tract and superficial evaporation in both passive and active situations. Our model sheds new light on previous explanations about the classic metabolism-ambient temperature U-shaped curve.


Introduction
Thermoregulation is a fundamental physiological process with implications not only for the organismal level but also, from an ecological standpoint, as a process in which organisms interact with environmental factors to define a species thermal niche (Hill, Wyse, & Anderson, 2012). Most mammals and birds maintain a relatively constant body temperature independent of environmental temperatures by means of homeostatic mechanisms that regulate rates of heat production and heat loss (Hill et al., 2012;Randall & Burggren, 2001). Therefore, mammals and birds are known as homeothermic endotherms. The uncoupling between body temperature and environmental conditions has allowed mammals and birds to colonize habitats that otherwise would not be suitable, thus expanding their geographical distribution over evolutionary time (Khaliq, Hof, Prinzinger, Böhning-Gaese, & Pfenninger, 2014;Kronfeld-Schor & Dayan, 2013).
The basic experimental facts concerning thermoregulation in homeothermic endotherms are known since the late eighteen hundreds (Scholander, Hock, Walters, Johnson, & Irving, 1950). An essential experimental result, which has been demonstrated many times for several species of birds and mammals, is the relationship between total internal heat produced due to metabolism, known as metabolic rate, and ambient temperature (Hill et al., 2012;Randall & Burggren, 2001). The shape of this relationship is captured in the socalled metabolism-ambient temperature U-shaped curve (Figure 1), which forms the basis for understanding the dynamics of the rate of heat transfer between an animal and its environment (Hill, Muhich, & Humphries, 2013;Hill et al., 2012;Randall & Burggren, 2001). This curve is usually obtained by measuring oxygen consumption and its shape is consistent across endotherms. Metabolic rate scales with body mass and ambient temperature according to a power function of the form Q = i 0 B 3/4 e −E/kT , where Q is the metabolic rate, i 0 is a normalization constant independent of body size and temperature, B is the body size, E is the activation energy, k = 1.38 × 10 −23 is Boltzmann constant, and T is the absolute temperature in K (Brown, Gillooly, Allen, Savage, & West, 2004;Darveau, Suarez, Andrews, & Hochachka, 2002;Gillooly, Charnov, West, Savage, & Brown, 2002;Heldmaier, Ortmann, & Elvert, 2004;Kleiber, 1932;Scholander et al., 1950). This empirically derived curve, remarkably continues to inspire conceptual and theoretical interest.
Despite the importance of this curve, its shape is commonly explained by a rather simple model that dates back to the early 1900s (Harris & Benedict, 1919;Kleiber, 1932). In these works, the authors applied a specially simple form of Fourier's law to heat flow in the body of homeothermic animals. The animal's body is represented by a core with a constant temperature T bc , surrounded by an insulating layer with heat conductivity k, thickness h, surface area S, and surface temperature T . According to the model (Kleiber, 1972), the rate of heat flow, dH/dt, may be written as 1 The term C ≡ k(S/h) is known as the average body condutance. Under the assumption that all the internal heat produced is exchanged with the ambient, i.e. Q = dH/dt, Equation (1) became the classical model for explaining the shape of the curve shown in Figure 1, which can be found nowadays in textbooks (e.g. Hill et al., 2012.) In the next section we review this classical model to expose its hidden simplifying hypotheses and its shortcomings. We conclude that this model is too simple to provide a satisfactory explanation for experimental data (Barnes & Buck, 2000;Brown & Lasiewski, 1972;Cooper, Withers, & Cruz-Neto, 2009;Gordon, 1990;Hart, 1957;Hill et al., 2012;Johnson, 1968;Scholander et al., 1950;Tattersall, 2016) that shapes the curve shown in Figure 1. In fact, the classical model assumes heat conduction as the only relevant heat transport mechanism in the animal's body, not taking into account heat transport due to blood flow. In addition, the model does not incorporate important nonlinear effects due to heat exchange by radiation, which is a heat exchange mechanism that accounts for at least 50% of total heat exchange of organisms (Gates, 2012, p. 24).

Figure 1.
Classic metabolism-ambient temperature U-shaped curve valid for mammals and birds. The thermoneutral zone (TNZ) is the range of ambient temperatures (T a ) within which an animal's resting metabolic rate is independent of ambient temperature and constant. The lower and upper limits of the TNZ are the lower (T LC ) and upper (T UC ) critical temperatures. Within the TNZ, body temperature (T bc ) is maintained by the modulation of body insulation, which is low near T UC and high near T LC , and by thermal control mechanism of vasomotor adjustment. Below the TNZ, body insulation is high and approximately constant, thus the body temperature is maintained with mechanisms that increase the rate of metabolic rate as ambient temperature falls. The absolute value of the slope of the curve at T a < T LC is the conductance C. Thus, for ambient temperatures below T UC we have basically dry heat exchange. Above the TNZ, body insulation is low and approximately constant, thus the body temperature is maintained by mechanisms of evaporative heat exchange that require an increase in the metabolic rate as ambient temperature rises.
In addition to the metabolism-ambient temperature curve, animal surface temperature as a function of ambient temperature is another commonly investigated relationship in thermoregulation studies. This relationship can be close to linear such as that observed for the bill surface temperature in sparrows (Greenberg, Cadena, Danner, & Tattersall, 2012), whereas this relationship can also be nonlinear such as that observed for the bill surface temperature in toucans (Tattersall, Andrade, & Abe, 2009). The linear dependence has positive slopes with values always smaller than one (Greenberg et al., 2012). The nonlinear dependencies have been explained as a consequence of vasomotor adjustment, which works as a thermoregulatory mechanism (Tattersall et al., 2009). However, to our knowledge, there is no clear theoretical explanation based on physics first principles for all these experimental results.
In this work, we introduce a four-layer model for heat production and exchange based on physics first principles. Our model considers the internal structure of the body, different types of internal heat transport, different mechanisms for heat exchange with the ambient (including nonlinear ones), as well as control mechanisms used for thermoregulation.
Our goal is to shed light on the roles played by specific thermal control mechanisms in the process of maintaining a prescribed temperature in the body core, which is defined as the set of the most important organs of an animal (Hill et al., 2012;Randall & Burggren, 2001). As for thermal control mechanisms, we consider (i) changes in metabolic rate, (ii) changes in insulation layer characteristics (e.g. raising of fur or feathers), (iii) changes in blood fluxes due to blood vasomotor adjustment in tissues near the body core, (iv) changes in blood fluxes due to blood vasomotor adjustment in tissues near the surface, (v) evaporation from the respiratory tract (e.g. panting), and (vi) superficial evaporation (e.g. sweating).
The thermal control mechanisms (i) and (ii) are considered in the classical model by means of dQ b /dt and h, respectively. To take (iii) and (iv) into account, we use an approach similar to Pennes' bio-heat transfer (BHT) model, which has been successfully used in studies of heat flux in vascularized tissues, such as hyperthermia and hypothermia treatments (Arkin, Xu, & Holmes, 1994;Okajima, Maruyama, Takeda, & Komiya, 2009). Mechanisms (v) and (vi) are important at high ambient temperatures. In particular for animals with no significant number of sweat glands, that rely mostly on evaporation from the respiratory tract, to take this type of evaporation into account, we assume its rate to be proportional to the volume of air breathed per unit of time and it increases by rapid shallow breathing, as seen typically in dogs and domestic fowl, and known as panting (Gates, 2012;Hill et al., 2012;Richards, 1970).
We model animals under experimental laboratory conditions, i.e. animals at rest, without direct exposure to sun light or reflected radiation, controlled ambient temperature and low humidity. Thus, for the sake of simplicity, we did not include control by postural changes. In Section 3.6, we provide further comments on how postural changes could be considered as an additional thermal control mechanism and how the effect of wind speed could be incorporated into our model.
Our model predicts how the metabolic rate and the body surface temperature depend on the ambient temperature under simplified conditions in which we can solve the model analytically. In particular, our model explains: • details of the classical U-shaped relationship between metabolic rate and ambient temperature shown in Figure 1; • season-dependent shifts in this relationship as a response to changes in the insulation layer properties (e.g. fur colour and thickness in mammals); • extension of the TNZ to high temperatures due to vasomotor adjustment; • linear and nonlinear responses of body surface temperature as a function of ambient temperatures as observed in infrared thermal assays.
Finally, we remark that our work concentrates in understanding qualitative aspects of thermoregulation. However, we believe that our model may be useful for further theoretical developments as well as for numerical applications.

The classical model revisited
According to Hill et al. (2012), the animal's metabolic rate is low and independent of T a inside the TNZ, where thermoregulation is achieved by autonomic modulation of body insulation. The TNZ is limited by its lower and upper critical temperatures T LC and T UC , respectively. Below the T LC body insulation is at its maximum value and heat production required for thermoregulation increases as T a falls. The rate of change of heat production is given by conductance, C. Therefore, for ambient temperatures below T UC the body temperature is maintained constant by dry heat exchange, which does not involve the evaporation or condensation of water. Above T UC , dry heat exchange is not enough to maintain body temperature constant and the metabolic rate increases with heat exchange through evaporation.
Mechanistic details of the animal-ambient heat exchange process are well-known from an experimental standpoint (Hill et al., 2012(Hill et al., , 2013Randall & Burggren, 2001;Rey et al., 2015). However, the theoretical model commonly used for the thermoregulation process relies heavily on arguments derived from Fourier's law of heat flow (Hill et al., 2013;Randall & Burggren, 2001), which only explains the shape of the metabolic-temperature curve for ambient temperatures below T UC .
If factors other than ambient temperature are held constant, the 'rate of dry heat transfer' is proportional to the difference T bc − T a , known as Fourier's law of heat flow (Hill et al., 2012). To explain the shape of the metabolism-temperature curve, Randall and Burggren (2001) write the Fourier's law of heat flow as As the ambient temperature drops within the TNZ, thermal insulation increases from its minimum value up to its maximum, for instance by raising fur or feathers. Therefore, thermal conductance C, given by the inverse of the animal's thermal insulation, I, is bounded between a minimum and a maximum value. Because the body temperature of a mammal or bird is typically higher than the animal's upper-critical temperature, we have that (T bc − T a ) > 0 and dry heat exchange carries heat away from the body when T a < T LC (Hill et al., 2012). To compensate for the heat loss, the organism makes heat metabolically at a rate that matches its rate of heat loss so that the body temperature can be kept constant. Equation (2) gives a straight line with negative slope −C for the relationship between Q and T a when T a < T LC . According to Randall and Burggren (2001), the thermal conductance determines the slope of the plot below the neutral zone. Therefore, low conductance implies a shallower slope and less metabolically heat production is required at low temperatures. When the ambient temperature is in the TNZ, Hill et al. (2012) suggest that a mammal or bird responds to a decrease in ambient temperature by decreasing their conductance, which increases the animal's resistance to heat loss, such that the animal's rate of heat loss and its metabolic rate remain constant.

Criticism of the classical explanation
All the previous explanations are interesting and revealing, but they deserve some criticisms as previously discussed by Kingma, Frijns, Shellen, and Van Marken Lichtenbelt (2014). First of all, a more precise approach should consider that heat exchange between an animal and the ambient depends directly on T (the body exchanges heat through its surface) and indirectly on T bc . This indirect dependence is due to heat transfer occurring inside the body, and it is at this point that the classical explanation has some implicit hypotheses that must be clarified. The main one is that the Fourier's law of heat flow is a linear law for one specific mechanism of heat transfer, that is, heat transfer by conduction. In other words, the classical explanation is based on the hypothesis that the only heat transfer mechanism inside the body is heat conduction. Figure 2. Schematic diagram of heat flow from the body core of an animal through internal tissues and coat (feathers or fur) to the ambient. The animal's body core is at temperature T bc , its surface is at temperature T , and the ambient temperature is denoted by T a . This system is described by Equation (3) in a simplified scenario that considers homogeneous media, steady-state situation, conduction as the only mechanism of heat transport inside the body, and convection as the only mechanism of heat exchange with the ambient.
Under this (over)simplifying hypothesis and assuming homogeneous conditions (constant coefficients) in steady state, and that heat is exchanged with the ambient through convection, Equation (2) can indeed be obtained. To do this based on first principles (see Figure 2), we start by observing that heat loss due to the convection mechanism is given by where T is the surface temperature and C s is the superficial convection coefficient. Assuming that conduction is the only heat transport mechanism inside the body, heat loss from the body core to the surface can be expressed by Q =C(T bc − T ), whereC is the mean body thermal conductance (see Equation (1)). Therefore, if all the heat produced in the body core is transferred to the surface and exchanged with the ambient, we can replace this expression in Equation (3) to obtain T = [C s /(C s +C)]T a + [C/(C s + C)]T bc . This expression can be replaced back in Equation (3) to obtain Fourier's law shown in Equation (2) with C = C sC /(C s +C). We conclude that only in this oversimplified scenario, Equation (2) is appropriate to describe heat exchange between the animal's body core and the ambient. However, the assumption that T is homogeneous in space is inadequate. In fact, although it is reasonable to assume that T bc is constant, surface temperature, T , is spatially heterogeneous because different body regions play different roles in the thermoregulation process. Therefore, a realistic thermoregulation model should consider different T associated with different body regions.
More important, conduction is not the only heat transport mechanism inside the body. Blood mass circulating through the body transports heat in an efficient way from the body core to the surface affecting T (Randall & Burggren, 2001). From the mathematical point of view this is a more complex mechanism than heat conduction. Therefore, it is not clear that T depends linearly on T a and T bc and, consequently, that Q can be expressed by Equation (2). Another point is that, in addition to convection, radiation is an important mechanism of heat exchange with the ambient that must be considered. Heat exchange by radiation is given by σ (T 4 − T 4 a ), where σ is the Stefan-Boltzmann constant and is the emissivity of the body surface, introducing a nonlinearity to the thermoregulation model. Moreover, heat exchange through evaporation, which is the most important mechanism of heat exchange at T a > T UC , is also not considered in the classical model. Therefore, the classical explanation (Equation (2)) for the shape of the metabolism-ambient temperature curve ( Figure 1) is too simple to be satisfactory.
We aim to develop and analyse a model for heat production and exchanges that considers the internal structure of the body, different types of internal heat transfer, different mechanisms for heat exchange with the ambient, and control mechanisms used for thermoregulation.

Definitions, hypotheses, relevant physical variables and parameters, derivation of the thermal balance equations
In this section, we describe the concepts, notations, relevant physical parameters, and variables used in our model. We also derive the relevant thermal balance equations.
We denote by a bounded domain corresponding to the stylized body of a thermoregulated animal under consideration. Different from previous works (Hill et al., 2012), which divided the animal's body into two layers (see Figure 3(a)), we divided the animal's body into four layers (see Figure 3(b)). Each layer corresponds to a group of organs or tissues with similar roles in thermoregulation (see Figure 3(b)). The animal body is also divided into J ∈ N regions represented by different colours in Figure 3(b), such that physical variables may change with time, but are assumed to be uniform across each region. Table A1 in Appendix 1 displays the relevant parameters and variables we use, and the next subsections detail the particularities of each layer. We use subscript and superscript indices to denote the layer and the region that a variable corresponds to, respectively. Subscript indices can be bc for body core, s for tissues surrounding the body core, p for perfused tissues near the surface, and c for coat. For instance, the variable T (i) p corresponds to the temperature distribution in the superficial perfused tissue of the ith region (see Figure 3).
In the next subsections we will derive the equations for variables and parameters of our model based on first principles, i.e. by balancing the thermal energy in each of the subdomains of the previously mentioned layers. These equations all have the following general, standard form (see Gates, 2012;Hill et al., 2012;Lienhard & Lienhard, 2001;Randall & Burggren, 2001): time rate of the thermal energy density = diffusion term + sum of sources of thermal energy densities − sum of sinks of thermal energy densities.
To obtain the equations with thermal energy balance, we have to identify the expressions for the terms of Equation (4) in each case. As mentioned previously, the heat transfer mechanisms considered in this work are (i) conduction inside and between layers of tissues, (ii) blood flow (mass transfer), (iii) evaporative cooling by the respiratory tract (panting), (iv) superficial evaporative cooling, (v) heat exchange by convection, and (vi) heat exchange due to radiation. All these mechanisms will be explained in detail in what follows.

Layer 1: body core
The innermost layer is the body core bc ⊂ , which includes all of the essential internal organs, such as heart, lungs, stomach, liver, brain, and so on. For proper functioning, The four-layer model of an animal's body introduced in this paper. Layer 1: the body core layer ( bc ) is formed by all of the essential internal organs, which, for proper functioning, are required to be maintained near an objective temperature (T obj ). Layer 2: The surrounding tissue layer ( s ) is the next layer and is formed by all blood perfused tissues that surrounds the body core and also take part in heat exchange. These surrounding tissues can be separated according to the layers that are on top of each of them. Layer 3: The next layer corresponds to the blood perfused tissues near the surface ( p ) that play a significant role in thermoregulation. Such tissues are on top of previous layers. Blood perfused tissues are distinguished according to the possibility of passive or active (via vasomotor adjustment) blood flux control. Layer 4: the coat layer ( c ) is formed by the subdomains protecting the previous subdomains. The coat is in direct contact with the ambient and exchanges heat with it. We assume that the animal exchanges heat with the external ambient by convection, by radiation, by evaporation via the respiratory tract and via superficial evaporation. Different colours represent different body regions (e.g. legs, head, trunk), which are denoted in our notation by the superscript index (i).
these organs are required to be maintained near a given fixed temperature (Hill et al., 2013;Randall & Burggren, 2001). For simplicity, we assume that the body core has homogeneous physical properties (we work with average values) and that its temperature may vary in time, but it is the same across the entire subdomain. For proper functioning, the body core is required to be maintained at the objective temperature T obj , which is the temperature that the thermoregulation mechanisms should attain (Hill et al., 2012;Randall & Burggren, 2001). The current body core temperature at time t is denoted by T bc (t), which may be different from T obj . We denote by m bc the known body core mass and by c bc its known (mean) specific heat. In such cases of spatial homogeneity it is more convenient to work with the integrated form of Equation (4). Then, integrating it on the body core volume, we obtain: time rate of the total thermal energy = total heat flux at the body core boundary + sum of sources of thermal energy.
− sum of sinks of thermal energy.
First, the term corresponding to the total thermal energy in the body core at time t is expressed as m bc c bc T bc (t) (see Lienhard & Lienhard, 2001, for instance). Second, due to the functioning of the organs in the body core (metabolism), there is heat production q bc per unit of volume and unit of time, which yields to a total heat contribution of q bc vol( bc ) per unit of time, where vol( bc ) is the known body core volume. The last two terms in Equation (5), are closely related to the heat transfer mechanisms. We observe that from the heat transfer mechanisms considered in this article, three of them directly contribute to the balance of thermal energy of the body core: (i) conduction to and from surrounding tissues, (ii) blood flow (mass transfer), and (iii) evaporative cooling by panting. The other two heat exchange mechanisms (superficial evaporative cooling and convection) only indirectly affect body core temperature through heat flow from their interfaces with the other layers of tissue.
Contributions due heat conduction: Thermal energy is exchanged between the body core and its surrounding tissues (i) s is its heat diffusion coefficient, ∇ denotes the gradient operator, n (i) s is the external unit normal to (i) s , and ∂T (i) s /∂n (i) s is the directional derivative of T (i) s in the normal direction. Therefore, the thermal energy flowing from s ∩¯ bc , and the total heat flux at the body core boundary is the sum of all the terms of this kind.
Contributions due blood flow: Blood leaves the body core to irrigate all the tissues in the organism and provides nutrients and oxygen (Hill et al., 2012;Randall & Burggren, 2001;Schmidt-Nielsen, 1997). In addition, blood leaves the body core at temperature T bc and carries thermal energy away. However, blood also returns from blood perfused tissues carrying some thermal energy back to the body core. We assume that blood flows from the body core to surrounding tissues regions with blood perfusion (i) s of Layer 2 and back through independent circulatory pathways. Similarly, blood flows from the body core to blood perfused regions (j) p of Layer 3 ( Figure 3). We assume that there is no direct blood flow between any surrounding tissues regions with blood perfusion (j) s of Layer 2 and any blood perfused region (j) p of Layer 3. This assumption prevents the possibility of counter current heat exchange. At the end of this section we discuss how counter current heat exchange could be incorporated in our model.

Concepts and notations
Let a generic blood perfused tissue be denoted by O, its temperature by T O . In addition, let the known physical parameters of this tissue be the mass density ρ O , the specific heat per unit of mass c O , the heat diffusion coefficient k O , and the heat production per unit volume q O .
We denote by φ O the normalized density of volumetric blood flux in O effectively used for heat exchange, according to the BHT approach (Pennes, 1948). Such φ O gives the blood vessels distribution used for heat exchange and their respective capacities. According to the mechanism used for blood flux control, the tissue region O can be classified either as passive or active blood flux region.
Passive blood flux subdomains: If O is less specialized for thermoregulation, blood vessels do not undergo significant constriction or dilation and constant volumetric blood flux is maintained (Arkin et al., 1994;Gordon, 1990;Hill et al., 2012;McCafferty et al., 2011;Savage & Brengelmann, 1996;Tattersall et al., 2009). Thus, we have: Active blood flux subdomains: If O is specialized for thermoregulation, it is provided with vasomotor adjustment mechanism and it is able to change the volumetric blood flux by contracting or dilating the blood vessels in response to the body temperature (Arkin et al., 1994;Gordon, 1990;Hill et al., 2012;McCafferty et al., 2011;Savage & Brengelmann, 1996;Tattersall et al., 2009). In this case, we have: unknown and may depend on time.
By using these previous notations, the amount of blood transported per unit of time from the body core to O is w O ρ bl when O is a passive blood flux domain and w O (t)ρ bl when O is an active blood flux domain. In these expressions ρ bl is the blood density, and w O and w O (t) are the volumetric blood flux per unit of time arriving at O, respectively, in the passive and active cases. Thus, the thermal energy per unit of time transferred from bc to O contributes with the following sink terms in Equation (5): where c bl is the blood specific heat (see Lienhard & Lienhard, 2001, for instance).
Blood also returns from blood perfused tissues, O, carrying thermal energy back to the body core. Because thermal balance is fast in the blood vessels of the perfused tissues (Arkin et al., 1994;Gordon, 1990;Kellogg, 2006), the blood returns to the body core at a temperature corresponding to the average temperature of the volumetric blood vessels distribution in O. Such averages temperatures are, respectively for active blood flux domains. Therefore, the thermal energy per unit of time reentering the body core contributes with the following source terms in Equation (5) is: The previous arguments will be applied to each of the blood perfused regions that are important for thermal regulation. Away from the body core, these regions lay basically in the the layer of surrounding tissues, where the muscles are located and in the third layer, where more superficial tissues are located. These latter tissues, by their position, facilitate heat exchange of the blood perfused in them with the ambient and thus play an important role in thermoregulation.
We denote the blood perfused regions in the layers of surrounding tissue by for active blood flux tissues. Similar arguments hold for the superficial blood perfused tissues. We denote by (i) p , i = 1, . . . , J p such regions. We consider the firstJ p (0 ≤J p ≤ J p ) of them to be passive blood flux domains and the last ones, active blood flux domains. We will also use such previous expressions when considering the balance of thermal energy in Layers 2 and 3. In blood perfused tissues near the surface (Layer 3), for instance, we will replace O by for active blood flux tissues. Similar arguments will be made for blood perfused tissues surrounding the body core (Layer 2).
Contributions of evaporative cooling from the respiratory tract (panting): By panting, water vapour leaves the body at body core carrying thermal energy away (Richards, 1970). We assume that the amount of heat loss depends on the respiratory ventilation rate, vr, which is the average volume of air per unit of time expelled by the animal. This ventilation rate is proportional to the respiratory frequency and depth of breathing (Richards, 1970). We are interested in the net ventilation rate that assumes values in the interval [v min , v max ], where v min is the basal ventilation rate at rest and v max is the maximum ventilation rate that depends on biological constraints of the organism. Notice that our phenomenological approach does not describe the biological details occurring in the respiratory tract during panting. Please see Fiala, Lomas, and Stohrer (1999) for more details. We assume the mass of water vapour expelled per unit of time in ventilation is written as where T bc is the body core temperature and T a is the external ambient temperature. The parameter α 0 (T bc , T a ) is the mass of water vapour expelled at basal ventilation rate, v min , andα(T bc , T a ) is a positive parameter. The increase of ventilation rate from v min to v is associated with an additional heat production per unit of volume, q v , due to extra work done by muscles of the respiratory tract. This additional heat production depends on the excess ventilation, v = vr − v min , which we assume to be given by whereq v (T bc , T a ) is a parameter that may depend on the body core and the ambient temperatures. Thus, the total energy per unit of time generated in the body core is (q bc + q v ) vol( bc ). Finally, the thermal energy per unit of time lost to the ambient due to the ventilation mechanism is the product of the mass of water expelled and the L is the latent heat of water vaporization per unit of mass (see Lienhard & Lienhard, 2001, for instance). Thus, This contributes with a sink term −α v L in Equation (5).

Balance of thermal energy in the body core
By using the previous notations and arguments in Equation (5) and rearranging the terms, we can write the following integro-differential equation for the balance of the thermal energy in the body core: The first term in the right-hand side corresponds to the total heat production in the body core per unit of time. The second term is the heat lost due to the evaporation by the respiratory tract. The third and fourth terms correspond to heat leaving the body core (at the body core temperature) due to passive and active blood flux per unit of time, respectively. The fifth and sixth terms are heat returning to the body core due to passive and active blood flux per unit of time, respectively, after temperature equalization due to heat exchange in the blood perfused tissues. The last term is heat exchanged per unit of time through thermal conduction among the body core and its surrounding tissues.

Layer 2: tissues surrounding the body core
The next innermost layer consists of muscles, bones and other tissues that surround the body core and that also take part in heat exchange in addition to other specific functions they may have. These surrounding tissues are denoted by (i) s ⊂ , i = 1, . . . , J s , and their temperature by T (i) s . The known physical variables of this layer are the mass density ρ (i) s , the volume, vol( s , according to the BHT approach (Pennes, 1948). According to the mechanism used for blood flux control we propose that each subdomains (i) s can be classified as passive or active blood flux subdomains.
Because the properties and temperatures in this layer are assumed not to be spatially uniform, we use the local form of the balance of thermal energy given by Equation (4) to describe the thermoregulation process in this layer (see Lienhard & Lienhard, 2001). The thermal energy density is now expressed as ρ (i) s c (i) s T (i) s . Next, the diffusion term is given by the divergence of the heat flux by conduction, that is, div(k As for the heat source terms densities, one of them is due to metabolism, q (i) s . Another source term is due to thermal energy carried by the blood flow coming from the body core, w s ρ bl c bl T bc , and then a sink term due to the loss of thermal energy carried by the blood flow returning to the body core, −w (i)

Balance of thermal energy in the tissues surrounding the body core
By using the expressions obtained in the last subsection in Equation (5), rearranging the terms, the balance of thermal energy per unit volume in (i) s is written as: where . . , J s . The first and the second terms on right-hand side correspond to heat conduction and the heat generated in the tissue per unit of time, respectively.
Notice that the previous treatments of the contribution of the blood flow to the heating are similar to BHT (Pennes, 1948). The distribution φ (i) s of blood vessels in (i) s appears as a product factor in the source-sink term of the original Pennes' model.

Layer 3: blood perfused tissues near the surface
The next layer in our model consists of blood perfused tissues that play a significant role in thermoregulation. Lying on top of each of the subdomains (i) p , i = 1, . . . , J, we have a corresponding (i) p ⊂ , i = 1, . . . , J, subdomain that corresponds to blood perfused tissue. Their temperature is denoted by T (i) p and the known physical variables of this layer are the mass density ρ (i) p , the specific heat per unit of mass c p , and the heat production per unit volume q (i) p . Again we use φ (i) p to describe the normalized density of volumetric blood flux in (i) p , according to the BHT approach (Pennes, 1948).

Balance of thermal energy in the blood perfused tissues near the surface
By proceeding exactly as we did for the regions of the previous layer, the balance of thermal energy per unit volume in (i) p is written as: where φ The first term in the previous right-hand side represents heat transport in the tissue due to diffusion by conduction. As before, the second term is the source-sink term associated with heat exchange among the tissue and the blood in the perfusing vessels. The third term is the heat production rate in the tissue.

Remark:
The previous two layers look quite similar since both are assumed to have similar nets of perfused blood vessels, and so one could think that they could be considered just one layer. However, we stress that these layers play different roles in thermoregulation, and it is important to consider them individually. In fact, blood perfused tissues near the surface (Layer 3) are always engaged in thermoregulation, even when the animal is at rest because near the surface they are more directly enrolled in heat exchange with the external ambient. The same is not true for the tissues surrounding the body core (Layer 2). These tissues display an important role in thermoregulation specially when the muscles are active during exercise. These situations require large increases of their metabolic rates q (i) s and lead to an increase of muscle temperatures T (i) s . The exceeding muscle temperature requires an increase of blood flow in the muscles to take away the heat excess and prevent muscle from damage.

Layer 4: coat
The outermost layer in our model is the coat subdomains denoted (i) c ⊂ , with i = 1, . . . , J c . The coat is in direct contact with the ambient and exchanges heat with it by convection and radiation (Hill et al., 2012;Randall & Burggren, 2001). We assume that inside the coat subdomains heat flows only by conduction and heat production inside the coat is negligible. An important aspect of this layer is that it can be used for active thermoregulation in certain animals, which is done by the raising of hairs or feathers (Hill et al., 2012;Randall & Burggren, 2001), for instance. We analyse the effects of this mechanism in the model by considering changes in coat width. The temperature of (i) The known physical variables of the coat are mass density ρ (i) c , the specific heat per unit of mass c c , the heat production per unit volume q (i) c (assumed to be zero), the convection coefficient d (i) > 0 at the interface with the ambient¯ (i) c ∩¯ a , and the emissivity (i) of the interface¯

Balance of thermal energy in the coat
The balance of thermal energy per unit of volume of the ith subdomain (i) where the term in the right-hand side represents heat transport in the tissue due to diffusion by conduction in the coat.

Remark:
Counter current heat exchange is related to the existence of a blood flow connection between the blood perfusion system of a region (j) s of Layer 2 and the blood perfusion system of the associated region (j) p of Layer 3. For instance, suppose we have a situation where blood flows from the body core to (1) s , then to (1) p and only then returns to the body core. Imagine no circulatory pathway from the body core to (1) p . This situation could be easily included in our model just by taking out or including appropriated terms of the equations that describe the balance of energy in the body core. In Equation (9), the indexes i in the first summation of the fifth line and the summation in the sixth line would start with i = 2 instead of 1. Equation (10) would be the same, and Equation (11) We will analyse this type of situation in future works.

Interface conditions
The previous balance equations holding for different body regions must be supplemented by conditions at their interfaces. In fact, heat is also exchanged through the interfaces between subdomains that are in contact with each other, and temperature and the thermal fluxes must be continuous at these interfaces. For instance, the interface condition associated to the interface¯ (i) s ∩¯ (i) p is such that

Heat exchange with the ambient and boundary conditions
The ambient is denoted by a = R 3 −¯ . The relevant physical parameters related to heat exchange between the animal and the ambient is the known air temperature T a and also its relative humidity r a (Gates, 2012;Hill et al., 2012;Randall & Burggren, 2001). Heat is exchanged with the ambient through heat fluxes due to convection, f conv , evaporation from the surface, f evap and radiation, f rad , 3 occurring at the body's boundary (Gates, 2012;Hill et al., 2012;Randall & Burggren, 2001), which is composed of the interfaces between the subdomains c i and the ambient. Thus, the boundary conditions of the problem are obtained by requiring that the total heat flux is the sum of the previous particular forms of the fluxes: To complete these conditions, following Cooney (1976), we next briefly describe each of these heat fluxes at the right-hand side.
Convection: a rather good approximation for the thermal convective flux at a part¯ c i ∩ a of the boundary between the body and the external ambient is given by the expression where d (i) is the convection coefficient associated with that part of the boundary, T (i) c is the temperature at the boundary surface, and T a the ambient temperature near the boundary.
Superficial evaporation: the details of the heat flux due to superficial evaporation are complex. Therefore, we propose a phenomenological approach and use a rather good approximation for the range of temperatures with which we are concerned: where T (i) c and T a are as before; r a is the relative humidity of the air; S r is the sweating rate per unit of area of the surface; L is the latent heat of water vaporization per unit of mass at ambient pressure (generally assumed as 1 atm); P(·) is the vapour pressure of water as a function of temperature in the condition of ambient pressure. Observe that r a P v (T a ) is the partial pressure of water vapour in the ambient air. P(·) decays in a complicated way with the increase of the temperature, but within the range we are considering, it may be approximated by an linear expression P(T) = g 1 T − g 2 for suitable positive coefficients g 1 and g 2 .
Moreover, the coefficient κ (i) evap may depend on several aspects such as the coat evaporative resistance and wind speed. Details of these dependences and values for the parameters involved in the present arguments can be found, for instance, in Cooney (1976), Yildirim (2005) and also Parsons (2013).
We also remark that the values of the sweating rate are restricted to the interval 0 ≤ S min r ≤ S r ≤ S max r < +∞. Here, S r = S min r corresponds to the situation of passive minimal superficial evaporation, without active sweating (latent heat loss), whereas S r > S min r corresponds to active sweating situations. Thus, both passive and active situations are included in the model. Radiation: the thermal flux due to radiation is given by where (i) is the emissivity coefficients (0 ≤ (i) ≤ 1) of the surface that are exchanging heat with the ambient, and σ is the Stefan-Boltzmann constant, 4 and the other symbols meanings are the same as before. Thus, the boundary condition associated with the part¯ c i ∩¯ a becomes: where d (i) 's and (i) 's are, respectively, the convection and emissivity coefficients (0 ≤ (i) ≤ 1) of the surfaces that are exchanging heat with the ambient, and σ is the Stefan-Boltzmann constant. 5

Remark:
The parameters in Equation (13) are assumed to be fixed known constants, as we are not considering postural changes or wind speed variation. However, they could be incorporated at the expense of having a more complex model, by using expressions known in specific cases. For instance, the convection coefficient d (i) depends on the position at the interface between the coat and the ambient with respect to its geometry and angle, and whether the convection is free or forced, which may be a function of wind speed (Gates, 2012).

Thermal feedback controls
The neurophysiological details of the control mechanisms involved in thermoregulation in endotherms are rather complex, and they are still important subjects of scientific investigation (Jessen, 2001). However, in order not to complicate too much the present model, we do not include the neurophysiological bases of the specific actions of the controls to be considered and instead take again a phenomenological approach. We exemplify this approach in the present context by describing how we consider the process of change in the thermal insulation in an specific part of the coat, say c i , by changing some of its specific physical characteristics, whose value at a time t is denoted by I i (t). For instance, under certain temperature limits, fast thermal insulation control can be obtained by the raising or the lowering of feathers or hairs. In this case, we could consider I i (t) as the width of the resulting coat at time t.
To obtain an equation for changes in I i (t), we observe that for a difference of the actual body temperature and the required objective temperature, the body must respond by increasing or decreasing the thermal insulation at a certain rate. This rate must be proportional to R c i (T bc (t) − T obj ), meaning that it must be a function of the difference between the actual body temperature and the objective temperature and its sign determines whether there will be an increase or a decrease of the thermal insulation.
Moreover, for biological reasons, there must be a minimum and a maximum value for I i (t). That is, through time the values of I i (t) must be in an interval [I min are, respectively, known minimum and maximum values for the characteristic variable that determines the thermal insulation. This is guaranteed if we assume the following feedback control equation for I i : We stress that for initial condition I min , Equation (14) forces I min for t > 0. Moreover, in this case, the thermal insulation decreases for R c i > 0 (which must be the case when T bc > T obj ) and increases for R c i < 0 (which must be the case when T bc < T obj ).
Thus, the dependence of such rate R c i on T bc − T obj must imply an automatic feedback control, which we describe as a proportional-integrated-differential (PID) controller (Bennett, 1996), with certain non-negative parameters (λ c i,P , λ c i,I , λ c i,D ) as follows: • As for the proportional part of the controller, when the current body core temperature is above the objective temperature it acts by decreasing the coat thermal insulation to loose heat to the ambient. On the contrary, when the body core temperature is below the objective temperature the proportional part of the controller increases the insulation to lower the heat lost to the ambient. i,f (τ ) = 0 for τ >τ .
• The differential part of the control acts by hampering sudden changes and oscillations in the thermal insulation. Therefore, The same arguments apply to the other control mechanisms we consider in the present work. The only possible difference is the − sign in the right-hand side of Equation (14), which, depending on physical/biological meaning of the control under consideration, should be changed to +. That would be, for instance, a situation for the vasomotor adjustment in the perfused tissues near the surface because it is necessary an increase of blood flow near surface to facilitate heat exchange when T bc (t) > T obj , in contrast with coat insulation which decreases under the same condition.
In the next subsection we summarize this approach in general terms and show the differences for each kind of control in Table 1.

General approach
We consider six thermal control mechanisms in animals: • Changes in the metabolic rate.
• Changes in the thermal insulation due to changes in properties of the coat.
• Changes in the blood fluxes due to vasomotor adjustment in the tissues surrounding the body core. • Changes in the blood fluxes due to vasomotor adjustment in the tissues near the surface.
• Changes in evaporation from the respiratory tract due to changes in ventilation rate.
• Changes in evaporation from the animal's surface due to changes in the sweating rate.
For each mechanism we assume that there is a variable z(t) that has to be controlled such that its value lies inside the interval z(t) ∈ [z min , z max ] at any given time t. To model and control the responses of z(t) as a function of changes in T bc we use a PID controller with non-negative parameters λ z P , λ z I and λ z D (Bennett, 1996). Such PID controller can be expressed as where κ = −1 when the thermal control mechanism being considered is based on vasomotor adjustment, and κ = 1 otherwise. 6 The proportional part of the controller guarantees that the feedback control acts by decreasing z(t) when T bc > T obj or increasing z(t) otherwise. The integral part of the controller acts similarly as the proportional part but with respect to the past values of T bc − T obj . The past values are modulated by a given positive forgetting kernel g z (·), which must be an integrable function on [0, +∞] and lim τ →+∞ g z (τ ) = 0. 7 Finally, the differential part of the control acts by hampering sudden changes and oscillations in z(t).

Interpretation of z(t) for different thermal control mechanisms
Minimum and maximum values of z(t) are usually related to physical or energetic constraints of the system. For instance, the minimum value of the metabolic rate corresponds to the value required to keep all the animal's essential functions working properly. On the other hand, the maximum value of metabolic rate corresponds to the maximum energy the animal is able to produce.
In the case of changes in the properties of the coat, the minimum and maximum values of thermal insulation may depend specifically of what property is changing. The thermal insulation can be modified at a fast time scale in certain temperature ranges by raising or lowering of feathers or hairs, which ultimately changes the effective width of the coat. In this scenario, minimum and maximum values of thermal insulation would represent the situation where the hair is totally flat and erect, respectively. Season-dependent coat properties may also change at a much slower time scale. For instance, the thickening of the fur during winter lowers the convection constant d (i) , and alterations in fur colour changes the associated coat emissivity (i) . Notice that both these variables have limits due to physical constraints. For the mechanism based on vasomotor adjustment, the minimum and maximum values of blood fluxes are a direct consequence of the minimum and maximum that vessels can be constricted and dilated, respectively. Minimum and maximum values of evaporation from respiratory tract correspond to the minimum and maximum frequencies the animal can breathe per unit of time, while minimum and maximum values of superficial evaporation correspond to the to minimum (passive evaporation) and maximum of the sweating rate. Table 1 shows the correspondence between z(t) and real physical variables for each thermal control mechanism.

Endurable temperatures and steady-state equations
We are interested in steady-state solutions of our model that correspond to a situation where thermal equilibrium between the animal and the ambient is achieved. To define these solutions more precisely we need to introduce definitions of endurable temperature and set of endurable temperatures.
Definition: endurable temperature is the ambient temperature T a such that when an organism is placed at an ambient at this temperature for long enough to reach the thermal equilibrium with the ambient, T bc remains at T obj .
Definition: set of endurable temperatures is the set of all endurable ambient temperatures T a . Therefore, our model admits a steady-state solution (T bc = T obj ) when the ambient is at an endurable temperature. From the physiological point of view, it means that the thermal control mechanisms incorporated into the model are capable of sustaining T bc = T obj as long as the ambient is maintained at temperature T a . The set of endurable temperatures is usually expected to be an interval and we then call it interval of endurable temperatures.
To determine the existence of endurable temperatures, we must analyse the existence of solutions of our model that do not depend on time and with T bc = T obj . This is done in the next section where we write the previous thermal balance equations as well as the feedback control equations with no time dependency.

Steady-state thermal balance equations
The steady-state version of Equation (9) that gives the thermal balance in the body core is where w The steady-state version of Equation (11) that gives the thermal balance in the blood perfused tissues is where, φ (i) = φ (i) (x) for i = 1, . . . ,J, and φ (i) = φ (i) (x, w (i) ) for i =J + 1, . . . , J. Finally, the steady-state version of Equation (12) that gives the thermal balance in the coat for each is

Steady-state equations for the feedback controls
The steady-state version of the PID controller expressed in Equation (16) can be written as Because we are looking for solutions where T bc = T obj , we have from Equation (22) that R c i ≡ 0 and Equation (21) is satisfied. In addition, as the system evolves towards the steady state according to Equation (16), we have that, at steady state, the controlled variable is always z min ≤ z ≤ z max for any initial condition inside this interval. These results can be summarized for each thermal control mechanism as follows: • The metabolic rate q in the body core can assume any value such that q min ≤ q ≤ q max given that T bc = T obj . • The thermal insulation I i associated with the coat properties of domain (i) c can assume any value such that I (i) min ≤ I (i) ≤ I (i) max given that T bc = T obj .
• The volumetric blood flux w (i) in the perfused tissue (i) p can assume any value such that w (i) min ≤ w (i) ≤ w (i) max given that T bc = T obj .
• The excess of ventilation rate v can assume any value such that v min ≤ v ≤ v max given that T bc = T obj .

Global steady-state balance
Integrating Equations (18)-(20), using the divergence theorem and the boundary conditions and plugging the result into Equation (17), we obtain The left side of Equation (23) represents the total dry heat exchange which does not involve evaporation of water. On the right side we have the total internal heat production due to metabolism, Q, given by minus the total thermal energy per unit of time lost to the ambient due to evaporative cooling. Equation (23) of energy balance shows that, at steady state, the total internal heat production must counterbalance the total heat exchange with the ambient. Most of the previous works on thermoregulation take this global balance for granted and rely on it to model experimental data. Instead, we show that the global energy balance is an emerging property of our model and indicates energetic consistency.

Steady-state solutions
In this section we present the analytical expressions for steady-state solutions of our fourlayer model in the case of a simplified system. Next, we will use such expressions to advance more complete explanations for some of the features of Figure 1.
The simplified system to be considered is composed of only one region. Although our general model may include the simultaneous action of all previously described six thermal control mechanisms, to simplify the derivation of the expressions, we will remove here the thermal controls due to vasomotor adjustments in Layer 2 and also due to superficial evaporation. The reason for this simplification is that the experimental measurements leading to Figure 1 are usually obtained with animals at rest, and vasomotor adjustments in Layer 2 and superficial evaporation are expected to be more engaged when the animal is in active situations, as discussed in Section 3.3.
We divide our analysis into two parts. First, we consider the case where only passive blood flux is present. We study the effects of the three remaining thermal control mechanisms individually. In the second case we incorporate vasomotor adjustment and we discuss the consequences. As we shall see, our findings for this interesting yet simplified system allow us to discuss and explain many qualitative aspects of thermoregulation phenomena.

The model
Our simplified system is composed of one region and layers stacked on top of each other as shown in Figure 4. We assume the system is in thermal equilibrium with the ambient so that we can use the steady-state formalism derived in Section 5. All known parameters are assumed to be constant. Notice that, because we have only one region, we omit the upper-script index i that runs over different regions from now on.
The subdomains of our system are expressed as bc =˜ × (−h bc , 0); where˜ ⊂ R 2 is given, as well as the positive constants used to specify the limits of each layer, h bc , h s , h p , and h c . For the sake of simplicity, we also assume that the heat diffusion coefficient is constant and equals to k in s , p , and c . Another assumption is that the metabolic rate per unit volume q is the same in bc , s , and p . Therefore, the total internal heat production due to metabolism given by Equation (24) can be written in this case as Finally, we consider only fast time scale changes in the animal's thermal insulation due to changes in the width of the coat layer. Therefore, in Sections 6.2 and 6.3 the thermal control variable I is equivalent to h c . In Section 6.3.2 we discuss the effect of slow time scale changes in the animal's thermal insulation and how they affect thermoregulation.

Passive blood flux
To recover the equations for the case where no vasomotor adjustment is present, we set J = J = 1 in the equations of Section 5. The unknowns of our problem are the temperatures at each layer. More precisely, we want to find the value of T bc as a function of T a , and the functions T s (z), T p (z) and T c (z). Under the assumptions described above, we can write the equations derived in Section 5 as together with the following conditions at the internal interfaces between layers and the boundary condition at the interface with the external ambient (35) First, we reduce Equation (23) to our simplified case as where a = σ/d, the surface temperature is defined as and Expression (36) is similar to Equation (3). However, in Equation (36) the thermal control mechanisms are explicitly shown and include a correction term that comes from the radiative heat exchange incorporated into our four-layer model. It is worth mentioning that the emissivity of biological tissues is known to be between 0.95 and 0.98 (Gates, 2012;Monteith & Unsworth, 2008), whereas the heat convection coefficient depends on the properties of the fluid, type of flow and the object shape. Values of the coefficient have been measured and tabulated for most common fluids and object shapes for both forced and free convection (Gates, 2012;Monteith & Unsworth, 2008). For instance, the heat convection coefficient of air lies in the interval [10 −3 , 10 4 ] W m −2 K −1 . Therefore, the constant a = σ/d is of order 10 −5 at most, which does not make the correction term negligible and reflects the importance of the terms introduced due to radiative heat exchange.
Next, we use Taylor expansion up to 2nd order in Equation (36) to find Next, by evaluating Equation (36) Finally, we combine Equations (26) and (27) to obtain an expression for T s (z) in terms of T p and T bc . Plugging the result into Equation (28) and using the boundary conditions (31) and (32), we find an expression for T p and T s in terms of known parameters and T bc . We use this result together with Equations (26), (33), and (34) to obtain an expression for T c in terms of T bc that can be written as where the coefficients G i , i = 1, . . . 8, depend only on T obj and the known parameters ρ bl , c bl , h bc , h s , h p , k and w. Please, refer to Appendix 2 for more details. Neglecting the quadratic terms in Equation (41) and substituting the expression for T in (42), we have As T a changes, the three thermal control mechanisms must change to satisfy Equation (43). To gain insights about the thermoregulation phenomena, we analyse this expression in the next section when thermal control mechanisms change individually.

Changes in thermal insulation of the coat
Increasing the metabolic rate is biologically expensive (Greenberg et al., 2012;Tattersall, 2016). Therefore, it is expected that the animal would keep, for as long as possible, the metabolic rate at its minimum, q min . Extra ventilation rate that would imply an additional metabolic rate should also be kept as v = 0 as long as possible. In fact, there are some ranges of T a that the resting animal is capable of maintaining a constant T bc with q = q min and v = 0. This is achieved by raising or lowering of feathers or hairs, which ultimately changes the coat width, h c , and therefore its insulation. Thus, assuming a constant q = q min and v = 0, Equation (43) gives the following relation between T a and h c where where M The interval [T TNZ 1 , T TNZ 2 ] corresponds to the TNZ, that is, the interval of comfort temperature attained with metabolic rate at its minimal value (Hill et al., 2012). This explains the presence of the flat region in the experimental data of the average internal heat production as functions of the ambient temperature (see Figure 5).

Changes in metabolic rate
Because coat width is biologically limited to the maximum value h max c , the mechanism of thermal regulation by changing the animal's insulation fails in maintaining T bc = T obj for ambient temperatures below T TNZ 1 . Therefore, to guarantee the existence of a steady-state solution at T a < T TNZ 1 , the metabolic rate q has to change as a function of T a according to Equation (43) with h c = h max c and v = 0 as From the definition of T TNZ 1 , we have Figure 5. Classic metabolism-ambient temperature curve according to our findings. (a) In the absence of vasomotor adjustment, we conclude that in the TNZ, represented by the B) interval, an endotherm maintains its body core temperature by changing its thermal insulation (see Equation (44)), while metabolic rate is kept at its minimum. Changes in thermal insulation are the result of changes in the width of the coat h c . At low temperatures, in the A) interval, the animal's thermal insulation is at its maximum value and the metabolic rate is required to increase to maintain the body core temperature according to the linear function expressed in Equation (52). At high temperatures, in the D) interval, the animal's thermal insulation is at its minimum and the mechanisms of evaporative cooling take place. As the ambient temperature rises, the mechanism of evaporative cooling requires an increase in the metabolic rate according to the linear function in Equation (57). (b) In the presence of vasomotor adjustment and with the animal's metabolic rate at its minimum value, we observed that the TNZ can either be expanded to the right or to the left, as shown by the C) interval, depending on the biological and physical parameters (see discussion at the end of Section 6.3).
Subtracting Equation (48) from Equation (47) we find which can be expressed as where with A > 0 for biological reasons. The metabolic rate given by Equation (25) can finally be written as where we use the fact that q v = 0 when v = 0, and V represents the animal's body volume.
Equation (52) indicates that Q increases linearly as the ambient temperature T a decreases below T TNZ 1 , as illustrated in Figure 5(a). In addition, as the ambient temperature decreases, the animal is capable of maintaining the required body core temperature by increasing the metabolic rate up to its maximum value Q max , which is defined by biological constraints. Figure 6. Schematic graph of metabolism-ambient temperature curve during summer and winter. Some animals change the thickness and brightness of their fur in a slow time scale seasonal-dependent process. Because thicker fur corresponds to larger values of h max c , our model suggests that the thicker the fur, the lower the value of T TNZ 1 according to Equation (46), which ultimately represents a shift of the dashed curve with respect to solid one. According to Equation (45), the same shift should be observed when the coat emissivity coefficient is lowered, which occurs when the fur becomes lighter. Our model also suggests that as either h max c increases or decreases, the slope of Q × T a , given by A in Equation (52), decreases.
The ambient temperature for which Q max is achieved is defined as and it represents a critical low temperature below which other behavioural mechanisms, like hibernation, would have to be used to avoid death. We highlight a fundamental difference between our model and the classical model used in the literature. The classical model serves as a mere expression for data fitting, which is of no much help in explaining the underlying physical and biological reasons behind certain behaviours, such as season-dependent slope in Figure 6. Just like the classical model, we also obtained the correct linear dependence between Q and T a given by Equation (52). However, this expression takes into consideration all the important physical mechanisms and biological parameters of the problem through the equation's coefficients. As we will se in Section 6.4, Equation (52) can be used to predict the behaviour depicted in Figure 6.

Changes in ventilation rate
At high ambient temperatures the coat width is at its minimum value h min c and the mechanism of thermal regulation by changing the animal's insulation fails. Under these circumstances, the animal has to increase the heat loss by evaporation from the respiratory tract to maintain T bc = T obj . This is done by increasing the ventilation rate through rapid shallow breathing, as seen typically in dogs and domestic fowl, and known as panting (Richards, 1970). To guarantee the existence of a steady-state solution, the excess of ventilation rate v changes as a function of T a according to Equation (43) with h c = h min c and q = q min as From the definition of T TNZ 2 , we have Subtracting Equation (55) from Equation (54) we find The last expression can be rewritten as where The previous expression holds for T TNZ 2 < T a ≤ T + , where Using Equations (7) and (57) in the definition of Q given by Equation (25), we find the relationship between the metabolic rate and the ambient temperature as We conclude that because B > 0 and with T TNZ 2 < T a ≤ T + , there is a steady-state solution in which T bc = T obj . Under these circumstances, Equation (60) indicates that Q changes linearly with positive slope as T a increases, as illustrated in Figure 5(a).

Active blood flux
To understand the role of the active blood flux in thermoregulation, we consider the same system illustrated in Figure 4, with absence of thermal control by changes in metabolic rate and animal's insulation. Therefore, the parameters corresponding to these thermal control mechanisms are assumed to be fixed. To recover the equations for the case where blood flux is controlled by vasomotor adjustment, we setJ = 0 and J = 1 in the equations of Section 5. In addition to the unknown temperatures of each layer, we also have now the volumetric blood flux w as an unknown variable. We assume that changes in w lead to changes in the superficial area˜ which is actively exchanging heat with the ambient. As a consequence, depends on w such that˜ =˜ (w).
If the blood flux uniformly fills out p , the blood volume in the perfused tissue must be equal to vol( p ) = area(˜ )(h p − h s ). Therefore, it should be true that Equations (26)-(35) derived in Section 6.2 are also valid under the assumptions described above. The reduced version of Equation (23) to our simple case with active blood flux is also given by Equation (36) but, because we have now a variable effective area for heat exchange according to Equation (61) and because we know vol( bc ) = area(˜ )h bc , we can rewrite Equation (36) as where Equation (62) can be analysed as it was done in Section 6.2 to obtain an expression similar to Equation (43) that includes now the unknown variable w and relates the ambient temperature and the different thermal control mechanisms. The resulting expression can now be analysed as thermal control mechanisms of changes in metabolic rate, insulation, and evaporation are considered individually. By doing so, we found a steady-state solution for ambient temperature below T TNZ 1 when we set h c = h max c , v = 0 and w = w min , and we let the metabolic rate change in the interval [q min , q max ]. We also found that there is a steady-state solution when q = q min , v = 0, w = w min and we let h c to change in the interval [h min c , h max c ]. This solution indicates that the animal is able to maintain the body core temperature constant when the ambient temperature changes in the interval [T TNZ 1 , T TNZ 2 ]. Finally, we found that above T TNZ 2 there is an steady-state solution that corresponds to changes in the excess ventilation rate in the interval [0, v max − v min ] while the other thermal control mechanisms are kept constant at h c = h min c , q = q min , w = w max . In summary, as T a drops and heat loss to the ambient increases, more heat needs to be produced by the animal's body so that the equilibrium between production and exchange is maintained. We have mathematically shown that the thermal control mechanism of increasing the metabolic rate linearly when T a < T TNZ 1 is due to the increase in heat production required to maintain T bc constant at these ranges of ambient temperature (Hill et al., 2012;Tattersall, 2016). The mechanism of changing the thermal insulation works for body temperature regulation within the TNZ (Hill et al., 2012;Kellogg, 2006;Savage & Brengelmann, 1996;Tattersall, 2016;Tattersall et al., 2009). Conversely, the mechanism of evaporative cooling acts when T a > T TNZ 2 (Dawson, 1982;Hill et al., 2012;McCafferty et al., 2011;Schmidt-Nielsen, 1997;Tattersall, 2016). Together, our findings explain the shape of the metabolism-ambient temperature U-shaped curve for the full range of ambient temperatures, as illustrated in Figure 5(a). We focus now on the thermal control mechanism of vasomotor adjustment. As we shall see in the next section, the action of this mechanism may extend the animal's TNZ beyond T TNZ 2 .

Changes in active blood flux
Similar computations to those performed to obtain Equation (42) can be done in the case of active blood flux. Assuming the case where q = q min and v = 0, we obtain the expression where the coefficientsG i , i = 1, . . . , 4 depend only on T obj and the known parameters ρ bl and c bl , h bc , h s , h p , and k. The dependence of T c (z) on w has been written explicitly for convenience. Using Equation (37), the surface temperature for the case of active blood flux can be expressed as Plugging this Equation (65) in Equation (62) we obtain an expression relating w and T a that can be approximately inverted and written as where the coefficients H 1 and H 2 depend on the same known parameters asG i plus the width of the coat layer, h c . As the ambient temperature changes, there must be a corresponding change in blood flux to satisfy Equation (66). Because we have H 2 > 0 for biological reasons, as the ambient temperature changes within the interval [T * 1 , T * 2 ], the active blood flux changes within the interval [w min , w max ] such that T * 1= H 1 − H 2 /w min and T * 2= H 1 − H 2 /w max . Therefore, for T a within the interval [T * 1 , T * 2 ], there is a correspondent steady-state solution with T bc = T obj . In addition, we conclude that the interval [T * 1 , T * 2 ] is contained within the animal's TNZ because the metabolic rate is at its minimum value.
The steady-state solutions of our four-layer model indicate that at least two thermal control mechanisms act in the animal's TNZ. Both changes in animal's thermal insulation and vasomotor adjustment can be used to keep the body core temperature constant in the TNZ. This can be achieved by having these two mechanisms totally overlapping each other in the same range of ambient temperature. However, we think that a more biologically reasonable possibility is that these mechanisms complement each other, such that the vasomotor adjustment is triggered in either region where the thermal insulation is close to its minimum or its maximum value.
Active blood flux at maximum insulation: if the biological and physical parameters are , then the TNZ is extended towards high ambient temperatures, as illustrated in Figure 5(b). This should be the case of animals that live in regions where T a is high, such as toucans or desert mammals (Greenberg et al., 2012;Tattersall, 2016). We remark that in such cases, the corresponding expression for the excess evaporation rate for ambient temperatures larger than T TNZ 3 is given by Active blood flux at minimum insulation: on the other hand, if the biological and physical parameters are such that T * 1 ≤ T TNZ 1 < T * 2 < T TNZ 2 , then the TNZ is extended towards low ambient temperatures. This should be the case of animals that live in regions where T a is high, such as toucans or desert mammals (Greenberg et al., 2012;Tattersall, 2016). This could be the case of certain aquatic mammals, like whales, which use shunt mechanisms in their vascular system for thermoregulation (Heyning, 2001;Parry, 1949).

Size and season-dependent changes
In addition to the fast time scale changes in the animal's thermal insulation caused for instance by raising fur or feathers, it is also known that slow time scale season-dependent mechanisms may affect thermal insulation. For instance, some animals tend to develop a thicker and/or lighter fur in winter (Hill et al., 2012) in response to the drop in T a .
To analyse the consequences of these slow time scale changes, we observe that thicker fur corresponds to larger values for h max c in Sec. 6.2. From the expression of T TNZ 1 given by Equation (46) On the other hand, the slope of the metabolism-ambient temperature curve when T a < T TNZ 1 , given by Equation (51), decreases as either h max c increases or decreases. Therefore, our four-layer model explains the seasonal changes in both the slope of the metabolism-ambient temperature curve and in the shift of T TNZ 1 to the left as illustrated in Figure 6. These changes are in agreement with experimental data of animals monitored in both winter and summer (Hart, 1957;Nilssen, Sundsfjord, & Blix, 1984;Scholander et al., 1950).

Changes in body size
Equation (51) implies that A decreases as h p + h bc increases and all the other parameters are held fixed (including the ambient temperature T a ). On the other hand, when T a ≤ T TNZ 1 , Equation (50) can be rewritten as q= A(T TNZ 1 − T a ) + q min , meaning that q decreases as h p + h bc increases and all the other parameters are held fixed. Since an increase h p + h bc corresponds to an increase in the animal's size (volume or weight), this argument leads to the conclusion that an increase in the animal size corresponds to a decrease in the metabolic rate when T a ≤ T TNZ 1 , which is in agreement to the experimental data (Heldmaier et al., 2004;Scholander et al., 1950).

The relationship between surface and ambient temperature
Distinct regions of the body of an endotherm are known to be at different temperatures at the same time for a given ambient temperature (Angilletta, Brandon, Matthew, & Justin, 2010;McCafferty et al., 2011). We use heuristic and approximation arguments together with previous reasoning to obtain an approximate expression for the surface temperature of different regions in terms of T a . We assume that in the short run each region equilibrates its temperature independently and is responsible for exchanging a fixed fraction f (i) of dry heat, i.e. total internal heat produced minus heat lost by evaporative cooling, such that The coefficients of convection d (i) and emissivity (i) are assumed to be constants. The surface temperature at¯ (i) c ∩¯ a is also constant and denoted by T (i) . Therefore, Equation (23) can be written for each region as . (69) We analyse the equation above as it was done to obtain Equation (41) from Equation (36). After some computation we find By taking into account the thermal control mechanisms studied in the previous sections and how they shape the metabolism-ambient temperature curve in absence of vasomotor adjustment, we can write the metabolic rate as a function of T a as where we assume that the metabolic rate cannot increase above Q max when T a < T − , and the positive constants and C 1 and C 2 are given by By using Equation (71) and Equation (8) in Equation (70), and recalling that the excess ventilation rate is v = 0 on [T − , T TNZ 2 ] and Equation (57), on [T TNZ 2 , T + ], we conclude that, where , , , .
In the absence of vasomotor adjustment, the effective area used for heat exchange area(¯ (i) c ∩¯ a ) does not depend on T a , and the coefficients D i , i = 1 . . . 5, are constants. Therefore, the surface temperature T (i) as a function of the ambient temperature T a consists of a piecewise linear function, as illustrated in Figure 7. Because the coefficients D i , i = 1 . . . 5, are positive, Equation (74) indicates that the slope of this function is smaller than one for T − < T a < T TNZ 1 , whereas the the slope is one in the TNZ. The slope of the linear function for T TNZ 2 < T a < T + , given by 1 + D 4 in Equation (74), depends on the term C 2 −αBL. This means that, depending on the efficiency of the evaporative cooling mechanism, the slope can be either smaller or greater than one. For this reason, the corresponding line shown in Figure 7 is merely illustrative.
When vasomotor adjustment is present, we obtain the same equations as before, but now with T TNZ , T + ]. In addition, the effective area used for heat exchange is altered due to contraction or dilation of blood vessels in response to the body temperature. According to Equation (74), changes in area(¯ (i) c ∩¯ a ) as a function of T a yields a nonlinear dependence between T (i) and T a . We conclude that the presence or absence of vasomotor adjustment leads to a striking change from linear to nonlinear behaviour of T × T a . This change explains why some authors have  (74). The slope of the linear function in the interval [T − , T TNZ 1 ] is smaller than one because D 1 > 0, while the slope in the interval [T TNZ 1 , T TNZ 2 ] is one. In the interval [T TNZ 2 , T + ], the slope can be either smaller or greater than one depending on physical and biological parameters of the animal being considered. In particular, it depends on the efficiency the evaporative cooling mechanism. Therefore, the line shown in the interval [T TNZ 2 , T + ] is merely illustrative.
used experimental data of T × T a acquired from thermal infrared assays, as an evidence of whether an animal relies on vasomotor adjustment for thermoregulation. A natural question in the present setting is whether steady-state solutions are stable solutions of the corresponding dynamical system. This is expected to be true, at least in a particular range of T a because we used feedback controls only. However, a mathematical proof of this statement is rather difficult and it would require a careful analysis of a nonlinear dynamical system in infinity dimensions (in our case we have an integro-partial differential equations model). A simpler analysis based on a linearization around a specific steady-state solution would require the study of the spectrum (eigenvalues) of the associated linearized integro-partial differential operator, which is also very hard. In any case, such analyses are out of the scope of the present paper.

Discussion and conclusions
We presented a dynamical four-layer model for thermoregulation of endotherms by taking into account heat production due to metabolic rate and heat exchange within the body given its internal structure. Our model assumes that heat exchange inside the body is due to conduction and heat transport by blood flow. In addition, heat is exchanged with the ambient through convection, radiation, and evaporation from the respiratory tract. In addition, four possible thermal control mechanisms are taken into account: changes in (i) thermal insulation, (ii) metabolic rate, (iii) changes in blood fluxes due to blood vasomotor adjustment in tissues near the body core, (iv) changes in blood fluxes due to blood vasomotor adjustment in tissues near the surface, (v) evaporation from the respiratory tract, and (vi) superficial evaporation.
Numerical simulations could be used to understand the evolution of the dynamical model and the interplay among the several thermal control mechanisms in rather realistic situations (Dudley, Bonazza, & Porter, 2013). Numerical simulations could also be helpful to elucidate the consequences of the interactions of several body regions with different thermal control mechanisms, and how the steady-state solutions are dynamically approached.
We analysed the steady-state solutions of our model analytically as an attempt to fully understand the metabolism-ambient temperature U-shaped curve. Our findings are able to explain qualitatively many aspects observed in experimental data, which are not properly explained in classical models of thermoregulation in endotherms. We explained (1) the linear increase in the metabolically rate away from the TNZ; (2) the role of vasomotor adjustment as a thermal control mechanism capable of expanding the range of the TNZ; (3) season and size-dependent changes of the metabolism-ambient temperature curve, and (4) the linear and nonlinear dependence of surface temperature as a function of T a .
We believe that from the physical first principles and the physiological concepts revised here and applied to build our model we can lay the foundations for a solid understanding of thermoregulation phenomena in endotherms. For instance, in its present form, our model considers all the possibilities of the changes in the metabolic rate (e.g. shivering and non-shivering mechanisms are implicitly considered) in a unified way. As it is, our model could contribute to highly important models that are currently applied to medical science (Fraguela, Matlalcuatzi, & Ramos, 2015), which aim to determine the temperature inside the incubator that would thermally stabilize a premature newborn in the shortest time possible.
Later on, a more detailed model, distinguishing different mechanisms which increase the metabolic rate can be developed. Another possibility of improving our model is to include transpiration as a mechanism for thermoregulation. Furthermore, field conditions with, for instance, direct exposition to the sun or wind can be incorporated into our model by including terms associated with the dependence of the convection coefficient on the wind speed, and so on. A probably more complex thermoregulation mechanism to be incorporated is postural change in response to ambient temperature.
The mathematical richness of our model also raises interesting questions for future work. For instance, steady-state solutions (with T bc = T obj ) are expected to be the attractors of our dynamical model. Additional investigation is required to test whether this is true and which steady states are verified for a given initial condition determined by experimental data. However, examining our arguments we observe that the steady-state equations admit other equilibrium states with T bc = T obj . However, such other equilibria cannot be attainable as asymptotic limits as t → +∞ starting from arbitrary initial data if they evolve according the stated evolution equations. Such equations include the control equations with coefficients R depending on T bc (t) − T obj (as in expression (15) and (16) for the special case of the time evolution of the thermal insulation of the coat).

Notes
1. Notice that we changed the original notation from Kleiber (1972). 2. A bar over the name of a set denotes its closure. 3. In the present work we do not consider heat exchange due to direct insolation (direct exposure to the sun rays) or reflected radiation. 4. All temperatures in this paper are measured in Celsius degrees. 5. All temperatures in this paper are measured in Celsius degrees. 6. Different from the other three mechanisms, thermal control based on the vasomotor adjustment requires the controlled variable to change in the same direction as T bc − T obj to facilitate heat exchange. 7. For example g z (τ ) = e −aτ , a > 0. scholarships from Conselho Nacional de Desenvolvimento Científico e Tecnológico (140773/2013-4). Work by J. L. Boldrini and S. F. dos Reis is supported by research fellowships from Conselho Nacional de Desenvolvimento Científico e Tecnológico (S. F. dos Reis -303544/2011-2 and J. L. Boldrini -306182/2014-9).
We use Equation (A2) to rewrite Equation (26) as Therefore, A 1 can be written as Solving Equation (28), we have T p (z) = B 1 exp(−az) + B 2 exp(az) + T obj + q wρ bl c bl , ( A 5 ) with a = wρ bl c bl k . ( A 6 ) Combining Equations (A5) and (A2) in Equation (31), we obtain exp(−ah s )B 1 + exp(ah s )B 2 = C 1 , ( A 7 ) with C 1 = − q wρ bl c bl − q 2k h s 2 + A 1 h s . ( A 8 ) Combining Equations (A5) and (A2) in Equation (32), we obtain − a exp(−ah s )B 1 + a exp(ah s )B 2 = C 2 , ( A 9 ) with Solving Equations (A7) and (A9) for B1 and B2, we obtain Using the definitions of B1 and B2, the last equation can be written as where Notice that a, D i , and E j for 1 ≤ i ≤ 2 and 1 ≤ j ≤ 4, are constants that depend only on known physical and biological parameters. Next, by replacing Equation (A17) in Equation (A4), we have which can be used in Equations (A11) and (A12) to determine B 1 and B 2 in terms of known parameters, the metabolic rate per unit volume q, the additional metabolic rate per unit volume due to evaporative cooling q v and the vapour mass lost by respiration α v .