Assessing cardiorespiratory interactions via lagged joint symbolic dynamics during spontaneous and controlled breathing

Introduction: Joint symbolic analysis (JSA) can be utilized to describe interactions between time series while accounting for time scales and nonlinear features. JSA is based on the computation of the rate of occurrence of joint patterns built after symbolization. Lagged JSA (LJSA) is obtained from the more classical JSA by introducing a delay/lead between patterns built over the two series and combined to form the joint scheme, thus monitoring coordinated patterns at different lags. Methods: In the present study, we applied LJSA for the assessment of cardiorespiratory coupling (CRC) from heart period (HP) variability and respiratory activity (R) in 19 healthy subjects (age: 27–35 years; 8 males, 11 females) during spontaneous breathing (SB) and controlled breathing (CB). The R rate of CB was selected to be indistinguishable from that of SB, namely, 15 breaths·minute−1 (CB15), or slower than SB, namely, 10 breaths·minute−1 (CB10), but in both cases, very rapid interactions between heart rate and R were known to be present. The ability of the LJSA approach to follow variations of the coupling strength was tested over a unidirectionally or bidirectionally coupled stochastic process and using surrogate data to test the null hypothesis of uncoupling. Results: We found that: i) the analysis of surrogate data proved that HP and R were significantly coupled in any experimental condition, and coupling was not more likely to occur at a specific time lag; ii) CB10 reduced CRC strength at the fastest time scales while increasing that at intermediate time scales, thus leaving the overall CRC strength unvaried; iii) despite exhibiting similar R rates and respiratory sinus arrhythmia, SB and CB15 induced different cardiorespiratory interactions; iv) no dominant temporal scheme was observed with relevant contributions of HP patterns either leading or lagging R. Discussion: LJSA is a useful methodology to explore HP–R dynamic interactions while accounting for time shifts and scales.


Introduction
Cardiorespiratory coupling (CRC) is defined as the interaction between cardiac and respiratory rhythms. CRC covers a variety of phenomena independent of each other (Penzel et al., 2016;Elstad et al., 2018). Respiratory sinus arrhythmia (RSA) is the most studied effect of CRC (Hirsch and Bishop, 1981;Eckberg, 1983;Hayano et al., 1996). RSA accounts for the periodical modification of heart period (HP) with respiratory activity (R) resulting from HP shortening during inspiration and HP lengthening during expiration under physiological conditions and breathing rates. RSA is a frequency-dependent phenomenon (Angelone and Coulter, 1964;Hirsch and Bishop, 1981;Eckberg, 1983;Brown et al., 1993), but it is not affected by voluntary control of breathing (Patwardhan et al., 1995;Pinna et al., 2006). Another phenomenon under the umbrella of CRC is cardiorespiratory phase synchronization, which takes the form of short intermittent epochs of stable occurrence of heartbeats at specific phases of the respiratory cycle (Schäfer et al., 1998;Schäfer et al., 1999;Bartsch et al., 2012;Kuhnhold et al., 2017;Mazzucco et al., 2017;Cairo et al., 2021;Ottolina et al., 2023). Cardiorespiratory phase synchronization has been observed in athletes (Schäfer et al., 1998;Cairo et al., 2021), healthy subjects (Schäfer et al., 1999;Bartsch et al., 2012), patients in intensive care units (Mazzucco et al., 2017;Ottolina et al., 2023), and post-infarction patients (Kuhnhold et al., 2017). Cardioventilatory coupling, defined as a stable temporal relationship between the onset of the inspiration and the last heartbeat preceding it (Galletly and Larsen, 1997;Tzeng et al., 2003;Larsen et al., 2010;Friedman et al., 2012), is another form of CRC independent of RSA and cardiorespiratory phase synchronization.
Therefore, it is clear that the cardiac and respiratory systems interact in complex and multifaceted ways, including various nonlinear types of interplay (Elstad et al., 2018;Abreu et al., 2023), and this influence might be mediated by the baroreflex as well (Dick et al., 2005;Abreu et al., 2020). As such, many bivariate methodologies that are capable of accounting for high-order statistical moments have been proposed in the context of CRC assessment (Moser et al., 1995;Galletly and Larsen, 1997;Schäfer et al., 1998;Penzel et al., 2016). Among them, joint symbolic analysis (JSA) was first proposed as an extension of univariate symbolic analysis (Baumert et al., 2002;Wessel et al., 2009;Kabir et al., 2011;Schulz et al., 2013a;Porta et al., 2015a;Porta et al., 2015b) and applied to beat-to-beat HP variability and R for the quantification of CRC in sports, health, and pathology (Kabir et al., 2011;Reulecke et al., 2012;Schulz et al., 2013b;Baumert et al., 2015;Schulz et al., 2015;Reulecke et al., 2018;Schulz et al., 2018;Abreu et al., 2019). The introduction of a time shift between patterns, thus leading to a lagged JSA (LJSA), allows monitoring CRC with different temporal schemes and latencies (Wessel et al., 2009;Suhrbier et al., 2010;Kabir et al., 2011;Baumert et al., 2015). Because some JSA tools can account for the time scales when analyzing interactions by classifying joint symbolic schemes composed of patterns featuring different frequency contents (Porta et al., 2015b;Bari et al., 2016;Porta et al., 2016), LJSA could be made scale-specific.
The temporal scheme in HP-R dynamic interactions remains poorly elucidated. Modeling approaches describe cardiorespiratory interactions because of the exogenous actions of R on HP (Baselli et al., 1994;Triedman et al., 1995;Porta et al., 2000a;Chen and Mukkamala, 2008;Porta et al., 2012a;Porta et al., 2013a). This hypothesis was followed even when applying JSA, given that in Kabir et al. (2011), the delay of HP to R was optimized, while advancements of HP to R were not considered. However, the leading action of the heart to the respiratory system has been identified (Galletly and Larsen, 1997;Galletly and Larsen, 1999) and confirmed in healthy individuals (Tzeng et al., 2003;Larsen et al., 2010;Friedman et al., 2012). This link was supported via crossspectral analysis, especially at breathing rates below 0.15 Hz, while in the range of respiratory rates between 0.15 to 0.25 Hz, heart rate fluctuations and R volume appear to be in phase (Saul et al., 1989). More sophisticated parametric modeling approaches corroborate the finding that HP might lead R (Yana et al., 1993;Perrott and Cohen, 1996;Porta et al., 2013b). This effect is compatible with the fastness of cardiac neural control compared to the slowness of the thoracic movements (Yana et al., 1993;Triedman et al., 1995;Perrott and Cohen, 1996;Porta et al., 2013b). However, some studies provided a less trivial interpretation involving the action of a latent confounder, such as the baroreflex (Abreu et al., 2020) triggering the onset of inspiration (Galletly and Larsen, 1999;Dick et al., 2005). We tested the hypothesis that the association between heartbeat and R could be higher at specific time shifts and that this result might depend on breathing rate.
The aim of the present study is the evaluation of the HP-R association in healthy subjects during SB and CB using an LJSA approach capable of differentiating interactions at different time scales. The rate of CB is selected to be indistinguishable from that of SB, namely, at 15 breaths·minute −1 (CB15), or significantly slower, namely, 10 breaths·minute −1 (CB10) but above the limit of 0.15 Hz. The protocol includes CB10 and CB15 because both paced R rates are in the range of R frequencies leading to very fast, in-phase modifications of heart rate and R volume (Saul et al., 1989), even though interactions occur at significantly different time scales that are slower in the case of CB10 than CB15. Accordingly, we tested a physiological range of lags, accurately chosen according to the type of collected signals, to explore variations of the HP-R association with time shift. The link of markers derived from LJSA and coupling strength was tested over simulations, and surrogate data were exploited to check the significance of interactions. Preliminary results were presented at the 12th conference of the European Study Group on Cardiovascular Oscillations (Cairo et al., 2022).

