Conceptual model of the respiratory system
Our modelling and simulation for reproducing airflow signals required simplification of airway structure in the lung. The geometrical and dimensional structures of the airways were proposed by several researchers, among which Horsefield et al. [10] and Weibel [31] dissected and measured the human lungs and airways, and more recently, and Tawhai et al. [11] and Kitaoka et al. [14] proposed algorithmic approaches for reconstructing the branching structure. We employed Weibel's morphometry of the lung, which provides the geometries and the dimensions based on the symmetric dichotomous structure of airway branching when the lung volume is assumed to be 75% of the total lung capacity (TLC). According to his morphometry, the trachea is defined as airway generation 0 and it is separated into two geometrically identical daughter branches. Each daughter is repeatedly branched up to 22 times, thus the lung is considered to have 24 (0–23) airway generations and the number of airway branches total ~1.7 × 107. In each airway generation, there are 2z(z is a generation number) identical branches whose lengths and radii are provided and the dimensions of airway branches in each generation differ from generation to generation. He also suggested that airway generations 0–16 comprise the conducting zone, whereas airway generations 17 – 23 make up the respiratory zone in which gases are exchanged.
Based on Weibel's morphometry of the lung, we simplified the geometry of the airways. For airway generation Z in conducting zone, a bundle of 2zidentical airway branches are considered as the big tube whose cross sectional area equals to 2ztimes the cross sectional area of a single airway branch. The tubes for each airway generation are represented by RCL T-networks shown in fig. 1 (b), and evaluation of the R's, C's, and L's in the RCL T-network are explained in equations (4) – (6). Meanwhile, the respiratory zone is considered as a big lump, the alveolar space, since alveolar ducts and sacs are scattered throughout the respiratory zone [8]. Although the upper airway is not presented in Weibel's morphometry, it is one of the chief sites for airway resistance.
Our conceptual model as a summary of the simplification, 'the upper airway + 17 conducting airway generations + the alveolar space', is demonstrated in fig. 1 (a), and circuit analogue converted from fig. 1 (a) is shown in fig. 1 (c).
Dimensions of airway branches
Since the dimensions of every single airway branch in the lung vary during respiration, it is crucial to track them to extract proper values of R's, C's, and L's. The radii of the airway branches in airway generation Z are determined by transmural pressure (P
tmZ
), pressure difference between the pleural cavity and inside airway generation Z. In order to determine P
tmZ
, we need to rely on the study about transpulmonary pressure (P
tp
), pressure difference between the pleural cavity and alveolar space. According to Salizar et al. [26], P
tp
is a function of lung volume (LV) and TLC. With given LV, P
tp
is obtained by,
P
tp
= -log (1 - LV / TLC) × 7.22. (1)
If LV is functional residual capacity (FRC) and no airflow exists in the airways, the lung is at rest. In this situation, air pressure at any site is the same and atmospheric (zero), therefore P
tmZ
and P
tp
are the same. Our simulation begins with assuming that the lung was at rest. The equation of Lambert et al., or Lambert's tube law [16], is a function of P
tmZ
, and gives the ratio of the cross sectional area (A
Z
) to the maximum cross sectional area of an airway branch (A
maxZ
) in airway generation Z. That is,
A
Z
/ Amax Z= 1.0 - (1.0 - α0)(1.0 - P
tmZ
/ P0)-N (2)
where α0, α0' and N are constants for each airway generation given from Lambert's tube law, and P0 = (α0-1)N/α0'. As mentioned earlier, the dimensions of the airway branches provided by Weibel's morphometry are based on 75 % of TLC. Let A
Z75
be Weibel's cross sectional areas in airway generation Z. By using
A
maxZ
for each airway generation can be found. Since A
maxZ
does not change whether or not the lung is at rest, the equation is valid to find out A
Z
during respiration.
The length of an airway branch is assumed to vary in away that;
l
z
/ l
FRC
= r
z
/ r
FRC
(3)
where l
z
and r
z
are the length and the radius of an airway branch in airway generation Z, respectively. l
FRC
and r
FRC
are the length and the radius on FRC, respectively.
Equivalent circuit elements modelling
The RCL T-network consists of two resistors (R
gZ
), two inductors (L
gZ
), a capacitor (C
gZ
), and a DC voltage source (P
tm_ini
). R
gZ
and L
gZ
are the airway resistance and the inertance of air in the entire airway generation Z, respectively. C
gZ
, the airway capacitance of the entire airway generation Z, represents inflation and deflation of the non-rigid airway wall during respiration. P
tm_ini
is the initial value of P
tmZ
. This initial pressure counterbalances P
pl
, which is negative when the lung is at rest.
The values of the circuit elements in the RCL T-network of airway generation Z can be obtained by totalling the 2zairway branches. Therefore, R
gZ
is given by
Similarly, L
g
z is given by
and C
gZ
is
C
gZ
= C
sZ
× N
z
[ml/cmH2O], (6)
where R
sZ
, L
sZ
, and C
sZ
indicate the values of a single airway branch in airway generation Z. N
z
is 2z, the number of airway branches in generation Z.
R
sZ
depends on the classic Poiseuille equation and the Zeta correction factor proposed by Pedley et al. [21]. The Zeta correction factor is given by
and
R
e
is Reynolds number 
, V is air velocity, ρ is the density of air, and η is the viscosity of air.
Because electric current flows at the speed of light in an electrical circuit, it is necessary to compensate for the incomparably slower behaviour of airflow in the airways by placing inductors. Inertance of air in a single airway branch (L
sZ
) depends on the dimension of the airway branch [5] and it is given by
To compute the values for C
sZ
in equation (2), relation between A
Z
and P
tmZ
in airway generation Z, is used again. That is,
The alveolar space is represented by a capacitor in the circuit analogue. Alveolar capacitance (C
as
) is the ratio of the alveolar volume (ΔAV) to ΔP
tp
. Lung volume (LV) consists of airway volume and alveolar volume, and airway volume is negligible compared to alveolar volume. Therefore,
which implies that
represents C
as
. It can be obtained by equation (1); however, the equation did not consider the hysteresis of P
tp
-volume curves that is caused by several proposed reasons [1]. To overcome this, we defined C
as
as 
*k (k<0), and k was continuously changed over the time course of the simulation. The values of k were empirically determined to achieve the acceptable shape of the hysteresis, which is shown in fig. 2.
A resistor and an inductor characterize the upper airway, from the nasal/oral cavity to larynx. The equations for the resistance of the upper airway (R
ua
) that Jackson et al. [12] validated are described below:
where F
ua
is the airflow rate in the upper airway. As Marchal et al. [18] estimated, the value of the inductor that represents inertance of air is 0.00003 [cmH2O·s2/ml].
Numerical methods for nonlinear circuit analogue
Although the number of elements is manageable and the structure of circuit is fairly simple, analysis of the circuit analogue in fig. 1 (c) is not trivial since the values of all the energy-storing (C's and L's) elements as well as resistive elements (R's) change at every sequence of the simulation time-step. To exemplify a general idea of numerical methods for the nonlinear circuit analogue, analysis of a simple nonlinear second-order system in fig. 3 was demonstrated. The system of fig. 3 is described by an equation:
where v
C
is the voltage difference between each end of the capacitor. Note that R(•) is a function of i and v
C
, and L(•) and C(•) are functions of v
C
just like in fig. 1 (c). Using backward Euler approximation [22], equation (14) is converted to a difference equation:
where Δt is time difference between sequence [n] and [n-1]. Equation (15) is the same as
By the definition of backward Euler approximation, the current of the circuit
. Then equation (16) can be restated as an equation;
In equation (17), the right hand side consists of all known values, and the left hand side is a function of v
C
[n]. Suppose that v
C
[n] is x, equation (17) can be expressed as
f(x) = c, (18)
where c is a constant. Equation (18) can be easily solved using an iteration method [24]. To do this computation, MATLAB (Mathworks, Natick, MA) codes were written.
Protocols for measurement airflow signals
In this study, The Vest™ (Advanced Respiratory, St.Paul, MN, USA; now named Hill-Rom Co.,Inc.), which delivers sine waveform compression pulses, was used for application of HFCC to a subject for measuring and recording the airflow signals at the mouth. The general usage of HFCC device was early described in the review of Hansen et al. [9] and the typical respiratory airflow during use of HFCC is shown in fig. 4.
The subject sat upright on a chair for measuring and recording the airflow signals at the mouth with an in-house built electronic spirometer. The subject worn a nose clip and breathed through a mouthpiece. After the HFCC device is properly set up, the compression pulses were applied and then the subject made several slow and large, but not to TLC, breaths. During the breaths the subject hold his glottis open until data collection was completed. This protocol was followed for the low (5 Hz), high (21 Hz), and medium frequencies (15 Hz) of HFCC pulses. Before each frequency recording of airflow signals, the subject rested for one minute. To ensure that the subject had adapted to HFCC pulses and had reached a steady state, only the last ten seconds of the one- minute breathing were recorded and analyzed for our study.
, V is air velocity, ρ is the density of air, and η is the viscosity of air.
