Fractal Physiology, Breath-to-Breath Variability and Respiratory Diseases: An Introduction to Complex Systems Theory Application in Pulmonary and Critical Care Medicine

1. Introduction

Physiologic data measured at the bedside often display fluctuations at scales spanning several orders of magnitude. These fluctuations are extremely inhomogeneous and appear irregular and complex whereas in the medical literature, they are often regarded as noise and are neglected. However, they may carry information about the underlying structure or function of the heart and lungs. Examples include fluctuations in heart rate, respiratory rate, lung volume and blood flow [1]. The central task of statistical physics is to study macroscopic phenomena that result from continuous microscopic interactions among many different components. Particularly, physiologic systems such as the cardiovascular and respiratory systems, are good candidates for such an approach, since they include multiple components and are affected by varying neuro-autonomic inputs, continuously over time [1].

Healthy state exhibits some degree of stochastic variability in physiologic variables, such as heart and respiratory rate. This variability is a measure of complexity that accompanies healthy systems and is responsible, according to Buchman, for their greater adaptability and functionality related to pathologic systems [2]. Loss of this variability has been shown to precede the onset of sepsis and multiple organ dysfunction syndrome (MODS) [3-6]. Studying physiological signals of critically ill patients, such as heart and respiratory rate can easily identify ‘hidden’ information concerning inherent dynamics and overall variability within time series [4]. Recognition that physiologic time series contain such information, related to an extraordinary complexity that characterizes physiologic systems, defies traditional mechanistic approaches based on conventional biostatistical methodologies and has fueled growing interest in applying techniques from statistical physics, for the study of living organisms [6]. Through those techniques different ‘physiomarkers’ reflecting variability of various biosignals (e.g., heart rate variability that is the variability of R-R interval in the electrocardiogram) can be estimated. These indices of healthy complexity seem to fulfill the requirements of contemporary critical care medicine for better and more accurate early warning signs, since they are based on high-frequency measurements (sampling rate at least 250 Hz). Different monitors sample original physiological signals at discrete sample intervals and the rate of sampling determines how well the signal is reconstructed. In this respect, and based on the Shannon-Nyquist theorem, accurate reproduction needs a sampling frequency at least two times the highest frequency component of a signal’s frequency spectrum, otherwise the signal is undersampled. On the contrary, conventional biomarkers, such as different cytokines, are typically measured once per day, exhibit marked pleiotropy and poorly reflect inherent dynamics of the system under study, leading finally to loss of information regarding real time changes in patient’s physiology. The combination of structural indices such as the left ventricular ejection fraction (LVEF) with autonomic function indices derived from heart rate variability analysis (HRV) has been recently proposed as the state-of-the-art method for risk assessment among patients with acute myocardial infarction or severe congestive heart failure [7].

In this respect, a few studies have explored indices derived from breathing pattern variability analysis in patients with pulmonary diseases or in critically ill patients during their stay in the Intensive Care Unit (ICU), for assessing readiness for liberation from mechanical ventilation, respectively. Reliable assessment of breathing variability involves a set of signal processing techniques that can be applied to various respiratory signals and can also extract different sets of information concerning intrinsic breathing dynamics.

2. Signal processing techniques

2.1. Linear methods

A considerable body of data suggests that healthy individuals exhibit breath-to-breath variability of breath components in a breath series [8,9]. Breath-to-breath variations have been traditionally treated as random uncorrelated white noise superimposed on the output of the respiratory controller [9,10]. According to Tobin, the random fraction aids respiratory system to perform tasks other than gas exchange, such as speaking [8,11]. Only simple statistics such as mean, variance and coefficient of variation (CV=standard deviation/mean) can estimate random variational fraction after averaging over many breathing cycles. Since variability in complex living systems is not only an artefact of biological noise but also an intrinsic property of various control mechanisms, different types of deterministic (non-random) variability have been described in the pattern of breathing. These types include correlated and oscillatory fractions [12-14].

Autocorrelation analysis calculates coefficients that quantify the fraction of variational activity that is correlated on a breath-to-breath basis. Observations through time may be correlated with a lagged version of themselves. A set of time series values is taken as the first set and the same set is taken as the second, except lagged. If the autocorrelation coefficient retreats from 1.00 as the lag increases and then returns to nearly 1.00, the studied time series behaves in a periodic way. A plot of autocorrelation coefficients on the vertical axis with different lag on the horizontal axis is termed a correlogram; whereas time values at which autocorrelation coefficients reproaches 1.00 indicate periodicities within the time series [13].

Periodic signals can be decomposed into a frequency spectrum of oscillating signals. The different frequency components can be estimated through the Fast Fourier Transformation (FFT) of a time series. The method is called power spectrum density (PSD) and displays in a plot the relative contribution (amplitude) of each frequency [15,16], whereas the area under the power spectral curve in a particular frequency band is considered to be a measure of variability (power) at that frequency. PSD of breathing signals can determine which fraction of variational activity is oscillatory at a particular frequency on a breath-to-breath basis. It has been proposed that the standard deviation (SD) for each breath component can be considered as a measure of gross breath-to-breath variability [13,17].

2.2. Non-linear methods

The above methods have been described as linear, easy to interpret and accessible. However, their application supposes stationary time series behaviour, meaning stability of statistical properties of signals along time. In these cases, any variation in measurements is considered to be random sampling error around a ‘true’ mean [18]. Furthermore, they present insensitivity to the orderliness of data and lack the ability of describing systems’ inherent dynamics. For instance, a time series can be variable but not complex. Conversely, a time series can be less variable but highly complex, therefore variability and complexity that better describes inherent dynamics of non-stationary signals, are two different and independent aspects of a time series [18]. For the above reasons, different nonlinear complexity assessment techniques have been studied as weaning descriptors in few human studies and for the estimation of breathing complexity in different experimental models. The methods mostly used include approximate entropy (ApEn), sample entropy (SampEn) and visual assessment through Poincaré plots [19-23].

