Homeostatic model of human thermoregulation with bi-stability.

All homoiothermic organisms are capable of maintaining a stable body temperature using various negative feedback mechanisms. However, current models cannot satisfactorily describe the thermal adaptation of homoiothermic living systems in a physiologically meaningful way. Previously, we introduced stress entropic load, a novel variable designed to quantify adaptation costs, i.e. the stress of the organism, using a thermodynamic approach. In this study, we use stress entropic load as a starting point for the construction of a novel dynamical model of human thermoregulation. This model exhibits bi-stable mechanisms, a physiologically plausible features which has thus far not been demonstrated using a mathematical model. This finding allows us to predict critical points at which a living system, in this case a human body, may proceed towards two stabilities, only one of which is compatible with being alive. In the future, this may allow us to quantify not only the direction but rather the extent of therapeutic intervention in critical care patients.

The discovery of homeostasis, the ongoing maintenance and defense of vital physiological variables such as blood pressure and blood sugar, is generally attributed to Walter Cannon, considered to be the primary author of this prominent principle underlying physiological regulation[1]. However, already in 1878, the French physiologist Claude Bernard described the processes of physiological control as the milieu interieur [2]. The evolution of the term “homeostasis” has a fascinating history, and generally covers a large array of physiological and behavioral responses which are elicited in a precisely coordinated way to facilitate the maximum adaptiveness of the living system to the surrounding environment. Although the term homeostasis was originally used to refer to processes within living organisms, it is currently widely used to describe autonomous control in any system and is generally understood as a way of maintaining the adaptiveness of any system, either a living one or a machine.

The core principle of the hypothesis of the homeostatic regulation of a biological parameter is that all responses act together in a highly coordinated fashion to adjust bodily parameters crucial for the survival and reproduction of the living system. This means that the homeostatic regulation is precise and targeted and requires constant fine-tuning, including multiple negative feedback mechanisms, as expressed by Walter Cannon in 1945[3].

Over the last decades, a large number of physiologically plausible models[4,5] have attempted to conceptualize homeostatic mechanisms using a systemic approach based on regulated variables, their set points, errors, and comparative information and their proxies. This is an extremely problematic approach, as there are countless physiological parameters that characterize the living organism at a given moment, and these countless parameters can further involve countless synergies or antagonisms. One may even argue that such an approach, based on atomizing the environmental influences in various categories and evaluating the influences of these drivers in different individuals, is reductionist to the point of no longer being useful. Therefore, while many explanations on the interaction of living organisms with the environment have been proposed, starting with research of Hans Selye[6,7], no completely satisfactory theory describing the interaction between organisms and the surrounding environment is currently available. Furthermore, no presently available methods feature a predictive function capable of making predictions about the thermal homeostasis of the body in given ambient temperatures with satisfactory precision.

In our previous study, we introduced stress entropic load (SEL), a novel variable designed to facilitate the objective physical measurement of the stress load of a living (human) body by monitoring energy and matter flows in the body using the proxy of entropy production[8]. Mathematically, SEL is a time-dependent function which may be calculated from different types of heat, temperature changes and gases exchange. In line with our initial definition, we further described an equation calculating approximate SEL for a short measurement interval. In this article, we present a novel approach to measure SEL development in humans, and we further test this approach on healthy male volunteers subjected to prolonged mental effort[9]. As the calculation itself is based on commonly measured variables, e.g. volumes of inspired/expired O2/CO2 and measurements of core/surface/ambient temperatures, the method is potentially usable in a practical setting, including under clinical conditions.

Our current goal is to link the concept of SEL with a model of human thermoregulation. Previous studies proposed several models which can be used to study mammal thermoregulation [4,10-12]. However, while simple phenomenological models[4] do not include all of the necessary variables for computing SEL, and thus cannot be used to model stress response during thermoregulation, complex simulation models [10-12] do include all of the necessary variables for computing SEL, but their analysis is limited by their complexity.

