**Alterations of Cardiovascular Complexity during Acute Exposure to High Altitude: A Multiscale Entropy Approach**

#### **Andrea Faini 1, Sergio Caravita 1,2, Gianfranco Parati 1,3,\* and Paolo Castiglioni <sup>4</sup>**


Received: 31 October 2019; Accepted: 10 December 2019; Published: 15 December 2019

**Abstract:** Stays at high altitude induce alterations in cardiovascular control and are a model of specific pathological cardiovascular derangements at sea level. However, high-altitude alterations of the complex cardiovascular dynamics remain an almost unexplored issue. Therefore, our aim is to describe the altered cardiovascular complexity at high altitude with a multiscale entropy (*MSE*) approach. We recorded the beat-by-beat series of systolic and diastolic blood pressure and heart rate in 20 participants for 15 min twice, at sea level and after arrival at 4554 m a.s.l. We estimated Sample Entropy and *MSE* at scales of up to 64 beats, deriving average *MSE* values over the scales corresponding to the high-frequency (*MSE*HF) and low-frequency (*MSE*LF) bands of heart-rate variability. We found a significant loss of complexity at heart-rate and blood-pressure scales complementary to each other, with the decrease with high altitude being concentrated at Sample Entropy and at *MSE*HF for heart rate and at *MSE*LF for blood pressure. These changes can be ascribed to the acutely increased chemoreflex sensitivity in hypoxia that causes sympathetic activation and hyperventilation. Considering high altitude as a model of pathological states like heart failure, our results suggest new ways for monitoring treatments and rehabilitation protocols.

**Keywords:** Sampen; cross-entropy; autonomic nervous system; heart rate; blood pressure; hypobaric hypoxia; rehabilitation medicine

#### **1. Introduction**

There is an increasing interest in the physiological adaptations that occur during exposures to high-altitude conditions, particularly in the alterations of autonomic cardiovascular control. This is due to the extraordinary development of mountain tourism that leads millions of people each year to stay in the high mountains for short periods of time. In addition to mountain tourism, millions of other people live permanently at high altitudes and are exposed to conditions that may cause episodes of mountain sickness [1]. The alterations to cardiovascular control caused by high altitudes are mainly due to hypobaric hypoxia, i.e., the low partial pressure of oxygen in the air, which produces an autonomic response by increasing the chemosensitivity, possibly altering the overall integrative autonomic regulation.

Interestingly, some of the cardiovascular changes that can be observed during exposure to high altitudes are similar to the autonomic alterations occurring in some diseased conditions, as in patients with heart failure, stroke or metabolic disorders [2,3]. Therefore, stays at high altitude can be viewed as an experimental model of some pathological conditions that affect autonomic cardiovascular regulation at sea level. In this regard, the study of cardiovascular control at high altitudes may help to better understand some pathophysiological mechanisms and may be beneficial for improving treatments and outcomes in rehabilitation medicine.

Most of the high-altitude studies in the literature describe the cardiovascular autonomic alterations with linear measures of heart rate variability, in relatively small groups of participants. By contrast, the literature reports very few nonlinear measures, which exclusively regard heart rate variability during sleep [4,5], when apneas induced by hypoxia may profoundly alter the heart rate dynamics. Furthermore, very few studies evaluate high-altitude alterations in beat-by-beat blood pressure variability [6,7] because of the logistical difficulties in performing such physiological recordings at high altitude. Therefore, the effect of a stay at high altitude on the complex dynamics of the cardiovascular system in the waking state remains an unexplored issue.

Our work contributes to filling this gap by assessing changes in cardiovascular complexity during a short-term stay at high altitude (4554 m a.s.l.). One of the most effective tools for extracting information on the complex dynamics of physiological systems is multiscale entropy, and our work is based on the multivariate and multiscale assessment of entropy and cross-entropy on beat-by-beat recordings of heart rate and arterial blood pressure in healthy volunteers. Our aim is to identify those aspects of the complex dynamics of the cardiovascular system that better describe the autonomic changes in response to the chemoreflex activation expected to occur in the waking state during exposure to hypobaric hypoxia. Given that high-altitude cardiovascular alterations in healthy individuals may be a model of some pathological conditions at sea level, the results of our study may indicate novel ways for monitoring the severity of a deteriorated autonomic cardiovascular control and the efficacy of treatment during rehabilitation programs.

#### **2. Materials and Methods**

#### *2.1. Subjects and Data Collection*

This study is based on data collected previously to evaluate the effectiveness of a drug (acetazolamide) for treating mountain sickness during acute exposure to high-altitude hypoxia [8]. In the present work, we consider the 20 Caucasian volunteers of the placebo group, who completed the hemodynamic recordings at sea level and at high altitude. The placebo group was composed of 10 males and 10 females with mean (SD) age of 37 (9.5) years old and body mass index of 22.3 (2.7) kg/m2. They were healthy lowlanders without known cardiovascular disease or chronic cardiovascular therapy, without a history of severe mountain sickness. None of them practiced professional sports, all lived in Milan or its surroundings, and had no recent exposure to altitudes above 2000 m. They took the placebo orally in pill form twice a day.

In all participants, the cardiovascular measurements were taken twice. The first recording session was performed at baseline, in the normobaric/normoxyc conditions of our laboratory in Milan (122 m a.s.l.). The second recording session was performed in the hypobaric/hypoxic conditions of the Margherita Hut (4554 m a.s.l. on Monte Rosa in the Italian Alps). The ascent from Milan to the Margherita Hut was completed in about 28 h, by car and cable car up to 3200 m a.s.l. and then by foot, spending one night for acclimatization in the Gniffetti hut (3647 m a.s.l.). Recordings at high altitude were performed 2 days after the arrival at the Margherita Hut.

The measurements were performed in a quiet environment with ambient temperature of about 20 ◦C. They consisted of the simultaneous recording of one-lead electrocardiogram, ECG (PowerLab, ADInstruments, Sydney, Australia at sea level; ECG100C, MP150 Biopac Systems, Santa Barbara, CA at high altitude) and of continuous arterial blood pressure at the finger artery (Nexfin, BMEYE, Amsterdam, The Netherlands at sea level; Portapres, Finapres Medical Systems, The Netherlands at high altitude) for 15 min. Brachial arterial blood pressure waveforms were reconstructed from the measured finger blood pressure waveforms through the transfer function method [9].

Hemoglobin oxygen saturation (NPB-295, Nellcor Puritan Bennett Inc., Plaseanton, CA, USA) and the respiratory activity with a spirometry device (Vmax SensorMedics 2200, Yorba Linda, CA, USA) were also measured during the recordings. Each subject rested in a semi-recumbent position and the recordings started after an adaptation period to the new posture of at least 5 min, ensuring that the participants felt comfortable with the setup and had no apparent urges that could influence their responses. A familiarization recording session was performed several days before the first session, at sea level, which allowed the participants to become accustomed to the experimental setup (the signals acquired during the familiarization sessions are not considered in this study). At high altitude, the Lake Louise Score Questionnaires was administered to evaluate the presence (score ≥ 5) or absence of acute mountain sickness [10].

The Ethical Committee of the Istituto Auxologico Italiano (Milan, Italy) approved the study protocol (procedure number CE Auxologico: 2010\_04\_13\_01); all the participants gave their written informed consent to the study procedures.

#### *2.2. Data Preprocessing and Spectral Analysis*

Recordings of finger blood pressure (BP) and of the ECG were digitized at 200 Hz and 12 bits. Each R peak of the ECG was identified by a derivative-and-threshold algorithm and a parabolic interpolation was used to refine the R wave fiducial point as suggested in [11]. The interval between consecutive R peaks, i.e., the R-R Interval (RRI), was calculated beat by beat for cardiac beats resulting from sinus node depolarization. The systolic (S) BP and diastolic (D) BP values were identified beat by beat from the finger BP signal. SBP and DBP values associated with premature beats (as identified from the ECG) or with calibrations of the device for recording finger BP were removed. The percentage of removed RRI beats was 0.5% on average, with 5% being the worst case; with respect to SBP/DBP, the percentage of removed beats increased to 2.8% on average, with 10% the worst case. The Pulse Interval (PI) was calculated beat by beat as the time interval between consecutive systolic peaks. The duration of each breathing cycle was identified on the respiratory signals as the interval between the start of consecutive inspiratory phases. The beat-by-beat series of RRI, SBP, and DBP were interpolated at 5 Hz before spectral analysis to obtain evenly sampled series and to linearly interpolate missing beats. Power spectra were calculated by the Welch periodogram with 80% overlapped Hann data windows of 240 s length and by integrating the periodogram over the very-low frequency (VLF, between 0.003 and 0.04 Hz), the low frequency (LF, between 0.04 and 0.15 Hz) and the high-frequency (HF, between 0.15 and 0.4 Hz) bands, as defined in international guidelines [11].

#### *2.3. Multiscale Entropy of RRI, SBP, and DBP*

We estimated the multiscale entropy (*MSE*) of RRI, SBP, and DBP with the method proposed in [12]. The method is based on the original approach to evaluate multiscale entropy as Sample Entropy, *SampEn*, of progressively coarse-grained sub-series [13] (where the coarse graining is obtained by decimation, taking one sample every *n* after a moving average of order *n*) with the same tolerance threshold at each coarse-graining order [14]. However, this includes the subsequent suggestion to replace the decimation with a delay *n* that increases the statistical consistency of the estimate [15]. Furthermore, our method substitutes the moving average with a Butterworth filter to improve the scale resolution, as shown in the case of decimation by others [16].

Briefly, given a time series of *N* samples **X** = {*x1,x2,* ... *,xN*}, let's call **Y***<sup>n</sup>* = {*y1 n,y2 n,* ... *,yN<sup>n</sup>*} the output of the zero-phase, 6th-order low-pass Butterworth filter with cut-off frequency *fc* = 0.5/*n* applied to **X**. The template vectors at a given embedding dimension *m* and scale *n* are built considering a delay of *n* samples between consecutive elements:

$$\mathbf{y}\_{i}^{m}(n) = \left[ y\_{i}^{n}, y\_{i+n}^{n}, \dots, y\_{i+(m-1)n}^{n} \right]^{T}, 1 \le I \le N-mn \tag{1}$$

The infinity norm distance between any couple of template vectors is

