1 Introduction

During pregnancy, the autonomic nervous system (ANS) undergoes significant changes to accommodate the physiological demands of the mother and support the developing fetus. The ANS is responsible for regulating various involuntary processes in the body, such as heart rate, blood pressure, digestion, and respiration (1). For instance, some of these changes are that heart rate tends to increase gradually during pregnancy, mainly due to the influence of hormonal changes and increased blood volume. Blood pressure also undergoes alterations, with a decrease in systemic vascular resistance and a subsequent reduction in blood pressure during the first two trimesters, followed by a gradual return to pre-pregnancy levels in the third trimester. One of the key adaptations in the ANS during pregnancy is an overall increase in sympathetic activity (2). This shift in autonomic balance helps to create an environment conducive to fetal growth and development. In addition, other systems related to the ANS, such as the respiratory system, are also affected during pregnancy, leading to higher minute ventilation (the total amount of air inhaled and exhaled per minute) (3). However, certain pregnancy complications, such as hypertensive disorders of pregnancy (HDP) or gestational diabetes mellitus (GDM), can disrupt the normal autonomic regulation, generating complications in both maternal and fetal health (4–7).

Heart rate variability (HRV) refers to the variation in time intervals between consecutive heartbeats (8), and it is considered a valuable measure of the ANS activity (9). Self-affinity, in the context of HRV, refers to the presence of statistical similarity or scaling properties across different time scales within the HRV signal. When HRV exhibits self-affinity, it means that the fluctuations in the time intervals between heartbeats display similar patterns or characteristics at different time scales. This property implies that the HRV signal possesses a certain level of complexity and organization (10, 11). For example, healthy individuals often exhibit statistical similarity patterns in their HRV, indicating a balanced and adaptable autonomic control of the heart rate. Reduced or altered self-affinity may indicate impaired autonomic regulation and increased risk of cardiovascular diseases, such as hypertension, heart failure, and arrhythmias. It can also be indicative of physiological stress, inflammation, or other disturbances in the body (12–15). The self-affinity analysis of HRV is commonly performed using techniques such as detrended fluctuation analysis (DFA) or fractal analysis (16). These methods quantify the presence and degree of self-affinity by examining how the fluctuations in HRV persist or decay across different time scales. Self-affinity in HRV has been studied in the context of pregnancy to understand the changes in ANS regulation in short-time recordings (17, 18), which suggest changes in the fractal-like behavior and complexity of HRV patterns during different stages of pregnancy. In general, the scaling exponents of HRV increase as gestation time advances, while the HRV complexity decreases (19). However, the reason is unclear since it may be influenced by hormonal changes, increased sympathetic activity, and other physiological adjustments related to pregnancy, given fetal development and growth.

As mentioned above, there are studies characterizing changes in HRV during pregnancy and the postpartum period, but they have been limited to short-term recordings of electrocardiographic (ECG) signals during wakefulness. Few studies have analyzed classic temporal and frequency indices of HRV during sleep at different stages of pregnancy (8, 20). However, sleep deserves more attention because it is a fundamental process that plays an essential role in maintaining overall health and well-being, and its importance is heightened during pregnancy. Sleep acts as a catalyst for immune function, facilitates tissue repair, and increases overall vitality. In addition, the emotional benefits of proper sleep are remarkable, acting as a shield against stress, anxiety, and mood swings - factors that are particularly prevalent during pregnancy (21). The effects of sleep extend to the growth and development of the fetus, as the body releases essential hormones that play a key role in the proper development of the fetus. An often underestimated aspect of sleep is its role in maintaining healthy blood pressure levels. This is particularly important during pregnancy when maintaining optimal blood pressure is essential to prevent complications such as pre-eclampsia. Sleep also plays a role in managing blood sugar levels and promoting healthy weight gain. Finally, sleep is a complex process involving many sub-processes during the different sleep stages, so it would be useful to establish indices to characterize it. Therefore, understanding the relationship between HRV, self-affinity, pregnancy, and sleep may provide insights into maternal cardiovascular health and sleep-related complications during this crucial period.

