Toward the modeling of mucus draining from human lung: role of airways deformation on air-mucus interaction

Chest physiotherapy is an empirical technique used to help secretions to get out of the lung whenever stagnation occurs. Although commonly used, little is known about the inner mechanisms of chest physiotherapy and controversies about its use are coming out regularly. Thus, a scientific validation of chest physiotherapy is needed to evaluate its effects on secretions. We setup a quasi-static numerical model of chest physiotherapy based on thorax and lung physiology and on their respective biophysics. We modeled the lung with an idealized deformable symmetric bifurcating tree. Bronchi and their inner fluids mechanics are assumed axisymmetric. Static data from the literature is used to build a model for the lung's mechanics. Secretions motion is the consequence of the shear constraints apply by the air flow. The input of the model is the pressure on the chest wall at each time, and the output is the bronchi geometry and air and secretions properties. In the limit of our model, we mimicked manual and mechanical chest physiotherapy techniques. We show that for secretions to move, air flow has to be high enough to overcome secretion resistance to motion. Moreover, the higher the pressure or the quicker it is applied, the higher is the air flow and thus the mobilization of secretions. However, pressures too high are efficient up to a point where airways compressions prevents air flow to increase any further. Generally, the first effects of manipulations is a decrease of the airway tree hydrodynamic resistance, thus improving ventilation even if secretions do not get out of the lungs. Also, some secretions might be pushed deeper into the lungs; this effect is stronger for high pressures and for mechanical chest physiotherapy. Finally, we propose and tested two a dimensional numbers that depend on lung properties and that allow to measure the efficiency and comfort of a manipulation.

1. Introduction

Lung secretions are used to protect the lungs against pollutants or infections, by capturing airborne particles or cells that may enter the lungs. Secretions form an heterogenous gel-like fluid consisting mostly in water and biopolymers (Lai et al., 2009). Secretions are motioned upwards thanks to two natural mechanisms. First, the cilia located on the bronchi walls, whose motion pushes the secretions upward the bronchial tree, this phenomenon is called mucociliary clearance (Grotberg, 2001; King, 2006; van der Schans, 2007). The second mechanism is cough, it involves high air flows into the bronchi (Basser et al., 1989; Zahm et al., 1991; Vandevenne, 1999). These mechanisms can however fail. For example, cough or cilia motion are not mature in young children, they can also been altered by pathologies or age. As a consequence, when secretions are too thick or overproduced, they can become stalled and bronchial infections may occur. Chest physiotherapy is used to treat patients whose secretions natural clearing mechanisms are altered. The techniques depend on the pathology and on the patient characteristics, typically her/his age. Bronchiolitis, asthma, chronic obstructive pulmonary disease or cystic fibrosis are typical diseases treated with dedicated chest physiotherapy techniques.

Chest physiotherapy is based on external maneuvers on the chest. It relies on a basic concept: the pressures applied by the practitioner on the chest transmit into the lungs and deforms the bronchi. Bronchi's volume changes create an air flow inside the lung and the mechanical constraints applied on the secretions by the air flow affect the secretions and have the potential to induce motion. Secretions are gel-like materials which are motioned only if they are submitted to a stress sufficiently strong. In manual chest physiotherapy, flow and lung volume are modulated by the practitioners during the manipulations in response to the feedback they feel from the patient. Mechanical devices based on this concept have also been developed (Arens et al., 1994; Langenderfer, 1998; Hristara-Papadopoulou et al., 2008). This concept represents however a very simplified view of the lung and thorax biomechanics and chest physiotherapy remains as of today a discipline which is mostly empirical. Maneuvers are efficient but the inner mechanisms involved in their efficiency remain either unknown or non-proved in a scientific way. Many chest physiotherapy techniques exist, and their definition and whether they should be used or not are decided during country wide consensus meetings. Consequently, the available pool of maneuvers can differ from one country to the other, mostly because their efficiency are difficult to compare (Pryor, 1999; Main et al., 2005).

Historically, most of the maneuvers were based on applying relatively strong pressures on the chest, bringing discomfort and, if ill performed, to risks to the patient, most particularly to baby and children. Thus, the use of manual chest physiotherapy for baby suffering bronchiolitis—a well-known and wide spread disease—has been questioned because of the possible risks induced by the manipulations (Postiaux and Lens, 1992; Beauvois et al., 2007). Recently, new techniques focusing on patient comfort and using low pressures on the chest, low air flows in the lungs and induced cough have emerged and have spread quickly in physiotherapists' community. The efficiency of these new low flow techniques is however as of today still controversial, because chest physiotherapy is based on the hypothesis that secretions can be mobilized (i.e., put to motion) only with high air flows.

Chest physiotherapy has reached a state where it needs to be supported by a strong scientific background. The different questions raised by physiotherapists would benefit from a better understanding of the coupled biomechanics of chest and lungs. An integrated model of chest physiotherapy that is able to uncover the effects of hand pressures on secretions mobilization would help the practitioners to chose and define knowledgeably the techniques they are using. Most particularly, the role of the bronchial tree geometry play an important role on the air fluid dynamics (Mauroy et al., 2004; Mauroy and Meunier, 2008; Bokov et al., 2010, 2014; Mauroy and Bokov, 2010) and may have an important role on secretions mobilization. Hence, in a previous paper (Mauroy et al., 2011), we studied the role of the interaction between air flow, secretion thickness and bronchial tree geometry assuming it was rigid. This first theoretical frame contained some of the core elements involved in chest physiotherapy and we showed that air flow needed to be high enough for secretions to be motioned and that applying the same given air flow in the lung leads to a stationary distribution of mucus, i.e., a constant flow manipulation might see its efficiency decrease over time. These first results indicate that air flow needs to be modulated toward higher air flows in order to be able to mobilize mucus all along the manipulation. Another method for increasing the air flow induced stress might be to decrease the volume of the bronchi while keeping the air flow rate constant. This is the idea used, for example, in low volume manipulations (Pryor, 1999) or in the mechanical chest physiotherapy technique of chest compression (Arens et al., 1994; Langenderfer, 1998). In order for our model to be able to capture such behaviors and to give some qualitative and comparative predictions, it was necessary to add to our previous model how lung, and its subunits bronchi and acini, deform during the manipulations. Moreover, integrating lung deformations would improve globally the quality of the predictions of our model.

Thus, in this work, we propose a new model that integrates lung deformations. The model input is an homogeneous pressure applied on the chest, and the model output is the air flow, the secretions motion and distribution and the lung configuration (airways sizes, tissue pressure, etc.). In Section 2, we describe the hypotheses of the model and in Section 3, we present and discuss the predictions of the model in the case of idealized manipulations of manual chest physiotherapy and of different devices of mechanical chest physiotherapy.

2. Methods

This study is based on three separate sub-models: one for the lung, one for the thorax and one to mimic chest physiotherapy manipulations. We build a minimal model, which includes only the core properties we identified as playing an important role in chest physiotherapy. This choice was mainly made in order to be able to interpret correctly the interactions between these core phenomena, it was also a way to limit the computation times of the numerical studies. As a consequence, our model can perform qualitative and comparative predictions, but not quantitative predictions.

This work assumes that lung tissue mechanics and lung generations are homogeneous. These choices were made to keep the model tractable in term of physical behavior. Thus, we chose to neglect phenomena such as hydrostatic pressures, inhomogeneous pressure distribution in the lung tissue or differences between airways or alveoli belonging to the same generations. The details of the model are given in the following subsections.

2.1. Terminology

Below is a list of definitions for some terms or variables used in this work.

FRC: functional residual volume, the volume of the lung at the end of a normal expiration.

Generation: the generation of an airway is an index accounting for the number of bifurcations in the path from the tree root (trachea) to that airway (tree root is zero, deepest alveoli are 23). We denote z(i) the generation index of airway i.

Lung volume: lung volume is denoted V_L it is the sum of the tracheobronchial tree volume V_tbt and of the acini volume V_ac.

Tree structure: for an airway indexed i, its parent branch is indexed p, and its daughter branches are indexed d₁ and d₂. Although p, d₁ and d₂ are function of i, we do not make i appear for the sake of notations simplification.

1 cmH₂O: 1 cmH₂O is equivalent to 98.07 Pa.

Φ^(z)_a: air flow in an airway in generation z.

Φ^(z)_{m, *}: mucus flow in an airway in generation z.

* = I for mucus inflow and * = O for mucus outflow.

