## Temporal orders and causal vector for physiological data analysis

Marcel Mły´ nczak 1

Abstract -In addition to the global parameterand time-series-based approaches, physiological analyses should constitute a local temporal one, particularly when analyzing data within protocol segments. Hence, we introduce the R package implementing the estimation of temporal orders with a causal vector (CV). It may use linear modeling or time series distance. The algorithm was tested on cardiorespiratory data comprising tidal volume and tachogram curves, obtained from elite athletes (supine and standing, in static conditions) and a control group (different rates and depths of breathing, while supine). We checked the relation between CV and body position or breathing style. The rate of breathing had a greater impact on the CV than does the depth. The tachogram curve preceded the tidal volume relatively more when breathing was slower.

The first two approaches (global time-independent, and global time-dependent) have already been checked in [6], [7]. In the former one, we tried to discover causal paths, and the results suggested body position is an essential factor to be taken into account. To bring a context, in a supine position, the values of a tidal volume seemed to cause heart activity variation, which affected average heart activity, and the chain was completed with respiratory timing parameter. For standing, in turn, the relation led from normalized respiratory activity variation to average heart activity [6].

Clinical relevance - The method can assess (1) relationships between two signals, having one of the two time-shifted in a given range, (2) curves of most optimal inter-signal shift, and (c) their stability across time; also for other physiological studies.

## I. INTRODUCTION

The network physiology concept rests on the assumption that several systems are combined in a network of dependencies with various loops, feedbacks, and delays in transmitting information (also influenced by many environmental, psychological, or demographic factors) [1], [2].

One of the combinations is that of respiration and heart activity. It appears multi-directional and very intricate. Respiratory sinus arrhythmia (RSA) is a phenomenon evident in heart rate modulation as the effect of successive inspirations and expirations [3]. On the other hand, cardiorespiratory coupling holds that cardiac activity triggers respiration, particularly because of increased sympathetic nervous activity and afferent baroreflex discharge [3], [4], [5].

Many mathematical approaches have been proposed to analyze relationships. We hypothesized (primarily from the sports medicine, normative perspective) that the parameterization of the cause-and-effect relationships would also be valuable. Therefore, we assumed that the analysis should comprise (1) global parameters examination without considering the impact of time; then (2) it should delve into global temporal (time series) relationships and causalities; and finally, (3) local analysis of narrow time intervals would be performed for the directionality, strength and stability assessments within each protocol segments (e.g., static supine as the first step of orthostatic maneuver).

1 M. Mły´ nczak is with the Warsaw University of Technology, Faculty of Mechatronics, Institute of Metrology and Biomedical Engineering, 8 Boboli Street, 02-525, Warsaw, Poland. marcel.mlynczak@pw.edu.pl

In the latter approach, temporal relationships were examined by (a) Granger causality frameworks (with extensions that consider zero-lag effects [8]) and (b) Time Series using Restricted Structural Equation Models (TiMINo) [9]. The most prominent 'global' combination appeared between the tachogram (RR intervals curve) and tidal volume signal [7], as a result of RSA. However, the results were weak and unstable, which is probably an effect of considering relatively long segments of data. But, on the other hand, temporal causality analyses require minimum lengths of data to work correctly, e.g., by definition and because of the stationarity criterion.

Therefore, there is a need for a framework to examine local, short segments of data and to explore the temporal orders (TO) between signals. The introduction and preliminary evaluation of such a technique is hence the main aim of this study.

## II. MATERIALS AND METHODS

A. Description and interpretation of the algorithm

The method is implemented as an R package, 'tempord', supplementing the paper. The 'installer', the CRANcompatible manual, and the scheme of the algorithm are available on the author's website [10]. The package uses external R packages: ggplot2 and TSdist .

The bivariate dataset should be loaded. The first signal is kept stationary during the analysis, while another is shifted in time (backward and forward, in the given range).

The most important parameters that define the entire setup (all are described in the R package manual), are:

- · method - determines which approach will be used (linear modeling or time series distance; further explained below);
- · scaling -states the standardization type: 0 no standardization, 1 - uniform standardization, or 2 - Gaussian standardization;
- · length of a signal segment to analyze per loop iteration (the same length for both signals; should not exceed the length of the entire signal); and
- · maximum shift in time - determines how far the second signal may be moved backward and forward; the range does not need to be symmetrical to zero shift.

