Human simulation model for predicting body temperature, blood pressure, and blood flow rate under various thermal exposures — Part 1: development of physical model for average Japanese male in their 20s

TL;DR

This paper develops a human body simulation model to predict body temperature, blood pressure, and blood flow rate in an average Japanese male in their 20s under thermal conditions.

Contribution

A new physical model combining thermal and cardiovascular simulations for predicting physiological responses in a specific demographic.

Findings

The physical model accurately reproduces experimental body temperature, blood flow rate, and blood pressure values under thermally neutral and supine conditions.

The model integrates a thermal network and cardiovascular system to simulate heat and blood transfer phenomena.

Abstract

In the mechanism of thermal health problems, not only body temperature but also blood pressure and blood flow rate are key factors. Human body simulation is a useful tool for predicting these physiological quantities under various conditions. Therefore, we aimed to develop a new human-body simulation model to predict body temperature, blood pressure, and blood flow rate. Our model consists of a “physical model” that simulates the phenomena of heat and blood transfer based on the physical properties of the human body, and an “autonomic regulation model” that simulates the phenomena in which the physical properties of the human body are regulated by the activity of the autonomic nervous system. The physical model consists of a “thermal network (TNW) model” that simulates heat flow and body temperature, and a “cardiovascular (CV) model” that simulates blood flow and blood pressure. In this study, we focused only on developing the physical model. In addition, based on experimental data from a recent study, parameter settings were derived to simulate an average Japanese male in their 20s. By applying the derived parameter settings to the physical model developed in this study, we confirmed that the calculated body temperature, blood flow rate, and blood pressure under a thermally neutral and supine condition reproduced the experimental values well.

The online version contains supplementary material available at 10.1007/s00484-025-03072-6.

Figures

Overall structure of the human simulation model

Schematic diagram of the thermal network model a: segmentation and decomposition of the body; b: configuration of the thermal network within an example segment (leg)

Schematic diagram of the cardiovascular model a: overview of the 1-D model; b: overview of the 0-D model, c: configuration of a typical peripheral circulation; d: configuration of the heart and pulmonary circulation

Simulation results and experimental data under the baseline condition (error bar: standard deviation of the measurements) a: body temperature; b: arterial blood flow rate and cardiac output; c: skin blood flow rate; d: blood pressure; e: time variation of simulated blood pressure (no experimental data for time variation of blood pressure)

Full text

Introduction

Inappropriate thermal environments can cause health problems. In hot environments, not only does the problem of hyperthermia arise, but problems of a decrease in orthostatic tolerance (Rowell 1986) and the occurrence of heat syncope (Ministry of the Environment of Japan 2022) also arise due to decreases in blood pressure (BP) and blood flow rate (BFR) in the brain. The decreases in BP and BFR in the brain are caused by vasodilation and increased BFR in the skin and extremities for thermoregulation. In cold environments, not only does the problem of hypothermia arise, but the problem of an increase in BP also arises due to vasoconstriction and reduced BFR in the skin and extremities for thermoregulation. The increase in BP can lead to the development of cerebrovascular disease in patients with hypertension (Umemura et al. 2019). In addition, rapid fluctuations in BP due to sudden changes in ambient temperature can cause cerebral stroke and heart attack (Tochihara 2012a, b). Thus, changes in body temperature (BT), BP, and BFR can directly trigger thermal health problems.

BP and BFR are also important for correctly assessing the risks of hyperthermia and hypothermia. It is well known that BFR plays an important role in maintaining the core temperature of the body; however, it is not solely related to BT regulation. BFR can change due to postural changes, for example, from a seated to a standing position, because they induce changes in BP distribution in the body (Rowell 1986). It can also change due to exercise, which requires an increase in oxygen supply to various parts of the body (Rowell 1986). Therefore, BFR is not necessarily controlled in the same manner under the same BT conditions, and the effects of posture and exercise on BP and BFR compete with BFR regulation to maintain BT. Consequently, the actual BFR is realized by compromising these factors.

Thus, BT, BP, and BFR are influenced not only by the thermal environment but also by individual behavior. Moreover, the physical constitution and predisposition of individuals are also important influencing factors. The risk of thermal health problems varies depending on the combination of these influencing factors. Human-body simulations are useful for assessing the risks brought about by various and complicated combinations of these factors.

Various human-body simulation models have been developed (Katic et al. 2016). The Stolwijk model (Stolwijk 1971) is the most representative model that considers the shape of the human body and can predict temperature distribution throughout the body. Autonomic thermoregulatory functions of the human body, such as vasomotion, sweating, and shivering heat production, have also been modeled. In recent years, advanced models, such as those proposed by Tanabe (Tanabe et al. 2001, 2006; Kobayashi et al. 2013; Takahashi et al. 2021), Huizenga and Zhang (Huizenga et al. 2001; Zhang et al. 2001), Fiala et al. 1999, 2001), have become well known. These models are called thermophysiological models, and are often used to evaluate human thermal stress in indoor and outdoor environments. However, these models primarily focus on predicting the BT, and their validity for predicting the BFR has not been sufficiently tested. This is because there was a paucity of available measurement data for the BFR. In addition, these models do not have the ability to predict BP.

