Sources of variability in expiratory flow profiles during sleep in healthy young children

Standard lung function tests are not feasible in young children, but recent studies show that the variability of expiratory tidal breathing flow-volume (TBFV) curves during sleep is a potential indirect marker of lower airway obstruction. However, the neurophysiological sources of the TBFV variability in normal subjects has not been established. We investigated sleep stages and body position changes as potential sources for the TBFV curve variability. Simultaneous impedance pneumography (IP), polysomnography (PSG) and video recordings were done in 20 children aged 1-7 years without significant respiratory disorders during sleep. The early part of expiratory TBFV curves are less variable between cycles of REM than NREM sleep. However, within individual sleep cycles, TBFV curves during N3 are the least variable. The differences in TBFV curve shapes between sleep stages are the main source of overnight variability in TBFV curves and the changes in body position have a lesser impact.


Introduction
Objective evidence of airway obstruction is essential for the diagnosis and monitoring of asthma (Global Initiative for Asthma, 2018). Usually lung function is tested by spirometry, but for subjects with limited cognitive or physical capabilities, such as patients with developmental disabilities, elderly or young children, the spirometric respiratory maneuvers are too demanding. However, tidal breathing flow profiles also contain information regarding the presence of airway obstruction.
Because tidal breathing variability analysis is still in its early phases as a field, studies so far have taken varied approaches leading to different measurement duration and technique, differing study populations (age and cognitive state, i.e. sleep or awake) and various extracted parameters from the tidal breathing flow-volume (TBFV) curves. Despite the differences in study methods, the overall finding has been that reduced tidal breathing variability is associated with airway obstruction in asthmatic adults (Veiga et al., 2011) and asthmatic children (Seppä et al., 2016), in adults with COPD (Dames et al., 2014;Motamedi-Fakhr et al., 2016;Niérat et al., 2017) and infants with bronchopulmonary dysplasia (Usemann et al., 2019).
Recently, we showed that such lack of variability distinguishes children with recurrent wheeze (asthma) from healthy children also during sleep at night when measured using impedance pneumography (IP) (Seppä et al., 2019). We also discovered that the obstruction-related reduction in tidal breathing variability is mostly associated with a specific part of the expiration, namely, 15-45% of expired volume. Moreover, obstruction-related changes in the TBFV curves were recently shown to be sleep stage-associated (Gracia-Tabuenca et al., 2019).
Based on the aforementioned findings, this IP measurement method may have a significant impact on the clinical practice of assessing lower airway obstructions (wheeze, asthma) in young children who cannot be easily assessed with any other method at the moment. Thus, it is important to explore the underlying neurophysiological mechanisms of a new method like this. The aim of this study is to investigate the sources of variability in expiratory TBFV curves during sleep to better understand the neurophysiological origins of the normal (healthy) TBFV curve variability. In this study, we investigated sleep stages and body position changes as sources for the normal TBFV curve variability. Firstly, we explored if there is temporal synchrony in how the different parts (volume ranges) of the expiratory TBFV curves change throughout the night. Secondly, we assessed how similar the TBFV curves are within and between different sleep stages, and lastly, we analyzed the potential effect of body position changes on the TBFV curve variability.

