# Multifractal-Multiscale Analysis of Cardiovascular Signals: A DFA-Based Characterization of Blood Pressure and Heart-Rate Complexity by Gender

Authors: Paolo Castiglioni; Davide Lazzeroni; Paolo Coruzzi; Andrea Faini
Journal: Complexity
 (2018)

---

## Abstract

Detrended Fluctuation Analysis (DFA) is a popular method for assessing the fractal characteristics of biosignals, recently adapted for evaluating the heart-rate multifractal and/or multiscale characteristics. However, the existing methods do not consider the beat-by-beat sampling of heart rate and have relatively low scale resolutions and were not applied to cardiovascular signals other than heart rate. Therefore, aim of this work is to present a DFA-based method for joint multifractal/multiscale analysis designed to address the above critical points and to provide the first description of the multifractal/multiscale structure of interbeat intervals (IBI), systolic blood pressure (SBP), and diastolic blood pressure (DBP) in male and female volunteers separately. The method optimizes data splitting in blocks to reduce the DFA estimation variance and to evaluate scale coefficients with Taylor’s expansion formulas and maps the scales from beat domains to temporal domains. Applied to cardiovascular signals recorded in 42 female and 42 male volunteers, it showed that scale coefficients and degree of multifractality depend on the temporal scale, with marked differences between IBI, SBP, and DBP and with significant sex differences. Results may be interpreted considering the distinct physiological mechanisms regulating heart-rate and blood-pressure dynamics and the different autonomic profile of males and females.


---

## Body

## 1. Introduction

Beat-by-beat measures of cardiovascular variables show an intrinsic variability, even when the cardiovascular system is observed in steady-state conditions. These spontaneous changes may reflect the processes underlying the cardiovascular homeostasis. Components of this variability show a fractal nature and in the last two decades different authors suggested that, at least for the heart rate, such components may be the output of a complex system that generates self-similar signals [1, 2]. In fact, the cardiovascular system, like several complex dynamical systems, is composed of interacting subsystems embedded in a fractal structure. In particular, fractal networks of vessels and of nervous and humoral pathways connect and hierarchically regulate local blood flows among several vascular beds. Accordingly, the cardiovascular system can be regarded as a dissipative system that preserves homeostasis evolving toward a self-organized state, not characterized by any intrinsic scale of time [3]. Properly assessing the fractal components of the spontaneous variability of cardiovascular signals is important, because it may help identifying early alterations in cardiovascular regulatory mechanisms and may contribute to stratifying more precisely the cardiovascular risk.The first descriptions of cardiovascular self-similarity were based on modeling the heart rate as a time series belonging to the families of fractional Gaussian noises or of fractional Brownian motions and on estimating the corresponding Hurst exponent [4]. Successive studies recognized that such an approach oversimplifies a more complex phenomenon, because the fractal characteristics of heart rate appear to depend on the scale of the observation [5]. This led some authors to propose multiscale approaches that quantify the cardiovascular complexity by a spectrum of self-similarity coefficients evaluated at different temporal scales [6–8]. Other authors provided evidence of the multifractal nature of heart rate [9–11], which means that the self-similar components of variability result from the superimposition of different fractal processes, interwoven at the same scales. This makes the inadequacy of methods based on monofractal models to describe the cardiovascular complexity even more apparent.For these reasons, the more recent research in the field of heart-rate variability is aimed at proposing methods that take into account both the scale dependency of self-similarity and its multifractal nature [12, 13]. Following this line of research, the aim of the present study is to describe the multifractal and multiscale characteristics of cardiovascular signals in healthy subjects under controlled conditions. This is done by adapting previously proposed methods of multifractal and multiscale analysis and by comparing three cardiovascular signals frequently recorded in physiological and clinical studies: interbeat interval (IBI, inverse of heart rate), systolic blood pressure (SBP), and diastolic blood pressure (DBP). We expect different fractal dynamics for these three signals because they are influenced by different cardiovascular effectors: DBP is mainly modulated by changes in vascular resistances, SBP by changes in cardiac output, and IBI by changes in cardiac outflows of the autonomic nervous system. Since males and females are characterized by a different autonomic profile [14], the analysis also focuses on sex differences in the multifractal and multiscale dynamics.


## 2. Methods

### 2.1. Multifractal-Multiscale DFA