$$d\_{ij}^{m}(n) = \left\| \mathbf{y}\_{i}^{m}(n) - \mathbf{y}\_{j}^{m}(n) \right\|\_{\infty'} \ 1 \le i, j \le N - mn, \ j > I + n \tag{2}$$

and the number of "paired-vectors" *np*(*m,n,r*), which are pairs of vectors with distance lower than a predefined tolerance threshold *r*, is calculated based on the infinity norm. Repeating these steps for the successive dimension *m* + *1*, the sample entropy of **Y***<sup>n</sup>* with delay *n* and tolerance *r* is:

$$\text{SampEn}(\mathbf{Y}^n, \mathbf{N}, m, n, r) = -\ln \frac{n\_p(m+1, n, r)}{n\_p(m, n, r)}\tag{3}$$

This leads to the definition of the *MSE* of **X**, which is a function of the scale *n*, at a given embedding dimension *m* and tolerance *r*, as

$$MSE(n) = SampleEn(\mathbf{Y}^{\text{dl}}, \mathbf{N}, m, n, r) \tag{4}$$

At *n* = 1, **Y**<sup>1</sup> = **X** and *MSE(1)* coincides with the *SampEn* of **X**.

For each RRI, SBP and DBP series, we calculated *MSE*(*n*) for 1 ≤ *n* ≤ 64 beats. The tolerance threshold is commonly set at *r* = 0.20 times the standard deviation of the time series in heart rate variability studies and we also adopted this choice. As to the embedding dimension, in addition to *m* = 2, traditionally used in *SampEn* analysis of heart rate variability, we also considered *m* = 1, because we previously showed that it provides *MSE(n)* profiles similar to *m* = 2 but with better statistical consistency [14]. To compare sea-level and high-altitude conditions over the same temporal scales, in seconds, we mapped the scale units from number of beats, *n*, to time *t*, in seconds, with the transformation:

$$t = n \times \text{} \tag{5}$$

where <*RRI*> is the mean RRI of each series, in seconds. We interpolated and resampled *MSE* to obtain 100 estimates at scales *t* exponentially distributed over the time axis between 1 s and 48 s. This range includes the scales associated with the traditional high-frequency (HF, 2.5 ≤ *t*< 6.7 s) and low-frequency (LF, 6.7 ≤ *t* < 25 s) bands of the heart rate variability spectra. As a concise way to represent the results, the *MSE* functions were averaged over the scales included in the HF and LF bands, obtaining the *MSEHF* and *MSELF* indices.

To evaluate the performance of our *MSE* estimator, we synthetized 10 series of white noise and 10 series of pink noise, each of *N* = 1000 samples. This length corresponds to 15' recordings at the average RRI of 900 ms. Then we calculated *MSE* over the scales associated with the *HF* and *LF* bands. The estimates in Figure 1 demonstrate the capability of our MSE method to faithfully describe the entropy structures of these random processes over the scales corresponding to the HF and LF bands.

**Figure 1.** (**a**) Profiles of Multiscale Entropy MSE calculated with embedding dimension *m* = 1: mean ± SD for 10 synthesized series of white noise and 10 synthesized series of pink noise, each of 1000 samples, simulating 15' beat-by-beat recordings with mean RRI equal to 900 ms. Gray bands show the ranges of scales corresponding to the HF and LF spectral bands. (**b**) MSE calculated with *m* = 2 for the same data of panel (**a**).

#### *2.4. Multiscale Cross-Entropy between SBP and PI*

To estimate the cross-entropy between blood pressure and heart rate, we used PI rather than RRI values to more easily couple blood pressure and heart rate series beat by beat (the number of valid beats of the RRI series may differ from those of SBP and PI due to the presence of calibration periods in the device measuring the finger arterial pressure). The multiscale cross-entropy, *XMSE*, between the SBP and PI series was defined extending the cross-sample entropy estimator, *XSampEn* [17], to multiple scales similarly to the way we defined *MSE* extending *SampEn* at multiple scales. However, the PI and SBP series are normalized to unit variance and zero mean before applying the Butterworth filters with cut-off frequency *fc* = 0.5/*n* to obtain the **P***<sup>n</sup>* = {*p1 n,p2 <sup>n</sup>* ... *,pN<sup>n</sup>*} and **S***n*={*s1 n,s2 <sup>n</sup>* ... *,sN<sup>n</sup>*} output series. The template vectors for the embedding dimension *m* at scale *n* are

$$\begin{aligned} \mathbf{p}\_{i}^{m}(n) &= \begin{bmatrix} p\_{i}^{n}, p\_{i+n'}^{n} \dots p\_{i+(m-1)n}^{n} \end{bmatrix}\_{\Gamma}^{T} \quad , \ 1 \le I \le N-mn \tag{6} \\\ \mathbf{s}\_{i}^{m}(n) &= \begin{bmatrix} s\_{i}^{n}, s\_{i+n'}^{n} \dots s\_{i+(m-1)n}^{n} \end{bmatrix}\_{\Gamma} \end{aligned} \tag{7}$$

Based on the distances between couples of vectors

$$d\_{ij}^{m}(n) = \left\| \mathbf{p}\_{i}^{m}(n) - \mathbf{s}\_{j}^{m}(n) \right\|\_{\text{op}'} \ 1 \le ij \le N - mn \tag{7}$$

the number of paired vectors with distance lower than a threshold *r*, *npx*(*m,n,r*), is calculated. Repeating these steps for *m* + *1*, the cross-*SampEn* between **P***<sup>n</sup>* and **S***<sup>n</sup>* with delay *n* is:

$$XSampleEn(\mathbf{P''}, \mathbf{S''}, N, m, n, r) = -\ln \frac{n\_{\text{px}}(m+1, n, r)}{n\_{\text{px}}(m, n, r)}\tag{8}$$

and the multiscale cross-entropy between SBP and PI is estimated as:

$$XMSE(n) = XSampleEn(\mathbf{P}^{\prime\prime}, \mathbf{S}^{\prime\prime}, \mathcal{N}\_r m\_r n\_r r) \tag{9}$$

*XMSE* between SBP and PI was assessed at the embedding dimensions *m* = 1 and *m* = 2, with *r* = 0.20. Clearly, at *n* = 1 *XMSE* coincides with the cross-SampEn. The scales were mapped from beats, *n*, into times, τ, according to Equation (5), then interpolated and resampled in a similar fashion between 1 and 48 s. Finally, the *XMSE* functions were averaged over the scales included in the HF and LF bands, obtaining the *XMSEHF* and *XMSELF* indices.

#### *2.5. Statistical Analysis*

All indices were compared between sea level and high altitude with paired t-tests whenever their distribution passed the Shapiro-Wilks Gaussianity test at p = 0.05, possibly after a Box-Cox transformation [18]. Otherwise, they were compared by the non-parametric Wilcoxon signed-rank test. As regards spectral analysis, powers were log-transformed and we considered only participants with an average breathing period shorter than 7 s because the HF band correctly reflects respiratory-driven modulations only in this case; for this reason, we discarded 5 participants from the statistical tests on spectral powers.

Furthermore, we calculated the Wilcoxon signed-rank test statistics of multiscale entropy at each scale separately to easily identify the scales better reflecting the alterations induced by high altitude in the *MSE*(*t*) and *XMSE*(*t*) profiles. With respect to the derived entropy indices (namely *SampEn*, *MSEHF,* and *MSELF* of RRI, SBP, and DBP; *XSampEn*, *XMSEHF*, and *XMSELF* between SBP and PI) we also calculated the 95% confidence intervals of the mean for the difference between high-altitude and sea-level conditions. The confidence intervals were obtained with the nonparametric bootstrap method by randomly sampling with replacement the original 20 measures 1000 times; this bootstrapping allows high-quality estimates of the confidence intervals to be obtained with a distribution-free approach that avoids making any assumption on the nature of the distributions.

The threshold for statistical significance was set at 5%. with a two-sided alternative hypothesis. All the statistical tests were performed with "R: A Language and Environment for Statistical Computing" software package (R Core Team, R Foundation for Statistical Computing, Vienna, Austria, 2018).

#### **3. Results**

The high-altitude conditions were characterized by faster and deeper breathing, by lower hemoglobin oxygen saturation, and by higher blood pressure and heart rate levels (Table 1). The spectral powers of RRI changed at high altitude, with decreased VLF, LF and HF powers and increased LF/HF powers ratio. By contrast, the high altitude did not change the LF and HF powers of SBP and DBP and only marginally decreased their VLF power (Table 1). Half of the participants (N = 10, 6 males) presented acute mountain sickness.


**Table 1.** General cardiorespiratory characteristics and spectral indices: mean (SD) at sea level and high altitude.

\* and \*\* indicate differences at 5% and 1% significance; 5 of the 20 participants are discarded from the statistics on spectral powers because their average breathing rate felt below the HF band.

#### *3.1. Multiscale Entropy*

Figure 2 compares the profiles of multiscale entropy. At sea level, the *MSE(t)* profile of RRI decreases with *t* from the maximum at 2 s, reaching a plateau at scales greater than 7 s. The high-altitude condition reduces *MSE* at the shorter end only, and thus the *MSE(t)* profile is flatter at high altitude. By contrast, SBP and DBP have higher *MSE* values at scales within the HF band and this pattern is more pronounced at high altitude due to the substantial reduction of *MSE* at scales within the LF band; in particular, at scales greater than 10 s, the entropy reduction is more significant for *m* = 1.

*SampEn*, which corresponds to *MSE* at the scale of 1 beat, decreases at high altitude for RRI only (Table 2), without changes for SBP or DBP. The decrease is more significant for *m* = 2.


**Table 2.** *SampEn*: mean (SD) over the group at sea level and high altitude.

*m* = embedding dimension; \*\* indicates differences at 1% significance.

**Figure 2.** (**a**) Profiles of Multiscale Entropy MSE at sea level (SL, blue lines) and high altitude (HA, red lines) for RRI calculated with embedding dimension *m* = 1: mean ± sem on 20 participants (gray bands show the ranges of scales corresponding to the HF and LF spectral bands); (**b**) MSE calculated as in (**a**) for SBP; (**c**) MSE calculated as in (**a**) for DBP; (**d**) MSE for RRI calculated as in (**a**) but with *m* = 2; (**e**) MSE calculated as in (**d**) for SBP; (**f**) MSE calculated as in (**d**) for DBP; (**g**) Wilcoxon signed-rank statistics *V* for the comparison between SL and HA for MSE of RRI; the red horizontal lines are the 5% (continuous) or 1% (dashed) percentiles of the distribution for the null hypothesis: when *V* is above these thresholds the hypothesis of similar entropies can be rejected at the corresponding significance level; (**h**) *V* statistics for the comparison between SL and HA for SBP MSE; (**i**) *V* statistics for the comparison between SL and HA for DBP MSE.

Similarly, we found a significant reduction of *MSEHF* for RRI and not for SBP or DBP; by contrast, *MSELF* decreased at high altitude for SBP and DBP but not for RRI, in these cases with greater statistical significance for *m* = 1 (Table 3).

The 95% confidence intervals of the differences between high altitude and sea level (Figure 3) confirm that exposure to high altitude reduces *SampEn* and *MSEHF* of RRI and reduces *MSELF* of SBP and DBP. The greater statistical power of the bootstrap method in Figure 3 compared to the paired *t*-test in Table 2 is reflected in the 95% confidence intervals of the *SampEn* variations which do not cross the zero line also for *m* = *1*.


**Table 3.** Multiscale entropy over the HF and LF band: mean (SD) at sea level and at high altitude.

*m* = embedding dimension; \* and \*\* indicate differences at 5% and 1% significance.

**Figure 3.** 95% confidence intervals of the difference between high-altitude and sea-level conditions of entropy indices. *From top to bottom*: Sample Entropy (*SampEn*), multiscale entropy over the HF (*MSEHF*) and over the LF (*MSELF*) bands; *m* is the embedding dimension.

#### *3.2. Multiscale Cross-Entropy*

The scale-by-scale profiles of *XMSE* between SBP and PI (Figure 4) show the highest values at scales within the HF band. The high altitude decreases cross-entropy at the shorter scales, the reduction being significant at 3–4 s.

However, the average decrease of SBP-PI cross-entropy in the HF band, *XMSEHF,* does not reach the statistical significance; furthermore, *XSampen* and *XMSELF* appear substantially similar in the two conditions (Table 4 and Figure 5).

**Figure 4.** (**a**) Profiles of multiscale cross-entropy XMSE between SBP and PI at sea level (SL, blue lines) and high altitude (HA, red lines): mean ± sem for embedding dimension *m* = 1; (**b**) XMSE from the same data of panel (**a**) calculated for embedding dimension *m* = 2; (**c**) Wilcoxon signed-rank statistics *V* for the comparison between SL and HA; the red horizontal lines are the 5% (continuous) or 1% (dashed) percentiles of the distribution for the null hypothesis: when *V* is above these thresholds, the hypothesis of similar entropies can be rejected at the corresponding significance level. Gray bands show the ranges of scales corresponding to the HF and LF spectral bands.

**Table 4.** SBP-PI cross-entropy indices: mean (SD) at sea level and at high altitude.


*m* = embedding dimension.

**Figure 5.** 95% confidence intervals of the difference between high-altitude and sea-level conditions of cross-entropy indices. *From top to bottom*: Cross sample entropy (*XSampEn*), cross multiscale entropy over the HF (*XMSEHF*) and the LF (*XMSELF*) bands; *m* is the embedding dimension.

#### **4. Discussion**

We described the effects of a short stay at high altitude on the complex dynamics of the cardiovascular system. Novelties of our study are that for the first time, the effects of high altitude are described with a multiscale entropy approach, and that the cardiovascular dynamics is evaluated considering not only the heart rate variability, as is usually done in studies on the autonomic cardiovascular control at high altitude, but also the beat-by-beat variability of arterial blood pressure. It should be noted that traditional spectral analysis and complexity analysis (when assessed as multiscale entropy) investigate very different aspects of the time series dynamics, with one being related to the amplitude of the fluctuations, and the other to their unpredictability/irregularity. The entropy approach (as originally proposed by Costa et al. [13,19] and which we follow even if using a different, statistically more consistent, estimator for short series) assumes that a proper decomposition of the entropic measure of irregularity at different temporal scales reveals the capability of dynamical systems to adapt to external perturbations and to environmental changes. Following this proposal, multiscale entropy has been assessed in several studies aimed at quantifying the loss of complexity in diseased states, as a way to assess the reduced adaptive capacity of the individual. Similarly, we applied multiscale entropy to quantify the degraded adaptive capacities of the cardiovascular system exposed to high altitude conditions. Interestingly, our results represent an example of the different information provided by spectral analysis and by entropy-based complexity analysis.

#### *4.1. Cardiorespiratoy Variables and Spectral Powers at High Altitude*

Our participants showed the marked decrease of hemoglobin oxygen saturation expected at such a high altitude and the hyperventilation triggered by the chemoreflex control of breathing in response to hypoxia [20]. Acute hypoxia also stimulates peripheral chemoreceptors, producing a sympathetic activation that increases blood pressure [21,22], explaining the increased blood pressure levels we observed at high altitude. The sympathetic activation could also explain the higher heart rate and LF/HF spectral measures, index of cardiac sympatho/vagal balance; furthermore, it could have induced a general decrease of the vagal modulations of heart rate, as quantified by the reduced HF and LF powers. These spectral changes are in line with those repeatedly reported on the acute autonomic effects after an ascent at high altitude [23–26] or at a simulated altitude of 3600 m asl in a hypobaric/hypoxic chamber [27,28]. By contrast, an increased LF power without HF power changes [5] was reported at 3600 m asl; however, differently from our study and from [23], recordings were performed during nighttime sleep when periodic breathing and apneas are likely to occur, making the HF power an unreliable index of the vagal respiratory modulations of the heart rate.

#### *4.2. Heart Rate Multiscale Entropy*

The multiscale entropy of heart rate is a measure of the complexity of the cardiovascular system [13]. The shortest possible scale at which the multiscale entropy can be calculated is the single beat, and in this case, *MSE* coincides with SampEn. Maneuvers eliciting the sympatho/vagal balance, like posture changes or pharmacological blockade [29,30], decrease the heart-rate SampEn, and thus the significant reduction of RRI SampEn at high altitude (Figure 3) could be a consequence of the sympathetic activation induced by the acute exposition to high-altitude hypoxia. Coherently with our result, a decrease of the heart-rate SampEn was observed during simulated high altitude [27], as well as in a real high-altitude environment [5,26]. The opposite trend was reported in participants suddenly exposed to hypoxia simulating the altitude of 8230 m asl in a hypobaric chamber [31], a result that could reflect an abrupt activation of a defensive autonomic response caused by the sudden exposure to such an extreme condition.

However, the novel results of our study are the *MSE* alterations at scales greater than 1 beat. We found a significant decrease in *MSE* over the scales of the HF band that does not extend to larger scales (Figures 2 and 3, Table 2). This suggests a loss of cardiovascular complexity that mainly affects the faster components, probably associated with ventilation, while the cardiac complexity at longer scales is preserved. For some aspects, like the sympathetic and ventilatory responses to hypercapnia [32], the high-altitude condition may represent a model of heart failure and it is worth noting that heart failure patients, compared to healthy subjects, have lower *MSE* values over a broad range of scales that includes the HF band [33].

#### *4.3. Blood Pressure Multiscale Entropy*

Another novel finding regards the alterations in the *MSE* of blood pressure. While exposure to high altitude reduces the *MSE* of RRI in the HF band and at shorter scales, it reduces the *MSELF* of SBP and DBP without affecting their SampEn or *MSEHF*. The lack of alterations in the faster components of blood pressure complexity could be related to the non-autonomic nature of the blood pressure dynamics at scales faster than the LF band. For instance, HF modulations of blood pressure are mainly due to the direct action of the respiratory mechanics, and not to autonomic modulations mediated by chemo- or baro-reflexes, as for RRI. Interestingly, an autonomic influence on blood pressure is expected over a range of larger scales that includes the LF band. At these scales, the sympathetic outflow that reaches the individual vascular districts modulates the arteriolar resistances in order to regulate the local supplies of blood. It could be possible that the high-altitude hypoxia induced overall sympathetic vasoconstrictions to make more oxygen available to the brain and the heart, and that the vasoconstriction substantially decreased the amplitude of local vasomodulations. Therefore, the observed loss of blood pressure complexity in the LF band might reflect an altered sympathetic control of local vascular districts.

#### *4.4. Blood Pressure-Heart Rate Multiscale Cross-Entropy*

We estimated cross-entropy with the *XMSE* estimator to evaluate the degree of asynchronicity between blood pressure and heart rate. We used PI values rather than RRI values to more easily couple the blood pressure and heart rate beat by beat. *XMSE* is based on the conditional probability that SBP-PI pairs of segments that are similar when observed over *m* beats remain similar when the segments are increased by one beat (the higher the probability, the lower the cross-entropy). *XMSE* can therefore provide a more general assessment of the synchronization between time series than other analysis tools, like the squared coherency spectrum, which may reflect the linear components only of the coupling between time series.

Even if *MSE* decreased substantially for RRI at the faster scales and for SBP at the lower scales, the *XMSE* between SBP and PI showed only a marginal decrease at scales around 5 s with a non-significant reduction of *XMSEHF*. This would indicate that the level of synchronization between the two series is preserved, if not even slightly increased in the HF band. A possible explanation for this trend is related to the mechanism that couples SBP with PI in the HF band, i.e., respiration. Each inspiratory phase increases the filling of the left ventricle and thus the stroke volume of the following beat, which in turn increases SBP. The increase in SBP is sensed by the baroreceptors and triggers a vagal baroreflex response that lengthens the following cardiac interval. In this way, a respiratory oscillation is mechanically generated in SBP and coupled to a baroreflex-mediated oscillation in PI, with the same frequency. The minute ventilation increased dramatically at high altitude (Table 1), and thus it may be responsible for SBP and PI coupled oscillations with larger amplitude, and thus for the increased SBP-PI synchronization in the HF band.

#### *4.5. Limitations and Conclusion*

Cardiorespiratory control adapts differently in males and females to short-term stays at high altitude [2,3], and although a recent study did not report an interaction with sex in the effect of high altitude on spectral indices of heart rate variability [34], it cannot be excluded that the alterations of multiscale entropy we described are gender dependent. Our participants were matched by sex (10 males and 10 females), so our results are not biased by the gender composition. However, a larger population is needed to stratify our results by sex.

Another factor possibly influencing the cardiac autonomic adaptation to high altitude is the presence of acute mountain sickness [23]. Our group was composed of 10 participants with acute mountain sickness (age 37.9 ± 8.9 years old) and 10 without acute mountain sickness (age 35.6 ± 9.6 years old), and the results should reflect those of a general population ascending to similar high altitudes [35]. However, a larger group of participants is required to evaluate the possible influence of acute mountain sickness on cardiovascular complexity.

Genetic factors [36,37] and acclimatization [38] may play a role in the capability of the autonomic cardiovascular control to adapt to high-altitude environments. Since all our participants were lowlanders belonging to the same Caucasian ethnicity, and with no recent exposure to high altitude, our study was not designed to investigate the possible influence of genetic factors or acclimatization. Larger groups are needed to evaluate whether similar alterations in multiscale entropy may also characterize different ethnicities or highlander populations. Cardiorespiratory diseases may produce alterations in cardiovascular control even during short-term exposure to moderate altitudes [39]. Since we included only healthy volunteers without known cardiovascular diseases or chronic therapies, future studies are needed to evaluate whether the alterations we observed may be exacerbated by diseased conditions.

We had to use different measuring devices at sea level and at high altitude due to organizational reasons, and in theory, this might have influenced our results. However, since sampling rates, digital resolution, and analog preprocessing filters were the same, we can exclude differences in the quality of the ECG recordings. As to the finger blood pressure, the measuring device at high altitude (Portapres, Finapres Medical Systems, The Netherlands), although specifically designed for portability, is based on the same physical principles and technologies of the laboratory device used at sea level (Nexfin, BMEYE, Amsterdam, The Netherlands). Therefore, we can exclude differences due to the quality of the measuring devices also for the BP measures.

In conclusion, we assessed the alterations induced by a short stay at high altitude in the complexity of the cardiovascular system with a multiscale entropy approach. The alterations indicate a loss of complexity at specific ranges of scales that differ between heart rate and blood pressure and are complementary to each other, being the complexity loss concentrated at the shorter scales for heart rate and at the longer scales for blood pressure. The changes can be ascribed to the increased chemoreflex sensitivity in hypoxia that causes sympathetic activation and hyperventilation. These results may contribute to understanding the physiological adaptations to high altitude; furthermore, considering high-altitude conditions as a model of pathological states like heart failure, they may also help to better understand the loss of cardiovascular complexity in patients, possibly suggesting effective ways to improve treatments or to monitor rehabilitation protocols.

**Author Contributions:** A.F., S.C., and G.P. designed the work. G.P. organized the experiments, A.F. and S.C. collected the data. A.F. and P.C. created the software, analyzed the recordings and interpreted the results. P.C. led the writing of the manuscript. A.F. made the figures. All the authors contributed to the interpretation of the data, drafted the work and revised it.

**Funding:** This research was supported by the IRCCS Istituto Auxologico Italiano and by the Italian Ministry of Health.

**Acknowledgments:** We express our gratitude to all HIGHCARE (High altitude cardiovascular research) investigators.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Hainsworth, R.; Drinkhill, M.J.; Rivera-Chira, M. The autonomic nervous system at high altitude. *Clin. Auton. Res.* **2007**, *17*, 13–19. [CrossRef] [PubMed]


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Zipf's Law of Vasovagal Heart Rate Variability Sequences**

#### **Jacques-Olivier Fortrat**

UMR CNRS 6015 Inserm 1083, Centre Hospitalier Universitaire Angers, 4 Rue Larrey CEDEX 9, 49933 Angers, France; jofortrat@chu-angers.fr

Received: 25 March 2020; Accepted: 1 April 2020; Published: 6 April 2020

**Abstract:** Cardiovascular self-organized criticality (SOC) has recently been demonstrated by studying vasovagal sequences. These sequences combine bradycardia and a decrease in blood pressure. Observing enough of these sparse events is a barrier that prevents a better understanding of cardiovascular SOC. Our primary aim was to verify whether SOC could be studied by solely observing bradycardias and by showing their distribution according to Zipf's law. We studied patients with vasovagal syncope. Twenty-four of them had a positive outcome to the head-up tilt table test, while matched patients had a negative outcome. Bradycardias were distributed according to Zipf's law in all of the patients. The slope of the distribution of vasovagal sequences and bradycardia are slightly but significantly correlated, but only in cases of bradycardias shorter than five beats, highlighting the link between the two methods (r = 0.32; p < 0.05). These two slopes did not differ in patients with positive and negative outcomes, whereas the distribution slopes of bradycardias longer than five beats were different between these two groups (−0.187 ± 0.004 and −0.213 ± 0.006, respectively; p < 0.01). Bradycardias are distributed according to Zipf's law, providing clear insight into cardiovascular SOC. Bradycardia distribution could provide an interesting diagnosis tool for some cardiovascular diseases.

**Keywords:** baroreflex; heart rate variability; self-organized criticality; vasovagal syncope; Zipf's law

#### **1. Introduction**

Complexity, the final frontier of the cardiovascular system, has emerged as a major topic over the last decade. Complexity was initially discovered from incidental findings when studying cardiovascular variability. The meaning and implications of these findings have long remained unclear. Cardiovascular complexity has since been more precisely described over time. It became more and more difficult to integrate it into the current view of cardiovascular physiology that is largely dominated by the deterministic homeostatic principle. According to this principle, physiological variables are regulated and maintained at their normal values thanks to negative feedback regulatory loops. Today, complexity challenges the homeostatic view on the cardiovascular system [1–3]. We recently demonstrated that at least part of this complexity is explained by the self-organization of a cardiovascular system poised at criticality [4,5]. We showed that occurrences of spontaneous vasovagal events are distributed according to Gutenberg Richter's law. This law has been initially described in earthquakes occurrences: the magnitude plotted against the total number of earthquakes of at least this magnitude draws a straight line on a log-log graph. This finding explained how vasovagal reaction may occur. Vasovagal reaction is a parallel bradycardia and decrease in blood pressure of varying intensity from self-limiting symptoms to loss of consciousness and prolonged postictal asthenia [6]. During vasovagal syncope, the blood pressure decrease is not compensated by an increase of the heart rate as expected due to blood pressure homeostatic regulatory mechanisms. Brain perfusion is compromised because of the blood pressure decrease, and loss of consciousness eventually occurs. The self-organized pathophysiology of vasovagal syncope is a major finding, but its implication for cardiovascular physiology in general remains limited. Self-organized criticality has, however, emerged as a major topic in the study of dynamical systems and as a unifying theory across science fields, including physics, chemistry, ecology, and biology [1,2,7,8]. A better understanding of its meaning and implications for the cardiovascular system is needed. The study of cardiovascular self-organized criticality through vasovagal events requires continuous beat-by-beat recordings of blood pressure and heart rate. These recordings are difficult to obtain in some environmental conditions and are limited by time. These difficulties are a barrier toward a better understanding of cardiovascular self-organized criticality. Zipf's law has initially been described based on word occurrence in a text: the frequency of any word in a text is inversely proportional to its rank of occurrence [7,9]. This law has been inscribed into beat-by-beat recordings of the heart rate (Heart Rate Variability, HRV). These recordings show a linear distribution of occurrence of non-specific consecutive heart rate sequences across several beats, these sequences being the "words" of the cardiovascular system "language" [10–12]. Zipf's law represents another argument for cardiovascular self-organized criticality but without physiological and medical implications, contrary to Gutenberg Richter's law [2,4,10]. Beat-by-beat recordings of the heart rate are easy to obtain by means of commercially available heart rate monitors, facilitating the study of its complex dynamics [13,14]. However, it is still unknown whether Gutenberg Richter's law, determined by blood pressure and heart rate recordings, and Zipf's law, determined only by heart rate recordings, provide the same information. The goals of our study were to check, first, whether Zipf's law is observed specifically in bradycardia sequences, and second whether the meaning of Zipf's and Gutenberg Richter's laws overlap.

#### **2. Materials and Methods**

#### *2.1. Patients*

This study focused on patients with a history of iterative vasovagal syncope. Patients and flow charts have previously been extensively described [4]. One hundred consecutive patients who came to our department for advice on their iterative loss of consciousness and who gave their informed consent were included (51 female, 43 ± 2 years, 1.67 ± 0.01 m, 68 ± 1 kg, mean ±Standard Error of the Mean, SEM). Thirty patients were excluded because their interview was not suggestive of vasovagal syncope or because of a history of heart disease. A detailed medical history is central to the diagnosis of vasovagal syncope, but the head-up tilt test may help in both diagnosis and management. The head-up tilt test identified three patients with an orthostatic hypotension and five patients with a postural tachycardia syndrome. These eight patients were excluded. From the remaining 62 patients, 34 had a positive outcome to the head-up tilt test with (near) syncope symptoms, and 24 of them could be paired in age and sex with patients with a negative outcome. The group of patients with a positive outcome was called T+ patients (16 female, 39 ± 3 years, 1.66 ± 0.01 m, 67 ± 2 kg). The group of patients with a negative outcome vas called T- patients (16 female, 39 ± 3 years, 1.69 ± 0.02 m, 69 ± 3 kg). Patients received a complete description of the experimental procedure before giving their written informed consent. The Comité Consultatif de Protection des Personnes dans la Recherche Biomédicales des Pays de la Loire (Regional Committee for the Protection of Persons, #00/08, May 30th, 2000), France, approved the experiment, which is in accordance with the declaration of Helsinki, Finland.

#### *2.2. Head-Up Tilt Table Test*

The head-up tilt table test was performed in a quiet room with a comfortable ambient temperature (22–24 ◦C). The patient was lying on the table for at least ten minutes of adaptation to the supine position. The test began after 10 min in the supine position, and was followed by a 45 min period in the head-up position at an inclination of 70◦ by means of a motorized inclination table (AkronA8622, Electro-Medical Equipment, Marietta, GA, USA). The head-up position was stopped before 45 min elapsed in the event of (pre) syncopal symptoms defining the positive outcome. Cardiovascular

monitoring was performed during the whole test by means of an electrocardiogram and a digital blood pressure monitor for medical purposes (MACvu, Marquette, Milwaukee, WI, USA; and Finometer, FMS system, Amsterdam, Netherlands).

#### *2.3. Signal Analysis*

We followed recommendations to obtain accurate measurements of RR-intervals to analyze the smallest heart rate fluctuations [15]. Lead 2 of the electrocardiogram was digitized with a sampling frequency set at 500 Hz (AT-MIO-16, 12bits, Labview5.1, National Instruments, Austin, TX, USA). Intervals between R-peaks of the electrocardiogram were determined off-line by means of a peak detection algorithm. Electrocardiograms and time series of RR-intervals were visually inspected to identify R-peak misdetections and ectopic beats, which were manually deleted. Bradycardia sequences were identified on the time series of each patient, taking care not to include the large bradycardia of the syncopal episode and the preceding 30 s in cases of positive outcomes. A bradycardia sequence was defined as successive RR intervals with an increasing value. The length of a bradycardia sequence was defined as the total number of beats involved in the sequence. For each time series, bradycardia sequences were classified according to their length and were counted. The rank of bradycardia sequences of a same length was determined by classifying them according to their frequency of occurrence. For each patient, a diagram was plotted with the natural log of the rank on the x-axis and the natural log of the length of the matching bradycardia sequences on the y-axis. A linear regression was performed for each diagram in order to obtain the correlation coefficient and the slope. A previous study showed a cardiovascular Zipf's distribution according to two straight lines with a tipping point [16]. In this study, the position of a tipping point was determined by the best linear fits for each diagram.

#### *2.4. Vasovagal Events*

The method to assess and quantify the vasovagal events has previously been extensively described [4]. Vasovagal events were defined as consecutive beats with a drop in the mean blood pressure and an increase in the RR-interval. We classified and counted these events according to their length in number of beats.

#### *2.5. Statistics*

Data are presented as the mean ± SEM. Statistics were performed by means of Prism 5.01 (GraphPad Software, San Diego, CA, USA). We considered that the distribution fitted a straight line when |r| > 0.95. Tests for normality were performed by means of d'Agostino-Pearson omnibus K2 tests. Spearman correlations between Zipf's and Gutenberg Richter's law parameters were performed thanks to the data of a previously published study on the same data set focusing on this former law [4]. Matched patients with and without (pre)syncopal symptoms during the head-up tilt test (T+ and T−) were compared by means of a paired t test. We set the statistical significance at p < 0.05.

#### **3. Results**

T+ and T− patients had comparable anthropomorphic characteristics, medical history, treatments, heart rate, and blood pressure, as previously reported (66 ± 1 vs. 67 ± 2 bpm, 131 ± 4 vs. 127 ± 4 mmHg, and 77 ± 4 vs. 74 ± 3 mmHg, heart rate, systolic, and diastolic blood pressure, respectively) [4]. Electrocardiography recordings were of a high quality, and the visual review identified only several false R peak detections. Ectopic beats were observed in only seven patients (three T+ and three T-) and were sparse (maximum of two per 5 min on two T- patients). The quality of the tachograms was good (Figure 1).

**Figure 1.** Tachogram of a patient obtained during a head-up tilt table test. The tachogram is the beat-by-beat heart rate plotted against time (y axis and x axis, respectively). The first part of the tachogram is obtained in the supine position and ends at the vertical dashed line (10 min). The second part of the tachogram is the head-up position. The head-up position was stopped at the (pre)syncope occurrence (arrow). The analyzed time series started at the beginning of the tachogram and ended before the (pre)syncope occurrence, so it was excluded. The heart rate is shown here in beats per minute for convenience but is measured as RR-intervals (in ms) on the electrocardiographic signal. (bpm: beats per minute).

Bradycardia sequences were very frequent with no difference between T+ and T− and involved a large number of beats (36.2 ± 1.3 and 37.2 ± 1.1 beats per minute, respectively; Figure 2). Their maximal length was 12.1 ± 0.6 and 10.6 ± 0.6 beats for T+ and T−, respectively, with no difference between groups.

**Figure 2.** One minute of a patient's heart rate over time (y and x axes, respectively). Each heart beat is indicated by a square. Each box indicates a bradycardia sequence. The heart rate is shown in beats per minute for convenience but is measured as RR-intervals (in ms) on the electrocardiographic signal. (bpm: beats per minute).

Bradycardia sequences were distributed according to their rank along a straight line in all of the patients (T+ and T−). More precisely, they were distributed along two straight lines: one for the bradycardia sequences of a maximum of five beats and a second one for the longer bradycardia sequences with a coefficient correlation superior to 0.95 in all patients (T+ and T−; Figure 3). The position of the tipping point was the same in all patients.

**Figure 3.** Distribution of the number of bradycardia sequences according to their rank in one patient (log-log plot in natural logarithm). The pattern is the same for all patients including the position of the tipping point.

The slope of the relationship was significantly different between T+ and T− in the case of the long bradycardia sequences (Figure 4) but not in the case of shorter ones (−0.82 ± 0.05 and −0.78 ± 0.07 for T+ and T−, respectively; no unit; p = 0.686).

**Figure 4.** Absolute value of the slope of Zipf's distribution of bradycardia sequences longer than five beats in a group of patients with a negative outcome to the head-up tilt test and a matched group of patients with a positive outcome (T− and T+ respectively; \*\*: p < 0.01). n.u.: no unit.

The link between Gutenberg Richter's and Zipf's distributions was determined by comparing the slope of the linear relationship drawn by these two distributions. Gutenberg Richter's distribution was assessed through the slope of the distribution of vasovagal sequences. Zipf's distribution was assessed through the slope of the short and long bradycardia sequences. The slope values of the linear relationships were not normally distributed in cases of vasovagal and short bradycardia sequences (p < 0.0001 and p < 0.01, respectively), but were normally distributed in cases of long bradycardia sequences (p = 0.06). The slopes of the linear relationship of vasovagal and short bradycardia sequences were slightly but significantly correlated (r = 0.324, p = 0.02, Figure 5). There was no correlation between the slopes of the vasovagal and long bradycardia sequences' relationship and between the slopes of short and long bradycardia sequences (r = 0.117, p = 0.441, and r = 0.111, p = 0.465, respectively, Figure 5, vasovagal sequence data are from [4]).

**Figure 5.** Link between Gutenberg Richter's and Zipf's distribution according to the slope of the linear relationship drawn by these two distributions. Gutenberg Richter's distribution was assessed through the slope of the distribution of vasovagal sequences (axis label: Vasovagal). Zipf's distribution was assessed through the slope of the short and long bradycardia sequences (axis label: Short bradycardia and Long bradycardia, respectively). The coefficients of correlation and significance are mentioned in the text.

#### **4. Discussion**

This study confirms Zipf's law of cardiovascular dynamics. It also shows that Zipf's and Gutenberg Richter's distributions provide complementary information despite a link between these two distributions.

Several authors have dealt with Zipf's law of cardiovascular dynamics by means of different approaches. To our knowledge, Kalda et al. were the first to demonstrate Zipf's law of cardiovascular dynamics [10]. These authors studied the statistical similarities of short series of RR-intervals. Yang et al. performed an analysis based on logic variations of consecutive heart beats, while Rodriguez et al. looked for the presence of mathematical properties of Zipf's series in RR-interval recordings [11,12]. The question remains whether these approaches used the most efficient way to assess Zipf's law because the natural language of the cardiovascular system remains unknown. In our study, we began with a physiological observation to attempt to align as far as possible with this unknown language. This physiological observation is a conundrum. It is possible to observe bradycardia episodes on consecutive heart beats when the heart rate spontaneously fluctuates in a normal subject [17]. These episodes contradict the deterministic homeostatic regulation of the cardiovascular function. A fast regulatory mechanism called the baroreflex should detect blood pressure fluctuations and compensate beat-by-beat for these fluctuations by adjusting the heart rate. A decrease in the heart rate decreases blood pressure, but the baroreflex should increase the heart rate to stop the blood pressure decrease. A bradycardia episode of several beats is not expected in cases of a well-working deterministic homeostatic baroreflex. A prolonged bradycardia episode could eventually lead to large arterial hypotension with compromised brain perfusion and loss of consciousness. The phenomenon is called vasovagal syncope and could paradoxically occur in any normal subject despite a well-working baroreflex [6]. In this study, we tried to stay close to the natural physiological language of the cardiovascular system. This approach allowed us to define a simple method with strong evidence of Zipf's law in the cardiovascular dynamics. However, further studies may help to better define the natural cardiovascular language to better characterize Zipf's law and the self-organized properties of the cardiovascular function.

Our approach also differs from the previous studies on Zipf's law in the cardiovascular system. Heart rate variability recordings were obtained by means of a Holter monitor in all of the three previous studies, while we studied quiet unmoving patients [10–12]. A Holter monitor is a device that records the heart rate during the normal daily life of the patient. The heart is constantly influenced by the various demands of daily life that include activities, stressors, and body position. On Holter recordings, heart rate variability is the result of daily life but is also intrinsic to regulatory mechanisms and their physiological delays. Daily life variability is totally absent in immobile and quiet patients, and solely the intrinsic variability remains. We previously demonstrated the influence of differences in experimental set-up in analyses of heart rate variability, including those with a focus on its complex dynamics [18,19].

The definitions of the sequences to study Gutenberg Richter's and Zipf's laws both included bradycardia episodes on consecutive heart beats. This point therefore identifies overlap between the meaning of Gutenberg Richter's and Zipf's distributions with these two close definitions. Thus, the slopes of Gutenberg Richter's and Zipf's distributions are correlated. However, only the short bradycardia sequences are correlated with vasovagal sequences and not the long ones. Moreover, the distribution of long bradycardia sequences differs between T− and T+ patients, contrary to Gutenberg Richter's distribution (Figure 4). This difference shows that Gutenberg Richter's and Zipf's laws provide complementary information about cardiovascular self-organized criticality.

Telling the difference between patients with and without a positive result in the diagnosis tool for vasovagal syncope is a challenge of cardiovascular medicine. Patients with vasovagal syncope are usually apparently healthy after a regular medical check-up [6]. The baroreflex is functioning well in these patients, who generally maintain their cardiovascular function well. Vasovagal syncope has remained a medical mystery for centuries [6]. Only recently, some studies focusing mainly on the complex dynamics of heart rate variability convincingly showed a difference between patients with a positive outcome of the diagnosis tool and patients with a negative one. Graff et al. demonstrated this difference by means of the entropy method, while Fortrat et al. achieved this by defining a marker of cardiovascular instability [20,21]. Questions remain about whether these two studies and Zipf's law finally focused on the same complex properties by means of different tools or whether the different methods provide complementary and unrelated information.

The main limitation of this study is the analysis of the supine and standing parts of the head-up tilt table test into single time series. The standing position requires major cardiovascular adaptions partly driven by the autonomic nervous system [22]. Cardiovascular dynamics and heart rate variability are different between these two positions, as demonstrated by means of signal analysis or only by looking at the time series (Figure 1). We previously demonstrated the influence of body position on the cardiovascular Zipf's distribution [16]. The analysis of a single time series for a whole head-up tilt table test was however necessary to collect enough of the sparse vasovagal events in order to perform Gutenberg Richter analysis [4]. Our study demonstrated that the focus can be placed on the heart rate to study cardiovascular self-organized criticality. Specifically, future studies should clarify the influence of body position adaptation on cardiovascular self-organized criticality and investigate the influence of autonomic nervous system adaptations on SOC. Future studies should also confirm that the tipping point of the bradycardia Zipf's distribution is linked to physiological influences like breathing and not to analysis artifacts like a finite size effect.

#### **5. Conclusions**

This study confirmed Zipf's law of cardiovascular dynamics and demonstrated that Zipf's and Gutenberg Richter's laws explore complementary aspects of the self-organized criticality of cardiovascular dynamics. Zipf's law provides an interesting and easy to implement tool to better characterize the self-organized criticality of cardiovascular dynamics. Zipf's law may also provide an interesting tool for the medical diagnosis of some cardiovascular diseases.

**Funding:** This research received no external funding.

**Acknowledgments:** We thank the patients for their cooperation. We thank the Centre Hospitalier Universitaire d'Angers (Angers University Hospital) and its Centre de Recherche Clinique (Clinical Research Centre) for their help with the administrative procedures for biomedical research projects.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**

1. Struzik, Z.R. Is heart rate variability dynamics poised at criticality? In Proceedings of the 8th Conference of the European Study Group on Cardiovascular Oscillations (ESGCO), Trento, Italy, 25–28 May 2014; IEEE: Piscataway, NJ, USA, 2014.


© 2020 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Day and Night Changes of Cardiovascular Complexity: A Multi-Fractal Multi-Scale Analysis**

#### **Paolo Castiglioni 1,\*, Stefano Omboni 2,3, Gianfranco Parati 4,5 and Andrea Faini <sup>5</sup>**


Received: 14 March 2020; Accepted: 16 April 2020; Published: 18 April 2020

**Abstract:** Recently, a multifractal-multiscale approach to detrended fluctuation analysis (DFA) was proposed to evaluate the cardiovascular fractal dynamics providing a surface of self-similarity coefficients α(*q*,τ), function of the scale τ, and moment order *q*. We hypothesize that this versatile DFA approach may reflect the cardiocirculatory adaptations in complexity and nonlinearity occurring during the day/night cycle. Our aim is, therefore, to quantify how α(*q*, τ) surfaces of cardiovascular series differ between daytime and night-time. We estimated α(*q*,τ) with −5 ≤ *q* ≤ 5 and 8 ≤ τ ≤ 2048 s for heart rate and blood pressure beat-to-beat series over periods of few hours during daytime wake and night-time sleep in 14 healthy participants. From the α(*q*,τ) surfaces, we estimated short-term (<16 s) and long-term (from 16 to 512 s) multifractal coefficients. Generating phase-shuffled surrogate series, we evaluated short-term and long-term indices of nonlinearity for each *q*. We found a long-term night/day modulation of α(*q*,τ) between 128 and 256 s affecting heart rate and blood pressure similarly, and multifractal short-term modulations at *q* < 0 for the heart rate and at *q* > 0 for the blood pressure. Consistent nonlinearity appeared at the shorter scales at night excluding *q* = 2. Long-term circadian modulations of the heart rate DFA were previously associated with the cardiac vulnerability period and our results may improve the risk stratification indicating the more relevant α(*q*,τ) area reflecting this rhythm. Furthermore, nonlinear components in the nocturnal α(*q*,τ) at *q* - 2 suggest that DFA may effectively integrate the linear spectral information with complexity-domain information, possibly improving the monitoring of cardiac interventions and protocols of rehabilitation medicine.

**Keywords:** multifractality; multiscale complexity; detrended fluctuation analysis; heart rate; blood pressure; self-similarity

#### **1. Introduction**

Time-series complexity is common in physiology. In fact, physiological systems often exhibit fractal geometries and are composed of several elements interacting nonlinearly, which are both typical features of a complex system [1]. The cardiovascular system, in particular, can be described as a complex, dynamical system because it is composed of a fractal network of branching tubes, the vasculature, connecting individual vascular beds that interact with each other to harmonize globally the local needs of blood supply. The overall cardiovascular regulation modulates the local blood flows thanks to the integrative control of the autonomic nervous system operating through effectors and feedbacks (the baro- and chemoreflexes) with nonlinear elements.

Complex dynamical systems are not characterized by an intrinsic time scale. This means that their derived time series may appear statistically self-similar when plotted at different scales. For this reason, the interest in methods that quantify self-similar (or fractal) properties of the cardiovascular dynamics is increasing. A very popular method is based on the detrended fluctuation analysis (DFA), which provides a self-similarity scale coefficient, α, directly related to the Hurst's exponent [2]. When DFA was originally proposed for the analysis of heart rate variability, it described a bi-scale fractal model providing a short-term coefficient (α1) for scales shorter than 16 beats and a long-term coefficient (α2) for longer scales [3]. The original bi-scale method was then extended in two ways. One way was to provide a multiscale spectrum of self-similarity coefficients, a function of the scale *n* in beats, α(*n*) [4–6]. Another way was to provide a multifractal spectrum of self-similarity coefficients, a function of the moment order *q*, α(*q*) [7,8]. The multifractal spectrum includes *q* = 2—the second-order moment used in the original DFA method for monofractal series—and allows detecting multifractality when α(*q*) differs substantially between positive and negative *q* orders. The multiscale and the multifractal methods were finally combined in the multifractal-multiscale DFA, a versatile approach that describes multifractal structures localized over specific scales and that provides a surface of scale coefficients, α(*q,n*) [9]. Recent works demonstrated the capability of the multifractal-multiscale DFA of heart rate variability to classify different types of cardiac patients [10] and to describe alterations in the heart rate complexity due to an impaired integrative autonomic control in paraplegic individuals [11].

It is less clear, however, whether complexity methods based on DFA can quantify nonlinear components. In this regard, theoretical analyses affirm that the information on the Hurst's exponent provided by the second-order moment DFA can be derived mathematically from the power spectrum, which is a linear method of analysis [12,13]. Actually, empirical quantifications of the degree of nonlinear information of the cardiovascular dynamics provided by the more advanced multifractal-multiscale approaches are missing.

In this work, we hypothesize that the versatile multifractal-multiscale DFA approach may reflect the cardio-circulatory adaptations in the overall complexity and, particularly, in the nonlinear dynamics of the cardiovascular time series that may occur during the day/night cycle. Circadian rhythms and differences in activity levels between daytime and night-time hours are expected to have a major influence on cardiovascular regulation. Knowing how this happens may help to better identify and interpret possible alterations associated with pathological conditions. In this regard, a description of the α(*q,n*) circadian modulations may be important in the rehabilitation medicine for correctly monitoring changes associated with treatments or the recovery from clinical interventions. To our knowledge, no studies addressed the quantification of the changes in the multifractal-multiscale DFA of heart rate variability associated with the day–night cycle. Furthermore, most of the studies on cardiovascular complexity are based on the analysis of heart rate variability only. This is due to the difficulty to better describe the status of the system by measuring other cardiovascular variables beat-by-beat in addition to the heart rate, as the systolic and the diastolic arterial blood pressure.

Therefore, our work aims to address the above open issues on cardiovascular complexity by quantifying the fractal dynamics of heart rate and blood pressure, the degree of nonlinearity, and possible night–day modulations of complexity. This will be done analyzing continuous 24-h blood pressure recordings and comparing self-similarity coefficients estimated by the multifractal-multiscale approach over daytime and night-time. In particular, we will define new indices of the degree of nonlinearity based on the multifractal-multiscale DFA to quantify the additional information provided by this complexity method compared to traditional spectral methods.

#### **2. Materials and Methods**

#### *2.1. Subjects and Data Collection*

The study is based on a historic database of 24-h ambulatory intra-arterial blood pressure recordings obtained at the University Hospital of Milan (Ospedale Maggiore Policlinico, Milan, Italy), for the diagnosis of hypertension [14]. Recordings were performed between the 1980s and the 1990s when intermittent noninvasive arm devices were not still in use in the clinical practice.

As inclusion criteria, we selected only adult (>18 yro) normotensive subjects in which the suspected hypertensive state was excluded after the clinical evaluation. Exclusion criteria were smoking; obesity; clinical or laboratory evidence of health abnormalities, like cardiovascular disease or diabetes; prior drug treatment for hypertension; any alteration in glucose metabolism or renal function; and administration of cardiovascular drugs in the 4 weeks preceding the recording. We also excluded blood pressure tracings of inadequate quality for a 24-h analysis. This led to selecting recordings of N = 14 normotensive subjects (3 females of which one in the childbearing age and 11 males) with age between 19 and 64 years.

Details of data collection are reported in [14]. Briefly, a catheter inserted into the radial artery of the non-dominant arm was connected to a transducing-perfusing unit secured to the thorax at the heart level. The blood pressure signal was stored on a magnetic tape recorder bound to the waist. During the recordings, the subjects were free to move within the hospital. Mealtimes and bedtimes were standardized. The blood pressure signal was digitized (170 Hz, 12 bits) and edited manually from movement artifacts, pulse pressure dampening, and premature beats. Each pulse wave was identified by a derivative-and-threshold algorithm [15]; systolic blood pressure (SBP) and diastolic blood pressure (DBP) were calculated for each pulse wave beat-by-beat. As suggested in [16], a parabolic interpolation refined the SBP fiducial point before calculating the inter-beat interval (IBI) as the interval between the times of occurrence of consecutive systolic peaks.

Two sub-periods were selected for the analysis after visual inspection of the tracings: the "Day" subperiod during daytime in the afternoon, when the subjects were not lying in bed and were free to perform normal daytime activities; the "Night" subperiod after 11 PM when the participants were asleep according to the schedule of the hospital. The selected segments had to be composed of at least 14,000 heartbeats, with a duration of at least 4 h during daytime and of at least 5 h during night-time, without evident nonstationarities.

The study was carried out after having obtained informed consent from the participants in accordance with the 1975 Declaration of Helsinki and following the recommendations of the ethical committee of the Ospedale Maggiore Policlinico (Milan, Italy).

#### *2.2. Multifractal-Multiscale Detrended Fluctuation Analysis*

We estimated the multifractal multiscale structure of the IBI, SBP, and DBP time series by the fast DFA algorithm available in [17]. Given the beat-by-beat series *xi* of length *L* beats, we calculated its cumulative sum, *yi*. We split *yi* into *M* maximally overlapped blocks of *n* beats (two consecutive blocks have *n*-1 beats in common). We detrended each block with least-square polynomial regression and calculated the variance of the residuals in each *k-th* block, σ<sup>2</sup> *<sup>n</sup>*(*k*). The variability function *Fq*(*n*) is the *q*-th moment of σ<sup>2</sup> *n* [7]:

$$\begin{cases} F\_q(n) &= \left(\frac{1}{M} \sum\_{k=1}^{M} \left(\sigma\_n^2(k)\right)^{q/2}\right)^{1/q} \text{ for } q \neq 0\\ &\quad F\_q(n) = \left. \varepsilon^{\frac{1}{2M} \sum\_{k=1}^{\infty} \ln\left(\sigma\_n^2(k)\right)} \right. & \text{ for } q = 0 \end{cases} \tag{1}$$

We evaluated Equation (1) for *q* between −5 and +5 and block sizes *n* between 6 and *L*/4 beats. We evaluated the multifractal multiscale coefficients as a function of the beat-scale *n*, αB(*q*,*n*), calculating the derivative of log *Fq*(*n*) vs. log *n* [17]. This was done for detrending polynomials of order 1 and 2 (see examples of the corresponding *Fq*(*n*) estimates in Figure 1. Previous empirical analyses suggested that the second-order polynomial overfits block sizes shorter than 12 beats, but at the same time, it appears to more efficiently remove long-term trends [17–19]. Therefore, we estimated a single αB(*q*,*n*) function combining the estimates after detrending of order 1 and 2 with a weighted average which weights more the order one at the shorter scales as proposed in [17].

**Figure 1.** Multifractal variability functions *Fq*(*n*) for inter-beat interval (IBI) with different orders of detrending polynomials: average over the group of participants. The *Fq*(*n*) functions are plotted in blue for *q* > 0, in black for *q* = 0, and in red for *q* < 0; the dashed line is *q* = 2, second-order moment of the traditional monofractal detrended fluctuation analysis (DFA). Upper panels: *Fq*(*n*) estimated with 1st order (linear) detrending during (**a**) Day and (**b**) Night. Lower panels: *Fq*(*n*) with 2nd-order (quadratic) detrending during (**c**) Day and (**d**) Night.

It should be noted that the parameter *q* in Equation (1) defines the moment order calculated for the variances of the residuals. In the traditional monofractal DFA, the variability function is defined as the root-mean-square of σ<sup>2</sup> *<sup>n</sup>*, which corresponds to the second-order moment, or *q* = 2. If the series is monofractal, all the moment orders *q* provide the same slope α. By contrast, for multifractal series positive moment order *q* weight more the contribution of the fractal components with greater amplitude, negative moment order *q* weight more the contribution of the fractal components with lower amplitude.

To compare Day and Night periods over the same temporal scales, in seconds, we mapped the scale units from number of beats, *n*, to time τ, in seconds, with the transformation

$$
\pi = n \times \mu\_{\text{IBI}} \tag{2}
$$

with

$$
\mu\_{\text{IBI}} = \frac{1}{L} \sum \mathbf{i}\_i^L = {}\_1 {}^1 {}\_2 \mathbf{i}\_i \tag{3}
$$

the mean of the IBI values for all the *L* beats composing the whole time series, in seconds. The obtained coefficients expressed as a function of the time scale, α*(q,*τ*)*, were spline-interpolated over τ and resampled (2500 points evenly spaced over the logarithmic τ axis from 6 s to 3600 s) to have estimates at the same temporal scales for each recording. Realigning the DFA coefficients in this way allows properly comparing the same time scales between conditions in which the cardiovascular signals are sampled at different heart rates. We considered scales between 8 and 2048 s. The largest scale (τ = 2048 s) is estimated on more than seven independent blocks of data even in the case of the recording with the shortest duration: this assures sufficient stability of the estimate as shown empirically in previous validations [5,20]. Scales shorter than τ = 8 s were not considered because at negative *q* orders high levels of estimation bias may be present [17].

We introduced multifractal short-term and long-term coefficients to concisely describe the multifractal multiscale structure. This was done by averaging α(*q*,τ) over short scales, with 8 ≤ τ ≤ 16 s, and over long scales, with 16 < τ ≤ 512 s, obtaining the multifractal short-term coefficient αS(*q*) and long-term coefficient αL(*q*).

#### *2.3. Nonlinearity Index*

For each series *j* we generated 100 Fourier phase-randomized series by shuffling the spectrum of the phases with the code available in [21]. This procedure removes possible nonlinear components in the dynamics of the original series, preserving its power spectrum and therefore the original firstand second-order moments [22]. Then, we calculated the multifractal multiscale coefficients of each of the 100 surrogates, α*i,j*(*q*,τ) with 1 ≤ *I* ≤ 100, to be compared with the coefficients of the original series *j*, αO,*<sup>j</sup>* (*q*,τ). For the comparison, we calculated π*<sup>j</sup>* (*q*,τ), defined at each *q* and τ as the percentile of the distribution of 100 surrogate α*i,j*(*q*,τ) coefficients in which was the original αO,*<sup>j</sup>* (*q*,τ) coefficient (to apply a 2-tail statistics, percentiles greater than 50% were transformed into their complement to 1 as in [23]). π*<sup>j</sup>* (*q*,τ) may range between 50% and 0%: the lower its value, the more significant the deviation of the original scale coefficient αO,*<sup>j</sup>* from the distribution of the 100 surrogate coefficients α*i*,*<sup>j</sup>* . Large deviations from the surrogates distribution are suggestive of nonlinear components in the original series. Therefore, we defined a short-term nonlinearity index at each moment order *q*, *NL*S(*q*), by calculating the percentage of scales in the range 8 ≤ τ ≤ 16 s, where π*<sup>j</sup>* (*q*,τ) was ≤1%. Similarly, we calculated the percentage of scales with π*<sup>j</sup> (q,*τ*)* ≤ 1% for 16 < τ ≤ 512 s to define the long-term nonlinearity index *NL*L(*q*). Both *NL*S(*q*) and *NL*L(*q*) may range between 0% and 100%. Their higher values indicate moment orders *q* that better detect the presence of nonlinear components.

#### *2.4. Spectral Analysis*

The IBI, SBP, and DBP beat-by-beat series were interpolated evenly at 5 Hz before spectral analysis. Power spectra were estimated by the Welch periodogram with 80% overlapped Hann data windows of 240 s length. The spectra were integrated over the very-low frequency (VLF, between 0.003 and 0.04 Hz), the low frequency (LF, between 0.04 and 0.15 Hz), and the high-frequency (HF, between 0.15 and 0.4 Hz) bands as indicated in the guidelines [16].

#### *2.5. Statistical Analysis*

The α(*q*,τ) coefficients of the *N* = 14 participants were compared between *Day* and *Night* at each τ and *q* by the Wilcoxon signed-rank test. The multifractal short- and long-term coefficients, αS(*q*) and αL(*q*), and nonlinearity indices, *NL*S(*q*) and *NL*L(*q*), were also compared between *Day* and *Night* at each *q* by the Wilcoxon signed-rank test. IBI, SBP, and DBP levels and power spectra were compared between *Day* and *Night* by the paired t-test, after log-transformation of the spectral indices to remove the skewness of their distribution [24]. The threshold for statistical significance was set at 5% with a two-sided alternative hypothesis. All the tests were performed with "R: A Language and Environment for Statistical Computing" software package (R Core Team, R Foundation for Statistical Computing, Vienna, Austria, 2019).

#### **3. Results**

#### *3.1. Day vs. Night*

The data segments selected for the analysis of the *Day* and *Night* periods were composed by a similar number of heartbeats: 20,779 (2744) beats during the *Day* and 20,329 (4059) beats during the *Night*, as average (SD) over the group. Means and spectral powers of the cardiovascular series are reported in Table 1. Because of the higher heart rate during the daytime, the segment duration was shorter in the *Day*, i.e., 4 h 30 (30 ), than in the *Night* period, i.e., 5 h 42 (36 ). For the same reason, the scale τ = 16 s that divides the αS(*q*) and αL(*q*) indices corresponds on average to 20.7 beats in the *Day* and 15.5 beats in the *Night* period, and the αL(*q*) upper scale at τ = 512 s corresponds to 661.2 beats and 495.3 beats in the *Day* and *Night* periods, respectively.


**Table 1.** Mean levels and spectral powers of cardiovascular series.

Values as mean (SD); *p* value after *T* test on log-transformed powers.

Figure 2 shows the α(*q*,τ) surfaces for IBI, SBP, and DBP, separately, during *Day* and *Night* periods (average over the group of patients). The figure suggests the presence of structural differences between heart rate and blood pressure in their complex dynamics: during the daytime, these differences appear particularly clear between 16 and 256 s, where IBI appears characterized by a relatively flat surface at all *q* orders while SBP and DBP show a dip around τ = 32 s for positive *q* orders.

**Figure 2.** Surfaces of multifractal multiscale DFA coefficients, α(*q*,τ), during Day and Night periods. Average over 14 participants, for scales τ between 8 and 2048 s and moment orders *q* between −5 and +5; IBI = inter-beat-interval; SBP = systolic blood pressure; DBP = diastolic blood pressure.

Even more obvious is the difference between daytime and night-time in each cardiovascular series. The difference is particularly evident for the IBI surface of scale coefficients, which shows a marked decrease of the α coefficients at scales between 128 and 256 s during the night, more pronounced at negative *q* orders. Similar deflections appear to also characterize the surfaces of DFA scale coefficients of SBP and DBP.

Figure 3 compares *Day* and *Night* cross sections of the α(*q*,τ) surfaces at each moment order *q*. As to IBI, differences at scales shorter than 16 s regard two distinct *q*-τ areas. In the main area, centered at *q* = −2, α is lower at night, while in the secondary narrower area, centered at *q* = 4, α is greater at night.

Remarkable *Day–Night* differences with lower α at night also regard scales between 128 and 256 s. They are evident at all the moment orders but extend over a larger range of scales τ for negative *q* values.

**Figure 3.** Day–Night comparison of cross sections of multifractal multiscale DFA coefficients. (**a**) Cross sections of α(*q*,τ) of IBI for scales τ between 8 and 2048 s and moment orders *q* between −5 and +5: average over the group of 14 participants in the *Day* subperiod; *q* < 0 in red, *q* > 0 in blue, *q* = 0 in black; the dotted line is α for *q* = 2 (second order moment of the monofractal DFA); (**b**) α(*q*,τ) of IBI as in panel (**a**) for the *Night* subperiod; (**c**) color map representing the statistical significance (*p* value) of the *Day* vs. *Night* comparison of IBI scale coefficients calculated at each τ and *q* after the Wilcoxon signed rank test; (**d**) α(*q*, τ) of SBP in the *Day* subperiod represented as in panel (**a**); (**e**) α(*q*, τ) of SBP in the *Night* subperiod represented as in panel (**a**); (**f**) color map of the *Day* vs. *Night* statistical significance for SBP scale coefficients; (**g**) α(*q*, τ) of DBP during *Day* represented as in panel (**a**); (**h**) α(*q*, τ) of DBP during *Night* represented as in panel (**a**); (**i**) color map of the *Day* vs. *Night* statistical significance for DBP coefficients.

Similarly to IBI, also the α(*q*,τ) coefficients of SBP and DBP show a significant decrease at *Night* for scales between 128 and 256 s for all the *q* orders. Significant *Day–Night* differences with greater α at night also appear in blood pressure at scales shorter than 32 s, but, differently from IBI, the changes are significant for positive *q* only.

The detailed representation of Figure 3 is summarized by the multifractal short- and long-term coefficients in Figure 4. The IBI multifractal short-term coefficient is significantly lower at night for −3 ≤ *q* ≤ 0, while the long-term coefficient is significantly lower at night for *q* ≤ 1. Moreover, the multifractal short-term coefficients of blood pressure are higher at night when *q* ≥ 2 (Figure 4c,e), while the long-term coefficients, as for IBI, are lower at night mainly for negative *q* (Figure 4d,f).

**Figure 4.** Day–Night comparison of multifractal short- and long-term coefficients. (**a**) Short-term coefficients, αS(*q*), for IBI in *Day* (open circles) and *Night* (solid circles) periods and for −5 ≤ *q* ≤ +5: median ±standard error of the median over N = 14 participants; the \* indicates *Day* vs. *Night* differences significant at *p* < 0.05; (**b**) long-term coefficients, αL(*q*), of IBI represented as in panel (**a**); (**c**) short-term coefficients of SBP and (**d**) long-term coefficients of SBP, represented as in panel (**a**); (**e**) short-term coefficients and (**f**) long-term coefficients of DBP, represented as in panel (**a**).

#### *3.2. Nonlinearity*

Figure 5 illustrates the degree of nonlinearity detected comparing α(*q*,τ) of the original and surrogate series during the daytime. A common feature to heart rate and blood pressure is the evidence

of nonlinear components at scales shorter than 64 s at all *q* but *q* = 2 (the moment order of the traditional monofractal DFA).

**Figure 5.** Assessment of nonlinearity during daytime. Upper panels refer to IBI: (**a**) α(*q*,τ) coefficients for the original series (average over N = 14 participants, see panel (**a**) for line colors); (**b**) α(*q*,τ) for the corresponding phase-randomized surrogate series; (**c**) color map of the percentile of the distribution of surrogate estimates in which is the original estimate (average over N = 14 participants). Mid panels refer to SBP: (**d**) α(*q*,τ) for the original series; (**e**) α(*q*,τ) for the corresponding surrogate series; (**f**) color map of percentiles. Lower panels refer to DBP: (**g**) α(*q*,τ) for the original series; (**h**) α(*q*,τ) for the corresponding surrogate series; (**i**) color map of percentiles.

At *Night* nonlinear components are more evident (Figure 6) and affect longer scales, particularly for IBI. Estimates at *q* = 2 appear to be linear also at night-time.

**Figure 6.** Assessment of nonlinearity during night-time. Upper panels refer to IBI: (**a**) α(*q*,τ) for the original series (average over N = 14 participants, see panel (**a**) for of line colors), (**b**) α(*q*,τ) for the corresponding phase-randomized surrogate series, and (**c**) color map of the percentile of the distribution of surrogate estimates in which is the original estimate (average over N = 14 participants). Mid panels refer to SBP: (**d**) α(*q*,τ) for the original series; (**e**) α(*q*,τ) for the surrogate series; (**f**) color map of percentiles. Lower panels refer to DBP: (**g**) α(*q*,τ) for the original series; (**h**) α(*q*,τ) for the surrogate series; (**i**) color map of percentiles.

Figure 7 summarizes these findings showing the short-term and the long-term nonlinearity indices, *NL*S(*q*) and *NL*L(*q*). The highest degree of nonlinearity is detected at *Night* by *NL*S(*q*), which is close to 100% for all the cardiovascular series between *q* = −2 and *q* = +4, with the notable exception of *q* = 2. In fact, at *q* = 2 *NL*<sup>S</sup> falls to 0% for all the signals. *NL*<sup>S</sup> tends to be higher at night with significant differences at some *q* < 0 for IBI and DBP and at *q* > 2 for IBI. Long-term nonlinear components are mainly present in IBI at night. In fact, *NL*L(*q*) of IBI is greater than 50% during night-time at all *q* but

*q* = 2. Furthermore, it is significantly greater at night for all *q* - 2. *NL*<sup>L</sup> too is close to 0% at *q* = 2, both during *Day* and *Night*, for heart rate and blood pressure.

**Figure 7.** Day vs. Night comparison of short-term and long-term indices of nonlinearity. (**a**) Short-term index, *NL*S(*q*), for IBI in *Day* (open circles) and *Night* (solid circles) periods and for −5 ≤ *q* ≤ +5: median ±standard error of the median over N = 14 participants; the \* indicates *Day* vs. *Night* differences significant at *p* < 0.05; (**b**) long-term nonlinearity index, *NL*L(*q*), of IBI; (**c**) short-term and (**d**) long-term nonlinearity index of SBP; (**e**) short-term and (**f**) long-term nonlinearity index of DBP.

#### **4. Discussion**

This work compared patterns of blood pressure and heart rate complexity between daytime and night-time as assessed by the multifractal multiscale DFA approach. To our knowledge, this is the first study addressing night–day changes of multifractality in different cardiovascular signals and on a continuum spectrum of temporal scales. Our work revealed specific scales τ and specific fractal components (as identified by *q*) where the cardiovascular complexity differs between wake at daytime and sleep at night. Furthermore, it introduced new indices of nonlinearity which highlight the areas of the α(*q*,τ) surface that better reflect the nonlinear dynamics. A brief discussion of these points follows.

#### *4.1. Day vs. Night*

Mean levels and spectral powers of heart and blood pressure (Table 1) reflect the day–night changes reported previously [25,26], i.e., lower heart rate and blood pressure at night due to the lower levels of physical activity and to the lying position, which are associated with a higher cardiac vagal tone (HF power of IBI), a lower cardiac sympatho/vagal balance (LF/HF powers ratio of IBI), and a lower vascular sympathetic tone (LF power of SBP and DBP [27,28]).

In addition to these known changes in heart rate and blood pressure mean levels and spectral powers, we reported clear changes in the α(*q*,τ) fractal structure. In IBI, the more evident change is the night decrease of coefficients around 128–256 s (Figure 4c). The decrease affects all the moment orders but it is amplified at negative *q* and thus the night/day modulation of the long-term multifractal index αL(*q*) is larger for *q* < 0 (Figure 4b). We may associate this night/day oscillation to an endogenous circadian rhythm previously described in the heart rate by a monofractal DFA exponent (i.e., for *q* = 2) estimated over scales between 20 and 400 beats [29]. This endogenous rhythm was hypothesized to contribute to the period of the cardiac vulnerability reported in epidemiological studies. Our work suggests that this night/day rhythm (1) is highlighted by a multifractal approach that assesses negative moment orders and (2) is better quantified in a narrower range of scales, between 128 s and 256 s. Therefore, our results may prove to be of clinical importance by allowing designing new tools for the complexity analysis of heart rate that better stratify the cardiovascular risk. Interestingly, our study also provides evidence that a night/day modulation with greater daytime values is present at the same scales in blood pressure too, suggesting that a common physiological mechanism is at the origin of the circadian oscillation in the heart rate and the blood pressure self-similarity coefficients.

By contrast, night–day changes at shorter scales affect heart rate and blood pressure differently. While short-term coefficients of blood pressure are greater at night for moment orders *q* ≥ 2, the main modulation of short-term scales of heart rate consists of lower values at night for −3 ≤ *q* ≤ 0 (Figure 4a). Further studies controlling the effects of posture and physical activity are needed to understand the nature of so different night/day changes between heart rate and blood pressure.

Night–day modulations of the heart rate self-similarity coefficients were also reported in a study on 24-h Holter's recordings performed on a large population of healthy subjects [30]. This study applied the bi-scale model as in [3], which originally defined a short-term coefficient α<sup>1</sup> for scales between 4 and 16 beats and a long-term coefficient α<sup>2</sup> for scales between 16 and 64 beats. The study in [30] used the scale *n* = 11 beats to separate α<sup>1</sup> from α<sup>2</sup> and reported a significant decrease in α<sup>1</sup> at night. We did not consider scales short as in this study because at τ < 8 s the multifractal estimates can be affected by large estimation bias for negative *q* orders. However, α<sup>1</sup> and the LF/HF powers ratio of the heart rate are correlated [13] and the reduction in the LF/HF powers ratio we reported at night in Table 1 is coherent with the night reduction of α<sup>1</sup> in [30]. These authors, however, also showed a significant increase of α<sup>2</sup> at night, which appears in contrast with the night decrease of the long-term scale coefficients reported both in [29] and in our work. To correctly interpret the results of the three studies, we should consider carefully the scale ranges where the coefficients are estimated. To illustrate this point, Figure 8 plots the coefficients we calculated as the derivative of log *Fq*(*n*) vs. log *n* in Equation (1), i.e., αB(*q*,*n*), for *q* = 2. The scale *n* is expressed in beats to facilitate the comparison with previous studies [29,30]. As the estimation bias is negligible for *q* = 2, α<sup>B</sup> is plotted from *n* = 6 beats. The night/day comparison shows a significant nocturnal decrease of α<sup>B</sup> at scales < 11 beats, in line with the α<sup>1</sup> results in [30], and greater night-time values at scales where α<sup>2</sup> was estimated in [30]. These greater values correspond to the small area of statistical significance that appears in our Figure 3 at scales τ ≤ 16 s and at orders *q* ≥ 2. The α coefficient calculated in [29] between 20 and 400 beats overlaps partially with α<sup>2</sup> but covers a much wider range of longer scales, which includes the band between 128 and 256 beats where we found a significant night decrease of αB. Therefore our study and the studies in [29,30] provide coherent results if the correct scale ranges are considered. The comparison of Figure 8 also highlights the importance to provide estimates of the scale coefficients as a

continuous function of the scale *n* to correctly identify phenomena which may occur in nearby scales with different characteristics.

**Figure 8.** Day–night comparison of the multiscale monofractal DFA coefficients of IBI plotted vs. the block size *n* in beats. (**a**) αB(*q*,*n*) calculated for *q* = 2 (second-order moment of the monofractal DFA) during daytime (red) and nighttime (blue): mean +/− sem over the group of N = 14 participants; the arrows indicate the scale ranges for estimating α<sup>1</sup> and α<sup>2</sup> as defined by Vanderput et al. in [31] and for estimating α as defined by Hu et al. in [30]. (**b**) W statistics for the day-night difference in αB(2,*n*); when W is above the red horizontal line, the difference at the corresponding scale is significant at *p* < 5%.

#### *4.2. Nonlinearity*

The comparison between original and Fourier-shuffled surrogates allowed us defining two concise indices of nonlinearity, *NL*S(*q*) and *NL*L(*q*), that indicate the moment orders and the scale ranges, where α(*q*,τ) provides information on nonlinear dynamics. These indices are close to 0% for *q* = 2, supporting previous theoretical speculations indicating that the monofractal DFA and the power spectrum provide similar information [12,13]. However, we also found clear nonlinear components for *q* between −2 and +4 at the short scales, more pronounced at night, both for the heart rate and the blood pressure. Furthermore, at night, substantial nonlinear components appear in heart rate at the longer scales. It

should be noted that higher nonlinear components at night have been previously demonstrated by a noise titration procedure applied to Volterra/Wiener models fitting 24-h heart rate series [30]. The similarity of results obtained with so different approaches supports the evidence that nonlinearity prevails at night.

The finding of important nonlinear components detected by the multifractal multiscale DFA method may help designing future clinical procedures aimed at better assessing the cardiovascular risk. Actually, a recent review of complexity-based methods for the analysis of heart rate variability reported that the prediction of cardiac events by the traditional short-term coefficient of the bi-scale monofractal DFA, α1, and by the standard spectral methods are correlated [31]. This inevitably reduces the additional prediction power of α<sup>1</sup> compared to the spectral method. The traditional bi-scale monofractal model is based on the second-order moment, *q* = 2. By showing that DFA coefficients evaluated for *q* - 2 provide information substantially different from that of the spectral powers, particularly at night, our study suggests that the multifractal multiscale DFA approach might effectively integrate the information of traditional spectral methods, possibly improving the clinical value of DFA.

Finally, an unexpected pattern in the Fourier-shuffled surrogate series of Figures 5 and 6 consists in systematically higher α values for positive than for negative *q* orders (blue lines above red lines) when α increases with τ and in the opposite pattern (blue lines below red lines) when α decreases with τ. As a possible explanation of this pattern, we may hypothesize that cross-over scales appear anticipated at shorter scales when *q* > 0 and delayed at larger scales when *q* < 0.

#### **5. Limitations and Conclusions**

Nowadays, the clinical practice replaced the continuous invasive measures with intermittent noninvasive blood pressure measures for monitoring free-moving subjects, limiting the number of recordings available for the present study. Thus, it was not possible to stratify our results by gender or age, factors possibly influencing the circadian profile of the cardiovascular complexity [30]. Future studies on cardiovascular complexity can make use of noninvasive instrumentation measuring arterial blood pressure at the finger site continuously for 24 h even in ambulant subjects [32]. However, the scale coefficients of SBP could be affected by the amplification of the Mayer waves when blood pressure is measured at the digital artery [24]. Furthermore, if IBI is derived as the series of intervals between consecutive R peaks of the electrocardiogram rather than between consecutive pulses of blood pressure, as in this study, results at the shortest scales might differ because of the different amplitude of the respiratory sinus arrhythmia [33].

Finally, this study represents the temporal scales in seconds of time and not in number of beats, a relatively new methodological aspect originally proposed for comparing conditions with markedly different heart rate levels after selective autonomic blockade [6]. We adopted the same approach here because of the day–night differences in the mean heart rate and because we expected the differences to involve neural/humoral mechanisms which depend on time delays in seconds and not in number of beats (let's think to the Mayer rhythm with a 10-s period due to the slow response of vascular resistances; or to the dynamics of removal of noradrenaline released by the sympathetic nerve endings, with a time constant of 1 min; or to long-term humoral fluctuations possibly responsible for the circadian component we found at scales of about 4 min). Mapping the temporal scales from beats *n* to time τ does not change the estimate of α, still based on the calculation of the derivative of log *Fq*(*n*) vs. log *n*. This axis transformation is similar to mapping the "cycles/beat" in "Equivalent Hz" in the spectral analysis of cardiovascular series [34]. However, if results obtained with scales expressed as τ in seconds are discussed in relation to other studies based on scales defined in number of beats, readers should be aware that discrepancies may arise because possible differences in the heart rate level between conditions or groups may change the ranges of scales that define short-term and long-term DFA coefficients.

In conclusion, the multifractal multiscale DFA provides a detailed description of the complexity features of the cardiovascular series and highlights circadian modulations occurring at specific scales

and affecting the individual fractal components differently. In perspective, by focusing on the more informative portions of the α(*q*,τ) surface it could be possible to design more powerful tools for assessing the cardiovascular risk. Furthermore, coefficients with *q* - 2 reflect well the nonlinear components during night-time sleep, suggesting that they may effectively integrate the spectral information with complexity-domain information. Therefore, the evaluation of the multifractal multiscale surface of scale coefficients during wake and sleep may improve the risk assessment in cardiovascular prevention, the evaluation of cardiovascular interventions as well as the monitoring of the efficacy of rehabilitation protocols.

**Author Contributions:** G.P. and S.O. selected the patients and organized the data collection. A.F. and P.C. designed the work, created the software, analyzed the recordings, and interpreted the results. P.C. led the writing of the manuscript. A.F. made the figures. All the authors contributed to the interpretation of the data, drafted the work, and revised it. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