Φ^(z)_alv: air flow in an alveolus connected to an airway in generation z.

l_z: length of the bronchi from generation index z.

max(a, b): maximum function, max(a, b) = a if a > b and max(a, b) = b otherwise.

n^(z)_alv: number of alveoli connected to an airway in generation z (0 for airways which are not alveolar ducts, 58 for alveolar ducts, i.e., whose generation is deeper than 17).

P_ext: external chest wall pressure. It is assumed constant everywhere on the chest.

P^(z)_air: air pressure at mid-airway length for airway i.

P_tissue(V_L): mean lung tissue pressure when the volume of the lung is V_L. This data comes from static lung volume-pressure relationships measured by Agostoni et al. (Agostoni and Hyatt, 2011).

r^(z)_a: air lumen area radius.

r^(z)_b: airway lumen area.

S^(z)_a: air lumen area of bronchus i. Air lumen area is assumed disk shaped.

S^(z)_b: airway i lumen area. Airway lumen area is assumed disk shaped.

V_{L, rs}(ΔP): Lung volume when it is submitted to a pressure difference ΔP between the thorax and the bronchi (same pressure into the whole lung, no air movement). This quantity is built from static lung volume-pressure relationships measured by Agostoni et al. (Agostoni and Hyatt, 2011).

z: generation index in the tree.

Note that the subscript (z) corresponds to the airway generation associated to the data. In the case of properties that are applicable to any airway, we may drop this subscript for the sake of formulas simplification.

2.2. Lung's Geometry

The geometrical model of the lung used in this work is scaled to have a total lung capacity (TLC) of 6.5 L, a vital capacity (VC) of 5 L, a residual volume (RV) of 1.5 L and a tidal volume at rest of 0.5 L (Weibel, 1984; Agostoni and Hyatt, 2011). According to Agostoni and Hyatt (2011), functional residual volume (FRC) is reached at 35% of VC, thus FRC = RV + 0.35 × VC is 3.25 L.

2.2.1. Bronchial Tree Model

The bronchial tree is modeled as an idealized airway tree with 23 generations. The generation index starts from generation zero (trachea) to generation 22 (deepest alveolar duct) and is incremented by one at each bifurcation. Although bronchial bifurcations are slightly asymmetrical (Tawhai et al., 2004; Mauroy and Bokov, 2010; Florens et al., 2011; Clément et al., 2012), we approximate the bronchial tree with a symmetric airway tree for our pioneering model to remain tractable. We assume that each airway bifurcates into two smaller identical airways, implying that all the airways belonging to a same generation are identical.

The sizes and lengths of the conducting airways (seventeen first generations) are chosen from Lambert's data (Lambert et al., 1982), the different parameters are recalled in Table 1. The airways are modeled as cylindrical tubes. We assume that mucus distribution in an airway is a constant thickness layer along the cylinder wall and that mucus cannot reach the alveoli.

TABLE 1

Table 1. Lung's data for the geometry of the conductive airways at FRC = 1.75 L (end expiration volume), from Lambert et al. (1982): z is the generation index; l is the airway length; d is the airway diameter.

2.2.2. Alveoli and Acini Geometrical Models

Alveoli are assumed to be spherical and to stem from airways walls from the seventeenth generation. The total number of alveoli N_a is documented in the literature (Ochs et al., 2004) and we used N_a = 480 10⁶. FRC volume can be decomposed as the volume of the bronchial tree plus the volume of the alveoli. Our model computations give a FRC volume for the bronchial tree of about 0.67 L and a FRC total alveoli volume of about 2.58 L. Consequently, the volume of an alveolus at FRC used in our model is v_a = 5.37 10⁶ μm³ (diameter of 217 μm). This value is slightly higher than the value measured in Ochs et al. (2004) (4.2 10⁶ μm³, i.e., a diameter of 200 μm), the discrepancy may arise from the measures being done in a lung's volume that is not FRC. The exchange surface in the lungs consists in the total alveoli surface and can be estimated in our model for FRC to ~ 70 m².

Alveoli are stemming from bronchi walls in the acini units, that corresponds to the last sixth generations of the lungs. In our model, each acinus has 2⁶ − 1 = 63 alveolar ducts. On the contrary of bronchi in the 17 first generations whose diameters and lengths decrease when generation index increases, the bronchi in the acini, called alveolar ducts, keep roughly similar size whatever their generation (Weibel, 1984). Thus, we assume that the number of alveoli is homogeneously distributed on the alveolar ducts walls, and to reach a total of 480 10⁶ alveoli, the amount of alveoli per alveolar ducts in the acini is 58, a bit higher than the range from 20 to 40 reported in Whimster (1970). It remains however fully in the range of physiological patterns, since there is quite large variation between individuals, as reported for example in Ochs et al. (2004).

2.3. Lung's Mechanics

