Structural changes in the COPD lung and related heterogeneity

This paper proposes a mathematical framework for understanding how the structural changes in the COPD lung reflect in model parameters. The core of the analysis is a correlation between the heterogeneity in the lung as COPD degree changes (GOLD II, III and IV) and the nonlinearity index evaluated using the forced oscillation technique. A low frequency evaluation of respiratory impedance models and nonlinearity degree is performed since changes in tissue mechanics are related to viscoelastic properties. Simulation analysis of our model indicates a good correlation to expected changes in heterogeneity and nonlinear effects. A total of 43 COPD diagnosed patients are evaluated, distributed as GOLD II (18), GOLD III (15) and GOLD IV (10). Experimental data supports the claims and indicate that the proposed model and index for nonlinearity is well-suited to capture COPD structural changes.


Introduction
Fractional calculus (FC) penetrated well the biomedical engineering field of research and numerous literature reviews witness its success [1,2]. It has been shown that both time domain and frequency domain analysis can be performed with various FC tools. Already several decades ago, modelling lung mechanics have been subject to fervent applications of constantphase element models, a successful tool emerging from fractional calculus [3][4][5]. Special attention has been given to viscoelastic properties, intrinsically requiring such models and providing relevant clinical information upon changes in lungs with disease [6][7][8].
Geometry of the lungs, either fractal or heterogeneous, associated with CT scans, has been employed to analyse images of lung tissue [9]. Since single fractal models of image processing algorithms for clinical analysis have been poor in capturing the high degree of variability in airway density, multi-fractal tools have been employed [9]. Same FC tools have been employed for image processing of other biological tissues as well [9,10].
In our previous work, we employed structure and mechanical properties to develop consistent impedance models for respiratory tree [11][12][13]. We have shown that such models converge to lumped impedance models with fractional order terms to account for frequency PLOS [13,14]. We expected that changes from airway remodelling with disease affecting morphology and structure would lead to changes in model parameters; this hypothesis has been verified in both simulations as well as in clinical studies. Changes in chronic obstructive pulmonary disease (COPD) are affecting mainly distal airways and lung parenchyma, involving mechanisms related to small airway disease [15]. Alveolar structure and morphology is drastically affected and thus plays an important role in determining the viscoelastic properties of lung parenchyma. Consequently, the diffusion mechanism is greatly dependent on these changes.
The goal of our study is to analyse what are the effects induced by these changes at alveolar levels upon impedance model parameters. We expect that since the structure plays a role in determining the appearance of a fractional order term in the impedance, this term may change as a result of changes due to disease. The core of the paper is to evaluate if there is a correlation between the heterogeneity in the lungs (as COPD degree progresses) and the non-linear index.
The paper is organized as follows: In Section II the theoretical background and the proposed theoretical framework is presented. In this section the case of one alveolus and the case of a network of alveoli is discussed. Next, the structural changes in patients with chronic obstructive disease are presented, as well as a short description of the device used to acquire the data. In section III the simulation and experimental results are shown and discussed followed by the conclusions.