Sleep is a cyclical process characterized by multiple stages that influence the oscillatory components of HRV. In particular, the high-frequency (between 0.15 Hz and 0.5 Hz) and low-frequency (between 0.04 Hz and 0.15 Hz) components contribute to the self-affinity of cardiac activity, which is often described by a monofractal structure with 1/f behavior (15, 22). Additionally, some studies have shown that HRV also shows some characteristics of multifractality during wakefulness at different pathologies (23, 24), in 24h recordings (25), during sleep (26, 27) and in different stages of sleep (28). However, during sleep, the dominance of the high-frequency component could obscure valuable insights into the multifractality of cardiovascular dynamics during sleep, since slow and small variations with valuable information could be hidden by this oscillatory rhythm (29, 30). To address this issue and accurately characterize a possible multifractal self-affinity inherent to the sleep process, it would be interesting to assess variations in HRV amplitude after mitigating the impact of high-frequency components.

Thus, we are interested in analyzing sleep HRV patterns during pregnancy using classical and Multifractal indices to better understand the dynamics of the ANS regulation during sleep in different physiological conditions (early-pregnancy, middle-pregnancy, and non-pregnancy). The main aim is a better understanding of the underlying mechanisms: How does the autonomic nervous system controlling heart rate is modified during the different stages of pregnancy and compared to non-pregnancy?

The objective of this research is to study the self-affinity patterns exhibited by the ANS throughout the HRV indices during sleep in three distinct women groups: non-pregnancy, early-pregnancy, and middle-pregnancy. The analysis of self-affinity behavior is carried out through the widely recognized Multifractal Detrended Fluctuation Analysis methodology. Additionally, the sleep period is segmented into four intervals, and classical temporal and spectral indexes of heart rate variability are computed and subsequently compared across these segments. Lastly, we have analyzed the original feature space of the HRV indexes and its transformed version using Principal Component Analysis (PCA) with the aim of highlighting differences between non-pregnancy, early-pregnancy, and middle-pregnancy individuals.

2 Methodology

Figure 1 illustrates the procedural workflow used in this study. The process begins with the ECG signal, in which specific segments relevant to our analysis were identified and selected. This is followed by a precise detection of the R peaks. From this point, the process is divided into two distinct analyses. The first analysis focuses on the multifractal assessment, which is based on the detection of maxima within the RR intervals. In parallel, the second analysis focuses on the more traditional HRV analysis together with maxima within the RR intervals, including classification (class separation) and statistical evaluations.

Figure 1. Schematic representation of the methodology process.

2.1 Polysomnography data