Our goal is to be able to compute the volume of the airways and alveoli as a function of two input pressures: the pressure applied on the thorax P_ext and the fluid pressures inside the airways P^(z)_air (z is the generation index). In our model, the pressure applied on the thorax may either come from muscles action, like in ventilation, or from chest physiotherapy manipulation; the pressure inside the airways is a consequence of the fluid mechanics in the tree. For the sake of simplification, we assume that the pressure in the lung tissue P_tissue to be homogeneous throughout the lung. We assume that lung tissue pressure P_tissue depends on the lung volume only and is computed from the static volume-pressure relationship from Agostoni et al. work (D'Angelo and Milic-Emili, 2005; Agostoni and Hyatt, 2011), see Figure 1 (dashed curve). Similarly, we assume that bronchi and alveoli mechanics do not depend on their location in the lung.

FIGURE 1

Figure 1. Data (static) from Agostoni et al. (Agostoni and Hyatt, 2011). Continuous curve: total respiratory lung volume-pressure relationship, used in our work to model acini mechanics. Dashed curve: Lung tissue pressure, used in our work to compute bronchi transmural pressure.

To model the lung mechanics, we divided the tree in two parts: the conductive airways and the acini. Mechanics of airways is based on the transmural pressure-airways section relationships proposed by Lambert et al. (1982). Mechanics of the acini is based on the total respiratory system lung volume-pressure relationship from Agostoni et al. (D'Angelo and Milic-Emili, 2005; Agostoni and Hyatt, 2011), see Figure 1 (continuous curve).

2.3.1. Conductive Airways Mechanics (Up to the 17th Generation)

To compute airways diameters, we used data from Lambert et al. (1982). Lambert et al. built data based static relationships of bronchi's airways lumen area vs. their transmural pressure for the seventeenth first generations. Their estimations are given per generation. For an airway belonging to generation z with a transmural pressure ΔP, they proposed that the airway lumen area S_z(Δp) is a sigmoid function of transmural pressure:

$\begin{array}{cc}{{𝒮}}_z(\Delta P)=\{\begin{array}{ll}{\alpha }_0(z){\left(1-\frac{\Delta P}{P_1(z)}\right)}^{-n_1(z)}{A}_m(z)\hfill & \text{when}\Delta P\le 0\hfill \\ \left(1-(1-{\alpha }_0(z)){(1-\frac{\Delta P}{P_2(z)})}^{n_2(z)}\right)\hfill & \hfill \\ A_m(z)\hfill & \text{when}\Delta P>0\hfill \end{array}& \left(1\right)\end{array}$

with $P_1(z)=\frac{{\alpha }_0(z)n_1(z)}{{{\alpha }^{\prime }}_0(z)}$ and $P_2(z)=-n_2(z)\frac{1-{\alpha }_0(z)}{{{\alpha }^{\prime }}_0(z)}$. A_m represents the maximal possible lumen area. The quantity α₀A_m is the airway area when transmural pressure is 0. At functional residual regime (FRC, end of normal expiration), the transmural pressure is about 500 Pa. The values used for the different parameters in the previous formulas are detailed in Table 2; the dependence of the bronchi surface area with transmural pressure are plotted on Figure 2.

TABLE 2

Table 2. Lung's data for the static mechanical behavior of the conductive airways, from Lambert et al. (1982).

FIGURE 2

Figure 2. Dependence of bronchi surface area (normalized) with bronchi transmural pressure. The numbers correspond to the generation index and the red dashed line represents bronchi state at FRC.

Airways transmural pressure is the difference between the pressure in the lung's tissues around the bronchi and the bronchi inner pressure due to fluid mechanics in the lung. The fluid mechanics in the airway tree is detailed below. We assume that tissue pressure depends on the volume of the lungs only; it is computed from the volume—lung pressure static relationship measured by Agostoni et al. (Agostoni and Hyatt, 2011). Using our terminology, the lumen area S^(z)_b of a bronchi in generation z (z ≤ 16) is

$\begin{array}{cc}{S}_b^{(z)}={{𝒮}}_z\left(P_{tissue}(V_L)-P_{air}^{(z)}\right)& \left(2\right)\end{array}$

We recall that P_tissue(V_L) is given by the volume-pressure relationship from Agostoni et al. (Agostoni and Hyatt, 2011), see Figure 1 (dashed curve). To simplify, we will assume that P^(z)_air is evaluated in the middle of the bronchus z. Moreover, the pressure—lumen area relationships assumes that airway walls are always at equilibrium with their corresponding transmural pressure, thus bronchi mechanics in our model is quasi-static. From these laws, we can derive the volume of the tracheobronchial tree,

$\begin{array}{cc}{V}_{tbt}(V_L,P_{air})={\displaystyle \sum _{z=0}^{16}{2}^z}{{𝒮}}_z\left(P_{tissue}(V_L)-P_{air}^{(z)}\right)L_b^{(z)}& \left(3\right)\end{array}$

2.3.2. Alveolar Ducts and Alveoli Mechanics (Last Six Generations)

The last sixth generations of the lung correspond to the acini. Each acinus in our model consists in 2⁶ − 1 = 63 alveolar ducts and each alveolar duct is connected to 58 alveoli (see model geometry section). An alveolar duct with its corresponding 58 alveoli will be referred to as an alveolar duct unit. There is 2¹⁷ × 63 alveolar duct units in our model. Lambert et al. data does not cover those generations, thus we used another approach to approximate the mechanical behavior of an alveolar duct unit.

The pressure difference between the inside of an alveolar duct unit in generation z and the body surface is P^(z)_air − P_ext. Consequently, the volume the lung would have under static condition with the same pressure difference in the whole tree is given by the respiratory system pressure—volume relationship from Agostoni and Hyatt (2011), see continuous curve on Figure 1. However, the lung volume measured by Agostoni et al. includes the tracheobronchial tree volume, thus to reach the volume of acini only, it is necessary to subtract from the lung volume the volume of the tracheobronchial tree as it is in the configuration of Agostoni et al. measurements: static and the glottis open. As a consequence, the reference pressure (= atmospheric pressure) spans the whole tree. Hence, the volume of an alveolar duct unit in generation index z in the tree is:

$\begin{array}{cc}\begin{array}{l}{v}_{adu}\left(P_{air}^{(z)},P_{ext}\right)=\frac{1}{2^{17}\times 63}(V_{L,rs}(P_{air}^{(z)}-P_{ext})\\ -V_{tbt}\left(V_{L,rs}(P_{air}^{(z)}-P_{ext}),0\right))\end{array}& \left(4\right)\end{array}$

Now, to compute alveolar ducts and alveoli volume, we assume their respective volume ratio in an alveolar duct unit remains unchanged whatever the lung volume. We assume the length of alveolar ducts to be L_ad = 0.7 mm (Haefeli-Bleuer and Weibel, 1988). The quantity α represents the volume ratio of an alveolar duct into an alveolar duct unit. Hence, for airways in a generation z in the acini (z > 16):

$\begin{array}{cc}\{\begin{array}{ll}{S}_b^{(z)}=\alpha v_{adu}\left(P_{air}^{(z)},P_{ext}\right)/L_{ad}\hfill & \text{(alveolar duct lumen area)}\hfill \\ v_{alv}^{(z)}=(1-\alpha )v_{adu}\left(P_{air}^{(z)},P_{ext}\right)\hfill & \text{(alveolus volume)}\hfill \end{array}& \left(5\right)\end{array}$

We compute α from FRC state and assume it is a constant: α = 0.17. Finally, we can compute the total volume of the acini in the tree,

$\begin{array}{cc}{V}_{ac}(P_{air},P_{ext})={\displaystyle \sum _{z=17}^{22}{2}^z}{v}_{adu}\left(P_{air}^{(z)},P_{ext}\right)& \left(6\right)\end{array}$

2.3.3. Computation of the Lung Volume V_L.

To compute the conductive airways lumen areas with Equation (3), we need to compute the volume of the lung V_L. Acini volume is known through Equation (4) and lung volume is the sum of acini volume and of conductive airways volume. Consequently, lung volume V_L depends on P_air and P_ext and is the solution of the equation:

$\begin{array}{cc}{V}_L=V_{ac}(P_{air},P_{ext})+V_{tbt}(V_L,P_{air})& \left(7\right)\end{array}$

Once computed, lung volume V_L is used to compute the tracheobronchial tree lumen areas using Equation (2). Consequently, the airway lumen area of generations z ≤ 16 expressed in Equation (2) can be reformulated as a function of P_ext and P_air: S^(z)_b = ${𝒮}$_z (P_tissue(V_L(P_ext, P_air)) − P^(z)_air). Merging with the equation on the first line of Equation (5), one can summarize the airway mechanics into a vectorial form:

$\begin{array}{cc}{𝒮}=H(P_{ext},P_{air})& \left(8\right)\end{array}$

with

$H{(P_{ext},P_{air})}^{(z)}=\{\begin{array}{ll}{{𝒮}}_z\left(P_{tissue}\left(V_L(P_{ext},P_{air})\right)-P_{air}^{(z)}\right)\hfill & \hfill \\ \text{ for }0\le z\le 16\hfill & \hfill \\ \alpha v_{adu}\left(P_{air}^{(z)},P_{ext}\right)/L_{ad}\hfill & \hfill \\ \text{ for }17\le z\le 22\hfill & \hfill \end{array}$

2.4. Air and Mucus Mechanics

In the previous section, we defined a model that mimics how airways and alveoli are behaving when they are submitted to an external pressure P_ext and an internal fluid pressure P_air. Whereas the external pressure P_ext is an input of the model, fluid pressures in the airways P_air are outputs. Fluid pressures and flows come from the interactions between airways deformations and fluid mechanics of air and mucus.

We will distinguish two lumen areas: the bronchi lumen area, whose surface is measured by the quantity S^(z)_b in generation z, and the air lumen area, whose surface is measured by the quantity S^(z)_a in generation z. The difference S^(z)_b − S^(z)_a represents the surface of the airway section filled with mucus.

2.4.1. Air Mechanics

In this section, we will derive the equations on the air lumen area surface S^(z)_a. Air is assumed incompressible and Newtonian. Its density is ρ_a = 1 kg.m⁻³ and its viscosity is μ_a = 1.8 10⁻⁵ Pa.s. We assume that air motion in a bronchus is at low regime with fully developed profiles as in Mauroy et al. (2011). The rate of air volume change in an airway is equal to the air flow getting in the airway—or getting out the airway, in which case we consider there is a negative airflow getting in the airway. Air flow getting in an airway in generation z comes from its two daughter branches standing in generation z + 1. The conservation law for air in generation z is:

$\begin{array}{cc}\frac{d{S}_a^{(z)}}{dt}{L}_b^{(z)}=-\left(2{\Phi }_a^{(z+1)}-{\Phi }_a^{(z)}\right)-n_{alv}^{(z)}{\Phi }_{alv}^{(z)}& \left(9\right)\end{array}$

The minus sign means that air flow coming in the tree is modeled as negative. Note that we could also have computed the balance with the air flow coming in from the parent branch instead of the daughter branch, however this would lead to the similar equations once rewritten at the tree level. The number n^(z)_alv corresponds to the number of alveoli connected to an airway in generation z, it is equal to 0 if the airway belongs to generations 0–16 and is equal to 58 if it belongs to a deeper generation, see geometry section upwards. Φ^(z)_alv corresponds to the part of the airway air flow that is induced by one alveolus belonging to generation z. Φ^(z)_alv is equal to the rate at which the volume v^(z)_alv changes in time. It is computed from Equation (5):

$\begin{array}{cc}\frac{d{v}_{alv}^{(z)}}{dt}=-{\Phi }_{alv}^{(z)}& \left(10\right)\end{array}$

Because of our hypotheses, the pressure drop per unit length in an airway, called $C=\frac{\partial p}{\partial z}$ depends on the generation index only. C^(z) in an airway of generation z can be computed from the air flow Φ^(z)_a in the bronchus and from air and mucus lumen areas, respectively denoted S^(z)_a and S^(z)_b:

$\begin{array}{cc}{C}^{(z)}=F({\Phi }_a^{(z)},S_a^{(z)},S_b^{(z)})& \left(11\right)\end{array}$

F computation is based on low regime and fully developed fluid dynamics hypotheses, as in our previous work (Mauroy et al., 2011). Details about how F is computed are recalled in Appendix A in Supplementary Material. Since pressure reference is that of atmospheric pressure, the pressure at the root of the tree (trachea) is zero. The pressure drop in an airway of generation z is C^(z)L^(z)_b and the pressure P^(z)_air in generation z is the sum of the pressure drops of the airways connecting the tree root (trachea) to a generation z airway:

$\begin{array}{cc}{P}_{air}^{(z)}=C^{(z)}\frac{L_b^{(z)}}{2}+{\displaystyle \sum _{g=0}^{z-1}{C}^{(g)}}{L}_b^{(g)}& \left(12\right)\end{array}$

The 1/2 in the first term means that P^(z)_air is air pressure at mid-bronchus length.

2.4.2. Mucus Mechanics

As in Mauroy et al. (2011), we assume that secretions are a Bingham fluid: it is solid if its inner shear stress σ is lower than a stress threshold σ₀ and liquid with viscosity μ_m otherwise. In this work, we assume σ₀ = 0.1 Pa and μ_m = 0.1 Pa.s, see (Mauroy et al., 2011). Secretions density is that of water, i.e., 1000 kg/m³. Mucus is motioned by air flow when the shear forces induced by air flow on mucus exceed the shear threshold σ₀. Mucus can have three states: either fully solid, fully liquid or both solid/liquid. The state of mucus depends on a radius r₀ = | 2σ₀ / C | where C is the pressure drop per unit length in the branch (Mauroy et al., 2011). Mucus is solid if r₀ > r_b, liquid if r₀ < r_a. If r_a < r₀ < r_b, then mucus is solid between radius r_a and r₀ and liquid between radius r₀ and r_b.

Mucus volume change in an airway in generation z is the balance between mucus input Φ^(z)_{m, I} and output Φ^(z)_{m, O}. Mucus volume in a generation z airway is the difference between airway total volume and the volume occupied by air: S^(z)_b L^(z)_b − S^(z)_a L^(z)_b. Since the airway length L_b is assumed constant, the time derivative of this equation leads to the conservation law for mucus:

$\begin{array}{cc}\left(\frac{d{S}_b^{(z)}}{dt}-\frac{d{S}_a^{(z)}}{dt}\right)L_b^{(z)}=|{\Phi }_{m,I}^{(z)}|-|{\Phi }_{m,O}^{(z)}|& \left(13\right)\end{array}$

The flow of mucus Φ^(z)_{m, O} that gets out of an airway through its interaction with air flow can be computed from air properties in the airway. Mucus flow Φ^(z)_{m, O} in an airway in generation z is a function of the pressure drop per unit length in the airway C^(z) and of the bronchi and air lumen areas of the airway:

$\begin{array}{cc}{\Phi }_{m,O}^{(z)}=G(C^{(z)},S_a^{(z)},S_b^{(z)})& \left(14\right)\end{array}$

Details on how to compute G are given in Appendix B in Supplementary Material, the method is the same as in Mauroy et al. (2011). Mucus flow in airways occurs only if the local shear stress induced by air flow on mucus overcomes the mucus yield stress σ₀. The flow of mucus that gets in a generation z airway comes either from the two daughter airways in generation z + 1 or from the parent airway in generation z − 1, depending on mucus flow direction. This is reflected by the sign of mucus flow:

$\begin{array}{cc}{\Phi }_{m,I}^{(z)}=\frac{1}{2}\max \left(0,{\Phi }_{m,O}^{(z-1)}\right)-2\min \left(0,{\Phi }_{m,O}^{(z+1)}\right)& \left(15\right)\end{array}$

If mucus flows Φ^{(z + 1)}_{m, O} in the two daughter airways (generation z + 1) are negative, then mucus is going up the tree and thus accounts twice for mucus inflow in an airway of generation z. If mucus in the parent airway (generation z − 1) is positive, then mucus in that airway is going down and accounts half for mucus inflow in an airway of generation z. If the flow in one of these neighboring airways are of the opposite sign as stated, then it does not account for the mucus inflow in an airway of generation z.

2.5. Solving Equations

The model consists in the set of Equations (8, 9, 10, 11, 12, 13, 14, 15). The system of equations is implemented numerically in C++ using a vectorial formulation based on Eigen 3.2.1 (http://eigen.tuxfamily.org/) and OpenMP (http://openmp.org/) multiprocessing. It is solved thanks to a full implicit first order method for differential Equations (9, 10, 13). A Newton method is used to solve Equation (7) in order to compute the function H, see Equation (8). Details of the mathematical and numerical implementations are given in Appendix C in Supplementary Material.

In order to be able to capture the effects of a 20 Hz oscillating external pressure, the typical time step in our simulations is 5.10⁻³ s. The time needed to compute one time step on 6 cores Xeon 5645 ranges from 0.1 to 0.2 s. In this study, the typical manipulation mimicked by our model covers about 200 s, with a total computation time of 4000–8000 s.

3. Results/Discussion

The model we describe in the previous sections allows to test how the changes in airway diameters affect air and mucus distributions in the tree. The model contains the minimal physics that we hypothesize to be main drivers of chest physiotherapy: the tree structure, linear air fluid mechanics, Bingham fluid mechanics to model secretions behavior and quasi-static deformable airways. In this section, we present the model predictions with pressure inputs mimicking chest physiotherapy. In Section 3.1, we define the initial state which is used in all the next simulations. This initial state is not representative of any specific pathology, it only reflects a lung with an excess of secretions near the eighth generation. In Section 3.2, we present how our model is responding to inputs that mimic manual chest physiotherapy. Next, in Section 3.3 we present how our model is responding to inputs that mimic high frequency chest wall oscillations. We study most particularly the role of an added static pressure on the chest. In order to be able to compare both the efficiency and comfort of these manipulations we define in Sections 3.4 and 3.5 two numbers that quantify manipulations efficiency and comfort.

3.1. Mimicking Chest Physiotherapy Manipulations

We will test different input pressures mimicking manual or mechanical chest physiotherapy. We recall that our model assumes the pressure to be homogeneous on the thorax. For each manipulation, we will focus on three quantity computed from the model outputs:

• the quantity of mucus expelled out of the tree v_out—the net reduction of mucus volume in the lung after a manipulation.

• the relative change in hydrodynamic resistance r during and after a manipulation—this quantity is correlated to the feeling for the patient to breath easily or not.

• the mean mucus position in the tree, that reflects how deep mucus stands in the lung. It is equal to

$mmp=\frac{-v_{out}+{\displaystyle {\sum }_{z=0}^N\left(z{2}^z(S_b^{(z)}-S_a^{(z)})L_b^{(z)}\right)}}{v_{out}+{\displaystyle {\sum }_{z=0}^N\left(2^z(S_b^{(z)}-S_a^{(z)})L_b^{(z)}\right)}}$

All numerical experiments are performed with a background chest motion that corresponds to an idealized normal ventilation of the patient. We model normal ventilation as a sinusoidal negative pressure applied on the thorax with a frequency of 0.2 Hz and an amplitude P_v = − 5 cmH₂O, see Figure 3A:

$P_{ext}^{ventil}(t)=P_v\frac{\left(1-\cos (2\pi t/5)\right)}{2}$

FIGURE 3

Figure 3. (A) External oscillating pressure to model normal rest ventilation. (B) Initial distribution of mucus in the tree: air is plotted in blue and mucus in green, bars height represents a proportion of airway surface at FRC. This distribution is not affected by normal rest ventilation.

With such data as an input, our model predicts a tidal volume of 0.5 L (6 L/min). This value for tidal volume is typical of rest regime ventilation (Weibel, 1984).

The initial mucus distribution in the tree is plotted on Figure 3B. It is homogenous in the first six generations and fills about ten percents of the bronchi lumen area. Then mucus thickness increases in the three next generations to reach its maximal thickness in generation eight, filling about fifty percents of the bronchi lumen area. Finally, mucus thickness decreases and reaches zero percent in generation sixteen. The initial configuration for mucus has been chosen in such a way mucus is almost not affected by normal rest ventilation. Other initial mucus distributions were tested and the results were not qualitatively affected by a change in the initial distribution.

3.2. Manual Physiotherapy

Manual chest physiotherapy is performed by applying pressures on the chest with hands. The pressures are applied during the expiration phase only and practitioners are guided by chest motion and reaction to their manipulations. Although the hand locations play an important role in chest physiotherapy, our model accounts only for pressures amplitude and time dependence: pressure is applied homogeneously onto the chest. We mimicked chest physiotherapy in our model by altering the ventilation signal during expiration, as shown on Figure 4A:

$P_{ext}^{manual}(t)=P_{ext}^{ventil}(t)+P_{cp}\max \left(\sin (2\pi t/5),0\right)$

FIGURE 4

Figure 4. (A) Chest pressure t → P^manual_ext(t) used as an input for our model to mimic manual chest physiotherapy. Pressure is assumed homogeneous on the thorax, case with P_cp = 20 cmH₂O. (B) Relative hydrodynamic resistance r variations during a manipulation with P_cp = 20 cmH₂O (blue). The red line corresponds to FRC relative hydrodynamic resistance before the manipulation (r_{|t = 0s} = 1.00), and the black line corresponds to FRC relative hydrodynamic resistance at the end of the manipulation (r_{|t = 230s} = 0.79).

Manipulation pressures are applied regularly during one session of 230 s, except for the 10 first and last seconds, where only normal ventilation occurs. The pressures applied are all the same during the whole manipulation. The 25 first seconds of pressure input are plotted on Figure 4A with P_cp = 20 cmH₂O. In manual chest physiotherapy, the pressure is not applied on the whole chest at once. The choice of the location for pressure application is made by the physiotherapist in order to optimize the pressure effect. Our model is for now not able to account for pressure spatial inhomogeneities and as a first approximation, we assume the hand pressure spreads over the whole chest homogeneously.

We investigated the role of pressure P_cp amplitude on mucus motion and distribution and on hydrodynamic resistance changes. Smaller chest pressures induce smaller air velocities and thus smaller shear stresses in the bronchi. As a consequence, the stress threshold that has to be overcome for secretions to be mobilized is only reached in highly obstructed bronchi. But the motion is quickly stopped when a new equilibrium between air flow and secretions distribution is reached. In particular, for chest pressures lower than 16.5 cmH₂O, we do not observe in the model any secretions going out of the lung, see Figure 5A. This does not mean however that secretions distribution has not been affected, as shown on Figure 5B. An important point is that by performing the manipulation, secretions are moving in such a way that their distribution always decreases the total hydrodynamic resistance of the tree, down to a value corresponding to an equilibrium between secretions and air flows. The first steps of the manipulations are the more efficient, as shown on Figure 4B. The higher the pressure, the lowest is the tree hydrodynamic resistance at the end of the manipulation. The first conclusion is then that the manipulations decrease the hydrodynamic resistance of the patient, which may lead to an improvement of patient breathing quality, even if no secretions get out of the tree.

FIGURE 5

Figure 5. (A) Volume of secretions expelled at the end of the manipulation as a function of the pressure amplitude. (B) Relative hydrodynamic resistance of the tree at the end of manipulation as a function of the pressure amplitude. (C) Mean secretions position in the tree as a function of the pressure amplitude. Mean position is expressed in generation index.

Secretions distribution in the tree are motioned from one generation to the next, and because manipulation pressures are applied during expiration, mucus tends to go up the tree, at least for moderate pressures amplitude, see Figure 5C. The mechanism is as follow: during pressure application, secretions are moving upward in the tree from one generation to the next, and thus fill the upper generations. Lung's recoil from the manipulation and lung opening due to inspiration are correlated to an inlet air flow in the lung, which potentially moves the secretions back into the lower generations. This backward secretions motion depends on how much secretions have accumulated during the pressure application. If the pressure is high, then secretions accumulate a lot in the upper bronchi, reducing drastically air lumen area, and by thus they are more affected by the inward air flow. This is why for high pressures, although there is less secretions in the tree at then end, they are globally deeper after the manipulation than before the manipulation. This phenomenon can be seen on Figure 6 where generations 11–13 have more secretions at the end of the manipulation.

FIGURE 6

Figure 6. Evolution of mucus distribution in the bronchi during the manipulation, with P_cp = 20 cmH₂O. The blue bars correspond to the surface of bronchi lumen area relatively to bronchi initial surface. The green bars correspond to the surface fraction of bronchi lumen area obstructed by secretions. At the end of manipulation, secretions have been globally motioned upward the tree, however some secretions have been pushed deeper in the lung, see generation 11–13.

3.3. Mechanical Physiotherapy

There exists a whole range of mechanical chest physiotherapy devices that are able to help secretions expectoration (Arens et al., 1994; Hristara-Papadopoulou et al., 2008; Mitchell et al., 2013). We will focus on devices that work by applying oscillatory positive pressures on the chest, a method that is referred to as high frequency chest wall oscillation (HFCWO). HFCWO covers different ways to apply the pressure on the thorax. Chest compression (CC) devices (Arens et al., 1994) apply a global static pressure on the chest, and low amplitude oscillations are applied on the whole thorax, thus working at low lung volume. Focused pulses technique (FPT) devices (Mitchell et al., 2013) do not apply a global static pressure but instead apply high amplitude oscillations on specific locations on the chest using pistons. In this section, we propose to mimic with our model the effects of these two types of devices. We chose an input external pressure that adds to normal ventilation a static pressure P_s and an oscillatory pressure with amplitude P_o and frequency f, see (Figure 7):

$P_{ext}^{oscill}(t)=P_{ext}^{ventil}(t)+P_s+\frac{P_o}{2}\sin (2\pi ft)$

FIGURE 7

Figure 7. Chest pressure used as an input for our model to mimic mechanical chest physiotherapy based on high frequency chest wall oscillations (20 Hz in this example). Pressure is assumed homogeneous on the thorax. (A) Pressure distribution on the thorax to mimic devices using chest compression (CC) using a static pressure. (B) Pressure distribution on the thorax to mimic devices using focused pulses technique (FPT) without static pressure (P_s = 0).

In the case of FPT, it is necessary to account for the pistons covering only a small fraction of the thorax. Thus, we use equivalent pressures P_o as inputs of the model. They represent the typical values applied by the pistons weighted by the surfaces ratio between pistons' heads and thorax. As for manual chest physiotherapy modeling, manipulation pressures are applied for one session of 230 s. The 10 first and last seconds are normal ventilation. The time range 10–15 s is used to start smoothly the manipulation by applying progressively the manipulation pressure. The time range 215–220 s is used to stop smoothly the manipulation by decreasing progressively the manipulation pressure.

3.3.1. An Example of High Frequency Chest Wall Oscillation

In this section, we study the predictions of the model with an input mimicking the effects of HFCWO with a static pressure P_s = 5.6 cmH₂O and an oscillatory pressure P_o = 1.2 cmH₂O at f = 20 Hz. Because of the static pressure, both ventilation and pressures oscillations are applied on a lung with low volume. Most particularly, this affects the net bronchi lumen area and the pleural pressure value and response (D'Angelo and Milic-Emili, 2005; Agostoni and Hyatt, 2011).

As for manual chest physiotherapy, our model predicts that the manipulation redistributes the secretions in the tree in such a way it reduces the hydrodynamic resistance of the tree by about 20 percents at the end of the manipulation, see Figure 8A. The manipulation efficiency is higher at the beginning and decreases with time, as shown on Figure 8B. Because the static pressure and the oscillations occurred during the whole ventilation cycle, the secretions spread a little more down the tree than for manual chest physiotherapy. Initially, mean mucus generation in the tree is 7.34, and it reaches 7.47 at the end of the manipulation, see Figure 9. With the parameters used in this section, our model predicts that no secretions are going out of the tree.

FIGURE 8

Figure 8. (A) Evolution of mucus distribution in the bronchi during the manipulation with P_s = 5.6 cmH₂O and P_o = 1.2 cmH₂O. (B) Relative hydrodynamic resistance changes during a manipulation with P_cp = 20 cmH2O (blue). The red line corresponds to FRC relative hydrodynamic resistance before the manipulation (r_{|t = 0s} = 1.00) while the black line corresponds to FRC relative hydrodynamic resistance at the end of the manipulation (r_{|t = 230s} = 0.79).

FIGURE 9

Figure 9. Role of static pressure P_s, with P_o = 1.2 cmH₂O and f = 20 Hz, data computed at the end of the manipulations for time t = 230s. (A) Volume of mucus getting out of the tree. (B) Relative hydrodynamic resistance of the tree. (C) Mean mucus position in the tree, in generation index. The red dashed line represents the value before the manipulation.

3.3.2. Role of Static Pressure P_s

We investigated the role of the static pressure P_s that controls the lung volume at which the manipulation is performed. Oscillating pressure P_o was fixed to 1.2 cmH₂O with a frequency f at 20 Hz. Surprisingly, we found that static pressure has only small effects on the final hydrodynamic resistance and on the mucus spreading in the tree. Increasing the static pressure leads to a slight improvement of the hydrodynamic resistance at the end of the manipulation, see Figure 9B. Larger static pressure also sends secretions a bit deeper in the tree, and increases slightly the mean mucus generation at the end of the manipulation, see Figure 9C. Finally, increasing static pressure reduces to amount of secretions getting out of the tree, see Figure 9A.

These results indicate that, in the frame of our model, increasing the static pressure may have no real effect. Indeed, increasing the static pressure goes with a constriction of the bronchi, which increases the hydrodynamic resistance of the bronchial tree, as shown in the previous example on Figure 8B. The air flow created by the oscillating pressure is then smaller if the static pressure is larger. Shear forces determine the motion of the secretions, and in a bronchi, they are proportional to the ratio between the air flow rate in the bronchi and the diameter of the bronchi at the power three. Our results suggest that the increase on shear stress gained by smaller airways diameters is almost compensated by the loss in term of air flow rate.

3.3.3. Influence of the Oscillating Pressure Amplitude P_o and Frequency f

Air flow in a bronchi determines the shear stress applied on the secretions and thus their motion. Air flow in a bronchi is partially related to velocity of the bronchi volume changes, see Equation (9). Thus, secretions motion is correlated to bronchi volume change velocity. This quantity is actually the ratio between the amount of volume change over the time during which this volume change occurs. The amount of volume change is given by oscillating pressure amplitude, while the time for these changes to occur is given by the oscillating pressure frequency. Thus, these two quantities are ought to play an important role on manipulation efficiency.

We investigated first the role of the oscillating pressure P_o. Because static pressure has only a small effect in our model, we assume P_s = 0.6 cmH₂O which makes the total pressure applied on the thorax oscillate between 0 and 1.2 cmH₂O. The frequency f is kept at 20 Hz. As predicted, the larger the amplitude of the oscillations, the more efficient is the manipulation, see Figure 10. As before, a more efficient mechanical manipulation tends to send more mucus down the tree, as shown on Figure 10C, because oscillations are symmetric and applied during the whole cycle of ventilation.

FIGURE 10

Figure 10. Role of oscillating pressure P_o, with P_s = 0.6 cmH₂O and f = 20 Hz, data computed at the end of the manipulations for time t = 230s. (A) Volume of mucus getting out of the tree. (B) Relative hydrodynamic resistance of the tree. (C) Mean mucus position in the tree, in generation index. The red dashed line represents the value before the manipulation.

The frequency f also plays an important role on the manipulation efficiency. There is a minimal frequency under which the manipulations have low efficiency. Once this minimal frequency is reached, the hydrodynamic resistance and the mean mucus position in the tree are only slightly affected by any supplementary increase of the frequency, as shown on Figure 11. Only the quantity of secretions that get out of the tree continue to increase regularly with the frequency, as shown on Figure 11A, mainly because the secretions in the root of the tree (mimicking the trachea) are all the more mobilized than the frequency is high.

FIGURE 11

Figure 11. Role of oscillating pressure frequency f, with P_s = 0.6 cmH₂O and P_o = 1.2 cmH₂O, data computed at the end of the manipulations for time t = 230s. (A) Volume of mucus getting out of the tree. (B) Relative hydrodynamic resistance of the tree. (C) Mean mucus position in the tree, in generation index. The red dashed line represents the value before the manipulation.

3.4. Shrek Number: A Measure for Chest Physiotherapy Efficiency

The mucus behavior in a generation i can be predicted with the ratio between the air constraints on the mucus measured by the air shear stress τ_i in an airway and the shear stress threshold σ₀ of the mucus. In an airway with radius r_i and length l_i, the air shear stress is τ_i = 4 μ_a Φ⁽ⁱ⁾_a / (π r³_i). We define the instantaneous Shrek number Sh as the average of the ratios τ_i / σ₀ over the tree generations:

Instantaneous Shrek (DreamWorks Animation SKG, Shrek, 2001) number:

$\begin{array}{cc}iSh=\frac{1}{N}{\displaystyle \sum _{i=0}^{N-1}\frac{4{\mu }_a{\Phi }_a^{(i)}}{\pi r_i^3{\sigma }_0}}& \left(16\right)\end{array}$

The instantaneous Shrek number can be approximated with a formula easier to evaluate, using lung hydrodynamic resistance R_aw and air flow at mouth level Φ_a, see details in Appendix E in Supplementary Material. This approximation is based on the hypothesis that the ratios of length over diameter of all the airways are equal to 3, which corresponds to the mean value measured (Tawhai et al., 2004):

$\begin{array}{cc}iSh~\frac{R_{aw}{\Phi }_a}{12N{\sigma }_0}& \left(17\right)\end{array}$

The larger the instantaneous Shrek number, the more secretions are mobilized. An instantaneous Shrek number smaller than 1 does not mean that mucus is not mobilized, since Shrek number represents an average over the generations. It means however that in many generations, the mucus is probably not mobilized. However, it is important to notice that airways impaired by mucus will have the highest Shrek numbers in the tree. Because chest physiotherapists are modulating the flow in their manipulation to keep mucus moving in the tree, they are indirectly and intuitively modulating the instantaneous Shrek number to keep it high enough to get a response from mucus. These modulations are paired with the real-time changes of lung hydrodynamic resistance.

Finally, the instantaneous Shrek number can be reformulated as a function of the homogeneous external pressure used in our model. The new formulation involves quantities related to lung function and to manipulation: the compliance of the lung C_L that relates the external pressure to the lung volume change, and the velocity with which changes in external pressure occurs ΔP_ext / Δt:

$\begin{array}{cc}iSh~\frac{C_L{R}_{aw}}{12N{\sigma }_0}\frac{\Delta P_{ext}}{\Delta t}& \left(18\right)\end{array}$

This formula comes from the simple fact that air flow in the lung equals lung volume change ΔV_L divided by the duration for this change Δt. Also the compliance of the lung checks C_L = ΔV_L / ΔP_ext. Combining these equations leads to the alternative formulation for instantaneous Shrek number iSh in Equation (18).

To be able to compare different chest physiotherapy manipulations in the hypotheses of our model, we compute the Shrek number which is the time average of the instantaneous Shrek number. It defines a global mean efficiency for the manipulation, with the limit that any average may have. If the manipulation has a duration T, then:

$\begin{array}{cc}Sh=\frac{1}{T}{\displaystyle {\int }_0^{T}i}Sh(t)dt& \left(19\right)\end{array}$

The mean Shrek number associated to normal ventilation in a normal lung is about 0.2.

It is important to highlight that Shrek number is only an indicator which aggregates many information, thus its use is correlated to an understanding of the physical phenomena involved in chest physiotherapy and to an understanding of its limitation due to the fact that it is an average.

On Figure 12, we compare the efficiency in term of hydrodynamic resistance decrease for all the manipulations we studied in the previous sections. Most particularly, we can see that two manipulations with the same Shrek number have globally the same effect on hydrodynamic resistance. As expected, modulating the static pressure in CC techniques leads to very similar Shrek numbers, while manual and FPTs are able to cover a wide range of Shrek number by modulating their respective pressure amplitude.

FIGURE 12

Figure 12. Relative reduction of lung hydrodynamic resistance after 230 s manipulations (model predictions) as a function of the time averaged Shrek number.

3.5. Comfort Number: A Measure for Chest Physiotherapy Comfort

The instantaneous comfort number of a manipulation is the lung tissue pressure shift relatively to normal ventilation. If at a given time of a manipulation, the volume of the lung is V_L while it would have been V^ventil_L without the manipulation, then the instantaneous comfort number is:

$\begin{array}{cc}iCom=\left|\frac{P_{tissue}(V_L)-P_{tissue}(V_L^{ventil})}{P_{tissue}(V_L^{ventil})}\right|& \left(20\right)\end{array}$

The instantaneous comfort number is zero if only ventilation occurs. The higher the number, the larger are the mechanical constraints in the lung tissues. Notice that this number only measures the mechanical stress inside the lung and is not directly related to its hydrodynamic resistance. For example if lung tissues are suffering high negative pressure because of an high inflation, the associated lung hydrodynamic resistance is probably low but mechanical stress is high.

The comfort number of a manipulation with duration T is the mean value computed over the whole manipulation,

$\begin{array}{cc}Com=\frac{1}{T}{\displaystyle {\int }_0^{T}i}Com(t)dt& \left(21\right)\end{array}$

As a consequence, a strong pressure applied for a short time can lead to similar comfort number as a moderate or even small pressure applied for a long time.

On Figure 13, the comfort number for different idealized chest physiotherapy techniques studied in this paper are plotted. CC techniques with moderate static pressure are associated to high comfort numbers, mainly because the static pressure is applied constantly during the manipulation. FPT are associated to low comfort numbers, because they are applying low pressures on chest. Notice that comfort number is only related to the net compressive forces induced into the lung by the manipulations, comfort number is not able to measure the role the frequency of pressure application may play on the patient comfort.

FIGURE 13

Figure 13. Comfort number computed from model predictions for the different chest physiotherapy techniques mimicked in this study. The higher the comfort number, the larger are the stress in the lung tissues. The comfort number for normal ventilation is zero.

We recall that in our model the pressure is applied onto the whole chest. Thus, in the case of manual chest physiotherapy, comfort number has to be adapted to reflect the fact that the pressure is applied only on a subpart of the chest. In a first approximation, it can be done by weighting the pressure used to compute comfort number with the ratio between the surface on which the pressure is applied over the whole chest area.

3.6. Model Limitations

Because the model described in this work contains the biophysical phenomena we identified as major in chest physiotherapy (“minimal” model), we are confident that it is able to predict main behaviors of chest physiotherapy. However, our predictions can only be qualitative or comparative since the biophysical phenomena included are simplified versions of the “real” phenomena and are based on static mechanics. For example, secretions are modeled by a simple Bingham fluid, or secretions thickness is assumed homogeneous in a bronchus. Moreover, quantitative predictions would need to include in the model some other phenomena which, although not strongly driving the system, might play a role on the quantitative data. Thus, any interpretation based on the results of our model has to be made in the frame of the model's hypotheses.

Our model is built from static mechanical data, such as Agostoni et al. pressure-volume curves (Agostoni and Hyatt, 2011) or (Lambert et al., 1982) mechanical behavior of the bronchi lumen areas. Our model shows that unsteady effects might play an important role in chest physiotherapy efficiency, so they will have to be accounted for in future work in order to improve the accuracy of our model predictions. Unfortunately, to our knowledge, no corresponding data is available in the literature.

Our model is intrinsically non-linear. Non linearity in the model comes from two different aspects. First, in the way we deal with the mechanics of the lung and thorax. For this aspect, our model does not rely on a modeling process but instead on -non linear- data about lung mechanical inner pressure and about lung volumes [pressure curves from Agostoni et al. (Agostoni and Hyatt, 2011)]. Secondly, non-linearity arises from the fluid mechanics of mucus which is based on Bingham fluid. Bingham fluid flows only if inner shear stress is above a threshold. Although non-linear, secretion mechanics remains a basic simplification of real secretion mechanics, which behavior may be highly dependent on the local physical and biological environment. Air fluid mechanics in our model is linear, and we ignored inertia and turbulence which may break the homogeneity of flows distribution in the tree, even if the tree has symmetric branching (Mauroy et al., 2003). We also chose to work with fully developed velocity profiles for three reasons: 1/ it is a reasonable first order approximation; 2/ it allowed to compute quasi-analytical expressions for otherwise costly to compute variables; 3/ it eased in the interpretation of mucus motion. However, flow oscillations in the bronchi and their interaction with deformable walls might lead to non-fully developed velocity profiles in bronchi where Womersley number is high enough (Womersley, 1955) and thus might affect locally mucus motion.

Asymmetric flow distribution can also be a consequence of an asymmetry in the geometry, either induced by asymmetric branching or by asymmetric secretions distribution in the tree. This last phenomenon may play a crucial role when a specific bronchus is more occluded by secretions than its neighbors: air flows might be redistributed toward the lowest resistive bronchial pathways and possibly miss the bronchi where secretions quantity is the highest. Chest physiotherapists are able to counteract this phenomena by applying pressures in a non-homogeneous manner on the thorax. By assuming in our model that pressure is applied homogeneously on the thorax, none of these effects can be predicted with our current model.

Chest physiotherapy is highly patient dependent: manipulation are based on real-time response to the praticien feeling (hand feeling, noises, mouth air flows, etc.). The data we used in this model for secretions properties (Lai et al., 2009), pressure-volume relationships from Agostoni et al. (Agostoni and Hyatt, 2011) and lung geometry and mechanics (Lambert et al., 1982) reflect either mean or specific behaviors and cannot be applied as such to a specific patient. It has to be clear that the predictions of our model are general predictions that cannot be mapped as such to any patient. Again, careful interpretations in the frame of the model hypotheses have to be made before reaching any specific conclusion.

We proposed two numbers, the Shrek number that measures an efficiency for the manipulation and the comfort number that measures a degree of compressive discomfort correlated to a manipulation. Both these numbers have been tested with our idealized model. They are however not validated in a medical way and there significance remains as of today an hypothesis. In this paper we propose these numbers as candidates for evaluating chest physiotherapy manipulations, but clinical measurements have to be made before any medical applications of these numbers are made. This leads to another question: how to evaluate quantitatively the values for these numbers in a clinical frame. Concerning Shrek number, we proposed a formula that allows to compute its value with typical medical data, see Equation (18) and Appendix E in Supplementary Material. Concerning comfort number, it is necessary to get information about patient pressure-volume lung curves, which may be difficult in the frame of everyday consultations. This may be solved by assuming that patients pressure-volume lung curves have globally the same shapes and can be fitted through the measurement of a few simple physiological parameters at the beginning of a manipulation. This aspect is however out of the scope of this paper since it can be made only in a clinical frame.

4. Conclusion

We developed a model of the thorax and lung that is able to link the pressure on the thorax with the air and secretions motion in the airways. The goal of this model is to give indications on how chest physiotherapy may interact with the biomechanics of the thorax and of the lung. The model we built in this work is based on core phenomena that are involved in chest physiotherapy and can predict qualitative behaviors and compare different techniques, in the limitations of the model's hypotheses.

Our results indicate that manipulations need to overcome secretions threshold in order to be able to mobilize secretions. Our model shows that manipulations main effect is to reduce the hydrodynamic resistance of the airway tree by secretions redistribution in the tree, eventually reaching a distribution that is not anymore affected by the manipulation. This effect means that patient breathing might be instantly improved by the manipulation. The second main effect is mucus expectoration, which is also another mean to improve the status of the patient. Mucus expectoration predicted by our model is non-negligible only for pressures high enough, and also, in the case of HFCWO, for frequencies high enough. These conclusions, essentially mechanistic, do not depend qualitatively on the initial mucus distribution.

All these phenomena are driven by the fact that the motor for secretions motion is the shear stress transmitted by air to the secretions. In a bronchus the shear stress is proportional to the air flow amplitude and to the inverse of the bronchus radius. Since air flow is created by applying pressure on the chest, the same method that is used to create the air flow will increase the hydrodynamic resistance of the airway tree, thus negatively affecting the air flow amplitude. Consequently, at some point, using a higher pressure will not really improve the efficiency of the manipulation, as shown on Figure 5. This phenomenon is well-known by chest physiotherapists and is also at the origin of forced expiration curve shape (Lambert, 2004).

Finally, we proposed two numbers to measure the efficiency of the manipulations: the Shrek number and the comfort number. They are both performing well in the case of our idealized manipulations. They need however to be validated with clinical studies.

Conflict of Interest Statement

We, the authors, declare that the financial support by the private company RespInnovation SAS (Seillans, France) did not affect in any way the results presented in this work. The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to thank RespInnovation SAS and the CNRS PEPS INSIS funding program Ecoulements Biologiques for supporting this work.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys.2015.00214

References

Agostoni, E., and Hyatt, R. E. (2011). “Static behavior of the respiratory system,” in Comprehensive Physiology, ed R. Terjung (Hoboken, NJ: John Wiley & Sons, Inc.), 113–130.

Arens, R., Gozal, D., Omlin, K. J., Vega, J., Boyd, K. P., Keens, T. G., et al. (1994). Comparison of high frequency chest compression and conventional chest physiotherapy in hospitalized patients with cystic fibrosis. Am. J. Respir. Crit. Care Med. 150, 1154–1157.

PubMed Abstract | Google Scholar

Basser, P. J., McMahon, T. A., and Griffith, P. (1989). The mechanism of mucus clearance in cough. J. Biomech. Eng. 111, 288–297.

PubMed Abstract | Google Scholar

Beauvois, E., Fouré, H., Gouilly, P., and Reychler, G. (2007). Kinésithérapie dans la bronchiolite: doute raisonné ou raison de douter? Kinésithérapie 7, 51–52. doi: 10.1016/S1779-0123(07)70357-1

CrossRef Full Text

Bokov, P., Mauroy, B., Mahut, B., Delclaux, C., and Flaud, P. (2014). Homothety ratio of airway diameters and site of airway resistance in healthy and COPD subjects. Respir. Physiol. Neurobiol. 191, 38–43. doi: 10.1016/j.resp.2013.10.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Bokov, P., Mauroy, B., Revel, M.-P., Brun, P.-A., Peiffer, C., Daniel, C., et al. (2010). Lumen areas and homothety factor influence airway resistance in COPD. Respir. Physiol. Neurobiol. 173, 1–10. doi: 10.1016/j.resp.2010.05.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Clément, R., Douady, S., and Mauroy, B. (2012). Branching geometry induced by lung self-regulated growth. Phys. Biol. 9:066006. doi: 10.1088/1478-3975/9/6/066006

PubMed Abstract | CrossRef Full Text | Google Scholar

D'Angelo, E., and Milic-Emili, J. (2005). “Statics of the respiratory system,” in Physiologic Basis of Respiratory Disease, eds Q. Hamid, J. Shannon, and J. Martin (Hamilton, ON: BC Decker Inc.,) 15–25.

Florens, M., Sapoval, B., and Filoche, M. (2011). Optimal branching asymmetry of hydrodynamic pulsatile trees. Phys. Rev. Lett. 106:178104. doi: 10.1103/PhysRevLett.106.178104

PubMed Abstract | CrossRef Full Text | Google Scholar

Grotberg, J. B. (2001). Respiratory fluid mechanics and transport processes. Annu. Rev. Biomed. Eng. 3, 421–457. doi: 10.1146/annurev.bioeng.3.1.421

PubMed Abstract | CrossRef Full Text | Google Scholar

Haefeli-Bleuer, B., and Weibel, E. R. (1988). Morphometry of the human pulmonary acinus. Anat. Rec. 220, 401–414.

PubMed Abstract

Hristara-Papadopoulou, A., Tsanakas, J., Diomou, G., and Papadopoulou, O. (2008). Current devices of respiratory physiotherapy. Hippokratia 12, 211–220.

PubMed Abstract | Google Scholar

King, M. (2006). Physiology of mucus clearance. Paediatr. Respir. Rev. 7, S212–S214. doi: 10.1016/j.prrv.2006.04.199

PubMed Abstract | CrossRef Full Text | Google Scholar

Lai, S. K., Wang, Y.-Y., Wirtz, D., and Hanes, J. (2009). Micro- and macrorheology of mucus. Adv. Drug Deliv. Rev. 61, 86–100. doi: 10.1016/j.addr.2008.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Lambert, R. K. (2004). Model of forced expiratory flows and airway geometry in infants. J. Appl. Physiol. 96, 688–692. doi: 10.1152/japplphysiol.00843.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Lambert, R. K., Wilson, T. A., Hyatt, R. E., and Rodarte, J. R. (1982). A computational model for expiratory flow. J. Appl. Physiol. 52, 44–56.

PubMed Abstract | Google Scholar

Langenderfer, B. (1998). Alternatives to percussion and postural drainage. A review of mucus clearance therapies: percussion and postural drainage, autogenic drainage, positive expiratory pressure, flutter valve, intrapulmonary percussive ventilation, and high-frequency chest compression with the ThAIRapy Vest. J. Cardiopulm. Rehabil. 18, 283–289.

PubMed Abstract | Google Scholar

Main, E., Prasad, A., and Schans, C. (2005). Conventional chest physiotherapy compared to other airway clearance techniques for cystic fibrosis. Cochrane Database Syst. Rev. 1, CD002011. doi: 10.1002/14651858.cd002011.pub2

PubMed Abstract | CrossRef Full Text | Google Scholar

Mauroy, B., and Bokov, P. (2010). The influence of variability on the optimal shape of an airway tree branching asymmetrically. Phys. Biol. 7:16007. doi: 10.1088/1478-3975/7/1/016007

PubMed Abstract | CrossRef Full Text | Google Scholar

Mauroy, B., Fausser, C., Pelca, D., Merckx, J., and Flaud, P. (2011). Toward the modeling of mucus draining from the human lung: role of the geometry of the airway tree. Phys. Biol. 8:056006. doi: 10.1088/1478-3975/8/5/056006

PubMed Abstract | CrossRef Full Text | Google Scholar

Mauroy, B., Filoche, M., Andrade, J. S., and Sapoval, B. (2003). Interplay between geometry and flow distribution in an airway tree. Phys. Rev. Lett. 90:148101. doi: 10.1103/PhysRevLett.90.148101

PubMed Abstract | CrossRef Full Text | Google Scholar

Mauroy, B., Filoche, M., Weibel, E. R., and Sapoval, B. (2004). An optimal bronchial tree may be dangerous. Nature 427, 633–636. doi: 10.1038/nature02287

PubMed Abstract | CrossRef Full Text | Google Scholar

Mauroy, B., and Meunier, N. (2008). Optimal Poiseuille flow in a finite elastic dyadic tree. Math. Model. Numerical Anal. 42, 507–533. doi: 10.1051/m2an:2008015

CrossRef Full Text

Mitchell, B., Lupi, F., Gianangeli, P., Poldi, G., Malerba, G., and Pontesilli, M. (2013). Hfcwo Vest. US Patent App. 13/521,495.

DreamWorks Animation SKG, Shrek. (2001). Shrek Use his Own Secretions as Core Material for Building Everyday Objects. Glendale, CA.

Ochs, M., Nyengaard, J. R., Jung, A., Knudsen, L., Voigt, M., Wahlers, T., et al. (2004). The number of alveoli in the human lung. Am. J. Respir. Crit. Care Med. 169, 120–124. doi: 10.1164/rccm.200308-1107OC

PubMed Abstract | CrossRef Full Text | Google Scholar

Postiaux, G., and Lens, E. (1992). De ladite ≪Accélération du Flux Expiratoire (AFE)≫: où Forced is… Fast (Expiration technique-FET). Ann. Kinésithér 19, 411–427.

Google Scholar

Pryor, J. A. (1999). Physiotherapy for airway clearance in adults. Eur. Respir. J. 14, 1418–1424.

PubMed Abstract | Google Scholar

Tawhai, M. H., Hunter, P., Tschirren, J., Reinhardt, J., McLennan, G., and Hoffman, E. A. (2004). CT-based geometry analysis and finite element models of the human and ovine bronchial tree. J. Appl. Physiol. 97, 2310–2321. doi: 10.1152/japplphysiol.00520.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Schans, C. P. (2007). Bronchial mucus transport. Respir. Care 52, 1150–1158.

PubMed Abstract | Google Scholar

Vandevenne, A. (1999). Rééducation Respiratoire: Bases Cliniques, Physiopathologiques et Résultats. Paris: Masson.

Weibel, E. R. (1984). The Pathway for Oxygen: Structure and Function in the Mammalian Respiratory System. Cambridge, MA: Harvard University Press.

Google Scholar

Whimster, W. F. (1970). The microanatomy of the alveolar duct system. Thorax 25, 141–149.

PubMed Abstract | Google Scholar

Womersley, J. R. (1955). Oscillatory motion of a viscous liquid in a thin-walled elastic tube—I: the linear approximation for long waves. Lond. Edinb. Dublin Philos. Mag. J. Sci. 46, 199–221.

Google Scholar

Zahm, J. M., King, M., Duvivier, C., Pierrot, D., Girod, S., and Puchelle, E. (1991). Role of simulated repetitive coughing in mucus clearance. Eur. Respir. J. 4, 311–315.

PubMed Abstract | Google Scholar