Theoretical background
In the context of electrical analogy, one usually considers voltage as the equivalent for respiratory pressure (P) and current as the equivalent for airflow (Q). The pressure in function of time is the accumulation of these components with the trans-pulmonary pressure at end-expiration (P 0 ): where elastance (1/C) (the reciprocal of the compliance) is directly related to volume (V), the resistance (R) is related to airflow and the inertia (I) is related to accelerating particles. This equation represents the global relationship between airflow and pressure in the lungs. Its equivalent frequency domain form is the respiratory impedance as a function of oscillatory frequency [16]. This principle has been applied by means of Womersley theory and Navier-Stokes equations to elastic tubes preserving structure and morphology of the airways [11]. According to the seminal work of Weibel and Mandelbrot [17][18][19], the structure of the lungs can be roughly approximated by a recurrent dichotomous tree of 24 levels, leading to a rough approximation of 2 24 alveolar ducts. At the end of each alveolar duct there is a porous bag called alveoli, which are grouped together like a lot of interlinked caves, rather than individual sacks.
A mathematical model of the respiratory tree up to 24th level has been proposed [11] and used in subsequent works in clinical studies [13]. The model is based on the geometrical parameters: radius and wall thickness (r, h), on the mechanical characteristics of the airway tube: complex elastic moduli (given by its modulus (|E|) and angle(φ E )) on the Poisson coefficient (ν P ), and on the air properties: viscosity and density (μ, ρ). Over the length ℓ of an airway tube, we have the corresponding definitions for electrical resistance: and for electrical capacitance where, M 1 is the modulus of the Bessel function of order 1, ε 1 is the phase angle of the complex Bessel function of order 1, and d ¼ ffiffiffiffiffiffiffiffiffiffiffi ffi or=m p the Womersley parameter. We have shown that the geometry of the lungs leads to a piecewise fractal structure, whose electrical equivalent leads to a recurrent ladder network. The admittance of such network has the form: with R e1 and C e1 denoting the resistance and compliance in the first airway (i.e. trachea), and λ the ratio of resistances per total levels, χ the ratio of compliances per total levels and s the Laplace operator. This continuous fraction expansion can be well-approximated for N ! 1 airway levels by the compact form of admittance [20]: with K(λ, χ) a gain factor depending on the values of the ratios and the fractional order n given by Notice that the impedance, Z N (s) will be the inverse of the admittance. In the frequency domain, the fractional order will lead to a constant-phase behaviour, i.e. a phase-locking in the frequency range given by the convergence conditions [13,20]. The lumped version of this model is the well-known constant-phase element model given by: where R (kPa s/L) and I (kPa s 2 /L) denote the central airway resistance and inertance respectively, whereas the last term consists of a constant-phase element which can be split into a real and imaginary part, denoting tissue damping and tissue elastance. These two latter properties are defined as averaged over the frequency range evaluated by the lung function test [3][4][5].
and their ratio denotes the degree of heterogeneity present in the tissue [21]. This ratio is in fact a characterization of sources for nonlinear contributions in the respiratory dynamical properties. If changes in R are small, then any changes in G r (kPa s (1−β) /L) will represent changes in parenchyma or in very small airways. The memory effects are introduced via s (1−β) [22]. These effects are observed in materials with viscoelastic properties. Chronic changes in H r (kPa s (1−β) /L) reflect changes in the intrinsic mechanical properties of the parenchyma [23]. This parametric model has been shown to reliably estimate airway and tissue properties [24], and its sensitivity to bronchodilatation in dogs at low frequencies has been assessed in [3]. The same model has been used to assess respiratory impedance in rats at low frequencies in [25]. The model has been used to evaluate respiratory impedance in healthy and sick patients [13,14,21]. The reliability of the parameter estimates from model fitting to raw data has shown that the inertance I and the resistance R may contain a high degree of uncertainty, while the constant-phase term delivers a reliable estimate of peripheral tissue properties. The study also showed that the steep dependency on frequency at low frequencies in the real part of impedance is consistently fitted by the constant-phase model from Eq (7) when compared to other models from literature. In order to quantify the non-linear contributions the following index has been introduced [13]: where each variable represents the sum of the absolute values of all the contributions in the pressure and flow signals a non-excited frequencies (even and odd) and the odd excited frequencies. The index expresses the relative ration of the contributions at non-excited frequencies with respect to the contributions at the excited frequencies. For a detailed overview on filtering techniques applied to extract the contributions see [13,26].

Proposed theoretical framework
Work of breathing in the lungs takes place with loss of kinetic energy and heat loss (production of heat through friction of air in airways and airductuli). To allow a systematic approach, we need to use an impedance model of one alveolus which will be later used for analysis of alveolar sac structures and their respective properties in presence of remodelling effects.