The main loop of the analysis is carried out for each point of the first signal with a shift applied to the second signal (according to the set resolutions). In the first step, the selection of the signal parts is taken for analysis during each iteration. Then follows scaling if selected.

The main parameters are: (1) the adjusted R-Squared measure, for linear modeling (LM) approach, or (2) the time series distance (TD) measure (e.g., Manhattan or Fourier types), for another, TD approach.

After all iterations, the matrix is filled and may be presented in the temporal orders graph. To facilitate assessing stabilities, the curve of maximum or minimum values (for all shifts considered) for each point in time (for linear modeling and time series distance, respectively) -hereafter called causal vector (CV), is calculated. If a threshold (the level at which the estimate of the parameter is excluded when determining CV) is chosen, only the values (respectively) above or below the threshold will produce those curves.

The setting of the right shift range is crucial. If the range is too wide, the impact of consecutive periods would be visible. Also, in general, the causal vector curve should be nonintermittent (when considering thresholds) and near-parallel to the X-axis. In such fragments of the signal, the temporal relation may be regarded as stable. Step changes in the course of the maximum values may reflect extrinsic changes in the study conditions.

## B. Study groups and the protocol

The participants created two groups (demographically described in Table I):

- · 10 elite athletes (A) -4 minutes of registrations for supine and standing body positions, with unconstrained and free breathing protocols (sample choice from the study group described in [6], [7]);
- · a control group of 10 students (B) - following a constrained breathing procedure comprising 12 breaths (6 normal, then 6 deep) each at rates of 6, 10, and 15 breaths per minute (BPM).

The experimental procedures involving human subjects described in this paper were performed according to the principles outlined in the Helsinki Declaration and the part regarding the athletes also approved by the Ethics Committee of Warsaw Medical University (permission AKBE/74/17).

TABLE I

THE DEMOGRAPHIC SUMMARY OF BOTH STUDY GROUPS.

| Group       | A                | B               |
|-------------|------------------|-----------------|
| Count       | 10               | 10              |
| Sex         | 5F & 5M          | 3F & 7M         |
| Weight [kg] | 68 . 5 ± 15 . 3  | 72 . 9 ± 13 . 2 |
| Height [cm] | 179 . 5 ± 12 . 6 | 176 . 2 ± 7 . 2 |

Single-lead ECG and impedance pneumography signals were acquired using the Pneumonitor 2 [11]. The tachogram was calculated from the ECG signal using the identification of R peaks based on the Pan-Tompkins algorithm. A tidal-volume-related impedance signal was obtained without calibration (only relative shape matters, as linear fitting, provides the best agreement between the impedance signal and the reference, pneumotachometry [12]). We reduced the sampling frequency from the original 250 Hz to 25 Hz for computational reasons.

## C. Cardiorespiratory analysis

The algorithm was applied to all the signals to find the graphs of temporal orders with the courses of a causal vector. Then, the CVs were analyzed concerning body positions and changes in breathing style, to determine how they affect the relationships.

During the analysis, average and standard deviation values of the estimated causal vector were calculated for supine/standing body positions in Group A, and for different depths and rates of breathing in Group B, using the main settings:

- · with Gaussian standardization;
- · segment length - 10 seconds; and
- · maximal shift in both directions - 2 seconds (based on average breathing frequency and the RSA phenomenon characteristics).

Then, to assess the stability more robustly, we checked the duration of the longest stable parts and the duration of CV to that of the entire signal. If we apply thresholds, 0 90 for . LM and 0 15 for TD approaches will be used arbitrarily. .

## III. RESULTS

Sample temporal order graphs for very regular signals recorded for a supine position are shown in Fig. 1, and for less regular signals obtained during standing - in Fig. 2. In both cases, we applied thresholds.

Fig. 1. Sample TO estimation for very regular signals acquired from Group A participant #1 (supine). For LM: navy blue means 1 (best) and red means 0 (worst); for TD: navy blue means the smallest possible value (best).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00123v2-figure-952732eb-4092-4345-b235-093a41409dc9.jpg)

Table II gathers the means and standard deviations of the causal vector in Group A for both supine and standing body positions without a threshold, for linear modeling, and TD Fourier kernels. The effect of respiratory sinus arrhythmia seems to be present, particularly for the supine body position. Alternative results obtained using considered methods come from different definitions, which may affect their compliance.

