Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Neural Mechanisms Underlying Breathing Complexity

  • Agathe Hess,

    Affiliations Laboratoire Matière et Systèmes complexes, UMR 7057, CNRS, Université Paris 7, Paris, France, Service de Radiologie, APHP, Hôpital Bichat-Claude Bernard, Paris, France

  • Lianchun Yu,

    Affiliations Laboratoire Matière et Systèmes complexes, UMR 7057, CNRS, Université Paris 7, Paris, France, Institute of Theoretical Physics, Lanzhou University, Lanzhou, China

  • Isabelle Klein,

    Affiliations Service de Radiologie, APHP, Hôpital Bichat-Claude Bernard, Paris, France, Unité Inserm 698, Université Paris 7, Paris, France

  • Marine De Mazancourt,

    Affiliations Laboratoire Matière et Systèmes complexes, UMR 7057, CNRS, Université Paris 7, Paris, France, Ecole Normale Supérieure, Paris, France

  • Gilles Jebrak,

    Affiliation Service de Pneumologie B, APHP, Hôpital Bichat-Claude Bernard, Paris, France

  • Hervé Mal,

    Affiliation Service de Pneumologie B, APHP, Hôpital Bichat-Claude Bernard, Paris, France

  • Olivier Brugière,

    Affiliation Service de Pneumologie B, APHP, Hôpital Bichat-Claude Bernard, Paris, France

  • Michel Fournier,

    Affiliation Service de Pneumologie B, APHP, Hôpital Bichat-Claude Bernard, Paris, France

  • Maurice Courbage,

    Affiliation Laboratoire Matière et Systèmes complexes, UMR 7057, CNRS, Université Paris 7, Paris, France

  • Gaelle Dauriat,

    Affiliation Service de Pneumologie B, APHP, Hôpital Bichat-Claude Bernard, Paris, France

  • Elisabeth Schouman-Clayes,

    Affiliation Service de Radiologie, APHP, Hôpital Bichat-Claude Bernard, Paris, France

  • Christine Clerici,

    Affiliations Département de Physiologie-Explorations fonctionnelles, APHP, Hôpital Bichat-Claude Bernard, Paris, France, Unité Inserm 700, Université Paris 7, Paris, France

  • Laurence Mangin

    laurence.mangin@bch.aphp.fr

    Affiliations Laboratoire Matière et Systèmes complexes, UMR 7057, CNRS, Université Paris 7, Paris, France, Département de Physiologie-Explorations fonctionnelles, APHP, Hôpital Bichat-Claude Bernard, Paris, France, Centre d’Investigation Clinique APHP, Hôpital Bichat, Paris, France

Abstract

Breathing is maintained and controlled by a network of automatic neurons in the brainstem that generate respiratory rhythm and receive regulatory inputs. Breathing complexity therefore arises from respiratory central pattern generators modulated by peripheral and supra-spinal inputs. Very little is known on the brainstem neural substrates underlying breathing complexity in humans. We used both experimental and theoretical approaches to decipher these mechanisms in healthy humans and patients with chronic obstructive pulmonary disease (COPD). COPD is the most frequent chronic lung disease in the general population mainly due to tobacco smoke. In patients, airflow obstruction associated with hyperinflation and respiratory muscles weakness are key factors contributing to load-capacity imbalance and hence increased respiratory drive. Unexpectedly, we found that the patients breathed with a higher level of complexity during inspiration and expiration than controls. Using functional magnetic resonance imaging (fMRI), we scanned the brain of the participants to analyze the activity of two small regions involved in respiratory rhythmogenesis, the rostral ventro-lateral (VL) medulla (pre-Bötzinger complex) and the caudal VL pons (parafacial group). fMRI revealed in controls higher activity of the VL medulla suggesting active inspiration, while in patients higher activity of the VL pons suggesting active expiration. COPD patients reactivate the parafacial to sustain ventilation. These findings may be involved in the onset of respiratory failure when the neural network becomes overwhelmed by respiratory overload We show that central neural activity correlates with airflow complexity in healthy subjects and COPD patients, at rest and during inspiratory loading. We finally used a theoretical approach of respiratory rhythmogenesis that reproduces the kernel activity of neurons involved in the automatic breathing. The model reveals how a chaotic activity in neurons can contribute to chaos in airflow and reproduces key experimental fMRI findings.

Introduction

Complexity is a universal phenomenon widely described in physics as well as in living organisms in biology and physiology. In the human brain, neural networks are complex [1] and communication between neurons occurs through a wild variety of codes such as bursting oscillations, which is a brief epoch of rapid firing. Such bursting behavior of the neuron oscillations may exhibit nonlinear deterministic chaos [2]. The human respiratory system displays several level of complexity: the bronchial tree has a fractal structure with various degrees of self-similarity and the airflow dynamics inside exhibits chaos during rhythmic breathing [3]. Why rhythmic breathing generates chaos in human airflow remains unknown. Breathing is maintained and controlled by a network of neurons in the brainstem that generate respiratory rhythm while receiving regulatory inputs. Pace-maker like neurons generating rhythmic breathing have been identified in 2 brainstem regions in rodents, one located in the rostral ventro-lateral (VL) medulla, the pre-Bötzinger complex [4][8], and the other close to this region, the parafacial respiratory group [9][13]. Recent evidence suggests that both groups of neurons are coupled oscillators that work in tandem to synchronize respiratory rhythm [9], [10], [13]. Moreover, these automatic neuronal groups have two important properties: they are capable of different synchronization regimes depending on the level of their excitabilities [13] and their dynamics exhibit chaotic spike-bursting oscillations in some circumstances [14]. Indeed, neural population activity recorded locally in the pre-Bötzinger complex of neonatal rat brainstem slices exhibit chaotic dynamics, when neuronal excitability is progressively elevated [14]. This is a strong argument to hypothesize that the chaos-like complexity of airflow in humans is an intrinsic property of central respiratory generators. In addition, both respiratory rhythm and airflow control have common genetic determinants [15]. However, breathing is also modulated by the state of airways [16], by the chest wall [17], the lung, by chemical afferents sensitive to hypercapnia, hypoxia or acidosis [3] and by mechanical afferents from the airway, lung, chest wall, respiratory muscles as well as by supra-pontine commands. A previous study has shown that the structural and mechanical properties of the bronchial tree, lung and chest wall in humans are not sufficient to generate chaos in airflow in the absence of a central neural drive [18]. Nevertheless, it is still unclear in humans to what extent the complex dynamics of the respiratory center contributes to airflow complexity.

We used both experimental and theoretical approaches to decipher the brainstem neural substrates of ventilatory complexity in humans. Complexity of airflow was estimated during inspiration and expiration at rest, and during an inspiratory effort with resistive load, used as an indirect neural stimulus. Brainstem regions of interest of the respiratory pacemakers were located with fMRI [19] in the rostral ventro-lateral medulla containing the pre-Bötzinger complex, and in the caudal ventro-lateral pons containing the parafacial group. Our goal was to evidence brainstem neural correlates of airflow complexity. We also analyzed airflow in a disease state in patients with chronic obstructive pulmonary disease (COPD). COPD is the most frequent chronic lung disease in the general population and is mainly due to tobacco smoke. Patients with COPD have an impaired lung function with an increased respiratory load due to small airways obstruction by inflammation and remodeling. Lung parenchyma destruction or emphysema is often associated with distal obstruction. Airflow obstruction associated with hyperinflation and respiratory muscles weakness are key factors contributing to load-capacity imbalance and hence increased respiratory drive [20]. At the end stage of the disease, the patients have respiratory insufficiency with home oxygen therapy while the neural respiratory drive is extremely high. We hypothesized that chaos in airflow should be altered in COPD patients but that such alterations should still correlate with the activity of the brainstem respiratory centers. Further, we developed a mathematical model of respiratory rhythmogenesis to reproduce the basic activity modes of neurons involved in the automatic breathing in healthy subjects and COPD patients. The model therefore reveals how a chaotic activity in neurons can contribute to chaos in airflow and reproduces key experimental fMRI findings.

Results

The characteristics of the whole population, healthy subjects and patients with chronic obstructive pulmonary disease (COPD), are shown in Tables 1 and S1. No difference was noted in end-tidal PCO2 (PETCO2) measurements between healthy subjects and COPD patients either during resistive load or during resting state fMRI (Figure 1).

thumbnail
Figure 1. Chaos characterization of airflow during inspiration (Vt/Ti) and expiration (Vt/Te) in the controls and COPD patients.

A: Noise Limit value (%), B: largest Lyapunov exponent, C: correlation dimension. The boxes encompass the interquartile range with indication of the median, the whiskers delimit the 95th percentile of the data distribution. Paired and unpaired Ttest.

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