LJSA
The two series HP HP i , i 1, . . . , N { } and R R i , { i 1, . . . , N}, where i is the cardiac beat counter and N is the series length, were separately converted into symbols via a uniform quantization procedure that covered the max-min range over ξ quantization bins (Porta et al., 2001;Cysarz et al., 2013). From the symbolic series HP ξ HP ξ i , i 1, . . . , N and R ξ R ξ i , i 1, . . . , N , we built the series of patterns of length L, HP ξ L HP ξ i,L , i L, . . . , N and are the patterns defined according to the technique of delay embedding. The parameters ξ = 6 and Frontiers in Network Physiology frontiersin.org L = 3 were set according to Porta et al. (2001) to ensure the consistency of the estimates of the symbolic pattern rates over short series of HP and R with N = 256. Each symbolic pattern was classified into one of four classes according to the variation between adjacent symbols (Porta et al., 2001): i) the no variation (0V) pattern when all the three symbols were equal; ii) the one variation (1V) pattern when two adjacent symbols were equal but the remaining one (i.e., the first or the third) was different; iii) the two like variation (2LV) pattern when the three symbols were different and showed a progressive increase, or decrease, from the first to the third one; iv) the two unlike variation (2UV) pattern when the two adjacent symbols were different with the second symbol of the pattern lower or higher than the first and the third one.
This approach allows the preservation of the amplitude features of the series according to the adopted coarse-graining procedure as well as the classification of patterns according to their temporal scale, given that the 0V pattern is the slowest one, because of its stable behavior, and the 2UV pattern is the fastest one, owing to its rapid changes. 1V and 2LV exhibited an intermediate range of time scales, with the 1V pattern slower than the 2LV one (Porta et al., 2001). The symbolization allowed the assessment of the autonomic control governing changes of HP (Guzzetti et al., 2005). The joint HP-R pattern was built by associating a pattern built over R, namely, R ξ i,L , with a pattern built over HP τ lags ahead, namely, HP ξ i+τ,L , where τ is an integer value, thus allowing the assessment of the association between HP and R with HP leading R for τ < 0, immediate HP-R association when τ 0 (i.e., within the current HP), and association between HP and R with HP lagging R for τ > 0.
We defined a joint coordinated (C) pattern associating one pattern built over HP and one pattern built over R when both belonged to the same family (Porta et al., 2015b). C patterns were labeled as 0V-0V, 1V-1V, 2LV-2LV, and 2UV-2UV. The percentage of C patterns (C%) over the total number of joint schemes (i.e., N − L − |τ| + 1) was computed. C% was deemed to be an index of the overall association between the HP and R series regardless of the time scales of the HP-R dynamic interactions. C% ranged from 0 (i.e., no coordinated behavior) to 100 (i.e., all patterns are coordinated) (Porta et al., 2015b). The percentages of each of the four joint pattern classes were computed over the total number of C joint patterns for each τ value so that their sum was always 100%. These markers were denoted as 0V-0V%, 1V-1V%, 2LV-2LV%, and 2UV-2UV%. These four indexes were deemed to assess the association between the HP and R series at a specific time scale of interactions set by the type of joint pattern: HP-R dynamic interactions occurred at the slowest and fastest scales in the presence of 0V-0V and 2UV-2UV joint patterns, respectively, and at time scales intermediate between those described by the 0V-0V and 2UV-2UV joint schemes when the 1V-1V and 2L-2LV joint patterns were detected. Indexes ranged from 0 (i.e., that specific scale did not contribute to the overall HP-R coordination) to 100 (i.e., that specific scale fully explained the overall HP-R coordination) (Porta et al., 2015b).
In-phase and out-of-phase patterns, usually referred to as symmetric and diametric patterns (Wessel et al., 2009), were merged within the same class because classification did not account for the sign of the variations (Porta et al., 2001). In the range of breathing rates considered in the present study, the phase between heart rate and R volume was expected to be about 0 (Angelone and Coulter, 1964;Eckberg, 1983;Saul et al., 1989). Therefore, because we recorded HP, being in phase opposition with heart rate, and a nasal flow, being in quadrature with the R volume, the expected phase shift between HP and our R signal was π/2. This phase shift at a rate of 10 and 15 breaths·minute −1 was equivalent to latencies of 1.5 s and 1 s, respectively, namely, less than two beats at the cardiac frequencies of the present study in any experimental condition. Therefore, the selected interval for τ values was −2 ≤ τ ≤ + 2.

