# Modeling of bio-heat transfers in lungs with fractional models

Enso Ndreko, Stéphane Victor, Jean-François Duhé, Pierre Melchior

► To cite this version:

Enso Ndreko, Stéphane Victor, Jean-François Duhé, Pierre Melchior. Modeling of bio-heat transfers in lungs with fractional models. Annual Reviews in Control, 2025, 60, pp.101010. 10.1016/j.arcontrol.2025.101010. hal-05311650

HAL Id: hal-05311650

https://hal.science/hal-05311650v1

Submitted on 13 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modeling of bio-heat transfers in lungs with fractional models

Enso Ndreko<sup>b</sup>, Stéphane Victor<sup>b,D</sup>,*, Jean-François Duhé<sup>a,D</sup>, Pierre Melchior<sup>b,D</sup>

## ARTICLE INFO

**Keywords:**

Lung thermal modeling

Bio-heat equation

Heat transfer

Thermal systems

Biological systems

Fractional calculus

Fractional systems

Model order reduction

Cardiopulmonary bypass

## ABSTRACT

In cardiac surgeries, when cardiopulmonary bypass (CPB) (or extracorporeal circulation (ECC)) is employed, the lungs are temporarily disconnected from the body. To minimize the risk of tissue damage or respiratory complications, the lungs are subjected to mild hypothermia. Incorporating dynamic heat transfer modeling offers the potential to enhance temperature regulation through a more advanced approach.

A complex thermal model, based on a thermal two-port network, offers a wide frequency range applicability, making it suitable for modeling the human breathing frequencies. This modeling approach can also be adapted to incorporate the influence of blood flow, which serves as a natural temperature regulator in the human body. This is accomplished by combining the thermal two-port network with the bio-heat equation.

The main contributions focus on introducing distinctive and simplified approximation models for the equivalent global impedance of thermal transfer within the lungs. These models, featuring minimal parameters, manifest comparable dynamic traits in the frequency domain, akin to the attributes of the two-port network model. This progress clears the way for broader utilization across various domains.

## 1. Introduction

Thermal modeling is essential in various fields where temperature plays a critical role. It finds application in areas such as air conditioning, industrial refrigeration, cooling of electronic devices (Künzi, 2015), building simulation (Danza et al., 2016; Parnis, 2012) and heat loss in biological tissues (Ismail et al., 2018; Jiang et al., 2004). The utilization of fractional calculus in bioengineering has revolutionized the study of diffusive phenomena and intricate geometries encountered in biology (Magin, 2004).

During cardiac surgeries, organs being exposed to the outside environment result on a temperature decrease. In fact, during normal conditions, the internal temperature of the human body is approximately 37 °C, and every minor shifts from this value can damage the biological tissues. Surgeons remove the lungs from the body of the patient during extracorporeal circulation (ECC) or cardiopulmonary bypass (CPB). Therefore, in order to lower the rate of tissue damage, mild hypothermia is employed (Giesinger & McNamara, 2018) by slowing down metabolism and O₂ consumption. However, this is frequently accomplished using empirical techniques. The integration of dynamic heat transfer models has the capacity to provide a more advanced strategy for temperature regulation.

Hildebrandt's works revolve around the modeling of lung elasticity, which he explored through experiments involving both balloons

and cats lungs (Hildebrandt, 1969, 1970). Within such modeling, the presence of fractional order models was indirectly indicated. In addition to Hildebrandt's contributions, there have been developments in constructing models that offer a simplified mechanical impedance equivalent for the entire pulmonary system (Diong et al., 2007; DuBois et al., 1956; Navajas et al., 1990).

The fractalist, Mandelbrot, examined the arborescent structure and recursive nature of the lungs (Weibel, 2005). Taking into account the geometry of the lungs, Ionescu and Oustaloup, suggested an equivalent circuit model (Ionescu et al., 2011). In their work, the model proved to be complex and therefore difficult to be implemented. To make up for this, Ionescu suggested a 4th-order fractional impedance model, suitable for analyzing experimental data (Ionescu et al., 2012).

A thermal model with fractional order characteristics is also documented in the literature (Pellet et al., 2011; Victor et al., 2020) where Havriliak-Negami transfer function is proposed as a black-box model and offline system identification. Inspired by Maillet et al. (2000), a thermal two-port network model (Fig. 2), with two series impedances Z₁ and Z₂ and a shunt impedance Z₃, is used and studied by J.-F. Duhé as a thermal modeling of a single lung branch (Duhé et al., 2022). These impedance expressions are complex and nonlinear in frequency which do not allow obtaining rational transfer functions.

<footnote><sup>a</sup> Universidad de Panamá, FIEC, Panamá, Panama</footnote>

<footnote><sup>b</sup> Univ. Bordeaux, CNRS, IMS UMR 5218, Bordeaux INP/ENSEIRB-MATMECA, 351 Cours de la Libération, 33405 Talence CEDEX, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: enso.ndreko@ims-bordeaux.fr (E. Ndreko), stephane.victor@ims-bordeaux.fr (S. Victor), jean-f.duhe-p@up.ac.pa (J.-F. Duhé), pierre.melchior@ims-bordeaux.fr (P. Melchior).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. 1D planar thermal wall system with section $S_w$, length $L$, thermal conductivity $k$, density $\rho$ and specific heat capacity $c$.</figcaption></figure>

When dealing with living tissues, the presence of blood circulation enables nutrient and gas exchange, which is essential for organic functionality, and at the same time, acts as a natural thermal regulator. Pennes (1948) came up with a bio-heat equation to take into account such phenomena by adding the blood perfusion term, the metabolic heat generation and an external heat source.

A two-port network model, which takes into account the blood perfusion term, has been studied and used in single thermal lung branch modeling in Duhé et al. (2021).

A global lung circuit model, derived from the heat equation, has been proposed in Duhé et al. (2023) (Fig. 4). Even though the complexity of this model is considerable, since it is derived from the heat-equation and takes into account the perfusion term, it is in fact a rough approximation with a very simple transfer function of the first kind.

The contributions of this paper are on proposing distinct and simplified approximation models for the equivalent global impedance of heat transfers in the lungs. Therefore, this will open the way for further application in system identification (see e.g. Victor et al., 2025), such as linear parameter varying system identification or real-time system identification.

The structure of this paper is as follows: Section 2 unveils the thermal two-port network by using the classic heat equation and the bio-heat one. Section 3 will study the equivalent global circuit model of the lung, and its analysis in the frequency domain. Then, Section 4 will propose several functional approximations in the frequency domain. Ultimately, Section 5 will provide conclusions, concluding remarks, and avenues for further research.

## 2. Two-port network formalism for a single lung branch level

A heat conduction on a simple plane wall, with length $L$ and the cross section $S_w$, is shown in its longitudinal direction $x$ in Fig. 1. The thermal notations $k$, $\rho$ and $c$ stand for the medium thermal conductivity, density and heat capacity respectively.

### 2.1. Classic two-port network

The classic case of the heat conduction is derived from the original heat equation in a homogeneous medium (a lung branch level in this case):

$$ \rho_a c_a \frac{\partial T}{\partial t} = k_a \nabla^2 T, $$

where $T$ defines the temperature and as conduction only goes on $x$-axis, the 1D heat equation can be written as:

$$ \rho_a c_a \frac{\partial T}{\partial t} = k_a \frac{\partial^2 T}{\partial x^2}. \qquad (1) $$

Applying the Laplace transform on (1), with null initial conditions, leads to:

$$ sT(x, s) = \frac{k_a}{\rho_a c_a} \frac{\partial^2 T(x, s)}{\partial x^2}, \qquad (2) $$

where $s$ denotes the Laplace variable. The solving of the general partial differential equation is expressed by

$$ T(x, s) = c_1 \cosh \left( \sqrt{\frac{s \rho_a c_a}{k_a}} x \right) + c_2 \sinh \left( \sqrt{\frac{s \rho_a c_a}{k_a}} x \right) \qquad (3) $$