Subjects
The study recruited children who were clinically referred to a polysomnography (PSG) study and healthy children without any relevant symptoms. The inclusion criteria for the referred group was age 1-7 years, both sexes, referral to PSG and a signed informed consent. The inclusion criteria for the healthy children were: age 2-7 years, both sexes, healthy at the time of inclusion based on history and clinical examination according to the investigator's judgment and the signed informed consent. All the children that met the inclusion criteria were in the range of normal regarding the apnea-hypopnea index, since subjects who had signs of obstruction or sleep apnea during the PSG were excluded.
The main exclusion criteria were history or risk of asthma, chronic conditions that may alter breathing pattern, continuous upper airway airflow limitation, or other significant clinical findings in the PSG (full list available at https://clinicaltrials.gov/ct2/show/NCT03408990).
One overnight measurement was done for each subject in a sleep laboratory at the Children's Hospital Srebrnjak, Zagreb, Croatia between January 2018 and December 2018. The study was approved by the local ethics committee, guardians of all subjects gave a written informed consent, and the Declaration of Helsinki and Good Clinical Practice were followed.

Measurement setup
The overnight measurement included a full PSG recording with video and a simultaneous TBFV curve shape measurement using IP equipment. IP is a method for indirect measurement of lung tissue aeration (volume) changes through skin electrodes. With correct electrode placement strategy (Seppä et al., 2013a) and signal processing (Seppä et al., 2011), the method can record TBFV curve shapes accurately in young children (Seppä et al., 2013b;Malmberg et al., 2017).
The PSG study was done by using the EEG-1200 (Nihon Kohden, CA, USA), the signals (6 channel electroencephalogram, electrooculogram, peripheral oxygen saturation, nasal pressure, abdominal and thoracic plethysmography) were analyzed using the Polysmith (Nihon Kohden, CA, USA) and Polaris (Nihon Kohden, CA, USA) softwares by an experienced sleep technician. The PSG-signals were analyzed in 30-second epochs. The IP measurement was done using a Ventica Recorder (Revenio Research Ltd., Finland), providing the following signals: IP, electrocardiogram (ECG) and accelerometer. Sleep staging was done according to the current AASM guideline (Berry et al., 2018). Apneic and hypopneic events were identified from the PSG signals. The changes in body position during the night were annotated from a video recording by a human operator.

Signal processing & calculating TBFV curves
IP signal processing to derive TBFV curves was done with a commercially available software (Ventica Analytics 2.0.1, Revenio Research Ltd., Finland) and the subsequent analysis of the curves was done using MATLAB (R2017b, MathWorks Inc., USA).
The Ventica Analytics software provides averaged TBFV curves, calculated by averaging raw TBFV curves in a moving 5-minute window, with 2.5-minute steps. Before generating the TBFV curves, distorted parts of the IP, ECG and accelerometer signals were discarded at the beginning of the signal processing. The resulting TBFV curves represent normal breathing, where sighs, coughs, crying etc. have been excluded.
We were only interested in the exhale part of the TBFV curves, based on the findings from previous studies, where it was shown that it is possible to differentiate between healthy subjects and asthmatics with the help of the exhale part (Gracia-Tabuenca et al., 2019;Seppä et al., 2019). Here we analyzed two ranges of interest, namely, 15-45% and 55-85% of expired volume. The earlier (15-45%) part is used in the expiratory variability index (EVI) parameter derived by the Ventica Analytics software and it has been shown to be the best at differentiating between wheezy and healthy children (Seppä et al., 2019). The latter (55-85%) part was included because it presumably represents mostly passive expiration (post-inspiratory muscle activity has ceased) and Gracia-Tabuenca et al. (2019) showed that parameters derived from that part show significant differences between asthma risk categories in young children.
To make sure that the analysis was accurate it had to be ensured that the Ventica-signals and the PSG-signals were synchronized so that the parameters derived and calculated from these signals are comparable. Both the PSG and Ventica measurements were started as simultaneously as possible, leading to at most a difference of some seconds in the start times. The temporal synchronization was done by manually comparing the breathing signals from both sources, finding multiple reference points seen in both signals (for ex. distortion caused by movement) and iteratively searching for the correct time shift for the PSG-signal so that the signals are aligned.

Averaged TBFV curves and their mutual correlation
To ensure that averaged TBFV curves only belong to one sleep stage, we rejected these averaged TBFV curves resulting from 5-minute windows situated partially in two different sleep stages. For each recording, then the Pearson correlation was calculated between all of the accepted averaged TBFV curves for the 15-45% and 55-85% ranges of the exhaled volume. This process results in two separate correlation matrices, which show how similar or different the TBFV curves are between different points in time during the night (Fig. 2). These matrices from all recordings were the base for the different analyses presented in Fig. 1.

Time-synchrony of change in different parts of TBFV curves
In order to assess whether the changes in the TBFV curve shape occur at the same time in different parts of the exhaled TBFV curve, we calculated the Pearson correlations between all the averaged TBFV curves of a measurement night within ranges of 30% of exhaled volume (0-30%, 5-35%, 10-40% all the way to 70-100%). This resulted in a correlation matrix for each range, where the axes are the sleep time and the correlation coefficients indicate how similar or different the averaged TBFV curves were at different points of the night. After establishing the correlations between the TBFV curves within each volume range, the cross-correlation between the correlation matrices of different ranges was calculated to find out if the changes occur synchronously in time in different parts of the TBFV curve. If the correlation between two ranges (i.e. 0-30% and 5-35%) is high, it means that the changes in the TBFV curve happens at similar points in time during the night. If the correlation is close to 0, it means that the TBFV curve changes occur at different time points in the two compared ranges of the curve.

TBFV curve variability within and between sleep stage cycles
To better understand if there are certain types of breathing patterns in the sleep stages and whether they differ from one sleep stage to another, the averaged TBFV curves of the whole night were divided according to sleep stages (wake, N1, N2, N3 and REM) and the crosscorrelations were calculated between the curves. These results show us if the TBFV curve variability differs in different sleep stages. We also extracted the averaged TBFV curve correlations within single sleep stage cycles to understand if the curve patterns behave differently within cycles in comparison to all sleep stages combined.

Effect of body position change on averaged TBFV curves
In addition to purely analyzing the averaged TBFV curves in regards of sleep stages, we also analyzed the effect of changes in body position on the averaged TBFV curves throughout the night. This analysis was broken down into sleep stages, to better see if body position change has different outcomes on breathing in different sleep stages. The analysis was done by taking the last averaged TBFV curve before a change in position (averaging window ends before change in position) and the first curve after (averaging window begins after change in position), and calculating the correlation between the two curves. If the averaged TBFV curve before the position change is in a different sleep stage than the first curve after the change, they were discarded from the analysis. The reason for this is to enable assessing solely body position change effect without simultaneous sleep stage change effect.

Statistical analysis
Pearson correlation was used in all correlation calculations between the averaged TBFV curves. This method was the same as used by Seppä et al. (2019) and was able to differentiate between healthy subjects and asthmatics. All tests of difference between multiple data groups were done using the Kruskal-Wallis test.   correlations) than the earlier part (panel A), 2) effect of sleep stage changes (especially stage 3 non-REM sleep; N3) on the curve shape is more prominent in the early part, and 3) the early part of the TBFV curve during N3 is very similar within one cycle (next to the matrix diagonal), but can be quite different between two cycles of N3 (away from the diagonal).