Simulations
To validate the ability of the LJSA approach to follow changes in the coupling strength, we simulated two stochastic processes that were interconnected unidirectionally or bidirectionally via a parameter modulating the coupling strength . The two processes Y 1 and Y 2 are defined as where W 1 and W 2 are Gaussian white noises with zero mean and variances assigned such that Y 1 and Y 2 exhibit unit variance. Three configurations are of value: i) full uncoupling between Y 1 and Y 2 with c 1 c 2 0 and Y 1 and Y 2 are two second-order autoregressive processes; ii) unidirectional interactions from Y 1 to Y 2 with c 1 0 and c 2 ≠ 0; iii) bidirectional interactions from Y 1 to Y 2 and vice versa with c 1 ≠ 0 and c 2 ≠ 0. In the uncoupling condition, Y 1 and Y 2 were set to feature the same dominant rhythm according to ρ 1 ρ 2 0.85 with phases φ 1 φ 2 ± 3π/10 corresponding to a dominant oscillation at the normalized frequency f 1 f 2 0.15 cycles · sample −1 , thus simulating the typical dynamics usually observed in HP and R series during SB and CB (i.e., 0.15 Hz with a mean HP equal to 1 s). In the unidirectional configuration from Y 1 to Y 2 , the coupling parameter c 2 was varied incrementally from 0 to 1.0 with 0.1 steps with c 1 0. In the bidirectional configuration from Y 1 to Y 2 and vice versa, c 1 c 2 c, and c was varied gradually from 0 to 1.0 with 0.1 steps.

Experimental protocol and data analysis Experimental protocol
The estimate of CRC via LJSA was performed on a historical database designed to quantify cardiorespiratory interactions while varying the breathing rate (Porta et al., 2000b;Porta et al., 2012a). The protocol was in keeping with the Declaration of Helsinki and approved by the Local Ethical Review Board of L. Sacco Hospital, Milan, Italy (protocol code: 1999-3; date of approval: 1/2/1999). Written signed informed consent was obtained from all subjects. We enrolled 19 healthy subjects (age: 27-35 years, median = 31 years; 8 males, 11 females). The health status of subjects was verified via physical examination, evaluation of the historical personal records, blood pressure measurement with a sphygmomanometer, and Frontiers in Network Physiology frontiersin.org standard assessment of the electrocardiographic trace by an expert cardiologist. During the experimental sessions, we acquired the electrocardiogram (ECG) from lead II via a bioamplifier (Marazza, Monza, Italy) and the R flow via a nasal thermistor (Marazza, Monza, Italy) at the sampling frequency of 300 Hz. The two signals were synchronized through a 12-bit analog-todigital board (National Instruments, Austin, TX, United States) plugged into a personal computer. Recordings were made at rest in a supine position during SB, CB10, and CB15. The SB session always preceded CB recordings, while the order of the CB10 and CB15 sessions was randomized. The subjects were trained to follow a computer-based metronome that provided the pace for starting both the inspiratory and expiratory phases. The inspiratory-to-expiratory ratio was set to a traditional 1:2 (De Maria et al., 2021) via the metronome sounds. The indication provided by the metronome was reinforced through verbal commands by the experimenter. Sessions lasted 10 min, and subjects were not allowed to talk for the entire duration of the protocol.