where $c_1$ and $c_2$ are determined according to the following boundary conditions, namely by considering a heat flux input at $x = 0$ and a heat flux output at $x = L$, it then comes:

$$ \begin{aligned} \phi_{in}(s) &= -k_S \frac{\partial T(x, s)}{\partial x} \Big|_{x=0} \\ \phi_{out}(s) &= -k_S \frac{\partial T(x, s)}{\partial x} \Big|_{x=L}, \\ T_{in}(s) &= T(x, s) \Big|_{x=0} \\ T_{out}(s) &= T(x, s) \Big|_{x=L} \end{aligned} \qquad (4) $$

or even, put under a matrix form (Maillet et al., 2000):

$$ \begin{bmatrix} T_{in}(s) \\ \phi_{in}(s) \end{bmatrix} = M \begin{bmatrix} T_{out}(s) \\ \phi_{out}(s) \end{bmatrix}, \qquad (5) $$

where

$$ M = \begin{bmatrix} \cosh(\delta L) & \frac{1}{k_a S_w \delta} \sinh(\delta L) \\ k_a S_w \delta \sinh(\delta L) & \cosh(\delta L) \end{bmatrix}, \qquad (6) $$

with

$$ \delta = \sqrt{a_a^s}, \qquad (7) $$

and thermal diffusivity $a_a = \frac{k_a}{\rho_a c_a}$.

Note that $\delta$ involves the presence of a fractional (non-integer) operator, as its time-domain representation is a half-order derivative.

### 2.2. Bio-heat two-port network

In living tissue, blood circulation plays a dual role. It allows vital nutrient and gas exchanges which are crucial for the organic function while serving as a natural temperature regulator. Pennes (1948) introduced a bio-heat equation to incorporate these processes. This was achieved by incorporating the perfusion term $\rho_b c_b \omega_b (T_a - T)$ into the heat equation ($\rho_b$ being the blood density and $c_b$ the blood heat capacity), along with metabolic heat generation $q_m$ and an external heat source $q_{ext}$.

Assuming no external heat source $q_{ext} = 0$ and that the effect of the metabolic heat is negligible, the bio-heat equation is expressed by:

$$ \rho_t c_t \frac{\partial T}{\partial t} = k_t \nabla^2 T + \rho_b c_b \omega_b (T_a - T), \qquad (8) $$

where $T_a$ is the arterial temperature, $\rho_t$ the tissue density, $k_t$ the tissue thermal conductivity and $c_t$ the tissue specific heat capacity. Or else, put under a one dimensional conduction, it then comes:

$$ \rho_t c_t \frac{\partial T}{\partial t} = k_t \frac{\partial^2 T}{\partial x^2} + \rho_b c_b \omega_b (T_a - T). \qquad (9) $$

It is assumed that the temperature of an internal organ remains constant at a core temperature $T_b$, which is close to the normal human body temperature of 37 °C. By introducing temperature variations of a tissue around this core temperature, denoted by $\Delta T = T - T_b$, the time and space partial derivatives for $\Delta T$ and $T$ are equivalent since $T_b$ is assumed to be constant. Additionally, under normal conditions, the arterial blood temperature $T_a$ approaches the core temperature $T_b$ ($T_a \approx T_b$). Therefore, (9) boils down to:

$$ \rho_t c_t \frac{\partial \Delta T}{\partial t} = k_t \frac{\partial^2 \Delta T}{\partial x^2} - \rho_b c_b \omega_b \Delta T. \qquad (10) $$

Denoting $a_t$ as the tissue thermal diffusivity and $h_b$ as the blood perfusion heat coefficient, they can be expressed by:

$$ a_t = \frac{k_t}{\rho_t c_t} \quad \text{and} \quad h_b = \frac{\rho_b c_b \omega_b}{k}, \qquad (11) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Thermal two-port network.</figcaption></figure>

and substituting them into (10), one gets the following simplified equation:

$$ \frac{1}{a_t} \frac{\partial \Delta T}{\partial t} = \frac{\partial^2 \Delta T}{\partial x^2} - h_b \Delta T. \qquad (12) $$

Again, applying the Laplace transform while assuming null initial conditions, one gets:

$$ \left[ \frac{s}{a_t} + h_b \right] \Delta T(x, s) = \frac{\partial^2 \Delta T(x, s)}{\partial x^2}. \qquad (13) $$

This results in a nearly identical two-port network configuration (5), by replacing $\delta$ with the following bio-heat one, $\delta_{bio}$:

$$ \delta_{bio} = \sqrt{\frac{s}{a_t} + h_b}. \qquad (14) $$

## 2.3. Equivalent circuit representation

Matrix M in both previous cases, gives rise to an electrical analogy as an equivalent simple circuit, based on the two-port network theory (Fig. 2).

The change from the matrix formulation (5) into an equivalent T circuit is described in Appendix A.

Indeed the thermal two-port network formalism illustrates two series impedances, $Z_1(s)$ and $Z_2(s)$, along with a shunt impedance, $Z_3(s)$ that are expressed by:

$$ \begin{cases} Z_1(s) = Z_2(s) = \frac{1}{kS_w\delta}[\coth(\delta L) - \operatorname{csch}(\delta L)] \\ Z_3(s) = \frac{1}{kS_w\delta}\operatorname{csch}(\delta L), \end{cases} \qquad (15) $$

where $\delta$ will be taken as $\sqrt{\frac{s}{a}}$ and $k = k_a$ for the classic heat conduction, and as $\sqrt{\frac{s}{a} + h_b}$ and $k = k_t$ for the bio-heat one.

# 3. Lung heat transfer model

## 3.1. Lung global heat transfer model

A lung representation is shown in Fig. 3, which highlights the temperature at each bifurcation point, as well as the corresponding length $L_n$ and cross-section $S_n$ at each branch level. Note that the branch level $n = 0$ corresponds to the trachea.

By using the two-port network formalism, as described in the previous Section 2, the global lung can be illustrated as an equivalent circuit model (Duhé et al., 2023), see Fig. 4.

Murray's law allows to establish a relationship between a main vessel that is bifurcated into two smaller vessels. If a main vessel has diameter $D_1$ and the two following vessels at the bifurcation have diameters $D_2, D_3$, and the work done by the fluid flow is to be minimized, the following relationship holds:

$$ D_1^3 = D_2^3 + D_3^3. $$

<figure><img src="image_4.png" /><figcaption>Fig. 3. Lung scheme.</figcaption></figure>

When both *children* vessels have the same diameters:

$$ D_1^3 = 2D_{children}^3, $$

the following relationship is obtained:

$$ \frac{D_{children}}{D_1} = 2^{-1/3}. $$

Consequently, by following Murray's law, Kuwahara et al. (2009) have shown that when considering two consecutive branches, the lengths and cross-sections can be related as follows:

$$ \frac{S_{n+1}}{S_n} = 2^{-\frac{2}{3}} \quad \text{and} \quad \frac{L_{n+1}}{L_n} = 2^{-\frac{1}{3}}. \qquad (16) $$

In human lungs, the maximal number of branch levels is $n = 23$, where the trachea level corresponds to $n = 0$ (Kuwahara et al., 2009). In practice, gas exchanges and the effect of blood perfusion starts from level $n = 14$. Therefore, in order to develop a comprehensive model of the lungs, a combination of the classic thermal two-port network and the bio-heat equation is suggested to account for these variations:

$$ \begin{cases} n < 14 & \text{classic} \\ n \ge 14 & \text{bio-heat}. \end{cases} \qquad (17) $$

Consequently, the impedances $Z_{1,n}$ (and therefore $Z_{2,n}$) and $Z_{3,n}$ at each level can be expressed as:

$$ \begin{cases} Z_{1,n}(s) = Z_{2,n}(s) \\ = \frac{1}{kS_02^{-\frac{n}{3}}\delta} \left[ \coth\left(\delta L_0 2^{-\frac{n}{3}}\right) - \operatorname{csch}\left(\delta L_0 2^{-\frac{n}{3}}\right) \right], \\ Z_{3,n}(s) = \frac{1}{kS_02^{-\frac{n}{3}}\delta} \operatorname{csch}(\delta L_0 2^{-\frac{n}{3}}). \end{cases} \qquad (18) $$

where $\delta$ is defined as in (7) and $k = k_a$ for the classic heat conduction case, or as in (14) and $k = k_t$ for the bio-heat one, where $h_b$ is the perfusion heat coefficient (Duhé et al., 2021).

A frequency domain analysis of the lung impedances is carried out from the trachea ($n = 0$) to a branch level. The physiological parameters of a human being for air, tissue, and blood are defined in Table 1.

## 3.2. Frequency domain analysis of $Z_{1,n}$ and $Z_{3,n}$

Fig. 5 illustrates the Bode gain diagrams of $Z_{1,n}$ and $Z_{3,n}$ for different lung level $n$.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Lung global equivalent circuit model.</figcaption></figure>

<table><caption>**Table 1**
Nominal simulation parameters for human beings (Coccarelli et al., 2016).</caption><thead><tr><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>k<sub>a</sub></td><td>Air thermal conductivity</td><td>0.025 W m<sup>-1</sup>K<sup>-1</sup></td></tr><tr><td>a<sub>a</sub></td><td>Air thermal diffusivity</td><td>2.2 &times; 10<sup>-5</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>k<sub>t</sub></td><td>Tissue thermal conductivity</td><td>0.28 W m<sup>-1</sup>K<sup>-1</sup></td></tr><tr><td>&rho;<sub>t</sub></td><td>Tissue density</td><td>550 kg m<sup>-3</sup></td></tr><tr><td>c<sub>t</sub></td><td>Tissue specific heat capacity</td><td>3718 J kg<sup>-1</sup>K<sup>-1</sup></td></tr><tr><td>a<sub>t</sub></td><td>Tissue thermal diffusivity</td><td>1.37 &times; 10<sup>-7</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>&rho;<sub>b</sub></td><td>Blood density</td><td>1060 kg m<sup>-3</sup></td></tr><tr><td>c<sub>b</sub></td><td>Blood specific heat capacity</td><td>3900 J kg<sup>-1</sup>K<sup>-1</sup></td></tr><tr><td>&omega;<sub>b</sub></td><td>Blood perfusion coefficient</td><td>0.0043 s<sup>-1</sup></td></tr><tr><td>L<sub>0</sub></td><td>Trachea length</td><td>0.1 m</td></tr><tr><td>d<sub>0</sub></td><td>Trachea diameter</td><td>1.4 cm</td></tr></tbody></table>

First, focusing on $Z_{1,n}$, it is interesting to note that it behaves like a low-pass filter with a bandwidth inside $[10^{-1}, 10^{1}]$ rad/s. Moreover, the higher the level $n$, the higher the static gain for all cases ($\delta$ for $n < 14$ and $\delta_{bio}$ otherwise), which directly comes from Murray's law (16). All frequency responses have a slope of $-20$ dB/dec in high frequency, which gives a fractional slope of $-10$ dB/dec. Consequently, there is a phase lock of $-45^{\circ}$ on the Bode phase diagram. These characteristics are natural when considering constant-phase elements, and also, as the expressions in (18) are given with the square root of the Laplace variable $s$.

The behavior of the shunt impedance $Z_{3,n}$ is more particular as the gain (and consequently the phase) goes down to $-\infty$ in high frequency. Again, the higher the level $n$, the higher the static gain for all cases ($\delta$ for $n < 14$ and $\delta_{bio}$ otherwise), which directly comes from Murray's law (16). For the first levels of $n < 14$, the gain diagrams of $Z_{3,n}$ appears to have an integrator behavior of $-20$ dB/dec in low frequency. For $n \ge 14$, $Z_{3,n}$ have a static gain that increases as the level $n$ increases.

### 3.3. Equivalent impedance from a branch level $k$ to a deeper branch level $n$

To analyze the lung equivalent impedance from a branch level $k$ to a deeper level $n$, denoted as $Z_{k-n}^{eq}$, we refer to an equivalent circuit model

<figure><img src="image_6.png" /><figcaption>Fig. 5. Bode diagram of $Z_1$ (top) and $Z_3$ (bottom) for different levels (the plain curves are plotted for $\delta$ and the dashed curves for $\delta_{bio}$).</figcaption></figure>

illustrated in Fig. 4. Given the expressions for $Z_{1,n}$ and $Z_{3,n}$ provided in (18), the equivalent impedance $Z_{k-n}^{eq}$ can be recursively expressed by:

$$Z_{k-n}^{eq} = Z_{1,k} + \frac{Z_{3,k} \left( Z_{(k+1)-n}^{*} + Z_{2,k} \right)}{Z_{3,k} + Z_{(k+1)-n}^{*} + Z_{2,k}}, \quad \text{for} \quad \begin{cases} k = 0, 1, \dots, n-1 \\ n = 1, 2, \dots, 23 \\ k \le n \end{cases} \qquad (19)$$

where $Z_{(k+1)-n}^{*} = \frac{Z_{(k+1)-n}^{eq}}{2}$, and $Z_{(k+1)-n}^{eq}$ is the equivalent impedance from a level $k+1$ to a deeper level $n$. It is assumed that the temperature and heat flux in a branch bifurcation of the same level are identical.
### 3.4. Frequency analysis of the equivalent impedance from trachea (level 0) to a branch level n

We are indeed interested on the equivalent impedance $Z_{0-n}^{eq}$, namely the equivalent impedance from trachea to level n.

Given the expressions of $Z_{1,n}$ and $Z_{3,n}$ provided in (18), the equivalent impedance $Z_{0-n}^{eq}$ from the trachea up to a branch level n, can be expressed by:

$$Z_{0-n}^{eq} = Z_{1,0} + \frac{Z_{3,0}(Z_{1-n}^{*} + Z_{2,0})}{Z_{3,0} + Z_{1-n}^{*} + Z_{2,0}}, \quad (20)$$

where $Z_{1-n}^{*}$ is defined by:

$$Z_{1-n}^{*} = \frac{(Z_{1-n}^{eq})^2}{2Z_{1-n}^{eq}} = \frac{Z_{1-n}^{eq}}{2}, \quad (21)$$

and $Z_{1-n}^{eq}$ is the equivalent impedance of one of the branches from level 1 up to branch level n.

Notice that the equivalent impedance from level 1 to level n, denoted $Z_{1-n}^{eq}$, appears in parallel with itself. This is because the lung branches are assumed to be symmetric, implying that the temperature profiles in branches at the same level are identical. As a result, the impedance $Z_{1-n}^{*}$ in Eq. (21) is calculated as the parallel combination of two identical impedances:

$$Z_{1-n}^{*} = \frac{Z_{1-n}^{eq}}{2}.$$

The same reasoning is applied to the other branches, namely $Z_{2-n}^{eq}$, $Z_{3-n}^{eq}$, and so forth as symmetrical branches are parallel to themselves, except for the trachea (branch level 0).

Following this, $Z_{1-n}^{*}$ is placed in series with $Z_{2,0}$. This combination is then in parallel with $Z_{3,0}$. Finally, the result of that parallel combination is connected in series with $Z_{1,0}$, forming the total equivalent impedance $Z_{0-n}^{eq}$. When it comes on calculating $Z_{1-n}^{eq}$, the same idea can be applied:

$$Z_{1-n}^{eq} = Z_{1,1} + \frac{Z_{3,1}(Z_{2-n}^{*} + Z_{2,1})}{Z_{3,1} + Z_{2-n}^{*} + Z_{2,1}}, \quad (22)$$

and so on for $Z_{2-n}^{eq}$, $Z_{3-n}^{eq}$, etc.