Fig. 2. Sample TO estimation for less regular signals acquired from another Group A participant #3 (standing). Colorbars description is the same as for Fig. 1.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00123v2-figure-f2c3137b-9a97-41e2-8e55-a7cc04573a61.jpg)

## TABLE II

SUMMARY OF THE MEAN ± SD OF CAUSAL VECTOR VALUE DEPENDING ON BODY POSITION, FOR ALL PARTICIPANTS IN GROUP A. A MINUS MEANS THAT THE TIDAL VOLUME CURVE IS AHEAD OF THE

TACHOGRAM. ALL VALUES ARE IN MILLISECONDS.

| Position   | Linear modeling   | TD Fourier-type   |
|------------|-------------------|-------------------|
| Supine     | 302 ± 466         | 469 ± 518         |
| Standing   | - 2 ± 244         | 239 ± 794         |

Table III collects the duration of the longest stable CV part and ratio of the duration to that of the signal) after applying a threshold only to the TD-Manhattan approach.

## TABLE III

SUMMARY OF THE MEAN ± STANDARD DEVIATIONS OF THE DURATION OF THE LONGEST STABLE CV PART, AND RATIO OF THE DURATION TO THE DURATION OF THE SIGNAL, FOR ALL PARTICIPANTS IN GROUP A, DEPENDING ON BODY POSITION. 'N' COUNTS, IN HOW MANY CASES

THE VALUE ABOVE THE THRESHOLD WAS FOUND.

| Position   |   N | Longest [sec]   | Ratio [%]       |
|------------|-----|-----------------|-----------------|
| Supine     |   7 | 15 . 1 ± 16 . 8 | 21 . 2 ± 20 . 1 |
| Standing   |   6 | 9 . 0 ± 4 . 0   | 11 . 0 ± 8 . 8  |

When applying arbitrary thresholds, there are relatively few and short fragments of constant causal vectors between signals. The mean separate causal vectors lasted 571 ms and 431 ms for supine and standing body positions, respectively. Next, sample visual results for the entire breathing protocol in Group B are shown in Fig. 3.

Similarly, as earlier, the mean and standard deviations were calculated for the TD-Manhattan approach (without thresholds) for different breathing rates and depths (in Group B) and are shown in Table IV.

The results suggest the rate of breathing has a more significant impact on the mean value of the causal vector than the depth of breathing. The comparison of the average CVs (calculated as a mean value independently of the depth) for Group B, along with the cycle durations, depending on the breathing rate, are stored in Table V.

Interestingly, the RR interval curve precedes the tidal volume more when breathing is slower. It is further discussed in the next section.

Fig. 3. Sample TO estimation for the whole protocol acquired from a Group B participant #7; using only the TD-Manhattan approach with a threshold, with shifts from -2 s to 5 s . For the last part of the signal (15 BPM breathing rate), the second period appears to be present for shifts close to 5 s . Colorbars description is the same as for Fig. 1.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00123v2-figure-46fa9ecb-e7e5-4844-9c2e-8fd97ae72200.jpg)

## TABLE IV

SUMMARY OF MEAN ± STANDARD DEVIATIONS OF CAUSAL VECTOR VALUES, FOR ALL GROUP B, FOR EVERY COMBINATION OF BREATHING DEPTH (NORMAL AND DEEP) AND BREATHING RATE (6, 10, AND 15 BREATHS PER MINUTE, BPM), FOR THE TD MANHATTAN APPROACH,

WITHOUT A THRESHOLD. ALL RESULTS ARE IN MILLISECONDS.

| Depth   | 6 BPM      | 10 BPM    | 15 BPM    |
|---------|------------|-----------|-----------|
| Normal  | 2469 ± 643 | 829 ± 280 | 323 ± 334 |
| Deep    | 2499 ± 380 | 854 ± 363 | 304 ± 317 |

## TABLE V

THE COMPARISON OF AVERAGE CAUSAL VECTOR VALUES, CALCULATED INDEPENDENTLY OF THE BREATHING DEPTH FOR ALL GROUP B, WITH THE CYCLE DURATION AND THE RATIO BETWEEN OPTIMAL SHIFT (AVERAGE CAUSAL VECTOR), DEPENDING ON THE BREATHING RATE (6, 10, AND 15 BREATHS PER MINUTE, BPM), FOR THE TD-MANHATTAN APPROACH, WITHOUT A THRESHOLD.