Beat-to-beat variability series extraction
The R-wave peaks were identified on the ECG through a thresholdbased algorithm applied to the ECG first derivative. The ith HP value was then estimated as the time interval occurring between the ith and (i + 1)th R-wave peaks. The R signal was sampled at the ith R-wave peak. Sequences of 256 consecutive values were selected within each experimental condition, thus focusing on short-term cardiac control (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996). Detections of the R-wave peak were visually checked and manually corrected if necessary. The effect of any isolated arrhythmic beat was mitigated through linear interpolation using the two most adjacent HPs computed between sinus beats. The 5% limit of corrections was never reached. From the HP series, we computed the mean, labeled as μ HP and expressed in ms.

Frequency domain analysis
Power spectral density was estimated separately over the HP and R series via a parametric approach based on the identification of the coefficients of an autoregressive model (Baselli et al., 1997). The least squares problem was solved via the Levinson-Durbin recursion (Kay and Marple, 1981). The model order was optimized via the Akaike information criterion between 10 and 16 (Akaike, 1974). The traditional method based on the decomposition of the power spectral density according to the residue theorem was utilized to compute the power associated to each spectral component (Baselli et al., 1997). The highfrequency band (HF, 0.15-0.4 Hz) was considered for the computation of the RSA from the HP series and the assessment of the respiratory frequency (f R ) from the R series. RSA was estimated by summing the power of all spectral components whose central frequency fell within the HF band, and this index was labeled as HF HP and expressed in ms 2 . From the R series, we extracted the f R as the central frequency of the dominant spectral component present in the HF band.
Power cross-spectral density was estimated over the HP and R series using a parametric approach based on the identification of the bivariate autoregressive model (Baselli et al., 1997). The least squares problem was solved via the Cholesky decomposition method (Porta et al., 2000a). The model order was fixed to 10. The power cross-spectral density was computed by taking the R series as the input and HP as the output. The phase φ HP−R (f) of the power cross-spectral density represented the phase shift from R to HP as a function of the frequency f, where negative values indicated that HP lagged R, and positive values indicated that HP led R. The phase was expressed in radians (rad) and ranged from −π to +π. The phase function was sampled at

Surrogate data approach
We utilized a surrogate data approach to test the null hypothesis of uncoupling between HP and R at a given time lag τ (H 0,τ ) and the null hypothesis of full uncoupling between HP and R regardless of the time lag (H 0 ). We constructed sets of surrogate pairs preserving the distribution and power spectral density of the original series (Schreiber and Schmitz, 1996) but fully adherent with H 0,τ and H 0 (Palus, 1997). We generated 100 realizations for each original pair of HP and R series in any experimental condition. The surrogate series were built to preserve the distribution and power spectral density of the original series, while phases were substituted with uniformly distributed random numbers ranging from 0 to 2π. Surrogate pairs were generated through an iteratively refined, amplitude-adjusted, Fourier transformbased procedure (Schreiber and Schmitz, 1996). The technique allowed the exact preservation of the original distribution, while the power spectral density was the best approximation of the initial power spectral density given 100 iterates. The use of two independent random phase sequences allowed the generation of fully uncoupled pairs (Palus, 1997). A fast Fourier transform procedure was applied to speed the construction of surrogates. The percentages of the C pattern families (i.e., 0V-0V%, 1V-1V%, 2LV-2LV%, and 2UV-2UV%) were computed over each set of surrogates, and the 95th percentile was extracted. At each time lag, if the percentage of the C pattern family computed over the original series was above the 95th percentile of the same index derived from surrogates and if this occurred for at least one of the 4°C pattern classes, the H 0,τ was rejected for the considered τ, and the alternative hypothesis, namely, the two series were significantly associated at the selected τ, was accepted. The H 0 was rejected if H 0,r was rejected in correspondence with at least one of the considered time lags (i.e., τ = −2, τ = −1, τ = 0, τ = +1, or τ = +2), and the alternative hypothesis, namely, the two series were significantly associated regardless of time shift, was accepted. The percentage of H 0 rejections was evaluated in any experimental condition (i.e., SB, CB10, and CB15).