ApEn was introduced by Pincus [19,20] as a quantification of regularity in data and compares each group of consecutive measurements over a predefined time window to every other group of measurements of the same time length. ApEn is a measure of the likelihood that patterns recur over specified time intervals. Regular signals are expected to have low ApEn, while complex ones take on higher ApEn values. Due to ApEn’s dependence on the record length an alternative statistic named sample entropy (SampEn) was introduced by Richmann and Moorman [21] with the benefit of reduced computational load.

Visual assessment methods include the application of Poincaré plots. For this plot analysis, each value of the original time series [e.g., respiratory rate (RR)] R-R_n is plotted against the value of the immediately following R-R_n+1 for a predetermined segment. The plot can be quantified by the two values of standard deviations, SD1 and SD2 as indicators of the dispersion of RR points [22,23]. SD2 is defined as the dispersion of points along the line-of-identity whereas SD1 is defined as the dispersion of points perpendicular to the line-of-identity through the centroid of the plot. Table 1 describes the mostly used methods for quantification of breathing variability and complexity.

Fractal analysis constitutes a subset of non-linear methods and will be discussed separately.

Method	Type	Characteristics	Limitations
Fast Fourier Transformation (FFT) [15,16]; Approximate entropy (ApEn) [19,20]; Sample entropy (SampEn); [21]; Poincaré plots for visual assessment [22,23]	Linear; Nonlinear; Nonlinear; Nonlinear	Suitable for analysis of stationary signals of long duration; Estimates regularity within non stationary signals. The higher its values the more easily predicted behavior of the signal and vice versa; Guidelines for parameter selection are provided. Low computational load. Better consistency than ApEn; Easy and fast estimation of signal dynamics (periodic, quasiperiodic, chaotic pattern)	Reduced capability of detecting adjacent peaks in signals with short duration; There are no guidelines for parameter selection and dependence on N (length of data) when short signals are processed; Reduced accuracy for short time series and wide confidence intervals for highly correlated signals; No automation in feature extraction

Table 1.

Summary of the most frequently used linear and nonlinear methods for estimation of breath-to-breath variability /complexity.

2.3. Fractals and power law: Basic concepts

Fluctuations of a variable can be characterised by its probability density distribution. A way of estimating its characteristics is the construction of a histogram after normalisation, so that the area under it will be equal to one. Often, this distribution N(x) of a variable x follows the so called power law form: N(x) = x^-d meaning that the relative frequency of a value x is proportional to x raised to the power of –d. If we plot the logarithms of this relationship we have a linear equation:

$$\left[N\left(x\right)=x^{-d}then,log\left(N\right)=-d*log\left(x\right)\right],$$E1

whereas d is the negative slope of a straight line fit to N. This slope is frequently called β slope or exponent [24].

Power law distribution behaves differently than Gaussian distributions. Its tails are very long (long-tail distribution), representing the relative frequency of occurrence of large events. This means that the probability of large or rare events is much higher compared with a Gaussian. Power laws describe dynamics that have a similar pattern of change at different scales and they are called ‘scale invariant’. On the contrary, Gaussians are characterised by typical values, such as those corresponding to their peaks [25]. Moreover, the power law describes a time series with many small variations and fewer and fewer large ones, whereas the pattern of variation is statistically similar regardless of its size. Magnifying or shrinking the scale of the signal reveals the same relationship, a property that has been called ‘self-similarity’ and is a fundamental characteristic of fractals [24,26]. Fractals are self similar objects because small parts of the structure at increasing magnification appear similar to the entire object. Akin to a coastline, fractals represent structures that have no fixed length, since it increases with increased magnification of measurement. This is why all fractals have noninteger dimensions, the so called fractal dimensions (FDs) [24,27].

The concept of fractals can be applied not only to structures that lack a characteristic length scale, but also to signals that lack a characteristic time scale. In this case, the relationship between the statistical properties of the fluctuations of the signal and the time window of observation follows the power law. The meaning of such behaviour is that future values in a time series are dependent on the past, displaying correlations over time, whereas the system that produces the signal exhibits a kind of memory [24,28].

In order to evaluate the power law of a signal it is necessary to compute the power spectrum. For that reason, a Fourier transformation is applied to the signal in order to decompose it to different frequency components that are included within the time series. Every time series can be considered as a sum of sinusoid oscillations with different frequencies. The fast Fourier transformation (FFT) that is a method for fast estimate of Fourier transform, transforms a signal to a sum of cosine and sine oscillations whose amplitudes determine their contribution to the whole signal. This frequency domain analysis displays the contribution of each sine wave as a function of its frequency, whereas its square is the power of that frequency in the whole spectrum of the signal [16]. The increased variability/complexity is a hallmark of health, whereas many large clinical studies in cardiovascular medicine have proven that loss of variability is associated with sudden cardiac death, post-myocardial infarction (MI) heart failure and ventricular fibrillation [29].

In the case of power law calculation, the plot of the log-log representation of the power spectrum (log power versus log frequency) gives rise to a straight line with a slope of approximately -1. As the frequency increases the size of variation drops by the same factor (scale invariance). The values of the β slope/exponent can reflect the inherent dynamics of a system. Values near 1 are supposed to reflect fractal-like behavior, whereas values lower than 0.5 represent a system without any correlations, lack of memory and finally chaotic-like and unpredictable evolution in time (white noise). On the contrary, values of β slope higher than 1 or even near 1.5 characterize strong correlations within the signal and a highly predictable and almost periodic evolution in time (brown noise) [26,27]. Goldbereger [30] has studied cardiovascular dynamics in health and disease and has found that both unpredictable (random-walk) and periodic behaviors represent loss of physiologic function and correlate with lack of fractal properties of heart rate signals in patients with cardiovascular diseases. Similar results have been found also in critically ill patients with severe sepsis and septic shock [26].