Twenty-one whole-night polysomnographies (PSG) from Healthy (n = 6) and Nocturnal Frontal Lobe Epilepsy (NFLE, n = 15) women in reproductive age were downloaded from the CAP database. The database consists of 108 PSG recordings in EDF format acquired at the Sleep Disorders Center of the Ospedale Maggiore of Parma, Italy. This database is freely available from Physionet (http://www.physionet.org). The sleep recordings present at least three electroencephalography (EEG) channels, electrooculography (EOG), electromyography (EMG) of the submentalis muscle, bilateral anterior tibial EMG, respiratory signal, and ECG (sampling frequencies varied between 100 Hz and 512 Hz). The sleep stages were scored by expert neurologists following the gold standard rules (31). The complete database comprises the following groups: healthy (n = 16), nocturnal frontal lobe epilepsy (n = 40), REM Behavior Disorder (n = 22), periodic leg movement (n = 10), insomniac (n = 9), narcoleptic (n = 5), sleep disorder breathing (n = 4), and bruxism (n = 2). We selected recordings of women from the healthy and NFLE groups since the sleep profile was similar and the HRV was not altered by pathological events such as apneas. Besides, this similarity between healthy and NFLE subjects has been previously documented in the literature (32).

Eighteen whole-night PSGs of early-pregnant and eighteen ones of middle-pregnant nulliparous women were downloaded from the National Sleep Research Resource (NSRR) repository (funded by the National Heart, Lung, and Blood Institute from USA). The recordings belong to the study, called The Nulliparous Pregnancy Outcomes Study: Monitoring Mothers-to-Be (nuMoM2b) (33, 34). The goal of the database was to provide valuable information to healthcare professionals and women who are pregnant or planning to become pregnant and to provide data and knowledge for future research. Most participants underwent two PSGs, the first between weeks 6 and 15 (3072 records for the first trimester), and the second between weeks 22 and 29 (2805 records for the second trimester). At each of the two recordings, data were collected through clinical measurements, personal interviews, and self-administered questionnaires. Participants ranged in age from 14 to 44 years. Each PSG has airflow (AF), snoring, nasal pressure, abdominal and thoracic effort, blood oxygen saturation (SpO2), and ECG (sampled at 200 Hz). Each PSG comes with an XML file with wake-sleep annotation and the time of occurrence of events: central apnea, obstructive apnea, hypopnea, SpO2 desaturation, and noise in SpO2. For the signal acquisition the Embletta Gold digital recorder was used and the technical details can be found in the manual of operation (33, 34). Table 1 shows the demographic information of the women who were recorded. Information is presented as median and min-max interval. NFLE and healthy recordings are shown as non-pregnancy.

Table 1. Median and interval of the demographic information of the women making up the study groups.

It is important to note that the data come from two different data sets, which could cause inhomogeneity in the data. However, the acquisition protocol follows the standards accepted for clinical and home monitoring of sleep, with sampling frequencies for the ECG sufficient to correctly detect the R-peaks (between 100 and 512 Hz). In addition, it is worth mentioning that the HRV signal extraction only depends on the identification of the R peak on the ECG, which is possible even in presence of noise or on different ECG leads, thus it presents the advantage of being one of the physiological signals more robust to the noise and device independent in different data sets. In addition we tested and compared the HRV metrics from the healthy and NFLE groups through the Mann-Whitney U test and Brown-Forsythe test: no one of the HRV metrics resulted different between the two groups confirming their homogeneity.

As a final note regarding the non-pregnant group, despite a thorough search of a database of sleep recordings from women, we did not find a free database with healthy, nulliparous, young women to match the nuMoM2b dataset. To overcome this situation for future research, we are in the process of formulating a protocol to acquire recordings for a control group, with the goal of complementing the nuMoM2b dataset.

2.2 Pre-processing of the ECG signal

The pre-processing stage helps to improve the quality of the ECG signal, reduces noise and artifacts, and prepares the signal for further analysis, such as HRV assessment or feature extraction for clinical diagnosis. However, motion artifacts contaminating the ECG signal and missing segments of ECG are common in sleep recordings, so their identification is necessary before any other procedure.

2.2.1 Non-useful segments and R-peak detection

From each PSG, the ECG signal was extracted and used to locate the inter-beat time by detecting the R-peaks. R-peak detection is an important step since self-affinity analysis will be based on the occurrence of R-peaks. Therefore, it is necessary to define the ECG segments in which R-peaks can be reliably detected.

The amplitude of the ECG signal depends on the subject and varies throughout the night due to body movements caused by changes in body position, apneas, etc. Therefore, it is necessary to normalize the ECG to mitigate the effect of ECG amplitude change on the development of non-useful segment detectors. For this purpose, non-overlapping 10 s segments of the ECG signal were normalized with respect to their variance. Then, to detect signal segments with motion artifacts (non-useful segments), the trend was removed from each segment, and its variance was calculated. If the variance was between a value of 0.04 and 0.2 (note that segments were previously normalized with respect to the variance), the segment was considered ECG; otherwise, it was annotated as a non-useful segment. The optimal threshold interval was based on the analysis of manually annotated segments from ten polysomnographic sleep ECG signals. For all the recordings, the first and last 20 min were considered non-useful segments.

We then used two different approaches for automatic R-peak detection: the widely recognized Pan-Tompkins algorithm (35) and a custom procedure based on wavelet decomposition. This dual-method strategy was adopted because of the variance in signal characteristics, which warranted the use of the most effective technique for detecting accurate R-peaks in each scenario.

Following this step, we proceeded to eliminate non-useful segments and wake epochs from the R-peak series. The segments remaining after filtering were then merged. It’s worth noting that this process has minimal impact on the assessment of the scaling behavior (36). To ensure the robustness of our analysis, we only considered recordings with a minimum duration of 5 h for the calculation of self-affinity, in line with maintaining at least 3 sleep NREM-REM cycles.

2.2.2 Heart rate variability fluctuations

After obtaining the R-peaks series from the previous step, the consecutive R-peak intervals, often referred to as RR intervals or inter-beat intervals, were calculated. These RR intervals were then subjected to cubic spline interpolation to create a regularly spaced time series with a 4 Hz sampling rate.

The RR intervals exhibit fluctuation around a mean value, representing the mean Heart Rate. For extended time periods like sleep, these RR intervals tend to form a Gaussian distribution. However, for shorter time spans, variations in the properties of the global Gaussian distribution can be observed. These variations might correlate with sleep stages or sequences of pathological events. Therefore, an appropriate method is needed to account for these local variations over time. Multifractal Detrended Fluctuation Analysis (MF-DFA) is particularly well-suited for this purpose because it extends beyond the limitations of the standard DFA by considering the scaling of qth-order moments at different time scales (37, 38). Unlike standard DFA, which specifically focuses on the second moment (q=2), MF-DFA considers different qth-order moments, providing a more comprehensive perspective on the multifractal nature of the data. The latter approach allows for identifying or detecting intricate correlations and patterns that may remain unnoticed using other methods (27).

However, during sleep, the prevalence of parasympathetic activity, combined with the significant effect of breathing on the cardiovascular system, enhances the oscillatory rhythm within the RR intervals. This rhythm contains remarkable information on the HRV dynamics, but could also prevent some techniques from revealing other interesting features in the RR intervals (29, 30). Therefore, a time series was constructed to counterbalance the influence of this oscillatory rhythm. This series includes only the maxima of the RR intervals (here called MM fluctuations). This time series could also be extracted with various analytical methods, such as the Hilbert transform (29), empirical mode decomposition, wavelets, filtering techniques, or adaptive filtering.

Figure 2 shows in the upper panel the RR intervals (black line) during the sleep of an early-pregnancy record and the overlapping amplitude of the MM fluctuations (gray line). As we can easily see, the MM fluctuations contain information related to the slow oscillations of the RR intervals. A zoom of the yellow window is shown in the lower panel, where the MM fluctuations are depicted with a dashed line. MM fluctuations were used to assess the multifractal characteristics of the cardiovascular system during the sleep phase using MF-DFA.

Figure 2. Example of the RR intervals during sleep and along with the superimposed MM fluctuations (gray line). The lower panel offers a zoomed-in view of the highlighted yellow window, showing the MM fluctuations (dashed line).

2.3 Multifractal detrended fluctuation analysis (MF-DFA)

MF-DFA is a technique used to assess the presence of fractal scaling behavior in physiological signals such as HRV. The MF-DFA method involves the following steps:

1. Obtain a time series: The time series can be a sequence of evenly spaced time points, n, representing measurements over time, x[n].

2. Calculate the cumulative sum of the time series: This step involves integrating the time series minus the mean time series,

y[n]=∑k=1n(x[k]−E(x[n])),n=1,…,N,(1)

where N is the number of data points and E(x[n]) is the mean.

3. Divide the cumulative sum, y[n], into non-overlapping segments: The time series is divided into smaller segments of equal length. Ns=[N/s], segment of length s.

4. Fit a polynomial function to each segment: A polynomial function, typically a straight line, is fitted to each segment of the cumulative sum, ys[n]).