Statistical analysis
The normality of data was checked via the Shapiro-Wilk test. The one-way repeated measures analysis of variance versus control (Dunnett's test for multiple comparisons), or the one-way Friedman repeated measures analysis of variance on ranks versus control (Dunn's test for multiple comparisons), when appropriate, was used to separately evaluate the effect of breathing modality (i.e., SB, CB10, Frontiers in Network Physiology frontiersin.org and CB15) and the effect of time shift (i.e., τ = −2, τ = −1, τ = 0, τ = +1, and τ = +2). When assessing the impact of the breathing modality on LJSA markers, data were pooled regardless of the time lag. When assessing the impact of the time shift on LJSA markers, data were pooled regardless of the experimental condition. The two-way repeated measures analysis of variance versus control (Holm-Sidak test for multiple comparisons) was utilized to check the impact of the experimental condition within the same time lag and the effect of time lag within the same experimental condition. In the case of both the one-way and the two-way analyses, the control conditions were SB and τ = 0. The χ 2 test (McNemar's test) was applied to the proportion of subjects featuring the rejection of H 0 to assess the impact of the CB versus SB. The level of significance of the test was lowered according to the number of comparisons (i.e., 2) to account for the multiple comparison issue. The same test was utilized to check the impact of CB versus SB within the same time lag and of the time lag versus τ = 0 within the same experimental condition on the proportions of subjects featuring the rejection of H 0,τ . In this case, the level of significance was lowered by a factor accounting for all the considered comparisons. Statistical analysis was carried out using a commercial statistical program (SigmaPlot, v.14.0, Systat Software, Inc., Chicago, IL, United States). A p < 0.05 was always considered significant.

Results over simulations
The vertical grouped error bar graphs of Figure 1 summarize the results of simulations of unidirectionally coupled ( Figures 1A, C, E,

FIGURE 1
The vertical grouped error bar graphs show C% (A,B), 0V-0V% (C,D), 1V-1V% (E,F), 2LV-2LV% (G,H), and 2UV-2UV% (I,J) computed over simulated processes as a function of c 2 . Results were obtained in the case of unidirectional interactions with c 1 0 (A,C,E,G,I) and bidirectional interactions with c 1 c 2 (B,D,F,H,J) from 20 realizations of Y 1 and Y 2 . Results were reported for τ −1 (black bars) and τ +1 (white bars). At c 2 0, Y 1 and Y 2 were uncoupled. The symbols § and * indicate a significant difference with p < 0.05 versus c 2 0 within the same time shift, being τ −1 and τ +1, respectively. The symbol # indicates a significant difference with p < 0.05 between the time lags within the same value of c 2 . Symbols indicating p < 0.05 were reported solely when the null hypothesis of uncoupling between Y 1 and Y 2 was rejected (i.e., the value of the markers was significantly above the one found at c 2 0). Data are reported as mean plus standard deviation. In the case of unidirectionally coupled processes, C%, 2LV-2LV %, and 2UV-2UV% computed at τ +1 were larger than the level set by the uncoupling condition (i.e., at c 2 0) when c 2 ≥ 0.4, c 2 ≥ 0.2, and c 2 ≥ 0.2, respectively, and these indexes were larger than the value obtained at τ −1 when c 2 ≥ 0.4, c 2 ≥ 0.2, and c 2 ≥ 0.5. Remarkably, 2UV-2UV% increased gradually with c 2 . 0V-0V% and 1V-1V% remained less than the value set by the uncoupling condition (i.e., at c 2 0), regardless of the time shift, because the link between the two processes occurred at fast time scales. 2UV-2UV% at τ −1 was larger than the value set by the uncoupled realizations as well, but it remained less than its value at τ +1, and the rate of the increase was less steep than that at τ +1.
In the case of bidirectionally coupled processes, C% and 2UV-2UV% computed at both τ ± 1 were larger than the level set by the uncoupling condition (i.e., at c 2 0) when c 2 ≥ 1.0 and c 2 ≥ 0.2, respectively. C% and 2UV-2UV% were alike when computed with τ ± 1. As in the case of the unidirectionally coupled processes, 2UV-2UV% increased gradually with c 2 . 0V-0V%, 1V-1V%, and 2LV-2LV% remained less than the value set by the uncoupling condition (i.e., at c 2 0), regardless of the time shift.

Impact of CB on time, spectral and crossspectral markers
The results of the time and frequency domain analysis of the HP variability and R signal are reported in Table 1. The μ HP and HF HP powers were similar across experimental conditions. f R did not change during CB15 compared to SB, while it was significantly slower during CB10.
The vertical box-and-whisker plots of Figure 2 show crossspectral markers, namely, φ HP−R (f R ) ( Figure 2A) and τ HP−R (f R ) ( Figure 2B), as a function of the experimental condition (i.e., SB, CB10, and CB15). The height of the box represents the distance between the first and third quartiles, with the median marked as a horizontal segment, and the whiskers denote the 5th and 95th percentiles. The phase is expressed in rad ( Figure 2A) and converted into a delay/advancement expressed in s in Figure 2B. Significance was tested against SB. Phases were more likely to be negative, thus indicating that, most frequently, HP lagged R. The median values were −2.15, −1.98, and −1.78 rad during SB, CB10, and CB15, respectively ( Figure 2A). After converting phase values into lags ( Figure 2B), the delay of HP to R was more pronounced during CB10 than SB and CB15, with median values of −0.71, −1.31, and −0.69 beats during SB, CB10, and CB15, respectively, and individual values being more negative than −2 beats in 5%, 32%, and 0% and never more positive than +2 beats during SB, CB10, and CB15, respectively.

Results of the surrogate data test
The simple bar graphs of Figure 3A show the percentage of H 0 rejections as a function of the experimental condition, while the grouped bar graph of Figure 3B summarizes the percentage of H 0,τ rejections as a function of the experimental condition at any time lag τ represented with bars of different filling color from light gray to black. Regardless of the experimental condition, HP and R series were found to be significantly coupled (above 79%), and the percentage of H 0 rejections peaked at 95% and 84% during CB10 and CB15, respectively ( Figure 3A). Remarkably, the percentage of H 0,τ rejections did not depend on either the time lag or experimental condition, with values ranging from 37% to 63% ( Figure 3B).

Impact of CB and time lag on LJSA markers
The simple error bar graph of Figure 4A shows C% as a function of the time lag after pooling all the data regardless of the experimental condition, while the simple error bar graph of Figure 4B shows C% as a function of the experimental condition after pooling all the data regardless of the time shift. The grouped error bar graph of Figure 4C summarizes C% as a function of the experimental condition when the time shift is represented by different filling colors of the bars from light gray to black. Significance was tested against τ = 0 and SB. C% decreased at τ = +1 compared to τ = 0 ( Figure 4A), did not vary with experimental condition (Figure 4B), and increased at τ = −1 during CB10 compared to SB ( Figure 4C). The simple error bar graphs of Figure 5 report 0V-0V% ( Figure 5A), 1V-1V% ( Figure 5B), 2LV-2LV% ( Figure 5C), and 2UV-2UV% ( Figure 5D) as a function of the time shift when data are pooled regardless of the experimental condition. Significance was tested against τ = 0. 0V-0V% decreased at τ = +1 ( Figure 5A), while 1V-1V%, 2LV-2LV%, and 2UV-2UV% did not vary with τ with respect to τ = 0 ( Figures 5B-D). Figure 6 has the same structure as Figure 5, but the percentages of pattern classes are given as a function of the experimental condition when data are pooled regardless of the time lag. Significance was tested against SB. 0V-0V% significantly decreased solely during CB15 ( Figure 6A), while the decrease of 2UV-2UV% was evident solely during CB10 ( Figure 6D). 2LV-2LV% increased both during CB10 and CB15 ( Figure 6C). CB did not affect 1V-1V% ( Figure 6B).

Discussion
The main findings of this study can be summarized as follows: ii) LJSA was useful to quantify the CRC strength while testing different The vertical box-and-whisker plots show φ HP−R (f R ) (A) and τ HP−R (f R ) (B) as a function of the experimental condition (i.e., SB, CB10, and CB15). The height of the box represents the distance between the first and third quartiles, with the median marked as a horizontal segment, and the whiskers denote the 5th and 95th percentiles. The symbol * indicates p < 0.05 versus SB.

FIGURE 3
The simple bar graphs show the percentage of H 0 rejections as a function of the experimental condition (i.e., SB, CB10, and CB15) (A), and the grouped bar graphs show the percentage of H 0,τ rejections as a function of the experimental condition with the time lag τ coded according to the filling color of the bar from light gray to black (B). No significant differences were detected across experimental conditions and lags.

FIGURE 4
The simple error bar graphs show C% as a function of the time lag (i.e., to τ = −2, τ = −1, τ = 0, τ = +1, and τ = +2) when data are pooled regardless of the experimental condition (A) and C% as a function of the experimental condition (i.e., SB, CB10, and CB15) when data are pooled regardless of the time lag (B). The grouped error bar graphs show C% as a function of the experimental condition with τ coded according to the filling color of the bar from light gray to black (C). The symbol # indicates a significant between-time lag difference versus τ = 0 with p < 0.05. Data are reported as mean plus standard deviation.

Frontiers in Network Physiology
frontiersin.org CRC from spontaneous fluctuations of HP and R according to the implementation proposed in Porta et al. (2015b). The characterization is grounded on the assessment of the probability of joint patterns, expressed as a percentage, describing coordinated behaviors between HP changes and R (i.e., C%). The C patterns were classified according to the frequency content of the features

FIGURE 6
The vertical grouped error bar graphs show 0V-0V% (A), 1V-1V% (B), 2LV-2LV% (C), and 2UV-2UV% (D) as a function of experimental condition (i.e., SB, CB10, and CB15). Data are pooled regardless of the time lag between HP and R series. The symbol * indicates a significant between-experimental condition difference versus SB with p < 0.05. Data are reported as mean plus standard deviation.
Frontiers in Network Physiology frontiersin.org comprising the joint behavior from the slowest (i.e., 0V-0V) to the most rapid (i.e., 2UV-2UV) passing through intermediate levels of rapidity (i.e., 1V-1V and 2LV-2LV), thus allowing the assessment of the contribution of time scales to the total coordination. The joint pattern is formed by imposing a time shift between the symbolic patterns, built separately over the two series and combined to form the joint pattern in such a way as to test different temporal schemes of cardiorespiratory interactions (Baumert et al., 2002;Wessel et al., 2009;Kabir et al., 2011): negative time shifts imply that HP leads R, no time shift implies HP-R immediate interaction, namely, associations between HP and R within the current HP, and positive time shifts imply that HP lags behind R. According to the definition of patterns built separately over the two series, both inphase and out-of-phase schemes are present within the same LJSA class (Wessel et al., 2009;Suhrbier et al., 2010). The selected approach has three potential strengths: i) being a model-free approach, it can describe nonlinear interactions between HP variability and R (Baumert et al., 2002;Wessel et al., 2009;Kabir et al., 2011;Porta et al., 2015a;Porta et al., 2015b); ii) different temporal schemes can be easily introduced by setting the delay/ advancement between patterns built separately over HP variability and R (Wessel et al., 2009;Suhrbier et al., 2010;Kabir et al., 2011;Baumert et al., 2015); iii) different time scales of interactions can be explored according to the definition of the joint patterns (Porta et al., 2015b;Bari et al., 2016;Porta et al., 2016). The range of lags explored in the present analysis was limited according to previous studies (Angelone and Coulter, 1964;Eckberg, 1983;Saul et al., 1989). Angelone and Coulter (1964) suggested that heart rate and thoracic movements are in-phase at f R = 10 breaths·minute −1 . Eckberg (1983) extended this observation by observing that the onset of HP shortening occurs progressively earlier during inspiration while slowing R such that it precedes the onset of inspiration at f R = 8 breaths·minute −1 . Through a broadband respiratory approach and cross-spectral analysis between heart rate and R volume, Saul et al. (1989) found that in the range of R frequencies between 0.15 and 0.25 Hz, heart rate and R volume interact with minimal lead (i.e., about zero phase), while robust phase leads are observed below 0.15 Hz. Therefore, because we recorded HP, being in phase opposition with heart rate, and an R flow signal, being in quadrature with the R volume, our two series are expected to be in quadrature. Cross-spectral analysis confirmed that the absolute value of the median phase at the R rate was smaller than π/2 in all the experimental conditions, and the values of lags expressed in samples are between −2 and +2 beats. Although the phase values are mostly negative, thus indicating that HP lagged R, the variability of phase values suggested the necessity of exploring even the positive phase, compatible with the observation that HP changes might precede R variations. The relevance of this exploration is supported by previous studies stressing the presence of a pathway from the heart to the respiratory system, even in healthy subjects (Tzeng et al., 2003;Larsen et al., 2010;Friedman et al., 2012). The relevance of this pathway was confirmed by modeling approaches as well (Yana et al., 1993;Perrott and Cohen, 1996;Porta et al., 2013b). The procedure adopted in the present study based on cross-spectral analysis allows the optimization of the range of lags compared to studies where a much wider interval of lags was considered (Galletly and Larsen, 1997), thus limiting the number of statistical comparisons and increasing the statistical power of the analysis. In addition, LJSA could provide a unique possibility of monitoring the evolution of the CRC strength by simultaneously accounting for nonlinearity, delays/ advancements, and time scales of the HP-R dynamic interactions. These characteristics result from the model-free nature of the method, the introduction of the time shift between HP and R patterns, and the types of joint schemes combining HP and R patterns that interact at different frequencies.

