Introduction

Most studies on physiological control systems have been based on measures of the average output over some period of time, and little attention has been paid to the mechanisms that determine how these control systems regulate their output in order to maintain a stable homeostatic internal environment. External fluctuations acting on the feedback loops within the control system often result in significant variabilities of the output, and it is becoming evident that variabilities carry useful information on the underlying structure and/or functioning of control systems [1].

One of the key life-support control systems of the body is the respiratory system. Many physiological variables associated with breathing such as tidal volume (V_T) or respiratory rate exhibit significant breath-to-breath variabilities [2–5]. While some studies correlated variability and irregularity in breathing pattern with the presence of obstructive lung disease [6] or the maturation in preterm infants [4,7], the origin of the variabilities is not well understood even in normal subjects.

In an effort to account for the correlated variability in V_T observed in babies, Cernelc et al. [4] proposed to add noise to the neural network model of the brain respiratory oscillator put forth by Botros and Bruce [8]. Their modeling results showed that the phrenic output of the model generated a cyclic pattern with breath-to-breath variations that mimicked the correlation properties of V_T [4]. However, the respiratory system is composed of several mechanical structures such as the lung and the chest wall and its abdominal and rib cage compartments. It is not known whether the centrally controlled signal from the phrenic nerve results in different variabilities of these compartments. Additionally, the role of the passive mechanical properties of the lung and chest wall in the breath-to-breath variations of respiratory parameters has not been determined.

In this study, we characterised the variability of resting breathing pattern in healthy adult humans by using a technique called Optoelectronic Plethysmography (OEP)[9–11] that is capable of measuring all ventilatory parameters, including long-term changes in end-expiratory volume (EEV) changes [12] on a breath-by-breath basis both for the total respiratory system and as well as for the different chest wall compartments such as pulmonary rib cage (RC,p), abdominal rib cage (RC,a) and abdomen (AB). This technique provides an ideal opportunity to assess variabilities of the total and compartmental ventilatory parameters as it does not require a mouthpiece and/or noseclips and it does not add dead space or mechanical load to the subjects. We thus quantified the variability and the long-range correlation properties in the breathing pattern and chest wall volumes using the detrended fluctuation analysis (DFA) [13] technique which can detect intrinsic correlation properties of complex time series. Moreover, to interpret our results, we also developed a simple model of the chest wall and the lung driven mechanically by the respiratory muscles which in turn receive input from the respiratory rhythm generator in the brain.

Methods

Experimental study

Sixteen young healthy subjects, 9 males and 7 females, age 28.3±4.6 years were studied while resting on a specially designed chair with their arm supported by armrests during a period of 35–40 minutes. The population was selected by considering healthy young subjects in order to characterize the variability pattern in a rather homogeneous group. Patients were studied in a seated position as we wanted to resume the standard posture in which FRC is normally evaluated. Chest wall volume (Vcw) was measured plethysmographically by OEP (Fig. 1A), a technique extensively validated in several conditions and postures [9–12,14,15]. Briefly, 89 reflective markers were placed over the thorax from the clavicles to the anterior superior iliac spines along pre-defined vertical and horizontal lines. The detailed position of the markers is described in detail elsewhere [9]. The three-dimensional position of each marker was measured by an automatic motion analyzer (ELITE System, BTS, Milano, Italy) at a frequency of 10 Hz using four special video cameras, two in front and two behind the subject. Chest wall volume was determined by approximating the chest wall surface by triangles connecting the markers and computing the volume enclosed by all these triangles. The use of a subset of the markers allowed us to measure separately the volume of three chest wall compartments subdividing Vcw into the pulmonary rib cage volume (V_RCp), defined as the part of the chest wall between the clavicular line and the horizontal line at the level of the xiphoid process, the abdominal rib cage volume (V_RCa), defined as the volume comprised between the xiphoid process and the costal margin and the abdominal volume (V_AB), defined as the volume comprised between the costal margin and the anterior superior iliac spine. This technique was proven to be reliable when used to measure total tidal volumes and all the other ventilatory parameters. Moreover, because it is based on the direct measurement of volumes, it can be used to track breath-to-breath changes in end-expiratory lung volume (EEV_L)[12]. In particular, it has been shown that the tidal volume of the total chest wall measured by OEP is extremely similar to the tidal volume measured by spirometry, with the coefficient of variations of the differences between techniques being below 5% [9], Therefore we assume that all breathing pattern parameters derived by the total chest wall volume measured by OEP accurately resembles the ones measured by spirometry or body plethismography. In this study, from the time course of Vcw all classical breathing pattern parameters were computed. In particular, we considered the total and compartmental chest wall tidal volumes (V_Tcw, V_Trc,p and V_Tab), the total and compartmental chest wall end expiratory volumes (EEVcw, EEVrc,p, EEVab), the total respiratory time (T_TOT) with its components for the inspiratory (T_I) and expiratory (T_E) phases, the duty cycle (T_I/T_TOT), the respiratory rate (RR) and the intra-breaths interval (IBI). The combination of these time-based parameters with volume changes allowed also the assessment of the minute ventilation (V_E) and mean inspiratory flow (V_I/ T_I).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. A): Experimental set-up for continuous recording of total and compartmental chest wall volumes with the patient in the seated position.