Impedance of one alveolus.
Let us assume a spherical alveolus characterized by an inner radius r 1 and an outer radius r 2 . We may assume without loss of generality that heat conduction is anisotropic, heat diffusion is unidirectional and spherical coordinates can be used for modelling. Following the line of thought presented in [22,27], approximating the spherical wall (i.e. difference between inner radius and outer radius) with many cells whose pressure is governed by: with R i and C i representing the resistance and compliance properties of the cell i and P i pressure profile assumed uniformly distributed along the cell face. The thickness of each cell is then given by the relation: h = r 2 − r 1 /N, with N the total number of cells and 0 i N. Introducing N 0 = r 1 /h, it follows that the radii of each cell is given by The surface of the cell becomes: and it follows that the resistance of each cell is given by: where λ is the thermal conductivity (W m −1˚C−1 ). The capacity per cell is defined as: with ρ mass density (kg m −3 ) and c specific heat (J kg −1˚C−1 ). The ladder network of series connected cells with resistance and capacity calculated in Eqs (12) and (13) in Laplace operator, provides the explicit impedance which can be evaluated over a range of frequencies. This can be either a ladder network made of recurrent elements, either one with same element values per cell. In our case, the latter applies since we assumed a homogeneous density of the alveolar wall. The lumped transfer impedance with heat loss of such a ladder network of this spherical alveolus is given by: with s the Laplace operator, a = λ(ρ Á c) thermal diffusivity (m 2 s −1 ). For specific values of lung tissue and frequencies (s ! jω) one may evaluate this lumped impedance. This lumped impedance has an operator at non-integer order, implying a constant phase of 45 degrees.

Network of alveoli.
Air in the lungs passes terminal bronchi, enters respiratory bronchi and consequently reaches the alveolar ducts and alveoli. This may be interpreted as penetration of air through a porous material, i.e. lung tissue, whose density increases with distance. Fig 1 illustrates the aforementioned transport phenomenon. The material is assumed to be homogeneous, sponge-like, whereas every alveolus has same properties. This assumption corresponds to healthy lung tissue.
Fractal properties have been applied to intergranullar materials for impedance models and sphere like geometry has been employed successfully to proof convergence to a fractional order impedance model [28,29]. The structure of such porous materials has been modelled with FC tools and led to fractal dimensions between dimension 2 (surface) and 3 (volumetric). Depending on the density of the material, the values for fractal dimension changed accordingly. In the context of lung pathology, changes in density occur in emphysematous lungs. The fractional order impedance term values reflects these changes and thus we speculate it might be of interest to study this aspect.
To illustrate the effect of changes in structure in tissue lungs, we will first discuss the nominal conditions of healthy lung tissue. The tools we employ are similar as for one alveolus, equivalent ladder network models. Porous material is not necessary permeable, therefore the pores we consider in our model are not linked to each other, thus similar to alveolar sacs (dead-end). To enter these pores, the airflow must counteract resistance with dissipative energy and must inflate/deflate with an elastic capacitance which decreases the peak pressure for the same given volume of air. The flow entering from the main stem we assume equally divided between the alveolar pores. Mandelbrot introduced the fractal character to porosity because of self-similarity principle [19]. In this context, we may employ recurrent ladder networks of parallel cells-each cell denoting an alveolus. Although Mandelbrot introduced fractality in porous materials as randomly distributed (si ze and location), due to the very compact structure of the lung ductuli we may consider this parallel arrangement of cells as deterministic [17]. Fig 2 depicts the recursive parallel arrangement of series RC cells.
In an idealized context, fundamental law of physics for air mass is given by: with v the air velocity, M the air mass and F the force. Air flow is given by Q(t) = v(t) Á S, with S the alveolar surface. Air pressure is then introduced as P(t) = F(t)/S. Notice that in case of disease, the breathing surface is decreased, and the same amount of pressure will require a higher force-which is indeed the case of increased work of breathing in patient with various degrees of respiratory obstruction. The arrangement from Fig 2 has an admittance of the form: r with their recurrent ratios. We have shown previously that such ladder networks converge to a lumped form [13,30]: with m ¼ r cþr the fractional order term value dependent exclusively on the ratios of the alveolar parameters. This relation defines the admittance of the air-alveolar interface.
Assuming initial conditions zero for both air flow and pressure at the interface (deviation values with respect to operating conditions), then Eq (15) becomes: and in Laplace transform: and which has similar form as the last term in Eq (7). Combining these relations leads to: where n = 1 + m with m (from Eq (17)) the non-integer order differentiator and t ¼ M The roots of the characteristic equation are: for k integer and These roots will be complex conjugated and determine the oscillatory nodes of the relaxation. The time constant depends on the mass M and when it varies, the roots move along a straight line with an angle ϕ determined exclusively by the fractional order n. Since this n does not depend on mass, the damping properties will not vary with the mass, but with the structural properties of the porous material-in our case the alveolar structure. The damping if calculated as: whereas the natural frequency of such a structure is given by In conclusion, we have shown that properties of a single alveolus, as well as those of a network of alveoli may be naturally characterized by fractional order impedance models.