Chaos in Airflow during Inspiration is Higher than during Expiration in Healthy Subjects

Linear and nonlinear measurements of the airflow.

The linear estimates (coefficient of variation (CV) and autocorrelation coefficient (AC)) of the airflow during inspiration and expiration are shown in Table S2. In the 25 healthy subjects, inspiratory flow yields higher variability (p<0.001) and lower value of the AC (p<0.001) than expiratory flow during unloaded breathing.

The number of time series that exhibits a positive noise limit value characterizing chaos in airflow is equivalent for inspiration and expiration (Table S3). In the time series with positive noise limit, chaos in airflow is increased during inspiration as compared with expiration (largest Lyapunov exponent (LLE) and the correlation dimension (CD), p<0.05) (Figure 1). The attractor of the airflow is reconstructed in the phase plane during inspiration with the corresponding time series in one healthy subject (Figure 2A).

thumbnail
Figure 2. The chaotic signatures of the airflow in one healthy subject and one COPD patient are evidenced.

The reconstructed attractors in the phase plane are shown on the left panel for one healthy subject during inspiration (A) and for one COPD patient during expiration (B). The corresponding time series are shown on the right panel.

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

Cerebral fMRI results.

In healthy subjects, we found that neural activity assessed in terms of the amplitude of low frequency oscillations (AlFO) of the BOLD signal located in the VL medulla is significantly higher than neural activity of the VL pons (p<0.001, n = 16) (Figure 3top, Figure 4). In COPD patients, the AlFO of the BOLD signal located in the VL pons, which contains the parafacial group, is significantly higher than the ALFO of the VL pons of healthy subjects (p<0.001, n = 16) (Figure 3top).

thumbnail
Figure 3. fMRI results of the brainstem respiratory centers at rest.

Top. Amplitude of the low frequency oscillations (AlFO) of the resting state BOLD signal computed in controls and COPD patients. In healthy subjects the AlFO of the rostral ventro-lateral (VL) medulla that contains the pre-Bötzinger complex is higher than the VL medulla of the patients. Conversely, the ALFO of the caudal (VL) pons, which contains the parafacial respiratory group is higher in patients than the VL pons of healthy subjects. Bottom. Simulated AlFO obtained after hemodynamic convolution of the theoretical neural states. For controls, the chosen network scheme is described in Figure 8B, while for COPD patients, the network scheme is described in Figure 8C. Of note the synchronization regime describe in Figure 8D gave the same results as 8C for the simulated AlFO of the BOLD signal.

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

thumbnail
Figure 4. Localization on fMRI images of the regions of interest, the rostral ventro-lateral (VL) medulla with the pre-Bötzinger complex and the caudal VL pons with the parafacial.

The regions of interest (brain mask with 4 cubes) were computed based on the recent article of Schwarzacher et al. (A) on individual standard images (sagittal and axial slices in B). Then the coordinates of the regions of interest were transformed from standard space to functional space (sagittal and axial slices in C). Two regions of interest of the right VL medulla and pons are shown in the sagittal slice (red and yellow) while two regions of interest of the VL medulla (red and blue) are shown on the axial slice. The region of interest of the left VL pons is not shown. Finally the mean time series were extracted for subsequent analyses (D). The corresponding mean time series of the VL medulla and VL pons are shown after extraction from the functional images, preprocessing analyses and regressing out with physiological covariates. The oscillations of the fMRI BOLD signal of the medulla in healthy subject (time series in red) are higher than those of the pons (time series in yellow). A: anterior, P: posterior, R: right, L: left.

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

COPD Patients Breathe with a Higher Level of Complexity during Expiration than Healthy Subjects

Linear and nonlinear measurements of the airflow.

The linear estimates (CV and AC) of the airflow during inspiration and expiration are shown in Table S2. In the 25 patients with COPD, expiratory flow yields higher variability (p = 0.06) and AC (p<0.001) than inspiratory flow during unloaded breathing. The number of time series that exhibits a positive noise limit value characterizing chaos in airflow is equivalent for inspiration and expiration in COPD patients (Table S3). However, the number of chaotic time series during expiration is higher in COPD patients than healthy subjects (p = 0.001, Table S3). The attractor of the airflow is reconstructed in the phase plane during expiration with the corresponding time series in one COPD patient (Figure 2B).

In the time series with positive noise limit, chaos in airflow is increased during expiration as compared with inspiration (NL values, p = 0.05, Figure 1A). Moreover, as compared with controls, the levels of airflow complexity of expiration (NL value, p<0.001; LLE p<0.001; CD, p<0.01) as well as inspiration (LLE, p<0.05; CD, p<0.05) is higher (Figure 1B–C).

COPD patients having hypoxia (n = 10) do not exhibit differences from those being normoxic (n = 15) in terms of ventilatory complexity (noise limit value, largest Lyapunov exponent and correlation dimension). Furthermore, when comparing the chaotic indexes in the control group (n = 25) and in the COPD patients being normoxic (n = 15), (PETCO2 being equivalent for both group), significant differences are evidenced with the noise limit value (NL controls: 5±7, NL COPD: 13±12 p<001) and the largest Lyapunov exponent (LLE controls: 0.15±0.08, LLE COPD: 0.27±0.1, p<0.001) of the expiratory flow.

Besides, COPD patients that exhibit severe dyspnea (Borg scale) have a significant higher level of expiratory flow chaos (correlation dimension) and AlFO of the VL pons than those with moderate and mild dyspnea (Figure S2).

Airflow complexity correlates with cerebral fMRI BOLD signal.

Univariate analysis in the whole population shows that the NL and the LLE values of the expiratory flow both positively correlates with AlFO of the VL pons (R2 = 0.4, p = 0.05 and R2 = 0.5, p = 0.04 for the NL and LLE respectively), the higher the complexity of expiration, the higher the neural activity of the VL pons. There is also an inverse relationship between the NL and the LLE of the expiratory flow and the pulmonary function index FEV1/FVC in the whole population (R2 = 0.45, p<0.05; R2 = 0.5, p<0.05, respectively, Figure S3). In healthy subjects, the chaotic level (NL) of airflow during inspiration strongly correlates with the neural activity of the VL medulla (R2 = 0.75, p = 0.01). In COPD patients, chaos (NL) during expiration correlates with the neural activity of VL pons (R2 = 0.4, p = 0.03). No correlation was evidenced between complexity of airflow and oxygen or carbon dioxide arterial pressures (PaO2, PaCO2). Multivariate analysis in the whole population showed that both neural activity of the VL pons and pulmonary function FEV1/FCV significantly predict the chaos of expiration (R2 = 0.4, F = 5.2 with p<0.01): the lower the pulmonary function, the higher the neural activity of the VL pons, the higher the chaotic level of expiration (Figure 5).

thumbnail
Figure 5. Central neural correlates of airflow dynamics in the whole population using multiple linear regression.

Both the amplitude of the low frequency oscillations (ALFO) and the pulmonary function index FEV1/FVC significantly predict airflow complexity in the whole population: the lower the pulmonary function, the higher the value of the AlFO of ventro-lateral pons and the higher the complexity of expiration.

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

Airflow complexity and cerebral fMRI results during inspiratory load.

Loading inspiration significantly increases the variability of the inspiratory flow, and the AC of the inspiratory as well as expiratory flows in both healthy subjects and patients with COPD (Table S2). In healthy subjects, inspiratory resistance significantly reduces airflow complexity during inspiration (Figure 6). Interestingly in COPD patients, loading inspiration leads to a diminution of complexity of inspiration (NL, LLE, CD) as well as expiration (NL, CD). Of note, loading inspiration did not change the PETCO2 (Figure S1) and saturation of both populations. Loading inspiration in healthy subjects and COPD patients also leads to a diminution of fMRI BOLD responses in the VL medulla (healthy subjects and COPD) and pons (COPD) (Figures 7 and S4). During inspiratory loading in the whole population, the mean negative BOLD signal of the VL pons correlates with the CD of the expiratory flow (R2 = 0.7, p<0.01) while the mean negative BOLD signal of the VL medulla correlates with the LLE of the inspiratory flow (R2 = 0.6, p = 0.05).

thumbnail
Figure 6. Loading inspiration leads to a diminution of complexity in airflow during inspiration in healthy subjects and during both inspiration and expiration and COPD patients.

Noise Limit value, largest Lyapunov exponent and correlation dimension values are given from top to bottom. Lo no inspiratory load; L20: loading inspiration with 20 cmH2O/L/sec. The boxes encompass the interquartile range with indication of the median, the whiskers delimit the 95th percentile of the data distribution. Paired and unpaired Ttests.

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