A specially-designed chair allowed the measurement of the three-dimensional position of the passive optical markers for OEP placed on the back of the subject even if the subject was laid on the seatback. B): Schematic representation of the model. The brain respiratory rhythm generator includes five neuron groups, expiratory (Exp), early inspiratory (E-I), ramp-inspiratory (R-I), late-inspiratory (L-I) and post-inspiratory (P-I). Each group receives inhibitory signals (represented by a thick line at the end of the link between two groups) and the L-I group receives an excitatory input (filled circle). Note that the R-I group includes self-excitation. Pulmonary stretch receptors (PSR) are represented as an additional group of neurons receiving inputs from the R-I and sending outputs to the L-I and E-I groups. To mimic the observed physiological fluctuations in respiratory parameters such as tidal volume (V_T) and respiratory times (T_TOT = T_I+T_E), the R-I group also receives a uniformly distributed neural noise and its output drives the respiratory muscles which, in turn, act on the passive structures of the respiratory system (lung and chest wall) represented by a viscoelastic Kelvin body to develop the tidal breathing. (See text for further explanation).

https://doi.org/10.1371/journal.pone.0116317.g001

In order to assess the repeatability of the correlation properties of various parameters in the same subject, the measurements were repeated in a subset of 12 subjects from 3 to 9 months after the first measurement. The protocol was approved by the institutional review board of Fondazione Don Carlo Gnocchi, Milano and written informed consent was obtained from each subject.

Data analysis

Long-range correlation properties by DFA.

First, from the time courses of Vcw, all end-expirartory (EE) and end-inspiratory (EI) points were identified and used to compute all the breathing pattern parameters and end-expiratory volumes for each breath and to build, for each parameter, a time series as a function of breath number. These time series were analyzed by the DFA as suggested by Peng et al. [13]. The root-mean-square fluctuation of the integrated and detrended time series of y(k) of a time series was computed as follows: To calculate F(n) the time series y(k) containing N data points, the original time series were first integrated: where x(i) is the ith value of time series of one of the respiratory variables and is the corresponding average of the time series. The time series y(k) was then divided into non-overlapping windows of equal length (n). A linear regression line was fit through the data points of y(k) in each window. The regression line y_n(k) established the local trend in that window. The time series y(k) was then detrended by subtracting the local trend, y_n(k), from the data in each window. The calculation of F(n) was repeated for different n and plotted as a function of n on a log-log plot. When F(n) shows a linear increase on the log-log plot, then F(n) is said to follow a power-law functional form: where α is the scaling exponent and A is the amplitude of the power-law fluctuation function. These parameters can be obtained as the slope and intercept, respectively, of a straight line fit through the data plotted on a double logarithmic graph.

For a random process, α has the value of 0.5. For a positively correlated signal (large fluctuations are likely to be followed by large fluctuations), α is >0.5, and for an anticorrelated signal (large fluctuations are likely to be followed by small fluctuations), α is between 0 and 0.5 [13]. If F(n) follows a power law over at least an order of magnitude time scale with an α different from 0.5, the corresponding variable is said to exhibit long-range correlations or scale-invariant behaviour. For example, α = 1 corresponds to 1/f noise and α = 1.5 corresponds to Brownian noise.