Our estimator of the multifractal-multiscale characteristics of beat-by-beat cardiovascular signals was based on detrended fluctuation analysis (DFA), a method originally proposed for calculating a scale exponent,α, strictly related to Hurst’s exponent of monofractal time series [5]. DFA has been successively extended to analyze multifractal time series, obtaining distributions of α coefficients that describe the superposition of different fractal processes [15]. Moreover, DFA has been also extended to provide multiscale evaluations, that is, a spectrum of α coefficients function of the observation scale [6, 16]. Therefore, DFA is a versatile technique easily adaptable for multifractal or multiscale analysis. In this regard, recently Gierałtowski et al. combined both the approaches: they proposed a multifractaland multiscale method for the DFA of heart-rate variability, exploiting the possibility of adapting the multifractal DFA algorithm in order to provide estimates separately at different scales [12]. This method was recently applied for modeling heart-rate variability during sleep and blood-pressure variability [17, 18]. In the present study, we followed a similar approach, introducing, however, important variants. These took into account specific properties of the beat-to-beat cardiovascular dynamics that regard the way local slopes are derived and their proper mapping in the time domain.Given a time series of the cardiovascular variablex(j) with mean μ, evaluated over N consecutive heart beats (1≤j≤N), its cumulative sum (1)yi=∑j=1ixj-μwas calculated for 1≤i≤N. Then, fixing a block sizen in number of beats, y(i) was split intoM blocks each containing a data segment of n beats. Often DFA is evaluated considering nonoverlapping consecutive blocks: in this case, M=int⁡(N/n) and a short segment ofN-nM data at the end of the series is not included in any of theM blocks ifN is not multiple ofn [19]. By contrast, we overlapped consecutive segments so that two successive blocks had n-1 beats in common (maximal overlapping),M was equal to (N-n+1) and all the data were included in at least one block for any sizen. As illustrated in Figure 1, maximal overlapping reduces the estimator variance substantially (this will allow evaluating the local slopes α with numerical differentiation formulas). Data were detrended in each of theM blocks with a least-square polynomial fitting of order 1. The standard deviation of the detrended data was calculated in each block k, σn(k), for 1≤k≤M.Figure 1
Examples of multifractal
F
q
(
n
)
functions plotted versus the block size n at different q values for nonoverlapped (a), 50% overlapped (b), and maximally overlapped blocks (c). Data from beat-to-beat IBI series of 2-hour duration recorded in a healthy volunteer sitting at rest.
(a)
(b)
(c)Data splitting was repeated for block sizesn between 6 and N/4 beats. Block sizes were selected as the closest integers to a distribution evenly spaced on a logarithmic scale, with density of about 13 samples per decade. For instance, for N=8400 beats, corresponding to a 2-hour recording at the heart rate of 70 bpm, we considered 34 block sizes n between 6 and 1827 beats. Because of the low size of the smallest block, we set the order of the fitting polynomial equal to 1, to avoid overfitting the data with a too high order, which might remove not only the trend but also the significant components of variability.According to the multifractal approach for DFA [15], a family of variability functions, Fq(n), which depend on the multifractal parameter q, are calculated for each block size n, as(2)Fqn=1M∑k=1Mσn2kq/21/qforq≠0Fqn=e1/2M∑k=1Mln⁡σn2kforq=0.If x(i) has a power law correlation—like fractional Gaussian noises or fractional Brownian motions—then Fq(n) increases as a power ofn for any choice of the parameter q:Fq(n)∝nα. For monofractal time series, the exponent α is associated with Hurst’s exponentH, being α=H for fractional Gaussian noises and α=H+1 for fractional Brownian motions. By contrast, if x(i) is a multifractal series with fractal components of different amplitude, α reflects the superposition of different power law correlations. In particular, α mainly reflects the fractal components with larger amplitude if q>0 and the fractal components with smaller amplitude if q<0. Therefore, α coefficients that depend onq are sign of multifractal dynamics.The multifractal exponentα(q) can be estimated from Fq(n) in (2) as the slope of a least-square linear regression between log⁡Fq(n) and log⁡n [20]. In this way, however, local deviations from the linear trend occurring at specific scalesn cannot be detected. The assessment of local deviations may reflect changes in the sympathetic and vagal cardiac control not otherwise visible [7, 16], revealing subtle alterations in the overall autonomic regulation of the cardiovascular system [21] and characterizing pathological conditions [6, 22]. To evaluate a local slope, that is, α as function ofn, methods with higher scale resolution are required. A simple way to obtain a multiscale representation is to calculate the least-square linear regression over a running window with constant width over the logn axis [12]. The estimated slope, α(q,n), is associated with the central scale n of the running window. However, the length of the running window influences the “smoothness” of the α(q,n) curves and limits the range of scales where α is estimated. Alternatively, α(q,n) could be estimated as derivative of log⁡Fq(n) versus log⁡n, as proposed for monofractal DFA [16]. We followed this approach and sincen was approximately spaced evenly on the logarithmic scale, we applied formula derived from Taylor’s expansion. Let us call {nl} with 1≤l≤lMAX the set of lMAX block sizes, where we calculated Fq(n) in (2). The 3-point expression of the derivative of log⁡Fq(n) versus log⁡n is(3)αBq,nl=log⁡Fqnl+1-log⁡Fqnl-1log⁡nl+1-log⁡nl-1.In (3), the pedixB of αB means that the scale coefficient is evaluated on the beat domain, n. For l=1 and l=lMAX (3) is not defined, and we used the expressions for right and left derivatives:(4)αBq,nl=-log⁡Fqnl+2+4log⁡Fqnl+1-3log⁡Fqnllog⁡nl+2-log⁡nlforl=1,αBq,nl=log⁡Fqnl-2-4log⁡Fqnl-1+3log⁡Fqnllog⁡nl-log⁡nl-2forl=lMAX.Equations (3) and (4) approximate the first derivative of log Fq(n) versus logn with errors proportional to the amplitude of the derivatives of order higher than 2. A better approximation is provided by the formula on 5 points with errors proportional to derivatives of order higher than 3. Therefore, instead of (3), for 2<l<lMAX-1 we used the following equation:(5)αBq,nl=8log⁡Fqnl+1-log⁡Fqnl-1-log⁡Fqnl+2-log⁡Fqnl-23log⁡nl+2-log⁡nl-2.Figure 2 illustrates how (3)–(5) derive αB(q,n) in a real case.Figure 2
Example of local slopes estimation. (a) shows the same Fq(n) function plotted for maximally overlapped blocks and q=-2 in Figure 1; the straight lines spanning over 5 points centered around the diamond symbols at n=20, n=40, n=161, and n=646 represent the local slopes α(n) evaluated by the 5-point derivative of (5). (b) shows the corresponding spectrum of local slopes.
(a)
(b)SinceαB is a function of the scale expressed in number of beats, when series with different mean heart rate are compared (e.g., a bradycardic versus a tachycardic subject or rest versus exercise conditions) the same scalesn, in beats, correspond to different temporal scales, τ, in seconds. Therefore, to associate each scale coefficient with its temporal scale, we mapped the beat domain into the time domain [7]. Given the αB(q,nl) coefficients evaluated on the beat domain for the {nl} set of scales, the coefficients evaluated on the corresponding set of temporal scales, {τl}, are(6)αq,τl=αBq,nlforτl=nl×μIBIwith μIBI being the mean IBI, in seconds.Theα(q,τl) coefficients were interpolated over the τ axis with a spline function to obtain estimates at the same temporal scales for each recording. On the basis of the analysis of synthesized series with known self-similarity structure (see Appendix A), for our application on healthy volunteers (see Section 2.3) we interpolated 256 points evenly spaced over the logarithmic τ axis, between τ=8 s and τ=512 s when q>-3 and between τ=10 s and τ=512 s when q≤-3. The largest scale (τ=512 s) corresponds to less than 10% the average duration of the recordings. Estimates were obtained forq between −5 and +5, with incremental step of 0.5. Figure 3 shows an example of estimated α(q,τ) coefficients.Figure 3
Example of multifractal and multiscale coefficients, α(q,τ), for IBI, SBP, and DBP. Data recorded in a male, normotensive volunteer. In this example, α(q,τ) of IBI decreases steeply with τ at scales shorter than 16 s, from values typical of fractional Brownian motions (α>1) to values typical of fractional Gaussian noises (α<1). Even if α is lower than 1 at larger scales, the α(q,τ) surface is not flat and a relative maximum appears at τ=32 s, more pronounced for positive rather than for negativeq values. SBP differs remarkably from IBI: α tends to be greater at scales > 64 s, and pronounced local maxima appear when q<0. DBP differs from IBI and SBP: α(q,τ) increases at τ>128 s for anyq, reaching values greater than 1; and no local maxima appear at τ<128 s.Finally, we defined a concise index of multifractality, function of the scaleτ:MFI(τ). For this purpose, fixing a parameter qr>0, we considered the range -qr≤q≤qr, symmetric around 0 with amplitude 2qr. For each temporal scale τ, we calculated the standard deviation of all α(q,τ) values estimated over the ±qr range, αSD(τ). The MFI(τ) index is defined as the ratio between αSD(τ) and the range of correspondingq values: (7)MFIτ=αSDτ2qr.Since αSD(τ)≥0, also MFI(τ)≥0, reaching values close to 0 if the series is monofractal at τ (Appendix A shows an example of MFI(τ) for synthesized monofractal series). In this study, MFI(τ) was calculated for τ between 10 s and 512 s, setting qr=5.