This study aims to i) demonstrate the utility of the stress entropic load for evaluating the adaptation strain on the organism; ii) use dynamical systems theory for estimating core/surface body temperature in a given ambient temperature; iii) describe a mechanism which facilitates the evaluation of stability in the homeostatic equilibrium in a wider range of ambient air temperatures.

Methods

In this article, we use dynamical modeling to describe human thermoregulation and compute SEL for each model’s possible trajectories. This approach enables us to simulate stress load during temperature changes in the human body. We propose a novel simple dynamical model describing the thermoregulation process. We perform standard equilibrium and bifurcation analysis of the dynamical model, specifically detecting a fold bifurcation[13]. A fold bifurcation may be associated with the hysteresis principle typical of living systems[14,15]. The analysis enables us to locate a homeostatic region[4] and an allostatic verge[8,16] in a dynamical system. A homeostatic region is a region defined by both parameters and initial conditions for state variables in which the homeostatic equilibrium is stable. An allostatic verge is a manifold at which a system enters or leaves a homeostatic region. We associate an allostatic verge with the hysteresis principle or, more generally, with catastrophe theory[14,15].

In this article, we cannot use the universal Golubitsky’s and Stewart’s theory[17] to analyze homeostatic equilibria in a dynamical system as a departure point because, in this theory, homeostasis is simplified. Instead of a process during which “some output variable remains approximately constant as an input parameter varies”[17], it presumes that the output variable remains constant. However, such a simplification does not work outside of idealized scenarios. In reality, no physiological variable is truly constant; there is always at least a small deviation from the optimal value. Therefore, we choose not use this simplification to study the efficiency of the homeostatic mechanism.

Results and conclusion

Model definition

A model of mammal thermoregulation is a typical example of a homeostatic system. A simple model of thermoregulation with heating saturation was proposed by Nijhout et al. 2014[4]. Their model centers on a basic assumption, i.e. that core body temperature remains approximately constant in a wide range of air temperatures [18]. Our model further develops this idea by assuming that core body temperature is affected by air temperature , mostly through an additional layer of skin. Hence, both air temperature , and core temperature affect skin temperature in our model. Skin temperature is thus a key variable for the thermodynamic modeling of entropy production[9,16]. In our model, the only direct heat transfer between the body core and the ambient air takes places through breathing. It is proportional to the volume of inhaled air, and, consequently, the volume of oxygen uptake . The volume of oxygen uptake is not homeostatically regulated, meaning no known sensor in the body is capable of directly measuring it[5]. In our conceptualization of the model, we use the volume of oxygen uptake as a parameter, as i) there is no direct sensor of oxygen uptake and ii) the mechanisms involved in the regulation of oxygen uptake seem to be too complex to analyze using current methods. Our mathematical model of mammal thermoregulation may thus be presented as follows:where Tables 1, 2 denote time-dependent state variables and parameters, and is a saturating function. Generally, when no heat loss occurs during the heat transfer between body core and skin, rate is equal to rate .

Table 1

State variables for model (1), (2).

State variable	Description
T_C	Core body temperature K
T_S	Skin temperature K

Table 2

Parameters for model (1), (2).

Parameter	Description	Default value
k_1	Rate of heat transfer between skin and ambient air min^-1	0.2
k_2	Rate of heat transfer between skin and core min^-1	0.2
k_3	Rate of heat transfer between core and skin min^-1	0.2
k_4	Rate of heat transfer between core and ambient air through breathing l^-1	0.06
T_A	Ambient air temperature K	300
V_O_2	Oxygen consumption in lmin^-1	0.3

Heat transfer rates , describe the interaction of the body and environment through skin temperature and core body temperature . Similarly, heat transfer rates , describe the interaction between core body and skin temperature. We assume that the rate is significantly lower than the other rates, i.e. , , . Saturating function describes the inner mechanism of the body thermoregulation regardless of the external environment. The term saturating function was introduced by Nijhout et al. 2014[4]. Heat transfer rates are illustrative and their estimation constitutes a task for future studies.