3. Breathing pattern dynamics: physiological and pathophysiological implications

Neurons in the brainstem govern respiratory rhythm through a network of coupled oscillators. Critical components of this network are located in a specialised region of the brainstem called the pre-Botzinger complex (pre-BotC). This complex system can generate a wide variety of rhythms [31]. Focal activation of this region has been shown to increase the frequency of inspiratory motor bursts. Del Negro and colleagues [32] showed that progressively elevating neuronal excitability of the pre-BotC of neonatal rats in vitro by increasing artificial cerebrospinal fluid K⁺ concentration causes periodic modulation of the inspiratory rhythm, reflected by cranial nerve XII motor discharge, in a well defined sequence of behavioural states, characterised by periodic oscillations, quasiperiodicity and ultimately disorganised aperiodic activity. In another experimental study with anesthetised adult cat models, Chen et al. [33] found that both focal hypoxia and chemical stimulation of pre-BotC can produce a marked excitation of phasic phrenic nerve discharge, characterized by high amplitude, rapid rate of rise, short duration bursts and reduced complexity, estimated with approximate entropy (low ApEn values).

The above studies support the hypothesis that central respiratory centers are responsible for different breathing patterns with various degrees of variability and complexity in different settings and levels of stimulation. In addition, they can also adapt ventilation to metabolic needs through integration of afferent information, since the respiratory pattern is created by integration of different inputs from chemoreceptors, chest wall and pulmonary receptors, the cerebrum, vagal afferents and non-respiratory central mechanisms [34,35].

Both hypercapnia and hypocapnia may alter ventilatory pattern. Jubran et al. [36] observed that hypercapnia in humans increases the gross variability of minute ventilation (MV) and V_T but decreases that of inspiratory (Ti) and expiratory (Te) times. In addition, hypercapnia makes ventilation more monotonous with increased autocorrelation between adjacent values. Fiamma et al. [37] studied breath-by-breath variability and complexity of different ventilatory signals (instantaneous ventilation, tidal volume, mean inspiratory flow, inspiratory and expiratory times) in eight healthy subjects with different CO₂ levels. Hypercapnia reduced breathing variability and increased all complexity indices whereas hypocapnia induced reciprocal changes. The authors concluded that chemoreceptors may exert a strong and inverse influence on ventilatory variability and complexity.

Cortical and subcortical effects upon breathing patterns dynamics are also influenced by slow-wave sleep that seems to reduce respiratory complexity [38] whereas panic-anxiety disorders may increase it [39]. Furthermore, Samon and Bruce reported that breathing complexity decreases with anesthesia and vagotomy [40].

Apart from chemoreceptor signalling, chest wall and pulmonary receptors may continuously affect central neural output, especially during resistive breathing. Brack and Tobin [11] measured breathing variability over one hour in ten patients with restrictive lung disease and in seven healthy subjects. They found that variability of Ti, Te and V_T, were significantly reduced in the patients group compared with the healthy group. Furthermore, autocorrelation coefficients were increased almost 3-fold in the patients group, indicating increased periodicity. According to the authors, the decreased breath-to-breath variability in restrictive lung disease patients is a compromise between increased effort and carbon dioxide clearance and arises from their voluntary control of ventilation, as they ‘choose’ to do it in order to reduce respiratory distress [11,41,42].

Of particular importance in the ICU setting is the potential impact of systemic inflammation on breath-to-breath dynamics as suggested by endotoxin response studies. In a clinical study of Preas and colleagues [43], 12 healthy subjects were randomized to receive endotoxin or saline. Administration of endotoxin after 3 to 4 hours increased RR, decreased Ti, produced dyspnea, augmented autocorrelation coefficients within RR time series and decreased random fraction of variational activity of frequency. These changes were related to changes in arterial carbon dioxide tension. The authors concluded that endotoxin has a direct effect on respiratory controller function whose increased output causes dyspnea. They suggested that decrease in random fraction of breath variability, meaning reduced freedom to vary the respiratory cycle, was attributed to a decrease in circulation time between the lung and the chemoreceptors, secondary to an increase in cardiac output. Since ibuprofen, a cyclooxygenase inhibitor, did not abolish dyspnea, something that seems to happen in healthy exercising subjects, the authors proposed that endotoxin augments respiratory center output through other alternative pathways [44, 45].

4. Fractals and power law in pulmonary physiology

Many organs in different biological systems have fractal structure. Fractal branching reduces the distances over which materials are transported, providing rapid and efficient delivery of nutrients [46]. The lung offers many examples of self-similarity properties. Weibel and Gomez [47] first measured the morphology of human airways and found an exponential relationship between the diameter and the generation number of the conducting airways. Mandelbrot [48], who was the first who introduced the term fractals, discovered a unifying scaling pattern of the branching in the lung. Its higher fractal dimension corresponds to a more complex branching, whereas a lower one reflects a more homogeneous structure. Moreover, regional pulmonary blood flow has been shown by Glenny to exhibit spatial and temporal fractal patterns [49]. The structure of alveolar surface has been also found to be well described by power laws, reflecting scale invariance [50]. The probability distribution of airway opening during inspiration behaves also according to the power law [51].

Another property of fractals and power laws in pulmonary physiology is error tolerance during development. In simulations of airway morphogenesis during lung development, West [52] compared a power law branching rule with an exponential decaying one and found that in the first case, the system was less susceptible to errors introduced into the branching process. These same properties suggest that living systems are capable to operate similarly at different scales, meaning that whenever environmental conditions change they can adapt more easily to their surroundings.