To ascertain that the correlations in breath-to-breath fluctuations of the parameters are real, we randomized (shuffled) the order of the original breath-to-breath time series. Such a rearrangement of the data results in an uncorrelated time series. Thus, although this procedure does not alter the distribution of the amplitudes in the time series, the correlated ordering should disappear, and hence α of the shuffled time series should be 0.5, allowing to establish the existence of long-range correlations for a given variable if the value of α is significantly different from that after shuffling.

To assess the sensitivity of α to the number of breaths analysed (i.e. the length of the recording), we also estimated α by considering different numbers of breaths. We produced 5 different time series for each of the EEV recording by including all the breaths from the first to the 2^N-th breath and computed α for each time series reporting α as a function of N.

Neuromechanical model of lung function fluctuations.

In order to account for the correlated behavior of the data, we modified the model of Cernelc et al. [4], which was able to account for the observed long-range correlations in breath-to-breath fluctuations of V_T, O₂ and CO₂ in babies. The model is based on the original neural oscillator model proposed by Botros and Bruce [8]. The neural oscillator consists of five coupled non-linear ordinary differential equations corresponding to the activities of five neuron groups in the brain respiratory centre (Fig. 1B). These neuron groups receive tonic inputs (TNI) from the periphery and other brain centres and the network, combined with the TNIs, generates a rhythmic output. Specifically, the ramp-inspiratory neuron group provides periodic outputs to the phrenic nerve which, in turn, activates the respiratory muscles. Additionally, we also implemented a feedback loop accounting for the pulmonary stretch receptors (Fig. 1B) that receive input from the ramp-inspiratory group and send excitatory and inhibitory signals to the late-inspiratory and early-inspiratory groups, respectively. Table 1 shows the parameters of the full model.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Parameter values defining the respiratory brain oscillator consisting of 5 neuron groups and the pulmonary stretch receptors.

https://doi.org/10.1371/journal.pone.0116317.t001

We first solved the differential equations of the network in the time domain using MATLAB (Mathworks Inc, Natick, MA, USA). Following a short transient, the solution of the network becomes strictly periodic without any irregularities. To mimic the observed irregularities, we added noise to the TNI of the first or ramp-inspiratory neuron group (TNI₁), based on considerations of Hoop et al. [16] suggesting that neural noise varies within the respiratory cycle most likely due to varying chemoreceptor and stretch receptor responses. Hoop et al. [16] found correlations in the neural noise itself. However, in order to test whether the dynamic nonlinear oscillator alone is able to generate long-range correlations, we added uncorrelated random noise to the input of the oscillator. The mean value of TNI₁ was 3 with a uniformly distributed noise between 1 and 5 (standard deviation, SD = 1.15) which changed on average 1.6 times within the respiratory cycle. Next, we simulated ~500 breaths by solving the 6 coupled nonlinear differential equations using MATLAB. The negative parts of the cycles corresponding to expiration were replaced by zeros, the signal was then resampled at equidistant time intervals and was further used to drive a mechanical model of the respiratory system.

The mechanical properties of the chest and respiratory system were modelled by using a Kelvin body which is a parallel combination of a spring (E₁ = 1) with an in series combination of another spring (E₀ = 1) and a dashpot (R₀ = 40). This model has been used to describe the mechanics of the respiratory system [17,18]. The output of the ramp inspiratory group was assumed to provide a pressure which was used to drive the Kelvin body. The volume displacement of the Kelvin body, which corresponds to tidal volume, was then solved in the time domain in MATLAB for several combinations of the parameters. When the driving pressure increases during inspiration, the volume displacement also increases with a certain time lag due to the viscoelastic time constant of the Kelvin body. Once inspiration ends, the Kelvin body passively relaxes and the volume displacement decreases. Depending on the neural noise reaching the brain neural oscillator within the respiratory cycle, the oscillator can slow down or accelerate even within a cycle. If the inspiration starts before the Kelvin body reaches its fully relaxed state, the next inspiration becomes superimposed on the relaxation curve. Consequently, the absolute displacement of the Kelvin body, which mimics EEV_L, also increases. Since the TNI₁ shows irregular fluctuations, the end-expiratory lung volume should also display apparently random cycle-by-cycle fluctuations.