Coupling strength and LJSA
Simulations support the ability of LJSA to follow variations of the coupling strength. Indeed, in situations of unidirectional or bidirectional coupling set between processes exhibiting a dominant and fast component, markers describing the association between realizations occurring at the fastest time scale were able to follow the progressive increase of the coupling strength. It is worth noting that an unspecific marker of association such as C% is much less powerful than a very specific index such as 2UV-2UV%, thus corroborating the relevance of following a specific class of joint pattern selected according to whether coupling could occur at faster or slower time scales. However, given the shortness of the considered patterns combined to form the joint scheme, the time scales of interactions remain roughly defined, but this characterization is sufficient to separate different classes of interactions. Differences between markers computed over different joint pattern classes support the relationship of the indexes with the temporal scheme of interactions and their ability to separate leading interactions from lagging ones (e.g., in the simulation of unidirectionally coupled processes). Because 0V-0V, 1V-1V, 2LV-2LV, and 2UV-2UV patterns belong to the class of C pattern, an increased percentage of one or more classes corresponds to a decrease of at least one of the others. This situation is evident in the simulation of bidirectionally coupled processes when 2UV-2UV% increased above the level set by the uncoupling situation and gradually increased with the coupling strength, while all the other classes became insignificant (i.e., less than the level of uncoupling). This effect might be particularly useful when coupling occurs at a dominant time scale, like in CRC, but it might be particularly limiting in situations where the association occurs simultaneously at different time scales.