Aging has been proven by Lipsitz and Goldbereger [30] to be significantly associated with loss of complexity of physiological signals, leading to decreased ability to adapt to different physiological insults. Using different algorithms for estimating fractal properties and power law behavior, these authors found that the β slope of different signals in elderly was either reduced (decreased lower than 1) or augmented (increased higher than 1) compared to younger adults, indicating chaotic or periodic behavior, respectively. Peng and co-workers [53] showed that aging was associated with a breakdown of fractal dynamics of respiratory signals via a decrease in β slope towards 0.5 (randomness). Concerning early stages of development in humans, one study [54] has found that ultrasonographic patterns for assessment lung maturity showed fractal properties with a power law behavior. In addition, the β slope increased with gestational age from 28 to 38 weeks. Szeto and co-workers [55], calculated β slopes of different respiratory signals in human fetus and showed its movement from randomness towards fractal behavior with gestational age. In conclusion, it seems that there is great variability in complexity with age in early life, after which complexity decreases with aging.

5. Fractal properties of the lung in disease: data from clinical and experimental studies

Alterations to fractal properties are related to different pathologies and could have clinical implications for diagnosis and treatment. Physiologic time series, such as heart and respiratory signals, show similar alterations in their power law behavior in different disease states. Mackey and Glass [56] have introduced the term ‘dynamic diseases’ to describe states with loss of fractal properties of organs and power law dynamics of signals that are produced form the above structures. For example, loss of heart rate variability that is the variability of the R-R in the electrocardiogram, has been found in patients with heart failure [57], atrial fibrillation, septic shock and multiple organ dysfunction syndrome [6,26]. In respiratory disorders, a classical example is the highly periodic variation in respiratory frequency, seen in Cheyne-Stokes respiration. Penzel [58] has observed loss of fractal properties of heart rate signals during episodes of obstructive sleep apnea.

Macklem [59] was the first who raised the question of whether airway function can be studied using tools from chaos theory and the paradigm of complex systems. Que and co-workers [60] studied the distribution of forced oscillatory resistance in asthmatics and demonstrated that lung function exhibit loss of fractal properties during severe asthma. Frey and colleagues [61] applied fractal methods to twice-daily peak expiratory flow (PEF) in asthmatic patients and showed that the β slope was reduced, whereas it become more regular with standard long-acting β2-agonist treatment and more random with short-acting β2-agonist treatment, respectively. Moreover, the authors were able to demonstrate that the higher the β exponent when a patient was not under any treatment, the larger the improvement of his condition upon administration of long-acting β2-agonist therapy.

In another study, Suki and co-workers [51] studied the dynamics of airway opening and crackles, using a simple mathematic model of the periphery of airway tree. Suki found that the time series of crackles emitted during airway opening follows a power law distribution. Additionally, as the crackles propagate up the tree, the sound amplitude is attenuated at successive bifurcations, whereas its distribution follows the power law. The same has been found for the time intervals of the ‘jumps’ by which airway resistance decreases upon lung inflation by a constant flow. In a study of Boser and colleagues [62], the fractal dimension of airways was computed using autopsy material from three groups: fatal asthma, non-fatal asthma and non-asthma controls. The authors were able to show that the average FDs of both fatal (1.72) and non-fatal asthma groups (1.76) were significantly lower than that of the third control group (1.83, p<0.05), whereas the lower fractal dimension correlated with a decreased overall structural complexity and pathologic severity of disease.

Venegas and colleagues [63], using positron emission tomography (PET) imaging and computer modeling showed that in cases of bronchoconstriction and when smooth muscle activation reaches a critical level, localized clusters of poorly ventilated lung regions can develop abruptly in discrete steps. These steps are called avalanches and can lead to new stable conditions. Because of the fractal structure of the airways, small initial heterogeneities that are always present and particularly in the diseased lung, can be amplified, leading to sudden patches of poorly ventilated lung regions. Another implication is that since airways are organized into a fractal network embedded in the elastic parenchyma, the constriction of one airway can propagate and cause an avalanche-like constriction in large parts of the lung. The same holds true for the opposite process, where opening of airways during inhalation takes place in discrete steps [63-65].

Suki [51] has also demonstrated that airway opening upon inflation occur in avalanches with power law distribution of both the size and time intervals between them. The significance of these results is that the probability of finding a large avalanche is much higher than it would be if the distribution were Gaussian or exponential, so both the magnitude and timing of pressure excursions applied at the airways (i.e., using mechanical ventilation) may be critical in triggering the avalanche process of alveolar recruitment [24].

In conclusion, these studies in asthma show that when the airways are likely to approach their critical closing threshold pressure, a small stimulus can provoke a catastrophic cascade of airway closure and for that reason, there is such poor correlation between the trigger and the outcome in asthmatic patients. Moreover, the history of symptom fluctuations seems related to the structural changes of the airway tree (power law distribution of airway diameter) [65].

Airway recruitment may affect alveolar recruitment as well. Sujeer and co-workers [66] have found in mathematical models that the recruited volumes upon inflation with constant flow are distributed according to a power law with a β slope equal to 2. From the above findings it can be supposed that since alveolar recruitment is influenced by airway structure, then the pressure-volume curve may carry information about the airway tree [24]. Whether such models have any value in acute lung injury (ALI) is unclear. In this syndrome, it has been found that the opening pressure distribution does not seem to be always Gaussian, something that is assumed to be the case in the avalanche model [67]. More studies are needed to investigate the pattern of recruitment in ALI, particularly in the case of gravity effect upon ventilation-perfusion mismatch at the level of alveoli [68].