Results
3.1. Patient characteristics and polysomnographic findings 20 subjects (8 female) aged 4.0 (1.4-6.9) (median and range) years were included in the study. 7 of the children had been referred to a polysomnography (PSG) study on clinical grounds and other were healthy volunteers. N1, N2, N3, REM and wake stages constituted 2.2 (1.7)% (mean, SD), 43.6 (5.4)%, 28.1 (4.5)%, 21.0 (5.1)%, 5.1 (5.5)% of sleep duration, respectively. These agree with previously reported values for proportions of sleep stages during sleep (Traeger et al., 2005). Because N1 constituted such a small part of the sleep time, it was not analyzed any further. Also the impact of apneic/hypopneic events were not analyzed because they were so rare and short in duration.

Time-synchrony of change in different parts of TBFV curves
As expected, the changes in the expiratory TBFV curve occurred more synchronously between overlapping, adjacent parts (above the black diagonal line in Fig. 3) and the synchronicity gradually reduced as the overlap reduced (towards the black line). More interesting is the synchronicity between non-overlapping curve parts. Unexpectedly, the highest degree of asynchrony was not found between the parts that are farthest away from each other (bottom right in Fig. 3), but between the range 15-45% and all the other parts. This means that 15-45% is the most independent range in terms of timing of changes in the curve. In fact, a change in the last parts of the curve (after 55% of exhaled volume) is more likely to be accompanied by a simultaneous change in the earliest parts (before 40%) than in the middle parts of the curve.