5. Calculate the q-th order fluctuation function: The q-th order deviations of the data points from the corresponding polynomial fit are calculated for each segment,

Fq(s)=(1N∑n=1N(y[n]−ys[n])q)1q,q∈R(2)

If the time series under analysis have fractal properties, the Fq(s) should follow a power-law of the form:

Fq(s)∼sH(q),(3)

where H(q) can be conceptualized as an extension of the Hurst exponent, with the equivalence H(2)≡α of the DFA exponent. If H(q) remains constant regardless of changes in q, the studied time series exhibits monofractal characteristics. Conversely, if H(q) varies with q, indicating a diverse scaling behavior, the time series possesses multifractal attributes. When q assumes positive values, H(q) represents the scaling behavior of the segments with large fluctuations. Conversely, for negative q values, H(q) characterizes the scaling dynamics of segments with small fluctuations.

The conventional approach to present the results of MF-DFA involves a two-step conversion process. First, H(q) is transformed into the q-order mass exponent (τ(q)), and subsequently, τ(q) is further transformed into the q-order singularity exponent (h) and the q-order singularity dimension (D(h)) (37, 38). This process finishes in the creation of a graphical representation known as the multifractal spectrum, where D(h) is plotted against h. This is:

τ(q)=qH(q)−1,(4)