The application of fractal analysis has also shed light into the morphology of the lung in cases of emphysema. Computed tomography (CT) is a sensitive method for assessing lung structure in different pathologies. In general, low attenuation area (LAA) clusters are depicted in pixels with density less than 950 Hounsfield units. These areas incorporate mostly air and assigned a value of 1, whereas pixels with a density higher than 950 include tissue with a value 0. Summing the number of pixels in a cluster gives the cluster size. In that way, a binary map of the lung can be constructed and a few studies have shown that in normal conditions, this map is highly heterogeneous [24]. Mishima and colleagues [69] have found that the probability distribution of LAA clusters follows a power law for both normal subjects and patients with chronic obstructive pulmonary disease (COPD). However, patients exhibited significantly smaller β slopes or exponents, which did not correlate with pulmonary function tests except for diffusion capacity of the lung. The authors suggested that the neighboring smaller LAA clusters tend to coalesce and form larger clusters as the weak elastic fibers separating them break under tension. This process does not change the % LAA but decreases the number of small clusters in favor of larger ones, which result in a reduction of the β slope. Another assumption derived from this study is that the likelihood of finding large LAA clusters is much higher in COPD patients than in normal controls [65].

Another possible application of fractals in pulmonary and critical care medicine seems to include the mechanical ventilation of critically ill patients. In an oleic acid injury animal model, Mutch and colleagues [70] introduced fluctuations according to an algorithm, to mechanical ventilation (biological variable tidal volume and respiratory frequency pro-portional to pre-defined minute ventilation values). Compared with conventional ventilation (with similar minute ventilation), this approach increased respiratory arrhythmia and oxygenation and decreased dead space. According to Suki [24,71], when fluctuations in the form of symmetrically distributed random noise is added to peak airway pressures (noisy ventilation), the mean does not change but isolated values can be augmented, leading to significant alveolar recruitment. In a mathematical model, the authors found that the recruited lung can be 200% larger in the case of biological variable ventilation than during conventional ventilation. Moreover, the standard deviation (SD) of the noise can be manipulated in order to achieve better oxygenation (system’s output), a phenomenon called ‘stochastic resonance’, which has already been confirmed in animal models of ALI [71].

6. Breathing variability and complexity indices as weaning predictors in mechanically ventilated patients

Engoren and colleagues [72] studied 10 control patients who had undergone cardiac surgery within an interval of 12 hours prior to this experiment and 21 patients who required prolonged (> 7 days) ventilatory support. The control group was studied during a weaning trial of 5 cm H₂O continuous positive airway pressure (CPAP). The patient group was studied during 60-to-120 min trials of spontaneous ventilation with 5 cm H₂O positive end-expiratory pressure (PEEP) with a constant (12.2 +/- 6.6 cm H₂O) level of pressure support (PS) for each trial. These patients passed 59 and failed 14 weaning trials. During spontaneous ventilation each breath’s instantaneous respiratory rate and tidal volume were recorded and their ApEn values were calculated for the terminal 1000 breaths, in a series of 100, 300 and 1000 breaths. Receiver operating characteristic (ROC) curves identified cut-point values of continuous variables that predicted a failed weaning outcome. While mean V_T did not vary between groups, mean RR and frequency-tidal volume ratio increased progressively from the control group to both successful weaning (SW) and failure weaning (FW) groups. Conversely, ApEn of RR did not vary between every pair groups, but entropy of V_T increased significantly from the control to FW group and from SW to group FW, with no difference between control and SW groups. The ROC curves for ApEn-V_T and frequency-tidal volume ratio showed similar sensitivity and specificity for predicting weaning failure whereas ApEn-V_T proved to be successful in separating SW from FW, with similar area under the curve (AUC) values of 0.74, 0.75 and 0.73 for 100, 300 and 1000 breaths respectively. Since 100 breaths during a weaning trial may take only 3 to 4 minutes compared with 30 to 40 minutes required for a 1000 breaths trial, ApEn-V_T assessment can be of significant value, permitting a faster determination of weaning tolerance.

According to the authors [72], increased irregularity of the FW group indicated enhanced external inputs to the respiratory controller whereas increased regularity of the SW group suggested greater component autonomy. The first case implies that enhanced external inputs from peripheral receptors that measure pH, Pco₂, Po₂ or from cortical areas that represent dyspnea or anxiety are responsible for increasing breathing complexity [19,20]. Conversely, many studies estimating heart rate variability have found inverse relations between mortality after myocardial infarction or endotoxemia and HRV [73,74]. It seems that patient’s voluntary control over breathing (which does not exist for heart rate) is responsible for different effects of cardiac and pulmonary disease on regularity.

El-Khatib et al. [75] studied 52 patients diagnosed mostly with lung diseases and under mechanical ventilation in two phases: 1. under synchronized intermittent mandatory ventilation (SIMV) with RR ≤ 4 breaths/min and with no PS for 60 minutes and 2. during a CPAP trial of 5 cm H₂O with no PS ventilation for other 60 minutes. Thirty-nine patients (75%) were successfully extubated and the remaining 13 patients failed weaning trial. In both patient groups, different breathing signals were collected (airway flow, volume, proximal airway pressure) and were divided into three intervals of 300 breaths. During SIMV trials, the spontaneous and mechanical peak flows (PF) versus tidal volume scattergrams were constructed and coefficients of variation for both variables were determined for each 300 breath interval separately and for the whole 900 breaths. During CPAP trials, the Kolmogorov entropy indicating unpredictability [76] within the flow-volume loops was estimated, along with the dimension of the respiratory breathing pattern, defined as the number of clusters or clouds of trajectories in the flow-volume loops space, for each data interval. The authors found that the CV of tidal volume and PF, the Kolmogorov entropy and the fractal dimension of the spontaneous breathing pattern were all significantly smaller in the successfully weaning group compared with the failure weaning group.

These results suggested that breath-to-breath variability during SIMV trials and complexity during CPAP trials were increased in patients who failed weaning trials. Patients from the two groups did not differ in terms of age, previous days on mechanical ventilation and blood gases at study entry, whereas the distribution of diagnosis was also similar between the two groups. Since the aim of this study was not to determine a cut-off threshold in terms of breathing pattern complexity for weaning prediction, no sensitivity or specificity were measured. Furthermore, the authors did not present any possible explanation of their results in terms of different pathophysiological mechanisms that may drive various levels of variability/ complexity.