CB affects CRC even in the absence of RSA modification
The invariable value of C% with the experimental condition might suggest that CRC strength did not vary during CB. This result disagrees with studies that suggested that CRC strength increases while slowing R (Porta et al., 2000b;Porta et al., 2023). However, when the analysis was made more specific by accounting for the time scales of the interactions governing CRC, we observed that 2LV-2LV% increased and 2UV-2UV% decreased during CB10, thus suggesting that slowing the f R diminishes the coordination between HP and R at the fastest time scale, but it improves it at intermediate time scale. It has Frontiers in Network Physiology frontiersin.org been proven that CB at the same f R as the SB did not alter RSA (Patwardhan et al., 1995;Pinna et al., 2006). The present study confirms this finding. However, SB and CB15 could not be considered indistinguishable. Indeed, 0V-0V% decreased and 2LV-2LV% increased during CB15 compared to SB as a likely result of the smaller impact of slow fluctuations on the HP series and a greater regularization of the HP fluctuations occurring at the f R (Porta et al., 2000a). Because the practice of CB is widely utilized to standardize the impact of f R on HP variability analysis (Patwardhan et al., 1995;Pinna et al., 2006), the present finding suggests that CB might have the wanted effects over some linear univariate parameters, such as RSA, but it might be ineffective, and even introduce additional confounding influences, over nonlinear bivariate markers.