thumbnail
Figure 7. Negative BOLD signal of the respiratory brainstem network during inspiratory resistive loading in healthy subjects and COPD patients.

Group analysis of healthy subjects (n = 16, left) and COPD patients (n = 17, right). Sagittal, coronal and axial slices are shown. In controls, negative BOLD signal is mainly evidenced in the ventro-lateral (VL) and dorsal medulla. In COPD patients, inhibition is located in the caudal lateral and dorsal pons, and in the lateral rostral medulla (color code in blue). A: anterior, P: posterior, R: right, L: left. Histograms showing the corresponding BOLD signal changes in the rostral medulla and caudal pons for controls and patients. C. The main coordinates (x,y,z) of the clusters that exhibit inhibitory BOLD signal are given in MNI space (Montreal neurological Institute).

https://doi.org/10.1371/journal.pone.0075740.g007

Of note, healthy subjects and COPD patients also exhibit positive BOLD signal in the activated brain regions known to be involved in the voluntary control of respiratory muscles, i.e. sensory-motor, premotor and supplementary motor cortex area (data not shown).

Comparison of the correlation dimension of the original time series with surrogates.

The correlation dimensions of the 137 experimental time series were compared with 5 surrogates (685 simulated time series) that match each original signal. Those surrogates were computed after assigning random phase. Significant differences were obtained between the original data paired with the corresponding average correlation dimension values from the matching surrogate (p<0.01, Wilcoxon signed-rank test), reinforcing the nonlinear features of the inspiratory and expiratory flows time series.

Mathematical Model of Respiratory Rhythmogenesis

The present model is the first attempt to reproduce respiratory rhythmogenesis in healthy humans and COPD patients with experimental data. The model considers two chaotic pacemakers, the inspiratory (Pre-Bötzinger) and expiratory (parafacial) generators that work together via chemical synaptic connection, either activated or inhibited, to synchronize the respiratory cycle. Different dynamics are evidenced depending on the excitability level of the neurons. In the model, the parameters J1 and J2 represent the excitability level of the parafacial and pre-Bötzinger respectively. Experimental results show that healthy subjects display more complexity during inspiration than expiration and that the low frequency oscillations of the BOLD signal located in the rostral VL medulla have higher amplitude than oscillations of the caudal VL pons. From this, we postulate that the pre-Bötzinger complex is highly likely more excitable than the parafacial group, and drives the respiratory rhythm (active inspiration). Simulation of this network scheme is shown in Figure 8 with two possible regimes depending on the parameter values J1 and J2. In the first regime (Figure 8A), the parafacial has a very low excitability and is entirely depressed with no action potential. This network scheme is similar to the one described in adult rats, the “no-handshake process” [13]. The corresponding attractor of this scheme entirely relies on the pre-Bötzinger dynamics (Figure 9A). In the second regime (Figure 8B), while the pre-Bötzinger is the dominant pacemaker still driving the respiratory cycle, the parafacial group is occasionally relieved by specific physiological conditions [21]. Experimental results show in COPD patients that airflow complexity is higher during expiration than inspiration and that the low frequency oscillations of the BOLD signal located in the VL pons have higher amplitude than the oscillations of the VL pons of healthy subjects. In patients, we therefore hypothesize that the expiratory neurons located in the VL pons are more excitable than the pre-Bötzinger and drive the respiratory cycle. In this network scheme (Figure 8C), the more excitable parafacial group triggers the pre-Bötzinger which in turn inhibits the parafacial with a post-inhibitory rebound burst. The parafacial then switch-off inspiration. This network scheme is similar to the “full-handshake process” described in neonatal rats [13]. The corresponding attractor of this synchronization process mainly relies on the parafacial neurons dynamics (Figure 9C). Another synchronization regime may coexist in the disease state, when the excitability level of the expiratory group is slightly lower: the “half-handshake” process in which the parafacial still triggers the pre-Bötzinger which in turn induces a delayed post-inhibitory rebound burst that triggers a new pre-Bötzinger activation (Figures 8D and 9D).

thumbnail
Figure 8. Simulations of different synchronization regimes in healthy subjects (A–B) and COPD patients (C–D) are depending on the excitability level of the parafacial repiratory group (J1) and the pre-Bötzinger complex (J2).

Other fixed parameter values of the model are: ε  = 0.005, d = 0.4, β = 0.4, a = 0.2, m0 = 0.864, m1 = 0.65, δ = 0.2, xth = −0.02 (threshold for calcium current), τ1 = 10 and τ2 = 2 for the parafacial while τ1 = 5, τ2 = 10 for the pre-Bötzinger. In COPD patients, the parafacial respiratory group of the brainstem has a higher excitability level than healthy subjects and drives the pre-Bötzinger (active inspiration and expiration). (See results for comments).

https://doi.org/10.1371/journal.pone.0075740.g008

thumbnail
Figure 9. Chaotic attractor of the 2 synchronized pacemakers for respiratory rhythmogenesis in healthy subjects (A–B) and COPD patients (C–D) after simulations.

Each attractor is given according to the different network regime presented in Figure 8. The figure reveals that the coupling between both neuronal pacemaker exhibit nonlinear deterministic chaos (9B–C–D).

https://doi.org/10.1371/journal.pone.0075740.g009

Modeling fMRI signal based on simulated neural activity in healthy subjects and COPD patients.

To confirm our hypotheses on respiratory rhythmogenesis in healthy subjects and COPD patients, we performed 5 runs of simulations (250 action potentials with chaotic bursting oscillations for both pacemakers) of the synchronization regimes shown in Figure 8B–C. fMRI signals can then be modeled as a result of the convolution of the obtained neural states with a hemodynamic response function and added noise (see methods). The amplitude of the low frequency oscillations of the fMRI signal is then computed and shown in Figure 3 (bottom). The model is able to replicate the experimental fMRI results in both healthy subjects and COPD patients.

Discussion

We are the first, to our knowledge, to identify and describe the brainstem neural substrates underlying breathing complexity in healthy humans and patients with lung disease. fMRI scans revealed neural activity in the rostral ventro-lateral medulla and caudal ventro-lateral pons fitting the neural dynamics of respiratory rhythmogenesis. We then provided evidence that these central neural activities significantly correlate with the dynamical characteristics of the inspiratory and expiratory airflow in healthy humans and COPD patients (Table S4). Further, we developed a mathematical model of chaotic pacemakers where different neuronal excitabilities entrain different synchronization regimes and complexities that replicate key fMRI findings in humans.

Source of Human Ventilatory Complexity

We decided to focus on the core automatic network generating respiratory rhythmogenesis [4], [5], [9][12], [22], [23] since previous experimental and clinical works highlighted its potential contribution to airway flow complexity [14], [18], [24]. In the present study, ventilatory complexity significantly correlates with the activity of the respiratory central pattern generators assess with cerebral fMRI: in COPD patients, the increase in airflow complexity during expiration comes along with the higher VL pons parafacial activity while healthy subjects exhibit higher VL medulla activity with greater complexity during inspiration. Such parallel changes underline the contribution of the respiratory pacemaker neurons in airflow complexity. Previous works analyzed the mechanisms modulating chaos in airflow but failed to decipher the brainstem neural contribution to airflow complexity in human. It was previously shown that mechanical loading conditions alter chaos with an increase complexity in circumstances improving the load capacity-balance of the respiratory system [25], that breathing complexity was impaired during carotid stenosis due to the effects of autonomic baroreflex impairment on breathing control [26], and finally that chemoreceptor stimulation of ventilation by hypercapnia led to a high level of complexity [3]. Interestingly, while Fiamma et al. [3] showed in one study that hypercapnia stimulated ventilation and increased airway flow chaos, Pattinson et al. [27] demonstrated in a neuroimaging work that carbon dioxide stimulus activates brainstem respiratory centers of the ventral pons, rostral pons and lateral medulla. Some of these activated area overlapped with our regions of interest during the block design paradigm. Besides, we used a theoretical approach of respiratory rhythmogenesis to reproduce the core activity modes of neurons involved in the automatic respiratory network scheme in humans with two synchronized chaotic pacemakers, one driving inspiration, the pre-Bötzinger complex and the other driving expiration, the parafacial group. We chose to develop a map-based model [28], [29] for its relative simplicity compared with Hodgkin-Huxley formalism, and for its ability to generate spontaneous chaotic bursting activity. The model was further refined to incorporate post-inhibitory rebound bursting behavior. The mathematical model we propose is in line with previous experimental and theoretical works [13], [14]. In addition, it is able to exhibit chaotic behavior depending on the parameter value J which is the excitability level of the neuron. Above all, it reveals how a chaotic activity in neurons (Figure 9) contributes to chaos in airflow (Figure 2). Through controlling the excitability levels of the pre-Bötzinger and parafacial neurons in the mathematical model, different synchronizations and level of complexity appear. The choice of the parameter values, among them J1 and J2, are motivated by 2 characteristics: the ability to exhibit chaotic spike bursting oscillations (J value between Jmin and Jmax) and the specific synchronization regimes. Finally we verified our hypotheses on respiratory rhythmogenesis in healthy human and COPD patients (re-activation of the parafacial) with the mathematical model of the full handshake process and we were able to mimic the experimentally fMRI signals of the brainstem ventro-lateral medulla and ventro-lateral pons (Figure 3).