The Bode diagrams of the equivalent impedance of different lung levels are shown in Fig. 6, which corresponds to the equivalent impedance $Z_{0-2}^{eq}$ up to $Z_{0-23}^{eq}$ with a level step of three.

Note that the gain Bode diagrams of $Z_{0-14}^{eq}$ up to $Z_{0-23}^{eq}$, which correspond to the bio-heat equation case, are plotted in Fig. 6. The gain diagrams exhibit a low-pass filter behavior, which are very similar to that observed in gain Bode diagrams of $Z_{1,n}$. Moreover, the gain diagrams are exactly the same.

Note how the classic heat conduction model (from impedance $Z_{0-1}^{eq}$ to $Z_{0-13}^{eq}$) behaves as a low-frequency capacitance, allowing sufficiently slow thermal disturbances to propagate fully through the lungs. For the deeper levels corresponding with the bio-heat transfers (from level [0 - 14] to [0 - 23]), a constant static gain appears in low frequency. The high-frequency behavior of all presented equivalent impedances are similar. For these cases, note that in high frequency, they behave like a constant-phase element with a phase lock of -45° and a slope of -10 dB/dec, very normal when we recall that the square root is present in the Laplace operator s. Since the equivalent impedances from $Z_{0-14}^{eq}$ up to $Z_{0-23}^{eq}$ are the same, from now on, model $Z_{0-14}^{eq}$ will be used instead of model $Z_{0-23}^{eq}$.

### 4. Model order reduction of the global lung impedance

The analytical expression of $Z_{0-14}^{eq}$ is defined with all intermediate expressions of $Z_{1,n}$ and $Z_{3,n}$ provided in (18). Consequently, in an objective to detect any variations of the global lung model, the analytical expressions of $Z_{0-14}^{eq}$ are too complex to be used in real-time. Therefore, given the frequency responses of $Z_{0-14}^{eq}$ and especially in an open-heart

<figure><img src="image_7.png" /><figcaption>Fig. 6. Bode diagram of the equivalent impedance $Z_{0-n}^{eq}$ from the trachea (level 0) up to a branch level.</figcaption></figure>

surgery with CPB, approximated compact models that are usable in real-time is sought in the following.

All calculations and simulations for fractional models have been carried out with the CRONE Toolbox developed in Matlab®. The CRONE toolbox has evolved into an object-oriented programming (Victor & Malti, 2023) and allows many enhancements such as overloading operators and scripts. It is now freely available for research and educational purposes: https://cronetoolbox.ims-bordeaux.fr.

Consider the Mean Square Error (MSE) criterion, involving the gain errors in dB, defined as:

$$MSE = \frac{1}{N} \sum_{i=1}^{N} e_i^2 = \frac{1}{N} e^T e, \quad (23)$$

where $e_i = |Z_{i(0)}| - |\hat{Z}_{i(0)}|$ is in the frequency range [10⁻⁷, 1] rad/s.

Even though this study focuses on gain optimization only, the criterion MSE may be modified in order to optimize the full complex values of the frequency response, which inherently includes gain and phase.

Under normal conditions, human breathing rates can approximately vary from 12 to 20 breaths per minute. This would mean a fluctuation between 1.26 rad/s to 2.09 rad/s. However, under surgery conditions (in cardiopulmonary bypass), breathing can go as low as 5 breaths per minute, which would be around 0.52 rad/s. Rogers et al. (2024) started by a 12 breaths per minute rate and during CPB it was decreased to 5, which justifies performing the optimization up to around 1 rad/s.

The minimization of the MSE is an optimization problem where the number of parameters may lead to local minima. To overcome this, the meta-heuristic optimization PSO (Particle Swarm Optimization) algorithm (Kennedy & Eberhart, 1995), inspired by swarm behavior observed in nature such as fish and bird schooling, is chosen. The PSO algorithm offers several advantages, including straightforward concepts, easy implementation, applicability to various problem types, efficiency in multi variable optimization, effectiveness in nonlinear, nonconvex environments and search spaces. The optimization process begins by generating a set of random candidate solutions within the search space boundaries. Through iterative steps, the algorithm improves these candidate solutions. Once the iterations are completed, the best candidate solution is identified and presented as the solution to the problem.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Gain diagram of $Z_{0-14}^{eq}$ and $Z_{asymp}$.</figcaption></figure>

<table><caption>**Table 2**
Estimated parameters of $Z_{asymp}$ and the resulting MSE.</caption><thead><tr><th>K</th><th>τ</th><th>MSE</th></tr></thead><tbody><tr><td>70 152.5115</td><td>57.5585</td><td>3.636</td></tr></tbody></table>

<table><caption>**Table 3**
Estimated parameters of $Z_{F-PZ1}$ and resulting MSE.</caption><thead><tr><th>K</th><th>τ<sub>1</sub></th><th>ν<sub>1</sub></th><th>τ<sub>2</sub></th><th>ν<sub>2</sub></th><th>MSE</th></tr></thead><tbody><tr><td>56 710.5385</td><td>338.5224</td><td>0.7184</td><td>6.1659</td><td>0.2540</td><td>1.615 × 10<sup>-1</sup></td></tr></tbody></table>

Other algorithms such as Flower Pollination Algorithm (FPA), Grey Wolf Optimizer (GWO), Differential Evolution (DE) and Genetic Algorithm (GA) were used as well. Nevertheless, during this work, it was noticed that PSO is the most effective when it comes to minimizing the MSE criterion. Indeed this is especially noticed in Section 4.4, when a large number of parameters is required.

In this paper, an initialization with PSO algorithm is done, and then, the optimization function fmincon is applied by using interior-point methods under Matlab®.

### 4.1. Asymptotic approximation

$Z_{0-14}^{eq}$ is a constant-phase element with a phase lock of $-45°$ and with a slope of $-10$ dB/dec when $\omega \to \infty$. As a result of that, a fractional integrator of order 0.5 should be present in high frequency.

An asymptotic model that corresponds to a first kind (which corresponds to the generalization of a first order system to a fractional one with a single $s^\nu$-pole), needs to fit the analytical global model $Z_{0-14}^{eq}$ in low and high frequency, and consequently is expressed as follows:

$$Z_{asymp}(s) = \frac{K}{1 + \tau s^{0.5}} \qquad (24)$$

The optimization provides the estimation of 2 parameters which are gathered in Table 2.

Fig. 7 shows both the gain Bode diagram of $Z_{0-14}^{eq}$ and of its asymptotic approximation $Z_{asymp}$ (in dashed red). As expected the middle frequencies are not well estimated, leading to a MSE of 3.636 that is quite high.

### 4.2. A single fractional pole-zero approximation

Inspired by the previous model (24), as the middle frequencies are not well estimated, a more complex model is sought with a single

<figure><img src="image_9.png" /><figcaption>Fig. 8. Gain diagram of $Z_{0-14}^{eq}$ and $Z_{F-PZ1}$.</figcaption></figure>

<table><caption>**Table 4**
Estimated parameters of $Z_{BW}$ approximation and the resulting MSE.</caption><thead><tr><th>a</th><th>b</th><th>c</th><th>d</th><th>α</th><th>MSE</th></tr></thead><tbody><tr><td>-454.5100</td><td>662.7800</td><td>-658.0100</td><td>552 459.0900</td><td>0.4985</td><td>9.641 × 10<sup>-1</sup></td></tr></tbody></table>

fractional cell of pole and zero approximation. The proposed model gets the following expression:

$$Z_{F-PZ1}(s) = K \frac{1 + \tau_2 s^{\nu_2}}{1 + \tau_1 s^{\nu_1}} \qquad (25)$$

The optimization provides the estimation of 5 parameters which are gathered in Table 3.

Fig. 8 shows both the gain Bode diagram of $Z_{0-14}^{eq}$ and of its fractional pole and zero approximation $Z_{F-PZ1}$ (in dotted yellow). As expected, the middle frequencies are better estimated leading to a MSE of 0.162 that is more than twenty times smaller than the asymptotic approximation.