Simulation models for BP and BFR have been developed to understand the physical phenomena in the cardiovascular system and support the diagnosis of cardiovascular diseases, such as arterial and aortic stenosis. The earliest studies were by Noordergraaf et al. (1963) and Westerhof et al. (1969). They proposed an analog computer model in which the systemic arterial tree was replaced with an electrical circuit. With reference to their model, an analytical lumped-parameter model (Avolio 1980) and one-dimensional (1-D) numerical fluid dynamics models (Schaarf and Abbrecht 1972; Stergiopulos et al. 1992; Olufsen et al. 2000) of the arterial tree were developed. Formaggia et al. (2006) and Mynard and Nithiarasu (2008) coupled a 1-D numerical fluid dynamics model of the arterial tree with a lumped-parameter model of the heart. Liang et al. (2009) completed the entire model of the closed-loop structure of the cardiovascular system. They coupled a 1-D numerical fluid dynamics model of the arterial tree with a lumped-parameter model of other elements, such as the peripheral circulations, venous system, and cardiopulmonary system. Thus, simulation technology for the cardiovascular system has already been developed outside thermophysiological modeling.

A few previous studies combined the aforementioned thermophysiological and cardiovascular models. Ghaddar and her colleagues (Salloum et al. 2007; Karaki et al. 2013; Rida et al. 2014) generated a model that combined the Stolwijk and Avolio models. Coccarelli et al. (2016, 2018) developed a model that combined the Fiala model and Mynard and Nithiarasu model. Liu et al. (2023) created a model that combined the Stolwijk model with a simple lumped-parameter cardiovascular model. However, the validity of their BFR predictions remains uncertain because their BFR predictions are essentially based on existing thermophysiological models. In addition, coupled cardiovascular models do not consider the venous system, or the cardiovascular model is overly simplified, making it impossible to make valid BP predictions for various conditions with different postures and thermal environments. Therefore, we aimed to develop a more realistic human-body simulation model that simultaneously predicts BT, BP, and BFR by coupling thermophysiological and cardiovascular models.

In human-body simulations, it is reasonable to treat the heat and blood transfer phenomena based on the physical properties of the human body separately from the phenomena in which these physical properties are regulated by the activities of the autonomic nervous system. This is a common concept in conventional thermophysiological models, which consist of a controlled system and a controlling system (Stolwijk 1971), or a passive system and an active system (Fiala et al. 1999, 2001). This study focused on developing a model that simulates the transfer of heat and blood based on the physical properties of the human body. The other model, which simulates phenomena in which the physical properties of the human body are regulated by the activity of the autonomic nervous system, will be addressed in the next paper. In addition, standard parameter settings were derived in this study. Although the ultimate goal of our present and future studies is to provide fine-tuned risk assessments according to differences in personal characteristics, it must be impossible to prepare tailor-made personalized parameter settings for every individual. Therefore, we aimed to prepare standard model parameter settings for personalizing later according to age, pre-existing diseases, etc. We chose an average Japanese male in their 20s as the standard person. This is because detailed experimental data on BFRs for this type of person were recently obtained (Goto et al. 2024) and could be used as reference data for deriving the parameter settings.

Structure and mathematical modeling of human simulation

Overall structure of the human simulation model

The overall structure of the proposed model is illustrated in Fig. 1. This model consists of a “physical model” that simulates the phenomena of heat and blood transfer based on the physical properties of the human body (i.e., including the phenomena of temperature and blood pressure propagation) and an “autonomic regulation model” that simulates the phenomena in which the physical properties of the human body are regulated by autonomic nervous system activity. The physical model consists of a “thermal network (TNW) model” that simulates heat flow and body temperature, and a “cardiovascular (CV) model” that simulates BFR and BP. The BFR calculated in the CV model was used as the BFR in the TNW model. This study only dealt with the physical model.Fig. 1. Overall structure of the human simulation model

Thermal network (TNW) model

The TNW model was developed based on the Stolwijk model (1971). In this model, the human body is divided into 16 segments, each of which is subdivided into four compartments: core, muscle, fat, and skin (Fig. 2a). However, the head is divided into the brain, skull, muscle, fat, and skin, and the chest is divided into the left heart, right heart, lung, core, muscle, fat, and skin. Arteries and veins are located in all segments, and superficial veins are located in the limb segments. In this study, each of the abovementioned compartments and blood vessels is referred to as a node.Fig. 2. Schematic diagram of the thermal network model a: segmentation and decomposition of the body; b: configuration of the thermal network within an example segment (leg)

Each compartment transfers heat to and from the adjacent compartments in the same segment by conduction (Fig. 2b). Blood flows into each compartment from the artery of the same segment, and blood flows into the vein of the same segment from the compartment. For the limb compartments, part of the blood from the compartments flows into the superficial vein. The heat balance equation for compartment I is expressed as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{c}C_{I}\frac{dT_{I}}{dt}=M_{I}+\sum_{J}\left\{\left({TC}_{I,J}+\rho c Q_{I,J}\right)\left(T_{J}-T_{I}\right)\right\}\\-C_{res}-E_{res}-E_{sk,I}+q_{I,ar}+q_{I,ve}\end{array}$$\end{document}

where QI,J is the BFR from node J to node I, and If QI,J is negative, it is set to zero. Sensible heat loss Cres and latent heat loss Eres due to respiration are included in the equation only when node I is the lung, and latent heat loss from the skin Esk,I is included in the equation only when node I is the skin. If node J is an artery or vein and is adjacent to compartment I, the thermal conductance TCI,J between I and J is set to zero, and qI,ar and qI,ve, described below, are included in the equation.