Because there was no feedback from the Kelvin body to the oscillator, we could simulate various long time series of the ramp-inspiratory group which were subsequently used as inputs to the Kelvin body. Since the solution for the displacement of the Kelvin body turned out to be highly sensitive to how the pressure input was sampled, we used a stiff differential equation solver. To test the numerical algorithm, first the noise in the neural network was set to zero which provided a strictly periodic output. This signal was used to drive the Kelvin body and the parameters of the solver were then adjusted until no variation was seen in the simulated EEV_L. Noise was then added to the neural oscillator and the simulations were repeated for short and long time constants of the Kelvin body. This procedure allowed us to examine the correlation properties of the time series of T_TOT, V_T as well as EEV_L using the DFA algorithm.

Statistical analysis.

Data are reported as mean ± SD and were compared using one-way or two-way ANOVA. When the normality or equal variance tests failed, either non-parametric tests were used or the data were transformed before the parametric test. Significance was accepted at the level of p<0.05.

Results

For each subject, we recorded approximately 500 breaths. Fig. 2 shows examples of the time course of the total chest wall volume (Vcw) and its pulmonary rib cage (Vrc,p) and abdominal (Vab) components for a representative subject. It can be seen that the fluctuations in volume are noticeably larger in the abdominal than in the rib cage compartment. The time series of several ventilatory parameters were then obtained by detecting EE and EI points and computing the value of each parameter on a breath-by-breath basis (Fig. 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Experimental tracings of total and compartmental chest wall volumes recorded continuously on a representative subject by Opto-Electronic Plethysmography (only 16 min of tracings are reported for clarity).

Vcw: total chest wall volume; Vab: abdominal volume; Vrc,p: volume of the pulmonary rib cage.

https://doi.org/10.1371/journal.pone.0116317.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Time series of the different breathing patter parameters computed from the volume recordings on a representative subject.

https://doi.org/10.1371/journal.pone.0116317.g003

The end-expiratory volume of the abdominal compartment (EEVab) showed a significantly larger variability than that of the pulmonary rib cage (EEVrc,p). Since also the magnitude of the changes of EEVab was larger than for the EEVrc,p, the fluctuations in the end-expiratory volume of the total chest wall (EEVcw) were very mostly dominated by those in EEVab. The tidal volumes of the two chest wall compartments and the total chest wall as well as the inspiratory flow, duty cycle and the minute ventilation all showed qualitatively similar fluctuations although somewhat different from EEVab.

The DFA plots corresponding to the time series in Fig. 3 are summarized in Fig. 4. It can be seen that except for minute ventilation (V_E), the fluctuation function of each ventilatory parameter in this subject increases linearly on a double logarithmic graph through about 2 decades of breath numbers. For this representative subject, the values of the exponent α obtained as the slope of the linear regression fits to the data were higher than 1 for all end-expiratory volumes. The values of α were between 0.58 and 0.65 for all other ventilatory parameters except V_E suggesting that EEV was considerably more correlated than the other parameters. With regard to V_E, when all 7 points on the DFA plot were included in the regression, α was 0.49 (implying a random process) with an r² = 0.972. However, it is clear that the last two points deviate from a straight line. By removing the last two points corresponding to the longest timescales from the regression, α became 0.63 with an r² = 0.99. Regarding all the subjects, we observed a similar trend (i.e. r²<0.98) in V_E in 4 and in V_T/Ti in 1 of the 16 subjects. The subject showing a low r² for the DFA fit of V_T/Ti was also one of the 4 subjects showing low r² in V_E. In these cases, we excluded the last 1 or 2 points from the linear regression to maximize r². In 2 of the 4 subjects, r² did not increase above 0.98 and hence the corresponding α values were excluded from the statistics.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Fluctuation (F) as a function of window size (n) for the time series from volume tracings of a representative subject.

Each panel refers to a specific variable of the breathing pattern. See text for abbreviations. The regression line determining α is also reported.

https://doi.org/10.1371/journal.pone.0116317.g004

Table 2 summarizes the mean±SD values of the average, the coefficient of variation (CV) and α of the ventilatory parameters EEVcw, EEVrc,p, EEVab, V_Tcw, V_Trc,p and V_Tab for all subjects. The CV of EEVab was statistically significantly larger than that of the EEVcw and EEVrc,p. The larger variability was also accompanied by stronger correlation properties of EEVab because the corresponding α was statistically significantly larger than the α of EEVcw or EEVrc,p. Interestingly, with regard to tidal volume, the largest variability was seen in the rib cage compartment which had a significantly larger CV than the total chest wall or the abdominal compartment. The largest α for the tidal volumes was observed in V_Tab which reached statistical significance compared to the α of V_Tcw. When the experiments were repeated on a second day 3 to 9 months later, no differences were found between any of the parameters.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2.

https://doi.org/10.1371/journal.pone.0116317.t002

Table 3 summarizes the mean±SD values of the average, CV and α of additional ventilatory parameters related to the respiratory period. The T_I/T_TOT showed considerably smaller variability than the other parameters. Indeed, the CV of T_I/T_TOT was statistically significantly smaller than that of the other parameters. The value of α for T_TOT was significantly larger than α for T_I/T_TOT or V_E. Again, on the second day, the parameter values were statistically the same as on the first day suggesting a good repeatability of the measurements.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3.

https://doi.org/10.1371/journal.pone.0116317.t003

To summarize, all ventilatory parameters showed significant fluctuations and had a correlation exponent higher than 0.5, implying the presence of long-range correlations in their fluctuations. While the end-expiratory lung volumes had the smallest variability with respect to the average, they also showed the strongest correlations. In particular, the abdominal compartment had the largest exponent implying the strongest control on the fluctuations over long time scales. On the other hand, the least amount of correlation (smallest value of α) was seen in the T_I/T_TOT while the parameters showing the largest coefficient of variation were the total and compartmental V_T.

The repeatability data are reported as Bland and Altman graphs in Fig. 5A. The average difference in α for the EEV_L time series is 0.037±0.107, with no evident dependence of the difference on the mean value. The sensitivity of α as a function of the length of the recording is also shown in Fig. 5B. The results demonstrate that α approaches a stable value already for time series of length 64.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. a) Bland-Altman analysis for the repeatability of the chest wall end-expiratory (EEVcw, upper panel) and tidal (V_T,cw, lower panel) volumes long term correlations.