In general, saturating function should have the following properties:

is smooth in .

There is no saturation at core temperature set point : .

When the core temperature is lower than set point , saturation increases the core temperature: ; when the core temperature is higher than set point , saturation reduces the core temperature: .

The saturation function is almost zero for a core temperature far away from the set point: as . Moreover, core temperature thresholds , exist, so that if or , then .

Function has only one maximum for , ; and one minimum for , .

An example of saturating function follows:

The proposed saturating function is similar to functions used in[4]; for computational reasons, we proposed a different formula to ensure is sufficiently smooth. The mechanism described by allows the body to bring the temperature back to a set-point [4]. The proposed saturating function is symmetric with respect to , but we do not require any symmetries in general. Parameters of function are denoted by Table 3.

Table 3

Parameters for saturating function (3).

Parameter	Description	Default value
p	Shape (peak) of saturating function f_S min^-1	2
v	Shape (reciprocal of variability) of saturating function f_S K^-2	0.5
TCSet	Core body temperature set-point [K]	310.15

To compare saturating functions of different shapes, we define the area below the saturation curve as

Therefore, the saturation functions of different shapes with the same area below the saturation curve exist. Figure 1 shows two different saturation functions with the same area below the saturation curve (). Example 1 shows a saturation curve of a higher variability () and smaller peak () in comparison to Example 2 (, ). Example 2 enables core temperature to stabilize closer to the set point than Example 1.

Figure 1

Different saturating functions with area below curve for parameters , (Example 1) and , (Example 2).

We can reduce the number of parameters of model (1), (2) using the following transformation , , , , :

The total number of parameters is thus reduced from nine to seven. In section Model analysis, the results are presented using the original parameters, i.e. before the transformation. The transformation aims to ensure that we only study relevant parameter dependencies.

Model analysis

The formal mathematical analysis of the model follows.

First, we analyze the equilibria of system (1), (2) for general sufficiently smooth function . By setting the right-hand sides of equations (1), (2) equal to zero and by then simplifying the resulting two equations, we obtain an equation for core body temperature equilibria .

The stability of solutions can be determined using the Jacobi matrix J of system (1), (2). Using stability conditions and , we obtain stability condition

For arbitrary , equation (4) has at least one solution with respect to . Moreover, if exists so that

it is then possible to establish for which equation (4) has three solutions with respect to . The solutions are either The condition (5) indicates that solutions 1 and 3 are stable. Depending on , typically, one of the following critical cases occurs: either solutions 1 and 2 collide, forming a saddle-node equilibrium , or , or solutions 2 and 3 collide, forming a saddle-node equilibrium , or .

solution 1: , solution 2: , and solution 3: ; for , and or

solution 1: , solution 2: , and solution 3: ; for , and .

The following paragraphs describe the consequences of this theorem, and present the results for as defined by equation (3). We analyzed system (1), (2), (3) with parameters from Tables 2 and 3 using equilibria and bifurcation detection and continuation methods in Matcont software[19,20], and using the DETools package in Maple software[21]. The goal of the analysis is to locate allostatic verges depending on parameters from 2 and 3. In different settings, allostatic verges can also depend on other quantities.

Depending on system parameters, two possible stable equilibria of system (1), (2), (3) exist. We can distinguish between them using core temperature at the equilibrium point. Core temperature at the first equilibrium point is close to the core temperature set-point . We thus call the equilibrium homeostatic. Core temperature at the second equilibrium point is close to ambient air temperature . We call this equilibrium non-homeostatic. We divided the entire parameter space into areas where either a homeostatic equilibrium exists and is stable (area 1), where a non-homeostatic equilibrium exists and is stable (area 2), or where both equilibria exist and are stable (area 3). By setting parameters to values from area 3, one additional saddle equilibrium exists, see Fig. 2. The separatrices of the saddle divide phase space, i.e. the space of state variables and , into two regions (region 1, region 2). Region 1 is a basin of attraction of the homeostatic equilibrium. Region 2 is a basin of attraction of the non-homeostatic equilibrium.