COPD Patients Breathe with a Higher Level of Complexity during Expiration than Controls

We found that patients with chronic obstructive pulmonary disease breathe with a higher level of complexity in airflow than healthy subjects. These unexpected findings cast doubt on the traditional view that complexity systematically degrades in disease state [30], [31]. Inspiratory and expiratory complexity changes parallel the activity of the VL medulla and VL pons, which contains the pre-Bötzinger and parafacial neurons respectively. It is therefore an in vivo estimate of the respiratory center function in humans as previously shown [18]. In healthy subjects, airflow complexity is higher during inspiration than expiration thus reflecting active inspiration while expiration is usually passive due to the elastic recoil of the lung. Conversely, patients with COPD have a higher level of complexity during expiration as compared with healthy subjects because they actively expire. In patients, fMRI revealed greater neuronal activity in the caudal VL pons region than in healthy subjects. Further studies are required to elucidate if patients having a high excitability of the caudal VL pons with the parafacial group, are those who effectively actively recruit their expiratory muscles, as suggested by Yan et al. [32]. We show that the excitability level of the neurons involved in respiratory rhythmogenesis in humans may vary depending on the physio-pathological conditions. These findings are in agreement with previous experiments in rats. In neonates, the parafacial expiratory group which has a high excitability level is dominant and drives the pre-Bötzinger [9], [10], [13], while in adults animals the parafacial is normally depressed and the pre-Bötzinger becomes dominant [5], [9], [10], [15]. Direct stimulation of parafacial neurons has been recently shown to promote active expiration in adult rats [33]. It is also possible to reactivate [34] the parafacial group during hypoxia [35]. Moreover, a previous study demonstrated that patients passively driven by a mechanical ventilator do not exhibit complexity in airflow whereas those with signs of active expiratory control displayed an increase complexity [18]. COPD patients have a forced expiratory flow limitation, which promotes the recruitment of abdominal muscles to sustain ventilation. The expiratory oscillator is probably turned on in patients to sustain ventilation in response to the increased respiratory load and hypoxia. Healthy subjects and COPD patients do differ in terms of PaO2. However, the contribution of O2 sensitive-chemoreceptors to the increase in airflow complexity in patients is weak since no difference between normoxic and hypoxic COPD patients is evidenced. Moreover, expiratory flow complexity differs between controls and normoxic COPD patients. From these results, we postulate that mechanical abnormalities due to disordered lung mechanics play a critical role in subsequent complexity alterations. Indeed, we found correlations between decrease pulmonary function and chaotic components in both univariate (Figure S3) and multivariate analyses (Figure 5). The increase in airflow complexity in patients is also related to systemic inflammation as shown during COPD [36]. A previous work in rats showed that brainstem cytokine level is high in a model of acute respiratory failure and this was strongly related to the increase in ventilatory complexity [24]. Finally, one additional explanation relies on the pathological narrowing of the bronchial tree and the direct “physical” consequences on the airflow: it is possible that some airflow turbulence due to local structural abnormalities and disordered lung mechanics directly contributes to increase airflow chaos, especially during expiration.

Interestingly, we could discriminate COPD patients with mild, moderate and severe dyspnea at rest according to expiratory flow complexity and the neural activity of the VL pons: patients with a severe dyspnea had a higher level of expiratory flow complexity and greater activity of the VL pons, as compared with patients having mild dyspnea. This difference was even less sensitive for the pulmonary function (Figure S2). Therefore, COPD patients having a severe dyspnea unexplained by a worsening of their pulmonary function, may exhibit an altered neuronal excitability of the VL pons, thereby reinforcing the central determinism of dyspnea.

Chaos in Airflow Decreases during Inspiratory Load, While Neural Activity of the Respiratory Centers Yields Negative BOLD Signals

Loading inspiration reduces airflow complexity with a parallel inhibition of the BOLD signal in the rostral medulla of healthy subjects. Our results differ from a previous study in 8 healthy subjects that did not find any effect of inspiratory loading on airflow chaos [37]. Differences in the experimental protocol may explain these discrepancies, i.e. the number of subjects included (25 healthy subjects in our study) and the duration of the load applied (15 minutes in our protocol). Furthermore, a previous work using fMRI found activation in the ventral pons of healthy subjects during inspiratory loading [38]. We point out that in the study of Gozal et al. [38] the protocol was different in terms of the load applied (30 cmH20/L/sec in their study), fMRI image acquisition and processing, specifically for the inclusion of confounding statistical regressors in the model. Moreover, negative BOLD signal changes were not specifically investigated [39]. Besides, it has been shown in 6 healthy subjects that voluntary hyperpnea targets the superior dorsal medulla of the brainstem [40]. In our study, the dorsal medulla showed significant de-activation during inspiratory resistive load. Differences in the stimulus applied (resistive load in our protocol) and in the characteristics of the healthy controls (16 controls in our study with older mean age 52±11) may explain these discrepancies. Of note, healthy subjects and COPD patients also exhibit activated brain regions known to be involved in the voluntary control of respiratory muscles, i.e. sensory-motor, premotor and supplementary motor cortex area (data not shown). The fact that the mechanical inspiratory load activates these cortical centers and de-activates in parallel the automatic network is physiologically relevant.

Loading inspiration in COPD patients leads to a diminution of airflow complexity of inspiration as well as expiration. These results are in line with the possible dual organization of respiratory rhythmogenesis in patients where reactivation of the parafacial occurred (Figure 8C–D). Once a stimulus is applied during inspiration it echoes on the other pacemaker due to the coupling characteristics. In addition, a diminution of fMRI BOLD responses in the two regions VL medulla and pons occurs in parallel in patients (Figure 7 and Table S3).

Study limitations.

A major challenge in application of fMRI to respiratory studies is the limited spatial and temporal resolutions of the BOLD signal [41], making it difficult to pinpoint precisely the specific brainstem respiratory related structures, which are generally rather small and heterogeneous with time-varying respiration related fluctuations. The pre-Bötzinger complex is a small structure and is bordered by other respiratory related nuclei including the Bötzinger complex. The parafacial respiratory group is a spread-out structure and contains both expiratory-related neurons and chemosensory neurons. We are however confident with our fMRI measurements for three reasons: (i) the first reason relies on the neuroanatomical paper recently published from Schwarzacher et al. [6]. The authors accurately identify in human brain autopsy the location of the pre-Bötzinger complex. The diameter of the complex is around 5–6 mm, in the ventro-lateral region of the rostral medulla, 9 mm from obex, below Fissura Pontomedullaris. For all participants of our fMRI protocol, we individually computed these coordinates in standard images. Then the regions of interest were centered on these coordinates and transformed from standard space to functional space for the extraction of the time series. The parafacial respiratory group is located near the pre-Bötzinger in the caudal ventro-lateral pons, ventro-laterally to the facial nerve nucleus VII [6], [7], above Fissura Pontomedullaris. (ii) The second reason relies on the de-activation regions evidenced during the block design paradigm with inspiratory resistive load. Theses inhibited regions overlapped the coordinates defined in the rest fMRI acquisition and we also found strong correlation between the mean negative BOLD signal and the chaotic component using the same coordinates than rest fMRI acquisition. (iii) The third reason concerns the theoretical part of the work. We modeled respiratory rhythmogenesis with two pacemakers that synchronously handshake one another, depending on their excitability level [13]. The resulting neural time series of the pre-Bötzinger and parafacial groups, convolved with a hemodynamic function plus noise replicate experimental fMRI signal in healthy subjects and COPD patients.

Furthermore, we cannot exclude the potential influence of emotion via the limbic system on the automatic network [42]. However before airflow recordings begin, the subjects were allowed to adapt for 5 minutes to the materials and were quiet. We also removed the first 2 minutes of recordings for subsequent analyses. Additionally, we took time to explain the fMRI protocol to both healthy subjects and COPD patients. For fMRI protocol, the participants were instructed to ‘keep their eyes closed and think of nothing in particular’. They were instructed to refrain from cognitive, language, and motor tasks. The participants knew that a physician was near the scanning room and they all had the possibility to stop the images acquisition if a problem arised. We therefore minimized as much as possible the possible influence of emotions on our experiments.