h=dτ(q)dq,(5)

D(h)=qh−τ(q).(6)

2.4 Time and frequency domain indices of HRV

We calculated temporal and frequency indices of HRV to establish a benchmark in line with current research standards. These indices serve as important reference points, as there is a well-established link between these indices and the function of the parasympathetic and sympathetic branches of the ANS. The sleep time was partitioned into four non-overlapping segments, named Q1, Q2, Q3 and Q4. This segmentation was done to effectively capture the variations that occur during different moments of sleep. We adopted this strategy in order to overcome the usual sleep staging based on EEG and provide parameters able to quantify sleep characteristics in a simpler way along the night. We then calculated the indices separately for each of these segments.

Time domain indices:

• meanRR, representing the average of the RR intervals.

• sdRR, which measures the deviation of the RR intervals with respect to meanRR.

• rmsRR, which calculates the root mean square of the time differences between consecutive normal heartbeats.

• meanMM, which represents the average MM fluctuation; it can be correlated with respiratory behavior.

• sdMM, which measures the deviation of the MM fluctuations with respect to meanMM.

Frequency domain indices:

The Power Spectral Density (PSD) was computed using the Welch method. Specifically, the signal was segmented into eight parts with a 50% overlap, and each segment was windowed using a Hamming window. The periodograms of these segments were then averaged to obtain the final PSD estimate, from which the following indices were calculated.

• Low frequency (LF), the power within the frequency band from 0.04 Hz to 0.15 Hz.

• High frequency (HF), the power within the frequency band from 0.15 Hz to 0.5 Hz; it is particularly associated with parasympathetic activity and respiration.

• LF/HF ratio, a traditional indicator of sympathovagal balance that provides insight into the balance between sympathetic and parasympathetic activity. However, the context of the recording must be considered for proper interpretation and, in particular the possible shift of the HF component towards the LF frequency range, may lead to erroneous interpretation; On the other hand, the LF/HF ratio was frequently used in sleep studies where the respiration frequency is more regular making the interpretation more robust (39, 40).

To ensure a meaningful comparison between different sleep segments and participants, we normalized the power of the PSD that belongs to the HF and LF bands with respect to the power of both the low and high-frequency bands. In essence, the summation of the power in the low and high-frequency bands approaches unity.

2.5 Multifractal analysis of the MM fluctuations