In another study, Bien and colleagues [77] investigated breathing pattern variability in 78 mechanically ventilated patients with systemic inflammatory response syndrome (SIRS) who had undergone abdominal surgery. The patients were divided in the successfully weaning group (n=57) and the failure weaning group (n=21). Within 1 hour before the measurement of breathing variability and under ventilation with PS with a pressure level between 10 and 20 cm H₂O, tidal volume, total breath duration (Ttot), Ti, Te and peak inspiratory flow were calculated. After obtaining these data, the ventilatory mode was switched to PS of 5 cm H₂O plus 5 cm H₂O PEEP for 30 minutes. The coefficients of variation and the SD1 and SD2 from the constructed Poincaré plots for the above variables were determined for at least 300 breaths. Moreover, AUC were computed for ROC curves and compared with airway occlusion pressure, its ratio to the maximal inspiratory pressure (P0.1/Pimax), RR, RR/V_T, and the product P0.1*( RR/V_T). The authors found that average values and the CVs of the 5 measured parameters were significantly lower in the failure weaning group than in the successful one. In addition, Poincaré plot analysis showed that both SD1 and SD2 of the 5 breathing parameters were significantly reduced in the weaning failure group whereas the SD1-SD2 ratio showed no statistical significance. Finally, AUCs of the CVs of the 5 measured variables and both SD1 and SD2 did not differ significantly from those of the 5 clinically used weaning predictors. All the AUC values were within the range of 0.73 to 0.80.

This was the first study reporting that WF patients exhibited a decreased breath-to-breath variability and complexity. According to the authors, the different results from the two previous studies may be due to different weaning protocols that were implemented and different patient characteristics. Their patients had a mean age of 68 years, no history of pulmonary disease and a short duration of ventilatory support (3 days) whereas patients in the El-Khatib study [75] had a mean age of 50 years with predominantly underlying lung diseases and a longer duration of mechanical ventilation (11 days). Finally, in the Engoren study [72] cardiac surgery patients were included and indices of cardiac function were not presented, especially when low cardiac output states are associated with a Cheyne-Stokes breathing pattern.

Recently, Wysocki and colleagues [78] examined breathing variability in 46 medical and surgical patients during 60 minutes spontaneous breathing trials (SBT) and during a 5-month period in four French tertiary university hospitals’ ICUs. There were 32 successful and 14 failure weaning cases. SBT was performed with complete disconnection from the ventilator, meaning none form of ventilatory support. The time series analysis was performed by two independent investigators, one of whom was blinded to clinical information. Furthermore, a signal processing technique for checking stationarity of breath time series was performed and nonstationary data were excluded from analysis. The authors reported significantly increased CVs of tidal volume, Ti, Te, Ttot, Ti/Ttot and V_T /Ti in the SW group, whereas the duration of autocorrelation (number of lags) was significantly shorter in the same group compared with patients who failed weaning trials, meaning that low number of breaths were significantly autocorrelated in cases of weaning success.

The different results from Engoren’s and El-Khatib studies [72,75] were attributed to the prevailing conditions during respiratory variability analysis. According to Caminal and Brochard, different levels of pressure support may alter breath-to-breath variability in an inverse way [79,80], something that could have played a role in the two previous studies. Furthermore, it was postulated the low breathing variability could cause the deterioration of respiratory mechanics rather than result from it, by promoting microatelectasis. In conclusion, Wysocki’s patients did not receive any ventilatory support during weaning trial, limiting accuracy of comparisons with populations from the first two studies.

In a recent study [81], we tried to investigate heart rate (HR) and respiratory rate complexity in patients with weaning failure or success, using both linear and nonlinear techniques. Forty-two surgical patients were enrolled in the study. There were 24 who passed and 18 who failed a weaning trial. Signals were analyzed for 10 minutes during two phases: 1. pressure support (PS) ventilation (15-20 cm H₂O) and 2. weaning trials with PS: 5 cm H₂O. Low and high frequency (LF, HF) components of HR signals, HR multiscale entropy (MSE), RR sample entropy, cross-sample entropy as an indication of coupling between cardiorespiratory signals and Poincaré plots were computed in all patients and during the two phases of PS. Multiscale entropy was recently introduced [82] for quantifying heart rate complexity. Briefly, for a given R-R time series, multiple ‘coarse-graining’ time series are constructed by averaging the data points within non-overlapping windows of increasing length τ, where τ represents the scale factor. Subsequently, sample entropy (SampEn) that represents the negative natural logarithm of the conditional probability that two sequences similar for m points remain similar at the next point with a tolerance r, is calculated for each time series and then plotted against the scale factor, giving rise to the MSE curve. The sum of SampEn over all scaling factors represents the MSE of a signal. For our MSE analysis, the parameter r was set at 15% of standard deviation (SD) of the time series, after normalization (SD=1) and the parameter m (embedding dimension), that is the length of sequences to be compared was set at 2 (data length ranging from 100 to 5000 data points) [82]. Furthermore, we extracted two more features from MSE curves after log transformation of SampEn and scale factor, the fast slopes for small time scales, i.e., those defined by SampEn values between scale 1 and 5 and the long slopes for higher scale values (Figures 1 & 2).

We found that weaning failure patients exhibited significantly decreased RR sample entropy, LF and HF components, compared with weaning success subjects (p<0.001). Their changes were opposite between the two phases, except for MSE that increased between and within groups (p<0.001). A new model including rapid shallow breathing index (RSBI), RR and cross sample entropies predicted better weaning outcome compared with RSBI, airway occlusion pressure at 0.1 sec (P_0.1) and RSBI* P_0.1 (conventional model, R² = 0.887 vs 0.463, p<0.001). Areas under the curve were 0.92 vs 0.86, respectively (p<0.005).

Figure 1.