Despite such a precision, this model does not satisfy the integral order in high frequency that should be of 0.5: the optimization has led to $\nu_1 - \nu_2 = 0.4644$. Even though being very close to the true order of 0.5, it is still not satisfactory.

### 4.3. Fractional Butterworth approximation

Inspired by Duhé et al. (2022), we can use another interesting approximation. The fractional Butterworth filter (Ali et al., 2013) is an approximation that is more compact:

$$Z_{BW}(s) = \frac{d}{s^{\alpha+\beta} + as^{\alpha} + bs^{\beta} + c} \qquad (26)$$

Given the high frequency behavior of $Z_{0-14}^{eq}$, $Z_{BW}(j\omega)$ should have:

$$\lim_{\omega \to \infty} Z_{BW}(j\omega) = \frac{d}{(j\omega)^{\alpha+\beta}} = \frac{d}{(j\omega)^{0.5}} \qquad (27)$$

with $\alpha, \beta \in [0, 2)$ for stability reasons. From this, by substituting $\beta$ with $0.5 - \alpha$, we lower the number of parameters by 1. Therefore we have:

$$Z_{BW} = \frac{d}{0.5 + as^{\alpha} + bs^{0.5-\alpha} + c} \qquad (28)$$

The optimization provides the estimation of 5 parameters which are gathered in Table 4.

Fig. 9 shows both the gain Bode diagram of $Z_{0-14}^{eq}$ and of its Butterworth approximation $Z_{BW}$ (in dashed-dotted purple). As expected, the middle frequencies are better estimated than the asymptotic approximation, however some oscillations appear in low frequency and the
<table><caption>**Table 5**
Estimated parameters of $Z_{\text{multi-PZ}}$ approximation and the resulting MSE.</caption><thead><tr><th>τ<sub>p1</sub></th><th>τ<sub>z1</sub></th><th>τ<sub>p2</sub></th><th>τ<sub>z2</sub></th><th>τ<sub>p3</sub></th><th>τ<sub>z3</sub></th><th>τ<sub>p4</sub></th><th>τ<sub>z4</sub></th></tr></thead><tbody><tr><td>3.3792</td><td>3.2540</td><td>94.4724</td><td>82.1647</td><td>1205.5523</td><td>847.3524</td><td>5948.0971</td><td>8583.1405</td></tr><tr><th>τ<sub>p5</sub></th><th>τ<sub>z5</sub></th><th>τ<sub>p6</sub></th><th>τ<sub>z6</sub></th><th>τ<sub>p7</sub></th><th>τ<sub>z7</sub></th><th>K</th><th>τ</th><th>MSE</th></tr><tr><td>46 767.1375</td><td>56 096.4034</td><td>327 738.1962</td><td>352 012.9613</td><td>4 917 213.9878</td><td>5 080 257.2074</td><td>70 470.7281</td><td>65.2489</td><td>9.33 × 10<sup>-4</sup></td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Gain diagram of $Z_{0-14}^{eq}$ and $Z_{BW}$.</figcaption></figure>

integral order is better in high frequency as compared to the fractional zero-pole approximation.

The optimization has led to a MSE of 0.964, that is almost four times smaller than the asymptotic approximation, but less precise than the single fractional pole-zero approximation, whilst having the same number of parameters.

### 4.4. Multiple pole-zero approximation

Inspired by the first reduced-order model (24), as the middle frequencies are not well estimated, a more complex model is sought by adding poles and zeros. The proposed model has the following expression:

$$Z_{\text{multi-PZ}} = \frac{K}{1 + \tau s^{0.5}} \prod_{i=1}^{N_{\text{cells}}} \frac{1 + \tau z_i s}{1 + \tau p_i s} \quad (29)$$

In order to determine the optimal number of cells $N_{\text{cells}}$, minimizing the MSE would lead to over-parameterize the system model. To limit to model complexity, let us use the Akaike Information Criterion (AIC) (Akaike, 1974; Ljung, 1999). Indeed the AIC is generally used in model order selection and will help fitting the goodness-of-fit of the model, while penalizing the model for over-parameterizing the number of parameter.

The AIC criterion is defined as (Ljung, 1999):

$$\text{AIC} = N \ln \left( \frac{1}{N} \sum_{i=1}^{N} e_i^2 \right) + 2n_p + N(\log(2\pi) + 1), \quad (30)$$

where N is the number of samples in the estimation data set, $e_i$ is the prediction errors, and $n_p$ is the number of parameters.

Fig. 10 illustrates the evolution of the AIC criterion according to the number of parameters, which clearly shows that the optimal number of pole-zero cells is $N_{\text{cells}} = 7$.

Based on this value, the optimization provides the estimation of 16 parameters which are gathered in Table 5.

Fig. 11 shows both the gain Bode diagram of $Z_{0-14}^{eq}$ and of its multiple pole and zero approximation $Z_{\text{multi-PZ}}$ (in diamond cyan). As

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** AIC versus the number of cell approximation $N_{\text{cells}}$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Gain diagram of $Z_{0-14}^{eq}$ and $Z_{\text{multi-PZ}}$.</figcaption></figure>

<table><caption>**Table 6**
Estimated parameters of $Z_{\text{imp}}$ approximation and the resulting MSE.</caption><thead><tr><th>K</th><th>τ<sub>1</sub></th><th>ν<sub>1</sub></th><th>τ<sub>2</sub></th><th>ν<sub>2</sub></th><th>MSE</th></tr></thead><tbody><tr><td>70 079.9721</td><td>2413.4335</td><td>1.0156</td><td>1406.6076</td><td>0.5156</td><td>1.096 × 10<sup>-3</sup></td></tr></tbody></table>

expected, the middle frequencies are extremely well estimated, leading to a MSE of $9.33 \times 10^{-4}$ that is almost four thousand times smaller as compared with the asymptotic model, but with a significative increase of parameter number.

Given the time constants of the pole and zero cells provided in Table 5, these frequencies are in the frequency range $[10^{-7}, 10^{-1}]$, which is exactly the range in which the asymptotic approximation $Z_{\text{asymp}}$ does not fit with the true global impedance $Z_{0-14}^{eq}$ (recall Fig. 7).
<figure><img src="image_13.png" /><figcaption>Fig. 12. Gain diagram of $Z_{0-14}^{eq}$ and $Z_{imp}$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Gain diagram of $Z_{0-14}^{eq}$ and $Z_{exp1}$.</figcaption></figure>

### 4.5. Implicit form approximation

In the previous Section 4.4, the multiple pole-zero approximation provides a very precise fit, however the number of parameters is significantly high (16 parameters), which shows that the model to be fitted in the middle frequency range could be a fractional order one. Such a number of parameter is too high and may lead to local minima in the optimization process, and especially in a real-time system identification.

Oustaloup's synthesis for fractional models (Oustaloup et al., 2000) helps approximating a fractional order operator in a frequency range $[\omega_b, \omega_h]$ as follows:

$$ 
\begin{aligned}
s_{[\omega_b, \omega_h]}^{\nu} &= C_0 \left( \frac{1 + \frac{s}{\omega_b}}{1 + \frac{s}{\omega_h}} \right)^{\nu} \\
&\approx \lim_{N \to \infty} C_0 \prod_{k=1}^{N} \frac{1 + \frac{s}{\omega_k}}{1 + \frac{s}{\omega_k}},
\end{aligned} 
\quad (31) $$

where

* $\nu \in \mathbb{R}$ defines the differential order

* $\omega'_1 = \eta^{1/2}\omega_b, \omega_1 = \alpha\omega'_1,$

* $\omega'_{k+1} = \alpha\eta\omega'_k,$

* $\omega_{k+1} = \alpha\eta\omega_k,$