No predominant HP-R temporal scheme is observed
As previously observed in Kabir et al. (2011) and Baumert et al. (2015), introducing a lag between HP and R might lead to observing higher values of CRC strength. However, the procedure adopted to maximize the value of CRC strength should also take into consideration that HP might lead R, in addition to lagged interactions of HP to R (Wessel et al., 2009;Suhrbier et al., 2010). Indeed, the pathway from HP to R is plausible because a constant latency between the inspiratory onset and the last heartbeat preceding it was found (Tzeng et al., 2003;Larsen et al., 2010;Friedman et al., 2012), even at physiological breathing rates (Eckberg, 1983;Saul et al., 1989). It is worth noting that the values of lag to be explored should be carefully chosen according to the type of the recorded signals, here HP and R flow, and with the help of specific analyses, here cross-spectral analysis, limiting the interval of values to be tested. The present study suggests that there is no dominant temporal scheme in the HP-R dynamic interactions, given that LJSA markers did not show a clear behavior at a specific time lag. Conversely, the presence of some peaks at both positive and negative time lags suggests the possible presence of closed-loop interactions, given that the open-loop condition should be characterized by the dominant presence of an assigned temporal scheme as suggested by simulations of unidirectionally coupled processes. Therefore, contrary to the warning raised in Larsen et al. (2010), our data do not support the statement that the CB could disrupt the anticipatory action of HP on R and highlight the importance of considering closed-loop models when describing the dynamic relationship between HP and R (Larsen et al., 2010;Porta et al., 2013b) instead of the more traditional approaches that exclusively model the action from R to HP (Baselli et al., 1994;Triedman et al., 1995;Porta et al., 2000b;Chen and Mukkamala, 2008;Porta et al., 2012a;Porta et al., 2012b;Porta et al., 2013a). While the pathway from R to HP is usually considered to be the result of the activity of the respiratory centers in the brain stem modulating vagal and sympathetic motoneuron firing and, in turn, sinus node activity (Gilbey et al., 1984;Eckberg, 2003;Skytioti and Elstad, 2022), fewer hypotheses have been advanced regarding potential physiological mechanisms underpinning interactions from the heart to the respiratory system. It has been proposed that HP changes at the f R precede modifications of the R signal due to the rapidity of vagal cardiac control compared to the slowness of modifications of the R volume (Yana et al., 1993;Triedman et al., 1995;Perrott and Cohen, 1996;Porta et al., 2013b). Less trivially, it has been suggested that modifications of the baroreceptor activity during systole occurring in the late expiration could contribute to starting the inspiratory phase at a preferred latency (Galletly and Larsen, 1999;Dick and Morris, 2004;Dick and Morris, 2004) and modulating the CRC strength (Abreu et al., 2020).

Conclusion
This study evaluates the effect of SB and CB on CRC strength in healthy subjects using LJSA over HP variability and R signal. LJSA accounts for temporal schemes with either HP leading or lagging R. Results show that CB has an impact on CRC strength not only when the breathing rate is slower than that of SB but also when the rates are similar. In addition, the study proves that in healthy subjects, no dominant temporal scheme governing HP-R dynamic interactions is present with relevant contributions of HP patterns either leading or lagging R. This result is compatible with closedloop dynamic interactions between the heart and the respiratory system. More experiments are needed to further elucidate the physiological meaning and mechanisms behind the observed closed-loop interactions and whether results might depend on the strategy adopted for LJSA computation, including symbolization and pattern construction procedures. In addition, future studies should compare the information derived from measures of association, like those proposed in the present study, with indexes more specifically designed to assess the directionality of coupling and causality (Bahraminasab et al., 2008;Staniek and Lehnertz, 2008;Chicarro and Andrzejak, 2009;Wessel et al., 2009;Suhrbier et al., 2010;Li et al., 2011;Martini et al., 2011;Dicken and Lehnertz, 2014).

Data availability statement
The dataset is relevant to the ESGCO 2022 challenge, "Characterizing cardiorespiratory interactions from spontaneous fluctuations of heart period and respiratory flow during controlled breathings," and it can be downloaded without any restriction at https://www.esgco.org/challenges.

Ethics statement
The protocols for the studies involving human participants were reviewed and approved by the Local Ethical Review Board of L. Sacco Hospital, Milan, Italy (protocol code: 1999-3; date of approval: 1/2/1999). The patients/participants provided their written informed consent to participate in this study.

Frontiers in Network Physiology
frontiersin.org