Heart rate multiscale entropy (MSE) curves between phase H (high PSV) and spontaneous breathing trial (SBT) or phase L (low PSV) in a weaning failure patient. SampEn and scale factor have been logarithmically trans-formed. The MSE curve during the SBT has increased SampEn values; however, the slope during small scales (fast slope) is decreased compared with the high PSV phase.

Figure 2.

Heart rate multiscale entropy (MSE) curves between phase H (high PSV) and spontaneous breathing trial (SBT) or phase L (low PSV) in a weaning success patient. The MSE curve during the SBT has increased SampEn values; moreover, the slope during small scales (fast slope) looks similar with the slope during the high PSV phase.

Different profiles of MSE curves on small time scales between the two groups and according to fast slopes may be due to the impact of respiratory rhythm on heart rate dynamics (respiratory sinus arrhythmia). In general, the lower the amplitude of respiratory modulation, the higher the entropy values tend to be [82]. Thus, in both groups but especially in weaning failure patients, sum MSE was significantly increased between different PS levels, whereas reduced fast slopes could reflect different time constants of respiratory effect upon heart rate control mechanisms.

In another similar study [83], we enrolled thirty-two mechanically ventilated patients undertaking a weaning trial. There were 22 who passed and 10 who failed a weaning trial. Tidal volume and mean inspiratory flow were analyzed for 10 minutes during two phases: 1. pressure support (PS) ventilation (15-20 cm H₂O) and 2. weaning trials with PS: 5 cm H₂O. Sample entropy (SampEn), fractal dimension (FD) and largest lyapunov exponents (LLE) of the two respiratory parameters were computed in all patients and during the two phases of PS. FD was computed according to the Higuchi method [84]. Briefly, FD was based on a measure of length L(k) of a time series, computed at different scales, by using a segment of k samples as a unit in each scale. The value of FD was calculated by a least-squares linear best-fitting procedure as the angular coefficient of the linear regression of the log-log graph of the mean of k values Lm(k) for m=1,2,3…k, with k being an interval time. The length Lm(k) originating from time m was calculated as the normalized sum of absolute differences between the values of point pairs that are ‘k samples distant’ and the length of curve of the time interval k, L(k) was calculated as the mean of the k values Lm(k). If the L(k) related to the scale used (k) linearly in a log-log plot with slope FD, then the curve was said to show fractal dimension.

LLE were computed as follows: Complex systems are considered sensitive to initial conditions and exhibit an exponential divergence in the phase space, which describes in a 3-dimensional axis their different states. Estimation of Lyapunov spectrum and largest Lyapunov exponents (LLE) can assess sensitivity to initial conditions. Briefly, if we consider two points in adjacent trajectories-states of the phase space with a distance between them d(0), after time t the average divergence (separation) will be:

$$[d\left(t\right)=d\left(0\right)*e^{LLE*(i}{}^{\Delta t)}]$$E2

(2)

whereas LLE is the largest Lyapunov exponent. In this study, we computed LLE of mean inspiratory flow and tidal volume signals, using the algorithm proposed by Rosenstein [85], which seems to be useful, particularly in small data sets. Values higher than 0 reflect an unstable and unpredictable system, where nearby points will diverge to any arbitrary separation. Increased LLEs reflect increased sensitivity to initial conditions and characterize unpredictable variations, whereas low values indicate regularity [85].

Weaning failure patients exhibited significantly decreased respiratory pattern complexity, reflected in reduced sample entropy and lyapunov exponents of respiratory flow time series, compared to weaning success subjects (p<0.001). In addition, their changes were opposite between the two phases of the weaning trials. A new model including rapid shallow breathing index (RSBI), its product with airway occlusion pressure at 0.1 sec (P_0.1), SampEn and LLE predicted better weaning outcome compared with RSBI, P_0.1 and RSBI* P_0.1 (conventional model, R² = 0.874 vs 0.643, p<0.001). Areas under the curve were 0.916 vs 0.831, respectively (p<0.05).

Vallverdu et al. [86] in another study with similar design compared with our two previous investigations examined heart rate and respiratory pattern complexity in 78 patients, during weaning trials and using information flow analysis, which describes the regularity of signals by estimating the auto- and mutual information functions. The authors were able to find reduced complexity and a more coupled nonlinear oscillator behavior in weaning failure subjects.

These results parallel those from Schmidt and colleagues [87] who reported increased LLE and Kolmogorov-Sinai entropy values of mean inspiratory flow signals in mechanically ventilated patients, after switching the ventilator from the pressure support mode to neurally adjusted ventilatory assist mode (NAVA). According to these authors, successful spontaneous breathing trials unmask underlying variability and complexity of central neural output, since inspiratory pressure inhibits the respiratory drive. This effect is nicely reflected through the increased complexity indices of flow and is responsible for better neuro-mechanical coupling.

Figure 3.

Phase space of minute ventilation (MV) in a weaning success patient. Different values x(i) of minute ventilation were plotted against the following ones in tree-dimensional axes: x, x+t and x+2t, giving rise to the phase space. Red dots represent data during high PS ventilation, whereas blue dots represent data during the performance of a SBT.

Figure 4.

Phase space of minute ventilation in a weaning failure patient. Red dots represent data during high PS ventilation, whereas blue dots represent data during the performance of a SBT. Scattering of data seem highly reduced, compared with findings in Figure 3, whereas ventilation values are positioned in different parts of the space.

Figure 5.

Distribution of LLE between the two groups of patients (S: successful weaning and U: unsuccessful), during high (H) and low (L) pressure support ventilation. It seems that SW group demonstrated an increase in LLE after shifting to spontaneous breathing, contrary to what have occurred to the FW group.

In another study, Mangin and colleagues [88] investigated ventilatory chaotic dynamics in 17 mechanically ventilated patients during switching the ventilator from the assist-control mode to pressure support mode. They were able to show that both fractal dimension and LLE were increased, particularly in 5 patients who were successfully extubated. Furthermore, the authors supposed that increased breathing complexity may also be attributed to higher vagal afferent feedback during unassisted breathing, as has already been shown by Sammon and Bruce [40].