* $\nu = \frac{\log \alpha}{\log(\alpha\eta)}, \alpha = \left(\frac{\omega_h}{\omega_b}\right)^{\frac{\nu}{N}}, \eta = \left(\frac{\omega_h}{\omega_b}\right)^{\frac{1-\nu}{N}}$

* $C_0 = \left| \frac{1 + j \frac{1}{\omega_h}}{1 + j \frac{1}{\omega_b}} \right|^\nu = \left( \frac{\omega_h}{\omega_b} \right)^\nu \left( \frac{1 + \omega_b^2}{1 + \omega_h^2} \right)^{\frac{1}{2}}$

From the previous Section 4.4 and given Oustaloup's synthesis, a relationship exists between the multiple pole-zero approximation and (31).

Since Oustaloup's synthesis converts an implicit form into a multiple pole and zero approximation, we can have a new model by combining $Z_{asymp}$ with (31) as follows:

$$ Z_{imp} = K \frac{(1 + \tau_2 s)^{\nu_2}}{1 + \tau_1 s^{\nu_1}}. \quad (32) $$

The optimization provides the estimation of 5 parameters which are gathered in Table 6.

Fig. 12 shows both the gain Bode diagram of $Z_{0-14}^{eq}$ and of its implicit approximation $Z_{imp}$ (in dashed star orange). As expected, the middle

frequencies are well estimated leading to a MSE of $1.096 \times 10^{-3}$ that is also very small compared to the asymptotic approximation.

Note that $\nu_1 - \nu_2 = 0.5$, a value that completely fits with the high frequency order constraint.

The resulting implicit model has a very good fit and a very small number of parameters, which makes it a serious candidate. Nevertheless, such an implicit form in the numerator is well suited for frequency approximation but in an objective of time domain simulation, and more specifically real-time domain system identification, an implicit form transfer function is yet not handled in the CRONE Toolbox. Therefore, such an approximation remains a very good fit, but further implementation developments are required.

### 4.6. Explicit forms and commensurate order approximations

Let us now consider a first explicit model:

$$ Z_{exp1} = K \frac{cs^{\gamma} + ds^{\delta} + 1}{as^{\alpha} + bs^{\beta} + 1}. \quad (33) $$

The optimization provides the estimation of 9 parameters which are gathered in Table 7.

Fig. 13 shows both the gain Bode diagram of $Z_{0-14}^{eq}$ and of the first fractional explicit form $Z_{exp1}$ (in star dotted purple). As expected, the middle frequencies are well estimated leading to a MSE of $3.099 \times 10^{-3}$ that is more than twenty times smaller than the asymptotic approximation.

As $d$ and $\delta$ are very close to 0, a simpler model can be found:

$$ Z_{exp2} = K \frac{cs^{\gamma} + 1}{as^{\alpha} + bs^{\beta} + 1}, \quad (34) $$

whose 7 optimized parameter are gathered in Table 8.

Fig. 14 shows both the gain Bode diagram of $Z_{0-14}^{eq}$ and of the second fractional explicit form $Z_{exp2}$ (in square dotted purple). The optimization result is even better than the first explicit model with two less parameters as the MSE is of $2.818 \times 10^{-3}$.

#### 4.6.1. Commensurate order model

In this latter explicit model, we notice that $\alpha$ is very close to $2\beta$ and $\gamma$ is very close to $\beta$, let us try to use a commensurate order model. In such a model, all the differentiation orders are integer multiple of a same order, called as the commensurate order. By doing so, the number of parameters gets even more reduced as only 5 parameters are estimated.

The fractional commensurate model has a compact form as follows:

$$ Z_{Comm} = K \frac{cs^{\alpha} + 1}{as^{2\alpha} + bs^{\alpha} + 1}. \quad (35) $$
<table><caption>**Table 7**
Estimated parameters of $Z_{exp1}$ approximation and the resulting MSE.</caption><thead><tr><th>K</th><th>a</th><th>b</th><th>c</th><th>d</th><th>α</th><th>β</th><th>γ</th><th>δ</th><th>MSE</th></tr></thead><tbody><tr><td>69 155.6216</td><td>7061.1382</td><td>82.1478</td><td>124.9752</td><td>0.0076</td><td>1.1308</td><td>0.5975</td><td>0.6304</td><td>-0.0409</td><td>3.099 × 10<sup>-3</sup></td></tr></tbody></table>

<table><caption>**Table 8**
Estimated parameters of $Z_{exp2}$ approximation and the resulting MSE.</caption><thead><tr><th>K</th><th>a</th><th>b</th><th>c</th><th>α</th><th>β</th><th>γ</th><th>MSE</th></tr></thead><tbody><tr><td>70 177.5267</td><td>7793.9494</td><td>81.2192</td><td>135.8739</td><td>1.1394</td><td>0.5936</td><td>0.6390</td><td>2.818 × 10<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Gain diagram of $Z_{0-14}^{eq}$ and $Z_{exp2}$.</figcaption></figure>

<table><caption>**Table 9**
Estimated parameters of $Z_{Comm}$ approximation and the resulting MSE.</caption><thead><tr><th>K</th><th>a</th><th>b</th><th>α</th><th>c</th><th>MSE</th></tr></thead><tbody><tr><td>69 379.9203</td><td>2897.3520</td><td>38.2663</td><td>0.4982</td><td>50.2583</td><td>1.986 × 10<sup>-2</sup></td></tr></tbody></table>

The optimization provides the estimation of 5 parameters which are gathered in Table 9.

Fig. 15 shows both the gain Bode diagram of $Z_{0-14}^{eq}$ and of its fractional commensurate model $Z_{Comm}$ (in dashed-dotted green). As expected, the middle frequencies are better estimated leading to a MSE of $1.986 \times 10^{-2}$ that is more than two hundred times smaller than the asymptotic model.

## 4.7. Summary and discussion

As the analytical global lung modeling of heat transfers is very complex. Indeed the equivalent global lung model $Z_{0-14}^{eq}$ (20) is expressed according to the impedances $Z_{1,n}$ and $Z_{3,n}$ of the 14 branch levels given in (18). In an objective of system identification in real-time, during cardiac surgery, the simulation of a more compact and a simple model is required.

The objective of this work is to propose an impedance approximation that behaves such as the true impedance $Z_{0-14}^{eq}$, in the perspective that the approximation model should be precise enough and implementable in real-time, more specifically in time-domain. This step of real-time system identification and time-domain simulation is proposed with the LMRPEM-2 algorithm in Victor et al. (2025) and is not detailed here.

To summarize the model order reduction, our 8 approximations are shown on Table 10. Since a compromise needs to be addressed between the model complexity and its accuracy, Table 10 presents the parameter number and the resulting MSE for each approximation. Moreover,

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Gain diagram of $Z_{0-14}^{eq}$ and $Z_{Comm}$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Absolute gain errors between the approximations and $Z_{0-14}^{eq}$.</figcaption></figure>

in order to highlight the accuracy differences, Fig. 16 compares the absolute gain errors to the analytical impedance $Z_{0-14}^{eq}$.

First of all, the asymptotic model is good enough in very low and very high frequency, but unfortunately its dynamical behavior is not suited in middle frequency. Indeed, it has the minimum requirements, namely a static gain and a fractional integral behavior in high frequency, which gives the lowest complexity. Even though it has the minimum number of parameters, it has also the highest MSE.

The most accurate model is the multiple pole and zero approximation that is based on the asymptotic model. It has the best MSE, but on the other hand, it has the maximal number of parameters, which is
<table><caption>**Table 10**
Summary of the approximations, the number of parameters and the resulting MSE.</caption><thead><tr><th>Approximation</th><th>n<sub>p</sub></th><th>MSE</th></tr></thead><tbody><tr><td>Z<sub>asyp</sub></td><td>2</td><td>3.636</td></tr><tr><td>Z<sub>F-PZ1</sub></td><td>5</td><td>1.615 × 10<sup>-1</sup></td></tr><tr><td>Z<sub>BW</sub></td><td>5</td><td>9.641 × 10<sup>-1</sup></td></tr><tr><td>Z<sub>multi-PZ</sub></td><td>16</td><td>9.33 × 10<sup>-4</sup></td></tr><tr><td>Z<sub>imp</sub></td><td>5</td><td>1.096 × 10<sup>-3</sup></td></tr><tr><td>Z<sub>exp2</sub></td><td>7</td><td>2.818 × 10<sup>-3</sup></td></tr><tr><td>Z<sub>Comm</sub></td><td>5</td><td>1.986 × 10<sup>-2</sup></td></tr></tbody></table>