The Stolwijk model ignores heat exchange through blood vessels between the heart and body segments, whereas the model in this study considers heat exchange through blood vessels. The artery node of each segment receives blood from the artery node of the neighboring body segment on the proximal side. However, the artery nodes of the abdomen, neck, and upper arm segments receive blood from the left heart compartment of the chest segment. The vein and superficial vein nodes of each segment receive blood not only from the compartments of the individual segments but also from the vein and superficial vein nodes of the neighboring body segments on the distal side. In addition, blood from the vein and superficial vein nodes of the thigh segments flows into the vein node of the abdomen, and blood from the vein nodes of the neck, abdomen, upper arm segments, and superficial vein nodes of the upper arm segments flows into the right heart compartment of the chest segment. Blood flows from the right heart to the left heart through the lung compartment.

Except for the head, hands, and feet segments, heat exchange in the artery and vein with adjacent compartments is considered as piston flow, and heat gains from adjacent compartment I in the artery and vein are calculated using the following equations:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{ar,I}=\rho c Q_{ar}\left\{1-\mathrm{exp}\left(-\frac{{TC}_{ar,I}}{\rho c Q_{ar}}\right)\right\}\left(T_{I}-T_{ar\_prox}\right)=-q_{I,ar}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{ve,I}=\rho c Q_{ve}\left\{1-\mathrm{exp}\left(-\frac{{TC}_{ve,I}}{\rho c Q_{ve}}\right)\right\}\left(T_{I}-T_{ve\_dist}\right)=-q_{I,ve}$$\end{document}

where Qar and Qve are the BFRs in the artery and vein, respectively. TCar,I and TCve,I are the composites of the convective heat transfer coefficient in the vessel and thermal conductance between the vessel and adjacent compartment.

Counterflow heat exchange between the artery and vein was also considered. Based on the same theory as for the heat exchangers in HVAC systems, this heat transfer can be calculated using the following equation:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{c}q_{ar,ve}=\rho c Q_{ve}\frac{1-\mathrm{exp}\left\{-\frac{{TC}_{ar,ve}}{\rho c Q_{ve}}\left(1-\frac{Q_{ve}}{Q_{ar}}\right)\right\}}{1-\frac{Q_{ve}}{Q_{ar}}\mathrm{exp}\left\{-\frac{{TC}_{ar,ve}}{\rho c Q_{ve}}\left(1-\frac{Q_{ve}}{Q_{ar}}\right)\right\}}\\\times\left(T_{ve\_dist}-T_{ar\_prox}\right)=-q_{ve,ar}\end{array}$$\end{document}

where TCar,ve is the composite of the convective heat transfer coefficients in the artery and vein as well as the thermal conductance of the tissue between the artery and vein.

If we neglect the time-derivative term and the interaction of heat transfer in Eqs. (2)–(4), the heat balance of the artery and vein can be expressed as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho c Q_{ar}\left(T_{ar\_prox}-T_{ar}\right)+q_{ar,adj}+q_{ar,ve}=0$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{c}\rho cQ_{ve}\left(T_{ve\_dist}-T_{ve}\right)+\sum_J\left\{\rho cQ_{ve,J}\left(T_J-T_{ve}\right)\right\}\\\:+q_{ve,adj}+q_{ve,ar}=0\end{array}$$\end{document}

where qar,adj and qve,adj represent the heat gains from the adjacent compartments of the artery and vein, respectively.

The artery and vein nodes of the head, hands, feet, and all superficial vein nodes were assumed to be perfect mixing blood pools. The heat balance equation for these vessels is the same as that in Eq. (1). The temperatures of the individual body nodes can be calculated by discretizing Eq. (1) with respect to time, and solving Eqs. (1), (5), and (6) simultaneously.

Cardiovascular (CV) model

Our CV model was developed based on the model proposed by Liang et al. (2009). This model consists of a one-dimensional model (1-D model) that reproduces the major arteries and a lumped-parameter model (0-D model) that reproduces the peripheral circulation, vena cava, heart, and pulmonary circulation.

The 1-D model consists of 77 major arteries throughout the body. The ID numbers and names of the arteries comprising the 1-D model are shown in Fig. 3a. The brachial branch (No.11, 29), intercostal (No.38), costal (No.40), and genicular (No.63, 78) arteries were simplified from multiple arteries that run in a complex manner and replaced by a single artery. The governing equations for the 1-D model are as follows:Fig. 3. Schematic diagram of the cardiovascular model a: overview of the 1-D model; b: overview of the 0-D model, c: configuration of a typical peripheral circulation; d: configuration of the heart and pulmonary circulation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}=0$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\partial Q}{\partial t}+\frac{\partial}{\partial x}\left(\frac{Q^2}{A}\right)+\frac{A}{\rho}\frac{\partial}{\partial x}\left(P+\rho g H\right)+\frac{2\pi r\nu}{\delta}\frac{Q}{A}=0$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P=\frac{Eh}{r_{0}\left(1-\sigma^{2}\right)}\left(\sqrt{\frac{A}{A_{0}}}-1\right)$$\end{document}