Figure 2

Phase portrait for area 3 for parameter values given by Tables 2 and 3.

The existence of bi-stability in area 3 is closely connected to a one-parameter dependent hysteresis phenomenon. Since in this model, hysteresis is a mechanism through which our system leaves or enters a homeostatic region, we analyze hysteresis to study allostatic verges. We identified the hysteresis phenomenon for air temperatures smaller than the core temperature set point . Similarly, it can be described for air temperatures larger than the core temperature set point . To identify the hysteresis phenomenon from either the trajectories of the dynamical system or the measured time series, we focus on two model scenarios of parameter change: (i) increasing and (ii) decreasing ambient air temperature . In the first scenario we start at a stable non-homeostatic equilibrium for below 307 K. We slowly increase . System (1), (2), (3) stabilizes on a non-homeostatic equilibrium until reaches an allostatic verge of 307 K. After crossing this allostatic verge, system (1), (2), (3) subsequently stabilizes on a homeostatic equilibrium until reaches the next allostatic verge at 321 K. In the second scenario we start at a homeostatic equilibrium for between 299 K and 321 K. We slowly decrease . System (1), (2), (3) stabilize on a homeostatic equilibrium until reaches an allostatic verge of 299 K. After crossing this allostatic verge, system (1), (2), (3) subsequently stabilizes on a non-homeostatic equilibrium. Allostatic verges are different for increasing and decreasing ambient air temperature scenarios (307 K and 299 K).

Moreover, allostatic verges are also dependent on other parameters of system (1), (2), (3). To describe the dependencies of allostatic verges on parameters, we studied several two-parameter sections of the parameter space. We provide the analysis in the following paragraphs.

We studied the equilibria of system (1), (2), (3) depending on ambient air temperature and oxygen consumption . For smaller area 3 gets wider, the change does not affect the size of area 1, see Fig. 3 and Table 4.

Figure 3

Transversal section through areas of stability of equilibria for parameter values given by Tables 2 and 3 and free parameters and .

Table 4

Estimated allostatic verges of ambient air temperature on the borders of areas 1, 2 and 3 for selected values of .

	Allostatic verges of T_AK on borders of areas
V_O_2	Area 2–3	Area 3–1	Area 1–3	Area 3–2
0.1	298	307	313.5	322.5
0.4	299.5	307	313.5	320.5
0.7	301	307	313.5	319.5

In addition to and , the areas of stability differ depending on the shape of saturation function . We study the dependence of ambient air temperature and the shape of . For the same area below the saturation curve, the higher peak p and the smaller variability enlarge area 3 and narrow area 1. Table 5 presents an example of such a dependency with the area below the saturation curve .

Table 5

Estimated allostatic verges of the ambient air temperature on borders of areas 1, 2 and 3 depending on the shape of the saturation curve (described by some values of parameters p, v) for the area below the saturation curve .

Shape		Allostatic verges of T_A[K] on borders of areas
v	p	Area 2–3	Area 3–1	Area 1–3	Area 3–2
1	4	295	307.5	312.5	325
0.9	3.6	296	307.5	312.5	324.5
0.8	3.2	296.5	307.5	313	323.5
0.7	2.8	297.5	307.5	313	323
0.6	2.4	298	307	313	322
0.5	2	299	307	313.5	321
0.4	1.6	300	306.5	314	320
0.3	1.2	301	306	314	319
0.2	0.8	302	305.5	315	318
0.1	0.4	303.5	304	316.5	317

The following Table 6 summarizes obtained results in terms of the homeostatic region. The necessity to characterize this homeostatic region by state-variables comes from the bi-stability between two possible equilibria for parameter area 3. The bi-stable mechanism enables the body to remain in a homeostatic equilibrium for more extreme temperatures during air temperature increase or decrease. On the other hand, the mechanism causes the destabilization of homeostatic equilibrium through even small perturbations.

Table 6