Variability of TBFV curve shapes in sleep stages
For both TBFV curve ranges, REM and wake stage showed higher correlations and lower variability (curves more similar to each other) than N2 or N3 (Fig. 4) when comparing TBFV curves from all the cycles of the same sleep stage. Curves from REM and wake stage were also mutually similar to each other. Slightly paradoxically, in the latter curve range, curves during N2 and N3 were better correlated with curves from REM than within N2 or N3 (Fig. 4, panel B). The variability in TBFV curves was significantly different between sleep stages in both curve ranges (p < 0.0001).
Interesting properties were discovered when the curve comparisons were done either only within individual sleep stage cycles or between different cycles of the same sleep stage (Fig. 5); The range 15-45% N3 showed even slightly higher correlations than REM within cycles, but, on the contrary, very low correlation between separate cycles. In practice this means that the 15-45% part of the TBFV curve is remarkably stable during each continuous cycle of N3 sleep, but can have significantly different shape between different cycles of N3 sleep. For the latter range 55-85 % such behaviour was not clearly present. Moreover, the REM and wake stages showed highest median correlations irrespective of analysis type (within/between cycles) for the latter range. Significant differences were found in both curve ranges, when comparing the correlation distributions shown in Fig. 5 (p < 0.0001).

Effect of body position change on TBFV curves
There were expected differences in the rate of movement in different sleep stages with N2, N3, REM and awake having 1.5, 0.4, 2.2, 0.5 movements per hour (mean), respectively.
The correlations in the 15-45% part of curves compared right before and after a movement event (Fig. 6) were lower than the correlations within each cycle of a certain sleep stage (Fig 5) which indicates that body position changes do have an effect on the TBFV curve shape. The effect was, however, smaller than what was caused by a change in sleep stage, especially for N3. Only the latter curve range showed statistical significance when comparing the effect of body position in different sleep stages (p < 0.0001), while in the early part of the curve there was no significance.
For the latter part of the curve the body position change seemed to have less effect on the curve shape (again comparing Figs. 6 and 5). If changes in body position would not have had any effect on the TBFV curve shape, the distributions in Fig. 6 should have been the same as the within sleep cycle distributions in Fig. 5. This analysis only included body position changes that occurred without a simultaneous sleep stage change.