Perspectives.

In this study, we decipher the brainstem neural substrates of airflow complexity in humans. We also shed new lights on the brainstem neural control of respiratory muscles in patients with COPD. The patients have an increased complexity of the airflow during expiration that correlates with the high activity of VL pons. COPD patients reactivate the parafacial neuronal group, as shown with the mathematical model and fMRI results, to sustain ventilation. These findings may be involved in the onset of respiratory failure when the neural network becomes overwhelmed by respiratory overload as suggested by previous works [43], [44]. Future works analyzing the relationships between automatic and cortical network from a theoretical and experimental viewpoint will help to clarify the mechanisms preceding acute respiratory failure. Moreover, we show that COPD patients having a severe dyspnea unexplained by a worsening of their pulmonary function, may exhibit an altered neuronal excitability of the VL pons, thereby reinforcing the central determinism of dyspnea. Identifying the activity of the respiratory pacemakers through both airflow complexity and functional imaging techniques opens new strategies to refine COPD patient phenotypes.

Methods

Participants and Protocol

Stable patients (n = 25) with COPD (no exacerbation for 4 weeks) were recruited from the Physiology and Respiratory disease departments of the Bichat University Hospital 2011–2012. Inclusion criteria were patients above eighteen having mild to severe COPD according to clinical and pulmonary function test criteria [45]. Exclusion criteria were home oxygen therapy, neurological disease, past history of stroke, psychiatric disorder, body mass index above 30 kg/m2, contraindication to cerebral functional magnetic resonance imaging. After given written informed consent, patients had a clinical examination and pulmonary function tests. In COPD patients, dyspnea was quantified at rest using Borg scale. An age-matched control group (n = 25) was recruited from the Centre d’ Investigation Clinique of the Bichat Hospital. The protocol was approved by the ethics committee Ile-de-France 1.

Subjects were comfortably seated and were asked to keep their eyes open. They wore a nose clip and breathed through a mouthpiece connected to a low resistance pneumotachograph (MLT1000L-AD Intruments) via a two-way non-rebreathing valve (Hans Rudolph 1410 Series). Ventilatory flow, digitized at 400-Hz sampling rate was recorded on a PC computer in the form of data files for subsequent analysis (Chart5, AD Instruments). Mouth pressure was measured at the mouthpiece and connected to a pressure transducer (MLT0699-AD Instruments). Ventilatory flow and mouth pressure were synchronously recorded on the PC computer via the PowerLab 4/25 (AD Instruments). End-tidal PCO2 (PETCO2), measured from a side port of the mouthpiece and finger oxygen saturation were connected to a portable Oxi-capnography (MD-660P Comdek) for continuous acquisition. Before recordings began, the subjects were allowed to adapt for 5 minutes to the materials and were quiet. Recordings were performed during 15 minutes at the same time of the day for all subjects. Two sets of measurements were performed in random order, one with subjects breathing spontaneously and one with subjects breathing during the continuous application of an inspiratory resistive load of 20 cmH20/L/sec (7100R20 Hans Rudolph). Reproducibility of our measurements was previously tested [26]. Ventilatory flow recordings will be available upon request to the corresponding author.

Linear and Nonlinear Analyses of Airflow

The first two minutes of recording were excluded from the analyses. Inspiratory (Vt/Ti) and expiratory (Vt/Te) flows were computed on a breath-by-breath basis during spontaneous breathing and during the inspiratory effort, i.e. during continuous application of the resistive load on the inspiratory phase of the respiratory cycle.

Analysis of Ventilation in the Time Domain and Autocorrelation Analyses

Fluctuations of the inspiratory and expiratory flows were first evaluated through their coefficients of variation (the ratio of the standard deviation to the mean). Autocorrelation of the flows was computed at a lag of one breath. It estimated the amount of correlated linear part of the flow [18], [26], [46].

Nonlinear Analyses

Chaos detection.

The noise titration technique [47] was used on the inspiratory and expiratory flow time series. This method has already been proven its accuracy to evidence the chaotic nature of human ventilation [3], [18], [26], [46]. It involved the simulation of time series with linear and nonlinear polynomial autoregressive model (Volterra-Wiener series) [48]. The best linear and nonlinear models are chosen according to the minimal information theoretic criterion. The null hypothesis, a stochastic time series with linear dynamics, is rejected if the best nonlinear model provided a significant better fit to the data than the best linear model using parametric (F-test) statistics at the 1% significance level. Once nonlinear determinism is indicated, white noise of increasing standard deviation is added to the data until nonlinearity can no longer be detected, i.e. the nonlinearity is ‘neutralized’. The noise limit (NL) is calculated as the percent of signal power added as noise to ‘titrate’ the data to the point of neutrality. Typically, an average NL value is obtained by repeating the titration procedure 5 times. Under this scheme, chaos is indicated by NL>0, and the value of NL provides a relative measure of chaos intensity. Conversely, if NL = 0, then it may be inferred that the series either is not chaotic or the chaotic component is already neutralized by the background noise (noise floor) in the data. We then estimated the largest Lyapunov exponent and the correlation dimension of the time series having a positive noise limit value.

Sensitivity to initial conditions.

Complex dynamical systems are sensitive to initial conditions, and exhibit an exponential divergence in the phase space. This can be quantified through the study of the Lyapunov exponents spectrum and the calculation of the largest Lyapunov exponent (λL: LLE). Consider two points on two nearby trajectories in the phase space, and assume the distance between them to be d(0). After time t, if the distance between the two trajectories becomes d(t), then the average divergence (separation after time t) can be written as:where λL is the LLE of the system. In the present study, we used the algorithm proposed by Rosenstein et al. that has been shown to be particularly useful for small data series [49].

Irregularity.

The correlation dimension is a fractal dimension reflecting the irregularity of the attractor of the system. It characterizes the “aperiodicity” of the system in the phase space. It is estimating by examining the scaling properties of the correlation sum [49]. From a time series , where N is the total number of points, the m dimensional vector in the phase space can be constructed by delay embedding:where, τ is the fixed time lag and m is the embedding dimension. Then the reconstructed trajectory of the actual dynamics can be written as where

The correlation dimension can be calculated from the correlation integral of the time series. The correlation integral can be computed as follows [49], [50]:where, r is scale length, and θ is the Heaviside step function. Scaling of the function C(r,m) can be written as:

The correlation dimension (Dcorr) can be defined byand for practical purpose, Dcorr can be obtained from the slope of ln C(r) vs ln r plot.

Time lag was first estimated by a drop of the autocorrelation to [49][51]. The optimal dimension was obtained after calculating the percentage of false nearest neighbors between points in phase space. A minimal number of false nearest neighbors was required [52]. The embedding dimension that adequately represents the system is the dimension that eliminates most of the false nearest neighbors allowing an adequate phase-space reconstruction of the underlying dynamic. An appropriate time lag and embedding dimension were estimated for each experimental time series.

Surrogate data.

In order to test the nonlinearity that governs the dynamics, we have applied surrogate test [53]. First the Fourier transform of the original time series is computed. The phase is replaced by random numbers and finally the inverse Fourier transform is applied. Power spectrum is thus preserved although the nonlinear structures are destroyed [51], [53]. Correlation dimension was estimated for both the original data and five surrogates that match each original signal. A global test was carried out by a Wilcoxon signed-rank test comparing the correlation dimension values computed on the original data paired with the corresponding average correlation dimension values form the matching surrogate. Significant Wilcoxon rank test between the original and surrogates implies the nonlinear dynamics of the original data [18], [26], [46], [53].

Cerebral Functional Magnetic Resonance Imaging

Protocol and image acquisition.

Participants were imaged while lying comfortably in the scanner. Three sets of images were performed: structural, resting state and block design paradigm. For the structural and functional resting state, the participants breathed spontaneously while during the block design paradigm, they breathed via a mouthpiece connected a two-way non-rebreathing valve (Hans Rudolph 1410 Series) with nose clip. A small plastic tube of one meter length was connected to the inspiratory limb of the T-valve for application of the resistive load (20 cm/L/sec).

Physiological monitoring synchronized with the images acquisition was performed for the resting state and block design paradigm. Chest expansion was measured with a pneumatic belt and electrocardiogram was acquired with chest electrodes [54]. Sampling rates were 10 ms and 1 ms respectively. Respiratory volume per time (RVT) was computed from the respiratory waveform (chest belt) [55]. Maximum minus minimum of the waveform was divided by the breathing period for each breath cycle and then interpolated to the imaging repeat time (RT). PETCO2 and saturation were also continuously recorded with 10 ms sampling rates. The RR cardiac interval, PETCO2 and saturation (maximum values per breath), were also interpolated to the imaging RT.