Homeostatic regions for mammal thermoregulation model (1) and (2).

	Parameters	State-variables
Homeostatic region 1	Area 1	Arbitrary T_C, T_S
Homeostatic region 2	Area 3	Region 1

For areas see Fig. 3 and Tables 4 and 5. For regions see Fig. 2.

Finally, having studied a dynamical system of mammal thermoregulation (1), (2), (3), we proceed to link our results to the stress entropic load concept. Assuming we are studying a short term event, we can use system (1), (2), (3) and the simplified formula published by Zlámal et al. 2018[9] to compute the immediate change of the stress entropic load ; see Supporting information for the formula. Our model describes a dynamic of core temperature and skin temperature depending on several parameters, including ambient air temperature and oxygen consumption . In the formula, the only free variable left is the CO liberation rate. As CO is the only free variable in the formula, the human body can thus optimize its stress entropic load by changing the CO liberation rate. The mechanism of regulating CO can therefore ensure that the following holds true in the homeostatic equilibrium: . Therefore, SEL production is zero at the homeostatic equilibrium, which is in perfect agreement with the underlying theory that SEL is directly related to stress levels recorded at a given point in time.

Discussion

In this paper, we propose a novel model of human thermoregulation using the concept of stress entropic load as a departure point. In doing so, we propose changes to the usual model of human thermoregulation in order to i) include all necessary variables for SEL computation into the model, such as core temperature , skin temperature , ambient temperature , oxygen consumption , and ii) model different adaptation strategies. Our results suggest that humans may use two types of adaptation strategies with respect to thermoregulation. The first adaptation strategy type takes into account adaptation costs measured by the area below the saturation curve while the second is based on CO production.

As for the first strategy: We studied the bi-stable area in system (1), (2) depending on the shape of the saturating function and oxygen consumption . The bi-stable mechanism enables the body to remain in a homeostatic equilibrium state for more extreme temperatures, but it leads to the possible destabilization of the homeostatic equilibrium through even small perturbations. With the decrease of oxygen consumption , the area of bi-stability (area 3) increases, but the homeostatic equilibrium area of stability remains approximately the same, see Fig. 3. Therefore, as oxygen consumption decreases, body can remain in homeostasis for more extreme temperatures. Similarly, depending on the shape of saturation sizes of both area 1 and area 3 changes. For the same area below the saturation curve, the saturating function with high peaks and small variability enlarges area 3 and narrows area 1; see Table 5. Therefore, body can remain in the homeostatic equilibrium for more extreme temperatures; however, the local stability is more fragile due to the limited resources needed to maintain homeostasis.

The thermoregulation mechanism described in our paper differs from classical chair-like relationships[4,17]. Golubitsky and Stewart[17] assumed that the derivative of the variable exhibiting homeostasis with respect to the parameter vanishes at homeostatic equilibrium. For the chair-like relationship, they assumed two additional conditions[17]. In our example of system (1), (2) with (3) the derivative of core temperature with respect to ambient temperature does not vanish for any saturating function shape (defined by parameters p and v). However, the implicit derivative of in the homeostatic equilibrium tends to zero as p tends to infinity. The theory is therefore not sufficient to describe this type of homeostatic system.

As for the second adaptation strategy which builds on CO production, we can use the dynamical systems of mammal thermoregulation to compute the change in stress entropic load (SEL) for short-term events [9]. Interestingly, our results show that the body can attain zero change in SEL at homeostatic equilibrium by manipulating the CO liberation rate. This is supportive of the notion that changes in SEL are representative of the current adaptation costs of the organism. In the light of our results, it is essential to mention that respiratory control in humans is primarily maintained via the partial pressure of CO in peripheral blood. The complete adaptation process consists of balancing both costs of adaptation, represented in the model by saturating function , and stress entropic load (SEL) through the CO liberation rate. Our findings may thus be interpreted to signify that manipulation with CO levels could reduce adaptation costs for thermoregulation, which may result in an increase of energy reserves for other purposes such as healing and regeneration.

Supplementary Information.