Discussion
Our findings display that the different parts of the TBFV curves show a degree of independence since they change shape asynchronously in time. An interesting result was that the degree of synchronicity is not the smallest between the earliest and latest parts of the curve. Instead, the last parts of the curve show a higher degree of synchronicity with the earliest parts than with the middle parts. This finding corroborates the special properties of the 15-45% part since it was also found to be the most susceptible for a reduction in its variability in the presence of lower airway obstruction (Seppä et al., 2019). This finding may be associated with the fact that the middle parts of the curve have been found to exhibit a higher absolute level of variability than the other parts (Seppä et al., 2019).
We found the curve shape variability to be considerably higher in the early part than in the late part of the curve which is in line with our previous findings (Seppä et al., 2019). This may be to some extent attributed to the approach we have taken in quantifying the variability (linear correlation between curves) and to the general shape of the expiratory TBFV curves. Namely, the flow peak of the curve occurs during the early section whereas the latter section is characterized by a monotonic decay of flow rate. This renders the earlier part of the curve more susceptible to variation (low correlation values between curves) because even if the curve shape remained the same around the peak flow region, small movements in the peak flow location (within the fixed window of 15-45%) cause a marked drop in the curve correlation Fig. 3. Mean time synchrony of changes between different parts in the averaged TBFV curves. High correlation means that the curve parts are likely to change shape at the same time during sleep. Low correlation means that they may change independently of each other. The area below the diagonal black line has only non-overlapping curve parts. The correlation coefficient values have been transformed logarithmically, to emphasize the differences. A. Hult, et al. Respiratory Physiology & Neurobiology 274 (2020) 103352 value. Such an effect is naturally less prominent in the late part of the curve due to the monotonic shape of the curve. This should, however, only affect the general level of the correlation values between the early and late parts and the interpretation of the effect of sleep stages or body position changes is still valid and comparable for both curve parts. The expiratory limb of the TBFV curves is shaped by a multitude of factors: passive, such as the combined recoil of lungs and thoracic cage and the caliber of lower and upper airways (Otis et al., 1950); and active, such as post-inspiratory inspiratory activity (PIIA), glottal breaking, and expiration interruption (Richter and Smith, 2014). These factors have a different contribution to the early and late part of the curves. The early part is most likely dominated by PIIA and glottal breaking, whereas the late part is mostly passive unless expiration is interrupted before reaching resting volume (typical in newborns) (van  der Ent et al., 1998;Hutten et al., 2008). Shee et al. (1985) showed that in awake healthy adults at 23% of the total expiratory time, 50% of the muscle activity is still present and that only at 79% of the expiratory time the activity ceases completely. Similar values were found by Mortola et al. (1984) for sleeping newborns. The decay time of PIIA means that the early section in our analysis is likely to be affected more by PIIA and the latter less if at all, suggesting that the latter curve shape would be more determined by the passive airway and chest wall mechanics. However, also the adductor and abductor muscles in the pharynx and larynx affect the upper airway patency and thus may also affect TBFV curves. They show both phasic (in concert with the respiratory cycles) (England et al., 1985) and tonic (continuous) activity and the activity level depends on sleep stage (Carroll and Donnelly, 2014). When comparing TBFV curves from all cycles of the same sleep stage throughout the night, we found both early and late parts to be generally less variable during REM than NREM sleep. An interesting feature is that, for the early part (15-45%), the lowest variability occurs within the same N3 sleep cycles, but at the same time, the highest variability appears between different N3 cycles of the night. NREM and in particular N3 is characterized by remarkably regular respiration and decreased muscle tone. Regular respiration in N3 is seen as a low variability in breath-to-breath amplitude and in frequency (Rostig et al., 2005), as well as in a low complexity in the short-term (3 minutes) flow signal (Burioka et al., 2003). Our results show that also the expiratory TBFV is regular during the same NREM cycle. This suggests that the above-mentioned factors shaping expiration remain stable within a NREM cycle. However, our results also show that expiration shape changes between NREM cycles. This suggest that there is more than one configuration of factors that produce a regular within cycle respiration. At different NREM cycles across the night, respiratory control may be attracted into different configurations of factors providing regular respiration (Donaldson, 1992), and for this reason, TBFV profiles between NREM cycles appear dissimilar.
The lower variability of the curve shape during REM in comparison to NREM might be unexpected considering that REM sleep is characterized by significantly higher variability in respiration rhythm (Willemen et al., 2014). The source of this variability is caused by REMdependent alterations in the motoneurons relaying the respiratory drive (Horner, 2010). Despite the pronounced breath-to-breath variations in duration and amplitude of the respiratory stimulus, our results show that the expiratory limb remains relatively similar within and between REM cycles. This may be explained by the decrease in tone in upper and lower respiratory musculature during REM, which leads to an increase in upper airways resistance (increased collapsibility), a decrease in the contribution of the rib-cage to respiration, and a lower functional residual capacity (Gaultier, 1995). These neurophysiological limitations during REM translate into a reduction in the dynamic range for the factors shaping the expiratory profile. Constricted factors increase the chances of profiles being similar within and between REM cycles across the night.
Because of the small sample size, it is difficult to analyze the results in the awake state. Putatively, in healthy individuals, muscle tone is high, most of the expiration is controlled (PIIA or glottal breaking) and the last part is passive (Morris et al., 1998). Moreover, the respiratory drive is modulated breath-to-breath to adapt to changing metabolic demand. Therefore, large variations in the early part within and between cycles may be due to a longer and variable active control of expiration. Whereas, low variation in the late part reflect a similar shorter passive exhale.
The interpretation of TBFV as a product of a given configuration of factors also serves to explain why a change in posture leads to higher differences during NREM than REM. During NREM, a postural change may perturb the stable factor configuration and displace respiratory control into another stable factor configuration. If the new configuration is different from the one before the perturbation, curve shapes will Fig. 6. Correlation between averaged TBFV curves right before and after a body position change. Results here should be compared to the correlation level presented in Fig. 5 within sleep stage cycle. Panels A and B represent TBFV curve parts of 15-45% and 55-85% of exhaled volume, respectively. A. Hult, et al. Respiratory Physiology & Neurobiology 274 (2020) 103352 present low correlation. During REM, the reduced dynamic range of the factors increases the chances that after a movement perturbation the curve shape will be similar to the preceding ones.
Since the averaged TBFV curves seemed to behave differently in different sleep stages, as demonstrated by the correlation distributions in Fig. 4, the TBFV curve variability of the whole night is affected by the relative amount of each sleep stage. As discovered in previous studies, the overall variability in tidal breathing patterns (from shorter recordings) was lowered during obstruction (Veiga et al., 2011;Dames et al., 2014;Frey and Bielicki, 2017;Fouzas et al., 2017;Hmeidi et al., 2018;Usemann et al., 2019). A similar lowering of the overall TBFV curve variability of the whole night could be achieved if the proportions of N2 or REM would grow in regards to the other sleep stages, since the averaged TBFV curves are least variable during these sleep stages. However, it is yet unclear whether the obstruction-related reduction in TBFV curve variability is more prominent in certain sleep stages or if it affects all stages equally. Hypothetically, following the above-mentioned interpretation, lower airway obstruction also reduces the dynamic range of some of the factors shaping TBFV. For instance, to ensure airways patency lungs are hyper-inflated by means of PIIA and glottal breaking (Pellegrino and Brusasco, 1997). Constricted factors would reduce the number of possible stable configurations between different NREM cycles, as well as lowering the chances of TBFV profiles being different within a cycle of any sleep stage or after motion perturbations.
It is also worth noting that the proportion of REM sleep decreases with age, taking approximately 30% of the sleep time in toddlers, while in adolescence and adulthood the time spent in REM sleep is around 20-25%. The proportion of REM sleep in the night decreases once more when reaching the age of 70 and onward. From birth to the age of 20, the percentage of time spent in N3 decreases rapidly from approximately 25% to around 15%. N3 sleep practically disappears after reaching the age of 60. (Carroll and Donnelly, 2014;Cardinali, 2018) In addition to natural change due to aging, there are also sleep disorders and other conditions which hinder normal sleep (Anders and Eiben, 1997;Moore et al., 2006;Parish, 2009). All of these things which can affect the amount of time spent in different sleep stages also indirectly affect the breathing patterns, which may result in increased or decreased variability of breathing for the night. This has to be taken into consideration when assessing the cause of decreased tidal breathing variability, since in addition to lower airway obstructions, there could be other factors in play as well.
This study has some limitations. The types of body position changes are not stratified in any way although it is apparent that, for instance, a movement of the leg is likely to have less impact on breathing than complete change in sleeping position from side to back. Additionally, our sample size is relatively small but there was no patient group with lower airway obstruction included, thus providing group homogeneity. To better understand the mechanisms that drive the known reduction in tidal breathing variability in presence of lower airway obstruction future sleep studies should include patients with uncontrolled asthma and monitor their respiratory muscle activity in addition to PSG.