Each data point represent the value of α measured in the same subject during two different experimental session with from 3 to 9 month time interval between the measurements. b) α values as a function of the number of breaths considered for the analysis for total and compartmental end-expiratory chest wall volumes.

https://doi.org/10.1371/journal.pone.0116317.g005

In order to better understand the origins of correlation properties of the breath-by-breath fluctuations in respiratory parameters, we simulated the fluctuations in V_T, T_TOT, and EEV_L using the Kelvin body driven by the respiratory neural oscillator network model. Examples of the corresponding time series are shown in Fig. 6A. Fluctuations in V_T of the model are most similar to those in V_Tcw (Fig. 3), whereas those in T_TOT and EEV_L are somewhat reminiscent of the corresponding fluctuations in the abdominal compartment. Fig. 6B shows the fluctuation functions for the model outputs on double logarithmic graphs. The model-based α for V_T is 0.571 overlapping with the experimental values of α for the chest wall (0.669 ± 0.13). The model correctly predicted both the α for T_TOT (0.694 vs 0.702 ± 0.1) and the α for EEV_L (1.045 vs 1.064 ± 0.136).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Left panels: time series derived from the simulated volume tracings provided by the model for tidal volume (V_T, upper panel), respiratory time (T_TOT, middle panel) and end-expiratory lung volume (EELV, lower panel).

Right panels: Fluctuation (F) as a function of window size (n) for the time series from volume tracings simulated by the model (V_T, upper panel, T_TOT, middle panel, EELV, lower panel).

https://doi.org/10.1371/journal.pone.0116317.g006

Discussion

The most important finding of this study is that the long-term correlation properties of EEV parameters are significantly different (greater) than those of all other respiratory indices such as V_T or T_TOT.

Even though several studies investigated the correlation properties of breathing pattern[2–5,7], only few addressed the correlation properties of EEV such as functional residual capacity (FRC) in humans, likely because of the difficulties in continuous monitoring. To our knowledge, the only study addressing this issue is by Hlastala et al. [19], in which the authors used a complicated set-up based on total body plethysmography in supine position to record FRC fluctuations and other respiratory parameters in four healthy subjects for 40 min. The experimental data reported, i.e. representative tracings of time series of breathing pattern parameters and FRC are quite similar to ours, showing much greater long time scale fluctuations for FRC compared to V_T and other respiratory parameters. In that study, however, the analysis of the fluctuations was limited to autocorrelation and spectral analysis, which showed major oscillations of FRC with periods from 8.3 and 28 min with amplitudes varying from 42 to 176 ml.