### 2.2. Power Spectral Analysis

The power spectrum of each cardiovascular series was also calculated. Beat-to-beat series were interpolated linearly at 10 Hz and resampled at 5 Hz. The Welch periodogram was estimated by splitting the resampled series in 50% overlapping Hann windows of 1638.4 s duration, by computing the FFT spectrum in each window and by averaging the spectra over all the windows. The final periodogram was smoothed with a broadband procedure [23].


### 2.3. Subjects and Experimental Protocol

We considered recordings previously collected in two studies aimed at evaluating the influence of sodium sensitivity on the cardiovascular control in normotensive [24] and hypertensive [25] healthy subjects. The original dataset in normotensive subjects consisted of recordings in 26 males and 45 females [24]. For the present analysis, we included all the 26 male participants and a subgroup of 26 female participants matched for age and body mass index. The original dataset in hypertensive subjects consisted of recordings in 30 males and 16 females, and, for the present analysis, we included all the 16 female participants and a subgroup of 16 male participants with age and body mass index matched with the female group. Table 1 summarizes by sex the general characteristics of the selected 84 participants.Table 1
General characteristics of participants by sex.
N
Age (years)
Body mass index (kg/m2)
Prevalence of hypertension
Females
42
34.3 (9.7)
22.6 (2.8)
38.1%
Males
42
34.4 (10.1)
23.6 (2.4)
38.1%
p
0.96
0.09Each participant was studied in a quiet environment in the morning, after 5 days of low-salt diet (30 mmol NaCl per day) to minimize the confounding effects of dietary sodium on cardiovascular variability. Continuous finger arterial blood pressure was recorded for about two hours, in sitting position at rest, by Portapres model-2 (Finapres Medical Systems B. V., Amsterdam, Netherlands). The finger cuff was placed on the mid finger of the left hand. SBP, DBP, and IBI (calculated as time interval between consecutive SBP values) were derived beat-by-beat for the whole duration of the recording. Brachial blood pressure was measured simultaneously with a cuff on the right arm every 15 minutes, and the SBP and DBP readings of the brachial device were used to calibrate beat-by-beat SBP and DBP values from the finger cuff.