16. In a context of real-time system identification, so many parameters could lead to local minima with a risk of no convergence.

A good approximation is the explicit model, and to reduce the number of parameters, the second one proposes to optimize only 7 parameters. In the scope of estimating all parameters in real-time, the most difficult part is to estimate efficiently the differentiation orders, however in the second explicit model, three differentiation orders need to be estimated which makes it more complex and time consuming in terms of convergence rate. The impact of estimating the differentiation orders is proposed in Duhé et al. (2024) and Victor et al. (2022) and shows that the coefficients will converge to the true parameters only once the differentiation orders have converged. So estimating three distinct differentiation orders in real-time might be very time consuming.

In-between, a set of 5 parameter approximation model have been found. The worst one is the Butterworth approximation, where some oscillations occur in low frequency, thus showing a resonance that is not present in the analytical model. Then, we enter in a class of models which proves that the analytical true impedance lung model also has a fractional behavior in middle frequency. The best model is the implicit approximation but such a model will be very difficult to simulate in time-domain and especially in real-time. The computation time may be over the sampling time (Victor et al., 2025). The single fractional pole-zero approximation is better than the Butterworth one, but it has the same drawback of having a resonance that is not present in the analytical model. The best compromise in terms of number of parameters and accuracy is the commensurate order approximation. Such a model is easily implementable and can be used in real-time by using the CRONE Toolbox in Matlab®.

## 5. Conclusion and final remarks

In this paper, the analogies between thermal conduction and electricity proved useful as a modeling tool by using the two-port network formalism. Such a formalism allows representing thermal diffusion as a circuit without restricting the frequency validity range. However, the exact analytical expressions of the impedances, both in series and shunt ones, resulted in a highly complex expression of the global lung impedance.

The two-port network formalism was used to obtain a global impedance model of bio-heat transfers in lungs. The frequency behavior of the thermal bio-heat transfers have enabled detecting the fractional behavior in middle and high frequency. Furthermore, the gain and phase Bode diagrams of the global lung models $Z_{0-n}^{eq}$ starting from branch level $n = 14$ up to the deeper levels (until $n = 23$) are very similar, therefore $Z_{0-14}^{eq}$ was used. Nevertheless, the complexity of this global model makes it difficult to be usable in real-time implementation and system identification.

On the other hand, one of the main aim of this study was to propose a compact (with few parameters) and simple global lung model which can be usable for real-time system identification. With this goal in mind, different approximations were proposed and optimized.

Even though an implicit form model was able to fulfill perfectly our accuracy and simplicity criteria, the current version of CRONE toolbox

is unable to support it for the present version. Nevertheless, such model could be used for future works.

The model recommended in this study is an explicit one with commensurate model structure. It gives pleasant results concerning accuracy and number of parameters. Such properties and simplicity makes it preferable in time-domain simulation and especially in real-time system identification.

For future works, the link with time-domain system identification should be explored. Implicit model time domain simulation should be also further studied. Finally, a study with real data set on sheep could be carried out to check the validity of the proposed analytical study and more specifically on the proposed approximations in frequency domain. Also, the optimization criterion could be improved by adding the phase error or by penalizing it with the number of parameters.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

S. Victor is an Editorial Board Member/Editor-in-Chief/Associate Editor/Guest Editor for this journal and was not involved in the editorial review or the decision to publish this article.

## Appendix A. Impedance expressions from the change of a matrix formulation to a two-port network circuit

For a T network as shown in Fig. 2, the input-output relationship may be expressed by means of a transmission matrix:

$$ \begin{bmatrix} T_{in}(s) \\ \phi_{in}(s) \end{bmatrix} = \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix} T_{out}(s) \\ \phi_{out}(s) \end{bmatrix}. \quad (A.1) $$

The transmission matrix can be expressed in function of the impedances $Z_1, Z_2$ and $Z_3$:

$$ \begin{bmatrix} A & B \\ C & D \end{bmatrix} = \begin{bmatrix} 1 + \frac{Z_1}{Z_3} & Z_1 + Z_2 + \frac{Z_1 Z_2}{Z_3} \\ \frac{1}{Z_3} & 1 + \frac{Z_2}{Z_3} \end{bmatrix}, \quad (A.2) $$

therefore, impedance expressions are obtained as the following:

$$ \begin{cases} Z_1 = \frac{A-1}{C}, \\ Z_2 = \frac{D-1}{C}, \\ Z_3 = \frac{1}{C}. \end{cases} \quad (A.3) $$

which lead to expressions in Eq. (15). It should be noted that $Z_1 = Z_2$ for this specific geometry as $A = D$, but it should not be taken as a general assumption.

## Appendix B. Effect of parameter variations on the global lung impedance model

In open-heart surgeries with ECC (or CPB), variations of the parameters provided in Table 1 can occur. Therefore, the influence of these parameters is now studied with a variation of ±10% of each parameter. The variations on $\rho_{tissue}$, $c_{tissue}$, $k_{tissue}$, $\rho_{blood}$, $c_{blood}$, $k_{blood}$ have no influence on the equivalent $Z_{0-14}^{eq}$ impedances. The only parameters affecting the most the frequency responses are those present in the early branch levels, when the effect of the bio-heat is not present. These parameters are the air diffusivity $a_{air}$ and the air thermal conductivity $k_{air}$ and their effect is shown in Figs. B.17 and B.18 respectively. The trachea impedance has a significant effect on the total impedance, consequently, these parameters also affect the global lung model.
<figure><img src="image_18.png" /><figcaption>Fig. B.17. Bode diagrams of $Z_{0-14}^{eq}$ with variations on $a_{air}$ of ±10% variation.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. B.18. Bode diagram of $Z_{0-14}^{eq}$ with variations on $k_{air}$ of ±10% variation.</figcaption></figure>

## References

Akaike, H. (1974). A new look at the statistical model identification. *IEEE Transactions on Automatic Control*, **19**, 716–723. http://dx.doi.org/10.1109/TAC.1974.1100705.

Ali, A., Radwan, A., & Soliman, A. (2013). Fractional order butterworth filter: Active and passive realizations. *IEEE Journal on Emerging and Selected Topics in Circuits and Systems*, **3**, 346–354. http://dx.doi.org/10.1109/JETCAS.2013.2266753.

Coccarelli, A., Boileau, E., Parthimos, D., & Nithiarasu, P. (2016). An advanced computational bioheat transfer model for a human body with an embedded systemic circulation. *Biomechanics and Modeling in Mechanobiology*, **15**, 1173–1190. http://dx.doi.org/10.1007/s10237-015-0751-4.

Danza, L., Belussi, L., Meroni, I., Salamone, F., Floreani, F., Piccinini, A., & Dabusti, A. (2016). A simplified thermal model to control the energy fluxes and to improve the performance of buildings. *Energy Procedia*, **101**, 97–104. http://dx.doi.org/10.1016/j.egypro.2016.11.013.

Diong, B., Nazeran, H., Nava, P., & Goldman, M. (2007). Modeling human respiratory impedance. *IEEE Engineering in Medicine and Biology Magazine*, **26**, 48–55. http://dx.doi.org/10.1109/memb.2007.289121.

DuBois, A., Brody, A., Lewis, D., & Burgess, B. (1956). Oscillation mechanics of lungs and chest in man. *Journal of Applied Physiology*, **8**, 587–594. http://dx.doi.org/10.1152/jappl.1956.8.6.587.