Structural changes in COPD
Chronic obstructive pulmonary disease (COPD) is often paired with emphysema and referred to as small airway disease [15,31]. It features a prolonged time constant for lung deflation, due to increased resistance of the small conducting airways and increased compliance as a result of emphysematous destruction. The latter implies breaking the alveolar walls and disruption of the adjacent connective walls, inflammation and tissue mass enlargement. Density of the lung is decreased, with airspace enlargement and CT scans have shown various degrees of heterogeneity in density of porous tissue with various degrees of obstruction [32].
Using a gross analogy to the Sierpinski triangle of various densities, as illustrated in Fig 3, we can enumerate the following changes in our impedance model as a result of airway remodelling in COPD: • the R and C parameters, along with their corresponding recurrent ratios, will vary in different respiratory zones; • the alveolar changes will affect the fractional order term value n; • within respiratory zones, different such fractional order term values will exist; • the porous structure of the network of alveoli will then contain multi-fractal spacial distribution and thus multi-fractal dynamics.
Fractal and multi-fractal analysis has been given lots of attention in medical applications [10]. Specifically for lung applications, dissimilarity principle has been applied to classify lung parenchyma based on image processing techniques [33], as well as pulmonary emphysema [34]. Airway remodelling has been correlated to fractal dimension in asthma [35], providing insight into space filling mechanism. Pulmonary embolism has been detected using multifractal analysis from lung scans and provided useful classification tool [9].
It is interesting to provide a correlation between changes in COPD at the alveolar level and fractional order impedance values evaluated in various degrees of obstruction. As illustrated in Fig 3, as COPD progresses, changes in the R, C parameters and structure of the alveolar network will affect the impedance and thus the lung function will change.
As the ratios of resistance and capacitance increase with various degrees of COPD obstruction (i.e. classified as GOLD2, GOLD3 and GOLD4, from mild to severe degree of obstruction), the fractional order term value n will vary. If the relative increase is the same, the value will remain unchanged. If resistance has absolute variance higher than capacitance, then the values of n will increase. Similarly, if capacitance has absolute variance higher than resistance, then the values of n will decrease.
The implications for damping are complex. Due to relative increase in compliance as a result of more empty spaces in the alveolar structure, damping will apparently decrease, but due to increased tissue density, the overall damping is increasing. Notice that if alveolar surface is decreased with pathology, the time constant defined in previous section for the network will also increase. This is in line with the pathology of chronic obstructive pulmonary disease (COPD) which has long deflation periods. Typically, the natural frequency in COPD is usually increased, implying that as the time constant is increased, the term in n from Eq (25) has to increase faster than the term in τ. In conclusion, local structural changes are more important than the network spacial distribution in terms of impedance.

Patients
The study includes 43 COPD diagnosed patients (! 60 years) who came for periodic evaluation of their lung function at Ghent University Hospital, Belgium. The study group included: 18 subjects diagnosed with COPD-GOLD II; 15 subjects diagnosed with COPD-GOLD III and 10 subjects diagnosed with COPD-GOLD IV. Biometric and spirometric variables are listed in Table 1. Written informed consent was obtained from all participants. This study and the

Lung function device
The device presented in this paper and used for measurements is a third generation of the prototype described in [36][37][38]. The system consist in a group of fans located on each extreme of the principal pipe (see Fig 4).
One group pushes the air into the tube, while the group on the other side of the tube, extracts the air. The controlled difference in speed between the two groups will generate a controllable pressure inside the pipe. The pipe has 2 inches diameter and has been filled with tubes with smaller diameter (straws) in order to preserve laminar flow. Despite the careful design, nonlinear effects are not completely avoided and their contribution may be accounted for. The end part of the device is formed a pneumotachograph and two pressure sensors. A mouthpiece is mounted on the device for every patient tested (single-use). Since the main goal is to design a suitable optimization technique in order to improve the signal-to-noise-ratio in the signal of interest, efforts have been made to develop a device able to perform low frequency measurements. By measuring at frequencies close to the breathing frequency of the patient, important distortions are present in the system coming from the interference of the breathing signal with the desired excitation signal (typically an optimized multisine). One of the approaches was to introduce open loop feedforward compensation (FF) and closed loop feedback compensation (FB). This has been achieved with online estimation of breathing period for FF and Proportional-Integral-Derivative controller design for FB.  Structural changes in the COPD lung and diffusion mechanism 3 Results and discussion These parameters are given below, with units as given above [39][40][41]:

Analysis of the model in simulation
It can be observed that the lumped impedance follows closely the explicit impedance in a limited frequency range. This interval will depend on the parameters Eq (26). Notice that changes with disease at alveolar level will affect the cell elements Eq (10), consequently will change the lumped impedance form and affect the fractional order operator. Using the impedance of an alveolus, we can now construct a network of alveoli. The origin of this impedance is not essential, as long as its form has the term in fractional order to help us understand changes with remodelling. Supposing that the flow is divided equally through the 5 alveolar sacs as in Fig 6, then all impedances are equal and receive the air-flow in parallel.
This allows us to write the following admittance equation: with Yt = 1/Zt the total admittance and Y i the admittances of each alveoli, with the form from Eq (14), with fractional order value n = 0.5. The total admittance will not change in phase, only in amplitude, dependent on the amount of alveoli we take into calculus of Eq (27). In pathology, as explained before, the mechanical properties are altered as a result of remodelling and lead to different fractional order values. Structure of alveolar sacs also modifies, leading to various consistencies of network of alveoli, as conceptually illustrated in Fig 3. So if we consider that we have the same admittance form as in Eq (27), but with the following elements: where Zs has the same form and values as in Eq (14) except the fractional order n. The phase will be altered, varying between −0.3 Ã 90˚and −0.5 Ã 90˚.
In COPD, we also have the breaking of alveolar walls, creating larger spaces. This will induce a different form of the admittance, as: with the same definitions as Eq (14) and n = 0.5 for i = 1, 2, 3 and n = 0.3 for i = 4, 5. This will also affect the phase of the total admittance, but the limit values will be the same as those from Eq (28). This mathematical framework provides the necessary support for the multi-fractal dynamics exhibited by the lung parenchyma and much discussed in literature under the resulting property of self-organized critically system Phase dynamics and thus also phase transitions are direct result of remodelling effects as a self-defence mechanism of the respiratory system towards disease. The result for the impedance obtained with Eq (27) is depicted in Fig 7 and from Eq (28) is depicted in Fig 8. It can be observed that the limit phase values are close to the expected values related from the fractional orders (−30˚, −45˚).
Finally, the result of the impedance with broken alveolar walls, as defined in Eq (29) is given in Fig 9. As a general remark about all these simulation results, the magnitude values are unrealistically high. In reality, the very high alveolar surface (2 24 alveolar sacs according to [17]) will The simulation studies performed here suggest that the convergence of the fractional order in model from Eq (7) has a structural origin and variations are directly related to structural changes. To our knowledge, such an analysis has not yet been reported in specialised literature.
Changes in gas exchange area paired with disease evolution are related to the mass of air which can be involved in the ventilation process. This implies that the term in Eq (21) has a  Structural changes in the COPD lung and diffusion mechanism time constant which is varying with the relative ratio between gas exchange surface and air mass. The damping factor of such a system is independent on this term, but the natural frequency is inversely related to it via relation Eq (25). The natural frequency will increase for same surface with lower air mass, as expected in COPD patients. This is then always reflected by the balance of inertial and capacitive properties in these patients which tend to be equal at frequencies around 15 Hz and higher (in healthy, the resonant frequency is around 8Hz). In terms of respiratory impedance, the resonant frequency is where the negative part of the imaginary part of impedance (reactance) is equal to the positive part (i.e. crosses the zero line) [13,16]. The damping factor is a parameter reflecting the capacity of the material for energy absorption. In the case of lung tissue damping is mostly characterize by viscoelasticity. The increased lung elastance in COPD due to more empty spaces in the alveolar structure will result in higher values of tissue damping than in healthy patients. The model from Eq (21) has been used to mimic healthy and COPD patient by changing morphological parameter values in Eqs (12) and (13). This impedance was then fitted to the model in Eq (7) to extract parameters from Eq (8). The results are given in Fig 10 for tissue damping (G r ).
It can be observed that as damping is increased, in COPD it relaxes as the alveolar walls are broken with COPD progress.