### 2.4. Statistics

We described statistical patterns inα(q,τ) estimates showing means and standard error of the means over the group, on the basis of previous observations reporting that DFA coefficients follow a normal distribution [26] and in MFI(τ) estimates showing median and standard error of the median, this latter estimated by bootstrapping using 100 bootstrap samples. Statistical inferences were performed with nonparametric tests for all the estimates not to make any assumption on the distribution of scale coefficients and multifractal indices at any τ. In particular, α(q,τ) coefficients were compared between signals (IBI versus SBP, IBI versus DBP and SBP versus DBP) by the paired Wilcoxon test; α(q,τ) and MFI(τ) were compared between males and females by the unpaired Mann–Whitney test. Power spectra were compared between genders by unpaired t-test after log-transformation, to obtain normal distributions of power spectra [27]. The analyses were performed with “R: A Language and Environment for Statistical Computing” software package (R Core Team, R Foundation for Statistical Computing, Vienna, Austria, 2017).


## 2.1. Multifractal-Multiscale DFA

Our estimator of the multifractal-multiscale characteristics of beat-by-beat cardiovascular signals was based on detrended fluctuation analysis (DFA), a method originally proposed for calculating a scale exponent,α, strictly related to Hurst’s exponent of monofractal time series [5]. DFA has been successively extended to analyze multifractal time series, obtaining distributions of α coefficients that describe the superposition of different fractal processes [15]. Moreover, DFA has been also extended to provide multiscale evaluations, that is, a spectrum of α coefficients function of the observation scale [6, 16]. Therefore, DFA is a versatile technique easily adaptable for multifractal or multiscale analysis. In this regard, recently Gierałtowski et al. combined both the approaches: they proposed a multifractaland multiscale method for the DFA of heart-rate variability, exploiting the possibility of adapting the multifractal DFA algorithm in order to provide estimates separately at different scales [12]. This method was recently applied for modeling heart-rate variability during sleep and blood-pressure variability [17, 18]. In the present study, we followed a similar approach, introducing, however, important variants. These took into account specific properties of the beat-to-beat cardiovascular dynamics that regard the way local slopes are derived and their proper mapping in the time domain.Given a time series of the cardiovascular variablex(j) with mean μ, evaluated over N consecutive heart beats (1≤j≤N), its cumulative sum (1)yi=∑j=1ixj-μwas calculated for 1≤i≤N. Then, fixing a block sizen in number of beats, y(i) was split intoM blocks each containing a data segment of n beats. Often DFA is evaluated considering nonoverlapping consecutive blocks: in this case, M=int⁡(N/n) and a short segment ofN-nM data at the end of the series is not included in any of theM blocks ifN is not multiple ofn [19]. By contrast, we overlapped consecutive segments so that two successive blocks had n-1 beats in common (maximal overlapping),M was equal to (N-n+1) and all the data were included in at least one block for any sizen. As illustrated in Figure 1, maximal overlapping reduces the estimator variance substantially (this will allow evaluating the local slopes α with numerical differentiation formulas). Data were detrended in each of theM blocks with a least-square polynomial fitting of order 1. The standard deviation of the detrended data was calculated in each block k, σn(k), for 1≤k≤M.Figure 1
Examples of multifractal
F
q
(
n
)
functions plotted versus the block size n at different q values for nonoverlapped (a), 50% overlapped (b), and maximally overlapped blocks (c). Data from beat-to-beat IBI series of 2-hour duration recorded in a healthy volunteer sitting at rest.
(a)
(b)
(c)Data splitting was repeated for block sizesn between 6 and N/4 beats. Block sizes were selected as the closest integers to a distribution evenly spaced on a logarithmic scale, with density of about 13 samples per decade. For instance, for N=8400 beats, corresponding to a 2-hour recording at the heart rate of 70 bpm, we considered 34 block sizes n between 6 and 1827 beats. Because of the low size of the smallest block, we set the order of the fitting polynomial equal to 1, to avoid overfitting the data with a too high order, which might remove not only the trend but also the significant components of variability.According to the multifractal approach for DFA [15], a family of variability functions, Fq(n), which depend on the multifractal parameter q, are calculated for each block size n, as(2)Fqn=1M∑k=1Mσn2kq/21/qforq≠0Fqn=e1/2M∑k=1Mln⁡σn2kforq=0.If x(i) has a power law correlation—like fractional Gaussian noises or fractional Brownian motions—then Fq(n) increases as a power ofn for any choice of the parameter q:Fq(n)∝nα. For monofractal time series, the exponent α is associated with Hurst’s exponentH, being α=H for fractional Gaussian noises and α=H+1 for fractional Brownian motions. By contrast, if x(i) is a multifractal series with fractal components of different amplitude, α reflects the superposition of different power law correlations. In particular, α mainly reflects the fractal components with larger amplitude if q>0 and the fractal components with smaller amplitude if q<0. Therefore, α coefficients that depend onq are sign of multifractal dynamics.The multifractal exponentα(q) can be estimated from Fq(n) in (2) as the slope of a least-square linear regression between log⁡Fq(n) and log⁡n [20]. In this way, however, local deviations from the linear trend occurring at specific scalesn cannot be detected. The assessment of local deviations may reflect changes in the sympathetic and vagal cardiac control not otherwise visible [7, 16], revealing subtle alterations in the overall autonomic regulation of the cardiovascular system [21] and characterizing pathological conditions [6, 22]. To evaluate a local slope, that is, α as function ofn, methods with higher scale resolution are required. A simple way to obtain a multiscale representation is to calculate the least-square linear regression over a running window with constant width over the logn axis [12]. The estimated slope, α(q,n), is associated with the central scale n of the running window. However, the length of the running window influences the “smoothness” of the α(q,n) curves and limits the range of scales where α is estimated. Alternatively, α(q,n) could be estimated as derivative of log⁡Fq(n) versus log⁡n, as proposed for monofractal DFA [16]. We followed this approach and sincen was approximately spaced evenly on the logarithmic scale, we applied formula derived from Taylor’s expansion. Let us call {nl} with 1≤l≤lMAX the set of lMAX block sizes, where we calculated Fq(n) in (2). The 3-point expression of the derivative of log⁡Fq(n) versus log⁡n is(3)αBq,nl=log⁡Fqnl+1-log⁡Fqnl-1log⁡nl+1-log⁡nl-1.In (3), the pedixB of αB means that the scale coefficient is evaluated on the beat domain, n. For l=1 and l=lMAX (3) is not defined, and we used the expressions for right and left derivatives:(4)αBq,nl=-log⁡Fqnl+2+4log⁡Fqnl+1-3log⁡Fqnllog⁡nl+2-log⁡nlforl=1,αBq,nl=log⁡Fqnl-2-4log⁡Fqnl-1+3log⁡Fqnllog⁡nl-log⁡nl-2forl=lMAX.Equations (3) and (4) approximate the first derivative of log Fq(n) versus logn with errors proportional to the amplitude of the derivatives of order higher than 2. A better approximation is provided by the formula on 5 points with errors proportional to derivatives of order higher than 3. Therefore, instead of (3), for 2<l<lMAX-1 we used the following equation:(5)αBq,nl=8log⁡Fqnl+1-log⁡Fqnl-1-log⁡Fqnl+2-log⁡Fqnl-23log⁡nl+2-log⁡nl-2.Figure 2 illustrates how (3)–(5) derive αB(q,n) in a real case.Figure 2
Example of local slopes estimation. (a) shows the same Fq(n) function plotted for maximally overlapped blocks and q=-2 in Figure 1; the straight lines spanning over 5 points centered around the diamond symbols at n=20, n=40, n=161, and n=646 represent the local slopes α(n) evaluated by the 5-point derivative of (5). (b) shows the corresponding spectrum of local slopes.
(a)
(b)SinceαB is a function of the scale expressed in number of beats, when series with different mean heart rate are compared (e.g., a bradycardic versus a tachycardic subject or rest versus exercise conditions) the same scalesn, in beats, correspond to different temporal scales, τ, in seconds. Therefore, to associate each scale coefficient with its temporal scale, we mapped the beat domain into the time domain [7]. Given the αB(q,nl) coefficients evaluated on the beat domain for the {nl} set of scales, the coefficients evaluated on the corresponding set of temporal scales, {τl}, are(6)αq,τl=αBq,nlforτl=nl×μIBIwith μIBI being the mean IBI, in seconds.Theα(q,τl) coefficients were interpolated over the τ axis with a spline function to obtain estimates at the same temporal scales for each recording. On the basis of the analysis of synthesized series with known self-similarity structure (see Appendix A), for our application on healthy volunteers (see Section 2.3) we interpolated 256 points evenly spaced over the logarithmic τ axis, between τ=8 s and τ=512 s when q>-3 and between τ=10 s and τ=512 s when q≤-3. The largest scale (τ=512 s) corresponds to less than 10% the average duration of the recordings. Estimates were obtained forq between −5 and +5, with incremental step of 0.5. Figure 3 shows an example of estimated α(q,τ) coefficients.Figure 3
Example of multifractal and multiscale coefficients, α(q,τ), for IBI, SBP, and DBP. Data recorded in a male, normotensive volunteer. In this example, α(q,τ) of IBI decreases steeply with τ at scales shorter than 16 s, from values typical of fractional Brownian motions (α>1) to values typical of fractional Gaussian noises (α<1). Even if α is lower than 1 at larger scales, the α(q,τ) surface is not flat and a relative maximum appears at τ=32 s, more pronounced for positive rather than for negativeq values. SBP differs remarkably from IBI: α tends to be greater at scales > 64 s, and pronounced local maxima appear when q<0. DBP differs from IBI and SBP: α(q,τ) increases at τ>128 s for anyq, reaching values greater than 1; and no local maxima appear at τ<128 s.Finally, we defined a concise index of multifractality, function of the scaleτ:MFI(τ). For this purpose, fixing a parameter qr>0, we considered the range -qr≤q≤qr, symmetric around 0 with amplitude 2qr. For each temporal scale τ, we calculated the standard deviation of all α(q,τ) values estimated over the ±qr range, αSD(τ). The MFI(τ) index is defined as the ratio between αSD(τ) and the range of correspondingq values: (7)MFIτ=αSDτ2qr.Since αSD(τ)≥0, also MFI(τ)≥0, reaching values close to 0 if the series is monofractal at τ (Appendix A shows an example of MFI(τ) for synthesized monofractal series). In this study, MFI(τ) was calculated for τ between 10 s and 512 s, setting qr=5.