Imaging was performed using a 3 Tesla MR scanner (General Electrics, USA) with a 64-channel head coil. T1-weight high resolution 3D volume covering the entire brain was acquired in controls (n = 16) and COPD patients (n = 17). Acquisition parameters were: 171 axial slices, 1.2 mm thickness with no gap echo time [Te] = 3.4 ms, repeat time [TR] = 8.6 ms, flip angle = 12°, matrix 256×256, field of view 240 mm×240 mm). The total acquisition time was 4 min 35 s.

T2-weighted echoplanar images were acquired for the resting state functional acquisition (52 axial slices, 4 mm thickness with no gap echo time [Te] = 19 ms, repeat time [TR] = 2000 ms, flip angle = 90°, matrix 64×64, field of view 240 mm×240 mm, and voxel dimension 3 mm3). Acquisition time was 10 min08 s, yielding 300 whole brain volume. For the resting state, the participants were instructed to ‘keep their eyes closed and think of nothing in particular’. They were instructed to refrain from cognitive, language, and motor tasks as much as possible, but not to fall asleep. Resting state fMRI scans will be available upon request to the corresponding author.

The second set of functional image was performed during a block design, which consists in 5 cycles of rest periods (36 sec), in alternate with active period (36 sec) during which the inspiratory resistive load (20 cmH20/L/sec) was applied on the breathing circuit. MRI parameters were: 52 axial slices, 4 mm thickness with no gap echo time [Te] = 33 ms, repeat time [TR] = 3000 ms, flip angle = 90°, matrix 64×64, field of view 240 mm×240 mm, and voxel dimension 3 mm3. Total acquisition time was 6 min12 s, yielding 120 whole brain volumes.

Image analyses.