In this study, the use of OEP allowed us to record accurate time series of FRC in a way that minimally perturbed natural breathing, as this technique does not require any connection to the airway opening and, therefore, no external load or dead space is added to the subject’s respiratory system. The accuracy of this technique for the long-term measurement of changes in EEV has been validated [12], showing the absence of any significant drift on long time scales. Therefore, the data presented in this study accurately reproduces natural variations in EEVs in conscious humans.

The long term correlations of EEV_L were studied by Zhang et al. [20] in anaesthetised tracheotomised rats exposed to either constant positive or negative airway pressure using the Hurst exponent (H). They found a value of 0.55 for H when positive airway pressure was applied which statistically significantly increased to 0.64 when negative airway pressure was applied. Interestingly, the difference in H between the positive and negative pressure conditions vanished if the rats were vagotomised, with H approaching 0.55 in all conditions. The authors concluded that although the exact mechanisms for the fractal behavior of EEV regulation are unknown, the results suggest that interaction of the pulmonary mechanoreceptor activity with the central respiratory pattern generator might contribute to the observed long-term correlation. Furthermore, in a subsequent study [21], the same authors found similar H values in rats with intact airways, supporting their previous conclusion. Since H is equal to α for fractional Gaussian noise, the nature of the long term correlations in anaesthetised rats appears to be much less correlated than in conscious humans in the present study. There are several possible reasons for the difference. For example, correlation properties of V_T and T_i have been shown to depend on sleep stage [5]. Additionally, the α of respiratory resistance at 8 Hz in human subjects decreased from wake to sleep [5] suggesting less active brain control of respiratory variables during sleep.

In order to interpret our experimental results, we extended the neural network model of Cernelc et al. [4] to include pulmonary stretch receptor feedback and the mechanical load of the respiratory system on the respiratory muscles. Without the pulmonary stretch receptor feedback (Fig. 1B), the model could predict α for V_T and EEV_L, but it provided a lower value for T_TOT (α = 0.534). With regard to the mechanical load, we hypothesized that while timing of the respiratory cycle and V_T are mostly determined by the respiratory oscillator and the inputs to the oscillator, EEV_L volume would also be influenced by the passive mechanical properties of the tissues. The reason is that, although some expiratory muscles are active during expiration and may also be active during inspiration for smoothing the shape of tidal volume or maintaining the posture [22], changes in lung volume are mostly determined by the passive mechanical properties of the respiratory tissues. While the lack of active expiration in the model is a limitation, the simulations are in reasonable agreement with the measured data. Thus, in the absence of active expiration, the tissues are stretched at end-inspiration and the stored elastic energy in the tissues lower lung volume via relaxation. Because the tissues are viscoelastic, this is similar to the discharge of a capacitor through a resistor. Thus, we used the simplest viscoelastic model, the Kelvin body, to mimic both inspiration and expiration. During inspiration, the Kelvin body simply acts as a filter and slightly shapes the inspiratory volume and therefore it has little effects on the correlations of V_T and T_TOT. Indeed, simulations without the Kelvin body changed α for both variables by less than 1.5%. Since during expiration the oscillator input to the Kelvin body is zero, the effect of the load on the α of EEV_L is substantial.

The assumptions used in and the limitations of the neural network model itself have been discussed in previous studies [4,7] and hence we focus here on the limitations of the mechanical model of the respiratory system. This model is a simple Kelvin body that has been used to fit the mechanical properties of the respiratory system in healthy human subjects [17]. The Kelvin body is a simple linear viscoelastic model with a dashpot and two springs. The dashpot with a spring in series mimics the viscoelastic properties of the tissues whereas the parallel spring gives the model a solid-like behaviour: following a step change in force (or pressure), the displacement (or volume) approaches a finite constant during creep and the time course of the creep is characterized by a single exponential [23]. However, frequency domain measurements of the respiratory tissues including the chest wall and the lung tissues have unequivocally proved that the tissues have a broad spectrum of time constants and this has led to the concept of the so called “constant phase” model [24]. The single exponential of the Kelvin body cannot account for the broad distribution of time constants of real tissues and hence the viscoelastic memory of this model can only describe tissue behaviour over a limited time scale. The constant phase model has a power law type of time constant distribution and its memory, characterized by the relaxation of the creep function, also has a power law form [25]. Thus, white noise passing through such a system generates long-range correlations and the corresponding DFA plot is a straight line over many decades of time on a double logarithmic graph. As seen in Fig. 5, however, the DFA plot for the EEV_L obtained from our model is curved because the memory function of the Kelvin body does not follow a power law and this is a limitation of the model. Indeed, for short time scale, a local α = 1.4 fits the first 4 points of simulated EEV_L time series whereas for long time scales a weaker correlation is seen with α = 0.7 as a best fit for the last 4 points. Furthermore, sensitivity analysis also showed that the longer the time constant of the Kelvin body, the higher the α of the simulated EEV_L. For example, a 10 times smaller time constant lowered α of the simulated EEV_L to 0.7.