Conclusions
The late part (55-85% of exhaled volume) of the TBFV curves is significantly more stable than the early part (15-45%). Changes in the curve shape show higher temporal synchronicity between the beginning and end of the curve than with the 15-45 % part. For the early part of the TBFV curves, cycles of REM sleep are less variable than those of NREM, when comparing all cycles of the same sleep stage throughout the night. However, when analyzing the variability only within individual cycles of the same sleep stage, N3 is the least variable. Similar behaviour is not observed for the late part of the curve where REM is less variable than the NREM stages both between and within sleep stage cycles. Body position changes affect predominantly the early part of the curve, but not as strongly as the sleep stage changes, which based on our findings are the main source of variability in the TBFV curves of the night in healthy subjects.

Conflict of interest statement
The authors declare that this study received funding from Revenio Research Ltd. The funder had the following involvement with the study: study design, data analysis, decision to publish and preparation of the manuscript. AH is an employee of Revenio Group Inc. that commercializes impedance pneumography technology. RGJ has nothing to disclose. JGT has nothing to disclose. MP has nothing to disclose. DP reports a grant and consultancy fees from Revenio Group during the conduct of the study, and grants from GlaxoSmithKline, Boehringer Ingelheim, Novartis and Chiesi, grants and personal fees from MSD, and personal fees from AbbVie, Sandoz, GlaxoSmithKline, Salveo, Menarini/Berlin Chemie and Salvus, outside the submitted work. VPS is an employee of Revenio Group Inc. that commercializes impedance pneumography technology and a shareholder in Tide Medical Ltd. that holds patents relating to impedance pneumography.

Author contributions
AH, VPS and DP conceived and designed the study. RGJ and DP collected the data. AH, VPS and MP designed and performed the analysis. AH, VPS and JGT wrote the paper. All authors contributed to manuscript revision, read and approved the submitted version.

Funding
The study was funded by Revenio Research Ltd.