## 2.2. Power Spectral Analysis

The power spectrum of each cardiovascular series was also calculated. Beat-to-beat series were interpolated linearly at 10 Hz and resampled at 5 Hz. The Welch periodogram was estimated by splitting the resampled series in 50% overlapping Hann windows of 1638.4 s duration, by computing the FFT spectrum in each window and by averaging the spectra over all the windows. The final periodogram was smoothed with a broadband procedure [23].


## 2.3. Subjects and Experimental Protocol

We considered recordings previously collected in two studies aimed at evaluating the influence of sodium sensitivity on the cardiovascular control in normotensive [24] and hypertensive [25] healthy subjects. The original dataset in normotensive subjects consisted of recordings in 26 males and 45 females [24]. For the present analysis, we included all the 26 male participants and a subgroup of 26 female participants matched for age and body mass index. The original dataset in hypertensive subjects consisted of recordings in 30 males and 16 females, and, for the present analysis, we included all the 16 female participants and a subgroup of 16 male participants with age and body mass index matched with the female group. Table 1 summarizes by sex the general characteristics of the selected 84 participants.Table 1
General characteristics of participants by sex.
N
Age (years)
Body mass index (kg/m2)
Prevalence of hypertension
Females
42
34.3 (9.7)
22.6 (2.8)
38.1%
Males
42
34.4 (10.1)
23.6 (2.4)
38.1%
p
0.96
0.09Each participant was studied in a quiet environment in the morning, after 5 days of low-salt diet (30 mmol NaCl per day) to minimize the confounding effects of dietary sodium on cardiovascular variability. Continuous finger arterial blood pressure was recorded for about two hours, in sitting position at rest, by Portapres model-2 (Finapres Medical Systems B. V., Amsterdam, Netherlands). The finger cuff was placed on the mid finger of the left hand. SBP, DBP, and IBI (calculated as time interval between consecutive SBP values) were derived beat-by-beat for the whole duration of the recording. Brachial blood pressure was measured simultaneously with a cuff on the right arm every 15 minutes, and the SBP and DBP readings of the brachial device were used to calibrate beat-by-beat SBP and DBP values from the finger cuff.