The multifractal spectrum was derived from the MM fluctuations of the RR intervals. This was done using the Multifractal Detrended Fluctuation Analysis (MF-DFA) method with linear detrending. The analysis included a spectrum of scales ranging from 16 to 2,048 and involved the evaluation of different q-orders, specifically [−5,−3,−1,0,1,3,5]. From the multifractal spectrum, we extracted the multifractal indices, including the position arg⁡maxhD(h) and the width of D(h) shown in Figure 4.

Generating surrogate data is a common technique used in time series analysis to test hypotheses and assess the significance of observed patterns by comparing them to randomized versions of the original data. There are several methods (random shuffle, phase randomization, iterative amplitude adjusted Fourier transform, etc.) for generating surrogate data, each designed to preserve certain characteristics of the original time series while destroying the ones to be tested. In this study, the random shuffle method was chosen, which consists of randomly shuffling the original data points to destroy any temporal correlation in the signal while preserving the distribution of the values (mean value, variance). Thus, this is useful to evaluate the presence of temporal structure or autocorrelation in the data, both in the short and in the long time scales. In general, the procedure involves: (a) taking the original time series, and (b) randomly permute the order of the data points. This process could be repeated many times to generate multiple surrogate data sets. When the surrogated data show lower correlation with respect to the original one we may conclude that measured correlations are a feature characterizing the signal.

Thus, to assess the type of multifractality, we calculated surrogate MM fluctuations. This involved analyzing the MM fluctuations that were computed from the same RR intervals but in different order. In the shuffling procedure, the RR interval values were randomly reordered, which effectively destroyed all correlations related to the sleep stages. The MM fluctuations were then calculated and resampled at 4 Hz.

2.6 Statistics and classification

We calculated temporal and spectral indices of HRV within different sleep segments, each corresponding to one-quarter of the total sleep duration (designated Q1, Q2, Q3, and Q4). We then performed a two-way mixed Analysis Of Variance (ANOVA) for each HRV index, considering PREGNANCY as a three-level between-subject factor (non-pregnancy, early-pregnancy, and middle-pregnancy) and SLEEPQ as a within-subject factor with four levels (Q1, Q2, Q3, and Q4). For each index, we first analyzed the interaction PREGNANCY*SLEEPQ; if the interaction was found significant, we analyzed simple main effects. Conversely, in the absence of a significant interaction, the main effects of PREGNANCY and SLEEPQ were evaluated separately. In case of significant differences from the above evaluations, pairwise multiple comparisons (i.e., post-hoc tests) were made using the Bonferroni correction to account for the increasing type-I error due to multiple testing. The effect of PREGNANCY was also analyzed on the multifractal indices derived from the MM fluctuations by means of a one-way ANOVA, followed by the appropriate Bonferroni-corrected post hoc analysis in case significant differences were found. The analyzed data distributions were checked for approximate normality and homogeneity of variance using the Shapiro-Wilk and Levene tests, respectively. The sphericity assumption for the mixed ANOVAs was tested through the Mauchly test; if not met, the Greenhouse-Geisser correction was applied to the degrees of freedom of the within-subject analysis. A p-value less than 0.05 was considered statistically significant for every analysis. The presence of outliers that may introduce severe bias in the results was established by inspecting boxplots of data grouped by PREGNANCY and SLEEPQ. In general, subjects were considered outliers if their data points appeared distant from the range [Qlow−1.5⋅IQR,Qup+1.5⋅IQR], with IQR, Qlow, and Qup indicating the interquartile range, the lower, and upper quartiles, respectively. Only the most distant outliers were excluded from the analysis, and no more than one subject per group of PREGNANCY (i.e., less than 6% of the original data). The statistical analysis was carried out in SPSS Statistics (IBM, U.S.A.).

2.6.1 Classification

In the field of biomedical engineering, it is of great importance to determine whether information derived from physiological signals contributes to the automatic detection of medical conditions. Given the pressing need to explore novel features that can improve the classification process across different states in women, it is useful to investigate whether HRV indices, even in the absence of sleep stage annotations, can serve as discriminative features.

To understand this, we evaluated the separability among groups by plotting the following indices from the different sleep quarters:

(a) three of the HRV indices that were statistically different among the three groups.

(b) the first three scores obtained from the PCA applied to the whole set of HRV indices.

PCA is a method used for dimensionality reduction. The primary objective is to transform the original variables into a new set of variables, the principal components, which are orthogonal and capture the maximum variance in the data.

Given a dataset X consisting of n observations and p variables, the steps involved in PCA are:

1. Compute the Covariance Matrix:

Σ=1n−1(X−X¯)T(X−X¯)(7)

Here, X¯ is the mean-centered data matrix.

2. Eigenvalue Decomposition: Solve the eigenvalue problem for Σ:

Σv=λv(8)

Where λ are the eigenvalues and v are the corresponding eigenvectors.

3. Order and Select Principal Components: Arrange the eigenvalues in descending order and choose the first k eigenvectors corresponding to the largest k eigenvalues to represent the k principal components.

4. Compute Scores: The scores represent the projection of the original data onto the principal components. For the ith observation xi, the scores zij for the jth principal component is given by:

zij=xiTvj(9)

Here, vj is the jth eigenvector.

5. Transform Data: Transform the original data X into the new coordinate system using the selected principal components:

Y=XVk(10)

Where Vk contains the first k eigenvectors.

In this context, the scores provide a representation of each observation in terms of the principal components. They capture the essential information from the original data in a reduced dimensionality space.

3 Results

This section presents the results of the statistical analysis conducted on heart rate fluctuations across the three PREGNANCY groups and the four sleep quarters (SLEEPQ). These results include the multifractal analysis combined with an examination of the surrogate MM time series. Subsequently, classical HRV indices and their temporal evolution during sleep are presented. Finally, the separability of the three groups in the original feature space of the HRV indices and after PCA is investigated.

3.1 Multifractal analysis of MM fluctuations

Figure 3 provides a representative illustration of the results obtained by applying the MF-DFA method, with q-orders = root-mean-square (RMS), to the RR intervals extracted from a full night recording of an early-pregnant subject (please see Figure 2 for the original RR intervals). These results (Figures 3A,B) are contrasted with those obtained after shuffling the same RR intervals (Figures 3C,D). When the actual sequence of RR intervals is considered, the regression lines of the different q-orders show a similar slope (Figure 3A, time series with a structure close to 1/f noise), which does not allow observing the multifractality of time series. Thus, as shown in Figure 3B, the multifractal spectrum is concentrated within a small interval of h close to one. Figure 3C shows the regression lines of the MF-DFA of the shuffled RR intervals also at different q-orders = RMS. Again, all slopes of the regression lines are close to h≈0.5, indicating a random white noise. The Figure 3D shows multifractal spectrum small D(h) width with singular exponents h concentrated around 0.5 and multifractal dimension close to 1.

Figure 3. Example of MF-DFA applied to RR intervals during sleep in an early-pregnancy recording. The figure includes plots of Fq(s) for different q-orders RMS and the corresponding regression for both (A) the original RR intervals and (C) the shuffled RR intervals. (B) shows the corresponding multifractal spectrum for the original RR intervals, while (D) shows the multifractal spectrum for the shuffled RR intervals.

Figure 4 shows the results of applying the MF-DFA method to the MM fluctuations obtained from the same sleep recording (see Figure 2 for the original MM fluctuations). Again, the MF-DFA was performed for q-orders = RMS, and the results are shown in Figure 4A. This time, the regression lines of the different q-orders do not show a similar slope, which implies that the scaling behavior is different depending on the q-order, around the 1/f noise. The multifractal spectrum spans over a wide interval of h (see Figure 4B). On the other hand, the regression lines of the MF-DFA of the surrogate MM intervals are shown in Figure 3C. The slopes of these regression lines are around h≈0.5, indicating that the surrogate series tends to have a random white noise behavior. As can be seen in Figure 3D, the multifractal spectrum is centered at 0.5, but the long-term correlations are not completely destroyed.