Duhé, J., Victor, S., Melchior, P., Abdelmoumen, Y., & Roubertie, F. (2021). Two-port network modeling for bio-heat transfers in lungs. In *11th IFAC symposium on biological and medical systems* (pp. 169–174). Ghent, Belgium: http://dx.doi.org/10.1016/j.ifacol.2021.10.250.

Duhé, J. F., Victor, S., Melchior, P., Abdelmoumen, Y., & Roubertie, F. (2022). Modeling thermal systems with fractional models: human bronchus application. *Nonlinear Dynamics*, **108**, 579–595. http://dx.doi.org/10.1007/s11071-022-07239-3.

Duhé, J. F., Victor, S., Melchior, P., Abdelmoumen, Y., & Roubertie, F. (2023). Online system identification of heat transfers in lungs with the LMRPEM-2 method. In *2023 international conference on fractional differentiation and its applications* (pp. 1–6). Ajman, United Arab Emirates: http://dx.doi.org/10.1109/ICFDA58234.2023.10153153.

Duhé, J. F., Victor, S., Melchior, P., Abdelmoumen, Y., & Roubertie, F. (2024). Recursive system identification for continuous-time fractional models. *IEEE Transactions on Control Systems Technology*, **1–13**. http://dx.doi.org/10.1109/TCST.2024.3407580.

Giesinger, R., & McNamara, P. (2018). The impact of therapeutic hypothermia on pulmonary hemodynamics of meconium aspiration syndrome. *American Journal of Respiratory and Critical Care Medicine*, **198**, http://dx.doi.org/10.1164/rccm.201802-0207LE.

Hildebrandt, J. (1969). Dynamic properties of air-filled excised cat lung determined by liquid plethysmograph. *Journal of Applied Physiology*, **27**, 246–250. http://dx.doi.org/10.1152/jappl.1969.27.2.246.

Hildebrandt, J. (1970). Pressure-volume data of cat lung interpreted by a plastoelastic, linear viscoelastic model. *Journal of Applied Physiology*, **28**, 365–372. http://dx.doi.org/10.1152/jappl.1970.28.3.365.

Ionescu, C., De Keyser, R., Sabatier, J., Oustaloup, A., & Lebron, F. (2011). Low frequency constant-phase behavior in the respiratory impedance. *Biomedical Signal Processing and Control*, **6**, 197–208. http://dx.doi.org/10.1016/j.bspc.2010.10.005, special Issue: The Advance of Signal Processing for Bioelectronics.

Ionescu, C. M., Hernandez, A., & De Keyser, R. (2012). Relationship between anatomical structure, fractional order models and fractal dimension in the respiratory impedance of healthy patients. *Vol. 45*, In *16th IFAC symposium on system identification* (pp. 1425–1430). http://dx.doi.org/10.3182/20120711-3-BE-2027.00170.

Ismail, N., Ghaddar, N., & Ghali, K. (2018). Electric circuit analogy of heat losses of clothed walking human body in windy environment. *International Journal of Thermal Sciences*, **127**, 105–116. http://dx.doi.org/10.1016/j.ijthermalsci.2018.01.025.

Jiang, G., Qu, T., Shang, Z., & Zhang, X. (2004). A circuit simulating method for heat transfer mechanism in human body. In *The 26th annual international conference of the IEEE engineering in medicine and biology society* (pp. 5274–5276). http://dx.doi.org/10.1109/IEMBS.2004.1404472.

Kennedy, J., & Eberhart, R. (1995). Particle swarm optimization. In *Proceedings of iCNN'95 - international conference on neural networks* (pp. 1942–1948). http://dx.doi.org/10.1109/ICNN.1995.488968.

Künzi, R. (2015). Thermal design of power electronic circuits. *ArXiv: Accelerator Physics*, **3**, 311–327. http://dx.doi.org/10.5170/CERN-2015-003.311.

Kuwahara, F., Sano, Y., Liu, J., & Nakayama, A. (2009). A porous media approach for bifurcating flow and mass transfer in a human lung. *Journal of Heat Transfer*, **131**, http://dx.doi.org/10.1115/1.3180699.

Ljung, L. (1999). *System identification – theory for the user* (2). Upper Saddle River, N.J. USA: Prentice-Hall.

Magin, R. (2004). Fractional calculus in bioengineering. *Critical Reviews in Biomedical Engineering*, **32**, 1–104. http://dx.doi.org/10.1615/critrevbiomedeng.v32.i1.10.

Maillet, D., André, S., Batsale, J., Degiovanni, A., & Moyne, C. (2000). *Loyola symposium series, Thermal quadrupoles: solving the heat equation through integral transforms*. John Wiley & Sons.

Navajas, D., Farre, R., Canet, J., Rotger, M., & Sanchis, J. (1990). Respiratory input impedance in anesthetized paralyzed patients. *Journal of Applied Physiology*, **69**, 1372–1379. http://dx.doi.org/10.1152/jappl.1990.69.4.1372.

Oustaloup, A., Lebron, F., Mathieu, B., & Nanot, F. (2000). Frequency-band complex noninteger differentiator: characterization and synthesis. *IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications*, **47**, 25–39. http://dx.doi.org/10.1109/81.817385.

Parnis, G. (2012). Building thermal modelling using electric circuit simulation. URL http://handle.unsw.edu.au/1959.4/51971.

Pellet, M., Melchior, P., Abdelmoumen, Y., & Oustaloup, A. (2011). Fractional thermal model of the lungs using havrioli-negami function. In *ASME 2011 international design engineering technical conferences and computers and information in engineering conference* (pp. 805–810). Washington DC USA: http://dx.doi.org/10.1115/DETC2011-48095.

Pennes, H. (1948). Analysis of tissue and arterial blood temperatures in the resting human forearm. *Journal of Applied Physiology*, **1**, 93–122. http://dx.doi.org/10.1152/jappl.1948.1.2.93.

Rogers, C., Mazza, G., Maishman, R., Thirard, R., Evans, J., de Jesus, S., Beard, C., Angelini, G., Millar, A., Jarad, N., Tomkins, S., Hillier, J., Suleiman, M., & Ascione, R. (2024). Low frequency ventilation during cardiopulmonary bypass to protect postoperative lung function in cardiac valvular surgery: The protection phase II randomized trial. *Journal of the American Heart Association*, **13**, Article e035011. http://dx.doi.org/10.1161/JAHA.124.035011.

Victor, S., Duhé, J., Melchior, P., Abdelmoumen, Y., & Roubertie, F. (2022). Long memory recursive prediction error method for identification of continuous-time fractional models. *Nonlinear Dynamics*, **110**, 635–648. http://dx.doi.org/10.1007/s11071-022-07628-8.
Victor, S. & Malti, R. (2023). Crone toolbox for system identification and modeling using fractional differentiation models. In *22nd IFAC world congress* (pp. 120–125). Yokohama, Japan: http://dx.doi.org/10.1016/j.ifacol.2023.10.1556.

Victor, S., Melchior, P., Pellet, M., & Oustaloup, A. (2020). Lung thermal transfer system identification with fractional models. *IEEE Transactions on Control Systems Technology*, *28*, 172–182. http://dx.doi.org/10.1109/TCST.2018.2877606.

Victor, S., Ndreko, E., Duhé, J., & Melchior, P. (2025). Real-time system identification of bio-heat transfers in lungs. *Nonlinear Dynamics*, *113*, 18645–18667. http://dx.doi.org/10.1007/s11071-025-11105-3.

Weibel, E. R. (2005). Mandelbrot's fractals and the geometry of life: A tribute to benoît mandelbrot on his 80th birthday. In G. A. Losa, D. Merlini, T. F. Nonnenmacher, & E. R. Weibel (Eds.), *Fractals in biology and medicine* (pp. 3–16). Basel: Birkhäuser Basel, http://dx.doi.org/10.1007/3-7643-7412-8_1.