## 2.4. Statistics

We described statistical patterns inα(q,τ) estimates showing means and standard error of the means over the group, on the basis of previous observations reporting that DFA coefficients follow a normal distribution [26] and in MFI(τ) estimates showing median and standard error of the median, this latter estimated by bootstrapping using 100 bootstrap samples. Statistical inferences were performed with nonparametric tests for all the estimates not to make any assumption on the distribution of scale coefficients and multifractal indices at any τ. In particular, α(q,τ) coefficients were compared between signals (IBI versus SBP, IBI versus DBP and SBP versus DBP) by the paired Wilcoxon test; α(q,τ) and MFI(τ) were compared between males and females by the unpaired Mann–Whitney test. Power spectra were compared between genders by unpaired t-test after log-transformation, to obtain normal distributions of power spectra [27]. The analyses were performed with “R: A Language and Environment for Statistical Computing” software package (R Core Team, R Foundation for Statistical Computing, Vienna, Austria, 2017).


## 3. Results

Figure4 shows α as a function of τ for specificq values over the whole group (for comparison, the traditional multiscale analysis corresponds to α values evaluated for q=2 only, and the traditional multifractal analysis corresponds to the generalized Hurst exponents shown in Appendix B). Figure 4 confirms patterns suggested in the example of Figure 3. IBI coefficients decrease with τ from values greater than 1 (as for fractional Brownian motions) at the shorter scales to values lower than 1 (as for fractional Gaussian noises) at the larger scales, with a minimum at τ around 250 s. A relative maximum appears at τ≅30 s for q≥2. Moreover, α increases asq decreases, at any τ. At the shorter scales also α of SBP decreases steeply with τ from values >1, with greater α estimates at lowerq values. However, unlike α of IBI, it remains stable around 1 (as for “1/f” processes) when τ>30 s. Similarly to IBI and to SBP, also α of DBP decreases with τ at the shorter scales. However, unlike IBI and SBP, at larger scales it shows an increasing trend with τ.Figure 4
Multifractal and multiscale coefficients,
α
(
q
,
τ
). Mean ± standard error of the mean over the group (N=84) as function of τ for fiveq values between −4 and +4. For clarity, the horizontal axis is plotted in a log scale.Figure5 compares scale coefficients among signals, at different q. Comparing IBI with SBP when q=0 and q=4, α is significantly greater for SBP almost at all the scales; however, this is not the case for scales τ between 16 and 35 s, where α of IBI shows a local maximum. When q=-4, differences between IBI and SBP are less significant, and their scale coefficients coincide at scale τ≤16 s.Figure 5
Comparison of
α
(
q
,
τ
)
among cardiovascular signals at different q values. (a) Mean ± standard error of the mean over the group (N=84) as function of τ. (b),from top to bottom: W statistics of Wilcoxon’s tests at each scale τ for comparing IBI and SBP, IBI and DBP, and SBP and DBP. Dashed horizontal lines represent the 5% threshold of statistical significance; W values above the threshold are marked by a red dot and indicate statistically significant differences at p<5%.
(a)
(b)Similarly, comparing IBI and DBP,α is greater for DBP almost at all the scales when q=0 and q=4. In this case, however, at the scales where α of IBI displays a relative maximum, α of DBP shows an absolute minimum, becoming as a result significantly lower than the IBI scale coefficients.Also SBP and DBP scale coefficients differ importantly. Forq=0 and q=4, α of SBP is significantly greater within 20<τ<100 s. Similar differences appear for q=-4 but in a slightly higher band (45<τ<240 s). At τ≤16 s, α is greater for DBP when q=-4 while when q=+4, SBP and DBP scale coefficients coincide.These results make it clear that the degree of multifractality is a function ofτ and of the type of cardiovascular signal. This is summarized by Figure 6, which shows the index of multifractality, MFI(τ), for IBI, SBP, and DBP separately. IBI reaches its highest degree of multifractality at τ=10 s; MFI(τ) of IBI decreases at larger scales up to a minimum at τ=30 s. A different pattern characterizes MFI(τ) of DBP and SBP: the highest degree of multifractality is not reached at the shortest scale (τ=10 s) as for IBI but between 16 and 32 s; and the lowest degree is reached at τ=64 s.Figure 6
Multifractal index
M
F
I
(
τ
)
for IBI, SBP, and DBP. Median ± standard error of the median over the group (N=84).Gender Differences. Figure 7 compares α(q,τ) for q=-4, q=0, and q=+4, in males and females. When q=+4 or q=0, α is greater in females at the larger scales (i.e., τ>60 s for IBI and SBP; τ>20 s for DBP). These differences vanish when q=-4. In addition, when q=0 and q=-4, α is lower in females at τ ≤ 10 s. The fact that differences between males and females depend onq suggests sex-related differences also in the level of multifractality: these are actually highlighted by Figure 8, which shows a higher degree of multifractality in males, at scales shorter than 16 s for IBI and at scales centered around 20 s for SBP and 32 s for DBP.Figure 7
Multifractal-multiscale coefficients,
α
(
q
,
τ
)
: comparison by sex. Mean ± standard error of the mean(upper panels) and W statistics of Mann–Whitney tests(lower panels) for comparing males (N=42) and females (N=42), as function of τ. Dashed horizontal lines represent the 5% significance threshold; W values above the threshold (red dots) indicate differences significant at p<5%. From top to bottom: q=-4; q=0; q=+4.Figure 8
Multifractal index
M
F
I
(
τ
)
: comparison by sex. (a) Median ± standard error of the median; (b) W statistics of Mann–Whitney tests for comparing males (N=42) and females (N=42), as function of τ. Dashed horizontal lines represent the 5% significance threshold; W values above the threshold (red dots) indicate differences significant at p<5%.
(a)
(b)Figure9 compares IBI, SBP, and DBP power spectra by gender. The three signals have common spectral patterns, all showing a peak around 0.10 Hz and a “1/f” component at frequencies lower than 0.03 Hz; moreover, IBI and SBP spectra also show a respiratory component at frequencies around 0.30 Hz. Although these patterns appear in both sexes, spectra differ significantly between males and females. IBI spectral powers are greater in males between 0.008 and 0.13 Hz. Also SBP spectral components are greater in males, but in a larger band including all frequencies higher than 0.004 Hz. By contrast, DBP spectra coincide in men and women, with exclusion of the spectral peak around 0.10 Hz, higher in men.Figure 9
Power spectra of IBI, SBP, and DBP: comparison by sex. (a) Mean ± standard error of the mean; (b) Student’s t statistics at each frequency for comparing males (N=42) versus females (N=42). Dashed horizontal lines represent the 5% significance threshold; t values above the threshold (red dots) indicate differences significant at p<5%.
(a)
(b)