Another limitation of the model is that the Kelvin body is linear, whereas the respiratory tissues exhibit non-linear mechanical properties [26]. An advantage of the Kelvin body model is that it is easy to simulate force-displacement data in the time domain with or without nonlinearities. Therefore, we tested the effects of mild nonlinearities on the output of the Kelvin body. The nonlinearity was an additional contribution (F_n) to the total force output of the Kelvin body in the form of a second order term proportional to the square of the displacement: F_n = E₂*x² where E₂ is the nonlinear spring and x is the displacement of the Kelvin body. We found that with E₂/E₁<0.2 there was no appreciable effect of nonlinearity on the simulated time series. For E₂/E₁ = 0.5, α of the simulated EEV_L decreased by 15% but the exponent of T_TOT and V_T did not change. Since the time constant R₀/E₀ and the nonlinearity had opposing effects on α of EEV_L, it was not possible to uniquely fit the Kelvin body to measured data. However, the purpose of the modelling was not to estimate parameters of individual subjects; rather, to examine whether the neural network together with a mechanical load is capable of qualitatively accounting for the measured correlations.

Despite the simplicity of the model, it provided correlated fluctuations with exponents that are in reasonable agreement with the experimentally measured exponents. These correlated fluctuations come from the brain but the input from the periphery is crucial for maintaining the oscillations [8]. On the other hand, without fluctuations in the inputs to the neural oscillator, the output of the oscillator would be a deterministic cyclic pattern, a limit cycle, with constant frequency and amplitude. One interesting implication of the modelling results is that the addition of the Kelvin body did not significantly affect the pattern including timing (T_TOT). This suggests that while the correlated properties of EEV_L are determined both by the brain oscillator and the viscoelastic memory of the tissues, breathing patterns originates almost exclusively from the brain and the neural noise.

There is another important implication of the fluctuations in V_T for cell function. Since the oscillator is a nonlinear network of neuron groups, any fluctuation coming from outside—including the lung and other brain regions—will generate cycle-by-cycle variations both in the rate and amplitude of the output which is assumed to be the force generated by the respiratory muscles. Consequently, all cells in the respiratory system are exposed to cycle-by-cycle fluctuations in stretch due to the variability in V_T. Since such fluctuations must have been present for hundreds of millions of years when the earliest mammals evolved, cells in the respiratory system must have adapted to such variabilities in stretch. One example of this is that surfactant secretion by alveolar type II epithelial cells is sensitive to stretch pattern and cells receiving variability in stretch secrete more surfactant both in vivo [27] and in vitro [28]. Therefore, we conclude that the origin of correlated fluctuations observed in breathing parameters is the brain respiratory center that receives inputs from other brain centers and the periphery. Furthermore, the cycle-by-cycle fluctuations in tidal volume and respiratory times likely have a significant impact on general functioning of stretch-sensitive adherent cells in the respiratory system.

In summary, we measured the cycle-by-cycle variations of several respiratory physiological variables and found that these fluctuations show long-range correlations that are highly reproducible in individual subjects. Using a novel neuromechanical model, we propose that the correlations in the timing and amplitude of the physiological variables originate from the brain with the exception of end-expiratory lung volume which shows the strongest correlations due to the contribution of the viscoelastic properties of the tissues. Finally, we also suggest that since cells in the respiratory system are exposed to fluctuation in cycle-by-cycle stretch related to variability in tidal volume, these finding may have implications on general cell function in the respiratory system.