Image processing was performed using FSL software (http://www.fmrib.ox.ac.uk/fsl, Oxford University).

Resting state fMRI.

Preprocessing steps included motion correction using MCFLIRT [56] slice timing corrections, non-brain removal using BET [57], spatial smoothing using a Gaussian kernel of full-width-half-maximum 6 mm, multiplicative mean intensity normalization of the volume at each time point. A brain mask was constructed with four regions of interest (cubes radii 6 mm) individually positioned on standard images over the brainstem in regions known to cover the respiratory generator nuclei in rostral ventro-lateral medulla oblongata and caudal ventro-lateral pons according to Schwarzacher et al. (Figure 4). These regions of interests were then transformed from standard space to functional space and the mean BOLD signal time series were then extracted. For all participants, the respiratory volume per time (RVT), the RR cardiac interval, PETCO2 [58], [59] and saturation were included in a multivariate regression linear model to account for significant influences on the BOLD signal. These covariates were then regress out.

Low frequency oscillations have been used in resting state fMRI in physiology and pathology to analyze the functional connectivity among brain regions [60][62]. Amplitude of the low frequency oscillations (AlFO) of the resulting BOLD signal time series is also a mean to assess neuronal activation with fMRI [63], [64]. BOLD time series were detrended and filtered between 0.01 and 0.08 Hz to remove the effects of very low-frequency drift and high-frequency noise. Fast Fourier transform (FFT) was applied and the power spectrum obtained. The average square root of the power spectral density was calculated across 0.01–0.08 Hz and this represents the AlFO. For normalization purposes, the AlFO of each regions of interest was divided by the global mean AlFO value of the whole brainstem. The standardized AlFO have a value about 1 and this procedure is analogous to that used in PET studies [65]. Finally the mean of the normalized AlFO of the 2 cubes of the medulla and the 2 cubes of the pons were averaged.

Block design fMRI.

Preprocessing step were the same as resting state fMRI with an additional high pass temporal filtering (Gaussian-weighted least-squares straight line fitting with sigma = 36 s). At the single level analysis we used a general linear model. Confounding regressors that potentially altered cerebral blood flow (RVT, PETCO2, RR cardiac interval, saturation) were included. Voxel-wise statistical analysis was extended to a second (group) level in a fixed-effects analysis. After analysis, statistical images were registered to high resolution structural and/or standard space images using FLIRT [66]. Registration from high resolution structural to standard was then further refined using FNIRT nonlinear registration.

Statistical analyses of airflow dynamics and fMRI data.

Matlab R2011a was used for statistical and signal processing analyses (Mathworks USA). Comparisons between clinical data among the groups were made using univariate analysis and χ2 test. The normality of the distributions of the discrete respiratory variables was ascertained using the Kolmogorov-Smirnov test. The occurrence of a positive noise limit in the airflow time series was compared using the χ2 test. Paired and unpaired T tests were used to study statistical differences of the linear and nonlinear measures of the inspiratory and expiratory flows among the groups. Pearson’s correlation coefficient was estimate for identifying significant relationships between airflow complexity and AlFO, FEV1/FVC, PaO2, PaCO2. Among the variables that had significant correlation with airflow complexity in univariate analysis, we then performed a multiple linear regression analysis to study the strength of its relation. During inspiratory load, the correlations were established in the whole population between airflow complexity and the mean negative BOLD signal.

Mathematical Model of Respiratory Rhythmogenesis

Two pacemaker-like neurons have been identified in mammals in the ventro-lateral column of the brainstem, the pre-Bötzinger complex inspiratory group and parafacial expiratory group respectively [5][13]. Previous works showed that the parafacial group exhibits pre-inspiratory activity [9], [35] as well as a rebound bursting after inspiration [35] while the dynamics of both pacemakers display chaotic spike-bursting oscillations [14]. We therefore chose to develop a map-based model for respiratory rhythmogenesis for its relative simplicity compared with Hodgkin-Huxley formalism, and for its ability to generate spontaneous chaotic bursting activity [28], [29]. The model is developed based on the discrete version of FitzHugh-Nagumo model by adding Heaviside step function H(x). Each pacemaker is modeled by the two dimensional original Courbage-Nekorkin map [28], [29] which is further refined to incorporate post-inhibitory rebound bursting behavior:(1)Where x qualitatively defines the dynamics of the membrane potential of the neuron and y is the common variable specifying the dynamics of all outwards ionic currents (recovery variable). β and d controls the threshold properties of the oscillation, ε is a positive parameter setting the time scale of the recovery variable y. J is associated with excitability properties of the neuron; F(x) is a piece-wise linear version of the cubic function in the FitzHugh-Nagumo model:

with , , and mo, m1>0

IT is a low-threshold calcium Ca2+ current [67] defines as:(2)Where k in equation (2) is a slow variable representing the inactivation of the low-threshold calcium conductance, which involves T-type Ca2+ calcium channels and produces a transmembrane current IT. δ represents the maximum conductance associated with IT. G(k) represents the dynamics of IT as follow:

(3)In this form the model is capable of post-inhibitory rebound bursting when xth is below the resting values of x. In equation (3), τ1 sets the duration of the burst and τ2 sets the duration of the hyperpolarization necessary to recruit a maximal post-inhibitory rebound response.

In equation (1), Isyn is the chemical synaptic coupling between the pre-Bötzinger complex and the parafacial group in the following form:Where K is the coupling strength which value is positive for excitatory synapse and negative for inhibitory synapse and rect is the rectangle function as described below:

Where ni is the step of the ith spike in the presynaptic neuron and τ is the duration of the postsynaptic current. A post-inspiratory inhibitory feedback is introduced from the pFRG with the same amplitude and duration of the rebound bursts for “inspiratory off-switch” to prevent the preBötC from reactivation.

Modeling fMRI signal based on simulated neural activity.

To confirm our hypotheses on different synchronization regimes of respiratory rhythmogenesis in healthy subjects and patients with respiratory failure, we performed 5 runs of simulations and then convolved the simulated neural states of the pre-Bötzinger complex and parafacial group with a hemodynamic response function. We used Statistical Parametrics Mapping software for the hemodynamic convolution: http://www.fil.ion.ucl.ac.uk/spm.

Under linear assumption, fMRI signals m(t) can then be modeled as a result of the convolution of neural states s(t) with a hemodynamic response function h(t), ε(t)is the noise.Where t is the time and ∶ denotes convolution, h(t) is the hemodynamic response function which is a mixture of two gamma functions. The parameter values of the hemodynamic response function are: delay of response relative to onset : 6 (s), delay of undershoot relative to onset = 16 (s), dispersion of response = 1 (s), dispersion of undershoot = 1 (s), ratio of response to undershoot 6 (s), onset = 0, length of kernel = 32 (s); ε(t) is the noise in the measurement assumed to be Gaussian white noise with mean zero and standard deviation 0.25. This value was chosen equal to the standard deviation of the mean BOLD time series. We model fMRI signal for 2 network schemes shown in figure 8 (B,C) and compute the AlFO of the modeled fMRI signal.

Supporting Information

Figure S1.

End-tidal PCO2 measurements during unloaded and inspiratory resistive load (ventilatory flow measurements) as well as during fMRI acquisition. Results are given for the 25 healthy subjects (A) and 25 COPD patients (B). C: End-tidal PCO2 measurements during resting state fMRI acquisition in healthy subjects (blue) and COPD patients (red). The means and standard deviations of the healthy subjects (n = 16) and the COPD patients (n = 17) are shown.

https://doi.org/10.1371/journal.pone.0075740.s001

(ZIP)

Figure S2.

Comparisons between COPD patients having mild, moderate and severe dyspnea (Borg scale) at rest according to expiratory flow complexity (A), the amplitude of the low frequency oscillations (AlFO) of the ventro-lateral (VL) pons (B), and the pulmonary function index (FEV1/FVC) (C). The patients with a severe dyspnea have a higher level of expiratory flow complexity and greater activity of the VL pons, as compared with patients having mild dyspnea. This difference is even less sensitive for the pulmonary function.

https://doi.org/10.1371/journal.pone.0075740.s002

(ZIP)

Figure S3.

Linear correlation between expiratory flow complexity (top: Noise limit, bottom: Largest Lyapunov exponent) and pulmonary function index (FEV1/FVC) in the whole population of healthy subjects and COPD patients. COPD patients are classified according to the diminution of their pulmonary function (GOLD classification).

https://doi.org/10.1371/journal.pone.0075740.s003

(ZIP)

Figure S4.

Negative BOLD signal of the cerebral fMRI during inspiratory resistive loading in healthy subjects (left) and COPD patients (right). Group analyses of the block design are given for the healthy subjects (left, n = 16) and COPD patients (right, n = 17). Sagittal and axial slices are shown on the top panel. Bottom: The corresponding mean time series of the ventro-lateral medulla of the 16 healthy subjects and 17 COPD patients are shown. The figures show the diminution of the BOLD signal during each application of the resistive load (5 cycles of rest (R: black line) and active task (A: red line) with resistive load).

https://doi.org/10.1371/journal.pone.0075740.s004

(ZIP)

Table S1.

Clinical characteristics of the 25 COPD patients.

https://doi.org/10.1371/journal.pone.0075740.s005

(DOCX)

Table S2.

Linear measures of the ventilatory variables (mean values, coefficients of variation and autocorrelation) during unloaded breathing and during inspiratory resistive load (20 cmH20/L/sec).

https://doi.org/10.1371/journal.pone.0075740.s006

(DOCX)

Table S3.

Number of time series that exhibit positive noise limit value for chaos characterization in the inspiratory and expiratory flow time series in controls and COPD patients.

https://doi.org/10.1371/journal.pone.0075740.s007

(DOCX)

Table S4.

Summary of the results concerning airflow complexity and brainstem respiratory centers activity in healthy subjects and patients with chronic obstructive pulmonary disease (COPD).

https://doi.org/10.1371/journal.pone.0075740.s008

(DOCX)

Acknowledgments

We thank the Unité de Recherche Clinique Paris Nord of the Bichat Hospital and Miss Naïma Beldjoudi for technical assistance. We also thank the staff of the physiology department for active participation in the study and the Critical Care department of the Bichat Hospital for the help in equipment. We thank Jean Champagnat for helpful discussion and Pr Poon from the MIT for providing the noise titration code. We are also grateful to the anonymous reviewers for their comments and analyses on the manuscript.

Author Contributions

Conceived and designed the experiments: LM. Performed the experiments: AH LY IK GJ HM OB MF GD ESC LM. Analyzed the data: AH LY MC MM CC LM. Wrote the paper: MC CC LM.

References

  1. 1. Chiaro DR (2010) Emergent complex neural dynamics Nature Physics. 6: 744–750.
  2. 2. Courbage M, Kazantsev VB, Nekorkin VI, Senneret M (2004) Emergence of chaotic attractor and anti-synchronization for two coupled monostable neurons. Chaos 14: 1148–1156.
  3. 3. Fiamma MN, Sraus C, Thibault S, Wysocki M, Baconnier P, et al. (2006) Effects of hypercapnia and hypocapnia on ventilatory variability and the chaotic dynamics of the ventilatory flow in humans. Am J Physiol 292: R1985–R1993.
  4. 4. Smith JC, Ellenberger HH, Ballanyi K, Richter DW, Feldman JL (1991) Pre-Bötzinger complex: a brainstem region that may generate respiratory rhythm in mammals. Science 254: 726–728.
  5. 5. Feldman JL, Del Negro CA (2006) Looking for inspiration: new perspectives on respiratory rhythm. Nature Rev Neurosc 7: 232–241.
  6. 6. Schwarzacher SW, Rüb U, Deller T (2011) Neuroanatomical characteristics of the human pre-Bötzinger complex and its involvement in neurodegenerative brainstem diseases. Brain 134: 24–35.
  7. 7. Lavezzi A, Matturri L (2008) Functional neuroanatomy of the human pre-Bötzinger complex with particular reference to sudden unexplained perinatal and infant death. Neuropathol 28: 10–16.
  8. 8. Gray P, Janczewski W, Mellen N, McCrimmon D, Feldman JL (2011) Normal breathing requires preBötzinger complex neurokinin-1 receptor-expressing neurons. Nature Neurosc 4: 927–930.
  9. 9. Onimaru H, Homma I (2003) A novel functional neuron group for respiratory rhythm generation in the ventral medulla. J Neurosci 23: 1478–1486.
  10. 10. Janczewski WA, Feldman JL (2006) Distinct rhythm generators for inspiration and expiration in the juvenile rat. J Physiol 57: 407–420.
  11. 11. Dubreuil V, Ramanantsoa N, Trochet D, Vaubourg V, Amiel J, et al. (2008) A human mutation in Phox2b causes lack of CO2 chemosensitivity, fatal central apnea, and specific loss of parafacial neurons. Proc. Natl. Acad. Sci. USA 105: 1067–1072.
  12. 12. Thoby-Brisson M, Karlen M, Charnay P, Champagnat J, Fortin G (2009) Genetic identification of an embryonic parafacial oscillator coupling to the preBötzinger complex Nature Neurosc. 12: 1028–1036.
  13. 13. Wittmeier S, Song G, Duffin J, Poon CS (2008) Pacemakers handshake synchronization mechanism of mammalian respiratory rhythmogenesis. Proc Natl Acad Sci USA 105, 18000–18008.
  14. 14. Del Negro CA, Wilson CG, Butera RJ, Rigatto H, Smith JC (2002) Periodicity, mixed mode oscillations, and quasi periodicity in a rhythm-generating neural network. Biophys J 82: 206–214.
  15. 15. Caubit X, Thoby-Brisson M, Voituron N, Filippi P, Bévengut M, et al. (2010) Teashirt 3 regulates development of neurons involved in both respiratory rhythm and airflow control. J Neurosci 30: 9465–9465.
  16. 16. Que CL, Kenyon CM, Olivenstein R, Macklem PT, Maksym GN (2001) Homeokinesis and short term variability of human airway caliber. J Appl Physiol 91: 1131–1141.
  17. 17. D’Angelo E, Robatto FM, Calderini E, Tavola M, Bono D, et al. (1991) Pulmonary and chest wall mechanics in anesthetized paralyzed humans. J Appl Physiol 70: 2602–2610.
  18. 18. Mangin L, Fiamma MN, Straus C, Derenne JP, Zelter M, et al. (2008) Source of human ventilatory chaos: lessons from switching controlled mechanical ventilation to inspiratory pressure support in critically ill patients. Resp Physiol Neurobiol 161: 189–196.
  19. 19. Logothetis NK, Pauls J, Augath M, Trinath T, Oeltermann A (2001) Neurophysiological investigation of the basis of the fMRI signal. Nature 412: 150–157.
  20. 20. Jolley CJ, Luo YM, Steier J, Reilly C, Seymour J, et al. (2009) Neural respiratory drive in healthy subjects and in COPD. Eur Resp J 33: 289–297.
  21. 21. Abraham KA, Feingold H, Fuller D, Jenkins M, Mateika J, et al. (2002) Respiratory-related activation of human abdominal muscles during exercice. J Physiol 541: 653–663.
  22. 22. Onimaru H, Homma I (2006) Point:Counterpoint: The parafacial respiratory group (pFRG)/pre-Bötzinger complex (preBötC) is the primary site of respiratory rhythm generation in the mammal. J Appl Physiol 100: 2094–2098.
  23. 23. Feldman JL, Del Negro C, Gray PA (2013) Understanding the rhythm of breathing: so near, yet so far. Ann Rev Physiol 5: 423–452.
  24. 24. Jacono PJ, Mayer CA, Hirsch YH, Wilson CG, Dick TE (2010) Lung brainstem cytokine levels are associated with breathing pattern changes in a rodent model of acute lung injury. Resp Physiol Neurobiol 178: 429–438.
  25. 25. Schmidt M, Demoule A, Cracco C, Gharbi A, Fiamma MN, et al. (2010) Neurally adjusted ventilatory assist increases respiratory variability and complexity in acute respiratory failure. Anesthesiology 112: 670–681.
  26. 26. Mangin L, Leseche G, Duprey A, Clerici C (2011) Ventilatory chaos is impaired in carotid atherosclerosis. PLoS ONE 6: e16297.
  27. 27. Pattinson KT, Mitsis GD, Harvey AK, Jbabdi S, Dirckx S, et al. (2009) Determination of the human brainstem respiratory control network and its cortical connections in vivo using functional and structural imaging. Neuroimage 44: 295–305.
  28. 28. Courbage M, Nekorkin VI, Vdovin LV (2007) Chaotic oscillations in a map-based model of neural activity. Chaos 17: 043109.
  29. 29. Courbage M, Nekorkin VI (2010) Map based models in neurodynamics. Int J Bif Chaos 20: 1631–1651.
  30. 30. Poon CS, Merrill CK (1997) Decrease of cardiac chaos in congestive heart failure. Nature 389: 492–495.
  31. 31. Goldberger AL, Amaral LA, Hausdorff J, Ivanov PC, Peng CK, et al. (2002) Fractal dynamics in physiology: Alterations with disease and aging. Proc Nat Acad Sci USA 99: 2466–2472.
  32. 32. Yan S, Sinderby C, Bielen P, Beck J, Comtois N, et al. (2000) Expiratory muscles pressure and breathing mechanics in chronic obstructive pulmonary disease. Eur Resp J 16: 684–690.
  33. 33. Pagliardini S, Janczewski W, Tan W, Dickson CT, Deisseroth K, et al. (2011) Active expiration induced by excitation of ventral medulla in adult anesthetized rats. J Neurosci 31: 2895–2905.
  34. 34. Milsom WK (2010) Adaptive trends in respiratory control: a comparative perspective. Am J Physiol 299: R1–R10.
  35. 35. Fortuna MG, West GH, Stornetta RL, Guyenet PG (2008) Bötzinger expiratory-augmenting neurons and the parafcial respiratory group. J Neurosc 28: 2506–2515.
  36. 36. Barnes PJ (2010) Chonic obstructive pulmonary disease: effects beyond the lungs. PLoS Med 7: e1000220.
  37. 37. Samara Z, Raux M, Fiamma MN, Gharbi A, Gottfried SB, et al. (2009) Effects of inspiratory loading on the chaotic dynamics of ventilatory flow in humans. Resp Physiol Neurobiol 165: 82–89.
  38. 38. Gozal D, Omidvar O, Kirlew KA, Hathout G, Hamilton R, et al. (1995) Identification of the human brain regions underlying responses to resistive inspiratory loading with functional magnetic resonance imaging. Proc Natl Acad Sci USA 92: 6607–6611.
  39. 39. Shmuel A, Augath M, Oeltermann A, Logothetis NK (2006) Negative functional MRI response correlates with decreases in neuronal activity in monkey visual area V1. Nature Neurosc 9: 569–577.
  40. 40. McKay LC, Evans KC, Frackowiak RS, Corfield DR (2003) Neural correlates of voluntary breathing in humans. J Appl Physiol 95: 1170–1178.
  41. 41. Logothetis NK (2008) What we can do and what we canot do with fMRI. Nature 453: 869–878.
  42. 42. Evans KC, Dougherty DD, Schmid AM, Scannell E, McCallister A, et al. (2009) Modulation of spontaneous breathing via limbic/paralimbic-bulbar circuitry: an event-related fMRI study. Neuroimage 47: 961–971.
  43. 43. Hopkinson NS, Sharshar T, Ross ET, Nickol AH, Dayer MJ, et al. (2011) Corticospinalcontrol of respiratory muscles in chronic obstructive pulmonary disease. Resp Physiol Neurobiol 141: 1–12.
  44. 44. Murphy PB, Kumar A, Reilly C, Jolley C, Walterspacher S, et al. (2011) Neural respiratory drive as a physiological biomarker to monitor change during acute exacerbations of COPD. Thorax 66: 602–608.
  45. 45. Pauwels RA, Buist AS, Calverley PM, Jenkis CR (2001) Hurd S. on behalf of the GOLD scientific committee (2001) Global strategy for the diagnosis, management and prevention of chronic obstructive pulmonary disease. Am J Resp Crit Care Med 163: 1256–1276.
  46. 46. Mangin L, Clerici C, Similowski T, Poon CS (2009) Chaotic dynamics of cardioventilatory coupling in humans: Effects of ventilatory modes. Am J Physiol 296: R1088–1097.
  47. 47. Poon CS, Barahona M. Titration of chaos with added noise (2001) Proc Natl Acad Sci USA. 98: 7001–7112.
  48. 48. Barahona M, Poon CS (1996) Detection of nonlinear dynamics in short, noisy time series. Nature 381: 215–217.
  49. 49. Rosenstein M, Collins J, De Luca C (1993) A practical method for calculating largest Lyapunov exponents for small data sets. Physica D 65: 117–134.
  50. 50. Grassberger P, Procacia I (1983) Characterization of strange attractors. Phys Rev Lett 50: 346–349.
  51. 51. Kantz H, Schreiber S (2004) Nonlinear Time series Analysis, 2nd ed. Cambridge University Press. Cambridge UK.
  52. 52. Lieber W, Pawelzik K, Schuster HG (1991) Optimal embeddings of chaotic attractor from topological considerations. Europhys Lett 14: 521–526.
  53. 53. Theiler J, Eubank S, Longtin A, Galdrikian B, Farmer JD (1992) Testing for nonlinearity in time series: the method of surrogate data. Physica D 58: 77–94.
  54. 54. Glover GH, Chang C (2009) Relationship between respiration, end-tidal CO2, and BOLD signals in resting-state fMRI. Neuroimage 47: 1381–1393 (2009)..
  55. 55. Birn R, Smith M, Jones T, Bandetti P (2008) The respiration response function: the temporal dynamics of fMRI signal fluctuations related to changes in respiration. Neuroimage 40: 644–654.
  56. 56. Jenkinson M, Banniste P, Brady M, Smith S (2002) Improved optimization for the robust and accurate linear registration and motion correction of brain images. Neuroimage 17: 825–841.
  57. 57. Smith S (2002) Fast Robust Automated Brain Extraction. Hum Brain Mapp 17: 143–155.
  58. 58. Wise RG, Ide K, Poulin MJ, Tracey I (2004) Resting fluctuations in arterial carbon dioxide induce significant low frequency variations in BOLD signal. Neuroimage 21: 1652–1664.
  59. 59. Peng T, Niazy R, Payne SJ, Wise RG (2013) The effects of respiratory CO2 fluctuations in the resting-state BOLD signal differ between eyes open and eyes closed. Magn Res Imaging in press.
  60. 60. Greicius M, Kranow B, Reiss A, Menon V (2003) Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. Proc Natl Acad Sci USA 100: 253–258.
  61. 61. Greicius M, Srivastava G, Reiss A, Menon V (2004) Default-mode network activity distinguishes Alzheimer’s disease from healthy aging: Evidence from functional MRI. Proc Natl Acad Sci USA 101: 4637–4642.
  62. 62. Fox M, Raichle M (2007) Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nature Rev Neurosc 8: 700–711.
  63. 63. Zhang Z Q, Lu G, Zhong Y, Tan Q, Chen H, et al. (2010) fMRI study of mesial temporal lobe epilepsy using amplitude of low frequency fluctuation analysis. Hum Brain Mapp 31: 1851–1861.
  64. 64. Jiang GH, Qiu YW, Zhang XL, Han LJ. Lv XF (2011) Amplitude low-frequency oscillations abnormalities in the heroin users: A resting state fMRI study. Neuroimage 57: 149–154.
  65. 65. Raichle ME, MacLeod AM, Snyder AZ, Powers WJ, Gusnard DA, et al. (2001) A default mode of brain function. Proc Natl Acad Sci USA 98: 676–682.
  66. 66. Jenkinson M, Smith S (2001) A global optimization method for robust affine registration of brain images. Med Image Anal 5: 143–156.
  67. 67. Smith GD, Cox CL, Sherman M, Rinzel J (2000) Fourier analysis of sinusoidally driven thlamocortical relay neurons and a minimal integrate-and-fire-or-burt model. J Neurophysiol 83: 588–610.