## 4. Discussion

We presented a novel algorithm for quantifying cardiovascular complexity based on previous researches that in various ways adapted DFA for assessing multifractal or/and multiscale aspects of heart-rate variability. By applying our method to data collected in healthy volunteers, we provided the first detailed description of differences in multifractal and multiscale features among those cardiovascular time series more often recorded in clinical settings or in physiological studies: IBI, SBP, and DBP.Three were the main results of our study. First, not only do self-similarity coefficients depend on the observational scale, but also the wayα changes with τ depends on the cardiovascular series (IBI, SBP, or DBP). Second, the degree of multifractality also depends on τ and on the type of cardiovascular signal. Third, at the scales where the signals show a multifractal nature (e.g., at τ<32 s), increasing the multifractal indexq from negative to positive values progressively decreases the estimate of α, thus indicating that fractal components with the lower α contribute more to the overall variability because positiveq values emphasize components with larger amplitude. In addition, we also showed significant gender differences in these complex self-similarity structures. Although this study was designed to provide a solid description of the multifractal and multiscale features of cardiovascular signals and not to find the mechanisms responsible for these features, in the following we may try to speculate on the possible origin and physiological meaning of our results.We showed thatα(q,τ) of IBI decreases with τ in healthy individuals sitting at rest, from values typical for fractional Brownian motions to values typical for fractional Gaussian noises. A similar behavior has been previously observed in other studies that applied DFA with traditional monofractal approaches (i.e., with q=2 only) to analyze the heart-rate variability in sitting volunteers at rest [16, 28]. A possible explanation for this trend has been already proposed, based on the hypothesis that the heart-rate dynamics depend on the superposition of two fractal processes simultaneously modulating the heart rate [7, 29]. One process, with relatively homogeneous fractal characteristics at all scales, resembling fractional Gaussian noise, depends on the cardiac vagal outflow. The other process depends on the cardiac sympathetic outflow: it appears as a fractional Gaussian noise at the longer scales and as a Brownian motion at the shortest scales, due to a low-pass filter with time constant of 66 s caused by the time of removal of noradrenaline released by the sympathetic endings. The mixture of these two processes may therefore explain the α decreasing trend with τ from fractional Brownian motion to fractional Gaussian noise and may explain why α, when evaluated at the shortest scales only, can be considered index of sympathovagal balance [30]. However, if this hypothesis holds, it should predict some specific features of the multifractal dynamics of heart rate. First, we should expect that the heart rate is multifractal at the shortest scales, where it depends on two processes with different dynamics (fractional Gaussian noise for the vagal modulations, Brownian motion for the sympathetic modulations) and with different amplitude (greater for vagal than for sympathetic modulations). Second, we should expect that multifractality decreases at scales τ larger and closer to the time constant of the low-pass filter on the sympathetic outflow, because amplitude and fractal dynamics of sympathetic heart-rate modulations should become more similar to amplitude and fractal dynamics of vagal modulations of heart rate. Third, we should also expect that, at scales shorter than the time constant of the low-pass filter modeling the sympathetic outflow dynamics (= 66 s), α increases when q decreases. In fact, negative q exponents in (2) amplify the contribution of the fractal components with lower amplitude and reduce the contribution of the fractal components with higher amplitude in multifractal dynamics. Since at the shorter scales the sympathetic modulations of heart rate have lower amplitude than the vagal heart-rate modulations, negativeq values should emphasize the Brownian motion contribution of the sympathetic outflow (with high α) rather than the fractal noise contribution of the vagal outflow (with low α). These three properties are actually demonstrated by our results on α(q,τ) of IBI, which therefore support our interpretative hypothesis.Interestingly, we found gender differences inα(q,τ) of IBI at τ≤12 s but only for q=0 and q=-4. When q=+4, in fact, α coefficients are the same in males and females over these scales (Figure 7). The IBI power spectra (Figure 9) indicate that males and females have the same amplitude for spectral components falling in the high-frequency band (>0.15 Hz) where only vagal modulations of heart rate are present and that males have higher amplitude of spectral components in the low-frequency band (around 0.1 Hz) where the heart rate is modulated importantly by both sympathetic and vagal outflows [31]. This means that, among our volunteers, males and females have a similar vagal tone, but males have a higher sympathovagal balance. In our interpretative hypothesis, α at the shorter scales represents mainly the vagal fractal modulation whenq is high, and whenq decreases the sympathetic contribution to α increases. Coherently, we observed that α at the shorter scales is greater in males, as expected for an index of sympathovagal balance, but only for low values ofq. No differences are found for q=+4, when we expect α to represent mainly the vagal component of the fractal dynamics, in line with the observation that the high-frequency spectral power associated with pure vagal modulations of heart rate are the same in males and females. The presence of similar vagal modulations of heart rate in males and females and of a higher sympathetic tone in males would also explain the higher IBI multifractal index in males at the shorter scales (Figure 8).We found marked differences betweenα(q,τ) of IBI and DBP (Figure 5). While α(q,τ) of IBI should reflect a mixture of fractal processes produced by the autonomic nervous system through vagal and sympathetic modulations of heart rate, α(q,τ) of DBP is expected to mainly reflect modulations of total peripheral resistances, which the autonomic nervous system modulates through vascular sympathetic outflows only. This would explain why the highest degree of multifractality occurs at the shortest scale for IBI, where the two branches of the cardiac autonomic nervous system modulate the heart rate with different fractal dynamics. By contrast, DBP shows the highest degree of multifractality between 16 and 32 s (Figure 6). A possible explanation is that the total peripheral resistance results from the combined effect of individual vascular resistances of several vascular beds, each with its own local regulation hierarchically controlled at the central level through vascular sympathetic outflows. Therefore, at scales around 16 and 32 s the DBP multifractal dynamics would reflect the superimposition of various local regulations associated with different vascular districts.A final consideration regards the different physiological information derivable from traditional spectral analyses and from a multifractal-multiscale analysis of cardiovascular signals. We found gender differences inα(q,τ) but not in power spectra at the corresponding frequencies, and vice versa we found power spectral differences without significant differences in the fractal structure at the corresponding scales. For instance, α(q,τ) of DBP is greater in females for τ between 32 and 256 s and q=0 or q=+4, but at the corresponding frequencies (between 0.031 and 0.004 Hz) the DBP spectra are exactly the same in males and females. On the other hand, α(q,τ) of SBP is very similar in males and females at τ<64 s, while the SBP spectra differ markedly by gender at frequencies higher than 0.01 Hz. This means that complexity methods and spectral methods are complementary approaches and that one method of analysis may reveal aspects of cardiovascular dynamics that go undetected with the other method. This is a promising perspective for designing new clinical tools based on multifractal-multiscale analysis that, used in addition to traditional frequency-domain and time-domain methods, might substantially increase the risk stratification in the healthy population or improve the prediction of adverse events in cardiac patients.



---

*Source: 4801924-2018-01-22.xml*