Experimental data
We do not have the possibility to measure non-invasively the impedance of the alveolar network. Instead, we employ the forced oscillation technique, a well-known lung function test broadly used to assess mechanical properties in lungs [16].
Figs 11 and 12 depict the boxplots for the tissue damping G r and tissue hysteresivity η r . Statistical analysis has been performed to identify if a significant difference between the group exists. Anova tests have been applied to the data and it has been observed a significant difference in tissue damping between the three groups: i.e. between: COPD-GOLD II and COPD-GOLD III; COPD-GOLD II and COPD-GOLD IV; COPD-GOLD III and COPD-GOLD IV.  The p values identified for this case is lower than 0.01. Same analysis has been performed also for tissue hysteresivity and a p value lower than 0.05 has been identified.
The damping factor is a material parameter reflecting the capacity for energy absorption. In materials similar to polymers, as lung tissue properties are very much alike polymers, damping is mostly caused by viscoelasticity, i.e. the strain response lagging behind the applied stresses [4,5]. The loss of lung parenchyma (empty spaced lung), consisting of collagen and elastin, both of which are responsible for characterizing lung elasticity, is the leading cause of increased elastance in COPD. Damping will increase due to increased tissue density in places where alveolar walls are not broken.
Since pathology of COPD involves significant variations between inspiratory and expiratory air-flow, an increase in the hysteresivity coefficient η r reflects increased inhomogeneities and structural changes in the lungs. In emphysematous lung, the caliber of small airways changes less than in the normal lung (defining compliant properties) and peripheral airway resistance may increase with increasing lung volume. Based on our proposed derivations to obtain Eqs (5) and (21) the model from Eq (7) has a strong theoretical basis linked to i) mechanical (resistance, elastance) properties of lung parenchyma, ii) structural layout (recurrent properties) and iii) tissue density (porous character). Alterations in these properties have an effect of the overall impedance parameters and these can now be linked to changes in small airways and alveolar levels. This in turn affects the degree of non-linearity in the tissue dynamics. The non-linear index T has been calculated for each group and the results obtained are presented in Fig 13. It can be noticed that there is a significant difference between the groups. From the results obtained we may conclude that the non-linear index T increases as COPD disease advances which is in line with prior rationale. Statistical analysis for the the non-linear T-index in the measured group has been performed and the results obtained are reported in Table 2.
In Fig 14 the evolution of the hysteresivity factor as a function of the non-linear index is presented.
The results obtained indicate that as COPD continues to deteriorate the alveolar walls the heterogeneity of the lungs will decrease as larger spaces occurs. However, the effects from changes at macromolecular level contribute to increased non-linear dynamics. The results indicate there exist a correlation between the non-linear index T and the heterogeneity factor (i.e. T increases ad the COPD disease advances). The added value of the T index is that it can be extracted directly from the measured signals without the necessity of Eq (7).

Conclusions
This paper presented the available tools emerging from fractional calculus to model changes in chronic obstructive pulmonary disease (COPD) and to understand their effect on model parameters values. The theoretical basis helps in understanding the interplay between all these structural and dynamical changes. The results obtained indicate that there is a direct relation between the non-linear index T and the heterogeneity in COPD lungs. A limitation of this study is that effects of decrease in vascularisation associated with broken alveolar walls are not investigated. This implies an augmentation of the model with complex mathematical  λ-ratio of resistances/total levels χ-ratio of compliances/total levels s-Laplace operator K(λ, X)-gain factor n-fractional order operator Z N (s)-impedance

Author Contributions
Conceptualization: DC RDK ED CI.

Data curation: DC RDK ED CI.
Formal analysis: DC RDK ED CI.
Funding acquisition: DC RDK ED CI.

Methodology: DC RDK ED CI.
Project administration: DC RDK ED CI.

Resources: DC RDK ED CI.
Software: DC RDK ED CI.