Eqs. (7) and (8) represent the mass and momentum conservation equations, respectively. Eq. (9) is a state equation that expresses the relationship between BP and the cross-sectional area of arteries. In Eq. (8), the momentum equation of the Liang model is modified so that the effect of gravity can be taken into account.

To calculate the BFR and BP for each artery, Eqs. (7) and (8) are solved using the finite difference method after substituting Eq. (9) into Eq. (8). In addition, the conservation of BFR and continuity of BP are satisfied at each connection of the arteries, that is, at each bifurcation in the arterial tree. The boundary conditions of the 1-D model, that is, the BFR at the root of the arterial tree and the BPs at the terminals of the arterial tree, can be provided by the calculation results of the 0-D model.

The 0-D model expresses the entire cardiovascular system, other than the arteries in the 1-D model, as electrical circuits. As shown in Fig. 3b, the terminal arteries of the 1-D model were connected to the peripheral circulations of the 0-D model. These peripheral circulations connect to the superior, inferior, or abdominal vena cavae and then to the cardiopulmonary system. Figure 3c shows the configuration of the typical peripheral circulation in our 0-D model. In the Liang model, there is only one route in the peripheral circulation. However, in our model, the peripheral circulation is divided into four routes (core, muscle, fat, and skin) for consistency with the TNW model. Each capillary and venule of each route and the vein into which the route rejoins are represented by a basic unit consisting of three elements (compliance C, viscosity resistance R, and inertial resistance L) that represent the vascular characteristics. However, for the arterioles, the four bifurcated units share a single compliance. The three venae cavae connecting downstream of the veins are represented by the same basic units as those of the veins. The BFR and BP at each site in the peripheral circulation and vena cava are calculated using the following equations:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{dV_{i}}{dt}=Q_{up}-Q_{i}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_{i}-P_{dw}=Q_{i}R_{i}+L_{i}\frac{dQ_{i}}{dt}+\rho g\left(H_{dw}-H_{i}\right)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_{i}=\frac{V_{i}}{C_{i}}$$\end{document}

Eqs. (10) and (11) represent the mass and momentum conservation equations, respectively. Eq. (12) is the state equation for the relationship between the BP and vascular volume.

The heart and pulmonary circulation of the 0-D model are illustrated in Fig. 3d. Here, elastance E is used instead of compliance C, and the viscoelastic resistance S and pressure-loss coefficient B of the heart valve are given as parameters. Elastance E refers to the ratio of the pressure change to the volume change (dP/dV), and elastance E and compliance C are the inverse of each other. The mass conservation equation in the heart is the same as that in peripheral circulation, Eq. (10), but the momentum and state equations are different as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_{i}-P_{dw}=Q_{i}R_{i}+{Q_{i}}^2 B_{i}+L_{i}\frac{dQ_{i}}{dt}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_{i}=\left(E_{A,i}e_{i}\left(t\right)+E_{B,i}\right)V_{i}+S_{i}\frac{dV_{i}}{dt}$$\end{document}