| Rate   | Optimal shift (OS)   | Cycle duration (CD)   | OS/CD   |
|--------|----------------------|-----------------------|---------|
| 6 BPM  | 2484 ms              | 10 s                  | 24 . 8% |
| 10 BPM | 842 ms               | 6 s                   | 14 . 0% |
| 15 BPM | 314 ms               | 4 s                   | 7 . 9%  |

## IV. DISCUSSION

The causal analysis appears to be a promising tool to extend the classical approaches of cardiorespiratory analysis and enable answers to new physiological questions. In a traditional approach (if considered at all), a graph of connections is rather taken as an input, 'prior information' based on medical knowledge, even if it is relatively hard to be established, specifically for an individual athlete. We would like to emphasize the opposite approach, in which the results of the causal-related analysis might be used as the input. We call this the 'bottom-up' strategy [6], [7].

It can be used by practitioners to (1) identify optimal training schedules to promote desirable progress and achieve optimal performance (based on the individual directions and strengths of cardiorespiratory relationships), or to (b) establish sufficient training loads; all considering interventions applied to the cause 'variables' [6], [7].

Moreover, the study protocols might be arranged in a way enabling testing of adaptation, recovery status, etc. Then, the possible levels of analysis might be supplemented with a local temporal order estimation. In our opinion, the presented concept may be used primarily to test the short temporal stability of relationships (this has already been discussed by Porta et al. [13], in the context of different physiological phenomena), but also to check their directionality and strength. The stability can be described as strong when the causal vector is almost constant over time. Protocol changes, like switching from supine to standing, may be compared regarding explored temporal orders. Directionality, in turn, can be indirectly assessed by the signal shift sign.

The results were that RR interval curve precedes the tidal volume one. It is not contradictory to RSA but it results from taking signals into consideration. The local maximum of RR interval occurs during the inspiratory phase, so based on the shape, it is earlier, still being caused by respiratory activity. Worth noted is that if tidal volume was replaced by airflow signal as derivative, temporal orders would change.

We then also observed (in the control group, where controlled breathing removes effects other than RSA) that the precedence of the tachogram curve before the tidal volume one is relatively greater when the breathing is forced to be slower. However, as the presented analysis is only an illustration, we do not treat it as a general medical conclusion. This issue certainly requires a further look. In our opinion, this may be connected to the slow breathing concept as an intervention that subjects may practice benefiting health (particularly cardiovascular) [14]. Here, the relation between the depth of breathing and causal vector values appeared insignificant; however, to maintain a lower breathing rate without disturbing respiratory homeostasis, a tidal volume must be increased (the relation between temporal and amplitude ventilation coefficients [14]).

We think various improvements may be considered:

- · introducing weights for different depths of shifts;
- · assessing temporal orders not only based on raw data, but also on signals' phases;
- · analyzing the series of sub-maximum values during a stable segment of data to indirectly evaluate the local strength of connections (not only a causal vector);
- · considering more sophisticated approaches;
- · transitioning an estimation to the multivariate case;
- · evaluating not only the signals' shapes but also their local extrema (enabling more to find phase relations and weaker and intermittent effects); or
- · enabling a more robust (perhaps even automatic) choice of the method's parameters, e.g., the thresholds for both LM and TD.

The method could be a useful complement, e.g., for transfer entropy (used in bivariate case); however, the limitations of the presented technique are (1) the relatively substantial dependence on the input parameter settings (in the presented setup they were manually fixed), and (2) the significance of the choice of analyzed signals. They should be carefully adapted or considered based on the problem.

## V. CONCLUSIONS

A technique for exploring temporal orders in physiological data was introduced using the example of cardiorespiratory data (tachogram and tidal volume curves) in two groups (elite athletes and control, made up of students). We proposed two different approaches. The former uses linear modeling, the latter time series distance calculation.

Respiratory sinus arrhythmia effect seems to be the most visible, particularly for the supine body position. The results also suggested that the rate of breathing has a greater impact on the mean value of the causal vector than does the depth of breathing. Interestingly, the RR interval curve precedes the tidal volume more when breathing is slower.

The method can be used as an initial step during causality analysis, as it may show the temporal relationships between signals, the stability of the temporal orders over time, highlighting the similar parts of the causal vector curve (which can be then compared with events in the study protocol).

## REFERENCES