These studies support our findings that transition between mechanical and spontaneous ventilation is associated with increased complexity of respiratory signals in weaning success patients, since duration of ventilation before the SBTs was similar between groups with different weaning outcome. Moreover, in a study of Burykin and Buchman [89] investigating cardiorespiratory dynamics and synchronization during controlled and unassisted breathing in 13 surgical patients, it was demonstrated that mechanical ventilation reduces significantly both heart and respiratory rate complexity whereas spontaneous respiration is more irregular with increased uncoupling of cardiorespiratory rhythms in weaning success patients.

7. Conclusions

According to Macklem [90], there is a continuum of thermodynamic systems that do not follow the 2^nd thermodynamic axiom (increase in entropy), since they exchange energy with their environment (open systems): from near to equilibrium, like crystals to far from equilibrium systems, like the weather. ‘The amount of energy dissipated determines where the system is situated along the continuum. Between the crystals and weather a sudden phase transition occurs over a small range of energy consumption. It is only there where life can flourish’. From the above statement it is concluded that either a decrease in energy consumption, i.e., during myocardial ischemia, or an increase may both contribute to a shift away from stable state. It has been shown by Que [60] that during asthma, the local increase in metabolic rate is associated with increased variability of respiratory impedance to flow. The same author has proposed the term ‘homeokinesis’ instead of homeostasis [60], as a basic property of living systems that ‘describes the ability of an organism to utilize external energy sources to maintain a highly organized internal environment fluctuating within acceptable limits in a far from equilibrium state’.

There is good basic science describing the physiological and fractal variability observed in respiratory dynamics and substantial evidence that conditions prevalent in the ICU impacts breath-to-breath variability. Different clinical studies have assessed the ability of markers of respiratory variability as predictors of successful separation from mechanical ventilation. However, these studies are not comparable since the authors implemented different weaning protocols in heterogeneous groups of patients and examined various indices of both breathing pattern variability (linear) and complexity (nonlinear) that by their nature capture different aspects of signal dynamics. It seems that the study of Wysocki [78] performed better, since it was conducted in 5 medical centres, stationarity tests were done and possible bias in the time series analysis was significantly controlled. Furthermore, lack of ventilatory support during weaning trial could have eliminated any possible impact of positive pressure ventilation upon breathing variability and could be responsible for the different findings between this and the first two studies. Unfortunately, the authors did not perform complexity nonlinear analysis. Since neither ApEn, nor any other complexity index was calculated, we cannot make assumptions regarding inherent breathing dynamics during discontinuation from mechanical ventilation, because variability and complexity may change in opposite directions.

We suggest that new studies are needed, for estimating both linear and non-linear indices of breathing variability and complexity for weaning outcome assessment, in different and homogeneous groups of patients, and comparing their diagnostic accuracy, along with other traditional and already used physiological predictors. Their association, in terms of ‘coupling’ or ‘uncoupling’ with other biosignals such as the electroencephalogram (EEG), with the aid of robust methods from signal processing theory (i.e., cross-ApEn, coherence etc) could uncover possible inter-relations between breathing and central nervous system dynamics during weaning trials.

However, we think that the development of a consensus between clinicians and modelers is urgently needed for the assessment of different methods used for breathing pattern variability prior to the performance of other studies. In the cardiology literature, the Task Force on HRV has standardized the use of a set of signal processing techniques, including target variables, frequency and duration of measurements and signal quality assessment, aiding in the development of better and more accurate diagnostic or prognostic indices in cardiovascular diseases [16]. Furthermore, many important questions remain unanswered such as the value of single or longitudinal variability measurements for clinically useful information or its prognostic capability for risk stratification between patients or within the same patient (trend analysis).

In conclusion, we believe that the most important aspect of research on breathing pattern dynamics remains the development of international databases of various respiratory signals, with free access from different investigators, such as the Web site Physionet (www.physionet.org), which has boosted research on heart rate dynamics in cardio-vascular medicine [91]. We also think that in addition, such efforts must also be undertaken by pulmonary physicians and those who treat patients with severe respiratory diseases. In conclusion, we suggest that implementation of fractal analysis and breathing pattern variability and complexity measurements in respiratory physiology will result in better understanding of different and dynamic changes during various pathologic states, whereas its impact upon prediction of severe complications will add value to such methods. Clinicians must begin to understand basic principles of complex systems theory and support an interdisciplinary approach for understanding pulmonary diseases, in order to treat their patients more effectively.

Abbreviations

ALI: acute lung injury

ApEn: approximate entropy

AUC: area under the curve

COPD: chronic obstructive pulmonary disease

CPAP: continuous positive airway pressure

CV: coefficient of variation

ECG: electrocardiogram

EEG: electroencephalogram

FD: fractal dimension

FFT: fast Fourier transformation

FW: failure weaning

HRV: heart rate variability

ICU: intensive care unit

LAA : low attenuation area

LLE: largest lyapunov exponents

LVEF: left ventricular ejection fraction

MODS: multiple organ dysfunction syndrome

MSE : multiscale entropy

MV: minute ventilation

PEEP: positive end-expiratory pressure

PEF: peak expiratory flow

PET: positron emission tomography

Pimax: maximal inspiratory pressure

P_0.1: airway occlusion pressure at 0.1 sec

Pre-Botz: pre-Botzinger complex

PSD: power spectrum density

PSV: pressure support ventilation

ROC: receiver operating characteristic

RR: respiratory rate

RSBI: rapid shallow breathing index

SampEn: sample entropy

SBT: spontaneous breathing trial

SD: standard deviation

SIMV: synchronized intermittent mandatory ventilation

SIRS: systemic inflammatory response syndrome

SW: successful weaning

Te: expiratory time

Ti: inspiratory time

Ttot: total breath duration

V_T: tidal volume