where ei(t) is a function of time representing the contraction and relaxation of the heart. In the case of the atria:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_{i}\left(t\right)=\left\{\begin{array}{l}0.5\left\{1+\mathrm{cos}\left(\pi \frac{T+1-\left(T_{cs,i}+T_{cp,i}\right)}{T_{rp,i}}\right)\right\}\quad\left(0<T\le T_{cs,i}+T_{cp,i}+T_{rp,i}-1\right)\\0\quad\left(T_{cs,i}+T_{cp,i}+T_{rp,i}-1<T\le T_{cs,i}\right)\\0.5\left\{1-\mathrm{cos}\left(\pi \frac{T-T_{cs,i}}{T_{cp,i}}\right)\right\}\quad\left(T_{cs,i}<T\le T_{cs,i}+T_{cp,i}\right)\\0.5\left\{1+\mathrm{cos}\left(\pi \frac{T-\left(T_{cs,i}+T_{cp,i}\right)}{T_{rp,i}}\right)\right\}\quad\left(T_{cs,i}+T_{cp,i}<T\le 1\right)\end{array}\right.$$\end{document}

In the case of the ventricles:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_{i}\left(t\right)=\left\{\begin{array}{l}0\quad\left(0<T\le T_{cs,i}\right)\\0.5\left\{1-\mathrm{cos}\left(\pi \frac{T-T_{cs,i}}{T_{cp,i}}\right)\right\}\quad\left(T_{cs,i}<T\le T_{cs,i}+T_{cp,i}\right)\\0.5\left\{1+\mathrm{cos}\left(\pi \frac{T-\left(T_{cs,i}+T_{cp,i}\right)}{T_{rp,i}}\right)\right\}\quad\left(T_{cs,i}+T_{cp,i}<T\le T_{cs,i}+T_{cp,i}+T_{rp,i}\right)\\0\quad\left(T_{cs,i}+T_{cp,i}+T_{rp,i}<T\le 1\right)\end{array}\right.$$\end{document}

For the pulmonary circulation, Eqs. (10) and (11) are used as the mass and momentum equations, respectively, and the following equation is used as the state equation:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_{i}=E_{0,i}\varPhi_{i}\mathrm{exp}\left(\frac{V_{i}}{\varPhi_{i}}\right)+S_{i}\frac{dV_{i}}{dt}$$\end{document}

The BFRs and BPs in the 0-D model can be calculated by substituting the state equation into the momentum conservation equation, discretizing the mass and momentum conservation equations with respect to time, and numerically solving them. The boundary conditions of the 0-D model, that is, the BFRs at the inlets of the peripheral circulations and BP at the outlet of the aortic valve, can be obtained from the calculation results of the 1-D model.

Derivation of standard parameter settings for the physical model

Parameter settings of TNW model

The parameters to be provided in the TNW model are the heat capacity and metabolic rate of each node and the thermal conductance between nodes. To reproduce the average Japanese male in their 20s according to the best available information, these values were derived based on the voxel model of Nagaoka et al. (2004). Nagaoka et al. provided voxel model data with 2 mm spatial resolution and 55 types of tissues and organs based on MRI data of volunteers, matching the average height and weight of adult men and women. In this study, a male model (TARO) was used to derive information on body configuration. The height and weight of the male model (TARO) are 172 cm and 65 kg, respectively. The surface area and volume of each segment and the volume of each compartment obtained using this voxel model are listed in Table S1 in the Supplementary Material.

The heat capacity and metabolic rate were calculated based on the volume of each component, as well as the volumetric specific heat, density, and basal metabolic rate per unit mass of each tissue or organ (Table S2, Supplementary Material). The calculated heat capacity and metabolic rate of each node are listed in Table 1. For the arteries and veins, the heat capacity was not set because the heat balance was calculated using Eqs. (5) and (6), neglecting the time-derivative term. The heat capacity of superficial veins was also ignored in this study.Table 1. Heat capacity, basal metabolic rate, and simplified geometry of each node in the TNW modelSegmentNodeHeat capacity [J/K]Basal metabolic rate [W]Segment length [m]Compartment or vessel radius [m]Segmental sum of vessel surface area [m^2^]HeadBrain5523.317.470–0.070Skull4208.61.1090.089Muscle2627.60.4860.097Fat1342.30.1270.103Skin958.90.1480.105Artery^*1^––––Vein^*1^––––NeckCore406.40.0950.0650.025Muscle1200.00.2220.047Fat282.30.0270.054Skin173.40.0270.056Artery^*2^––0.0040.016Vein^*3^––0.0040.019ChestRight heart801.44.5960.354–^*4^Left heart801.44.596–^*4^Lung7383.11.1560.058Core8924.92.0400.079Muscle14050.12.6010.099Fat3266.90.3090.106Skin1815.20.2800.108AbdomenCore26520.128.2860.5300.075Muscle27446.65.0820.101Fat7009.40.6630.111Skin2829.20.4360.113Artery^*2^––0.0070.046Vein^*3^––0.0080.055Upper arm (right or left)Core587.00.2140.2880.016Muscle3892.20.7210.038Fat1006.30.0950.045Skin640.80.0990.047Artery^*2^––0.0030.010Vein^*3^––0.0040.024Lower arm (right or left)Core305.00.1110.1760.015Muscle1929.30.3570.034Fat371.00.0350.039Skin343.70.0530.041Artery^*2^––0.0020.006Vein^*3^––0.0020.014Hand (right or left)Core193.20.0700.2570.010Muscle449.10.0830.016Fat270.40.0260.020Skin293.40.0450.022Artery^*2^––4.5 × 10^− 4^0.003Vein^*3^––5.4 × 10^− 4^0.004Thigh (right or left)Core1000.40.3650.2540.023Muscle11748.22.1750.068Fat2879.30.2720.080Skin987.10.1520.082Artery^*2^––0.0030.011Vein^*3^––0.0030.013Leg (right or left)Core1035.60.3780.2860.022Muscle5176.20.9580.045Fat1244.40.1180.053Skin740.70.1140.055Artery^*2^––0.0020.012Vein^*3^––0.0020.028Foot (right or left)Core703.30.2570.2880.018Muscle843.40.1560.023Fat660.10.0620.030Skin484.40.0750.032Artery^*2^––4.4 × 10^− 4^0.007Vein^3^––5.3 × 10^− 4^0.0081 The actual arteries and veins of the head run complicatedly inside and outside the skull and are difficult to simplify. Therefore, the radius and surface area were not set, and the thermal conductance between the artery/vein and other nodes was neglected.*2 The radius and surface area of the artery in each segment were given as the average radius and sum of the surface areas of the arteries of the CV model considered to be in that segment. For the hands and feet, the radius and surface area of the arteries were provisionally set to 1/4 of the adjacent segment (lower arms and legs) and 20% of the surface area of the core (bone), respectively, because no quantitative reference information could be found.*3 The radius and surface area of the vein were set to 1.2 times those of the artery, referring to Levick (2010). However, the surface area in the upper arms, lower arms, and legs was assumed to be 2.4 times that of the artery because two veins run parallel to one artery in these body regions.*4 The radius of the heart was not set because only the heat transfer between the heart and lung by blood flow was considered, and the heat transfer between the heart and lung by conduction was neglected.

To determine the thermal conductance between the nodes, we simplified the geometries of the segments, such as the head to a sphere, the hand and foot to a combined cylinder-hemisphere, and the other parts to cylinders. Table 1 presents the simplified geometries (length and radius of each segment and node) obtained so that the surface area and volume of each segment and the volume of each compartment were consistent with those in Table S1 in the Supplementary Material. Tables S3 and S4 in the Supplementary Material provide the specific calculation methods. The calculated thermal conductance values between the nodes are listed in Table S5 in the Supplementary Material.

Parameter settings of CV model

The parameters to be provided in the 1-D model are the length, radius at the reference state of each artery, and thickness and Young’s modulus of the wall of each artery. These parameters were set based on previous studies (Table 2). The length of the arteries was given with reference to Westerhof et al. (1969), which is the basis for many cardiovascular models, and similar values have been used in several recent studies (Olufsen et al. 2000; Wang and Parker 2004; Liang et al. 2009). In Westerhof, a human with a height of 175 cm and weight of 75 kg was assumed, which is almost the same as the voxel model of Nagaoka et al. (2004) with respect to height, but slightly heavier with respect to weight. For the radii of arteries, we referred to a recent study by Liang et al. (2009). Different researchers have provided different values for vascular wall thickness h and Young’s modulus E. These values are often adjusted by examining the calculated BP waveforms, and there is no clear rationale for setting these values. Therefore, the relative thickness h/r0 and Young’s modulus E of the vessel wall were given with reference to Westerhof et al. (1969).Table 2. Characteristics of each artery in the CV modelID No.Name of arteryLength l [m]^*1^Proximal radius r0,prox [m]^*2^Distal radius r0,dist [m] ^*2^Relative wall thickness h/r0^*3^Young’s modulus E [Pa] ^3^1Ascending aorta0.0400.015250.014050.112400,0002Brachiocephalic0.0340.006500.006200.137400,0003R. common carotid0.1770.004000.003700.167400,0004, 21Internal carotid0.1770.003000.002750.192800,0005, 22External carotid0.1770.002000.002000.232800,0006, 24Subclavian I0.0340.004250.004070.161400,0007, 25Vertebral0.1480.002000.002000.232800,0008, 26Subclavian II0.0150.004070.003990.163400,0009, 27Thyrocervical trunk0.1000.001000.001000.342800,00010, 28Subclavian III0.1700.003990.003170.173400,00011, 29Brachial branch0.1500.001200.001200.308800,00012, 30Brachial0.2350.003170.002300.198400,00013, 31Radial0.2350.001750.001400.266800,00015, 33Ulnar I0.0670.002150.002150.223800,00016, 34Interosseous0.0790.001000.001000.3421,600,00017, 35Ulnar II0.1710.002030.001800.238800,00019Aortic arch I0.0200.014050.013490.113400,00020L. common carotid0.2080.004000.003700.167400,00023Aortic arch II0.0390.013490.012460.114400,00037Thoracic I0.0800.012460.010580.117400,00038Intercostal0.0730.003000.003000.187400,00039Thoracic II0.0600.010580.009360.120400,00040Costal0.0730.002330.002330.214400,00041Thoracic III0.0600.009360.008280.124400,00042Celiac I0.0200.003500.003000.181400,00043Common hepatic0.0660.002750.002500.201400,00044Celiac II0.0200.003000.002500.196400,00045Gastric0.0710.001750.001500.261400,00046Splenic0.0630.002000.002000.232400,00047Abdominal aorta I0.0450.008280.007550.128400,00048Superior mesenteric0.0590.004000.003500.169400,00049Abdominal aorta II0.0150.007550.007320.130400,00050, 52Renal0.0320.002750.002750.196400,00051Abdominal aorta III0.0150.007320.007100.131400,00053Abdominal aorta IV0.0650.007100.006220.135400,00054Inferior mesenteric0.0500.002000.001750.241400,00055Abdominal aorta V0.0600.006220.005500.141400,00056, 71Common iliac0.0580.004000.003700.167800,00057, 72Internal iliac0.0450.002000.002000.232800,00058, 73External iliac0.1440.003700.003140.176800,00059, 74Deep femoral0.1260.002550.001860.2211,600,00060, 75Femoral0.2540.003140.002750.189800,00061, 76Descending genicular0.1350.001970.001440.2551,600,00062, 77Popliteal I0.1370.002750.002750.1961,600,00063, 78Genicular0.0500.001190.001190.3111,600,00064, 79Popliteal II0.0520.002750.002750.1961,600,00065, 80Anterior tibial0.3430.001750.001750.2501,600,00067, 82Posterior tibial I0.0250.002500.002500.2061,600,00068, 83Peroneal0.3180.001300.001300.2951,600,00069, 84Posterior tibial II0.3220.001750.001750.2501,600,0001 The length of the arteries was given with reference to Westerhof et al. (1969). For arteries not addressed by Westerhof, the length was estimated by referring to other studies (Avolio 1980; Liang et al. 2009) and human anatomy charts (Otani and Kawahara 2012; Tortora and Nielsen 2009). For the brachial branch (No.11, 29), intercostal (No.38), costal (No.40), and genicular (No.63, 78) arteries, the length was given expediently, because these were modeled as substitutes for multiple arteries running in a complex manner and could not be given accurately.*2 The radii of the arteries were given with reference to Liang et al. (2009). For arteries not addressed by Liang et al., the radii were estimated by referring to Avolio (1980) or by actual measurements (Lorbeer et al. 2018; Szpinda 2005). Expedient values were set for the radii of the brachial branch, costal, and genicular arteries for the same reason as the length.*3 The relative thickness and Young’s modulus of the vessel wall were given with reference to Westerhof et al. (1969). For arteries not addressed by Westerhof, the values were set by referring to similar arteries addressed by Westerhof.

The 0-D model requires parameters for the peripheral circulations, venae cavae, heart, and pulmonary circulation. These parameter values, particularly for the peripheral circulations and heart, should vary owing to autonomic regulation in response to thermal and postural conditions. This study considered only the state of the cardiovascular system under the condition of minimal autonomic regulation because autonomic regulation was beyond the scope of this study. We attempted to derive the parameter settings to reproduce the BFR and BP measured under the thermally neutral and supine condition (baseline condition) (Goto et al. 2024) by modifying the parameter values of Liang’s model. In addition, we primarily modified the parameter values for the peripheral circulations, because the BFR distribution in the body is mainly determined by those. The parameter values for the venae cavae, heart, and pulmonary circulation were not modified (Tables S6 and S7, Supplementary Material), but only the heart rate was slightly modified (from 60 to 62 beats/min).

Regarding the parameters of peripheral circulations, Liang et al. (2009) set the magnitude of viscosity resistance R of the distal arterial end, arteriole, capillary, venule, and vein connected to each artery at ratios of 17%, 54%, 21%, 6.5%, and 1.5% to the sum of those viscosity resistances (total viscosity resistance RT), respectively. In addition, they set the magnitude of the inertial resistance L of the arteriole, capillary, venule, and vein to be proportional to the square root of the total viscosity resistance RT, and the magnitude of the compliance C of the arteriole, capillary, venule, and vein to be inversely proportional to the total viscosity resistance RT. Therefore, we adjusted the RLC values while maintaining these relationships and reproducing the arterial BFRs measured under the baseline condition (Goto et al. 2024). However, we changed the ratio of R such that the arteriole had 90% of the RT, and the remaining 10% was allocated to the distal arterial end, capillary, venule, and vein. This is because the ratio of the viscosity resistance of the arteriole to the RT was estimated to be much larger, considering that the BFR is mainly controlled by the vasomotor activity of the arteriole, that skin BFR in hot environments can be ten times greater than that in a neutral environment, and that muscle BFR during exercise can be ten times greater than that at rest (Rowell 1986).

Furthermore, the peripheral circulation in our model differed from that in the Liang model – our model has four routes, whereas Liang’s model has only a single route, as described in the "Cardiovascular (CV) model" section. Therefore, we decomposed the single-route RLC circuit into a four-route RLC circuit. The decomposed RLC elements of the peripheral circulation were adjusted to reproduce a realistic BFR distribution, which was determined based on the measurements (Table S8, Supplementary Material), while maintaining the BPs of the capillary and venule before and after the decomposition. Consequently, the peripheral circulation parameter values were determined (Table S9, Supplementary Material).

Simulation results by the developed physical model

Using the model parameter settings derived as described above, we performed a simulation to reproduce the experiment by Goto et al. (2024) under the baseline condition. The TNW model was solved using an implicit method with a computational time interval of 10 s. The CV model was solved with a time interval of 0.0001 s, the 1-D model was solved using the two-step Lax–Wendroff method, and the 0-D model was solved using the Runge–Kutta method with a fourth-order accuracy. The average BFR per heartbeat, calculated by the CV model, was passed to the TNW model every 10 s and used as the BFR in the TNW model.

The ambient temperature and relative humidity were set at 28 °C and 40%, respectively. Other boundary conditions, such as heat transfer coefficients, were given based on our measurements using a thermal manikin (Table S10, Supplementary Material).

A comparison between the simulation results and experimental data is shown in Fig. 4a–d. The arterial and skin BFRs were reproduced by adjusting the 0-D model. The BTs and BPs simulated with these BFR values were in good agreement with the experimental values. Figure 4e shows the time variation of the simulated BP at several sites. The systolic/diastolic BPs at the upper arm and ankle shown in Fig. 4d correspond to the maximum/minimum BPs at the brachial and posterior tibial arteries shown in Fig. 4e, respectively.Fig. 4. Simulation results and experimental data under the baseline condition (error bar: standard deviation of the measurements) a: body temperature; b: arterial blood flow rate and cardiac output; c: skin blood flow rate; d: blood pressure; e: time variation of simulated blood pressure (no experimental data for time variation of blood pressure)

Discussion

The novelty of this study lies in the coupling of a closed-loop cardiovascular model, which handles the entire pathway of an actual circulatory system, with a thermophysiological model. In previous models by Ghaddar and her colleagues (Salloum et al. 2007; Karaki et al. 2013; Rida et al. 2014) and Coccarelli et al. (2016, 2018), the cardiovascular model did not include the venous system. Consequently, these models cannot account for effects such as gravity-induced venous expansion in the lower limbs when people are standing, the accompanying decrease in cardiac filling pressure and stroke volume (SV), or changes in venous volume and cardiac filling pressure caused by alterations in peripheral vascular resistance due to thermally induced vasodilation and vasoconstriction. Therefore, it is not possible to fairly simulate the changes in BP caused by variations in posture or thermal environments. Liu et al. (2023) employed a closed-loop cardiovascular model in their model but excessively simplified the structure of the arterial system. Therefore, it is impossible to accurately simulate the complex pressure wave propagation phenomena, including the reflection of pressure waves from the peripheral ends of the arterial tree. It is also impossible to accurately simulate situations in which thermophysiological responses, baroreflexes, or physical activity cause the BFR to increase or decrease non-uniformly within the body because the circulatory routes were not divided sufficiently. In this study, the CV model represented the arterial system in a form close to its actual structure, and the division of organs and tissues was sufficiently detailed. Therefore, it is possible to simulate BP propagation and various nonuniform blood flow regulation responses in a realistic manner. However, this study did not address the proposal of an autonomic regulation model for reproducing thermoregulatory responses and its coupling with the physical model.

Another advantage of our model is that it is based on actual BFR measurements throughout the body. Table 3 compares the BFR simulated by our model with the basal BFRs of the Stolwijk (1971) and Gordon (1976) models, which are frequently referenced in other conventional thermophysiological models. In our model, the sum of the BFRs in the head and limb segments was based on the measurements obtained using the ultrasound Doppler technique, whereas the BFRs of the segments in the Stolwijk and Gordon models were estimated by adding the BFR per unit mass of each tissue. The BFR of each segment of our model was roughly similar to that of the Gordon model but differed from that of the Stolwijk model. The BFR to the head segment in our model was greater than that of Stolwijk and Gordon, likely because neither Stolwijk nor Gordon considered the BFR to organs in the head other than the brain. Comparing the BFRs of body tissues between our model and Gordon’s, there were apparent differences in the BFRs for skin and fat; however, the sum of the BFRs for skin and fat was roughly similar. Gordon estimated a very low BFR of fat based on the fat metabolic rate. However, past studies (Larsen et al. 1966; Ardilouze et al. 2004) measured abdominal subcutaneous tissue BFR at 26 or 40 mL/min/kg (2.6 or 4.0 mL/min/100 g), although the room temperature and clothing conditions were unspecified. These values are much larger than the fat BFR of the Gordon model, and we believe that our model better represents the actual subcutaneous BFR conditions. Meanwhile, we found a fact when we simulated using our model with replacement of the BFR values from Goto’s to Stolwijk’s and Gordon’s. The fact was that the accuracy of the simulated BT distribution deteriorated with Stolwijk’s BFRs but slightly with Gordon’s BFRs (Fig. S1, Supplementary Material). This suggests that our BFR values are at least more valid than those of Stolwijk.Table 3. Comparison of basal blood flow rates of the body compartmentsSegmentCompartmentBFR [mL/min]BFR per unit mass [mL/min/kg]GotoStolwijkGordonGotoStolwijkGordonHeadBrain904.0750.0814.0602.0419.0579.0Skull371.00.00.0194.00.00.0Muscle23.22.010.030.05.431.5Fat46.52.20.179.95.90.3Skin20.324.096.279.988.9509.0Segmental sum1370.0778.0920.0272.0194.0224.0NeckCore4.5–0.027.7–0.0Muscle10.6–18.230.0–31.3Fat4.5–0.037.1–0.3Skin1.7–5.637.1–174.0Segmental sum21.3–23.831.2–25.9Chest and abdomenCore2070.03500.02830.0127.0238.0155.0Muscle366.0100.0377.030.05.631.4Fat75.142.72.416.96.00.3Skin21.435.094.917.425.977.4Segmental sum2540.03680.03300.074.189.781.7Upper arm and lower arm (right or left)Core1.97.00.03.36.30.0Muscle104.09.577.560.65.631.4Fat13.81.70.123.03.40.3Skin6.04.223.423.017.483.9Segmental sum125.022.3101.040.16.326.2Hand (right or left)Core0.40.80.03.36.40.0Muscle16.22.04.4122.057.131.4Fat13.00.30.0111.04.40.3Skin8.616.718.5111.0175.0201.0Segmental sum38.119.822.985.159.241.9Thigh and leg (right or left)Core4.222.40.03.36.40.0Muscle254.028.6115.051.15.631.4Fat41.74.30.423.33.60.3Skin10.623.835.023.239.671.1Segmental sum311.079.1151.036.67.619.0Foot (right or left)Core1.51.30.03.36.20.0Muscle12.70.26.851.14.831.3Fat10.40.40.136.33.80.3Skin4.725.015.736.3208.0123.0Segmental sum29.226.922.526.456.120.4Total4930.04750.04840.074.463.866.9

Despite its advantages, this study has some limitations. First, the model was developed for an average Japanese male in their 20s, and it cannot predict physiological quantities for individuals with different characteristics, such as size, age, sex, heredity, or pre-existing disorders. To resolve this issue, it is necessary to adjust the model parameter settings to suit individual characteristics, and further studies are required. Second, the distribution of compartment BFRs in each segment was set based on experimental measurements to the extent possible; however, certain assumptions were still necessary. Therefore, there is still room for consideration regarding the validity of these settings. Finally, this study was limited to the development of the physical model that could simulate BT, BP, and BFR only under the baseline condition. To simulate various conditions, an autonomic regulation model must be developed and coupled with the physical model. We plan to report on this issue in our next paper.

Conclusions

In this study, a new physical model that couples a TNW model for simulating BT with a CV model for simulating BFR and BP was developed. An average Japanese male in their 20s was set as the standard person, and the parameter settings for the person were derived based on recent experimental data and the literature. By applying the parameter settings to the model, it was confirmed that the simulated BT, BFR, and BP values under a thermally neutral and supine condition reproduced the experimental values well.

Supplementary Information

Below is the link to the electronic supplementary material.

Supplementary Material 1 (DOCX 192 KB)