- [1] R. P. Bartsch, K. K. L. Liu, A. Bashan and P. C. Ivanov, Network physiology: How organ systems dynamically interact, PLOS ONE, vol. 10, pp. 1-36, 2015.
- [2] J. Fatisson, V. Oswald and F. Lalonde, Influence diagram of physiological and environmental factors affecting heart rate variability: an extended literature overview, Heart International, vol. 11, pp. 5000232, 2016.
- [3] P. Larsen, Y. Tzeng, P. Sin and D. Galletly, Respiratory sinus arrhythmia in conscious humans during spontaneous respiration, Respiratory Physiology and Neurobiology, vol. 174, pp. 111-118, 2010.
- [4] G. A. Reyes del Paso, W. Langewitz, L. J. Mulder, A. Van Roon and S. Duschek, The utility of low frequency heart rate variability as an index of sympathetic cardiac tone: a review with emphasis on a reanalysis of previous studies, Psychophysiology, vol. 50, pp. 477-487, 2013.
- [5] T. Penzel, J. W. Kantelhardt, R. P. Bartsch, M. Riedl, J. F. Kraemer, N. Wessel, C. Garcia, M. Glos, I. Fietze and C. Schobel, Modulations of heart rate, ecg, and cardio-respiratory coupling observed in polysomnography, Frontiers in Physiology vol. 7, no. 460, 2016.
- [6] M. Mły´ nczak and H. Krysztofiak, Discovery of causal paths in cardiorespiratory parameters: a time-independent approach in elite athletes, Frontiers in Physiology, vol. 9, no. 1455, 2018.
- [7] M. Mły´ nczak and H. Krysztofiak, Cardiorespiratory temporal causal links and the differences by sport or lack thereof, Frontiers in Physiology, col. 10, no. 45, 2019.
- [8] L. Schiatti, G. Nollo, G. Rossato and L. Faes, Extended granger causality: a new tool to identify the structure of physiological networks, Physiological Measurement vol. 36, pp. 827, 2015.
- [9] J. Peters, D. Janzing and B. Scholkopf, Causal inference on time series using restricted structural equation models, In Advances in Neural Information Processing Systems, pp. 154-162, 2013.
- [10] http://marcelmlynczak.com/#software
- [11] M. Mły´ nczak, W. Niewiadomski, M. Zyli´ nski and G. Cybulski, Am-˙ bulatory devices measuring cardiorespiratory activity with motion, in Proceedings of the 10th International Joint Conference on Biomedical Engineering Systems and Technologies (Porto: SciTePress), pp. 91-97, 2017.
- [12] M. Mły´ nczak, W. Niewiadomski, M. Zyli´ nski and G. Cybulski, Assess-˙ ment of calibration methods on impedance pneumography accuracy, Biomedical Engineering/Biomedizinische Technik, vol. 61, no. 6, pp. 587-593, 2015.
- [13] A. Porta, V. Bari, B. De Maria, B. Cairo, E. Vaini, N.M. Perseguini, J. Milan-Mattos, P. Rehder-Santos, V. Minatel, A.C.M. Takahashi and A.M. Catai, Comparison between probabilistic and Wiener-Granger causality in assessing modifications of the cardiac baroreflex control with age, Physiological Measurement, vol. 39, no. 10, 104004, pp. 1-10, 2018.
- [14] M.A. Russo, D.M. Santarelli and D. O'Rourke, The physiological effects of slow breathing in the healthy human, Breathe, vol. 13, no. 4, pp. 298-309, 2017.
Marcel Młyńczak
Temporal orders and causal vector for physiological data analysis
In addition to the global parameter- and time-series-based approaches,
physiological analyses should constitute a local temporal one, particularly
when analyzing data within protocol segments. Hence, we introduce the R package
implementing the estimation of temporal orders with a causal vector (CV). It
may use linear modeling or time series distance. The algorithm was tested on
cardiorespiratory data comprising tidal volume and tachogram curves, obtained
from elite athletes (supine and standing, in static conditions) and a control
group (different rates and depths of breathing, while supine). We checked the
relation between CV and body position or breathing style. The rate of breathing
had a greater impact on the CV than does the depth. The tachogram curve
preceded the tidal volume relatively more when breathing was slower. Clinical
relevance - The method can assess (1) relationships between two signals, having
one of the two time-shifted in a given range, (2) curves of most optimal
inter-signal shift, and (c) their stability across time; also for other
physiological studies.